[Rsiena-commits] r289 - in pkg: RSiena RSiena/R RSiena/data RSiena/inst/doc RSiena/man RSiena/src/model/effects RSiena/src/model/effects/generic RSiena/src/model/tables RSiena/tests RSienaTest RSienaTest/R RSienaTest/data RSienaTest/inst/doc RSienaTest/man RSienaTest/src/model/effects RSienaTest/src/model/effects/generic RSienaTest/src/model/tables RSienaTest/tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Sep 10 21:43:32 CEST 2015
Author: tomsnijders
Date: 2015-09-10 21:43:31 +0200 (Thu, 10 Sep 2015)
New Revision: 289
Added:
pkg/RSiena/src/model/effects/AlterCovariateActivityEffect.cpp
pkg/RSiena/src/model/effects/AlterCovariateActivityEffect.h
pkg/RSiena/src/model/effects/DoubleInPopEffect.cpp
pkg/RSiena/src/model/effects/DoubleInPopEffect.h
pkg/RSiena/src/model/effects/DoubleOutActEffect.cpp
pkg/RSiena/src/model/effects/DoubleOutActEffect.h
pkg/RSiena/src/model/effects/generic/DoubleOutActFunction.cpp
pkg/RSiena/src/model/effects/generic/DoubleOutActFunction.h
pkg/RSiena/src/model/effects/generic/InJaccardFunction.cpp
pkg/RSiena/src/model/effects/generic/InJaccardFunction.h
pkg/RSiena/src/model/effects/generic/MixedOutStarFunction.cpp
pkg/RSiena/src/model/effects/generic/MixedOutStarFunction.h
pkg/RSiena/src/model/effects/generic/OutJaccardFunction.cpp
pkg/RSiena/src/model/effects/generic/OutJaccardFunction.h
pkg/RSienaTest/man/extract.sienaBayes.Rd
pkg/RSienaTest/src/model/effects/AlterCovariateActivityEffect.cpp
pkg/RSienaTest/src/model/effects/AlterCovariateActivityEffect.h
pkg/RSienaTest/src/model/effects/DoubleInPopEffect.cpp
pkg/RSienaTest/src/model/effects/DoubleInPopEffect.h
pkg/RSienaTest/src/model/effects/DoubleOutActEffect.cpp
pkg/RSienaTest/src/model/effects/DoubleOutActEffect.h
pkg/RSienaTest/src/model/effects/generic/DoubleOutActFunction.cpp
pkg/RSienaTest/src/model/effects/generic/DoubleOutActFunction.h
pkg/RSienaTest/src/model/effects/generic/InJaccardFunction.cpp
pkg/RSienaTest/src/model/effects/generic/InJaccardFunction.h
pkg/RSienaTest/src/model/effects/generic/MixedOutStarFunction.cpp
pkg/RSienaTest/src/model/effects/generic/MixedOutStarFunction.h
pkg/RSienaTest/src/model/effects/generic/OutJaccardFunction.cpp
pkg/RSienaTest/src/model/effects/generic/OutJaccardFunction.h
Modified:
pkg/RSiena/ChangeLog
pkg/RSiena/DESCRIPTION
pkg/RSiena/R/Sienatest.r
pkg/RSiena/R/effects.r
pkg/RSiena/R/globals.r
pkg/RSiena/R/initializeFRAN.r
pkg/RSiena/R/phase1.r
pkg/RSiena/R/phase3.r
pkg/RSiena/R/print01Report.r
pkg/RSiena/R/sienaDataCreate.r
pkg/RSiena/R/sienaGOF.r
pkg/RSiena/R/sienaModelCreate.r
pkg/RSiena/R/sienaRI.r
pkg/RSiena/R/sienaTimeTest.r
pkg/RSiena/R/sienaprint.r
pkg/RSiena/data/allEffects.csv
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/includeInteraction.Rd
pkg/RSiena/man/siena07.Rd
pkg/RSiena/man/sienaAlgorithmCreate.Rd
pkg/RSiena/man/sienaGOF.Rd
pkg/RSiena/man/sienaRI.Rd
pkg/RSiena/src/model/effects/AllEffects.h
pkg/RSiena/src/model/effects/EffectFactory.cpp
pkg/RSiena/src/model/effects/MixedNetworkEffect.cpp
pkg/RSiena/src/model/effects/MixedNetworkEffect.h
pkg/RSiena/src/model/effects/generic/CovariateMixedNetworkAlterFunction.h
pkg/RSiena/src/model/effects/generic/MixedNetworkAlterFunction.cpp
pkg/RSiena/src/model/effects/generic/MixedNetworkAlterFunction.h
pkg/RSiena/src/model/tables/TwoNetworkCache.cpp
pkg/RSiena/src/model/tables/TwoNetworkCache.h
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/effects.r
pkg/RSienaTest/R/globals.r
pkg/RSienaTest/R/initializeFRAN.r
pkg/RSienaTest/R/phase1.r
pkg/RSienaTest/R/phase3.r
pkg/RSienaTest/R/print01Report.r
pkg/RSienaTest/R/sienaBayes.r
pkg/RSienaTest/R/sienaDataCreate.r
pkg/RSienaTest/R/sienaGOF.r
pkg/RSienaTest/R/sienaModelCreate.r
pkg/RSienaTest/R/sienaRI.r
pkg/RSienaTest/R/sienaTimeTest.r
pkg/RSienaTest/R/sienaprint.r
pkg/RSienaTest/data/allEffects.csv
pkg/RSienaTest/inst/doc/RSiena.bib
pkg/RSienaTest/inst/doc/RSiena_Manual.pdf
pkg/RSienaTest/inst/doc/RSiena_Manual.tex
pkg/RSienaTest/man/RSiena-package.Rd
pkg/RSienaTest/man/includeInteraction.Rd
pkg/RSienaTest/man/siena07.Rd
pkg/RSienaTest/man/sienaAlgorithmCreate.Rd
pkg/RSienaTest/man/sienaGOF.Rd
pkg/RSienaTest/man/sienaRI.Rd
pkg/RSienaTest/src/model/effects/AllEffects.h
pkg/RSienaTest/src/model/effects/EffectFactory.cpp
pkg/RSienaTest/src/model/effects/MixedNetworkEffect.cpp
pkg/RSienaTest/src/model/effects/MixedNetworkEffect.h
pkg/RSienaTest/src/model/effects/generic/CovariateMixedNetworkAlterFunction.h
pkg/RSienaTest/src/model/effects/generic/MixedNetworkAlterFunction.cpp
pkg/RSienaTest/src/model/effects/generic/MixedNetworkAlterFunction.h
pkg/RSienaTest/src/model/tables/TwoNetworkCache.cpp
pkg/RSienaTest/src/model/tables/TwoNetworkCache.h
pkg/RSienaTest/tests/parallel.Rout.save
Log:
New versions 1.1-289. Various minor changes, new effects.
Modified: pkg/RSiena/ChangeLog
===================================================================
--- pkg/RSiena/ChangeLog 2015-07-18 14:54:37 UTC (rev 288)
+++ pkg/RSiena/ChangeLog 2015-09-10 19:43:31 UTC (rev 289)
@@ -1,27 +1,63 @@
-2015-07-18 R-Forge Revision 288
+2015-09-10 R-Forge Revision 289
Changes in RSiena and RSienaTest:
- * New effects homXOutAct, FFDeg, BBDeg, RFDeg, diffXTransTrip.
+ * New defaults for siena07() in sienaAlgorithmCreate():
+ doubleAveraging=0, diagonalize=0.2 (for MoM).
+ * Improved one-step approximations to expected Mahalanobis
+ distances in sienaGOF() (control variates for score function).
+ * Permit 3-way interactions with one ego and two dyadic effects
+ (initializeFRAN.r) (this was erroneously not allowed).
+ * New effects Jin, Jout, JinMix, JoutMix, altXOutAct,
+ doubleInPop, doubleOutAct.
+ * print01Report() now reports in-degrees also for two-mode networks.
+ * Better error handling for sienaTimeTest and scoreTest
+ (EvaluateTestStatistic() in Sienatest.r).
+ * inOutAss is dyadic.
+ * Check for positive derivative matrix at the end of phase 1
+ omitted for effects with fixed parameters
+ (phase1.r, CalculateDerivative()).
+ * Corrected effectName and functionName of inPopIntn, outPopIntn,
+ inActIntn, and outActIntn ('in' and 'out' were missing).
+ * Extended error message for "invalid effect" in initializeFRAN().
+ * Console message when allowOnly=TRUE is active in sienaDependent().
+ * Some change in includeInteraction.Rd.
+Changes in RSiena:
+ * New effects homXOutAct, FFDeg, BBDeg, RFDeg, diffXTransTrip
+ (ported from RSienaTest).
* sameXInPop and diffXInPop also added for two-mode networks;
but they are not dyadic!
* In names of behavior effects and statistics dropped the
- (redundant) parts "behavior" and "beh.".
- * plot.sienaRI: new parameter "actors";
+ (redundant) parts "behavior" and "beh.".
+Changes in RSienaTest:
+ * New function extract.sienaBayes() (bayesTest.r).
+ * sienaBayes(): options diagonalize=0.2, doubleAveraging=0 for estimation
+ of initial models in initialization phase; save initial results in case of
+ divergence during initialization phase;
+ check for initialEstimates done only for non-rate parameters.
+
+2015-07-18 R-Forge Revision 288
+Changes in RSiena and RSienaTest:
+ * plot.sienaRI: new parameter "actors";
proportions with piechart improved; function of parameter radius changed.
* fra2 constructed only if findiff (phase3.r).
* siena.table does no more produce the double minus sign in html output.
Changes in RSiena:
* Correction of error for two-mode networks in sienaGOF.
Changes in RSienaTest:
- * sienaBayes: new parameter nSampRates;
- correction in use of prevBayes (with extra createStores());
- more efficient calculation of multivariate normal density.
+ * New effects homXOutAct, FFDeg, BBDeg, RFDeg, diffXTransTrip.
+ * sameXInPop and diffXInPop also added for two-mode networks;
+ but they are not dyadic!
+ * In names of behavior effects and statistics dropped the
+ (redundant) parts "behavior" and "beh.".
+ * sienaBayes: new parameter nSampRates;
+ correction in use of prevBayes (with extra createStores());
+ more efficient calculation of multivariate normal density.
* Small changes in HowToCommit.tex.
2015-06-02 R-Forge Revision 286
(Version number was kept the same...)
Changes in RSienaTest:
* Correction of error for two-mode networks in sienaGOF.
-
+
2015-05-21 R-Forge Revision 286
Changes in RSiena and RSienaTest:
* test11 dropped from \tests\parallel.R.
@@ -45,11 +81,11 @@
diagonalize < 1 (phase2.r).
* When there are missings in constant or changing monadic covariates,
and centered=FALSE for their creation by coCovar() or
- varCovar(), the mean will be imputed (used to be 0, which was an error).
- For changing covariates, this is the global mean.
+ varCovar(), the mean will be imputed (used to be 0, which was an error).
+ For changing covariates, this is the global mean.
* In coCovar() and varCovar() there is a new argument imputationValues,
which are used (if given) for imputation of missing values.
- Like all missings, they are not used for the calculation
+ Like all missings, they are not used for the calculation
of the target statistics in the Method of Moments,
but only for the simulations.
* For non-centered covariates, averageAlter and averageInAlter
@@ -61,9 +97,9 @@
value(ego)=value(alter) are now set appropriately at 0.5 (was 0).
* Decent error message when there are (almost) all NA in the dependent
behavioural variable (effects.r, function getBehaviorStartingVals()).
- * The centering within effects for similarity variables at distance 2
+ * The centering within effects for similarity variables at distance 2
now is done by the same similarity means just as for the simX effect
- (attr(mydata$cCovars$mycov, "simMean"), etc.).
+ (attr(mydata$cCovars$mycov, "simMean"), etc.).
(Background: this was done earlier by "simMeans",
the implementation of which contained an error;
but the differentiation between "simMean" and "simMeans" is not important,
@@ -135,7 +171,7 @@
* sienaBayes: the stop caused by singularity of the precision
matrix after the multi-group estimation now is circumvented;
still a warning is printed to the screen. This allows the use of
- more than one elementary effects having the same target statistic.
+ more than one elementary effects having the same target statistic.
* sienaBayes: option priorRatesFromData changed to values 0-1-2,
with 0 = former FALSE, 1 = former TRUE, 2 = robust estimation of prior
for rate parameters from estimates at the end of initialization phase.
@@ -159,24 +195,24 @@
* Effect AltsAvAlt renamed to avXAlt (with error message).
* effects object no longer used as argument for print01Report
(with error message).
- * New effects sameXInPop, transRecTrip2,
- totAlt, avInAlt, totInAlt, totRecAlt,
- totXAlt, avXInAlt, totXInAlt,
+ * New effects sameXInPop, transRecTrip2,
+ totAlt, avInAlt, totInAlt, totRecAlt,
+ totXAlt, avXInAlt, totXInAlt,
avAltDist2, totAltDist2, avTAltDist2, totAAltDist2,
avXAltDist, totXAltDist2, avTXAltDist2, totAXAltDist2,
avInAltDist2, totInAltDist2, avTInAltDist2, totAInAltDist2,
avXInAltDist2, totXInAltDist2, avTXInAltDist2, totAXInAltDist2,
- XWX1, XWX2, cl.XWX1, cl.XWX2.
+ XWX1, XWX2, cl.XWX1, cl.XWX2.
* gwesp.. effects get endowment=TRUE in allEffects.csv, and thereby
obtain endowment and creation effects.
* In effect group behaviorBipartiteObjective, the following meaningless
effects were dropped: avAlt, avSim, totSim, avSimPopEgo, avSimPopAlt,
- behDenseTriads, simDenseTriads.
+ behDenseTriads, simDenseTriads.
* Effect class covarBehaviorBipartiteObjective split into
covarABehaviorBipartiteObjective (for covariates on first mode)
- and covarBBehaviorBipartiteObjective (for covariates on second mode)
- (allEffects.csv, effects.r).
- * Added effect class covarABehaviorBipartiteObjective
+ and covarBBehaviorBipartiteObjective (for covariates on second mode)
+ (allEffects.csv, effects.r).
+ * Added effect class covarABehaviorBipartiteObjective
and covarBBehaviorBipartiteObjective to effectsDocumentation.R.
* Undid duplication of covarNetNetObjective.
* For non-invertible covariance matrices, give diagnostic for the
@@ -185,8 +221,8 @@
* Correction of igraphNetworkExtraction() in sienaGOF-auxiliary.Rd
(the earlier version dropped isolated nodes from simulated networks).
* In sienaGOF-auxiliary.Rd, the example of constraint is replaced
- by the example of eigenvector centrality (because constraint
- is undefined for isolated nodes, leading to computational problems).
+ by the example of eigenvector centrality (because constraint
+ is undefined for isolated nodes, leading to computational problems).
* sienaRIDynamics restored, after corrections.
* "file" parameter of sienaRI dropped (implied platform dependence).
* Parameter showAll added to descriptives.sienaGOF.
@@ -209,10 +245,10 @@
Changes in RSienaTest:
* sienaBayes(): new parameters nImproveMH and priorRatesFromData;
truncate initial rate parameters depending on prior. Small change
- in example in help page.
+ in example in help page.
* glueBayes() corrected so that it can be applied sequentially.
* Small improvement in print.summary.sienaBayesFit().
- * multipleBayesTest now allows matrix parameter to test
+ * multipleBayesTest now allows matrix parameter to test
linear combinations.
* Improved plot.multipleBayesTest (shows truncation at 0).
@@ -227,7 +263,7 @@
and z$regrCor in phase1.r and phase3.r.
* print01Report() improved for changing dyadic covariates,
and some layout improvement for long variable names,
- and bug corrected for reporting upOnly/downOnly cases.
+ and bug corrected for reporting upOnly/downOnly cases.
Changes in RSienaTest:
* Adapted sienaBayes and glueBayes and print.sienaBayes to allow
inclusion of interaction effects without the corresponding
Modified: pkg/RSiena/DESCRIPTION
===================================================================
--- pkg/RSiena/DESCRIPTION 2015-07-18 14:54:37 UTC (rev 288)
+++ pkg/RSiena/DESCRIPTION 2015-09-10 19:43:31 UTC (rev 289)
@@ -1,8 +1,8 @@
Package: RSiena
Type: Package
Title: Siena - Simulation Investigation for Empirical Network Analysis
-Version: 1.1-288
-Date: 2015-07-18
+Version: 1.1-289
+Date: 2015-09-10
Author: Ruth Ripley, Krists Boitmanis, Tom A.B. Snijders
Depends: R (>= 2.15.0), utils
Imports: Matrix, tcltk, lattice, parallel, MASS
Modified: pkg/RSiena/R/Sienatest.r
===================================================================
--- pkg/RSiena/R/Sienatest.r 2015-07-18 14:54:37 UTC (rev 288)
+++ pkg/RSiena/R/Sienatest.r 2015-09-10 19:43:31 UTC (rev 289)
@@ -247,7 +247,7 @@
z2 <- fra[test]
if (inherits(try(id11 <- solve(d11), silent=TRUE), "try-error"))
{
- Report('Error message for inversion of d11: \n', cf)
+ cat('Error for inversion of d11 \n')
oneSided <- NA
v9 <- d22
v9[] <- NA
@@ -275,7 +275,7 @@
if (inherits(try(vav <- solve(v9), silent=TRUE), "try-error"))
## vav is the inverse variance matrix of ov
{
- Report('Error message for inversion of v9: \n', cf)
+ cat('Error for inversion of v9\n')
vav <- v9
vav[] <- NA
cvalue <- NA
Modified: pkg/RSiena/R/effects.r
===================================================================
--- pkg/RSiena/R/effects.r 2015-07-18 14:54:37 UTC (rev 288)
+++ pkg/RSiena/R/effects.r 2015-09-10 19:43:31 UTC (rev 289)
@@ -986,7 +986,7 @@
# restrict to covariates on second node set
covObjEffects <-
covObjEffects[covObjEffects$shortName %in%
- c("altX", "altSqX", "homXOutAct"), ]
+ c("altX", "altSqX", "homXOutAct", "altXOutAct"), ]
if (!moreThan2)
{
covObjEffects <-
Modified: pkg/RSiena/R/globals.r
===================================================================
--- pkg/RSiena/R/globals.r 2015-07-18 14:54:37 UTC (rev 288)
+++ pkg/RSiena/R/globals.r 2015-09-10 19:43:31 UTC (rev 289)
@@ -144,14 +144,14 @@
ch <- c("=", "-", " ")[level]
if (missing(dest))
{
- Report(c("@", level, "\n", text, "\n"), sep="", fill=fill)
+ Report(c("\n", "@", level, "\n", text, "\n"), sep="", fill=fill)
Report(rep(ch, sum(nchar(text)) + 3), sep="", fill=fill)
Report("\n\n")
}
else
{
dest <- deparse(substitute(dest))
- Report(c("@", level, "\n", text, "\n"), hdest=dest, sep="", fill=fill)
+ Report(c("\n", "@", level, "\n", text, "\n"), hdest=dest, sep="", fill=fill)
Report(rep(ch, sum(nchar(text))), hdest=dest, sep="", fill=fill)
if (level < 3)
{
Modified: pkg/RSiena/R/initializeFRAN.r
===================================================================
--- pkg/RSiena/R/initializeFRAN.r 2015-07-18 14:54:37 UTC (rev 288)
+++ pkg/RSiena/R/initializeFRAN.r 2015-09-10 19:43:31 UTC (rev 289)
@@ -56,7 +56,8 @@
{
bad <- which(!(userlist %in% deflist))
print(userlist[bad])
- stop("invalid effect requested: see above ")
+ cat("invalid effect requested: see above; \n")
+ stop("there seems to be a mismatch between data set and effects object.")
}
}
if (!inherits(effects, "data.frame"))
@@ -1853,11 +1854,11 @@
}
else
{
- if (egoCount < 2 && dyadCount != 3)
+ if (egoCount < 2 && (egoCount + dyadCount < 3))
{
- stop("invalid network interaction specification: ",
- "must be at least two ego or all dyadic ",
- "effects")
+ stop("invalid network 3-way interaction specification: ",
+ "must be at least two ego effects ",
+ "or all ego or dyadic effects")
}
}
## construct a name
Modified: pkg/RSiena/R/phase1.r
===================================================================
--- pkg/RSiena/R/phase1.r 2015-07-18 14:54:37 UTC (rev 288)
+++ pkg/RSiena/R/phase1.r 2015-09-10 19:43:31 UTC (rev 289)
@@ -21,10 +21,10 @@
z$SomeFixed <- FALSE
z$Phase <- 1
f <- FRANstore()
- int <- z$int
+ int <- z$int # used for multiple processes: number of processes for MoM
z <- AnnouncePhase(z, x)
z$phase1Its <- 0
- ## fix up iteration numbers if using multiple processors
+ ## fix up iteration numbers if using multiple processes
if (10 %% int == 0)
{
firstNit <- 10
@@ -130,7 +130,7 @@
##finish phase 1 iterations and do end-of-phase processing
zsmall <- makeZsmall(z)
xsmall <- NULL
- int <- z$int
+ int <- z$int # used for multiple processes: number of processes for MoM
if (z$n1 > z$phase1Its)
{
nits <- seq((z$phase1Its+1), z$n1, int)
@@ -302,12 +302,13 @@
dfra <- derivativeFromScoresAndDeviations(z$ssc, z$sf2)
+ fromBayes <- !is.null(x$fromBayes)
z$jacobianwarn1 <- rep(FALSE, z$pp)
- if (any(diag(dfra) <= 0))
+ if ((any(diag(dfra)[!z$fixed] <= 0)) && (!fromBayes))
{
for (i in 1 : z$pp)
{
- if (dfra[i, i] < 0)
+ if ((dfra[i, i] < 0) && (!z$fixed[i]))
{
##browser()
z$jacobianwarn1[i] <- TRUE
@@ -344,6 +345,14 @@
someFixed <- FALSE
if (any(diag(dfra) <= 0 & !z$fixed))
{
+ if (fromBayes)
+ {
+ # Use the values for diag(dfra) from startupGlobal computed in sienaBayes
+ neg <- which((diag(dfra) <= 0 & !z$fixed)[!z$effects$basicRate])
+ diag(dfra)[!z$effects$basicRate][neg] <- x$ddfra[neg]
+ }
+ else
+ {
neg<- which(diag(dfra) <= 0 & !z$fixed)
z$fixed[neg] <- TRUE
z$newFixed[neg] <- TRUE
@@ -352,6 +361,7 @@
' is/are fixed.\n'), cf)
Report(c('Estimation problem with parameter(s)', neg,
'this/these parameter(s) is/are fixed.\n'), outf)
+ }
}
z$dfra <- dfra
z$cdSomeFixed <- someFixed
Modified: pkg/RSiena/R/phase3.r
===================================================================
--- pkg/RSiena/R/phase3.r 2015-07-18 14:54:37 UTC (rev 288)
+++ pkg/RSiena/R/phase3.r 2015-09-10 19:43:31 UTC (rev 289)
@@ -295,6 +295,8 @@
sep="", collapse=" + ")
Report(thetext, outf)
Report('\n',outf)
+ if (is.null(x$fromBayes))
+ {
cat('*** Warning: Covariance matrix not positive definite *** \n')
cat('*** Standard errors not reliable ***\n')
cat('The following is approximately a linear combination \n')
@@ -311,6 +313,7 @@
Report('Do not use any reported standard errors.\n', outf)
errorMessage.cov <- '*** Warning: Noninvertible estimated covariance matrix ***'
}
+ }
z$msfinv <- NULL
cov.est <- NA * z$msfc
}
Modified: pkg/RSiena/R/print01Report.r
===================================================================
--- pkg/RSiena/R/print01Report.r 2015-07-18 14:54:37 UTC (rev 288)
+++ pkg/RSiena/R/print01Report.r 2015-09-10 19:43:31 UTC (rev 289)
@@ -178,15 +178,19 @@
tmpdepvar <- depvar[, , k]
use <- tmpdepvar %in% c(10, 11)
tmpdepvar[use] <- tmpdepvar[use] - 10
+ if (attr(depvar, "type") != "bipartite")
+ {
+ diag(tmpdepvar) <- 0
+ }
outdeg <- rowSums(tmpdepvar, na.rm=TRUE)
indeg <- colSums(tmpdepvar, na.rm=TRUE)
- diag(tmpdepvar) <- 0
missrow <- rowSums(is.na(tmpdepvar))
misscol <- colSums(is.na(tmpdepvar))
}
if (attr(depvar, "type") == "bipartite")
{
tmp <- format(cbind(1:atts$netdims[1], outdeg))
+ tmp2 <- format(cbind(1:atts$netdims[2], indeg))
}
else
{
@@ -197,9 +201,14 @@
Report(tmp[, 1], fill=60, outf)
Report("out-degrees\n", outf)
Report(tmp[, 2], fill=60, outf)
- if (attr(depvar, "type") != "bipartite")
+ if (attr(depvar, "type") == "bipartite")
{
Report("in-degrees\n", outf)
+ Report(tmp2[, 2], fill=60, outf)
+ }
+ else
+ {
+ Report("in-degrees\n", outf)
Report(tmp[, 3], fill=60, outf)
}
## report structural values
@@ -290,20 +299,8 @@
", number of missing values ",
"are:\n"),
sep="", outf)
- if (attr(depvar, "type") != "bipartite")
+ if (attr(depvar, "type") == "bipartite")
{
- Report("Nodes\n", outf)
- tmp <- format(cbind(1:atts$netdims[1],
- missrow, misscol))
- Report(tmp[, 1], fill=60, outf)
- Report("missing in rows\n", outf)
- Report(tmp[, 2], fill=60, outf)
- Report("missing in columns\n", outf)
- Report(tmp[, 3], fill=60, outf)
- mult <- atts$netdims[1] - 1
- }
- else
- {
Report("Senders\n", outf)
tmp <- format(cbind(1:atts$netdims[1],
missrow))
@@ -318,6 +315,18 @@
Report(tmp[, 2], fill=60, outf)
mult <- atts$netdims[2]
}
+ else
+ {
+ Report("Nodes\n", outf)
+ tmp <- format(cbind(1:atts$netdims[1],
+ missrow, misscol))
+ Report(tmp[, 1], fill=60, outf)
+ Report("missing in rows\n", outf)
+ Report(tmp[, 2], fill=60, outf)
+ Report("missing in columns\n", outf)
+ Report(tmp[, 3], fill=60, outf)
+ mult <- atts$netdims[1] - 1
+ }
Report(c("Total number of missing data: ",
sum(missrow),
", corresponding to a fraction of ",
@@ -1277,15 +1286,15 @@
length(atts$vCovars) > 0)
{
netnames <- atts$netnames
- vCovarSim2 <-
- lapply(data, function(x)
- lapply(x$vCovars, function(y) attr(y, "simMeans")))
- behSim2 <-
- lapply(data, function(x)
- lapply(x$depvars, function(y) attr(y, "simMeans")))
- cCovarSim2 <-
- lapply(data, function(x)
- lapply(x$cCovars, function(y) attr(y, "simMeans")))
+# vCovarSim2 <-
+# lapply(data, function(x)
+# lapply(x$vCovars, function(y) attr(y, "simMeans")))
+# behSim2 <-
+# lapply(data, function(x)
+# lapply(x$depvars, function(y) attr(y, "simMeans")))
+# cCovarSim2 <-
+# lapply(data, function(x)
+# lapply(x$cCovars, function(y) attr(y, "simMeans")))
if (nData > 1)
{
vCovarSim <-
Modified: pkg/RSiena/R/sienaDataCreate.r
===================================================================
--- pkg/RSiena/R/sienaDataCreate.r 2015-07-18 14:54:37 UTC (rev 288)
+++ pkg/RSiena/R/sienaDataCreate.r 2015-09-10 19:43:31 UTC (rev 289)
@@ -604,6 +604,7 @@
attr(depvars[[i]], 'noMissing') <- rep(0, observations)
attr(depvars[[i]], 'noMissingEither') <- rep(0, observations - 1)
attr(depvars[[i]], 'nonMissingEither') <- rep(0, observations - 1)
+ someOnly <- FALSE
if (type == 'behavior')
{
attr(depvars[[i]], 'noMissing') <- FALSE
@@ -631,9 +632,15 @@
if (attr(depvars[[i]], 'allowOnly'))
{
if (all(mydiff >= 0, na.rm=TRUE))
+ {
attr(depvars[[i]], 'downonly')[j] <- TRUE
+ someOnly <- TRUE
+ }
if (all(mydiff <= 0, na.rm=TRUE))
+ {
attr(depvars[[i]], 'uponly')[j] <- TRUE
+ someOnly <- TRUE
+ }
}
}
rr <- range(depvars[[i]], na.rm=TRUE)
@@ -700,9 +707,15 @@
if (attr(depvars[[i]], 'allowOnly'))
{
if (all(mydiff at x >= 0, na.rm=TRUE))
+ {
attr(depvars[[i]], 'uponly')[j] <- TRUE
+ someOnly <- TRUE
+ }
if (all(mydiff at x <= 0, na.rm=TRUE))
+ {
attr(depvars[[i]], 'downonly')[j] <- TRUE
+ someOnly <- TRUE
+ }
}
}
else
@@ -735,9 +748,15 @@
if (attr(depvars[[i]], 'allowOnly'))
{
if (all(mydiff >= 0, na.rm=TRUE))
+ {
attr(depvars[[i]], 'uponly')[j] <- TRUE
+ someOnly <- TRUE
+ }
if (all(mydiff <= 0, na.rm=TRUE))
+ {
attr(depvars[[i]], 'downonly')[j] <- TRUE
+ someOnly <- TRUE
+ }
}
}
}
@@ -936,6 +955,12 @@
}
attr(depvars[[i]], 'name') <- names(depvars)[i]
}
+ if (someOnly)
+ {
+cat('For some variables, in some periods, there are only increases, or only decreases.\n')
+cat('This will be respected in the simulations.\n')
+cat('If this is not desired, use allowOnly=FALSE when creating the dependent variables.\n')
+ }
## create the object
z <- NULL
z$nodeSets <- nodeSets
Modified: pkg/RSiena/R/sienaGOF.r
===================================================================
--- pkg/RSiena/R/sienaGOF.r 2015-07-18 14:54:37 UTC (rev 288)
+++ pkg/RSiena/R/sienaGOF.r 2015-09-10 19:43:31 UTC (rev 289)
@@ -123,6 +123,11 @@
{
ttcSimulation <- system.time( simStatsByPeriod <- lapply(period,
function (j) {
+ if (verbose)
+ {
+ cat(" Period ", j, "\n")
+ flush.console()
+ }
simStatsByPeriod <- sapply(1:iterations, function (i)
{
if (verbose && (i %% 100 == 0) )
@@ -237,25 +242,21 @@
mhdTemplate <- rep(0, sum(sienaFitObject$test))
names(mhdTemplate) <- rep(0, sum(sienaFitObject$test))
+ JoinedOneStepMHD_old <- mhdTemplate
+ OneStepMHD_old <- lapply(period, function(i) (mhdTemplate))
JoinedOneStepMHD <- mhdTemplate
- JoinedPartialOneStepMHD <- mhdTemplate
- JoinedGmmMhdValue <- mhdTemplate
-
OneStepMHD <- lapply(period, function(i) (mhdTemplate))
- PartialOneStepMHD <- lapply(period, function(i) (mhdTemplate))
- GmmMhdValue <- lapply(period, function(i) (mhdTemplate))
obsMhd <- NULL
- ExpStat <- lapply(period, function(i) {
- colMeans(simStatsByPeriod[[i]])
- })
+ ExpStat <-
+ lapply(period, function(i) {colMeans(simStatsByPeriod[[i]])})
+ simStatsByPeriod_tilde <-
+ lapply(period, function(i) {
+ t(apply(simStatsByPeriod[[i]],1, function(x){x - ExpStat[[i]]}))})
+
OneStepSpecs <- matrix(0, ncol=sum(sienaFitObject$test),
nrow=length(sienaFitObject$theta))
- PartialOneStepSpecs <- matrix(0, ncol=sum(sienaFitObject$test),
- nrow=length(sienaFitObject$theta))
- GmmOneStepSpecs <- matrix(0, ncol=sum(sienaFitObject$test),
- nrow=length(sienaFitObject$theta))
if (robust) {
covInvByPeriod <- lapply(period, function(i) ginv(
cov.rob(simStatsByPeriod[[i]]) ))
@@ -276,29 +277,18 @@
effectsObject <- sienaFitObject$requestedEffects
nSims <- sienaFitObject$Phase3nits
for (i in period) {
+ names(OneStepMHD_old[[i]]) <-
+ effectsObject$effectName[sienaFitObject$test]
names(OneStepMHD[[i]]) <-
effectsObject$effectName[sienaFitObject$test]
- names(PartialOneStepMHD[[i]]) <-
+ }
+ names(JoinedOneStepMHD_old) <-
effectsObject$effectName[sienaFitObject$test]
- names(GmmMhdValue[[i]]) <-
- effectsObject$effectName[sienaFitObject$test]
- }
names(JoinedOneStepMHD) <-
effectsObject$effectName[sienaFitObject$test]
- names(JoinedPartialOneStepMHD) <-
- effectsObject$effectName[sienaFitObject$test]
- names(JoinedGmmMhdValue) <-
- effectsObject$effectName[sienaFitObject$test]
rownames(OneStepSpecs) <- effectsObject$effectName
colnames(OneStepSpecs) <- effectsObject$effectName[sienaFitObject$test]
- rownames(PartialOneStepSpecs) <- effectsObject$effectName
- colnames(PartialOneStepSpecs) <-
- effectsObject$effectName[sienaFitObject$test]
- rownames(GmmOneStepSpecs) <- effectsObject$effectName
- colnames(GmmOneStepSpecs) <-
- effectsObject$effectName[sienaFitObject$test]
-
counterTestEffects <- 0
for(index in which(sienaFitObject$test)) {
if (verbose) {
@@ -337,68 +327,91 @@
mmThetaDelta <- as.numeric(ScoreTest(length(doTests), D,
sigma, fra, doTests,
maxlike=sienaFitObject$maxlike)$oneStep )
- mmPartialThetaDelta <- rep(0,length(theta0))
- mmPartialThetaDelta[length(theta0)] <- mmThetaDelta[length(theta0)]
# \mu'_\theta(X)
+ JacobianExpStat_old <- lapply(period, function (i) {
+ t(SF[,i,]) %*% simStatsByPeriod[[i]]/ nSims })
JacobianExpStat <- lapply(period, function (i) {
- t(SF[,i,]) %*% simStatsByPeriod[[i]]/ nSims })
+ t(SF[,i,]) %*% simStatsByPeriod_tilde[[i]]/ nSims })
# List structure: Period, effect index
thetaIndices <- 1:sum(effectsToInclude)
- # \Gamma_i(\theta)
- ExpStatCovar <- lapply(period, function (i) {
+ # \Gamma_i(\theta) i=period, j=parameter, k=replication
+ ExpStatCovar_old <- lapply(period, function (i) {
lapply(thetaIndices, function(j){
Reduce("+", lapply(1:nSims,function(k){
simStatsByPeriod[[i]][k,] %*% t(simStatsByPeriod[[i]][k,]) * SF[k,i,j]
})) / nSims
- - JacobianExpStat[[i]][j,] %*% t(ExpStat[[i]]) - ExpStat[[i]] %*% t(JacobianExpStat[[i]][j,])
+ - JacobianExpStat[[i]][j,] %*%
+ t(ExpStat[[i]]) - ExpStat[[i]] %*% t(JacobianExpStat[[i]][j,])
})
})
+ ExpStatCovar <- lapply(period, function (i) {
+ lapply(thetaIndices, function(j){
+ Reduce("+", lapply(1:nSims,function(k){
+ simStatsByPeriod_tilde[[i]][k,] %*%
+ t(simStatsByPeriod_tilde[[i]][k,]) * SF[k,i,j] })) / nSims})})
# \Xi_i(\theta)
+ JacobianCovar_old <- lapply(period, function (i) {
+ lapply(thetaIndices, function(j){
+ -1 * covInvByPeriod[[i]] %*% ExpStatCovar_old[[i]][[j]] %*%
+ covInvByPeriod[[i]] })
+ })
JacobianCovar <- lapply(period, function (i) {
lapply(thetaIndices, function(j){
- -1 * covInvByPeriod[[i]] %*% ExpStatCovar[[i]][[j]] %*% covInvByPeriod[[i]] })
+ -1 * covInvByPeriod[[i]] %*% ExpStatCovar[[i]][[j]] %*%
+ covInvByPeriod[[i]] })
})
+ Gradient_old <- lapply(period, function(i) {
+ sapply(thetaIndices, function(j){
+ ( obsStatsByPeriod[[i]] - ExpStat[[i]] ) %*%
+ JacobianCovar_old[[i]][[j]] %*%
+ t( obsStatsByPeriod[[i]] - ExpStat[[i]] )
+ })
+ -2 * JacobianExpStat_old[[i]] %*% covInvByPeriod[[i]] %*%
+ t( obsStatsByPeriod[[i]] - ExpStat[[i]] )
+ })
Gradient <- lapply(period, function(i) {
sapply(thetaIndices, function(j){
( obsStatsByPeriod[[i]] - ExpStat[[i]] ) %*%
JacobianCovar[[i]][[j]] %*%
t( obsStatsByPeriod[[i]] - ExpStat[[i]] )
})
- -2 * JacobianExpStat[[i]] %*%
- covInvByPeriod[[i]] %*%
+ -2 * JacobianExpStat[[i]] %*% covInvByPeriod[[i]] %*%
t( obsStatsByPeriod[[i]] - ExpStat[[i]] )
})
- OneStepSpecs[effectsToInclude,counterTestEffects] <- theta0 + mmThetaDelta
- PartialOneStepSpecs[effectsToInclude,counterTestEffects] <- theta0 + mmPartialThetaDelta
+ OneStepSpecs[effectsToInclude,counterTestEffects] <-
+ theta0 + mmThetaDelta
for (i in 1:length(obsMhd)) {
- OneStepMHD[[i]][counterTestEffects] <- as.numeric(obsMhd[i] + mmThetaDelta %*% Gradient[[i]] )
- PartialOneStepMHD[[i]][counterTestEffects] <- as.numeric( obsMhd[i] + mmPartialThetaDelta %*% Gradient[[i]] )
+ OneStepMHD_old[[i]][counterTestEffects] <-
+ as.numeric(obsMhd[i] + mmThetaDelta %*% Gradient_old[[i]] )
}
- JoinedOneStepMHD[counterTestEffects] <- Reduce("+",OneStepMHD)[counterTestEffects]
- JoinedPartialOneStepMHD[counterTestEffects] <- Reduce("+",PartialOneStepMHD)[counterTestEffects]
+ for (i in 1:length(obsMhd)) {
+ OneStepMHD[[i]][counterTestEffects] <-
+ as.numeric(obsMhd[i] + mmThetaDelta %*% Gradient[[i]] )
}
+ JoinedOneStepMHD_old[counterTestEffects] <-
+ Reduce("+",OneStepMHD_old)[counterTestEffects]
+ JoinedOneStepMHD[counterTestEffects] <-
+ Reduce("+",OneStepMHD)[counterTestEffects]
+ } # end 'for index'
}
names(res) <- names(obsStats)
class(res) <- "sienaGOF"
attr(res, "scoreTest") <- (sum(sienaFitObject$test) > 0)
attr(res, "originalMahalanobisDistances") <- obsMhd
+ attr(res, "oneStepMahalanobisDistances") <- OneStepMHD
attr(res, "joinedOneStepMahalanobisDistances") <-
JoinedOneStepMHD
+ attr(res, "oneStepMahalanobisDistances_old") <- OneStepMHD_old
+ attr(res, "joinedOneStepMahalanobisDistances_old") <-
+ JoinedOneStepMHD_old
attr(res, "oneStepSpecs") <- OneStepSpecs
- attr(res, "partialOneStepMahalanobisDistances") <- PartialOneStepMHD
- attr(res, "joinedPartialOneStepMahalanobisDistances") <-
- JoinedPartialOneStepMHD
- attr(res, "partialOneStepSpecs") <- PartialOneStepSpecs
- attr(res, "gmmOneStepSpecs") <- GmmOneStepSpecs
- attr(res, "gmmOneStepMahalanobisDistances") <- GmmMhdValue
- attr(res, "joinedGmmOneStepMahalanobisDistances") <- JoinedGmmMhdValue
attr(res,"auxiliaryStatisticName") <-
attr(obsStats,"auxiliaryStatisticName")
attr(res, "simTime") <- attr(simStats,"time")
@@ -457,13 +470,13 @@
originalMhd <- attr(x, "originalMahalanobisDistances")
if (attr(x, "joined")) {
cat("-----\nCalculated joint MHD = (",
- sum(originalMhd),") for current model.\n")
+ round(sum(originalMhd),2),") for current model.\n")
}
else
{
for (j in 1:length(originalMhd)) {
cat("-----\nCalculated period ", j, " MHD = (",
- originalMhd[j],") for current model.\n")
+ round(originalMhd[j],2),") for current model.\n")
}
}
invisible(x)
@@ -476,17 +489,9 @@
if (attr(x, "scoreTest")) {
oneStepSpecs <- attr(x, "oneStepSpecs")
oneStepMhd <- attr(x, "oneStepMahalanobisDistances")
-# Tom deleted the following statements, because they lead to warnings
-# when checking. In the current code they are superfluous.
-# gmmMhd <- attr(x, "gmmOneStepMahalanobisDistances")
-# gmmOneStepSpecs <- attr(x, "gmmOneStepSpecs")
-# partialOneStepSpecs <- attr(x, "partialOneStepSpecs")
-# partialOneStepMhd <- attr(x, "partialOneStepMahalanobisDistances")
-# joinedPartialOneStepMhd <-
-# attr(x, "joinedPartialOneStepMahalanobisDistances")
joinedOneStepMhd <- attr(x, "joinedOneStepMahalanobisDistances")
-# joinedGmmMhd <- attr(x, "joinedGmmOneStepMahalanobisDistances")
- cat("\nOne-step estimates for modified models.\n")
+ cat("\nOne-step estimates and predicted Mahalanobis distances")
+ cat(" for modified models.\n")
if (attr(x, "joined")) {
for (i in 1:ncol(oneStepSpecs)) {
a <- cbind(oneStepSpecs[,i, drop=FALSE] )
@@ -494,8 +499,8 @@
rownames(b) <- c("MHD")
a <- rbind(a, b)
a <- round(a, 3)
- cat("\n**Model", colnames(a)[1], "\n")
- colnames(a) <- "MMD"
+ cat("\n**Model including", colnames(a)[1], "\n")
+ colnames(a) <- "one-step"
print(a)
}
}
@@ -508,8 +513,8 @@
rownames(b) <- c("MHD")
a <- rbind(a, b)
a <- round(a, 3)
- cat("\n**Model", colnames(a)[1], "\n")
- colnames(a) <- c("MMD")
+ cat("\n**Model including", colnames(a)[1], "\n")
+ colnames(a) <- c("one-step")
print(a)
}
}
Modified: pkg/RSiena/R/sienaModelCreate.r
===================================================================
--- pkg/RSiena/R/sienaModelCreate.r 2015-07-18 14:54:37 UTC (rev 288)
+++ pkg/RSiena/R/sienaModelCreate.r 2015-09-10 19:43:31 UTC (rev 289)
@@ -21,14 +21,14 @@
function(fn,
projname="Siena", MaxDegree=0, useStdInits=FALSE,
n3=1000, nsub=4, n2start = NULL, dolby=TRUE,
- maxlike=FALSE, diagonalize=1.0*!maxlike,
+maxlike=FALSE, diagonalize=0.2*!maxlike,
condvarno=0, condname='',
firstg=0.2, reduceg=0.5, cond=NA, findiff=FALSE, seed=NULL,
pridg=0.05, prcdg=0.05, prper=0.2, pripr=0.3, prdpr=0.3,
prirms=0.05, prdrms=0.05, maximumPermutationLength=40,
minimumPermutationLength=2, initialPermutationLength=20,
modelType=1, mult=5, simOnly=FALSE, localML=FALSE,
-truncation=5, doubleAveraging=nsub, standardizeVar=(diagonalize<1))
+truncation=5, doubleAveraging=0, standardizeVar=(diagonalize<1))
{
model <- NULL
model$projname <- projname
Modified: pkg/RSiena/R/sienaRI.r
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 289
More information about the Rsiena-commits
mailing list