[Rsiena-commits] r267 - in pkg: RSiena RSiena/R RSiena/src/model RSienaTest RSienaTest/src/model
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Mar 30 23:14:42 CEST 2014
Author: natalie
Date: 2014-03-30 23:14:41 +0200 (Sun, 30 Mar 2014)
New Revision: 267
Added:
pkg/RSiena/R/sienaRI.r
pkg/RSiena/R/sienaRIDynamics.r
Modified:
pkg/RSiena/DESCRIPTION
pkg/RSiena/changeLog
pkg/RSiena/src/model/Model.cpp
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/changeLog
pkg/RSienaTest/src/model/Model.cpp
Log:
Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION 2014-03-18 18:21:38 UTC (rev 266)
+++ pkg/RSiena/DESCRIPTION 2014-03-30 21:14:41 UTC (rev 267)
@@ -1,8 +1,8 @@
Package: RSiena
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.1-266
-Date: 2014-03-17
+Version: 1.1-267
+Date: 2014-03-30
Author: Ruth Ripley, Krists Boitmanis, Tom A.B. Snijders
Depends: R (>= 2.15.0)
Imports: Matrix
Added: pkg/RSiena/R/sienaRI.r
===================================================================
--- pkg/RSiena/R/sienaRI.r (rev 0)
+++ pkg/RSiena/R/sienaRI.r 2014-03-30 21:14:41 UTC (rev 267)
@@ -0,0 +1,551 @@
+#/******************************************************************************
+# * SIENA: Simulation Investigation for Empirical Network Analysis
+# *
+# * Web: http://www.stats.ox.ac.uk/~snijders/siena/
+# *
+# * File: sienaRI.r
+# *
+# * Description: Used to determine, print, and plots relative importances of effects
+# * in for potential desicions of actors at observation moments.
+# *****************************************************************************/
+
+##@sienaRI. Use as RSiena:::sienaRI()
+sienaRI <- function(data, ans=NULL, theta=NULL, algorithm=NULL, effects=NULL)
+{
+ if (!inherits(data, "siena"))
+ {
+ stop("no a legitimate Siena data specification")
+ }
+ if(!is.null(ans))
+ {
+ if (!inherits(ans, "sienaFit"))
+ {
+ stop(paste("ans is not a legitimate Siena fit object", sep=""))
+ }
+ if(!is.null(algorithm)||!is.null(theta)||!is.null(effects))
+ {
+ warning(paste("some information are multiply defined \n results will be based on 'theta', 'algorithm', and 'effects' stored in 'ans' (as 'ans$theta', 'ans$x', 'ans$effects')", sep=""))
+ }
+ if(sum(ans$effects$include==TRUE & (ans$effects$type =="endow"|ans$effects$type =="creation")) > 0){
+ stop("sienaRI does not yet work for models that contain endowment or creation effects")
+ }
+ contributions <- getChangeContributions(algorithm = ans$x, data = data, effects = ans$effects)
+ RI <- expectedRelativeImportance(conts = contributions, effects = ans$effects, theta =ans$theta)
+ }else{
+ if (!inherits(algorithm, "sienaAlgorithm"))
+ {
+ stop(paste("algorithm is not a legitimate Siena algorithm specification", sep=""))
+ }
+ algo <- algorithm
+ if (!inherits(effects, "sienaEffects"))
+ {
+ stop(paste("effects is not a legitimate Siena effects object", sep=""))
+ }
+ if(sum(effects$include==TRUE & (effects$type =="endow"|effects$type =="creation")) > 0)
+ {
+ stop("sienaRI does not yet work for models containinf endowment or creation effects")
+ }
+ effs <- effects
+ if (!is.numeric(theta))
+ {
+ stop("theta is not a legitimate parameter vector")
+ }
+ if(length(theta) != sum(effs$include==TRUE & effs$type!="rate"))
+ {
+ if(length(theta) != sum(effs$include==TRUE))
+ {
+ stop("theta is not a legitimate parameter vector \n number of parameters has to match number of effects")
+ }
+ warning("length of theta does not match the number of objective function effects\n theta is treated as if containing rate parameters")
+ paras <- theta
+ ## all necessary information available
+ ## call getChangeContributions
+ contributions <- getChangeContributions(algorithm = algo, data = data, effects = effs)
+ RI <- expectedRelativeImportance(conts = contributions, effects = effs, theta = paras)
+ }else{
+ paras <- theta
+ ## all necessary information available
+ ## call getChangeContributions
+ contributions <- getChangeContributions(algorithm = algo, data = data, effects = effs)
+ RI <- expectedRelativeImportance(conts = contributions, effects = effs, theta = paras)
+ }
+ }
+ RI
+}
+
+##@getChangeContributions. Use as RSiena:::getChangeContributions
+getChangeContributions <- function(algorithm, data, effects)
+{
+ ## The following initializations data, effects, and model
+ ## for calling "getTargets" in "siena07.setup.h"
+ ## is more or less copied from "getTargets" in "getTargets.r".
+ ## However, some modifications have been necessary to get it to work.
+ f <- unpackData(data,algorithm)
+
+ effects <- effects[effects$include,]
+ if (!is.null(algorithm$settings))
+ {
+ effects <- addSettingsEffects(effects, algorithm)
+ }
+ else
+ {
+ effects$setting <- rep("", nrow(effects))
+ }
+ pData <- .Call('setupData', PACKAGE=pkgname,
+ list(as.integer(f$observations)),
+ list(f$nodeSets))
+ ## register a finalizer
+ ans <- reg.finalizer(pData, clearData, onexit = FALSE)
+ ans<- .Call('OneMode', PACKAGE=pkgname,
+ pData, list(f$nets))
+ ans<- .Call('Behavior', PACKAGE=pkgname, pData,
+ list(f$behavs))
+ ans<-.Call('ConstantCovariates', PACKAGE=pkgname,
+ pData, list(f$cCovars))
+ ans<-.Call('ChangingCovariates',PACKAGE=pkgname,
+ pData,list(f$vCovars))
+ ans<-.Call('DyadicCovariates',PACKAGE=pkgname,
+ pData,list(f$dycCovars))
+ ans<-.Call('ChangingDyadicCovariates',PACKAGE=pkgname,
+ pData, list(f$dyvCovars))
+
+ storage.mode(effects$parm) <- 'integer'
+ storage.mode(effects$group) <- 'integer'
+ storage.mode(effects$period) <- 'integer'
+
+ effects$effectPtr <- rep(NA, nrow(effects))
+ depvarnames <- names(data$depvars)
+ tmpeffects <- split(effects, effects$name)
+ myeffectsOrder <- match(depvarnames, names(tmpeffects))
+ ans <- .Call("effects", PACKAGE=pkgname, pData, tmpeffects)
+ pModel <- ans[[1]][[1]]
+ for (i in 1:length(ans[[2]]))
+ {
+ effectPtr <- ans[[2]][[i]]
+ tmpeffects[[i]]$effectPtr <- effectPtr
+ }
+ myeffects <- tmpeffects
+ for(i in 1:length(myeffectsOrder)){
+ myeffects[[i]]<-tmpeffects[[myeffectsOrder[i]]]
+ }
+ ans <- .Call("getTargets", PACKAGE=pkgname, pData, pModel, myeffects, parallelrun=TRUE, returnActorStatistics=FALSE, returnStaticChangeContributions=TRUE)
+ ans
+}
+
+expectedRelativeImportance <- function(conts, effects, theta, effectNames = NULL)
+{
+ waves <- length(conts[[1]])
+ effects <- effects[effects$include == TRUE,]
+ noRate <- effects$type != "rate"
+ effects <- effects[noRate,]
+ if(sum(noRate)!=length(theta))
+ {
+ theta <- theta[noRate]
+ }
+ effectNa <- attr(conts,"effectNames")
+ effectTypes <- attr(conts,"effectTypes")
+ networkNames <- attr(conts,"networkNames")
+ networkTypes <- attr(conts,"networkTypes")
+ networkInteraction <- effects$interaction1
+ effectIds <- paste(effectNa,effectTypes,networkInteraction, sep = ".")
+
+ currentDepName <- ""
+ depNumber <- 0
+ for(eff in 1:length(effectIds))
+ {
+ if(networkNames[eff] != currentDepName)
+ {
+ currentDepName <- networkNames[eff]
+ actors <- length(conts[[1]][[1]][[1]])
+ if(networkTypes[eff] == "oneMode")
+ {
+ choices <- actors
+ }else if(networkTypes[eff] == "behavior"){
+ choices <- 3
+ }else{
+ stop("so far, sienaRI works only for dependent variables of type 'oneMode' or 'behavior'")
+ }
+ depNumber <- depNumber + 1
+ currentDepEffs <- effects$name == currentDepName
+ effNumber <- sum(currentDepEffs)
+
+ RIs <- data.frame(row.names = effectIds[currentDepEffs])
+ RIs <- cbind(RIs, matrix(0, nrow=effNumber, ncol = actors))
+ entropies <- vector(mode="numeric", length = actors)
+
+ currentDepObjEffsNames <- paste(effects$shortName[currentDepEffs],effects$type[currentDepEffs],effects$interaction1[currentDepEffs],sep=".")
+ otherObjEffsNames <- paste(effects$shortName[!currentDepEffs],effects$type[!currentDepEffs],effects$interaction1[!currentDepEffs],sep=".")
+
+ expectedRI <- list()
+ RIActors <- list()
+ absoluteSumActors <- list()
+ entropyActors <-list()
+ for(w in 1:waves)
+ {
+ currentDepEffectContributions <- conts[[1]][[w]][currentDepEffs]
+ currentDepEffectContributions <- sapply(lapply(currentDepEffectContributions, unlist), matrix, nrow=actors, ncol=choices, byrow=TRUE, simplify="array")
+
+ distributions <- apply(apply(currentDepEffectContributions, c(2,1), as.matrix), 3, calculateDistributions, theta[which(currentDepEffs)])
+ distributions <- lapply(apply(distributions, 2, list), function(x){matrix(x[[1]], nrow=effNumber+1, ncol=choices, byrow=F)})
+
+ entropy_vector <- unlist(lapply(distributions,function(x){entropy(x[1,])}))
+ ## If one wishes another measure than the L^1-difference between distributions,
+ ## here is the right place to call some new function instead of "L1D".
+ RIs_list <- lapply(distributions,function(x){L1D(x[1,], x[2:dim(x)[1],])})
+ RIs_matrix <-(matrix(unlist(RIs_list),nrow=effNumber, ncol=actors, byrow=F))
+
+ RIs <- RIs_matrix
+ entropies <- entropy_vector
+
+ RIActors[[w]] <- apply(RIs, 2, function(x){x/sum(x)})
+ absoluteSumActors[[w]] <- colSums(RIs)
+ entropyActors[[w]] <- entropies
+ expectedRI[[w]] <- rowSums(RIActors[[w]] )/dim(RIActors[[w]])[2]
+ }
+ RItmp <- NULL
+ RItmp$dependentVariable <- currentDepName
+ RItmp$expectedRI <- expectedRI
+ RItmp$RIActors <- RIActors
+ RItmp$absoluteSumActors <- absoluteSumActors
+ RItmp$entropyActors <- entropyActors
+ if(!is.null(effectNames))
+ {
+ RItmp$effectNames <- effectNames[currentDepEffs]
+ }else{
+ RItmp$effectNames <- paste(effectTypes[currentDepEffs], " ", effects$effectName[currentDepEffs], sep="")
+ }
+ class(RItmp) <- "sienaRI"
+ if(depNumber == 1){
+ RI <- RItmp
+ }else if(depNumber == 2){
+ RItmp1 <- RI
+ RI <- list()
+ RI[[1]]<-RItmp1
+ RI[[2]]<-RItmp
+ }else{
+ RI[[depNumber]]<-RItmp
+ }
+ }
+ }
+ if(depNumber>1)
+ {
+ warning("more than one dependent variable\n return value is therefore not of class 'sienaRI'\n but a list of objects of class 'sienaRI'")
+ }
+ RI
+}
+
+calculateDistributions <- function(effectContributions = NULL, theta = NULL)
+{
+ effects <- dim(effectContributions)[1]
+ choices <- dim(effectContributions)[2]
+ effectContributions[effectContributions=="NaN"]<-0
+ distributions <- array(dim = c(effects+1,choices))
+ distributions[1,] <- exp(colSums(theta*effectContributions))/sum(exp(colSums(theta*effectContributions)))
+ for(eff in 1:effects)
+ {
+ t <- theta
+ t[eff] <- 0
+ distributions[eff+1,] <- exp(colSums(t*effectContributions))/sum(exp(colSums(t*effectContributions)))
+ }
+ distributions
+}
+
+entropy <- function(distribution = NULL)
+{
+ entropy <- -1*(distribution %*% log(distribution)/log(length(distribution)))
+ certainty <- 1-entropy
+ certainty
+}
+
+KLD <- function(referenz = NULL, distributions = NULL)
+{
+ if(is.vector(distributions))
+ {
+ kld <- (referenz %*% (log(referenz)-log(distributions)))/log(length(referenz))
+ }
+ else
+ {
+ kld <- colSums(referenz *(log(referenz)-t(log(distributions))))/log(length(referenz))
+ }
+ kld
+}
+
+## calculates the L^1-differenz between distribution "reference" (which is a vector of length n)
+## and each row of distributions (which is a matrix with n columns)
+L1D <- function(referenz = NULL, distributions = NULL)
+{
+ if(is.vector(distributions))
+ {
+ l1d <- sum(abs(referenz-distributions))
+ }
+ else
+ {
+ l1d <- colSums(abs(referenz-t(distributions)))
+ }
+ l1d
+}
+
+##@print.sienaRI Methods
+print.sienaRI <- function(x, ...){
+ if (!inherits(x, "sienaRI"))
+ {
+ stop("not a legitimate Siena relative importance of effects object")
+ }
+ cat(paste("\n Expected relative importance of effects for dependent variable '", x$dependentVariable,"' at observation moments:\n\n\n", sep=""))
+ waves <- length(x$expectedRI)
+ effs <- length(x$effectNames)
+ colNames = paste("wave ", 1:waves, sep="")
+ line1 <- paste(format("", width =63), sep="")
+ line2 <- paste(format(1:effs,width=3), '. ', format(x$effectNames, width = 56),sep="")
+ for(w in 1:length(colNames))
+ {
+ line1 <- paste(line1, format(colNames[w], width=8)," ", sep = "")
+ line2 <- paste(line2, format(round(x$expectedRI[[w]], 4), width=8, nsmall=4)," ",sep="")
+ }
+ line2 <- paste(line2, rep('\n',effs), sep="")
+ cat(as.matrix(line1),'\n \n', sep='')
+ cat(as.matrix(line2),'\n', sep='')
+ invisible(x)
+}
+
+##@summary.sienaRI Methods
+summary.sienaRI <- function(object, ...)
+{
+ if (!inherits(x, "sienaRI"))
+ {
+ stop("not a legitimate Siena relative importance of effects object")
+ }
+ class(x) <- c("summary.sienaRI", class(x))
+ x
+}
+##@print.summary.sienaRI Methods
+print.summary.sienaRI <- function(x, ...)
+{
+ if (!inherits(x, "summary.sienaRI"))
+ {
+ stop("not a legitimate summary of a Siena relative importance of effects object")
+ }
+ print.sienaRI(x)
+ invisible(x)
+}
+
+
+##@plot.sienaRI Methods
+plot.sienaRI <- function(x, file = NULL, col = NULL, addPieChart = FALSE, radius = 1, width = NULL, height = NULL, legend = TRUE, legendColumns = NULL, legendHeight = NULL, cex.legend = NULL, cex.names = NULL, ...)
+{
+ if (!inherits(x, "sienaRI"))
+ {
+ stop("not a legitimate Siena relative importance of effects object")
+ }
+ waves <- length(x$expectedRI)
+ actors <- dim(x$RIActors[[1]])[2]
+ if(legend)
+ {
+ if(!is.null(legendColumns))
+ {
+ if(is.numeric(legendColumns))
+ {
+ legendColumns <- as.integer(legendColumns)
+ }else{
+ legendColumns <- NULL
+ warning("legendColumns has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(legendColumns))
+ {
+ legendColumns <-floor((actors+2)/11)
+ }
+ if(!is.null(legendHeight))
+ {
+ if(is.numeric(legendHeight))
+ {
+ legendHeight <- legendHeight
+ }else{
+ legendHeight <- NULL
+ warning("legendHeight has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(legendHeight))
+ {
+ legendHeight <- max(0.6,ceiling(length(x$effectNames)/legendColumns)*0.2)
+ }
+ }
+ if(!is.null(height))
+ {
+ if(is.numeric(height))
+ {
+ height <- height
+ }else{
+ height <- NULL
+ warning("height has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(height))
+ {
+ if(legend)
+ {
+ height <- (2*waves+2*legendHeight)*1
+ }else{
+ height <- (2*waves)*1
+ }
+ }
+
+ if(!is.null(width))
+ {
+ if(is.numeric(width))
+ {
+ width <- width
+ }else{
+ width <- NULL
+ warning("width has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(width))
+ {
+ if(addPieChart)
+ {
+ width = (actors/3+2)*1
+ }else{
+ width = (actors/3+1)*1
+ }
+ }
+
+ if(!is.null(cex.legend))
+ {
+ if(is.numeric(cex.legend))
+ {
+ cex.legend <- cex.legend
+ }else{
+ cex.legend <- NULL
+ warning("cex.legend has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(cex.legend))
+ {
+ cex.legend <- 1.3
+ }
+
+ if(!is.null(cex.names))
+ {
+ if(is.numeric(cex.names))
+ {
+ cex.names <- cex.names
+ }else{
+ cex.names <- NULL
+ warning("cex.names has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(cex.names))
+ {
+ cex.names <- 1
+ }
+
+ if(!is.null(radius))
+ {
+ if(is.numeric(radius))
+ {
+ rad <- radius
+ }else{
+ rad <- NULL
+ warning("radius has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(radius))
+ {
+ rad <- 1
+ }
+
+ createPdf = FALSE
+ if(!is.null(file))
+ {
+ if (is.character(file)){
+ pdf(file = file, width = width, height = height)
+ createPdf = TRUE
+ }else{
+ file = NULL
+ warning("file has to be a of type 'character' \n could not create pdf")
+ }
+ }
+ if(is.null(file))
+ {
+ windows(width = width, height = height)
+ }
+
+ if(!is.null(col))
+ {
+ cl <- col
+ }else{
+ alph <- 175
+ green <- rgb(127, 201, 127,alph, maxColorValue = 255)
+ lila <-rgb(190, 174, 212,alph, maxColorValue = 255)
+ orange <- rgb(253, 192, 134,alph, maxColorValue = 255)
+ yellow <- rgb(255, 255, 153,alph, maxColorValue = 255)
+ blue <- rgb(56, 108, 176,alph, maxColorValue = 255)
+ lightgray <- rgb(184,184,184,alph, maxColorValue = 255)
+ darkgray <- rgb(56,56,56,alph, maxColorValue = 255)
+ gray <- rgb(120,120,120,alph, maxColorValue = 255)
+ pink <- rgb(240,2,127,alph, maxColorValue = 255)
+ brown <- rgb(191,91,23,alph, maxColorValue = 255)
+ cl <- c(green,lila,orange,yellow,blue,lightgray,darkgray,gray,pink,brown)
+ while(length(cl)<length(x$effectNames)){
+ alph <- (alph+75)%%255
+ green <- rgb(127, 201, 127,alph, maxColorValue = 255)
+ lila <-rgb(190, 174, 212,alph, maxColorValue = 255)
+ orange <- rgb(253, 192, 134,alph, maxColorValue = 255)
+ yellow <- rgb(255, 255, 153,alph, maxColorValue = 255)
+ blue <- rgb(56, 108, 176,alph, maxColorValue = 255)
+ lightgray <- rgb(184,184,184,alph, maxColorValue = 255)
+ darkgray <- rgb(56,56,56,alph, maxColorValue = 255)
+ gray <- rgb(120,120,120,alph, maxColorValue = 255)
+ pink <- rgb(240,2,127,alph, maxColorValue = 255)
+ brown <- rgb(191,91,23,alph, maxColorValue = 255)
+ cl <- c(cl,green,lila,orange,yellow,blue,lightgray,darkgray,gray,pink,brown)
+ }
+ }
+ bordergrey <-"gray25"
+
+ if(addPieChart)
+ {
+ if(legend)
+ {
+ layoutMatrix <- matrix(c(1:(2*waves+1),(2*waves+1)), byrow= TRUE, ncol=2, nrow=(waves+1))
+ layout(layoutMatrix,widths=c((actors/3),2),heights=c(rep(1,waves),legendHeight))
+ }else{
+ layoutMatrix <- matrix(c(1:(2*waves)), byrow= TRUE, ncol=2, nrow=waves)
+ layout(layoutMatrix,widths=c((actors/3),2),heights=rep(1,waves))
+ }
+ par( oma = c( 0, 0, 2, 0 ),mar = par()$mar+c(-1,0,-3,-2), xpd=T , cex = 0.75, no.readonly = TRUE )
+ }else{
+ if(legend)
+ {
+ layoutMatrix <- matrix(c(1:(waves+1)), byrow= TRUE, ncol=1, nrow=(waves+1))
+ layout(layoutMatrix,widths=c((actors/3)),heights=c(rep(1,waves),legendHeight))
+ }else{
+ layoutMatrix <- matrix(c(1:waves), byrow= TRUE, ncol=1, nrow=waves)
+ layout(layoutMatrix,widths=c((actors/3)),heights=rep(1,waves))
+ }
+ par( oma = c( 0, 0, 2, 0 ),mar = par()$mar+c(-1,0,-3,1), xpd=T , cex = 0.75, no.readonly = TRUE )
+ }
+ for(w in 1:waves)
+ {
+ barplot(cbind(x$RIActors[[w]], x$expectedRI[[w]]),space=c(rep(0.1,actors),1.5),width=c(rep(1,actors),1), beside =FALSE, yaxt = "n", xlab="Actor", cex.names = cex.names, ylab=paste("wave ", w, sep=""),border=bordergrey, col = cl, names.arg=c(1:actors,"exp. rel. imp."))
+ axis(2, at=c(0,0.25,0.5,0.75,1),labels=c("0","","0.5","","1"))
+ axis(4, at=c(0,0.25,0.5,0.75,1),labels=c("0","","0.5","","1"))
+ if(addPieChart)
+ {
+ pie(x$expectedRI[[w]], col = cl, labels=NA, border = bordergrey, radius = rad)
+ mtext("exp. rel. imp.",side = 1, line = 1, cex=cex.names*0.75)
+ }
+ }
+ if(legend)
+ {
+ plot(c(0,1), c(0,1), col=rgb(0,0,0,0),axes=FALSE, ylab = "", xlab = "")
+ legend(0, 1, x$effectNames, fill=cl, ncol = legendColumns, bty = "n", cex=cex.legend)
+ }
+ if(createPdf)
+ {
+ dev.off()
+ }
+ invisible(cl)
+}
+
Added: pkg/RSiena/R/sienaRIDynamics.r
===================================================================
--- pkg/RSiena/R/sienaRIDynamics.r (rev 0)
+++ pkg/RSiena/R/sienaRIDynamics.r 2014-03-30 21:14:41 UTC (rev 267)
@@ -0,0 +1,503 @@
+#/******************************************************************************
+# * SIENA: Simulation Investigation for Empirical Network Analysis
+# *
+# * Web: http://www.stats.ox.ac.uk/~snijders/siena/
+# *
+# * File: sienaRIDynamics.r
+# *
+# * Description: Used to determine, print, and plots relative importances of effects
+# * in sequences of simulated micro-steps
+# *****************************************************************************/
+
+##@sienaRIDynamics. Use as RSiena:::sienaRIDynamics()
+sienaRIDynamics <- function(data, ans=NULL, algorithm=NULL, effects=NULL, theta=NULL , depvar=NULL, intervalsPerPeriod=NULL)
+{
+ if(length(data$depvars)>1){
+ if(is.null(depvar)){
+ warning("If the models contains more than one dependent variables, \n it should be specified by variable 'depvar' for which dependent variable the relative importances should be calculated. \n\n")
+ warning(paste("As 'depvar = NULL', relative importances are calculated for variable ", currentNetName, sep=""))
+ currentNetName <- attr(data$depvars,"names")[1]
+ }else if(!(depvar %in% attr(data$depvars,"names"))){
+ stop("'depvar' is not a name of a dependent variable")
+ }else{
+ currentNetName <- depvar
+ }
+ }else{
+ currentNetName <- depvar
+ }
+ if (!inherits(data, "siena"))
+ {
+ stop("data is not a legitimate Siena data specification")
+ }
+ if(!is.null(ans))
+ {
+ if (!inherits(ans, "sienaFit"))
+ {
+ stop("ans is not a legitimate Siena fit object")
+ }
+ if(!is.null(algorithm)||!is.null(theta)||!is.null(effects))
+ {
+ warning("some information are multiply defined \n results will be based on 'theta', 'algorithm', and 'effects' stored in 'ans' (as 'ans$theta', 'ans$x', 'ans$effects')")
+ }
+ if(!is.null(intervalsPerPeriod))
+ {
+ if(is.numeric(intervalsPerPeriod))
+ {
+ intervalsPerPeriod <- as.integer(intervalsPerPeriods)
+ }else{
+ intervalsPerPeriod <- NULL
+ warning("'intervalsPerPeriod' has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(intervalsPerPeriod))
+ {
+ intervalsPerPeriod <- 10
+ }
+ RIValues <- calculateRIDynamics(data = data, theta= ans$theta, algorithm = ans$x, effects = ans$effects, depvar = currentNetName, intervalsPerPeriod=intervalsPerPeriod)
+ }else{
+ if (!inherits(algorithm, "sienaAlgorithm"))
+ {
+ stop("algorithm is not a legitimate Siena algorithm specification")
+ }
+ algo <- algorithm
+ if (!inherits(effects, "sienaEffects"))
+ {
+ stop("effects is not a legitimate Siena effects object")
+ }
+ effs <- effects
+ if(!is.numeric(theta))
+ {
+ stop("theta is not a legitimate parameter vector")
+ }
+ if(length(theta) != sum(effs$include==TRUE))
+ {
+ if(length(theta) == sum(effs$include==TRUE & effs$type!="rate"))
+ {
+ stop("vector of model parameters has wrong dimension, maybe rate parameters are missing")
+ }
+ stop("theta is not a legitimate parameter vector \n number of parameters has to match number of effects")
+ }
+ paras <- theta
+ if(!is.null(intervalsPerPeriod))
+ {
+ if(is.numeric(intervalsPerPeriod))
+ {
+ intervalsPerPeriod <- as.integer(intervalsPerPeriods)
+ }else{
+ intervalsPerPeriod <- NULL
+ warning("'intervalsPerPeriod' has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(intervalsPerPeriod))
+ {
+ intervalsPerPeriod <- 10
+ }
+ ## all necessary information available
+ RIValues <- calculateRIDynamics(data = data, theta= paras, algorithm = algo, effects = effs, depvar = currentNetName, intervalsPerPeriod=intervalsPerPeriod)
+ }
+ RIValues
+}
+
+##@calculateRIDynamics calculateRIDynamics simulates sequences of micro-steps, and aggregates the relative importances of effects over micro-steps of same time intervals
+calculateRIDynamics <- function(data=NULL, theta=NULL, algorithm=NULL, effects=NULL, depvar=NULL, intervalsPerPeriod=10, returnActorStatistics=NULL)
+{
+ x <- algorithm
+ currentNetName <- depvar
+ z <- NULL
+ z$FRAN <- getFromNamespace(x$FRANname, pkgname)
+ x$cconditional <- FALSE
+ z$print <- FALSE
+ z$Phase <- 3
+ z <- initializeFRAN(z, x, data, effects, prevAns=NULL, initC=FALSE, returnDeps=FALSE)
+ z$returnChangeContributions <- TRUE
+ z$theta <- theta
+ if (!is.null(x$randomSeed))
+ {
+ set.seed(x$randomSeed, kind="default")
+ }
+ else
+ {
+ if (exists(".Random.seed"))
+ {
+ rm(.Random.seed, pos=1)
+ RNGkind(kind="default")
+ }
+ }
+ chains <- x$n3
+ periods <- data$observation-1
+ effects <- effects[effects$include==TRUE,]
+ noRate <- effects$type != "rate"
+ thetaNoRate <- theta[noRate]
+ effectNames <- effects$shortName[noRate]
+ effectTypes <- effects$type[noRate]
+ networkName <- effects$name[noRate]
+ networkInteraction <- effects$interaction1[noRate]
+ effectIds <- paste(effectNames,effectTypes,networkInteraction, sep = ".")
+ currentNetObjEffs <- effects$name[noRate] == currentNetName
+ RIintervalValues <- list()
+ for(period in 1:periods)
+ {
+ RIintervalValues[[period]] <- data.frame(row.names = effectIds[currentNetObjEffs])
+ }
+
+ for (chain in (1:chains))
+ {
+ ans <- z$FRAN(z, x)
+ for(period in 1:periods)
+ {
+ microSteps <- length(ans$changeContributions[[1]][[period]])
+ stepsPerInterval <- microSteps/intervalsPerPeriod
+
+ RItmp <- data.frame(row.names = effectIds[currentNetObjEffs])
+ interval <- 1
+ stepsInIntervalCounter <- 0
+ for(microStep in 1:microSteps)
+ {
+ #currentNetName <- attr(ans$changeContributions[[1]][[period]][[microStep]],"networkName")
+ if(attr(ans$changeContributions[[1]][[period]][[microStep]],"networkName")==currentNetName){
+ stepsInIntervalCounter <- stepsInIntervalCounter + 1
+ ## distributions[1,] contains the probabilities of the available choices in this micr-step
+ ## distributions[2,],distributions[3,], ...contains the probabilities of the available choices
+ ## if the parameter of the second, third, ... effects is set to zero.
+ distributions <- calculateDistributions(ans$changeContributions[[1]][[period]][[microStep]], thetaNoRate[currentNetObjEffs])
+ ## If one wishes another measure than the L^1-difference between distributions,
+ ## here is the right place to call some new function instead of "L1D".
+ RItmp[,stepsInIntervalCounter] <- standardize(L1D(distributions[1,], distributions[2:dim(distributions)[1],]))
+ }
+ if (microStep > interval * stepsPerInterval || microStep == microSteps)
+ {
+ if(chain == 1)
+ {
+ RIintervalValues[[period]][,interval] <- rowSums(RItmp)/length(RItmp)
+ }
+ else
+ {
+ RIintervalValues[[period]][,interval] <- RIintervalValues[[period]][,interval] + rowSums(RItmp)/length(RItmp)
+ }
+ interval <- interval + 1
+ stepsInIntervalCounter <- 0
+ }
+ }
+ }
+ }
+ for(period in 1:periods)
+ {
+ RIintervalValues[[period]] <- RIintervalValues[[period]]/chains
+ }
+ RIDynamics <- NULL
+ RIDynamics$intervalValues <-RIintervalValues
+ RIDynamics$dependentVariable <- currentNetName
+ RIDynamics$effectNames <- paste(effectTypes[currentNetObjEffs], " ", (effects$effectName[noRate])[currentNetObjEffs], sep="")
+ class(RIDynamics) <- "sienaRIDynamics"
+ RIDynamics
+}
+
+standardize <- function(values = NULL)
+{
+ newValues <- values/sum(values)
+}
+
+
+##@print.sienaRIDynamics Methods
+print.sienaRIDynamics <- function(x, ...){
+ if (!inherits(x, "sienaRIDynamics"))
+ {
+ stop("x is not of class 'sienaRIDynamics' ")
+ }
+ periods <- length(x$intervalValues)
+ intervals <- length(x$intervalValues[[1]])
+ effs <- length(x$effectNames)
+ cat(paste("\n Relative importance of effects in micro-steps of dependent variable '", x$dependentVariable,"'. \n \n",sep=""))
+ cat(paste(" Periods between observations are devided into ", intervals, " intervals. \n \n",sep=""))
+ cat(paste(" Displayed results are aggregations over intervals:\n", sep=""))
+ for(p in 1:periods){
+ cat(paste("\n\nPeriod ",p,":\n", sep=""))
+ colNames = paste("int ", 1:intervals, sep="")
+ line1 <- paste(format("", width =52), sep="")
+ line2 <- paste(format(1:effs,width=3), '. ', format(x$effectNames, width = 45),sep="")
+ col <- 0
+ for(i in 1:length(colNames))
+ {
+ col <- col + 1
+ line1 <- paste(line1, format(colNames[i], width=8)," ", sep = "")
+ line2 <- paste(line2, format(round(x$intervalValues[[p]][[i]], 4), width=8, nsmall=4)," ",sep="")
+ if(col == 5)
+ {
+ line2 <- paste(line2, rep('\n',effs), sep="")
+ cat(as.matrix(line1),'\n \n', sep='')
+ cat(as.matrix(line2),'\n', sep='')
+ line1 <- paste(format("", width =52), sep="")
+ line2 <- paste(format(1:effs,width=3), '. ', format(x$effectNames, width = 45),sep="")
+ col<-0
+ }
+ }
+ if(col>0)
+ {
+ line2 <- paste(line2, rep('\n',effs), sep="")
+ cat(as.matrix(line1),'\n \n', sep='')
+ cat(as.matrix(line2),'\n', sep='')
+ }
+ }
+ invisible(x)
+}
+#
+##@summary.sienaRIDynamics Methods
+summary.sienaRIDynamics <- function(object, ...)
+{
+ if (!inherits(x, "sienaRIDynamics"))
+ {
+ stop("x is not of class 'sienaRIDynamics' ")
+ }
+ class(x) <- c("summary.sienaRIDynamics", class(x))
+ x
+}
+##@print.summary.sienaRIDynamics Methods
+print.summary.sienaRIDynamics <- function(x, ...)
+{
+ if (!inherits(x, "summary.sienaRIDynamics"))
+ {
+ stop("not a legitimate summary of a 'sienaRIDynamics' object")
+ }
+ print.sienaRIDynamics(x)
+ invisible(x)
+}
+
+
+##@plot.sienaRIDynamics Methods
+plot.sienaRIDynamics <- function(x, staticValues = NULL, file = NULL, col = NULL, ylim=NULL, width = NULL, height = NULL, legend = TRUE, legendColumns = NULL, legendHeight = NULL, cex.scale = NULL, cex.legend = NULL, cex.axis = NULL, cex.names = NULL, ylab = "", xlab = "", ...)
+{
+ if (!inherits(x, "sienaRIDynamics"))
+ {
+ stop("x is not of class 'sienaRIDynamics' ")
+ }
+ if(legend)
+ {
+ if(!is.null(legendColumns))
+ {
+ if(is.numeric(legendColumns))
+ {
+ legendColumns <- as.integer(legendColumns)
+ }else{
+ legendColumns <- NULL
+ warning("legendColumns has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(legendColumns))
+ {
+ legendColumns <- 3
+ }
+ if(!is.null(legendHeight))
+ {
+ if(is.numeric(legendHeight))
+ {
+ legendHeight <- legendHeight
+ }else{
+ legendHeight <- NULL
+ warning("legendHeight has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(legendHeight))
+ {
+ legendHeight <- 1
+ }
+ }
+ if(!is.null(height))
+ {
+ if(is.numeric(height))
+ {
+ height <- height
+ }else{
+ height <- NULL
+ warning("height has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(height))
+ {
+ if(legend)
+ {
+ height <- 4
+ }else{
+ height <- 3
+ }
+ }
+
+ if(!is.null(width))
+ {
+ if(is.numeric(width))
+ {
+ width <- width
+ }else{
+ width <- NULL
+ warning("width has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(width))
+ {
+ width <- 8
+ }
+
+ if(!is.null(cex.scale))
+ {
+ if(is.numeric(cex.scale))
+ {
+ cex.scale <- cex.scale
+ }else{
+ cex.scale <- NULL
+ warning("cex.scale has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(cex.scale))
+ {
+ cex.scale <- 1
+ }
+
+ if(!is.null(cex.legend))
+ {
+ if(is.numeric(cex.legend))
+ {
+ cex.legend <- cex.legend
+ }else{
+ cex.legend <- NULL
+ warning("cex.legend has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(cex.legend))
+ {
+ cex.legend <- 1
+ }
+
+ if(!is.null(cex.names))
+ {
+ if(is.numeric(cex.names))
+ {
+ cex.names <- cex.names
+ }else{
+ cex.names <- NULL
+ warning("cex.names has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(cex.names))
+ {
+ cex.names <- 1
+ }
+
+ if(!is.null(cex.axis))
+ {
+ if(is.numeric(cex.axis))
+ {
+ cex.axis <- cex.axis
+ }else{
+ cex.axis <- NULL
+ warning("cex.axis has to be of type 'numeric' \n used default settings")
+ }
+ }
+ if(is.null(cex.axis))
+ {
+ cex.axis <- 1
+ }
+
+ createPdf = FALSE
+ if(!is.null(file))
+ {
+ if (is.character(file)){
+ pdf(file = file, width = width, height = height)
+ createPdf = TRUE
+ }else{
+ file = NULL
+ warning("file has to be a of type 'character' \n could not create pdf")
+ }
+ }
+ if(is.null(file))
+ {
+ windows(width = width, height = height)
+ }
+
+ if(!is.null(col))
+ {
+ cl <- col
+ }
+ if(is.null(col))
+ {
+ alph <- 175
+ green <- rgb(127, 201, 127,alph, maxColorValue = 255)
+ lila <-rgb(190, 174, 212,alph, maxColorValue = 255)
+ orange <- rgb(253, 192, 134,alph, maxColorValue = 255)
+ yellow <- rgb(255, 255, 153,alph, maxColorValue = 255)
+ blue <- rgb(56, 108, 176,alph, maxColorValue = 255)
+ lightgray <- rgb(184,184,184,alph, maxColorValue = 255)
+ darkgray <- rgb(56,56,56,alph, maxColorValue = 255)
+ gray <- rgb(120,120,120,alph, maxColorValue = 255)
+ pink <- rgb(240,2,127,alph, maxColorValue = 255)
+ brown <- rgb(191,91,23,alph, maxColorValue = 255)
+ cl <- c(green,lila,orange,yellow,blue,lightgray,darkgray,gray,pink,brown)
+ while(length(cl)<length(x$effectNames)){
+ alph <- (alph+75)%%255
+ green <- rgb(127, 201, 127,alph, maxColorValue = 255)
+ lila <-rgb(190, 174, 212,alph, maxColorValue = 255)
+ orange <- rgb(253, 192, 134,alph, maxColorValue = 255)
+ yellow <- rgb(255, 255, 153,alph, maxColorValue = 255)
+ blue <- rgb(56, 108, 176,alph, maxColorValue = 255)
+ lightgray <- rgb(184,184,184,alph, maxColorValue = 255)
+ darkgray <- rgb(56,56,56,alph, maxColorValue = 255)
+ gray <- rgb(120,120,120,alph, maxColorValue = 255)
+ pink <- rgb(240,2,127,alph, maxColorValue = 255)
+ brown <- rgb(191,91,23,alph, maxColorValue = 255)
+ cl <- c(cl,green,lila,orange,yellow,blue,lightgray,darkgray,gray,pink,brown)
+ }
+ }
+ bordergrey <-"gray25"
+ values <- x$intervalValues
+ periods <- length(values)
+ effectNames <- x$effectNames
+ effectNumber <- length(effectNames)
+ lineTypes <- rep("solid",effectNumber)
+ legendNames <- effectNames
+ if(is.null(ylim))
+ {
+ ylim <-c(0,ceiling(max(unlist(lapply(values, max)))*10)*0.1)
+ }else{
+ ylim <-ylim
+ }
+ if(legend)
+ {
+ layout(rbind(1:periods, rep(periods+1,periods)),widths=rep(4, periods),heights=c(3,1))
+ }else{
+ layout(rbind(1:periods),widths=rep(4, periods),heights=c(3))
+ }
+ par( oma = c( 1, 3, 1, 3 ),mar = par()$mar+c(-5,-4.1,-4,-2.1), xpd=T )
+ for(period in 1:periods){
+ timeseries<-ts(t(values[[period]]))
+ plot.ts(timeseries, plot.type = "single", col = cl, lty = lineTypes, lwd = rep(1.5,effectNumber), bty = "n",xaxt = "n",yaxt = "n", ylab = ylab, xlab = xlab, ylim = ylim)
+ for(eff in 1:effectNumber)
+ {
+ points(ts(t(values[[period]]))[,eff], col = cl[eff], type = "p", pch = 20)
+ if(!is.null(staticValues))
+ {
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 267
More information about the Rsiena-commits
mailing list