[Rsiena-commits] r269 - in pkg: RSiena RSiena/R RSiena/inst/doc RSiena/man RSiena/tests RSienaTest RSienaTest/R RSienaTest/inst/doc RSienaTest/man RSienaTest/tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Apr 13 19:50:19 CEST 2014
Author: tomsnijders
Date: 2014-04-13 19:50:17 +0200 (Sun, 13 Apr 2014)
New Revision: 269
Added:
pkg/RSienaTest/R/sienaRI.r
pkg/RSienaTest/R/sienaRIDynamics.r
pkg/RSienaTest/man/sienaRI.Rd
Modified:
pkg/RSiena/DESCRIPTION
pkg/RSiena/NAMESPACE
pkg/RSiena/R/.Rhistory
pkg/RSiena/R/effectsDocumentation.r
pkg/RSiena/R/print01Report.r
pkg/RSiena/R/siena08.r
pkg/RSiena/R/sienaDataCreate.r
pkg/RSiena/R/sienaRI.r
pkg/RSiena/R/sienaRIDynamics.r
pkg/RSiena/R/sienaprint.r
pkg/RSiena/R/sienautils.r
pkg/RSiena/changeLog
pkg/RSiena/inst/doc/RSiena.bib
pkg/RSiena/inst/doc/RSiena_Manual.pdf
pkg/RSiena/inst/doc/RSiena_Manual.tex
pkg/RSiena/man/RSiena-package.Rd
pkg/RSiena/man/coDyadCovar.Rd
pkg/RSiena/man/sienaDependent.Rd
pkg/RSiena/man/sienaGroupCreate.Rd
pkg/RSiena/man/sienaRI.Rd
pkg/RSiena/man/varDyadCovar.Rd
pkg/RSiena/tests/parallel.Rout.save
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/NAMESPACE
pkg/RSienaTest/R/effectsDocumentation.r
pkg/RSienaTest/R/print01Report.r
pkg/RSienaTest/R/siena08.r
pkg/RSienaTest/R/sienaBayes.r
pkg/RSienaTest/R/sienaDataCreate.r
pkg/RSienaTest/R/sienaprint.r
pkg/RSienaTest/R/sienautils.r
pkg/RSienaTest/changeLog
pkg/RSienaTest/inst/doc/RSiena.bib
pkg/RSienaTest/inst/doc/RSiena_Manual.pdf
pkg/RSienaTest/inst/doc/RSiena_Manual.tex
pkg/RSienaTest/man/RSiena-package.Rd
pkg/RSienaTest/man/coDyadCovar.Rd
pkg/RSienaTest/man/sienaDependent.Rd
pkg/RSienaTest/man/sienaGroupCreate.Rd
pkg/RSienaTest/man/varDyadCovar.Rd
pkg/RSienaTest/tests/parallel.Rout.save
Log:
New version to make sienaRI.Rd() operational (I hope!). Also incorporation of some other changes: to sienaBayes, to have the option of centering dyadic covariates, and some other small things - see the changeLog.
Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION 2014-04-08 20:21:33 UTC (rev 268)
+++ pkg/RSiena/DESCRIPTION 2014-04-13 17:50:17 UTC (rev 269)
@@ -1,8 +1,8 @@
Package: RSiena
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.1-268
-Date: 2014-04-08
+Version: 1.1-269
+Date: 2014-04-13
Author: Ruth Ripley, Krists Boitmanis, Tom A.B. Snijders
Depends: R (>= 2.15.0)
Imports: Matrix
Modified: pkg/RSiena/NAMESPACE
===================================================================
--- pkg/RSiena/NAMESPACE 2014-04-08 20:21:33 UTC (rev 268)
+++ pkg/RSiena/NAMESPACE 2014-04-13 17:50:17 UTC (rev 269)
@@ -7,7 +7,7 @@
varCovar, varDyadCovar, setEffect, includeEffects, includeInteraction,
effectsDocumentation, sienaDataConstraint, print.xtable.sienaFit,
installGui, siena08, iwlsm, sienaTimeTest, includeTimeDummy,
- sienaGOF, descriptives.sienaGOF, sienaRI, sienaRIDynamics, sparseMatrixExtraction,
+ sienaGOF, descriptives.sienaGOF, sienaRI, sparseMatrixExtraction,
networkExtraction, behaviorExtraction,
OutdegreeDistribution, IndegreeDistribution, BehaviorDistribution,
siena.table, xtable,
@@ -52,7 +52,4 @@
S3method(summary, sienaRI)
S3method(print, sienaRI)
S3method(plot, sienaRI)
-S3method(summary, sienaRIDynamics)
-S3method(print, sienaRIDynamics)
-S3method(plot, sienaRIDynamics)
S3method(print, chains.data.frame)
Modified: pkg/RSiena/R/.Rhistory
===================================================================
--- pkg/RSiena/R/.Rhistory 2014-04-08 20:21:33 UTC (rev 268)
+++ pkg/RSiena/R/.Rhistory 2014-04-13 17:50:17 UTC (rev 269)
@@ -1 +1 @@
-print("Hi Natalie")
+print("Hi Natalie, best wishes from Tom")
Modified: pkg/RSiena/R/effectsDocumentation.r
===================================================================
--- pkg/RSiena/R/effectsDocumentation.r 2014-04-08 20:21:33 UTC (rev 268)
+++ pkg/RSiena/R/effectsDocumentation.r 2014-04-13 17:50:17 UTC (rev 269)
@@ -10,7 +10,8 @@
# *****************************************************************************/
##@effectsDocumentation Documentation
-effectsDocumentation <- function(effects= NULL, type="html", display=(type=="html"),
+effectsDocumentation <- function(effects= NULL, type="html",
+ display=(type=="html"),
filename=ifelse(is.null(effects), "effects", deparse(substitute(effects))))
{
if (is.null(effects))
@@ -21,8 +22,13 @@
}
else
{
- x <- as.data.frame(effects[, c("name", "effectName", "shortName", "type",
- "interaction1", "interaction2",
+ if (!inherits(effects,"sienaEffects"))
+ {
+ stop(paste(deparse(substitute(effects)),
+ "is not a sienaEffects object."))
+ }
+ x <- as.data.frame(effects[, c("name", "effectName", "shortName",
+ "type", "interaction1", "interaction2",
"parm", "interactionType")])
}
storage.mode(x$parm) <- "integer"
Modified: pkg/RSiena/R/print01Report.r
===================================================================
--- pkg/RSiena/R/print01Report.r 2014-04-08 20:21:33 UTC (rev 268)
+++ pkg/RSiena/R/print01Report.r 2014-04-13 17:50:17 UTC (rev 269)
@@ -713,26 +713,57 @@
(length(myvar) - nrow(myvar)), 1),
width=3, nsmall=1), '%)\n'), outf)
}
- Report("\nMeans of covariates:\n", outf)
- Report(c(format("minimum maximum mean", width=38,
+ Report("\nInformation about dyadic covariates:\n", outf)
+ Report(c(format("minimum maximum mean centered", width=48,
justify="right"), "\n"), outf)
+ any.cent <- 0
+ any.noncent <- 0
for (i in seq(along=covars))
{
atts <- attributes(x$dycCovars[[i]])
+ if (atts$centered)
+ {
+ cent <- " Y"
+ any.cent <- any.cent+1
+ }
+ else
+ {
+ cent <- " N"
+ any.noncent <- any.noncent+1
+ }
Report(c(format(covars[i], width=10),
format(round(atts$range2[1], 1),
nsmall=1, width=8),
format(round(atts$range2[2], 1),
nsmall=1, width=7),
format(round(atts$mean, 3),
- nsmall=3, width=10), "\n"), outf)
+ nsmall=3, width=10), cent, "\n"), outf)
}
Report('\n', outf)
- Report(c(ifelse(nCovars == 1,
- "This mean value is ", "These mean values are "),
- "subtracted from the dyadic covariate",
+
+ s.plural <- ifelse((any.cent >= 2),"s","")
+ if (any.noncent >= 1)
+ {
+ Report(c('The <mean> listed for the non-centered variable',
+ s.plural, ' is the attribute, not the observed mean.',
+ '\n'), sep="", outf)
+ }
+ if (any.noncent <= 0)
+ {
+ Report(c("The mean value", ifelse(nCovars == 1, " is", "s are"),
+ " subtracted from the",
+ ifelse(nCovars == 1, " centered", ""), " covariate",
ifelse(nCovars == 1, ".\n\n", "s.\n\n")), sep="", outf)
}
+ else if (any.cent >= 1)
+ {
+ Report(c("For the centered variable", s.plural,
+ ", the mean value", ifelse(any.cent == 1, " is", "s are"),
+ " subtracted from the covariate", s.plural,
+ ".\n"), sep="", outf)
+ }
+ }
+
##@reportChangingDyadicCovariates internal print01Report
reportChangingDyadicCovariates <- function()
{
@@ -775,21 +806,56 @@
width=3), '%)\n'), outf)
}
}
- Report("\nMeans of covariates:\n", outf)
+ Report("\nInformation about changing dyadic covariates:\n", outf)
+ Report(c(format("minimum maximum mean centered", width=48,
+ justify="right"), "\n"), outf)
+ any.cent <- 0
+ any.noncent <- 0
for (i in seq(along=covars))
{
atts <- attributes(x$dyvCovars[[i]])
+ if (atts$centered)
+ {
+ cent <- " Y"
+ any.cent <- any.cent+1
+ }
+ else
+ {
+ cent <- " N"
+ any.noncent <- any.noncent+1
+ }
Report(c(format(covars[i], width=10),
format(round(atts$mean, 3),
- nsmall=3, width=10), "\n"), outf)
+ nsmall=3, width=10), "\n"), cent, outf)
}
Report('\n', outf)
- Report(c(ifelse(nCovars == 1, "This ", "These "),
- "global mean value",
+ s.plural <- ifelse((any.cent >= 2),"s","")
+ if (any.noncent >= 1)
+ {
+ Report(c('The <mean> listed for the non-centered variable',
+ s.plural, ' is the attribute, not the observed mean.',
+ '\n'), sep="", outf)
+ }
+ if (nCovars >= 1)
+ {
+ if (any.noncent <= 0)
+ {
+ Report(c("The mean value",
ifelse(nCovars == 1, " is ", "s are "),
- "subtracted from the changing dyadic covariate",
+ " subtracted from the",
+ ifelse(nCovars == 1, " centered", ""), " covariate",
ifelse(nCovars ==1, ".\n\n", "s.\n\n")), sep="", outf)
}
+ else if (any.cent >= 1)
+ {
+ Report(c("For the centered variable", s.plural,
+ ", the mean value", ifelse(any.cent == 1, " is", "s are"),
+ " subtracted from the covariate", s.plural,
+ ".\n"), sep="", outf)
+ }
+ }
+ }
+
##@reportCompositionChange internal print01Report
reportCompositionChange <- function()
{
@@ -885,7 +951,7 @@
{
revision <- paste(" R-forge revision: ", rforgeRevision, " ", sep="")
}
- Report(c("SIENA version ", packageValues[[1]], " (",
+ Report(c("RSiena version ", packageValues[[1]], " (",
format(as.Date(packageValues[[2]]), "%d %m %Y"), ")",
revision, "\n\n"), sep="", outf)
Modified: pkg/RSiena/R/siena08.r
===================================================================
--- pkg/RSiena/R/siena08.r 2014-04-08 20:21:33 UTC (rev 268)
+++ pkg/RSiena/R/siena08.r 2014-04-13 17:50:17 UTC (rev 269)
@@ -208,6 +208,8 @@
## estimates
+ Report(c("\nTests for mean parameters use a t-distribution with N-1 d.f.,\n" ,
+ "where N = number of included groups.\n"), sep="", outf)
Report(c("\nUpper bound used for standard error is",
format(round(x$bound, 4), width=9, nsmall=2), ".\n"), sep="", outf)
@@ -472,7 +474,7 @@
revision <- paste(" R-forge revision: ", rforgeRevision,
" ", sep="")
}
- Report(c("SIENA version ", packageValues[[1]], " (",
+ Report(c("RSiena version ", packageValues[[1]], " (",
format(as.Date(packageValues[[2]]), "%d %m %Y"), ")",
revision, "\n\n"), sep="", outf)
}
Modified: pkg/RSiena/R/sienaDataCreate.r
===================================================================
--- pkg/RSiena/R/sienaDataCreate.r 2014-04-08 20:21:33 UTC (rev 268)
+++ pkg/RSiena/R/sienaDataCreate.r 2014-04-13 17:50:17 UTC (rev 269)
@@ -115,7 +115,7 @@
rr <- range(x, na.rm=TRUE)
nonMissingCount <- sum(!is.na(x))
}
- attr(x,'mean') <- varmean
+ attr(x,'mean') <- ifelse(attr(x,'centered'), varmean, 0)
attr(x,'range') <- rr[2] - rr[1]
storage.mode(attr(x, 'range')) <- 'double'
attr(x,'range2') <- rr
@@ -178,7 +178,7 @@
attr(x, "meanp") <- colMeans(x, dims=2, na.rm=TRUE)
nonMissingCounts <- colSums(!is.na(x), dims=2)
}
- attr(x, "mean") <- varmean
+ attr(x, "mean") <- ifelse(attr(x, "centered"), varmean, 0)
attr(x, "range") <- rr[2] - rr[1]
storage.mode(attr(x, "range")) <- "double"
attr(x, "name") <- name
@@ -1161,7 +1161,7 @@
j <- match(atts$cCovars[covar], names(group[[i]]$cCovars))
if (is.na(j))
{
- stop("inconsistent covariate names")
+ stop("inconsistent actor covariate names")
}
thisrange[, i] <- range(group[[i]]$cCovars[[j]],
na.rm=TRUE)
@@ -1206,7 +1206,7 @@
j <- match(atts$vCovars[covar], names(group[[i]]$vCovars))
if (is.na(j))
{
- stop("inconsistent covariate names")
+ stop("inconsistent actor covariate names")
}
vartotal <- vartotal + attr(group[[i]]$vCovars[[j]], "vartotal")
nonMissingCount <- nonMissingCount +
@@ -1225,7 +1225,7 @@
j <- match(atts$vCovars[covar], names(group[[i]]$vCovars))
if (is.na(j))
{
- stop("inconsistent covariate names")
+ stop("inconsistent actor covariate names")
}
group[[i]]$vCovars[[j]] <- group[[i]]$vCovars[[j]] -
@@ -1244,7 +1244,7 @@
j <- match(atts$vCovars[covar], names(group[[i]]$vCovars))
if (is.na(j))
{
- stop("inconsistent covariate names")
+ stop("inconsistent actor covariate names")
}
thisrange[, i] <- range(group[[i]]$vCovars[[j]],
na.rm=TRUE)
@@ -1287,7 +1287,7 @@
j <- match(atts$dycCovars[covar], names(group[[1]]$dycCovars))
if (is.na(j))
{
- stop("inconsistent covariate names")
+ stop("inconsistent dyadic covariate names")
}
dycCovarMean[covar] <- attr(group[[1]]$dycCovars[[j]], "mean")
dycCovarRange[covar] <- attr(group[[1]]$dycCovars[[j]], "range")
@@ -1305,7 +1305,7 @@
j <- match(atts$dyvCovars[covar], names(group[[i]]$dyvCovars))
if (is.na(j))
{
- stop("inconsistent covariate names")
+ stop("inconsistent dyadic covariate names")
}
sparse <- attr(group[[i]]$dyvCovars[[j]], "sparse")
vardims <- attr(group[[i]]$dyvCovars[[j]], "vardims")
@@ -1330,8 +1330,14 @@
thisrange[, i] <- range(group[[i]]$dyvCovars[[j]],
na.rm=TRUE)
}
+ centered <- attr(group[[i]]$dyvCovars[[j]], "centered")
+ if (i <= 1) {centered1 <- centered}
+ if (centered != centered1)
+ {
+ stop("inconsistent centering of dyadic covariates")
+ }
}
- dyvCovarMean[covar] <- vartotal / nonMissingCount
+ dyvCovarMean[covar] <- ifelse(centered, vartotal / nonMissingCount, 0)
rr <- range(thisrange, na.rm=TRUE)
dyvCovarRange[covar] <- rr[2] - rr[1]
}
@@ -1569,7 +1575,7 @@
covarsub <- match(varname, cCovars)
if (is.na(covarsub))
{
- stop('covariate names inconsistent')
+ stop('actor covariate names inconsistent')
}
attribs <- attributes(objlist[[i]]$cCovars[[j]])
if (is.na(ccnodeSets[covarsub]))
@@ -1606,7 +1612,7 @@
covarsub <- match(varname, dycCovars)
if (is.na(covarsub))
{
- stop('covariate names inconsistent')
+ stop('dyadic covariate names inconsistent')
}
attribs <- attributes(objlist[[i]]$dycCovars[[j]])
if (is.null(dycnodeSets[[covarsub]]))
@@ -1632,7 +1638,7 @@
covarsub <- match(varname, dyvCovars)
if (is.na(covarsub))
{
- stop('covariate names inconsistent')
+ stop('dyadic covariate names inconsistent')
}
attribs <- attributes(objlist[[i]]$dyvCovars[[j]])
if (is.null(dyvnodeSets[[covarsub]]))
@@ -1698,6 +1704,7 @@
attr(newcovar, "nonMissingCount") <-
attr(const[[j]], "nonMissingCount")
attr(newcovar, "mean") <- attr(const[[j]], "mean")
+ attr(newcovar, "centered") <- attr(const[[j]], "centered")
attr(newcovar, "range") <- attr(const[[j]], "range")
attr(newcovar, "rangep") <- rep(attr(const[[j]], "range"),
dim3)
Modified: pkg/RSiena/R/sienaRI.r
===================================================================
--- pkg/RSiena/R/sienaRI.r 2014-04-08 20:21:33 UTC (rev 268)
+++ pkg/RSiena/R/sienaRI.r 2014-04-13 17:50:17 UTC (rev 269)
@@ -1,551 +1,552 @@
-#/******************************************************************************
-# * 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)
-}
-
+#/******************************************************************************
+# * 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{
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 269
More information about the Rsiena-commits
mailing list