[Rsiena-commits] r96 - in pkg: RSiena RSiena/R RSiena/inst/doc RSiena/man RSienaTest RSienaTest/doc RSienaTest/inst/doc
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jun 4 16:12:29 CEST 2010
Author: ripleyrm
Date: 2010-06-04 16:12:27 +0200 (Fri, 04 Jun 2010)
New Revision: 96
Added:
pkg/RSiena/man/includeTimeDummy.Rd
pkg/RSiena/man/plot.sienaTimeTest.Rd
Modified:
pkg/RSiena/R/effects.r
pkg/RSiena/R/effectsMethods.r
pkg/RSiena/R/sienaTimeTest.r
pkg/RSiena/R/sienaeffects.r
pkg/RSiena/R/sienaprint.r
pkg/RSiena/R/simstatsc.r
pkg/RSiena/changeLog
pkg/RSiena/inst/doc/s_man400.pdf
pkg/RSiena/man/includeEffects.Rd
pkg/RSiena/man/includeInteraction.Rd
pkg/RSiena/man/setEffect.Rd
pkg/RSiena/man/sienaTimeTest.Rd
pkg/RSienaTest/changeLog
pkg/RSienaTest/doc/s_man400.tex
pkg/RSienaTest/inst/doc/s_man400.pdf
Log:
Bug fixes, extensions to sienaTimeTest
Modified: pkg/RSiena/R/effects.r
===================================================================
--- pkg/RSiena/R/effects.r 2010-06-04 13:27:18 UTC (rev 95)
+++ pkg/RSiena/R/effects.r 2010-06-04 14:12:27 UTC (rev 96)
@@ -837,7 +837,8 @@
## add starting values for the other objects
if (groupx && length(x) > 1)
{
- period <- xx$observations ### periods used so far
+ period <- xx$observations ##periods used so far
+
for (group in 2:length(x))
{
xx <- x[[group]]
@@ -992,6 +993,7 @@
},
stop('error type'))
}
+ period <- period + xx$observations ##periods used so far
}
}
effects <- do.call(rbind, effects)
@@ -1035,7 +1037,7 @@
}, z = depvar, y = dif)
startRate <- tmp[1, ]
##tendency
- tmp <- rowSums(tmp[-1, ]) + 2
+ tmp <- rowSums(tmp[-1, , drop=FALSE]) + 2
tendency <- log((tmp[2] * (tmp[3] + tmp[4])) /
(tmp[4] * (tmp[1] + tmp[2])))
untrimmed <- tendency
Modified: pkg/RSiena/R/effectsMethods.r
===================================================================
--- pkg/RSiena/R/effectsMethods.r 2010-06-04 13:27:18 UTC (rev 95)
+++ pkg/RSiena/R/effectsMethods.r 2010-06-04 14:12:27 UTC (rev 96)
@@ -15,43 +15,51 @@
if (!inherits(x, "sienaEffects"))
stop("not a legitimate Siena effects object")
- nDependents <- length(unique(x$name))
- userSpecifieds <- x$shortName[x$include] %in% c("unspInt", "behUnspInt")
- endowments <- !x$type[x$include] %in% c("rate", "eval")
-
- specs <- x[x$include, c("name", "effectName", "include", "fix", "test",
- "initialValue", "parm")]
- if (nDependents == 1)
+ if (typeof(fileName)=="character")
{
- specs <- specs[, -1]
+ sink(fileName, split=TRUE)
}
- if (any(endowments))
+
+ if (nrow(x) > 0)
{
- specs <- cbind(specs, type=x[x$include, "type"])
- }
- if (any(userSpecifieds))
- {
- specs <- cbind(specs, x[x$include, c("effect1", "effect2")])
- if (any (x$effect3[x$include] > 0))
+ nDependents <- length(unique(x$name))
+ userSpecifieds <- x$shortName[x$include] %in% c("unspInt", "behUnspInt")
+ endowments <- !x$type[x$include] %in% c("rate", "eval")
+
+ specs <- x[x$include, c("name", "effectName", "include", "fix", "test",
+ "initialValue", "parm")]
+ if (nDependents == 1)
{
- specs <- cbind(specs, effect3=x[x$include, "effect3"])
+ specs <- specs[, -1]
}
+ if (any(endowments))
+ {
+ specs <- cbind(specs, type=x[x$include, "type"])
+ }
+ if (any(userSpecifieds))
+ {
+ specs <- cbind(specs, x[x$include, c("effect1", "effect2")])
+ if (any (x$effect3[x$include] > 0))
+ {
+ specs <- cbind(specs, effect3=x[x$include, "effect3"])
+ }
+ }
+ specs[, "initialValue"] <- format(round(specs$initialValue,digits=5),
+ width=10)
+ row.names(specs) <- 1:nrow(specs)
+
+ print(as.matrix(specs), quote=FALSE)
+
}
- specs[, "initialValue"] <- format(round(specs$initialValue,digits=5),
- width=10)
- row.names(specs) <- 1:nrow(specs)
-
- if (typeof(fileName)=="character")
+ else
{
- sink(fileName, split=TRUE)
+ print.data.frame(x)
}
-
- print(as.matrix(specs), quote=FALSE)
-
if (typeof(fileName)=="character")
{
sink()
}
+
invisible(x)
}
Modified: pkg/RSiena/R/sienaTimeTest.r
===================================================================
--- pkg/RSiena/R/sienaTimeTest.r 2010-06-04 13:27:18 UTC (rev 95)
+++ pkg/RSiena/R/sienaTimeTest.r 2010-06-04 14:12:27 UTC (rev 96)
@@ -1,7 +1,7 @@
##*****************************************************************************
## * SIENA: Simulation Investigation for Empirical Network Analysis
## *
-## * Web: http://stat.gamma.rug.nl/siena.html
+## * Web: http://www.stats.ox.ac.uk/~snidjers/siena
## *
## * File: sienaTimeTest.r
## *
@@ -10,90 +10,161 @@
## * sienaTimeFix, which is called to set up time dummy interacted effects.
## ****************************************************************************/
##@sienaTimeTest siena07 Does test for time homogeneity of effects
-sienaTimeTest <- function (sienaFit)
+sienaTimeTest <- function (sienaFit, effects=NULL, condition=FALSE)
{
observations <- sienaFit$f$observations
+ # There must be more than 2 observations to do a time test!
if (observations <=2)
{
stop("You must have at least three time periods to test
for non-heterogeneity across time.")
}
-## Find out which dummies have already been included in the model:
- indRateEffects <- which(sienaFit$effects$shortName=="Rate")
- indDummiedEffects <- grep("Dummy", sienaFit$effects$effectName)
- indBaseEffects <- setdiff(1:nrow(sienaFit$effects), c(indRateEffects,
+ ## Screen out the undesired effects
+ if (!is.null(effects)) {
+ escreen = setdiff(1:nrow(sienaFit$effects), effects)
+ } else {
+ escreen = 99999
+ }
+ # Identify which effects are rate parameters
+ indRateEffects <- which(sienaFit$effects$shortName[-escreen]=="Rate")
+ # Identify which effects are estimated dummy terms
+ indDummiedEffects <- grep("Dummy", sienaFit$effects$effectName[-escreen])
+ # The effects which will be tested are stored here. Take all of the
+ # effects, and take out the rate and dummy effects. These indices
+ # will have to be changed for the moments, scores, etc. after we
+ # screen the sienaFit ingredients.
+ indBaseEffects <- setdiff(1:nrow(sienaFit$effects[-escreen,]), c(indRateEffects,
indDummiedEffects))
+ baseNames=sienaFit$effects[-escreen,]$effectName[indBaseEffects]
+ # toTest will hold booleans for which of the time dummies have not
+ # been estimated and thus are candidates for time tests.
toTest <- array(TRUE, dim=c(length(indBaseEffects),
sienaFit$f$observations - 2))
+ rownames(toTest) <- sienaFit$effects[-escreen,]$effectNumber[indBaseEffects]
+ colnames(toTest) <- 2:(sienaFit$f$observations - 1)
+ # dummyByEffect gets passed from sienaTimeTest so that other functions
+ # know which dummies belong to which base effects.
dummyByEffect <- array(0, dim=c(length(indBaseEffects),
sienaFit$f$observations - 2))
- rownames(toTest) <- sienaFit$effects$effectNumber[indBaseEffects]
- colnames(toTest) <- 2:(sienaFit$f$observations - 1)
dimnames(dummyByEffect) <- dimnames(toTest)
-## Must be able to screen out DummyX ego effects that are fixed:
- dscreen <- which(sienaFit$effects$shortName=='egoX' & sienaFit$effects$fix
- & length(grep("Dummy", sienaFit$effects$effectName)) > 0)
+ # dscreen is the first important screening vector, which will determine
+ # which egoX dummies are fixed, so that we do not consider them as included
+ dscreen <- which(sienaFit$effects[-escreen,]$shortName=='egoX' &
+ sienaFit$effects[-escreen,]$fix
+ & length(grep("Dummy",
+ sienaFit$effects[-escreen,]$effectName)) > 0)
if (length(dscreen)==0)
{
dscreen <- 99999
}
- for (i in sienaFit$effects$effectNumber[sienaFit$effects$timeDummy != ',']){
- tmp <- toString(sienaFit$effects$timeDummy[
- sienaFit$effects$effectNumber == i])
+ ## If the estimation was unconditional, the rate parameters will have scores
+ ## and moments which must also be screened out. Ruth, is there a simple way to
+ ## check conditioning? I tried $conditional and it doesn't seem to do what I
+ ## intuitively expected. For now, I just check the dimensionality of the scores,
+ ## as it will match the number of included "effects" on dimension 3 if uncond.
+ ## estimation was used.
+ if (dim(sienaFit$sf2[,,-escreen])[3] == dim(sienaFit$effects[,,-escreen])[1]) {
+ rscreen <- indRateEffects
+ } else {
+ rscreen <- 99999
+ }
+ ## Go through each effect which had a time dummy included, and incorporate this
+ ## information into the toTest vector. i.e. if a time dummy was estimated, set
+ ## its element in toTest equal to FALSE so that we do not time test it
+ for (i in sienaFit$effects[-escreen,]$effectNumber
+ [sienaFit$effects[-escreen,]$timeDummy != ',']){
+ tmp <- toString(sienaFit$effects[-escreen,]$timeDummy[
+ sienaFit$effects[-escreen,]$effectNumber == i])
tmp <- strsplit(tmp, split=",", fixed=TRUE)[[1]]
if (length(which(!tmp == '')) > 0)
{
- if (tmp[1]=='isDummy' & !(i %in% sienaFit$effects$
+ ## The effect we are looking at is a time dummy.
+ if (tmp[1]=='isDummy' & !(i %in% sienaFit$effects[-escreen,]$
effectNumber[dscreen]))
{
+ ## Dont test this dummy...
toTest[rownames(toTest)==as.numeric(tmp[3]),
- colnames(toTest)==as.numeric(tmp[2])] <- FALSE
+ colnames(toTest)==as.numeric(tmp[2])] <- FALSE
+ ## We want to be able to reference this effect given an
+ ## index for the base effect and a time period, so store
+ ## this information in dummyByEffect -- this is used
+ ## extensively in plot.sienaTimeTest
dummyByEffect[rownames(toTest)==as.numeric(tmp[3]),
- colnames(toTest)==as.numeric(tmp[2])] <-
- which(sienaFit$effects$effectNumber[-dscreen]==i)
+ colnames(toTest)==as.numeric(tmp[2])] <-
+ which(sienaFit$effects[-escreen,]$
+ effectNumber[-c(rscreen,dscreen)]==i)
}
}
else
{
+ ## The effect we are looking at had a time dummy,
+ ## nothing required for now.
next
}
}
+ ## nEffects, nSims, nameslist, nDummies convert commonly used ingredients
+ ## from sienaFit into an easily accessed form based on the screens
+ ## set up above
nEffects <- length(indBaseEffects) + sum(!toTest)
- nSims <- sienaFit$n3
+ ## With the use of multiple nodes, sometimes the sienaFit object comes back
+ ## with the wrong number of iterations!! Fixing it by looking elsewhere:
+ ## Used to be: nSims <- sienaFit$n3
+ nSims <- dim(sienaFit$sf2[,,-escreen])[1]
nameslist <- list(
Iteration=paste("it", 1:nSims, sep=""),
Wave=paste("Wave", 1:(observations - 1), sep=""),
- Effect=sienaFit$effects$effectName
+ Effect=sienaFit$effects[-escreen,]$effectName[-c(dscreen,rscreen)]
)
nDummies <- sum(toTest)
nTotalEffects <- nDummies + nEffects
- obsStats <- t(sienaFit$targets2[-dscreen, ])
- moment <- sienaFit$sf2[, , -dscreen] - rep(obsStats, each=nSims)
+ ## obsStats, moment, scores are the crucial ingredients from sienaFit which
+ ## screen for the base effects and make the rest of the code clean
+ obsStats <- t(sienaFit$targets2[-c(dscreen,rscreen,escreen), ])
+ moment <- sienaFit$sf2[, , -c(dscreen,rscreen,escreen)] - rep(obsStats, each=nSims)
+ scores <- sienaFit$ssc[ , , -c(dscreen,rscreen,escreen)]
+ ## Because the sienaFit object does not have a strict class definition,
+ ## the $sf2 and $targets2 arrays cannot be expected to always have the
+ ## proper format. The best we can do is therefore to die gracefully if
+ ## the arrays do not line up:
+ G <- array(0, dim=c(nSims, observations - 1, nEffects + nDummies))
+ SF <- array(0, dim=c(nSims, observations - 1, nEffects + nDummies))
+ if (sum(dim(G[, , 1:nEffects]) != dim(moment))+
+ sum(dim(SF[, , 1:nEffects]) != dim(scores))>0) {
+ stop("The moments and scores in your sienaFit have unexpected dimensions.\n
+ It is possible that your model specifications are not yet implemented\n
+ in sienaTimeTest. Please contact the developers.\n\nDid you include
+ the base effect?\n")
+ }
+ ## Will be used to construct the dummy names for output
dummyNames <- rep("", nDummies)
- G <- array(0, dim=c(nSims, observations - 1, nEffects + nDummies))
+ ## Set the base effects G equal to the moments from sienaFit
G[, , 1:nEffects] <- moment
+ ## inc used for incrementing through the dummies
inc <- nEffects
for (i in 1:nrow(toTest))
{
for (j in 1:ncol(toTest))
{
+ ## Go through each dummy to be tested
if (toTest[i, j])
{
inc <- inc + 1
+ ## And add scores and moments for the specific time period j+1
G[, j + 1, inc] <- moment[, j + 1, i]
dummyNames[inc-nEffects] <- paste("(*)Dummy", j + 1, ":",
nameslist$Effect[i], sep="")
}
}
}
+ ## Put names onto G for easy reference
dimnames(G) <- list(nameslist$Iteration, nameslist$Wave,
- c(nameslist$Effect[-dscreen], dummyNames))
+ c(nameslist$Effect, dummyNames))
+ ## Make the covariance matrix for the new moments
sigma <- cov(apply(G, c(1, 3), sum))
- SF <- array(0, dim=c(nSims, observations - 1, nEffects + nDummies))
- SF[, , 1:nEffects] <- sienaFit$ssc[ , , -dscreen]
+ ## Basically repeat this process for the scores:
+ SF[, , 1:nEffects] <- scores
inc <- nEffects
-## Add any dummies that have been previously estimated:
dummyProps <- list()
for (i in 1:nrow(toTest))
{
@@ -102,34 +173,48 @@
if (toTest[i, j])
{
inc <- inc + 1
- SF[, j + 1, inc] <- sienaFit$ssc[, j + 1, i]
- dummyNames[inc-nEffects] <- paste("(*)Dummy", j + 1, ":",
- nameslist$Effect[i], sep="")
+ SF[, j + 1, inc] <- scores[, j + 1, i]
+ ## Save some information on these dummies for later;
+ ## these operations dont relate directly to the scores
dummyByEffect[i, j]=inc
- dummyProps$shortName[inc] <- sienaFit$effects$shortName[i]
- dummyProps$type[inc] <- sienaFit$effects$type[i]
+ dummyProps$shortName[inc] <- sienaFit$effects[-escreen,]$shortName[i]
+ dummyProps$interaction1[inc] <- sienaFit$effects[-escreen,]$interaction1[i]
+ dummyProps$type[inc] <- sienaFit$effects[-escreen,]$type[i]
dummyProps$period[inc] <- j + 1
}
}
}
- dimnames(SF) <- list(nameslist$Iteration, nameslist$Wave,
- c(nameslist$Effect[-dscreen], dummyNames))
+ ## Copy the dimnames for G
+ dimnames(SF) <- dimnames(G)
+ ## We have now set up all of the ingredients properly, so we may proceed with
+ ## the score type test of Schweinberger (2007)
D <- derivativeFromScoresAndDeviations(SF, G)
fra <- apply(G, 3, sum) / nSims
doTests <- c(rep(FALSE, nEffects), rep(TRUE, nDummies))
jointTest <- ScoreTest(nTotalEffects, D, sigma, fra, doTests, maxlike=FALSE)
jointTestP <- 1 - pchisq(jointTest$testresOverall, nDummies)
- individualTest <- jointTest$testresulto[1:nDummies]
+ if (! condition) {
+ individualTest <- jointTest$testresulto[1:nDummies]
+ } else {
+ individualTest <- sapply(1:nDummies, function (i)
+ { doTests <- rep(FALSE, nEffects + nDummies)
+ doTests[nDummies+i] <- TRUE
+ test <- ScoreTest(nTotalEffects, D, sigma, fra, doTests, FALSE)
+ test$testresulto[1]
+ })
+ }
individualTestP <- 2 * (1-pnorm(abs(individualTest))[1:nDummies])
rownames(jointTestP) <- c("Joint Significant Test")
colnames(jointTestP) <- c("p-Val")
- thetaOneStep <- c(sienaFit$theta[-dscreen], rep(0, nDummies)) +
- jointTest$oneStep
+ thetaOneStep <- c(sienaFit$theta[-c(dscreen,rscreen,escreen)], rep(0, nDummies)) +
+ jointTest$oneStep
effectTest <- sapply(1:length(indBaseEffects), function (i)
{
doTests <- rep(FALSE, nEffects + nDummies)
tmp <- which(dummyProps$shortName ==
- sienaFit$effects$shortName[i])
+ sienaFit$effects[-escreen,]$shortName[i] &
+ dummyProps$interaction1 ==
+ sienaFit$effects[-escreen,]$interaction1[i])
if (length(tmp) > 0)
{
doTests[tmp] <- TRUE
@@ -142,17 +227,19 @@
NA
}
})
+
dim(effectTest) <- c(length(indBaseEffects), 1)
effectTestP <- round(1 - pchisq(effectTest, apply(toTest, 1, sum)), 5)
- rownames(effectTestP) <- nameslist$Effect[indBaseEffects]
+ rownames(effectTestP) <- baseNames
colnames(effectTestP) <- c("p-Val")
- thetaStar <- cbind(c(sienaFit$theta[-dscreen], rep(0, nDummies)),
+ thetaStar <- cbind(c(sienaFit$theta[-c(dscreen,rscreen,escreen)], rep(0, nDummies)),
thetaOneStep,
- round(c(2-2 * pnorm(abs(sienaFit$theta[-dscreen]/
- sqrt(diag(sienaFit$covtheta)[-dscreen]))),
+ round(c(2-2 * pnorm(abs(sienaFit$theta[-c(dscreen,rscreen,escreen)]/
+ sqrt(diag(sienaFit$covtheta)[-c(dscreen,
+ rscreen,escreen)]))),
individualTestP), 5))
colnames(thetaStar) <- c("Initial Est.", "One Step Est.", "p-Value")
- rownames(thetaStar) <- dimnames(SF)[[3]]
+ rownames(thetaStar) <- dimnames(G)[[3]]
returnObj <- list(
JointTest=jointTestP,
EffectTest=effectTestP,
@@ -169,10 +256,12 @@
DummyIndexByEffect=dummyByEffect,
DummyStdErr=sqrt(diag(jointTest$covMatrix)),
OriginalEffects=nEffects,
- OriginalThetaStderr=sqrt(diag(sienaFit$covtheta))[-dscreen],
+ OriginalThetaStderr=sqrt(diag(sienaFit$covtheta))[-c(dscreen,
+ rscreen,escreen)],
SienaFit=sienaFit,
DummyProps=dummyProps,
- ToTest=toTest
+ ToTest=toTest,
+ ScreenedEffects=setdiff(c(rscreen,escreen),99999)
)
class(returnObj) <- "sienaTimeTest"
returnObj
@@ -207,18 +296,11 @@
}
tmp <- paste(" (", 1:length(rownames(x$IndividualTest)), ") ",
rownames(x$IndividualTest), "\n", sep="")
- cat("\nUse the following indices for plotting the pairwise moment
- correlations:\n", tmp)
+ cat("\nUse the following indices for plotting:\n", tmp)
tmp <- paste(" (", 1:length(x$NonRateIndices), ") ",
rownames(x$IndividualTest)[x$NonRateIndices], "\n", sep="")
- cat("\nUse the following indices for plotting the effect-wise
- fitted parameters:\n", tmp)
- effectNames <- rownames(x$IndividualTest)
- dummies <- grepl("Dummy", effectNames)
- dummyIndex <- paste(" (", 1:sum(dummies), ") ", effectNames[dummies],
- "\n", sep="")
cat("\nIf you would like to fit time dummies to your model, use the
- following indices:\n", dummyIndex)
+ timeDummy column in your effects object.")
cat("\nType \"?sienaTimeTest\" for more information on this output.\n")
invisible(x)
}
@@ -321,14 +403,14 @@
if (length(effects)==1)
{
yaxis <- timetest$IndividualTest[as.vector(c(effects,
- timetest$DummyIndexByEffect[effects, ])), 2]
+ timetest$DummyIndexByEffect[effects, ])), 2]
dim(yaxis) <- c(1, timetest$Waves)
}
else
{
yaxis <- timetest$IndividualTest[as.vector(t(cbind(effects,
- timetest$DummyIndexByEffect[effects, ]))), 2]
+ timetest$DummyIndexByEffect[effects, ]))), 2]
yaxis <- matrix(yaxis, nrow=length(effects), ncol=timetest$Waves,
byrow=TRUE)
}
@@ -339,7 +421,7 @@
basevals[, 1] <- 0
yaxis <- yaxis + basevals
pvals <- timetest$IndividualTest[c(effects, as.vector(
- t(timetest$DummyIndexByEffect[effects, ]))), 3]
+ t(timetest$DummyIndexByEffect[effects, ]))), 3]
dummysd <- abs(c(vals / qnorm(1 - pvals / 2)))
dummysd[effects] <- timetest$OriginalThetaStderr[effects]
dim(dummysd) <- c(length(effects), timetest$Waves)
@@ -353,8 +435,7 @@
ymin=min(yaxis[i, ] - scale * abs(yaxis[i, ]))
ymax=max(yaxis[i, ] + scale * abs(yaxis[i, ]))
xyplot(yaxis[i, ] ~ xaxis,
- type = "p", main = rownames(timetest$IndividualTest)
- [timetest$NonRateIndices[effects[i]]],
+ type = "p", main = rownames(timetest$EffectTest)[effects[i]] ,
sub=paste("p=", timetest$EffectTest[effects[i]]), bty="n",
xlab="Wave", ylab="Parameter Value", auto.key=TRUE,
ylim=c(ymin, ymax), xlim=c(0, length(xaxis) + 1),
@@ -485,18 +566,19 @@
warning("Time dummy not implemented for multi-group projects")
effects$timeDummy <- ","
}
- covar <- effects$interaction1 != ""
- if (any(effects$timeDummy[covar] != ","))
- {
- warning("Time dummy not implemented for covariate effects")
- effects$timeDummy[covar] <- ","
- }
- eval <- effects$type =="eval"
- if (any(effects$timeDummy[!eval] !=','))
+ # Josh tested these covariate effects, they work as-is for sienaTimeFix.
+ # covar <- effects$interaction1 != ""
+ # if (any(effects$timeDummy[covar] != ","))
+ # {
+ # warning("Time dummy not implemented for covariate effects")
+ # effects$timeDummy[covar] <- ","
+ # }
+ implemented <- (effects$type == "eval" | effects$shortName == "RateX")
+ if (any(effects$timeDummy[!implemented] !=','))
{
warning("Time dummy effects are only implemented",
- " for one mode network effects of type eval.")
- effects$timeDummy[!eval] <- ","
+ " for one mode network effects of type eval or for RateX.")
+ effects$timeDummy[!implemented] <- ","
}
if (all(effects$timeDummy==',') )
{
@@ -505,6 +587,7 @@
}
else
{
+## One mode, eval effects, or RateX effects:
alreadyDummied <- grep("isDummy", effects$timeDummy)
effects$timeDummy[effects$timeDummy=="all"] <-
paste(2:(data$observations-1), collapse = ",")
@@ -515,8 +598,9 @@
## all of the previous dummied effects within the column.
effects <- effects[-alreadyDummied, ]
}
- dummiedEffects <- effects$effectNumber[effects$timeDummy != ',']
+ dummiedEffects <- effects$effectNumber[effects$timeDummy != ',' & (effects$type=='eval' | effects$shortName=='RateX')]
covToAdd <- NULL
+ rateCovToAdd <- NULL
dummyCombos <- list()
ctr=1
## This might need to be changed for sienaGroup:
@@ -543,9 +627,49 @@
}
if (length(tmp) > 0)
{
- dummyCombos[[ctr]]=list(effectNumber=i, periods=tmp)
- ctr=ctr + 1
- covToAdd <- unique(c(covToAdd, tmp))
+ if (effects$type[effects$effectNumber==i]=='eval') {
+ dummyCombos[[ctr]]=list(effectNumber=i, periods=tmp)
+ ctr=ctr + 1
+ covToAdd <- unique(c(covToAdd, tmp))
+ } else if (effects$shortName[effects$effectNumber==i]=='RateX') {
+ ## RateX effect, has to be dealt with differently. Just add them now:
+ for (p in tmp) {
+ dname <- paste(effects$interaction1[effects$effectNumber==i],
+ "Dummy",p,sep="")
+ base <- matrix(0,nact,nper-1)
+ ## Figure out the base values:
+ dvind <- which(names(data$cCovars) ==
+ effects$interaction1[effects$effectNumber==i])
+ ## Stick them into the right time spot
+ base[,p] <- data$cCovars[[dvind]]
+ ## make a new varCovar:
+ base <- varCovar(base)
+ base <- addAttributes.varCovar(base, name=dname)
+ data$vCovars[[length(data$vCovars)+1]] <- base
+ names(data$vCovars)[length(data$vCovars)] <- dname
+ ## Now add the rate term:
+ tmprow <- allEffects[allEffects$functionName==
+ 'Amount of change x xxxxxx' & allEffects$type=='rate'
+ & allEffects$effectGroup=='covarNonSymmetricRate', ]
+ tmprow$name <- effects$name[effects$shortName=='RateX' &
+ effects$type=='rate'][1]
+ tmprow$effectFn <- 'NULL'
+ tmprow$statisticFn <- 'NULL'
+ tmprow$netType <- 'oneMode'
+ tmprow$groupName <- 'Group1'
+ tmprow$group <- 1
+ tmprow$fix <- FALSE
+ tmprow$include <- TRUE
+ tmprow$effectNumber <- max(effects$effectNumber) + 1
+ tmprow <- tmprow[, colnames(effects)]
+ tmprow$effectName <- gsub('xxxxxx', dname, tmprow$effectName)
+ tmprow$functionName <- gsub('xxxxxx', dname, tmprow$functionName)
+ tmprow$interaction1 <- dname
+ tmprow$timeDummy <- paste('isDummy', p, i, sep=',')
+ rownames(tmprow) <- dname
+ effects <- rbind(effects, tmprow)
+ }
+ }
}
}
## Add the required covariate effects to the effect objects
@@ -583,10 +707,10 @@
rownames(tmprow) <- dname
effects <- rbind(effects, tmprow)
}
- for (i in 1:length(dummyCombos))
+ for (i in seq(along=dummyCombos))
{
baseNum=dummyCombos[[i]]$effectNumber
- for (j in 1:length(dummyCombos[[i]]$periods))
+ for (j in seq(along=dummyCombos[[i]]$periods))
{
dname <- paste("Dummy", dummyCombos[[i]]$periods[j], sep="")
dummyNum <- effects$effectNumber[rownames(effects)==dname]
@@ -626,4 +750,39 @@
list(effects=effects, data=data)
}
}
+##@includeTimeDummy DataCreate
+includeTimeDummy <- function(myeff, ..., timeDummy="all", name=myeff$name[1],
+ type="eval", interaction1="", interaction2="", include=TRUE,
+ character=FALSE)
+{
+ if (character)
+ {
+ dots <- sapply(list(...), function(x)x)
+ }
+ else
+ {
+ dots <- substitute(list(...))[-1] ##first entry is the word 'list'
+ }
+ if (length(dots) == 0)
+ {
+ stop("need some effect short names")
+ }
+ if (!character)
+ {
+ effectNames <- sapply(dots, function(x)deparse(x))
+ }
+ else
+ {
+ effectNames <- dots
+ }
+ use <- myeff$shortName %in% effectNames &
+ myeff$type==type &
+ myeff$name==name &
+ myeff$interaction1 == interaction1 &
+ myeff$interaction2 == interaction2
+ myeff[use, "timeDummy"] <- timeDummy
+ myeff[use, "include"] <- include
+ print.data.frame(myeff[use,])
+ myeff
+}
Modified: pkg/RSiena/R/sienaeffects.r
===================================================================
--- pkg/RSiena/R/sienaeffects.r 2010-06-04 13:27:18 UTC (rev 95)
+++ pkg/RSiena/R/sienaeffects.r 2010-06-04 14:12:27 UTC (rev 96)
@@ -75,17 +75,20 @@
{
shortNames <- dots
}
- ## check that we have a spare row
- ints <- myeff[myeff$name == name & myeff$shortName %in%
- c("unspInt", "behUnspInt") &
- (is.na(myeff$effect1) | myeff$effect1 == 0)&
- myeff$type == type, ]
- if (nrow(ints) == 0)
+ ## if want to include, check that we have a spare row
+ if (include)
{
- stop("No more interactions available:",
- "recreate the effects object requesting more interactions")
+ ints <- myeff[myeff$name == name & myeff$shortName %in%
+ c("unspInt", "behUnspInt") &
+ (is.na(myeff$effect1) | myeff$effect1 == 0)&
+ myeff$type == type, ]
+ if (nrow(ints) == 0)
+ {
+ stop("No more interactions available:",
+ "recreate the effects object requesting more interactions")
+ }
+ ints <- ints[1, ]
}
- ints <- ints[1, ]
## find the first underlying effect
shortName <- shortNames[1]
interact1 <- interaction1[1]
@@ -122,7 +125,7 @@
stop("Second effect not unique")
}
effect2 <- myeff[use, "effectNumber"]
- ## find the third underlying effect, if any
+ ## find the third underlying effect, if any
if (length(shortNames) > 2)
{
@@ -148,11 +151,22 @@
{
effect3 <- 0
}
- intn <- myeff$effectNumber == ints$effectNumber
- myeff[intn, "include"] <- include
- myeff[intn, c("effect1", "effect2", "effect3")] <-
- c(effect1, effect2, effect3)
-
+ if (include)
+ {
+ intn <- myeff$effectNumber == ints$effectNumber
+ myeff[intn, "include"] <- include
+ myeff[intn, c("effect1", "effect2", "effect3")] <-
+ c(effect1, effect2, effect3)
+ }
+ else
+ {
+ intn <- (myeff$effect1 == effect1) & (myeff$effect2 == effect2)
+ if (effect3 > 0)
+ {
+ intn <- intn & (myeff$effect3 == effect3)
+ }
+ myeff[intn, "include"] <- FALSE
+ }
print.data.frame(myeff[intn, c("name", "shortName", "type", "interaction1",
"interaction2", "include", "effect1", "effect2",
"effect3")])
@@ -162,6 +176,7 @@
##@setEffect DataCreate
setEffect <- function(myeff, shortName, parameter=0,
fix=FALSE, test=FALSE, initialValue=0,
+ timeDummy=",",
include=TRUE, name=myeff$name[1],
type="eval", interaction1="", interaction2="",
character=FALSE)
@@ -188,8 +203,9 @@
myeff[use, "fix"] <- fix
myeff[use, "test"] <- test
myeff[use, "initialValue"] <- initialValue
+ myeff[use, "timeDummy"] <- timeDummy
print.data.frame(myeff[use, c("name", "shortName", "type", "interaction1",
"interaction2", "include", "parm", "fix", "test",
- "initialValue")])
+ "initialValue", "timeDummy")])
myeff
}
Modified: pkg/RSiena/R/sienaprint.r
===================================================================
--- pkg/RSiena/R/sienaprint.r 2010-06-04 13:27:18 UTC (rev 95)
+++ pkg/RSiena/R/sienaprint.r 2010-06-04 14:12:27 UTC (rev 96)
@@ -115,9 +115,9 @@
{
for (j in 1:length(addtorow$command))
{
- ii <- grep(i-1, addtorow$pos[[j]])
- if (length(ii))
- if (i == 1 | addtorow$command[j] == 'Network Dynamics')
+ ii <- match(i-1, addtorow$pos[[j]])
+ if (!is.na(ii))
+ if (i == 2 | addtorow$command[j] == 'Network Dynamics')
cat( addtorow$command[j], '\n')
else
cat('\n', addtorow$command[j], '\n', sep='')
@@ -208,8 +208,16 @@
##@sienaFitThetaTable Miscellaneous
sienaFitThetaTable <- function(x, tstat=FALSE)
{
+ effects <- x$requestedEffects
pp <- x$pp
- nrates <- length(x$rate)
+ if (x$cconditional)
+ {
+ nrates <- length(x$rate)
+ }
+ else
+ {
+ nrates <- 0
+ }
pp <- pp + nrates
## mydf stores the data before formatting
mydf <- data.frame(dummy=rep(" ", pp),
@@ -285,7 +293,7 @@
if (nBehavs > 0)
{
- behEffects <- x$effects[x$effects$netType == 'behavior',]
+ behEffects <- effects[effects$netType == 'behavior',]
behNames <- unique(behEffects$name)
}
if (nBehavs > 1)
@@ -295,12 +303,12 @@
behNames)],
'> ', behEffects$effectName,
sep='')
- x$effects$effectName[x$effects$netType=='behavior'] <-
+ effects$effectName[effects$netType=='behavior'] <-
behEffects$effectName
}
mydf[nrates + (1:x$pp), 'row'] <- 1:x$pp
- mydf[nrates + (1:x$pp), 'type' ] <- x$effects$type
- mydf[nrates + (1:x$pp), 'text' ] <- x$effects$effectName
+ mydf[nrates + (1:x$pp), 'type' ] <- effects$type
+ mydf[nrates + (1:x$pp), 'text' ] <- effects$effectName
mydf[nrates + (1:x$pp), 'value' ] <- theta
mydf[nrates + (1:x$pp), 'se' ] <- ses
if (!is.null(x$tstat))
@@ -311,7 +319,7 @@
if (nBehavs > 0 && nOneModes > 0)
{
- nOneModeEff <- nrow(x$effects) - nrow(behEffects)
+ nOneModeEff <- nrow(effects) - nrow(behEffects)
addtorow$command[addsub] <-
'Behavior Dynamics'
addtorow$pos[[addsub]] <- nrates + 2 + nOneModeEff
Modified: pkg/RSiena/R/simstatsc.r
===================================================================
--- pkg/RSiena/R/simstatsc.r 2010-06-04 13:27:18 UTC (rev 95)
+++ pkg/RSiena/R/simstatsc.r 2010-06-04 14:12:27 UTC (rev 96)
@@ -679,6 +679,20 @@
networks <- vector('list', observations)
actorSet <- attr(depvar, "nodeSet")
compActorSets <- sapply(compositionChange, function(x)attr(x, "nodeSet"))
+ thisComp <- match(actorSet, compActorSets)
+ compChange <- any(!is.na(thisComp))
+ if (compChange)
+ {
+ stop("Composition change is not yet implemented for bipartite",
+ "networks")
+ action <- attr(compositionChange[[thisComp]], "action")
+ ccOption <- attr(compositionChange[[thisComp]], "ccOption")
+ }
+ else
+ {
+ ccOption <- 0
+ action <- matrix(0, nrow=attr(depvar, "netdims")[1], ncol=observations)
+ }
sparse <- attr(depvar, 'sparse')
if (sparse)
{
@@ -700,23 +714,20 @@
}
## extract missing entries
netmiss[[i]] <- netmat[is.na(netmat[,3]), , drop = FALSE]
- netmiss[[i]] <-
- netmiss[[i]][netmiss[[i]][, 1] != netmiss[[i]][, 2], ,
- drop=FALSE]
## carry forward missing values if any
- for (j in seq(along=netmiss[[i]][,1]))
+ if (i == 1) # set missings to zero
{
- if (i == 1) # set missings to zero
- {
- networks[[i]][netmiss[[i]][j, 1],
- netmiss[[i]][j, 2]] <- 0
- }
- else
- {
- networks[[i]][netmiss[[i]][j, 1], netmiss[[i]][j, 2]] <-
- networks[[i-1]][netmiss[[i]][j, 1], netmiss[[i]][j, 2]]
- }
+ netmat <- netmat[!is.na(netmat[,3]), ]
+ networks[[i]] <- spMatrix(nActors, nActors, netmat[, 1],
+ netmat[, 2], netmat[,3])
}
+ else
+ {
+ netmiss1 <- netmiss[[i]][, 1:2]
+ storage.mode(netmiss1) <- 'integer'
+ networks[[i]][netmiss1[, 1:2]] <-
+ networks[[i-1]][netmiss1[, 1:2]]
+ }
}
for (i in 1:observations)
{
@@ -753,8 +764,6 @@
x2[x2 %in% c(10, 11)] <- NA
mymat1 at x <- x1
mymat2 at x <- x2
- diag(mymat1) <- 0
- diag(mymat2) <- 0
mydiff <- mymat2 - mymat1
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 96
More information about the Rsiena-commits
mailing list