[Rsiena-commits] r339 - in pkg: RSiena RSiena/R RSiena/data RSiena/inst/scripts RSiena/man RSiena/src RSiena/src/model/effects RSiena/src/model/ml RSiena/tests RSienaTest RSienaTest/R RSienaTest/data RSienaTest/doc RSienaTest/inst/scripts RSienaTest/man RSienaTest/src RSienaTest/src/model/effects RSienaTest/tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 26 20:35:39 CET 2019
Author: tomsnijders
Date: 2019-02-26 20:35:38 +0100 (Tue, 26 Feb 2019)
New Revision: 339
Added:
pkg/RSiena/man/print.sienaTest.Rd
pkg/RSiena/src/model/effects/AltersCovariateMaximumEffect.cpp
pkg/RSiena/src/model/effects/AltersCovariateMaximumEffect.h
pkg/RSiena/src/model/effects/AltersCovariateMinimumEffect.cpp
pkg/RSiena/src/model/effects/AltersCovariateMinimumEffect.h
pkg/RSienaTest/src/model/effects/AltersCovariateMaximumEffect.cpp
pkg/RSienaTest/src/model/effects/AltersCovariateMaximumEffect.h
pkg/RSienaTest/src/model/effects/AltersCovariateMinimumEffect.cpp
pkg/RSienaTest/src/model/effects/AltersCovariateMinimumEffect.h
Modified:
pkg/RSiena/ChangeLog
pkg/RSiena/DESCRIPTION
pkg/RSiena/NAMESPACE
pkg/RSiena/R/Sienatest.r
pkg/RSiena/R/phase1.r
pkg/RSiena/R/phase3.r
pkg/RSiena/R/print07Report.r
pkg/RSiena/R/robmon.r
pkg/RSiena/R/siena07.r
pkg/RSiena/R/sienaRI.r
pkg/RSiena/R/sienaeffects.r
pkg/RSiena/R/sienaprint.r
pkg/RSiena/R/sienatable.r
pkg/RSiena/data/allEffects.csv
pkg/RSiena/inst/scripts/RSienaSNADescriptives.R
pkg/RSiena/inst/scripts/Rscript01DataFormat.R
pkg/RSiena/inst/scripts/Rscript02SienaVariableFormat.R
pkg/RSiena/inst/scripts/Rscript03SienaRunModel.R
pkg/RSiena/inst/scripts/Rscript04SienaBehaviour.R
pkg/RSiena/man/RSiena-package.Rd
pkg/RSiena/man/Wald.Rd
pkg/RSiena/man/getEffects.Rd
pkg/RSiena/man/includeEffects.Rd
pkg/RSiena/man/print.sienaEffects.Rd
pkg/RSiena/man/print.sienaMeta.Rd
pkg/RSiena/man/setEffect.Rd
pkg/RSiena/man/siena07.Rd
pkg/RSiena/man/siena08.Rd
pkg/RSiena/man/sienaAlgorithmCreate.Rd
pkg/RSiena/man/sienaFit.Rd
pkg/RSiena/man/sienaGOF.Rd
pkg/RSiena/man/sienaTimeTest.Rd
pkg/RSiena/man/xtable.Rd
pkg/RSiena/src/model/effects/AllEffects.h
pkg/RSiena/src/model/effects/CovariateAndNetworkBehaviorEffect.cpp
pkg/RSiena/src/model/effects/CovariateAndNetworkBehaviorEffect.h
pkg/RSiena/src/model/effects/EffectFactory.cpp
pkg/RSiena/src/model/effects/IsolateEffect.cpp
pkg/RSiena/src/model/effects/IsolateEffect.h
pkg/RSiena/src/model/effects/SameCovariateActivityEffect.cpp
pkg/RSiena/src/model/effects/SameCovariateActivityEffect.h
pkg/RSiena/src/model/ml/MLSimulation.cpp
pkg/RSiena/src/sources.list
pkg/RSiena/tests/parallel.R
pkg/RSiena/tests/parallel.Rout.save
pkg/RSienaTest/ChangeLog
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/NAMESPACE
pkg/RSienaTest/R/Sienatest.r
pkg/RSienaTest/R/bayesTest.r
pkg/RSienaTest/R/phase1.r
pkg/RSienaTest/R/phase3.r
pkg/RSienaTest/R/print07Report.r
pkg/RSienaTest/R/robmon.r
pkg/RSienaTest/R/siena07.r
pkg/RSienaTest/R/sienaBayes.r
pkg/RSienaTest/R/sienaRI.r
pkg/RSienaTest/R/sienaeffects.r
pkg/RSienaTest/R/sienaprint.r
pkg/RSienaTest/R/sienatable.r
pkg/RSienaTest/data/allEffects.csv
pkg/RSienaTest/doc/HowToCommit.pdf
pkg/RSienaTest/doc/HowToCommit.tex
pkg/RSienaTest/doc/RSiena.bib
pkg/RSienaTest/doc/RSiena_Manual.pdf
pkg/RSienaTest/doc/RSiena_Manual.tex
pkg/RSienaTest/inst/scripts/RSienaSNADescriptives.R
pkg/RSienaTest/inst/scripts/Rscript01DataFormat.R
pkg/RSienaTest/inst/scripts/Rscript02SienaVariableFormat.R
pkg/RSienaTest/inst/scripts/Rscript03SienaRunModel.R
pkg/RSienaTest/inst/scripts/Rscript04SienaBehaviour.R
pkg/RSienaTest/man/RSiena-package.Rd
pkg/RSienaTest/man/Wald.Rd
pkg/RSienaTest/man/extract.sienaBayes.Rd
pkg/RSienaTest/man/includeEffects.Rd
pkg/RSienaTest/man/print.sienaBayesFit.Rd
pkg/RSienaTest/man/print.sienaEffects.Rd
pkg/RSienaTest/man/print.sienaMeta.Rd
pkg/RSienaTest/man/setEffect.Rd
pkg/RSienaTest/man/siena07.Rd
pkg/RSienaTest/man/siena08.Rd
pkg/RSienaTest/man/sienaBayes.Rd
pkg/RSienaTest/man/sienaFit.Rd
pkg/RSienaTest/man/sienaGOF.Rd
pkg/RSienaTest/man/sienaTimeTest.Rd
pkg/RSienaTest/man/xtable.Rd
pkg/RSienaTest/src/model/effects/AllEffects.h
pkg/RSienaTest/src/model/effects/CovariateAndNetworkBehaviorEffect.cpp
pkg/RSienaTest/src/model/effects/CovariateAndNetworkBehaviorEffect.h
pkg/RSienaTest/src/model/effects/EffectFactory.cpp
pkg/RSienaTest/src/model/effects/IsolateEffect.cpp
pkg/RSienaTest/src/model/effects/IsolateEffect.h
pkg/RSienaTest/src/sources.list
pkg/RSienaTest/tests/parallel.R
pkg/RSienaTest/tests/parallel.Rout.save
Log:
Version 1.2-16 of both packages
Modified: pkg/RSiena/ChangeLog
===================================================================
--- pkg/RSiena/ChangeLog 2019-01-15 16:54:20 UTC (rev 338)
+++ pkg/RSiena/ChangeLog 2019-02-26 19:35:38 UTC (rev 339)
@@ -1,3 +1,49 @@
+2019-02-26 R-Forge Revision 339, packages version 1.2-16.
+Changes in RSiena and RSienaTest:
+ * New effects maxAAlt, minXAlt (thanks to Per Block, Marion Hoffman,
+ Isabel Raabe, and Timon Elmer), outIsolate.
+ * Effect name of isolate effect (not shortName) changed to in-isolate.
+ * New parameter thetaValues in siena07, allowing simulations in Phase 3
+ with varying parameters according to matrix input.
+ Corresponding changes in print.sienaFit.
+ * sienaRI: earlier check for bipartite dependent variables; if so, stop.
+ * Argument "verbose" added to includeEffects and setEffect.
+ (It existed already for includeInteraction.)
+ * Wald.RSiena, Multipar.RSiena and score.Test also produce one-sided
+ test statistic if df=1.
+ * objects produced by Wald.RSiena, Multipar.RSiena and score.Test
+ have class "sienaTest".
+ * print method added for objects of class sienaTest.
+ * Changed/added some keywords and concepts in .Rd files.
+Changes in RSienaTest:
+ * Allow arbitrary number of interactions in sienaBayes (via getEffects).
+ * simpleBayesTest and multipleBayesTest corrected; they were wrong
+ for models with user-defined interactions.
+ * glueBayes extended (to allow good use for models with fixed parameters)
+ and check for prior rates relaxed (to allow priorRatesFromData=2
+ leading to different prior rates.)
+ * New default (prevBayes$nwarm >= 1) for newProposalFromPrev in sienaBayes.
+ * siena.table for sienaBayes results extended with between-groups s.d.
+ * Argument "verbose" added to extract.posteriorMeans.
+ * sienaBayes: check that prevAns is a sienaFit object;
+ totally skip initial multigroup estimation if prevAns is given with
+ prevAns$n3 >= 500 and (!usePrevOnly).
+
+
+2019-01-02 R-Forge Revision 338, packages version 1.2-15.
+Changes in RSiena and RSienaTest:
+ * New effects sameXReciAct and diffXReciAct.
+ * Increased length of MLSimulation::preburnin.
+ * Changes in sienaFitThetaTable, averageTheta.last, sdTheta.last in sienaprint.r
+ to enable change in extract.posteriorMeans (see below).
+ * siena.table adapted for sienaBayes results.
+Changes in RSienaTest:
+ * New functions shortBayesResults and plotPostMeansMDS.
+ * extract.posteriorMeans works also for an object saved as PartialBayesResult
+ (which has NA results for as yet unfinished runs)
+ (this was achieved by changes in sienaprint.r, see above).
+ * Examples added to help page for print.sienaBayesFit (in dontrun...).
+
2018-12-05 R-Forge Revision 337, packages version 1.2-14.
Changes in RSiena and RSienaTest:
* Check for length of parameter mult for ML estimation (initializeFRAN),
@@ -12,7 +58,9 @@
prior mean is used for initial parameter values;
prewarming phase introduced before improveMH;
in case prevBayes is used, also parameter nImproveMH in the function
- call of sienaBayes supersedes this parameter in the prevBayes object.
+ call of sienaBayes supersedes this parameter in the prevBayes object;
+ checks for dimensions of prior mean and covariance matrix
+ put earlier in the initialization phase.
2018-10-29 R-Forge Revision 336, packages version 1.2-13.
Changes in RSiena and RSienaTest:
Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION 2019-01-15 16:54:20 UTC (rev 338)
+++ pkg/RSiena/DESCRIPTION 2019-02-26 19:35:38 UTC (rev 339)
@@ -2,8 +2,8 @@
Package: RSiena
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.2-14
-Date: 2018-12-05
+Version: 1.2-16
+Date: 2019-02-26
Author: Ruth Ripley, Krists Boitmanis, Tom A.B. Snijders, Felix Schoenenberger
Depends: R (>= 2.15.0)
Imports: Matrix, tcltk, lattice, parallel, MASS, methods
Modified: pkg/RSiena/NAMESPACE
===================================================================
--- pkg/RSiena/NAMESPACE 2019-01-15 16:54:20 UTC (rev 338)
+++ pkg/RSiena/NAMESPACE 2019-02-26 19:35:38 UTC (rev 339)
@@ -57,6 +57,7 @@
S3method(addAttributes, varCovar)
S3method(addAttributes, coDyadCovar)
S3method(addAttributes, varDyadCovar)
+S3method(print, sienaTest)
S3method(print, sienaTimeTest)
S3method(summary, sienaTimeTest)
S3method(print, summary.sienaTimeTest)
Modified: pkg/RSiena/R/Sienatest.r
===================================================================
--- pkg/RSiena/R/Sienatest.r 2019-01-15 16:54:20 UTC (rev 338)
+++ pkg/RSiena/R/Sienatest.r 2019-02-26 19:35:38 UTC (rev 339)
@@ -325,10 +325,21 @@
if (any(test & (!ans$fix))) cat('Warning: some tested parameters were not fixed; do you know what you are doing??? \n')
fra <- colMeans(ans$sf, na.rm=TRUE)
redundant <- (ans$fix & (!test))
- teststat <- EvaluateTestStatistic(ans$maxlike, test, redundant, ans$dfra, ans$msf, fra)$cvalue
+ tests <- EvaluateTestStatistic(ans$maxlike, test, redundant, ans$dfra, ans$msf, fra)
+ teststat <- tests$cvalue
df <- sum(test)
+ if (df == 1)
+ {
+ onesided <- tests$oneSided
+ }
+ else
+ {
+ onesided <- NULL
+ }
pval <- 1 - pchisq(teststat, df)
- list(chisquare = teststat, df = df, pvalue = pval)
+ t.ans <- list(chisquare=teststat, df=df, pvalue=pval, onesided=onesided)
+ class(t.ans) <- "sienaTest"
+ t.ans
}
##@Wald.RSiena Calculate Wald test statistics
@@ -351,9 +362,19 @@
th <- A %*% ans$theta
covmat <- A %*% sigma %*% t(A)
chisq <- drop(t(th) %*% solve(covmat) %*% th)
- d.f. <- nrow(A)
- pval <- 1 - pchisq(chisq, d.f.)
- list(chisquare = chisq, df = d.f., pvalue = pval)
+ df <- nrow(A)
+ pval <- 1 - pchisq(chisq, df)
+ if (df == 1)
+ {
+ onesided <- sign(th) * sqrt(chisq)
+ }
+ else
+ {
+ onesided <- NULL
+ }
+ t.ans <- list(chisquare=chisq, df=df, pvalue=pval, onesided=onesided)
+ class(t.ans) <- "sienaTest"
+ t.ans
}
##@Multipar.RSiena Calculate Wald test statistic for hypothesis that subvector = 0.
@@ -365,3 +386,30 @@
A[cbind(1:k,c(...))] <- 1
Wald.RSiena(A, ans)
}
+
+##@print.sienaTest Methods
+print.sienaTest <- function(x, ...)
+{
+ if (!inherits(x, "sienaTest"))
+ {
+ stop("not a legitimate sienaTest object")
+ }
+ cat(paste('chi-squared = ',
+ format(round(x$chisquare, digits=2), nsmall = 2),
+ ', d.f. = ', x$df, '; ', sep=''))
+ if ((x$df == 1) & (!is.null(x$onesided)))
+ {
+ cat(paste('one-sided Z = ',
+ format(round(x$onesided, digits=2), nsmall = 2), '; ', sep=''))
+ }
+ if (x$pvalue < 0.001)
+ {
+ cat(' p < 0.001. \n')
+ }
+ else
+ {
+ cat(paste(' p = ',
+ format(round(x$pvalue, digits=3), nsmall = 2), '. \n', sep=''))
+ }
+ invisible(x)
+}
Modified: pkg/RSiena/R/phase1.r
===================================================================
--- pkg/RSiena/R/phase1.r 2019-01-15 16:54:20 UTC (rev 338)
+++ pkg/RSiena/R/phase1.r 2019-02-26 19:35:38 UTC (rev 339)
@@ -554,5 +554,9 @@
z$ntim <- matrix(NA, nrow=nIterations, ncol=f$observations - 1)
}
z$sims <- vector("list", nIterations)
+ if (z$thetaFromFile)
+ {
+ z$thetaUsed <- matrix(NA, nrow = nIterations , ncol = z$pp)
+ }
z
}
Modified: pkg/RSiena/R/phase3.r
===================================================================
--- pkg/RSiena/R/phase3.r 2019-01-15 16:54:20 UTC (rev 338)
+++ pkg/RSiena/R/phase3.r 2019-02-26 19:35:38 UTC (rev 339)
@@ -589,6 +589,14 @@
}
}
## iteration proper
+
+ if (z$thetaFromFile)
+ {
+ if (nit <= dim(z$thetaValues)[1])
+ {
+ zsmall$theta <- z$thetaValues[nit,]
+ }
+ }
if (z$int == 1)
{
zz <- x$FRAN(zsmall, xsmall)
@@ -634,6 +642,10 @@
z$sims[[z$nit]] <- zz$sims
z$chain[[z$nit]] <- zz$chain
fra <- fra + z$targets
+ if (z$thetaFromFile)
+ {
+ z$thetaUsed[z$nit, ] <- zsmall$theta
+ }
}
else
{
@@ -651,6 +663,10 @@
z$sf2s <- z$sf2s + zz[[i]]$fra
}
z$sims[[z$nit + (i - 1)]] <- zz[[i]]$sims
+ if (z$thetaFromFile)
+ {
+ z$thetaUsed[z$nit + (i - 1), ] <- zsmall$theta
+ }
}
if (z$FinDiff.method)
{
Modified: pkg/RSiena/R/print07Report.r
===================================================================
--- pkg/RSiena/R/print07Report.r 2019-01-15 16:54:20 UTC (rev 338)
+++ pkg/RSiena/R/print07Report.r 2019-02-26 19:35:38 UTC (rev 339)
@@ -173,6 +173,11 @@
if (x$simOnly)
{
Heading(3, outf, "Simulated statistics")
+ if (z$thetaFromFile)
+ {
+ Report('Note: Parameters in Phase 3 were varying. \n', bof)
+ Report('Note: Parameters in Phase 3 were varying. \n', outf)
+ }
Report('Estimated means and standard deviations, standard errors of the mean \n', bof)
Report('Estimated means and standard deviations, standard errors of the mean \n', outf)
dmsf <- diag(z$msf)
@@ -181,7 +186,7 @@
mean.stats <- colMeans(z$sf) + z$targets
# cov.dev <- z$msf
sem <- sqrt(dmsf/dim(z$sf)[1])
- if (x$dolby)
+ if ((x$dolby) & (!z$thetaFromFile))
{
scores <- apply(z$ssc, c(1,3), sum) # z$nit by z$pp matrix
mean.scores <- colMeans(scores)
@@ -195,7 +200,7 @@
format(round(sem, 4) ,width=8, nsmall=4), sep='')
PrtOutMat(as.matrix(mymess1), outf)
PrtOutMat(as.matrix(mymess1), bof)
- if (x$dolby)
+ if ((x$dolby) & (!z$thetaFromFile))
{
Report('Standard errors of the mean are less than s.d./sqrt(n) \n', outf)
Report('because of regression on scores (Dolby option). \n', outf)
Modified: pkg/RSiena/R/robmon.r
===================================================================
--- pkg/RSiena/R/robmon.r 2019-01-15 16:54:20 UTC (rev 338)
+++ pkg/RSiena/R/robmon.r 2019-02-26 19:35:38 UTC (rev 339)
@@ -120,6 +120,16 @@
{
z$cl <- NULL
}
+
+ # check dimensionality of thetaValues
+ if (z$thetaFromFile)
+ {
+ if (dim(z$thetaValues)[2] != z$pp)
+ {
+ stop(paste('Matrix thetaValues, if given, should have', z$pp, 'columns.'))
+ }
+ }
+
z$newFixed <- rep(FALSE, z$pp)
z$AllNowFixed <- FALSE
if (!z$haveDfra)
@@ -353,8 +363,23 @@
z <- terminateFRAN(z, x)
## #####################################################
## call to FRAN changes covariance matrix for conditional estimation
- if (!x$simOnly)
+ if (x$simOnly)
{
+ if (z$thetaFromFile)
+ {
+ z$theta <- colMeans(z$thetaUsed)
+ z$covtheta <- cov(z$thetaUsed)
+ z$se <- sqrt(diag(z$covtheta))
+ z$thetaValues <- NULL # not needed any longer, superseded by z$thetaUsed
+ }
+ else
+ {
+ z$covtheta <- matrix(0, z$pp, z$pp )
+ z$se <- rep(0, z$pp)
+ }
+ }
+ else
+ {
z$diver<- (z$fixed | z$diver | diag(z$covtheta) < 1e-9) & (!z$AllUserFixed)
z$covtheta[z$diver, ] <- NA # was Root(diag(z$covtheta)) * 33
##not sure this does not use very small vals
Modified: pkg/RSiena/R/siena07.r
===================================================================
--- pkg/RSiena/R/siena07.r 2019-01-15 16:54:20 UTC (rev 338)
+++ pkg/RSiena/R/siena07.r 2019-02-26 19:35:38 UTC (rev 339)
@@ -12,7 +12,9 @@
##@siena07 siena07
siena07 <- function(x, batch = FALSE, verbose = FALSE, silent=FALSE,
- useCluster = FALSE, nbrNodes = 2, initC=TRUE,
+ useCluster = FALSE, nbrNodes = 2,
+ thetaValues = NULL,
+ initC=TRUE,
clusterString=rep("localhost", nbrNodes), tt=NULL,
parallelTesting=FALSE, clusterIter=!x$maxlike,
clusterType=c("PSOCK", "FORK"), cl=NULL, ...)
@@ -169,6 +171,29 @@
z$pb <- list(pb=NULL, pbval=0, pbmax=1)
}
+ ## create theta values for phase 3, if necessary
+ if (is.null(thetaValues))
+ {
+ z$thetaValues <- NULL
+ z$thetaFromFile <- FALSE
+ }
+ else
+ {
+ if (!x$simOnly)
+ {
+ cat('The thetaValues parameter was given\n')
+ cat('but simOnly was not specified in the algorithm object.\n')
+ cat('This is inconsistent.\n')
+ stop('Inconsistent combination simOnly - thetaValues')
+ }
+ if (!is.matrix(thetaValues))
+ {
+ stop('thetaValues should be a matrix.')
+ }
+ z$thetaValues <- thetaValues
+ z$thetaFromFile <- TRUE
+ }
+
z <- robmon(z, x, useCluster, nbrNodes, initC, clusterString,
clusterIter, clusterType, cl, ...)
Modified: pkg/RSiena/R/sienaRI.r
===================================================================
--- pkg/RSiena/R/sienaRI.r 2019-01-15 16:54:20 UTC (rev 338)
+++ pkg/RSiena/R/sienaRI.r 2019-02-26 19:35:38 UTC (rev 339)
@@ -17,6 +17,11 @@
{
stop("no a legitimate Siena data specification")
}
+ datatypes <- sapply(data$depvars, function(x){attr(x,"type")})
+ if (any(datatypes == "bipartite"))
+ {
+ stop("sienaRI works only for dependent variables of type 'oneMode' or 'behavior'")
+ }
if(!is.null(ans))
{
if (!inherits(ans, "sienaFit"))
@@ -25,7 +30,7 @@
}
if(!is.null(algorithm)||!is.null(theta)||!is.null(effects))
{
- warning(paste("some information are multiply defined \n",
+ warning(paste("some informations are multiply defined \n",
"results will be based on 'theta', 'algorithm', and 'effects'\n",
"stored in 'ans' (as 'ans$theta', 'ans$x', 'ans$effects')\n", sep=""))
}
Modified: pkg/RSiena/R/sienaeffects.r
===================================================================
--- pkg/RSiena/R/sienaeffects.r 2019-01-15 16:54:20 UTC (rev 338)
+++ pkg/RSiena/R/sienaeffects.r 2019-02-26 19:35:38 UTC (rev 339)
@@ -11,7 +11,7 @@
includeEffects <- function(myeff, ..., include=TRUE, name=myeff$name[1],
type="eval", interaction1="", interaction2="",
fix=FALSE, test=FALSE,
- character=FALSE)
+ character=FALSE, verbose=TRUE)
{
if (!inherits(myeff, 'sienaEffects'))
{
@@ -63,7 +63,11 @@
{
# print.data.frame(myeff[use, c("name", "shortName", "type",
# "interaction1", "interaction2", "include")])
- print.sienaEffects(myeff[use,])
+ if (verbose)
+ {
+ print.sienaEffects(myeff[use,], includeRandoms =
+ any(myeff$random & (myeff$shortName != 'density')))
+ }
}
if (hasArg('initialValue'))
{
@@ -233,7 +237,8 @@
}
if (verbose)
{
- print.sienaEffects(myeff[intn,], includeRandoms = random)
+ print.sienaEffects(myeff[intn,], includeRandoms =
+ any(myeff$random & (myeff$shortName != 'density')))
}
myeff
}
@@ -246,7 +251,7 @@
include=TRUE, name=myeff$name[1],
type="eval", interaction1="", interaction2="",
effect1=0, effect2=0, effect3=0,
- period=1, group=1, character=FALSE)
+ period=1, group=1, character=FALSE, verbose=TRUE)
{
if (!character)
{
@@ -298,6 +303,9 @@
# print.data.frame(myeff[use, c("name", "shortName", "type", "interaction1",
# "interaction2", "include", "parm", "fix", "test",
# "initialValue", "timeDummy", "period", "group")])
- print.sienaEffects(myeff[use,])
+ if (verbose)
+ {
+ print.sienaEffects(myeff[use,], includeRandoms = random)
+ }
myeff
}
Modified: pkg/RSiena/R/sienaprint.r
===================================================================
--- pkg/RSiena/R/sienaprint.r 2019-01-15 16:54:20 UTC (rev 338)
+++ pkg/RSiena/R/sienaprint.r 2019-02-26 19:35:38 UTC (rev 339)
@@ -235,6 +235,7 @@
##@print.sienaFit Methods
print.sienaFit <- function(x, tstat=TRUE, ...)
{
+ objectName <- deparse(substitute(x))
if (!inherits(x, "sienaFit"))
{
stop("not a legitimate Siena model fit")
@@ -249,9 +250,13 @@
}
else
{
+ if (is.null(x$thetaFromFile))
+ {
+ x$thetaFromFile <- FALSE
+ }
if(x$x$simOnly)
{
- cat("Parameter values\n\n")
+ cat("Parameter values used for simulations\n\n")
}
else
{
@@ -274,15 +279,14 @@
mymat[, 'value'] <- format(round(mydf$value, digits=4))
if(x$x$simOnly)
{
- mymat[, 'se'] <- ' '
mymat[, 'tstat'] <- ' '
}
else
{
- mymat[, 'se'] <- format(round(mydf$se, digits=4))
mymat[, 'tstat'] <- paste(" ", format(round(mydf$tstat, digits=4)))
mymat[is.na(mydf$tstat), 'tstat'] <- ' '
}
+ mymat[, 'se'] <- format(round(mydf$se, digits=4))
mymat[, 'type'] <- format(mymat[, 'type'])
mymat[, 'text'] <- format(mymat[, 'text'])
mymat[mydf$row < 1, 'row'] <-
@@ -291,9 +295,18 @@
paste(format(mydf[mydf$row >= 1, 'row']), '.', sep='')
if(x$x$simOnly)
{
- mymat <- rbind(c(rep("", 4), "Parameter", "", "", "",
+ if (x$thetaFromFile)
+ {
+ mymat <- rbind(c(rep("", 4), "Mean", "", "Standard", "",
""),
- c(rep("", 4), " value",rep("", 4)), mymat)
+ c(rep("", 4), " value", "", "Deviation", rep("", 2)), mymat)
+ }
+ else
+ {
+ mymat <- rbind(c(rep("", 4), "Parameter", "", "Standard", "",
+ ""),
+ c(rep("", 4), " value", "", "Deviation", rep("", 2)), mymat)
+ }
}
else
{
@@ -323,14 +336,22 @@
if(x$x$simOnly)
{
- cat('\nEstimated means and standard deviations, standard errors of the mean \n')
+ cat('\nSimulated means and standard deviations')
+ if (x$thetaFromFile)
+ {
+ cat('\n')
+ }
+ else
+ {
+ cat(', standard errors of the mean \n')
+ }
if (x$x$maxlike)
{
cat('Note: statistics for ML simulations are score functions.\n')
}
dmsf <- diag(x$msf)
mean.stats <- colMeans(x$sf) + x$targets
-# cov.dev may be dropped - just for now (07-10-13) I keep it in
+# cov.dev may be dropped
# cov.dev <- x$msf
sem <- sqrt(dmsf/dim(x$sf)[1])
if (x$x$dolby)
@@ -355,10 +376,29 @@
mymess1 <- paste(format(1:x$qq,width=3), '. ',
format(x$requestedEffects$functionName[selection], width = 56),
format(round(mean.stats, 3), width=8, nsmall=3), ' ',
- format(round(sqrt(dmsf), 3) ,width=8, nsmall=3), ' ',
+ format(round(sqrt(dmsf), 3) ,width=8, nsmall=3), ' ', sep='')
+ if (x$thetaFromFile)
+ {
+ mymess1 <- paste(mymess1, rep('\n',x$qq), sep='')
+ }
+ else
+ {
+ mymess1 <- paste(mymess1,
format(round(sem, 4) ,width=8, nsmall=4),
rep('\n',x$qq), sep='')
+ }
cat(as.matrix(mymess1),'\n', sep='')
+ cat("\nSimulated statistics are in ...$sf")
+# paste(objectName,'$sf',sep=""))
+ if (x$returnDeps)
+ {
+ cat("\nand simulated dependent variables in ...$sims.\n")
+# paste(objectName,'$sims',sep=""), ".\n")
+ }
+ else
+ {
+ cat(".\n")
+ }
}
else
{
@@ -380,40 +420,52 @@
cat(names(x$x$MaxDegree)[i],':',x$x$MaxDegree[i],'\n')}
cat('\n')
}
- if (any(x$x$UniversalOffset > 0)) {
+ if (any(x$x$UniversalOffset > 0))
+ {
cat(' Offsets for symmetric networks, and for settings model (if any): \n')
- for (i in 1:length(x$x$UniversalOffset)){
- cat(names(x$x$UniversalOffset)[i],':',x$x$UniversalOffset[i],'\n')}
+ for (i in 1:length(x$x$UniversalOffset))
+ {
+ cat(names(x$x$UniversalOffset)[i],':',x$x$UniversalOffset[i],'\n')
+ }
}
- if (any(x$x$modelType != 1)) {
+ if (any(x$x$modelType != 1))
+ {
cat('\n Model Type:\n')
- for (i in 1:length(x$x$modelType)){
+ for (i in 1:length(x$x$modelType))
+ {
cat(names(x$x$modelType)[i],':',
- ModelTypeStrings(x$x$modelType[i]),'\n')}
+ ModelTypeStrings(x$x$modelType[i]),'\n')
+ }
cat('\n')
}
- if (any(x$x$behModelType != 1)) {
+ if (any(x$x$behModelType != 1))
+ {
cat('\n Behavioral Model Type:\n')
- for (i in 1:length(x$x$behModelType)){
+ for (i in 1:length(x$x$behModelType))
+ {
cat(names(x$x$behModelType)[i],':',
- BehaviorModelTypeStrings(x$x$behModelType[i]),'\n')}
+ BehaviorModelTypeStrings(x$x$behModelType[i]),'\n')
+ }
cat('\n')
}
try(if (x$errorMessage.cov > '')
{cat('\nWarning:', x$errorMessage.cov, '\n')}, silent=TRUE)
- # "Try" for downward compatibility
+ # "Try" for compatibility with previous versions
cat("\nTotal of", x$n, "iteration steps.\n\n")
- if ((x$x$dolby)&(x$x$simOnly))
+
+ if ((x$x$dolby) & (x$x$simOnly) & (!x$thetaFromFile))
{
- cat('(Standard errors of the mean are less than s.d./', sqrt(x$n),' \n',
- sep='')
+ cat('(Standard errors of the mean are less than s.d./',
+ sqrt(x$n),' \n', sep='')
cat('because of regression on scores (Dolby option).) \n')
}
if (x$termination == "UserInterrupt")
+ {
cat(" \n*** Warning ***",
"Estimation terminated early at user request.\n")
}
+ }
invisible(x)
}
@@ -703,7 +755,7 @@
averageTheta.last <- function(z, groupOnly=0, nfirst=z$nwarm+1)
{
ntot <- sum(!is.na(z$ThinPosteriorMu[,1]))
- ntott <- dim(z$ThinParameters)[1]
+ ntott <- sum(!is.na(z$ThinParameters[,1,1]))
if (nfirst > ntot)
{
@@ -722,7 +774,8 @@
!z$generalParametersInGroup, drop=FALSE], na.rm=TRUE)
postVarMean[z$ratePositions[[group]]] <- apply(
z$ThinParameters[nfirst:ntott,
- group, !z$generalParametersInGroup, drop=FALSE], 3, var, na.rm=TRUE)
+ group, !z$generalParametersInGroup, drop=FALSE], 3,
+ var, na.rm=TRUE)
}
if (is.null(z$priorRatesFromData))
@@ -733,7 +786,7 @@
if (groupOnly != 0)
{
thetaMean[(z$set1)&(!z$basicRate)] <- colMeans(
- z$ThinParameters[(nfirst):ntott, groupOnly,
+ z$ThinParameters[nfirst:ntott, groupOnly,
z$varyingGeneralParametersInGroup, drop=FALSE], na.rm=TRUE)
}
else
@@ -741,24 +794,24 @@
if ((z$priorRatesFromData <0) | z$incidentalBasicRates)
{
thetaMean[(z$set1)&(!z$basicRate)] <- colMeans(
- z$ThinPosteriorMu[(nfirst):dim(z$ThinPosteriorMu)[1],
- , drop=FALSE], na.rm=TRUE)
- postVarMean[z$varyingObjectiveParameters] <- sapply(1:(dim(z$ThinPosteriorSigma)[2]),
- function(i){mean(z$ThinPosteriorSigma[nfirst:dim(z$ThinPosteriorSigma)[1],i,i], na.rm=TRUE)})
+ z$ThinPosteriorMu[nfirst:ntot, , drop=FALSE], na.rm=TRUE)
+ postVarMean[z$varyingObjectiveParameters] <-
+ sapply(1:(dim(z$ThinPosteriorSigma)[2]),
+ function(i){mean(z$ThinPosteriorSigma[nfirst:ntot,i,i], na.rm=TRUE)})
}
else
{
thetaMean[(z$set1)&(!z$basicRate)] <- colMeans(
- z$ThinPosteriorMu[(nfirst):dim(z$ThinPosteriorMu)[1],
+ z$ThinPosteriorMu[nfirst:ntot,
z$objectiveInVarying, drop=FALSE], na.rm=TRUE)
- postVarMean[z$varyingObjectiveParameters] <- sapply(1:(dim(z$ThinPosteriorSigma)[2]),
- function(i){mean(z$ThinPosteriorSigma[nfirst:dim(z$ThinPosteriorSigma)[1],i,i], na.rm=TRUE)}
+ postVarMean[z$varyingObjectiveParameters] <-
+ sapply(1:(dim(z$ThinPosteriorSigma)[2]),
+ function(i){mean(z$ThinPosteriorSigma[nfirst:ntot,i,i], na.rm=TRUE)}
)[z$objectiveInVarying]
}
}
thetaMean[z$set2] <-
- colMeans(z$ThinPosteriorEta[nfirst:dim(z$ThinPosteriorEta)[1],, drop=FALSE],
- na.rm=TRUE)
+ colMeans(z$ThinPosteriorEta[nfirst:ntot,, drop=FALSE], na.rm=TRUE)
thetaMean[z$fix & (!z$basicRate)] <- z$thetaMat[1,z$fix & (!z$basicRate)]
list(thetaMean, postVarMean)
}
@@ -767,7 +820,7 @@
sdTheta.last <- function(z, groupOnly=0, nfirst=z$nwarm+1)
{
ntot <- sum(!is.na(z$ThinPosteriorMu[,1]))
- ntott <- dim(z$ThinParameters)[1]
+ ntott <- sum(!is.na(z$ThinParameters[,1,1]))
if (nfirst >= ntot-1)
{
stop('Sample did not come beyond warming')
@@ -877,6 +930,10 @@
theEffects <-theEffects[which(theEffects$type!='gmm'),]
}
pp <- dim(theEffects)[1]
+ if (is.null(x$thetaFromFile))
+ {
+ x$thetaFromFile <- FALSE
+ }
if (x$cconditional)
{
nrates <- length(x$rate)
@@ -947,16 +1004,36 @@
{
mydf[1, 'text'] <- 'Rate parameter of conditioning variable'
}
+ if (x$x$simOnly)
+ {
+ mydf[1, 'se'] <- NA
+ mydf[1, 'value'] <- NA
+ }
+ else
+ {
mydf[1, 'value'] <- x$rate[1]
mydf[1, 'se'] <- x$vrate[1]
}
+ addtorow$command[addsub] <-
+ 'Other parameters: '
+ addtorow$pos[[addsub]] <- 3
+ addsub <- addsub + 1
+ }
else ## observations > 2
{
nn <- length(x$rate)
nnstr <- as.numeric(paste('0.', as.character(1:nn), sep=""))
mydf[1:nn, 'row'] <- nnstr
+ if (x$x$simOnly)
+ {
+ mydf[1:nn, 'se'] <- rep(NA, nn)
+ mydf[1:nn, 'value'] <- rep(NA, nn)
+ }
+ else
+ {
mydf[1:nn, 'value'] <- x$rate
mydf[1:nn, 'se'] <- x$vrate
+ }
if (length(x$f$types) == 1)
{
mydf[1:nn, 'text'] <- paste('Rate parameter period', 1:nn)
@@ -1260,7 +1337,15 @@
}
else
{
- cat("Note: this summary does not contain a convergence check.\n\n")
+ cat("Note: this summary does not contain a convergence check.\n")
+ if (is.null(nfirst))
+ {
+ cat("Note: the print function for sienaBayesFit objects")
+ cat(" can also use a parameter nfirst,\n")
+ cat(" indicating the first run")
+ cat(" from which convergence is assumed.\n")
+ }
+ cat("\n")
if (length(x$f$groupNames) > 1)
{
cat("Groups:\n")
@@ -1561,3 +1646,31 @@
invisible(x)
}
+##@shortBayesResult abbreviated sienaBayesFit results
+shortBayesResults <- function(x, nfirst=NULL){
+ if (!inherits(x, "sienaBayesFit"))
+ {
+ stop('x must be a sienaBayesFit object')
+ }
+ if (is.null(nfirst))
+ {
+ nfirst <- x$nwarm+1
+ }
+ df1 <- sienaFitThetaTable(x, fromBayes=TRUE, nfirst=nfirst)[[1]][,
+ c("text", "value", "se", "cFrom", "cTo", "postSd", "cSdFrom", "cSdTo" )]
+ df1$postSd[is.na(df1$cSdFrom)] <- NA
+ df1$postSd <- as.numeric(df1$postSd)
+ df1$cSdFrom <- as.numeric(df1$cSdFrom)
+ df1$cSdTo <- as.numeric(df1$cSdTo)
+ df2 <- as.data.frame(x$requestedEffects[,c("name","shortName", "interaction1", "interaction2",
+ "type", "randomEffects", "fix", "parm", "period", "effect1", "effect2", "effect3", "group")])
+ df2$period <- as.numeric(df2$period)
+ replace1 <- function(x){ifelse(x=="text", "effectName", x)}
+ replace2 <- function(x){ifelse(x=="value", "postMeanGlobal", x)}
+ replace3 <- function(x){ifelse(x=="se", "postSdGlobal", x)}
+ replace4 <- function(x){ifelse(x=="postSd", "postSdBetween", x)}
+ dfs <- cbind(df2, df1)
+ dfr <- dfs
+ names(dfr) <- replace1(replace2(replace3(replace4(names(dfs)))))
+ dfr
+}
Modified: pkg/RSiena/R/sienatable.r
===================================================================
--- pkg/RSiena/R/sienatable.r 2019-01-15 16:54:20 UTC (rev 338)
+++ pkg/RSiena/R/sienatable.r 2019-02-26 19:35:38 UTC (rev 339)
@@ -6,19 +6,46 @@
## * File: sienatable.r
## *
## * Description: This file contains the code to save a latex or html table of
-## * estimates for a sienaFit object
+## * estimates for a sienaFit or sienaBayesFit object
## * Written by Charlotte Greenan; small modifications by Tom Snijders.
## *
## ***************************************************************************/
##@siena.table siena07 Saves latex or html table of estimates
-## for a sienaFit object
+## for a sienaFit or sienaBayesFit object
siena.table <- function(x, type='tex',
file=paste(deparse(substitute(x)),'.',type,sep=""),
vertLine=TRUE, tstatPrint=FALSE,
- sig=FALSE, d=3)
+ sig=FALSE, d=3, nfirst=NULL)
{
+ objectName <- deparse(substitute(x))
+ fromBayes <- FALSE
+ xkind.string <- "sienaFit"
tstat <- tstatPrint
+ if (!inherits(x, "sienaFit"))
+ {
+ if (inherits(x, "sienaBayesFit"))
+ {
+ fromBayes <- TRUE
+ sig <- FALSE
+ tstat <- TRUE # for sienaBayesFit, the place for t-stats is used for between-groups sd.
+ if (is.null(nfirst))
+ {
+ nfirst <- x$nwarm + 1
+ cat("Note: the print function for sienaBayesFit objects")
+ cat(" can also use a parameter nfirst,\n")
+ cat(" indicating the first run")
+ cat(" from which convergence is assumed.\n")
+ cat(" The default value used now is nfirst =",
+ x$nwarm + 1, ".\n")
+ }
+ xkind.string <- "sienaBayesFit"
+ }
+ else
+ {
+ stop('x must be a sienaFit or sienaBayesFit object')
+ }
+ }
effects <- x$requestedEffects
p <- x$pp
condrates <- 0
@@ -29,6 +56,17 @@
condrates <- length(x$rate)
}
+ if (fromBayes)
+ {
+ xx <- shortBayesResults(x, nfirst=nfirst)
+ theta <- xx$postMeanGlobal
+ ses <- xx$postSdGlobal
+ sd.between <- xx$postSdBetween
+ se.string <- '(psd)'
+ sdb.string <- 'psd-between'
+ }
+ else
+ {
theta <- x$theta
theta[diag(x$covtheta) < 0.0 | x$fixed] <- NA
ses <- sqrt(diag(x$covtheta))
@@ -46,6 +84,8 @@
{
maxlincomb.t <- maxlincomb.t + 10^{-dd} #needs to be rounded up
}
+ se.string <- '(s.e.)'
+ }
if (length(x$condvarno) == 0)
{
condvarno <- 0
@@ -66,7 +106,7 @@
max.ses.width <- max.width(ses)
max.theta.width <- max.width(theta)
- max.tstat.width <- max.width(theta/ses)
+ max.tstat.width <- ifelse(fromBayes, max.width(sd.between), max.width(theta/ses))
## signif converts t values into daggers and asterisks
@@ -148,8 +188,15 @@
}
else
{
+ if (fromBayes)
+ {
+ tcsplit <- c("","")
+ }
+ else
+ {
tcsplit <- c("N","A.")
}
+ }
if (int.width>0)
{
@@ -214,7 +261,7 @@
{
if (type == "html")
{
- if (tstat == TRUE)
+ if (tstat)
{
amp5 <- rep("</TD><TD align=\"right\">",pp)
amp6 <- rep(".",pp)
@@ -291,10 +338,17 @@
if (type == "html")
{
- if (tstat == TRUE)
+ if (tstat)
{
+ if (fromBayes)
+ {
+ start.tstat <- "<TD>betw. sd.</TD>"
+ }
+ else
+ {
start.tstat <- "<TD>t stat.</TD>"
}
+ }
else
{
start.tstat <- ""
@@ -302,12 +356,23 @@
startTable <- tableSection(c("<TABLE border=1 rules=none frame = hsides>",
paste("<TR><TD>Effect</TD><TD> par.</TD><TD></TD>
- <TD> (s.e.)</TD>",start.tstat,"</TR>")))
+ <TD> ",se.string,"</TD>",start.tstat,"</TR>")))
midTable <- tableSection(c("",""))
indentTable <- tableSection("")
ruleTable <- tableSection("<TR> <TD colspan='9'><HR/></TD> </TR>")
footnoteStart <- c("</TABLE>","<TABLE border=1 rules=none frame = below>")
+ if (fromBayes)
+ {
+ footnote <- c(paste("<TR> <TD colspan=9 align=left> par=posterior mean;
+ psd = posterior standard deviation;
+ </TD> </TR> "),
+ paste("<TR> <TD colspan=9 align=left>
+ betw. sd = posterior between-groups stand. deviation.
+ </TD> </TR> <TR> </TR>"), "</TABLE>")
+ }
+ else
+ {
footnote <- c(paste(" <TR> <TD colspan=9 align=left>
all convergence t ratios < ",
max.t,".</TD> </TR> <TR> </TR>",
@@ -315,6 +380,7 @@
Overall maximum convergence ratio ",
maxlincomb.t,".</TD> </TR> <TR> </TR>",
sep="",collapse=""),"</TABLE>")
+ }
if (sig)
{
footnote <- c(footnoteStart, "<TR> <TD colspan=4 align=left> † p < 0.1;
@@ -328,25 +394,32 @@
}
else
{
- if (tstat == TRUE)
+ if (vertLine)
{
- start.tstat <- "r@{.}l"
- start.tstat2 <- "&\\multicolumn{2}{c}{$t$ stat.}"
+ linesep="|"
}
else
{
- start.tstat <- ""
- start.tstat2 <- ""
+ linesep=""
}
- if (vertLine)
+ if (tstat)
{
- linesep="|"
+ start.tstat <- "r@{.}l"
+ if (fromBayes)
+ {
+ start.tstat2 <- paste("&\\multicolumn{2}{c", linesep, "}{betw. sd}")
}
else
{
- linesep=""
+ start.tstat2 <- paste("&\\multicolumn{2}{c", linesep, "}{$t$ stat.}")
}
- startTable <- tableSection(c(paste("% Table based on sienaFit object",
+ }
+ else
+ {
+ start.tstat <- ""
+ start.tstat2 <- ""
+ }
+ startTable <- tableSection(c(paste("% Table based on", xkind.string, "object",
deparse(substitute(x)), ',', date()),
paste("\\begin{tabular}{l",
linesep,
@@ -355,14 +428,24 @@
"\\hline",
"\\rule{0pt}{2ex}\\relax",
paste("Effect &\\multicolumn{2}{c}{par.}&\\multicolumn{2}{c",
- linesep,
- "}{(s.e.)}",
+ "}{",se.string,"}",
start.tstat2,"\\\\[0.5ex]"),
"\\hline"))
midTable <- tableSection(c("\\hline",
"\\rule{0pt}{2ex}\\relax"))
indentTable <- tableSection("\\rule{0pt}{2ex}\\relax")
ruleTable <- tableSection("\\hline")
+ if (fromBayes)
+ {
+ footnote <- c(paste("\\multicolumn{7}{l}\n ",
+ "{\\footnotesize{par = posterior mean; psd = posterior standard deviation;}}\\\\\n",
+ "\\multicolumn{7}{l}\n ",
+ "{\\footnotesize{betw. sd = posterior between-groups stand. deviation.}}\\\\\n",
+ sep="",collapse=""),
+ "\\end{tabular}")
+ }
+ else
+ {
footnote <- c(paste("\\multicolumn{5}{l}\n ",
"{\\footnotesize{convergence $t$ ratios all $<$ ", max.t,
".}}\\\\\n",
@@ -379,6 +462,7 @@
$^{\\ast\\ast\\ast}$ $p$ $<$ 0.001;}}\\\\" ,footnote)
}
}
+ }
endTable <- tableSection(footnote)
@@ -441,8 +525,27 @@
remove <- unique(c(basicRates,fixed.2))
- if (tstat == TRUE)
+ if (tstat)
{
+ if (fromBayes)
+ {
+ remove <- is.na(sd.between)
+ mainTable$tstat1[-mid][-remove] <-
+ sapply(sd.between[rows],mystr,max.tstat.width)[1,][-remove]
+ mainTable$tstat2[-mid][-remove] <-
+ sapply(sd.between[rows],mystr)[2,][-remove]
+
+ if (any(remove))
+ {
+ if (type=='tex')
+ {
+ mainTable$tstat1[-mid][remove] <- "\\omit"
+ mainTable$tstat2[-mid][remove] <- "-"
+ }
+ }
+ }
+ else
+ {
mainTable$tstat1[-mid][-remove] <-
sapply(theta[rows]/ses[rows],mystr,max.tstat.width)[1,][-remove]
mainTable$tstat2[-mid][-remove] <-
@@ -458,7 +561,9 @@
}
}
- if (sig == TRUE)
+ }
+
+ if (sig)
{
mainTable$signif[-mid][-remove] <- sapply(theta[rows]/ses[rows],signif)[-remove]
}
@@ -474,7 +579,7 @@
rateTable$se1 <- sapply(x$vrate,mystr,max.vrate.width)[1,]
rateTable$se2 <- sapply(x$vrate,mystr)[2,]
- if (tstat==TRUE)
+ if (tstat)
{
if (type=='tex')
{
Modified: pkg/RSiena/data/allEffects.csv
===================================================================
--- pkg/RSiena/data/allEffects.csv 2019-01-15 16:54:20 UTC (rev 338)
+++ pkg/RSiena/data/allEffects.csv 2019-02-26 19:35:38 UTC (rev 339)
@@ -3,7 +3,8 @@
behaviorOneModeObjective,xxxxxx total similarity,xxxxxx total similarity,totSim,TRUE,yyyyyy,,eval,FALSE,FALSE,FALSE,FALSE,FALSE,",",0,0,objective,NA,NA,0,0,0,0,,FALSE
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 339
More information about the Rsiena-commits
mailing list