[Rsiena-commits] r276 - in pkg: RSiena RSiena/R RSiena/inst/doc RSiena/man RSienaTest RSienaTest/R RSienaTest/doc RSienaTest/inst/doc RSienaTest/man RSienaTest/tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jun 18 17:07:02 CEST 2014
Author: tomsnijders
Date: 2014-06-18 17:07:01 +0200 (Wed, 18 Jun 2014)
New Revision: 276
Added:
pkg/RSienaTest/R/bayesTest.r
pkg/RSienaTest/doc/HowToCommit.pdf
pkg/RSienaTest/doc/HowToCommit.tex
pkg/RSienaTest/man/bayesTest.Rd
Modified:
pkg/RSiena/ChangeLog
pkg/RSiena/DESCRIPTION
pkg/RSiena/R/sienaRI.r
pkg/RSiena/R/sienaRIDynamics.r
pkg/RSiena/R/sienaprint.r
pkg/RSiena/inst/doc/RSiena_Manual.pdf
pkg/RSiena/inst/doc/RSiena_Manual.tex
pkg/RSiena/man/RSiena-package.Rd
pkg/RSiena/man/siena07.Rd
pkg/RSiena/man/sienaGroupCreate.Rd
pkg/RSienaTest/ChangeLog
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/NAMESPACE
pkg/RSienaTest/R/sienaBayes.r
pkg/RSienaTest/R/sienaRI.r
pkg/RSienaTest/R/sienaRIDynamics.r
pkg/RSienaTest/R/sienaprint.r
pkg/RSienaTest/inst/doc/RSiena_Manual.pdf
pkg/RSienaTest/inst/doc/RSiena_Manual.tex
pkg/RSienaTest/man/RSiena-package.Rd
pkg/RSienaTest/man/print.sienaBayesFit.Rd
pkg/RSienaTest/man/siena07.Rd
pkg/RSienaTest/man/sienaBayes.Rd
pkg/RSienaTest/tests/parallel.R
Log:
Modified: pkg/RSiena/ChangeLog
===================================================================
--- pkg/RSiena/ChangeLog 2014-06-10 10:07:55 UTC (rev 275)
+++ pkg/RSiena/ChangeLog 2014-06-18 15:07:01 UTC (rev 276)
@@ -1,20 +1,31 @@
-2014-06-02 R-Forge Revision 276
+2014-06-17 R-Forge Revision 276
Changes in RSiena and RSienaTest:
- * Higher writefreq for batch operation (phase2.r, phase3.r)
+ * Some documentation of returnChains in siena07.Rd.
+Changes in RSienaTest:
+ * Added function simpleBayesTest() (bayesTest.r, bayesTest.Rd).
+ * Added function multipleBayesTest() with print and plot
+ methods (bayesTest.r, bayesTest.Rd).
+ * Small changes in various functions in sienaBayes.r.
+ * Corrected starting printing sienaBayesFit at nfirst (sienaprint.r).
+
+2014-06-02 R-Forge Revision 275, but version predated 1.1-276
+Changes in RSiena and RSienaTest:
+ * Higher writefreq for batch operation (phase2.r, phase3.r)
* Error check for name of dependent variable to condition upon (initializeFRAN.r)
- * Some corrections to account for behavior dependent variables that are not
+ * Some corrections to allow behavior dependent variables that are not
coded as integer (in initializeFRAN and unpackBehavior, initializeFRAN.r)
* Parameters fix and test now used in includeEffects() (sienaeffects.r)
* Small bug fix in CalculateDerivative3 (phase3.r)
+ * Avoidance of ":::" in siena01.r and document.r by small changes.
Changes in RSienaTest:
- * Changes in sienaBayes.r and its print and summary methods;
+ * Changes in sienaBayes.r and its print and summary methods;
new function glueBayes.
2014-05-05 R-Forge Revision 275
Changes in RSiena and RSienaTest:
* Error message and stop in case of mismatch variable names
in algorithm object and in data (initializeFRAN.r).
- * Force behavior and its attribute "change" to be integer in function
+ * Force behavior and its attribute "change" to be integer in function
unpackBehavior() (initializeFRAN.r).
* Some tabs in .Rd files change to spaces.
@@ -22,7 +33,7 @@
Changes in RSiena and RSienaTest:
* Replaced (!strcmp...) by (!(strcmp...)) one more time
in sienainternals07.cpp.
-
+
2014-04-25 R-Forge Revision 273
Changes in RSiena and RSienaTest:
* NAMESPACE: removed entries sienaRIDynamics with the corresponding
@@ -36,7 +47,7 @@
* Small change in print01Report.r to improve reporting
two files of composition change.
* sienaRI.r: require that <<file>> argument is not NULL for non-Windows
- operating systems.
+ operating systems.
Changes in RSienaTest:
* Minor changes in sienaBayes.r.
Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION 2014-06-10 10:07:55 UTC (rev 275)
+++ pkg/RSiena/DESCRIPTION 2014-06-18 15:07:01 UTC (rev 276)
@@ -2,7 +2,7 @@
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
Version: 1.1-276
-Date: 2014-06-02
+Date: 2014-06-18
Author: Ruth Ripley, Krists Boitmanis, Tom A.B. Snijders
Depends: R (>= 2.15.0)
Imports: Matrix
Modified: pkg/RSiena/R/sienaRI.r
===================================================================
--- pkg/RSiena/R/sienaRI.r 2014-06-10 10:07:55 UTC (rev 275)
+++ pkg/RSiena/R/sienaRI.r 2014-06-18 15:07:01 UTC (rev 276)
@@ -26,9 +26,11 @@
{
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){
+ 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{
@@ -128,7 +130,9 @@
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 <- .Call("getTargets", PACKAGE=pkgname, pData, pModel, myeffects,
+ parallelrun=TRUE, returnActorStatistics=FALSE,
+ returnStaticChangeContributions=TRUE)
ans
}
@@ -169,12 +173,14 @@
currentDepEffs <- effects$name == currentDepName
effNumber <- sum(currentDepEffs)
- RIs <- data.frame(row.names = effectIds[currentDepEffs])
- RIs <- cbind(RIs, matrix(0, nrow=effNumber, ncol = actors))
+# 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=".")
+# 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()
@@ -183,10 +189,16 @@
for(w in 1:waves)
{
currentDepEffectContributions <- conts[[1]][[w]][currentDepEffs]
- currentDepEffectContributions <- sapply(lapply(currentDepEffectContributions, unlist), matrix, nrow=actors, ncol=choices, byrow=TRUE, simplify="array")
+ 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)})
+ 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,
@@ -194,11 +206,11 @@
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
+# RIs <- RIs_matrix
entropies <- entropy_vector
- RIActors[[w]] <- apply(RIs, 2, function(x){x/sum(x)})
- absoluteSumActors[[w]] <- colSums(RIs)
+ RIActors[[w]] <- apply(RIs_matrix, 2, function(x){x/sum(x)})
+ absoluteSumActors[[w]] <- colSums(RIs_matrix)
entropyActors[[w]] <- entropies
expectedRI[[w]] <- rowSums(RIActors[[w]] )/dim(RIActors[[w]])[2]
}
@@ -212,7 +224,9 @@
{
RItmp$effectNames <- effectNames[currentDepEffs]
}else{
- RItmp$effectNames <- paste(effectTypes[currentDepEffs], " ", effects$effectName[currentDepEffs], sep="")
+ RItmp$effectNames <-
+ paste(effectTypes[currentDepEffs], " ",
+ effects$effectName[currentDepEffs], sep="")
}
class(RItmp) <- "sienaRI"
if(depNumber == 1){
@@ -240,12 +254,14 @@
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)))
+ 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[eff+1,] <-
+ exp(colSums(t*effectContributions))/sum(exp(colSums(t*effectContributions)))
}
distributions
}
@@ -332,7 +348,10 @@
##@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, ...)
+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"))
{
@@ -405,7 +424,7 @@
{
if(addPieChart)
{
- width = (actors/3+2)*1
+ width = (actors/3+2)*0.9
}else{
width = (actors/3+1)*1
}
@@ -522,7 +541,9 @@
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 )
+ 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)
{
@@ -532,11 +553,16 @@
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 )
+ 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."))
+ 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)
Modified: pkg/RSiena/R/sienaRIDynamics.r
===================================================================
--- pkg/RSiena/R/sienaRIDynamics.r 2014-06-10 10:07:55 UTC (rev 275)
+++ pkg/RSiena/R/sienaRIDynamics.r 2014-06-18 15:07:01 UTC (rev 276)
@@ -6,7 +6,7 @@
# * File: sienaRIDynamics.r
# *
# * Description: Used to determine, print, and plots relative importances of effects
-# * in sequences of simulated micro-steps
+# * in sequences of simulated micro-steps
# *****************************************************************************/
##@sienaRIDynamics. Use as RSiena:::sienaRIDynamics()
@@ -51,7 +51,7 @@
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')")
}
effs <- ans$effects
- }
+ }
if(!is.null(intervalsPerPeriod))
{
if(is.numeric(intervalsPerPeriod))
@@ -66,7 +66,7 @@
{
intervalsPerPeriod <- 10
}
- RIValues <- calculateRIDynamics(data = data, theta= c(ans$rate,ans$theta), algorithm = ans$x, effects = effs, depvar = currentNetName, intervalsPerPeriod=intervalsPerPeriod)
+ RIValues <- calculateRIDynamics(data = data, theta= c(ans$rate,ans$theta), algorithm = ans$x, effects = effs, depvar = currentNetName, intervalsPerPeriod=intervalsPerPeriod)
}else{
if (!inherits(algorithm, "sienaAlgorithm"))
{
@@ -76,7 +76,7 @@
if (!inherits(effects, "sienaEffects"))
{
stop("effects is not a legitimate Siena effects object")
- }
+ }
effs <- effects
if(!is.numeric(theta))
{
@@ -86,9 +86,9 @@
{
if(length(theta) == sum(effs$include==TRUE & effs$type!="rate"))
{
- stop("vector of model parameters has wrong dimension, maybe rate parameters are missing")
+ 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")
+ stop("theta is not a legitimate parameter vector \n number of parameters has to match number of effects")
}
paras <- theta
if(!is.null(intervalsPerPeriod))
@@ -122,7 +122,7 @@
z$print <- FALSE
z$Phase <- 3
z <- initializeFRAN(z, x, data, effects, prevAns=NULL, initC=FALSE, returnDeps=FALSE)
- z$returnChangeContributions <- TRUE
+ z$returnChangeContributions <- TRUE
z$theta <- theta
if (!is.null(x$randomSeed))
{
@@ -137,38 +137,41 @@
}
}
chains <- x$n3
- periods <- data$observation-1
+ 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]
+# networkName <- effects$name[noRate]
networkInteraction <- effects$interaction1[noRate]
- effectIds <- paste(effectNames,effectTypes,networkInteraction, sep = ".")
+ effectIds <- paste(effectNames,effectTypes,networkInteraction, sep = ".")
currentNetObjEffs <- effects$name[noRate] == currentNetName
-# The old code leading to a crash:
-# RIintervalValues <- list()
-# for(period in 1:periods)
-# {
-# RIintervalValues[[period]] <- data.frame(row.names = effectIds[currentNetObjEffs])
-# }
+# The old code leading to a crash...
+ RIintervalValues <- list()
+# ...threfore, the following line is added
+ blanks <- matrix(0, sum(currentNetObjEffs),intervalsPerPeriod)
+ for(period in 1:periods)
+ {
+ RIintervalValues[[period]] <- data.frame(blanks, row.names = effectIds[currentNetObjEffs])
+ }
for (chain in (1:chains))
{
-cat("The following line leads to an error\n")
-browser()
- ans <- z$FRAN(z, x)
+#cat("The following line leads to an error\n")
+#browser()
+ ans <- z$FRAN(z, x)
for(period in 1:periods)
{
microSteps <- length(ans$changeContributions[[1]][[period]])
stepsPerInterval <- microSteps/intervalsPerPeriod
-# new code TS (7 lines):
- blanks <- matrix(NA, sum(currentNetObjEffs),
- ceiling(microSteps/stepsPerInterval))
- RIintervalValues <-
- as.list(rep(data.frame(blanks,
- row.names = effectIds[currentNetObjEffs]), periods))
- blanks <- matrix(NA, sum(currentNetObjEffs), microSteps)
+# Tom's new code TS (7 lines).... Natalie slightly modified the next 4 lines
+# and moved them one loop higher (out of the "chains"-loop)
+# blanks <- matrix(NA, sum(currentNetObjEffs),
+# ceiling(microSteps/stepsPerInterval))
+# RIintervalValues <-
+# as.list(rep(data.frame(blanks,
+# row.names = effectIds[currentNetObjEffs]), periods))
+ blanks <- matrix(0, sum(currentNetObjEffs), intervalsPerPeriod)
RItmp <- data.frame(blanks, row.names = effectIds[currentNetObjEffs])
interval <- 1
stepsInIntervalCounter <- 0
@@ -178,53 +181,53 @@
if(attr(ans$changeContributions[[1]][[period]][[microStep]],
"networkName")==currentNetName)
{
- stepsInIntervalCounter <- stepsInIntervalCounter + 1
- ## distributions[1,] contains
+ stepsInIntervalCounter <- stepsInIntervalCounter + 1
+ ## distributions[1,] contains
## the probabilities of the available choices in this micro-step
- ## distributions[2,],distributions[3,], ...contains
+ ## distributions[2,],distributions[3,],distributions[4,], ...contains
## the probabilities of the available choices
- ## if the parameter of the second, third, ... effects is set to zero.
+ ## if the parameter of the first, second, third ... effects is set to zero.
distributions <- calculateDistributions(
- ans$changeContributions[[1]][[period]][[microStep]],
+ ans$changeContributions[[1]][[period]][[microStep]],
thetaNoRate[currentNetObjEffs])
- ## If one wishes another measure
+ ## If one wishes another measure
## than the L^1-difference between distributions,
- ## here is the right place
+ ## here is the right place
## to call some new function instead of "L1D".
- RItmp[,stepsInIntervalCounter] <-
- standardize(L1D(distributions[1,],
- distributions[2:dim(distributions)[1],]))
+ RItmp[,stepsInIntervalCounter] <-
+ standardize(L1D(distributions[1,],
+ distributions[2:dim(distributions)[1],]))
}
- if (microStep >
- interval * stepsPerInterval || microStep == microSteps)
+ if (microStep >
+ interval * stepsPerInterval || microStep == microSteps)
{
if(chain == 1)
{
- RIintervalValues[[period]][,interval] <-
+ RIintervalValues[[period]][,interval] <-
rowSums(RItmp)/length(RItmp)
}
else
{
- RIintervalValues[[period]][,interval] <-
- RIintervalValues[[period]][,interval] +
+ 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], " ",
+ RIDynamics$effectNames <- paste(effectTypes[currentNetObjEffs], " ",
(effects$effectName[noRate])[currentNetObjEffs], sep="")
- class(RIDynamics) <- "sienaRIDynamics"
+ class(RIDynamics) <- "sienaRIDynamics"
RIDynamics
}
@@ -232,7 +235,7 @@
{
newValues <- values/sum(values)
newValues
-}
+}
##@print.sienaRIDynamics Methods
@@ -261,8 +264,8 @@
if(col == 5)
{
line2 <- paste(line2, rep('\n',effs), sep="")
- cat(as.matrix(line1),'\n \n', sep='')
- cat(as.matrix(line2),'\n', 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
@@ -271,8 +274,8 @@
if(col>0)
{
line2 <- paste(line2, rep('\n',effs), sep="")
- cat(as.matrix(line1),'\n \n', sep='')
- cat(as.matrix(line2),'\n', sep='')
+ cat(as.matrix(line1),'\n \n', sep='')
+ cat(as.matrix(line2),'\n', sep='')
}
}
invisible(x)
@@ -301,8 +304,10 @@
##@plot.sienaRIDynamics Methods
-plot.sienaRIDynamics <- function(x, staticRI = NULL, file = NULL, col = NULL, ylim=NULL, width = NULL, height = NULL, legend = TRUE, legendColumns = NULL, legendHeight = NULL, cex.legend = NULL, ...)
-{
+plot.sienaRIDynamics <- function(x, staticRI = NULL, file = NULL, col = NULL,
+ ylim=NULL, width = NULL, height = NULL, legend = TRUE,
+ legendColumns = NULL, legendHeight = NULL, cex.legend = NULL, ...)
+{
if(!inherits(x, "sienaRIDynamics"))
{
stop("x is not of class 'sienaRIDynamics' ")
@@ -313,7 +318,9 @@
{
warning("staticRI is not of class 'sienaRI' and is therefore ignored")
staticValues <- NULL
- }else{
+ }
+ else
+ {
if(staticRI$dependentVariable != x$dependentVariable)
{
warning("staticRI does not correspond to x and is therefore ignored,\n staticRI and x do not refer to the same dependent variable")
@@ -368,7 +375,7 @@
if(is.null(legendHeight))
{
legendHeight <- 1
- }
+ }
}
if(!is.null(height))
{
@@ -389,7 +396,7 @@
height <- 3
}
}
-
+
if(!is.null(width))
{
if(is.numeric(width))
@@ -404,7 +411,7 @@
{
width <- 8
}
-
+
if(!is.null(cex.legend))
{
if(is.numeric(cex.legend))
@@ -419,7 +426,7 @@
{
cex.legend <- 1
}
-
+
createPdf = FALSE
if(!is.null(file))
{
@@ -435,7 +442,7 @@
{
windows(width = width, height = height)
}
-
+
if(!is.null(col))
{
cl <- col
@@ -467,10 +474,10 @@
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
+# bordergrey <-"gray25"
+ values <- x$intervalValues
periods <- length(values)
effectNames <- x$effectNames
effectNumber <- length(effectNames)
@@ -484,20 +491,24 @@
}
if(legend)
{
- layout(rbind(1:periods, rep(periods+1,periods)),widths=rep(4, periods),heights=c(3,legendHeight))
+ layout(rbind(1:periods, rep(periods+1,periods)),
+ widths=rep(4, periods),heights=c(3,legendHeight))
}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 )
+ 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 ="", xlab = "", ylim = ylim)
+ plot.ts(timeseries, plot.type = "single", col = cl,
+ lty = lineTypes, lwd = rep(1.5,effectNumber), bty = "n",
+ xaxt = "n",yaxt = "n", ylab ="", 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))
{
- points(xy.coords(1,staticValues[[period]][eff]),col = cl[eff], type = "p", pch = 1, cex = 1.75)
+ points(xy.coords(1,staticValues[[period]][eff]),col = cl[eff],
+ type = "p", pch = 1, cex = 1.75)
}
}
ax <- ((ylim[1]*10):(ylim[2]*10))/10
@@ -519,8 +530,9 @@
}
if(legend)
{
- plot(c(0,1), c(0,1), col=rgb(0,0,0,0),axes=FALSE, ylab = "", xlab = "")
- legend(0, 1, legendNames, col = cl[1:effectNumber], lwd = 2, lty = lineTypes, ncol = legendColumns, bty = "n",cex=cex.legend)
+ plot(c(0,1), c(0,1), col=rgb(0,0,0,0),axes=FALSE, ylab = "", xlab = "")
+ legend(0, 1, legendNames, col = cl[1:effectNumber], lwd = 2,
+ lty = lineTypes, ncol = legendColumns, bty = "n",cex=cex.legend)
}
if(createPdf)
{
Modified: pkg/RSiena/R/sienaprint.r
===================================================================
--- pkg/RSiena/R/sienaprint.r 2014-06-10 10:07:55 UTC (rev 275)
+++ pkg/RSiena/R/sienaprint.r 2014-06-18 15:07:01 UTC (rev 276)
@@ -518,98 +518,114 @@
}
##@averageTheta.last Miscellaneous
-averageTheta.last <- function(z, groupOnly=0)
+averageTheta.last <- function(z, groupOnly=0, nfirst=z$nwarm+1)
{
ntot <- sum(!is.na(z$ThinPosteriorMu[,1]))
+ if (nfirst > ntot)
+ {
+ stop('Sample did not come beyond warming')
+ }
thetaMean <- rep(NA, z$pp)
for (group in 1:z$nGroup)
{
thetaMean[z$ratePositions[[group]]] <- colMeans(
- z$ThinParameters[(z$nwarm+1):ntot,
+ z$ThinParameters[nfirst:ntot,
group, !z$generalParametersInGroup, drop=FALSE], na.rm=TRUE)
}
if (groupOnly != 0)
{
thetaMean[(z$set1)&(!z$basicRate)] <- colMeans(
- z$ThinParameters[(z$nwarm+1):ntot, groupOnly,
+ z$ThinParameters[(nfirst):ntot, groupOnly,
z$varyingGeneralParametersInGroup, drop=FALSE], na.rm=TRUE)
}
else
{
thetaMean[(z$set1)&(!z$basicRate)] <- colMeans(
- z$ThinPosteriorMu[(z$nwarm+1):ntot,
+ z$ThinPosteriorMu[(nfirst):ntot,
z$objectiveInVarying, drop=FALSE], na.rm=TRUE)
}
if (any(z$set2))
{
thetaMean[z$set2] <-
- colMeans(z$ThinPosteriorEta[(z$nwarm+1):ntot,, drop=FALSE]
+ colMeans(z$ThinPosteriorEta[nfirst:ntot,, drop=FALSE]
, na.rm=TRUE)
}
thetaMean
}
##@sdTheta.last Miscellaneous
-sdTheta.last <- function(z, groupOnly=0)
+sdTheta.last <- function(z, groupOnly=0, nfirst=z$nwarm+1)
{
ntot <- sum(!is.na(z$ThinPosteriorMu[,1]))
+ if (nfirst >= ntot-1)
+ {
+ stop('Sample did not come beyond warming')
+ }
sdTheta<- rep(NA, z$pp)
for (group in 1:z$nGroup)
{
sdTheta[z$ratePositions[[group]]] <- apply(
- z$ThinParameters[(z$nwarm+1):ntot,
+ z$ThinParameters[nfirst:ntot,
group, !z$generalParametersInGroup, drop=FALSE], 3, sd)
}
if (groupOnly != 0)
{
sdTheta[(z$set1)&(!z$basicRate)] <- apply(
- z$ThinParameters[(z$nwarm+1):ntot, groupOnly,
+ z$ThinParameters[nfirst:ntot, groupOnly,
z$varyingGeneralParametersInGroup, drop=FALSE], 3, sd)
}
else
{
sdTheta[(z$set1)&(!z$basicRate)] <- apply(
- z$ThinPosteriorMu[(z$nwarm+1):ntot,
+ z$ThinPosteriorMu[nfirst:ntot,
z$objectiveInVarying, drop=FALSE], 2, sd)
}
if (any(z$set2))
{
sdTheta[z$set2] <-
- apply(z$ThinPosteriorEta[(z$nwarm+1):ntot,, drop=FALSE],
+ apply(z$ThinPosteriorEta[nfirst:ntot,, drop=FALSE],
2, sd)
}
sdTheta
}
##@credValues Miscellaneous
-credValues <- function(z, groupOnly=0)
+credValues <- function(z, theProbs = c(0.025,0.975), tested = 0,
+ groupOnly=0, nfirst=z$nwarm+1)
{
ntot <- sum(!is.na(z$ThinPosteriorMu[,1]))
+ if (nfirst >= ntot-1)
+ {
+ stop('Sample did not come beyond warming')
+ }
credVals <- matrix(NA, length(z$set1), 3)
# But for the rate parameters NAs will remain.
- cvp <- function(x){c(quantile(x, probs = c(0.025,0.975)), 1-(ecdf(x))(0))}
+ cvp <- function(x, test0){c(quantile(x, probs=theProbs), 1-(ecdf(x))(test0))}
if (groupOnly != 0)
{
credVals[(z$set1)&(!z$basicRate), ] <- t(apply(
- z$ThinParameters[(z$nwarm+1):ntot, groupOnly,
- z$varyingGeneralParametersInGroup, drop=FALSE], 3, cvp))
+ z$ThinParameters[nfirst:ntot, groupOnly,
+ z$varyingGeneralParametersInGroup, drop=FALSE], 3,
+ cvp, test0 = tested))
}
else
{
- credVals[(z$set1)&(!z$basicRate), ] <- t(apply(
- z$ThinPosteriorMu[(z$nwarm+1):ntot,
- z$objectiveInVarying, drop=FALSE], 2, cvp))
+ credVals[(z$set1)&(!z$basicRate), ] <-
+ t(apply(
+ z$ThinPosteriorMu[nfirst:ntot,
+ z$objectiveInVarying, drop=FALSE], 2, cvp, test0 = tested))
}
if (any(z$set2))
{
credVals[z$set2, ] <-
- t(apply(z$ThinPosteriorEta[(z$nwarm+1):ntot,, drop=FALSE], 2, cvp))
+ t(apply(z$ThinPosteriorEta[nfirst:ntot,, drop=FALSE], 2,
+ cvp, test0 = tested))
}
credVals
}
##@sienaFitThetaTable Miscellaneous
-sienaFitThetaTable <- function(x, fromBayes=FALSE, tstat=FALSE, groupOnly=0)
+sienaFitThetaTable <- function(x, fromBayes=FALSE, tstat=FALSE, groupOnly=0, nfirst)
{
theEffects <- x$effects
# abcd here it was not OK when interactions had creation/endowment effects
@@ -724,7 +740,7 @@
}
if (fromBayes)
{
- theta <- averageTheta.last(x, groupOnly)
+ theta <- averageTheta.last(x, groupOnly, nfirst = nfirst)
}
else
{
@@ -758,9 +774,9 @@
if (fromBayes)
{
- mydf[nrates + (1:xp), 'se' ] <- sdTheta.last(x, groupOnly)
+ mydf[nrates + (1:xp), 'se' ] <- sdTheta.last(x, groupOnly, nfirst=nfirst)
mydf[nrates + (1:xp), 'random' ] <- NA
- credVal <- credValues(x, groupOnly)
+ credVal <- credValues(x, groupOnly = groupOnly, nfirst = nfirst)
mydf[nrates + (1:xp), 'cFrom' ] <- credVal[,1]
mydf[nrates + (1:xp), 'cTo' ] <- credVal[,2]
mydf[nrates + (1:xp), 'p' ] <- credVal[,3]
@@ -892,7 +908,7 @@
}
##@maketemp Methods
-makeTemp <- function(x, groupOnly=0, ...)
+makeTemp <- function(x, groupOnly=0, nfirst, ...)
{
# This reproduces a part of sienaFit but with Bayesian headings;
# for use in print.sienaBayesFit and summary.sienaBayesFit
@@ -901,16 +917,19 @@
x$requestedEffects[x$requestedEffects$groupName==name1,]
x$pp <- length(x$theta)
x$fixed <- x$fixed[x$effects$groupName==name1]
- tmp <- sienaFitThetaTable(x, fromBayes=TRUE, groupOnly=groupOnly)
+ if (is.null(nfirst))
+ {
+ nfirst <- x$nwarm + 1
+ }
+ tmp <- sienaFitThetaTable(x, fromBayes=TRUE, groupOnly=groupOnly, nfirst=nfirst)
mydf <- tmp$mydf
mymat <- as.matrix(mydf[,names(mydf)!="tstat"])
- mynames <- colnames(mymat)
+# mynames <- colnames(mymat)
# mymat <- cbind(mymat, rep.int(NA, dim(mymat)[1]))
# mymat <- cbind(mymat, matrix(NA, dim(mymat)[1], 4))
# colnames(mymat) <- c(mynames, 'random', 'credFrom','credTo', 'p')
mymat[, 'value'] <- format(round(mydf$value, digits=4))
mymat[, 'se'] <- format(round(mydf$se, digits=4))
-# x$basicRate
mymat[x$set1, 'random'] <- " + "
mymat[x$set2, 'random'] <- " - "
mymat[x$basicRate, 'random'] <- " "
@@ -1022,7 +1041,7 @@
stop("This object did not come beyond the warming phase.\n")
}
}
- tmps <- makeTemp(x)
+ tmps <- makeTemp(x, nfirst=nfirst)
tmp <- tmps[[1]]
tmp1 <- tmps[[2]]
addtorow <- tmp$addtorow
@@ -1141,7 +1160,11 @@
cat("\nBasic rates parameters are treated as incidental parameters.\n\n")
}
cat("\nAlgorithm specifications were nwarm =",x$nwarm,", nmain =", x$nmain,
- ", nrunMHBatches =", x$nrunMHBatches, ".\n")
+ ", nrunMHBatches =", x$nrunMHBatches, ", mult =", x$mult, ".\n")
+ if (!is.null(nfirst))
+ {
+ cat("For the results, nwarm is superseded by nfirst = ", nfirst, ".")
+ }
if (ntot < x$nwarm + x$nmain)
{
cat("This object resulted from an intermediate save, after",
@@ -1176,7 +1199,7 @@
for (h in 1:length(x$f$groupNames))
{
cat("\n", x$f$groupNames[h], "\n")
- tmps <- makeTemp(x, groupOnly=h)
+ tmps <- makeTemp(x, groupOnly=h, nfirst=nfirst)
tmp <- tmps[[1]]
tmp1 <- tmps[[2]]
addtorow <- tmp$addtorow
Modified: pkg/RSiena/inst/doc/RSiena_Manual.pdf
===================================================================
--- pkg/RSiena/inst/doc/RSiena_Manual.pdf 2014-06-10 10:07:55 UTC (rev 275)
+++ pkg/RSiena/inst/doc/RSiena_Manual.pdf 2014-06-18 15:07:01 UTC (rev 276)
@@ -50,18 +50,21 @@
endstream
endobj
487 0 obj <<
-/Length 1321
+/Length 1325
/Filter /FlateDecode
>>
stream
-xÚWKsÛ6¾ëW°7j&DˆA²=9Ï&3NÓXí¡I0 Ëh(RáéNýëÝÅ‚ªh+“I¦ãv?|ûœFÛ(^®Ò¯ŒO6«Ç/x Åò"çÑæ&©f©â‘Ö)SѦŽÞÇ—¦] O¦Y'2ñMׯÿܼÍâT“sÍ4,¥¤õîÊÙÖ J#®˜TZ `¢K¦E “œñi¼x8ç’á˜Çïܾ±ˆ’‚e¨ @AV›u!ânG²ì #í«ÖýUÛ~ðz ¯ób©÷Çð!rÊ4î>¤\BxÒjCf,/ùRñ¢û£ªHí÷©‡±(ãnÍc::á)KgÊ’Þ®AÃLk
-‘}ÛÛÊ™ºûÏKœ•YF^âð¡°—ÓYF¿µîQÀB7âä@HÝ
-¿ü
-Ñ©\'Jfñ3»7ý¸³-
-Ž§‚y|5šÑy¶ ‡= »jø‰¶ÞLàm›š4žvMc·Ö[Å%pÒ蘒)0Ïsz=µ–DÅ#BàW”pÏ7+04âyÆ4@ˆœ3]fQµ[}ZI¿I¿~a– _ìxô¬[ý
-÷±gÑdFNN }¦CZ—¬ÔBûDÏ! Ìàœe³)×ÃØ›j<K Ä @µôNH¸€Ü1¾zõüÍų8F¿¥|®i– «hÊK& þAîÊí¦âÔµàW¥ãWí F·=.åPŒkQÄ=í?ßí]ï*_¨ðùÆúùŒ$ö‘6.ZÓ 7<ä)€h)NxrÒp gè³êvûi´áØ}ß%”ô¢,Ë%&/&2¡l{ƒeZèx¼5#Í*Ó÷Î"h‘ÇÝ4Îû–&f'XJ¦Àš½›ÍF•ZßaÙq
- l1“JiA/Qí÷D3Í8¡ @x ò‡jl¼%éÒ„=YOö áUd±©*(6×nésìpTt .Ô‡Öì\E«U]Ÿt`u‹(£§[ÓÞl…mèXùdi ý4ü*õÛèÁÚ»—«è½?î´õ “Öï)Ÿq߈ Uœò{Xñ£ï!@Ù} Lo»“£¹‰,Ul1…JˆIÃîÛÛßi~žv¢8\mýÎ@&§´³@[ý¿´yjÎñ†RC–mn]8ÜßQÀÂÎøܤ› 6Î÷( ÷6OËsý¢•,çÇV2;ú„Ì
-V¾ hÕQeôîz Þ>\Å{SáâÇu±Ý"O.}¥á&‡É²_àÂaíîlƒÍ™€·ÄLþ!oÍ”Tó>#¸W#Ø÷» %<àuÇ‹øÚ÷/àUS û¾ÒmÔz»oLeCL%Óz¶]Ó9Á ßÚf³ ,ß ÖPì¡M`’*ä×â½ÐGÃn¦¶Â68ü€o…2X¨!‡ŽÎšÃa\h,„º$ø¹wãh[’±¦oœí …G¼àÚXÆ™¬É‹`tž/ó6Îç]¡X¦‘ûòE›ó£/|ÃÇ·’ï÷€,1dðP»Æ~iü½»?On .AæÒùÔ»56к¢ïãM`§}MAï·h¶× ÍÖLKóŸÞö˜ªI 1»mÂIŒ YŠ•PÍ ¥P´ò^ÊPÄ0ÞôöÓ³9ÐÚ´Z('¤&ñ ŽÛ]K3»–VL8õ>F0õ†òV!± ”*Oè,=|T˜ ÁgLémkÉ}5 ~>Ú™d¹úzR –Ê/„ ¡çñ˜'xÕŠ\/[,Á#½LKÿïK©¡# lM¬”êòúƒÇë¿ìzEÂ
+xÚWKsÛ6¾ëW°7j&DˆA°=9Ï&3NÓXí¡I0Ëh(R!)§:õ¯w ª¢L&™ŽgØýðípžl’<y¹È¿2>Y- ¿àe"+MÉ“ÕM"rÍrÅs¦’Õ:yŸ^Úv)tº·Í2“¹Hoº~ùçê5hšSMÎ5Ó°”“Ö»+ïZK‚*O¸bRi‚™®˜e`Lr Å÷ã-ÀÃ9— Ç2}çw; @’V ‚ Y–F¤Ý–d/ØFÚWÿkíú!è¼.Í\ïáC.äRUyÚ}ȹ²„ð¤[f0¬- X°²âsÅ‹vÝ Uí@j¿ÿCH=Œ¦J»%OéèŒç,Ÿ(KBx»
+»o|k‰ìÛÞÕÞ®»ÿ¼ÄYUä% J{i˜.
+‚øõwˆú'BênhüåoˆÎúÇe¦d‘>s;Û[×¢àx*X¦W£ }`›qØ㩯‡ŸhëÍ ü¢]³&§]Ó¸ Vq‘
+̪˜ó§×ûÖ‘(7B€cQáAÆ=_-8LÀÒ„—Ó€!JÎtU$õvñi!Ã&ý†…I:.<~±åɳnñ+üÝÇžD³ 9;©y]±J 2½„44˜Â%+&[.®‡±·õx–< ˆ€Œk¼qÉ ƒ|õêù›
+¨9–ÌœOÅñ!/òcÍ@yÅ$@”»òÛ}êZð«Òé«öÎ
+£ß —J¨Æ¥0iOûÏ·;ßû:T*|¾q!G>£‰}¤‹Ö6‡Áy
+ Z‰žœ4ü@ÇYú¬»ín?ºxì®ï2Êza+J‰Ù‹™L(›ÞbŽ·v¤YmûÞ;5eÚíÇißÑdÀôKÉX@³·“Ù¨rCë[¬[“®!ƒ ¦rC9"èà9ª» žhöN¬ ¨þ¡ ›`I>7¡EOšèÉþ#œ L‘Úº†jóí†>ÇGEçáÂúÐÚ¯i²ªë³¬ne t×´7YáúV!Y^AC¿J%ý&y°öîå"y Ž;í}Ǥ
+{*dÜ7bBçü Vúè{ˆPq SëÁîÅähn&+•:L¡
+bÒ°ûöãöwÚ_ž§)w@»þVg “SÚE¤þ_Ú<·çxé!ËV·> î£ï¨@`akCnÒUç{”„‹›çÕ¹›~ÖJ +ù±•LN£>! ÃŒàó€Ö UFﯡ·àõÃUº³5.~\Û
+òä2Tnq˜Ìû. †ÑmÏ6Ø’ xLLä òÖLÁ÷Á½iľ׸m,á¯;nÒëп€×šê8ô5n»¨Ö»]ck{d.™Ö“íšÎ‰–èôÖ5;œ•`ù†€tle°†bm“”‘_‹tv£†ÝìÛÛàð>ªh¡†D ::k
+‡õm¤1 ê’àçÞ£kIÆÙ¾ñ®_AðP€kc g²¦4Ñ貜ç lœÏ;£X¡‘ûòE[ò£/BÃÇÇRè÷€,1dðR»Æ~iý»?ïý \¢Ì¥©wk]¤uEßÇ›Àù@ûš‚ÞoÐì [š™0sóŸÞö˜ªI /1·iâIŒ^YŠUPÍ3¥X´ò^ÊXÄ0ÞôîÓ>&fs µý.¾j¡œšÄ s8nw-ÍTê[Z±ñ ߆Á4jÉ[Fb=(UЙ{ø¨0½‚Ϙһ֑ûÖtøùh’•êëI-X.¿n€†žÇSžáU+J=o±¯ô*¯Âÿ/•†Ž$°5±JFªóë ¯ÿ*bF
endstream
endobj
477 0 obj <<
@@ -342,16 +345,25 @@
endstream
endobj
698 0 obj <<
-/Length 1523
+/Length 1525
/Filter /FlateDecode
>>
stream
-xÚÍWKsÛ6¾ûWðHÍD0@ð™žâ6ñ$ÓL¬¶3Mz IHB͇† ìøÔ¿Þ],HKu&éE"—ûÂî~» î<î]_p÷µ¸¸|JO„L†qà-–^š1ÎS/á!2ö¥÷Åÿ4“ÂÏ»~6—Iä¿Ÿýµø`ÅÎâ8ˆPŒ{ó dw7ºÑu^9åû® …`Yƒ`ʸÌHP01› Dþµ~ÐÍŠ4tj©:ÕÊŠÔ ŒebÐÀY”…¤áµjfóÇþÖ ²ìJÄh>A+pûþíÇ7z%œŽ'×+Ò©{üO|mè=ßl:Uè¼W%úuÞÓÓÓ,HüvK/öŽ£%ýzPRçÍ,Hý-Ì*mJÇÒ¡m”{èèCÝv
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 276
More information about the Rsiena-commits
mailing list