[Rsiena-commits] r250 - in pkg: RSiena RSiena/R RSiena/inst/doc RSiena/man RSiena/tests RSienaTest RSienaTest/R RSienaTest/inst/doc RSienaTest/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 4 17:15:34 CET 2013
Author: tomsnijders
Date: 2013-12-04 17:15:33 +0100 (Wed, 04 Dec 2013)
New Revision: 250
Added:
pkg/RSiena/man/Wald.Rd
pkg/RSienaTest/man/Wald.Rd
Modified:
pkg/RSiena/DESCRIPTION
pkg/RSiena/NAMESPACE
pkg/RSiena/R/Sienatest.r
pkg/RSiena/R/initializeFRAN.r
pkg/RSiena/R/phase2.r
pkg/RSiena/R/phase3.r
pkg/RSiena/R/print01Report.r
pkg/RSiena/R/print07Report.r
pkg/RSiena/R/printInitialDescription.r
pkg/RSiena/R/robmon.r
pkg/RSiena/R/siena08.r
pkg/RSiena/R/sienaDataCreate.r
pkg/RSiena/R/sienaGOF.r
pkg/RSiena/R/sienaeffects.r
pkg/RSiena/R/sienaprint.r
pkg/RSiena/R/sienatable.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/coCovar.Rd
pkg/RSiena/man/getEffects.Rd
pkg/RSiena/man/includeInteraction.Rd
pkg/RSiena/man/setEffect.Rd
pkg/RSiena/man/sienaAlgorithmCreate.Rd
pkg/RSiena/man/sienaCompositionChange.Rd
pkg/RSiena/man/sienaGOF-auxiliary.Rd
pkg/RSiena/man/sienaGOF.Rd
pkg/RSiena/man/varCovar.Rd
pkg/RSiena/tests/parallel.Rout.save
pkg/RSienaTest/DESCRIPTION
pkg/RSienaTest/NAMESPACE
pkg/RSienaTest/R/Sienatest.r
pkg/RSienaTest/R/initializeFRAN.r
pkg/RSienaTest/R/phase3.r
pkg/RSienaTest/R/print01Report.r
pkg/RSienaTest/R/sienaDataCreate.r
pkg/RSienaTest/R/sienaGOF.r
pkg/RSienaTest/R/sienaeffects.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/inst/doc/effects.pdf
pkg/RSienaTest/man/RSiena-package.Rd
pkg/RSienaTest/man/coCovar.Rd
pkg/RSienaTest/man/includeInteraction.Rd
pkg/RSienaTest/man/setEffect.Rd
pkg/RSienaTest/man/varCovar.Rd
Log:
Small changes: Wald tests; setEffect and siena07(prevAns) now can handle interaction effects better; actor covariates do not need to be centered; and some other small things.
Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/DESCRIPTION 2013-12-04 16:15:33 UTC (rev 250)
@@ -1,8 +1,8 @@
Package: RSiena
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.1-249
-Date: 2013-10-31
+Version: 1.1-250
+Date: 2013-12-04
Author: Ruth Ripley, Krists Boitmanis, Tom A.B. Snijders
Depends: R (>= 2.15.0)
Imports: Matrix
Modified: pkg/RSiena/NAMESPACE
===================================================================
--- pkg/RSiena/NAMESPACE 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/NAMESPACE 2013-12-04 16:15:33 UTC (rev 250)
@@ -7,9 +7,11 @@
varCovar, varDyadCovar, setEffect, includeEffects, includeInteraction,
effectsDocumentation, sienaDataConstraint, print.xtable.sienaFit,
installGui, siena08, iwlsm, sienaTimeTest, includeTimeDummy,
- sienaGOF, sparseMatrixExtraction, networkExtraction, behaviorExtraction,
+ sienaGOF, descriptives.sienaGOF, sparseMatrixExtraction,
+ networkExtraction, behaviorExtraction,
OutdegreeDistribution, IndegreeDistribution, BehaviorDistribution,
- siena.table, xtable)
+ siena.table, xtable,
+ Wald.RSiena, Multipar.RSiena)
import(Matrix)
Modified: pkg/RSiena/R/Sienatest.r
===================================================================
--- pkg/RSiena/R/Sienatest.r 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/R/Sienatest.r 2013-12-04 16:15:33 UTC (rev 250)
@@ -260,7 +260,11 @@
Report('Error message for inversion of v9: \n', cf)
vav <- v9
vav[] <- NA
+ cvalue <- NA
+ oneSided <- NA
}
+ else
+ {
cvalue <- t(ov) %*% vav %*% ov
if (cvalue < 0) cvalue <- 0
if (sum(test) == 1)
@@ -277,5 +281,40 @@
oneSided <- 0
}
}
+ }
list(cvalue=cvalue, oneSided=oneSided, covMatrix=v9)
}
+
+##@Wald.RSiena Calculate Wald test statistics
+Wald.RSiena <- function(A, ans)
+{
+ if (is.vector(A)){A <- matrix(A, nrow=1)}
+ if (dim(A)[2] != ans$pp){stop(paste('A must have', ans$pp, 'columns.'))}
+ sigma <- ans$covtheta
+ if (any(is.na(sigma))){ # happens when some coordinates were fixed
+ # in the call of siena07 leading to ans;
+ # then the non-used part of sigma,
+ # which partially consists of NA,
+ # is replaced by the identity matrix.
+ zero.cols <- apply(A, 2, function(colum){all(colum==0)})
+ sigma[zero.cols, ] <- 0
+ sigma[, zero.cols] <- 0
+ diag(sigma)[zero.cols] <- 1
+ }
+ 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)
+}
+
+##@Multipar.RSiena Calculate Wald test statistic for hypothesis that subvector = 0.
+Multipar.RSiena <- function(ans, ...)
+{
+ p <- length(ans$theta)
+ k <- length(c(...))
+ A <- matrix(0, nrow=k, ncol=p)
+ A[cbind(1:k,c(...))] <- 1
+ Wald.RSiena(A, ans)
+}
\ No newline at end of file
Modified: pkg/RSiena/R/initializeFRAN.r
===================================================================
--- pkg/RSiena/R/initializeFRAN.r 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/R/initializeFRAN.r 2013-12-04 16:15:33 UTC (rev 250)
@@ -266,8 +266,8 @@
## z$FinDiffBecauseSymmetric <- TRUE
if (x$maxlike)
{
- stop("Maximum likelihood method not implemented",
- "for symmetric networks")
+# stop("Maximum likelihood method not implemented",
+# "for symmetric networks")
}
if (x$modelType == 1)
{
@@ -1971,13 +1971,13 @@
paste(x[c("name", "shortName",
"type", "groupName",
"interaction1", "interaction2",
- "period")],
+ "period", "effect1", "effect2", "effect3")],
collapse="|"))
efflist <- apply(effects, 1, function(x)
paste(x[c("name", "shortName",
"type", "groupName",
"interaction1", "interaction2",
- "period")],
+ "period", "effect1", "effect2", "effect3")],
collapse="|"))
use <- efflist %in% oldlist
effects$initialValue[use] <-
Modified: pkg/RSiena/R/phase2.r
===================================================================
--- pkg/RSiena/R/phase2.r 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/R/phase2.r 2013-12-04 16:15:33 UTC (rev 250)
@@ -323,12 +323,13 @@
}
}
## limit change. Reporting is delayed to end of phase.
-# The old version had a numerical ifelse condition:
-# maxrat<- max(ifelse(z$sd, abs(fra)/ z$sd, 1.0))
-# But ifelse with a numerical first argument only tests the value 0,
-# which here has probability 0... So this really has no effect.
-# I (TS) do not understand this.
-# Therefore I changed it.
+# The truncation has been different from version 1.1-227 to 1.1-243,
+# due to a misunderstanding.
+# In version 1.1-244 it was changed back to the old Siena 3 way.
+# Except now the threshold is 5 instead of 10.
+# For the case !x$diagg it might be better
+# to base truncation on some multivariate norm of z$dfra.
+ maxRatio <- max(ifelse(z$fixed, 1.0, abs(fra)/ z$sd))
if (x$diagg)
{
changestep <- fra / diag(z$dfra)
@@ -338,10 +339,9 @@
changestep <- as.vector(fra %*% z$dinvv)
}
changestep[z$fixed] <- 0.0
- maxratt <- max(abs(changestep))
- if (maxratt > 5)
+ if (maxRatio > 5)
{
- changestep <- 5*changestep/maxratt
+ changestep <- 5*changestep/maxRatio
z$truncated[z$nit] <- TRUE
}
fchange <- as.vector(z$gain * changestep)
Modified: pkg/RSiena/R/phase3.r
===================================================================
--- pkg/RSiena/R/phase3.r 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/R/phase3.r 2013-12-04 16:15:33 UTC (rev 250)
@@ -156,7 +156,7 @@
tstat <- rep(NA, z$pp)
tstat[!use]<- sf[!use] / sqrt(dmsf[!use])
tstat[use & use2] <- 0
- tstat[use & !use2] <- 999
+ tstat[use & !use2] <- NA # was 999
z$tstat <- tstat
# tconv.max = Maximum value of t-ratio for convergence,
# for any linear combination.
@@ -267,7 +267,7 @@
if (inherits(try(cov <- solve(dfrac), silent=TRUE),"try-error"))
{
Report('Noninvertible estimated covariance matrix : \n', outf)
- errorMessage.cov <- '***Warning: Noninvertible estimated covariance matrix ***'
+ errorMessage.cov <- '***Warning: linear dependencies between statistics ***'
cov <- NULL
}
}
@@ -300,11 +300,13 @@
{
z$diver <- (z$fixed | z$diver | diag(cov) < 1e-9) & (!z$AllUserFixed)
## beware: recycling works for one direction but not the other
- diag(cov)[z$diver] <- 99 * 99
- cov[z$diver, ] <- rep(Root(diag(cov)), each=sum(z$diver)) * 33
- diag(cov)[z$diver] <- 99 * 99
- cov[, z$diver] <- rep(Root(diag(cov)), sum(z$diver)) * 33
- diag(cov)[z$diver] <- 99 * 99
+ diag(cov)[z$diver] <- NA
+# cov[z$diver, ] <- rep(Root(diag(cov)), each=sum(z$diver)) * 33
+# diag(cov)[z$diver] <- 99 * 99
+# cov[, z$diver] <- rep(Root(diag(cov)), sum(z$diver)) * 33
+ cov[z$diver, ] <- NA
+ cov[, z$diver] <- NA
+ diag(cov)[z$diver] <- NA
}
z$covtheta <- cov
}
@@ -396,7 +398,7 @@
z$errorMessage.cov <- '*** Warning: Noninvertible derivative matrix ***'
fchange <- 0
z$dinv <- z$dfrac
- z$dinv[,] <- 999
+ z$dinv[,] <- NA # 999
}
else
{
Modified: pkg/RSiena/R/print01Report.r
===================================================================
--- pkg/RSiena/R/print01Report.r 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/R/print01Report.r 2013-12-04 16:15:33 UTC (rev 250)
@@ -511,25 +511,50 @@
width=3, nsmall=1), '%)\n'), outf)
}
Report("\nInformation about covariates:\n", outf)
- Report(c(format("minimum maximum mean", width=38,
+ 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$cCovars[[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)
}
if (nData <= 1)
{
+ if (any.noncent <= 0)
+ {
Report(c("The mean value", ifelse(nCovars == 1, " is", "s are"),
- " subtracted from the covariate",
+ " subtracted from the",
+ ifelse(nCovars == 1, " centered", ""), " covariate",
ifelse(nCovars == 1, ".\n\n", "s.\n\n")), sep="", outf)
}
+ else if (any.cent >= 1)
+ {
+ s.plural <- ""
+ if (any.cent >= 2){s.plural <- "s"}
+ 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)
+ }
+ }
}
##@reportChangingCovariates internal print01Report
reportChangingCovariates <- function()
@@ -600,13 +625,25 @@
}
}
Report("\nInformation about changing covariates:\n", outf)
- Report(c(format("minimum maximum mean", width=38,
+ 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))
{
if (use[i])
{
atts <- attributes(x$vCovars[[i]])
+ if (atts$centered)
+ {
+ cent <- " Y"
+ any.cent <- any.cent+1
+ }
+ else
+ {
+ cent <- " N"
+ any.noncent <- any.noncent+1
+ }
Report(c(covars[i], '\n'), outf) # name
for (j in 1:(ncol(x$vCovars[[i]])))
{
@@ -621,15 +658,27 @@
}
Report(c(format("Overall", width=29),
format(round(atts$mean, 3), width=10, nsmall=3),
- "\n"), outf)
-
+ cent, "\n"), outf)
}
}
if (nData <= 1)
{
- Report("\nThe overall mean value", outf)
- Report(c(ifelse(nCovars == 1, " is", "s are"),
- "subtracted from the covariate.\n\n"), 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)
+ {
+ s.plural <- ""
+ if (any.cent >= 2){s.plural <- "s"}
+ 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)
+ }
}
}
##@reportConstantDyadicCovariates internal print01Report
@@ -1061,7 +1110,6 @@
netnames[i], ".\n"), sep = "", outf)
}
}
- Report("\n", outf)
if (sum(atts$types == 'oneMode') > 0)
{
Modified: pkg/RSiena/R/print07Report.r
===================================================================
--- pkg/RSiena/R/print07Report.r 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/R/print07Report.r 2013-12-04 16:15:33 UTC (rev 250)
@@ -176,6 +176,7 @@
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)
+# sf and cov.dev may be dropped - just for now (07-10-13) I keep them in
# sf <- colMeans(z$sf)
mean.stats <- colMeans(z$sf) + z$targets
# cov.dev <- z$msf
@@ -196,7 +197,7 @@
PrtOutMat(as.matrix(mymess1), bof)
if (x$dolby)
{
- Report('Standard errors of the mean are less than s.d./n \n', outf)
+ 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/printInitialDescription.r
===================================================================
--- pkg/RSiena/R/printInitialDescription.r 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/R/printInitialDescription.r 2013-12-04 16:15:33 UTC (rev 250)
@@ -182,21 +182,21 @@
Report(c(format(per + periodFromStart, width=3),
" ==> ",
- format(per + 1 + periodFromStart, width=2),
+ format(per + 1 + periodFromStart, width=3),
format(matchange[, per], width=10),
format(attr(depvar, "distance")[per],
width=10),
- jaccard, format(misd, width=6), " (",
+ jaccard, format(misd, width=10), " (",
round(100 * misd/(ntot + misd)), "%)\n"),
sep="", outf)
}
else
{
Report(c(per + periodFromStart, " ==> ",
- format(per + 1 + periodFromStart, width=2),
+ format(per + 1 + periodFromStart, width=3),
format(matchange[, per], width=10),
format(attr(depvar,"distance")[per], width=10),
- format(misd, width=7), " (",
+ format(misd, width=10), " (",
round(100 * misd/(ntot + misd)), "%)\n"),
sep="", outf)
Modified: pkg/RSiena/R/robmon.r
===================================================================
--- pkg/RSiena/R/robmon.r 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/R/robmon.r 2013-12-04 16:15:33 UTC (rev 250)
@@ -6,7 +6,7 @@
# * File: robmon.r
# *
# * Description: This module contains the function robmon which controls
-# * the phases of the robbins-munro stochastic approximation algorithm.
+# * the phases of the robbins-monro stochastic approximation algorithm.
# *****************************************************************************/
##args:x: model object - intended to be read only
## z: model fitting object
@@ -333,10 +333,10 @@
if (!x$simOnly)
{
z$diver<- (z$fixed | z$diver | diag(z$covtheta) < 1e-9) & (!z$AllUserFixed)
- z$covtheta[z$diver, ] <- Root(diag(z$covtheta)) * 33
+ z$covtheta[z$diver, ] <- NA # was Root(diag(z$covtheta)) * 33
##not sure this does not use very small vals
- z$covtheta[, z$diver] <- Root(diag(z$covtheta)) * 33
- diag(z$covtheta)[z$diver] <- 999
+ z$covtheta[, z$diver] <- NA # was Root(diag(z$covtheta)) * 33
+ diag(z$covtheta)[z$diver] <- NA # was 999
}
z$termination <- 'OK'
z
Modified: pkg/RSiena/R/siena08.r
===================================================================
--- pkg/RSiena/R/siena08.r 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/R/siena08.r 2013-12-04 16:15:33 UTC (rev 250)
@@ -670,11 +670,10 @@
width=9),
" (d.f. = ", 2 * y$ns, "), ",
reportp(y$scoreplusp, 3), "\n"), sep="", outf)
- Report("Combination of left one-sided p-values:\n", outf)
- Report(c("Chi-squared = ",
- format(round(y$scoreminus, 4), width=9),
- " (d.f. = ", 2 * y$ns, "), ",
- reportp(y$scoreminusp, 3), "\n"), sep="", outf)
+ Report("Bonferroni combination of left and right one-sided p-values:\n",
+ outf)
+ Report(c(reportp(2*min(y$scoreminusp, y$scoreplusp), 3), "\n"),
+ sep="", outf)
}
}, y=x))
}
Modified: pkg/RSiena/R/sienaDataCreate.r
===================================================================
--- pkg/RSiena/R/sienaDataCreate.r 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/R/sienaDataCreate.r 2013-12-04 16:15:33 UTC (rev 250)
@@ -14,18 +14,30 @@
##@addAttributes.coCovar DataCreate
addAttributes.coCovar <- function(x, name, ...)
{
+ storage.mode(x) <- 'double'
varmean <- mean(x, na.rm=TRUE)
range2 <- range(x, na.rm=TRUE)
attr(x, 'moreThan2') <- length(table(x)) > 2
vartotal <- sum(x, na.rm=TRUE)
nonMissingCount <- sum(!is.na(x))
+ if (attr(x, "centered"))
+ {
x <- x - varmean
+ }
+ else
+ {
+ x <- x - 0.0
+ }
attr(x, 'mean') <- varmean
rr <- rangeAndSimilarity(x, range2)
if (rr$range[2] == rr$range[1] && !any(is.na(x)))
+ {
attr(x, 'poszvar') <- FALSE
+ }
else
+ {
attr(x, 'poszvar') <- TRUE
+ }
attr(x, 'range') <- rr$range[2] - rr$range[1]
storage.mode(attr(x, 'range')) <- 'double'
attr(x, 'range2') <- range2
@@ -51,12 +63,23 @@
attr(x, 'range') <- cr[2] - cr[1]
storage.mode(attr(x, 'range')) <- 'double'
attr(x, 'mean') <- varmean
+ if (attr(x, "centered"))
+ {
x <- x - varmean
+ }
+ else
+ {
+ x <- x - 0.0
+ }
rr <- rangeAndSimilarity(tmpmat, cr)
if (rr$range[2] == rr$range[1] && !any(is.na(tmpmat)))
+ {
attr(x, 'poszvar') <- FALSE
+ }
else
+ {
attr(x, 'poszvar') <- TRUE
+ }
attr(x, 'simMean') <- rr$simMean
attr(x, 'moreThan2') <- length(unique(x)) > 2
attr(x, 'name') <- name
@@ -1193,6 +1216,10 @@
attr(group[[i]]$vCovars[[j]], "nonMissingCount")
}
varmean <- vartotal / nonMissingCount
+#browser() # Hier kijken hoe je moet centreren in de groep.
+ j <- match(atts$vCovars[covar], names(group[[1]]$vCovars))
+ if (attr(group[[1]]$vCovars[[j]],"centered"))
+ {
for (i in 1:length(group))
{
j <- match(atts$vCovars[covar], names(group[[i]]$vCovars))
@@ -1204,6 +1231,7 @@
group[[i]]$vCovars[[j]] <- group[[i]]$vCovars[[j]] -
varmean
}
+ }
simTotal <- 0
simCnt <- 0
anyMissing <- FALSE
@@ -1356,6 +1384,7 @@
attr(x, "meanp") <- rep(atts$mean, ncol(x))
attr(x, "range") <- atts$range
attr(x, 'mean') <- atts$mean
+ attr(x, 'centered') <- atts$centered
attr(x, 'vartotal') <- atts$vartotal
attr(x, 'nonMissingCount') <- atts$nonMissingCount
attr(x, 'simMeans') <- atts$simMeans
Modified: pkg/RSiena/R/sienaGOF.r
===================================================================
--- pkg/RSiena/R/sienaGOF.r 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/R/sienaGOF.r 2013-12-04 16:15:33 UTC (rev 250)
@@ -465,6 +465,7 @@
originalMhd[j],") for current model.\n")
}
}
+ invisible(x)
}
##@summary.sienaGOF siena07 Summary method for sienaGOF
@@ -516,6 +517,7 @@
}
cat("\nComputation time for auxiliary statistic calculations on simulations: ",
attr(x, "simTime")["elapsed"] , "seconds.\n")
+ invisible(x)
}
##@plot.sienaGOF siena07 Plot method for sienaGOF
@@ -636,11 +638,12 @@
}
else
{
- key <- attr(x,"key")
+ key <- attr(x,"key")[screen]
}
}
else
{
+ key <- key[screen] ## added 1.1-244
if (length(key) != ncol(obs))
{
stop("Key length does not match the number of variates.")
@@ -689,6 +692,109 @@
}
+##@descriptives.sienaGOF siena07 Gives numerical values in the plot.
+descriptives.sienaGOF <- function (x, center=FALSE, scale=FALSE,
+ perc=.05, key=NULL, period=1)
+{
+# adapted excerpt from plot.sienaGOF
+ if (attr(x,"joined"))
+ {
+ x <- x[[1]]
+ }
+ else
+ {
+ x <- x[[period]]
+ }
+
+ sims <- x$Simulations
+ obs <- x$Observations
+ itns <- nrow(sims)
+
+ screen <- sapply(1:ncol(obs),function(i){
+ (sum(is.nan(rbind(sims,obs)[,i])) == 0) }) &
+ (diag(var(rbind(sims,obs)))!=0)
+ sims <- sims[,screen, drop=FALSE]
+ obs <- obs[,screen, drop=FALSE]
+ ## Need to check for useless statistics here:
+ n.obs <- nrow(obs)
+
+ if (is.null(key))
+ {
+ if (is.null(attr(x, "key")))
+ {
+ key=(1:sum(screen))
+ }
+ else
+ {
+ key <- attr(x,"key")[screen]
+ }
+ }
+ else
+ {
+ if (length(key) != ncol(obs))
+ {
+ stop("Key length does not match the number of variates.")
+ }
+ key <- key[screen]
+ }
+
+ sims.themin <- apply(sims, 2, min)
+ sims.themax <- apply(sims, 2, max)
+ sims.min <- pmin(sims.themin, obs)
+ sims.max <- pmax(sims.themax, obs)
+
+ if (center)
+ {
+ sims.median <- apply(sims, 2, median)
+ sims <- sapply(1:ncol(sims), function(i)
+ (sims[,i] - sims.median[i]) )
+ obs <- matrix(sapply(1:ncol(sims), function(i)
+ (obs[,i] - sims.median[i])), nrow=n.obs )
+ sims.min <- sims.min - sims.median
+ sims.max <- sims.max - sims.median
+ }
+
+ if (scale)
+ {
+ sims.range <- sims.max - sims.min + 1e-6
+ sims <- sapply(1:ncol(sims), function(i) sims[,i]/(sims.range[i]))
+ obs <- matrix(sapply(1:ncol(sims), function(i) obs[,i]/(sims.range[i]))
+ , nrow=n.obs )
+ sims.min <- sims.min/sims.range
+ sims.max <- sims.max/sims.range
+ }
+
+# ymin <- 1.05*min(sims.min) - 0.05*max(sims.max)
+# ymax <- -0.05*min(sims.min) + 1.05*max(sims.max)
+ screen <- sapply(1:ncol(obs),function(i){
+ (sum(is.nan(rbind(sims,obs)[,i])) == 0) }) &
+ (diag(var(rbind(sims,obs)))!=0)
+ sims <- sims[,screen, drop=FALSE]
+ obs <- obs[,screen, drop=FALSE]
+ sims.themin <- sims.themin[screen, drop=FALSE]
+ sims.themax <- sims.themax[screen, drop=FALSE]
+
+ ind.lower = max( round(itns * perc/2), 1)
+ ind.upper = round(itns * (1-perc/2))
+ ind.median = round(itns * 0.5)
+ yperc.mid = sapply(1:ncol(sims), function(i)
+ sort(sims[,i])[ind.median])
+ yperc.lower = sapply(1:ncol(sims), function(i)
+ sort(sims[,i])[ind.lower] )
+ yperc.upper = sapply(1:ncol(sims), function(i)
+ sort(sims[,i])[ind.upper] )
+ violins <- matrix(NA, 6, ncol(sims))
+ violins[1,] <- sims.themax
+ violins[2,] <- yperc.upper
+ violins[3,] <- yperc.mid
+ violins[4,] <- yperc.lower
+ violins[5,] <- sims.themin
+ violins[6,] <- obs
+ rownames(violins) <- c('max', 'perc.upper', 'median', 'perc.lower', 'min', 'obs')
+ colnames(violins) <- key
+ violins
+}
+
##@changeToStructural sienaGOF Utility to change
# values in X to structural values in S
# X must have values 0, 1.
@@ -840,6 +946,8 @@
dimsOfDepVar<- attr(obsData[[groupName]]$depvars[[varName]], "netdims")
isbipartite <- (attr(obsData[[groupName]]$depvars[[varName]], "type")
=="bipartite")
+# sparseData may be dropped - if that's OK
+# sparseData <- (attr(obsData[[groupName]]$depvars[[varName]], "sparse"))
# For bipartite networks in package <network>,
# the number of nodes is equal to
# the number of actors (rows) plus the number of events (columns)
Modified: pkg/RSiena/R/sienaeffects.r
===================================================================
--- pkg/RSiena/R/sienaeffects.r 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/R/sienaeffects.r 2013-12-04 16:15:33 UTC (rev 250)
@@ -46,8 +46,9 @@
{
cat(paste("There is no effect with short name "))
cat(paste(effectNames,", \n", sep=""))
- cat(paste("and with interaction1 = <",interaction1,">", sep=""))
- cat(paste(" and interaction2 = <",interaction2,">,\n", sep=""))
+ cat(paste("and with interaction1 = <",interaction1,">, ", sep=""))
+ cat(paste("interaction2 = <",interaction2,">, ", sep=""))
+ cat(paste("and type = <",type,">, \n", sep=""))
cat(paste("for dependent variable",name,".\n"))
}
else
@@ -61,9 +62,8 @@
##@includeInteraction DataCreate
includeInteraction <- function(myeff, ...,
include=TRUE, name=myeff$name[1],
- type="eval", interaction1=rep("", 3),
- interaction2=rep("", 3), character=FALSE,
- verbose=TRUE)
+ type="eval", interaction1=rep("", 3), interaction2=rep("", 3),
+ character=FALSE, verbose=TRUE)
{
if (character)
{
@@ -222,6 +222,7 @@
timeDummy=",",
include=TRUE, name=myeff$name[1],
type="eval", interaction1="", interaction2="",
+ effect1=0, effect2=0, effect3=0,
period=1, group=1, character=FALSE)
{
if (!character)
@@ -235,8 +236,25 @@
myeff$interaction2 == interaction2 &
(is.na(myeff$period) | myeff$period == period) &
myeff$group == group
+ if (shortName == "unspInt")
+ {
+ use <- use & (myeff$include) & (myeff$effect1 == effect1) &
+ (myeff$effect2 == effect2) & (myeff$effect3 == effect3)
+ }
if (sum(use) == 0)
{
+ cat(paste("There is no effect with short name "))
+ cat(paste(shortName,", \n", sep=""))
+ cat(paste("and with interaction1 = <",interaction1,">, ", sep=""))
+ cat(paste("interaction2 = <",interaction2,">, ", sep=""))
+ cat(paste("type = <",type,">, ", sep=""))
+ cat(paste("period = <",period,">, ", sep=""))
+ if (shortName == "unspInt")
+ {
+ cat(paste("effects1-2-3 = <",effect1, effect2, effect3,">,", sep=" "))
+ }
+ cat(paste("and group = <",group,">, \n ", sep=""))
+ cat(paste("for dependent variable",name,".\n"))
stop("Effect not found")
}
if (sum(use) > 1)
Modified: pkg/RSiena/R/sienaprint.r
===================================================================
--- pkg/RSiena/R/sienaprint.r 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/R/sienaprint.r 2013-12-04 16:15:33 UTC (rev 250)
@@ -48,6 +48,10 @@
}
cat('Dependent variables: ', paste(names(x$depvars), collapse=", "), "\n")
cat('Number of observations:', x$observations, "\n\n")
+ if (length(x$compositionChange) > 0)
+ {
+ cat('With composition change.\n')
+ }
if (!is.null(x$nodeSet))
{
tmp <- sapply(x$nodeSet, length)
@@ -69,6 +73,10 @@
mymat <- matrix("", 7, 2)
mymat[1,] <- c("Dependent variable", attr(xj, "name"))
mymat[2,] <- c("Type", attr(xj, "type"))
+ if (attr(xj,"symmetric"))
+ {
+ mymat[2,2] <- c(mymat[2,2],", symmetric")
+ }
mymat[3,] <- c("Observations", attr(xj, "netdims")[3])
if (attr(xj, "type") == "bipartite")
{
@@ -268,6 +276,7 @@
cat('\nEstimated means and standard deviations, standard errors of the mean \n')
dmsf <- diag(x$msf)
mean.stats <- colMeans(x$sf) + x$targets
+# cov.dev may be droppec - just for now (07-10-13) I keep it in
# cov.dev <- x$msf
sem <- sqrt(dmsf/dim(x$sf)[1])
if (x$x$dolby)
@@ -292,7 +301,7 @@
cat("\nTotal of", x$n, "iteration steps.\n\n")
if ((x$x$dolby)&(x$x$simOnly))
{
- cat('(Standard errors of the mean are less than s.d./', x$n,' \n',
+ 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')
}
@@ -323,6 +332,16 @@
print.sienaFit(x)
cat(c("Overall maximum convergence ratio: ",
sprintf("%8.4f", x$tconv.max), "\n\n"))
+
+ if (x$maxlike)
+ {
+ cat('Autocorrelations during phase 3 : \n')
+ cat(paste(format(1:x$pp, width=3), '. ',
+ format(x$sfl, width=8, digits=4),
+ '\n', collapse="", sep=""))
+ cat ('\n')
+ }
+
if (sum(x$test) > 0) ## we have some score tests
{
testn <- sum(x$test)
@@ -440,10 +459,12 @@
##@print.sienaAlgorithm Methods
print.sienaAlgorithm <- function(x, ...)
{
+ cat(' Siena Algorithm specification.\n')
cat(' Project name:', x$projname, '\n')
cat(' Use standard initial values:', x$useStdInits, '\n')
cat(' Random seed:', x$randomSeed,'\n')
cat(' Starting value of gain parameter:', x$firstg, '\n')
+ cat(' Reduction factor for gain parameter:', x$reduceg, '\n')
cat(' Diagonalization parameter:', x$diagonalize, '\n')
cat(' Dolby noise reduction:', x$dolby, '\n')
if (any(x$MaxDegree > 0))
Modified: pkg/RSiena/R/sienatable.r
===================================================================
--- pkg/RSiena/R/sienatable.r 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/R/sienatable.r 2013-12-04 16:15:33 UTC (rev 250)
@@ -1,4 +1,4 @@
-## /*****************************************************************************
+## /***************************************************************************
## * SIENA: Simulation Investigation for Empirical Network Analysis
## *
## * Web: http://www.stats.ox.ac.uk/~snijders/siena
@@ -7,10 +7,12 @@
## *
## * Description: This file contains the code to save a latex or html table of
## * estimates for a sienaFit object
+## * Written by Charlotte Greenan.
## *
-## ****************************************************************************/
+## ***************************************************************************/
-##@siena.table siena07 Saves latex or html table of estimates for a sienaFit object
+##@siena.table siena07 Saves latex or html table of estimates
+## for a sienaFit object
siena.table <- function(x,type='tex',
file=paste(deparse(substitute(x)),'.',type,sep=""),
vertLine=TRUE,tstatPrint=FALSE,
@@ -37,12 +39,12 @@
{
max.t <- max.t + 10^{-d} #needs to be rounded up
}
- maxlincomb.t1 <- x$tconv.max
- maxlincomb.t <- round(maxlincomb.t1, digits = d)
- if (maxlincomb.t < maxlincomb.t1)
- {
- maxlincomb.t <- maxlincomb.t + 10^{-d} #needs to be rounded up
- }
+# maxlincomb.t1 <- x$tconv.max
+# maxlincomb.t <- round(maxlincomb.t1, digits = d)
+# if (maxlincomb.t < maxlincomb.t1)
+# {
+# maxlincomb.t <- maxlincomb.t + 10^{-d} #needs to be rounded up
+# }
if (length(x$condvarno) == 0)
{
condvarno <- 0
@@ -294,7 +296,7 @@
linesep=""
}
startTable <- tableSection(c(paste("% Table based on sienaFit object",
- deparse(substitute(x))),
+ deparse(substitute(x)), ',', date()),
paste("\\begin{tabular}{l",
linesep,
"r@{.}l r@{.}l",start.tstat,
@@ -312,9 +314,11 @@
ruleTable <- tableSection("\\hline")
footnote <- c(paste("\\multicolumn{5}{l}\n ",
"{\\footnotesize{convergence $t$ ratios all $<$ ", max.t,
- ",}}\\\\\n", "\\multicolumn{5}{l}",
- "{\\footnotesize{overall maximum convergence ratio ",
- maxlincomb.t,".}}", sep="",collapse=""),
+ ".}}\\\\\n",
+# "\\multicolumn{5}{l}",
+# "{\\footnotesize{Overall maximum convergence ratio ",
+# maxlincomb.t,".}}",
+ sep="",collapse=""),
"\\end{tabular}")
if (sig == TRUE)
Modified: pkg/RSiena/R/sienautils.r
===================================================================
--- pkg/RSiena/R/sienautils.r 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/R/sienautils.r 2013-12-04 16:15:33 UTC (rev 250)
@@ -122,7 +122,7 @@
}
##@coCovar Create
-coCovar <- function(val, nodeSet="Actors")
+coCovar <- function(val, centered=TRUE, nodeSet="Actors")
{
##vector, numeric or factor
if (!is.vector(val))
@@ -135,11 +135,12 @@
}
out <- val
class(out) <- "coCovar"
+ attr(out, "centered") <- centered
attr(out, "nodeSet") <- nodeSet
out
}
##@varCovar Create
-varCovar<- function(val, nodeSet="Actors")
+varCovar<- function(val, centered=TRUE, nodeSet="Actors")
{
##matrix, numeric or factor, nrow = nactors and cols = observations-1
if (!is.matrix(val))
@@ -152,6 +153,7 @@
}
out <- val
class(out) <- "varCovar"
+ attr(out, "centered") <- centered
attr(out, "nodeSet") <- nodeSet
out
}
Modified: pkg/RSiena/changeLog
===================================================================
--- pkg/RSiena/changeLog 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/changeLog 2013-12-04 16:15:33 UTC (rev 250)
@@ -1,3 +1,26 @@
+2013-12-04 R-Forge Revision 250
+Changes in RSiena and RSienaTest:
+ * New option "centered" in coCovar and varCovar (sienautils.r,
+ sienaDataCreate.r, print01Report.r, sienaprint.r).
+ * setEffect, updateTheta, and prevAns in siena07() now also cater
+ for unspecified interactions (sienaeffects.r, initializeFRAN.r).
+ * Wald.RSiena and Multipar.RSiena added (Sienatest.r).
+ * Error occurrence with message about cvalue in EvaluateTestStatistic
+ corrected (Sienatest.r).
+ * Divergent parameters in siena07() get NA for their rows and columns
+ in the covariance matrix (function phase3.2 in phase3.r).
+The following changes in revision 244 were ported from RSienaTest to RSiena:
+ * In siena08, also report Bonferroni combination
+ of the two Fisher combinations.
+ * In phase2.r, rolled back change in truncation from version 1.1-227
+ to the earlier procedure.
+ * descriptives.sienaGOF added.
+ * Minor changes of output in siena.table, print07Report.r, print.siena,
+ printInitialDescription.r, and in error message for includeEffects.
+ * Change artificial results from 999 to NA (robmon.r, phase3.r)
+ * For ML estimation: added autocorrelations during phase 3
+ to print.summary.sienaFit (sienaprint.r)
+
2013-11-29 R-Forge Revision 249
Changes in RSiena and RSienaTest:
* Fix warning in EpochSimulation regarding String-plus-int.
Modified: pkg/RSiena/inst/doc/RSiena.bib
===================================================================
--- pkg/RSiena/inst/doc/RSiena.bib 2013-11-29 12:44:07 UTC (rev 249)
+++ pkg/RSiena/inst/doc/RSiena.bib 2013-12-04 16:15:33 UTC (rev 250)
@@ -237,7 +237,19 @@
Address = {Chicago},
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 250
More information about the Rsiena-commits
mailing list