Nothing
##*****************************************************************************
## * SIENA: Simulation Investigation for Empirical Network Analysis
## *
## * Web: http://www.stats.ox.ac.uk/~snijders/siena
## *
## * File: sienaTimeTest.r
## *
## * Description: This file contains the sienaTimeTest code for testing the
## * significance of additional time dummies interacted with effects, and
## * sienaTimeFix, which is called to set up time dummy interacted effects.
## ****************************************************************************/
##@sienaTimeTest siena07 Does test for time homogeneity of effects
sienaTimeTest <- function (sienaFit, effects=NULL, excludedEffects=NULL,
condition=FALSE)
{
if (!inherits(sienaFit, "sienaFit"))
{
stop("not a legitimate sienaFit object")
}
if (!sienaFit$OK)
{
stop("Error end of estimation algorithm")
}
else if (sienaFit$termination == "UserInterrupt")
{
stop("User interrupted run, object possibly incomplete")
}
if (!is.null(sienaFit$gmm))
{
if (sienaFit$gmm)
{
stop("sienaTimeTest cannot be applied to results of GMoM estimation")
}
}
if (!is.null(sienaFit$sf2.byIteration))
{
if (!sienaFit$sf2.byIteration)
{
stop("sienaTimeTest needs sf2 by iterations (use lessMem=FALSE)")
}
}
waveNumbers <- attr(sienaFit$f, "periodNos")
nWaves <- length(waveNumbers)
fitEffects <- sienaFit$requestedEffects
## There must be more than 2 periods (more than 1 wave)
## to do a time test!
if (nWaves < 2)
{
stop("You must have at least two time periods (three waves) to test ",
"for non-heterogeneity across time.")
}
## if this is a maximum likelihood fit or finite differences we need
## the optional objects sdf2
if (sienaFit$maxlike || sienaFit$FinDiff.method)
{
if (is.null(sienaFit$sdf2[[1]][[1]]))
{
stop("rerun siena07 with the byWave option TRUE")
}
}
## get the desired effects
if (!is.null(effects))
{
## a little validation helps ensure we have at least one effect we want
if (any(!is.numeric(effects)))
{
stop("non numeric effect number requested")
}
if (any(!effects %in% 1:nrow(fitEffects)))
{
stop("Some effect number requested not found")
}
use <- 1:nrow(fitEffects) %in% effects
if (any(fitEffects$basicRate[use]))
{
stop("siena time tests are inappropriate for basic rates")
}
# if (any(grepl("Dummy", fitEffects$effectName[use])))
# {
# stop("siena time tests are inappropriate for time dummy effects")
# }
}
else
{
use <- !fitEffects$basicRate
}
## if (sienaFit$maxlike || sienaFit$FinDiff.method)
## {
## stop("Not yet implemented for finite differences or maxlike")
## }
## Identify the effects which will potentially be tested
baseInFit <- use & !grepl("Dummy", fitEffects$effectName)
nBaseEffects <- sum(baseInFit)
if (nBaseEffects == 0)
{
stop("No effects available to test")
}
fixedDummies <- fitEffects$shortName=='egoX' &
fitEffects$fix & grepl("Dummy", fitEffects$effectName)
## establish effects for top left of derivative matrix D
estimatedInFit <- use & !fixedDummies
topleftEffectNumbers <- fitEffects$effectNumber[estimatedInFit]
nRowsToTest <- nBaseEffects * nWaves
## construct the data frame that will be the basis for all calculations
toTest <- data.frame(baseEffect=rep(1:nBaseEffects, each=nWaves),
effectNumber=rep(fitEffects$effectNumber[baseInFit],
each=nWaves),
type=rep(fitEffects$type[baseInFit], each=nWaves),
period = rep(waveNumbers, nBaseEffects),
period1 = rep(1:nWaves, nBaseEffects),
toTest=rep(TRUE, nRowsToTest),
baseRowInD=rep(0, nRowsToTest),
effectName=rep(fitEffects$effectName[baseInFit],
each=nWaves),
rowInD=rep(0, nRowsToTest))
toTest$toTest[toTest$period == 1] <- FALSE
toTest$rowInD <- match(toTest$effectNumber,
topleftEffectNumbers, nomatch=0)
toTest$baseRowInD <- toTest$rowInD
toTest$rowInD[toTest$period > 1] <- 0
## Go through each effect which is an estimated time dummy, and
## incorporate this information into the toTest data frame. 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 which(grepl("isDummy", fitEffects$timeDummy) & estimatedInFit))
{
tmp <- toString(fitEffects$timeDummy[i])
tmp <- strsplit(tmp, split=",", fixed=TRUE)[[1]]
if (any(tmp != ""))
{
## Dont test the dummy for the corresponding effect
thisRow <- toTest$effectNumber == tmp[3] & toTest$period == tmp[2]
toTest[thisRow, "toTest"] <- 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 rowInD -- this is used
## extensively in plot.sienaTimeTest
toTest[thisRow, "rowInD"] <-
match(fitEffects$effectNumber[i],
topleftEffectNumbers)
}
}
## nEffects, nSims, nameslist, nDummies convert commonly used ingredients
## from sienaFit into an easily accessed form based on the screens
## set up above
nEffects <- sum(!toTest$toTest)
## add row indices in D to the rows to be tested
toTest$rowInD[toTest$toTest] <- nEffects + 1:sum(toTest$toTest)
toTest <- toTest[order(toTest$rowInD), ]
nSims <- sienaFit$Phase3nits
nameslist <- list(Iteration=paste("it", 1:nSims, sep=""),
Wave=paste("Wave", waveNumbers, sep=""),
Effect=fitEffects$effectName[estimatedInFit]
)
nDummies <- sum(toTest$toTest)
ttt <- toTest$toTest
type <- ifelse(toTest$type == "eval", "",
paste(" (", toTest$type, ")", sep=""))
toTest$dummyNames[!ttt] <- paste(fitEffects$effectName[estimatedInFit],
type[!ttt], sep="")
toTest$dummyNames[ttt] <- paste("(*)Dummy", toTest$period[ttt], ":",
toTest$effectName[ttt], type[ttt], sep="")
## obsStats, moment, scores are the crucial ingredients from sienaFit which
## screen for the base effects and make the rest of the code clean
## be careful as the original logical vectors cannot be used once rows
## have been extracted.
obsStats <- t(sienaFit$targets2[estimatedInFit, , drop=FALSE])
moment <- sienaFit$sf2[, , estimatedInFit, drop=FALSE] -
rep(obsStats, each=nSims)
G <- array(0, dim=c(nSims, nWaves, nEffects + nDummies))
## Set the base effects G equal to the moments from sienaFit
G[, , 1:nEffects] <- moment
## copy over the others
subs1 <- cbind(rep(1:nSims, nDummies), rep(toTest$period1[ttt], each=nSims),
rep(toTest$rowInD[ttt], each=nSims))
subs2 <- cbind(rep(1:nSims, nDummies), rep(toTest$period1[ttt], each=nSims),
rep(toTest$baseRowInD[ttt], each=nSims))
G[subs1] <- moment[subs2]
## Put names onto G for easy? reference
## Use dimnames(G) <- list(NULL, NULL) to view the data as a matrix!
dimnames(G) <- list(nameslist$Iteration, nameslist$Wave, toTest$dummyNames)
## Make the covariance matrix for the new moments
sigma <- cov(apply(G, c(1, 3), sum))
if (length(sigma) == 1)
{
sigma <- matrix(sigma,1,1)
}
if (!(sienaFit$maxlike || sienaFit$FinDiff.method))
{
scores <- sienaFit$ssc[ , , estimatedInFit, drop=FALSE]
SF <- array(0, dim=c(nSims, nWaves, nEffects + nDummies))
SF[, , 1:nEffects] <- scores
SF[subs1] <- scores[subs2]
## Put names onto SF for easy reference
dimnames(SF) <- dimnames(G)
D <- derivativeFromScoresAndDeviations(SF, G, , , , TRUE, )
}
else
{
if (sienaFit$FinDiff.method)
{
##sienaFit$sdf2 is a list of arrays, one for each iteration
derivs <- t(sapply(sienaFit$sdf2, function(x)
x[, estimatedInFit, estimatedInFit, drop=FALSE]))
}
else
{
##sienaFit$sdf2 is a list of lists of Matrix objects
derivs <-
t(sapply(sienaFit$sdf2, function(xx)
array(t(sapply(xx, function(x)
as.matrix(x[estimatedInFit,
estimatedInFit,
drop=FALSE]))),
dim=c(nWaves, nEffects, nEffects))))
}
derivs <- array(derivs, dim=c(nSims, nWaves, nEffects, nEffects))
DF <- array(0, dim=c(nSims, nWaves, nEffects + nDummies,
nEffects + nDummies))
DF[, , 1:nEffects, 1:nEffects] <- derivs
for (wave in 2:nWaves)
{
thisWave <- toTest$period == wave & toTest$toTest
subs1 <- (1: (nEffects + nDummies))[thisWave]
subs2 <- toTest$baseRowInD[thisWave]
DF[, wave, subs1, subs1] <- derivs[, wave, subs2, subs2]
DF[, wave, subs1, 1:nEffects] <- derivs[, wave, subs2, 1:nEffects]
DF[, wave, 1:nEffects, subs1] <- derivs[, wave, 1:nEffects, subs2]
}
D <- t(apply(DF, c(3, 4), mean))
}
## The following parts about excluding effects were added later
## (preceding code is mainly from JL; parts about excluding are from TS).
## This leads to some duplication in functionality;
## but since the preceding construction works well,
## I (TS) wanted to keep it mostly as it was/is.
## Check for fixed effects and exclude them from testing.
if (is.null(excludedEffects))
{
excludedEffects <- which(fitEffects$fix)
}
else
{
excludedEffects <-
unique(c(excludedEffects, which(fitEffects$fix)))
}
## In case some effects are excluded from testing, this is done now:
if (!is.null(excludedEffects))
{
useRows <- rep(TRUE, dim(toTest)[1])
for (i in seq_along(excludedEffects))
{
useRows[toTest$baseRowInD==excludedEffects[i]] <- FALSE
}
useRows[1:nBaseEffects] <- TRUE
toTest <- toTest[useRows,]
D <- D[useRows, useRows]
sigma <- sigma[useRows, useRows, drop=FALSE]
nDummies <- sum(toTest$toTest)
G <- G[,,useRows]
nRowsToTest <- dim(toTest)[1]
exclusions <- 1:nBaseEffects %in% excludedEffects
}
else
{
exclusions <- rep(FALSE, nBaseEffects)
}
## Check if there are any collinearities in the effects and dummies;
## note that time dummies added by sienaTimeFix were already excluded.
rankSigma <- qr(sigma)$rank
extraExclusions <- rep(FALSE, nBaseEffects)
if (rankSigma < dim(sigma)[1])
{
es <- ifelse((nRowsToTest - nDummies == 1), '', 's')
message("TimeTest constructed a null hypothesis with ",
nRowsToTest - nDummies, " estimated parameter", es)
message("and ", nDummies," dummy variables to be tested.")
message("However, there are ", dim(sigma)[1] - rankSigma,
" linear dependencies between these.\n")
message("This may be because some of the parameters are already")
message("interactions with time dummies or other time variables.")
subset0 <- rep(FALSE, dim(sigma)[1])
subset0[1:nBaseEffects] <- TRUE
for (i in 1:nBaseEffects)
{
if (!(i %in% excludedEffects))
## Check if the dummies for this effect
## are dependent on the base effects.
{
subset <- subset0 | (toTest$baseRowInD==i)
if (qr(sigma[subset,subset])$rank
< dim(sigma[subset,subset])[1])
{
extraExclusions[i] <- TRUE
}
}
}
if (any(extraExclusions))
{
useRows <- rep(TRUE, dim(toTest)[1])
for (i in 1:nBaseEffects)
{
if (extraExclusions[i])
{
useRows[toTest$baseRowInD==i] <- FALSE
}
useRows[1:nBaseEffects] <- TRUE
}
toTest <- toTest[useRows,, drop=FALSE]
D <- D[useRows, useRows, drop=FALSE]
sigma <- sigma[useRows, useRows, drop=FALSE]
nDummies <- sum(toTest$toTest)
G <- G[,,useRows, drop=FALSE]
nRowsToTest <- dim(toTest)[1]
excludedNumber <- sum(extraExclusions)
es <- ifelse((excludedNumber == 1), '', 's')
message("Automatic discovery of dependencies yielded the exclusion of effect", es)
if (excludedNumber <= 1)
{
message(paste(row(fitEffects)[baseInFit,drop=FALSE][extraExclusions,drop=FALSE][1],'.',
fitEffects$effectName[baseInFit,drop=FALSE][which(extraExclusions)],"\n "))
}
else
{
message(paste(row(fitEffects)[baseInFit,drop=FALSE][extraExclusions,drop=FALSE][,1],'.',
fitEffects$effectName[baseInFit,drop=FALSE][which(extraExclusions)],"\n "))
}
rankSigma <- qr(sigma)$rank
if (rankSigma < dim(sigma)[1])
{
cat("After these exclusions, there still are linear dependencies.\n")
cat("Advice: use sienaTimeTest with a smaller set of effects.\n\n")
stop("Please rerun sienaTimeTest with appropriate excluded effects.")
}
}
else
{
cat("Automatic discovery of dependencies had no effect.\n")
cat("Advice: use sienaTimeTest with a smaller set of effects.\n\n")
stop("Please rerun sienaTimeTest with appropriate excluded effects.")
}
}
## We have now set up all of the ingredients properly, so we may proceed
## with the score type test of Schweinberger (2007-2012)
fra <- apply(G, 3, sum) / nSims
doTests <- toTest$toTest
redundant <- rep(FALSE, length(doTests))
jointTest <- ScoreTest(nrow(toTest), D, sigma, fra, doTests, redundant,
maxlike=sienaFit$maxlike)
jointTestP <- 1 - pchisq(jointTest$testresOverall, nDummies)
if (! condition)
{
individualTest <- jointTest$testresulto[1:nDummies]
}
else
{
individualTest <- sapply(1:nDummies, function (i)
{
doTests <- rep(FALSE, nEffects + nDummies)
doTests[nEffects + i] <- TRUE
test <- ScoreTest(nrow(toTest), D, sigma,
fra, doTests, redundant, maxlike=sienaFit$maxlike)
test$testresulto[1]
}
)
}
individualTestP <- 2 * (1- pnorm(abs(individualTest)))
if (!is.na(jointTestP))
{
rownames(jointTestP) <- "Joint Significant Test"
colnames(jointTestP) <- "p-Val"
}
thetaOneStep <- c(sienaFit$theta[estimatedInFit, drop=FALSE], rep(0, nDummies)) +
jointTest$oneStep
## Define the null hypothesis for the later tests (effect-wise and group-wise)
## as the null hypothesis tested in the overall test:
nullHyp <- !toTest$toTest
effectTest <- as.vector(by(toTest, toTest$baseEffect, function (x)
{
doTests <- rep(FALSE, nEffects + nDummies)
if (any(x$toTest))
{
if (!condition)
{
## Control for all estimated effects.
doTests[toTest$baseEffect == (x$baseEffect[1]) &
toTest$toTest] <- TRUE
testresult <- partial.scoreTest(D, sigma, fra,
doTests, nullHyp, maxlike=sienaFit$maxlike)
testresult
}
else
{
## Control for everything.
doTests[toTest$baseEffect == (x$baseEffect[1]) &
toTest$toTest] <- TRUE
redundant <- rep(FALSE, length(doTests))
test <- ScoreTest(nEffects + nDummies, D, sigma,
fra, doTests, redundant, maxlike=sienaFit$maxlike)
test$testresOverall
}
}
else
{
NA
}
}
))
dim(effectTest) <- c(length(effectTest), 1)
effectTestP <- matrix(NA, dim(effectTest)[1], 3)
effectTestP[,1] <- round(effectTest[,1],2)
effectTestP[,2] <- tapply(toTest$toTest, toTest$baseEffect, sum)
effectTestP[,3] <- round(1 - pchisq(effectTestP[,1], effectTestP[,2]), 3)
rownames(effectTestP) <- toTest$dummyNames[toTest$period == 1]
colnames(effectTestP) <- c("chi-sq.", "df", "p-value")
pvalues <-
round(c(2 * (1 -
pnorm(abs(sienaFit$theta[estimatedInFit] /
sqrt(diag(sienaFit$covtheta)[estimatedInFit])))),
individualTestP), 3)
thetaStar <- cbind(c(sienaFit$theta[estimatedInFit], rep(0, nDummies)),
thetaOneStep, pvalues)
colnames(thetaStar) <- c("Initial Est.", "One Step Est.", "p-Value")
rownames(thetaStar) <- dimnames(G)[[3]]
## put things on toTest to make plot easier
toTest[, c("InitialEst", "OneStepEst", "p.value")] <- thetaStar
toTest$effectTest <- NA
toTest$effectTest <- effectTestP[toTest$baseEffect, 1]
type <- ifelse(toTest$type =="eval", "", paste(" (",
as.character(toTest$type), ")", sep=""))
toTest$effectName <-
factor(paste(toTest$effectName, type,
" \n(chisq=", toTest$effectTest, ")", sep=""))
toTest$valsplus <- toTest$OneStepEst +
ifelse(toTest$period == 1, 0, toTest$OneStepEst[toTest$baseEffect])
toTest$dummysd <- abs(toTest$OneStepEst / qnorm(1 - toTest$p.value / 2))
toTest$dummysd[toTest$period == 1] <-
sqrt(diag(sienaFit$covtheta))[estimatedInFit][toTest$period == 1]
# Added by Tom: groupwise tests.
cef <- attr(sienaFit$f,"condEffects")
if (inherits(cef,"sienaGroupEffects"))
{
# This is a multi-group result
groupsKind <- "Group"
groupSizes <- attr(sienaFit$f,"groupPeriods") - 1
nGroups <- sienaFit$f$nGroup
if (nGroups != length(groupSizes)) {stop("No good group sizes 1.")}
groups <- rep(1:nGroups, groupSizes)
if (nWaves != length(groups)) {stop("No good group sizes 2.")}
toTest$whichGroup <- groups[toTest$period1]
}
else
{
groupsKind <- "Period"
nGroups <- sienaFit$observations-1
toTest$whichGroup <- toTest$period
}
groupTest <-
as.vector(by(toTest, toTest$whichGroup, function (x)
{
doTests <- rep(FALSE, nEffects + nDummies)
if (any(x$toTest))
{
doTests[toTest$whichGroup == (x$whichGroup[1]) &
toTest$toTest] <- TRUE
testresult <- partial.scoreTest(D, sigma, fra, doTests,
nullHyp, maxlike=sienaFit$maxlike)
testresult
}
else
{
NA
}
}
))
# The preceding does not yield the test for the first group/period,
# because it was used as the reference category.
# This is a test for the sum of the dummies for the other groups.
# Therefore matrix A is created as the design matrix for this test.
remainingEffects <- unique(toTest$baseEffect[toTest$toTest])
nRemainingEffects <- length(remainingEffects)
if (nRemainingEffects >= 1)
{
A <- matrix(0, nRemainingEffects, nRowsToTest)
for (i in seq_along(remainingEffects))
{
A[i, which(toTest$baseEffect==remainingEffects[i])[2:nWaves]] <- 1
}
groupTest[1] <- transformed.scoreTest(D, sigma, fra, A,
nullHyp, maxlike=sienaFit$maxlike)
dim(groupTest) <- c(nGroups, 1)
groupTestP <- matrix(NA, nGroups, 3)
groupTestP[,1] <- round(groupTest[,1], 2)
groupTestP[,2] <- tapply(toTest$toTest, toTest$whichGroup, sum)
groupTestP[1,2] <- nRemainingEffects
groupTestP[,3] <- round(1 - pchisq(groupTestP[,1], groupTestP[,2]), 3)
rownames(groupTestP) <- paste(groupsKind, 1:nGroups)
colnames(groupTestP) <- c("chi-sq.", "df", "p-value")
}
else
{
groupTest[1] <- NA
groupTestP <- NA
cat('Nothing remains...\n')
}
returnObj <- list(
sienaFitName=deparse(substitute(sienaFit)),
JointTest=jointTestP,
EffectTest=effectTestP,
IndividualTest=thetaStar,
JointTestStatistics=jointTest,
EffectTestStatistics=effectTest,
IndividualTestStatistics=individualTest,
GroupTest=groupTestP,
GroupTestStatistics=groupTest,
GroupsKind=groupsKind,
CovDummyEst=jointTest$covMatrix,
Moments=G,
Deriv=D,
BaseRowInD=match(which(baseInFit),
which(estimatedInFit)),
Waves=dim(G)[2],
Sims=dim(G)[1],
Effects=dim(G)[3],
DummyStdErr=sqrt(diag(jointTest$covMatrix)),
OriginalEffects=nEffects,
OriginalThetaStderr=
sqrt(diag(sienaFit$covtheta))[estimatedInFit],
ToTest=toTest,
ScreenedEffects=which(!use),
RequestedExclusions=exclusions,
ExtraExclusions=extraExclusions,
WaveNumbers=waveNumbers,
IndividualTestsOrthogonalized=condition
)
class(returnObj) <- "sienaTimeTest"
returnObj
}
##@transformed.scoreTest siena07 score test for part of the additional effects
transformed.scoreTest <- function(dfra, msf, fra, A, nullHyp, maxlike)
# score test of A %*% fra;
# Boolean vector nullHyp defines vector of estimated effects
# non-zero columns of A should not overlap with nullHyp
{
test <- colSums(abs(A)) > 0
if (any(test & nullHyp))
{
stop("Null hyp and test in transformed score test should not overlap")
}
if (dim(A)[1] + sum(nullHyp) < 2)
{
stop("Meaningless transformed score test")
}
A1 <- matrix(0, sum(nullHyp), dim(dfra)[1])
A1[cbind(1:sum(nullHyp), which(nullHyp))] <- 1
AA <- rbind(A1,A)
D <- AA %*% dfra %*% t(AA)
sig <- AA %*% msf %*% t(AA)
g <- AA %*% fra
tested <- rep(FALSE, dim(AA)[1])
tested[(sum(nullHyp)+1):dim(AA)[1]] <- TRUE
redundant <- rep(FALSE, dim(AA)[1])
EvaluateTestStatistic(maxlike, tested, redundant, D, sig, g)$cvalue
}
##@partial.scoreTest siena07 score test for part of the additional effects
partial.scoreTest <- function(dfra, msf, fra, test, nullHyp, maxlike)
# test gives vector of tested effects,
# nullHyp gives vector of estimated effects
{
if (any(test & nullHyp))
{
stop("Null hypothesis and test in partial score test should not overlap")
}
if (!any(test)){stop("Nothing to be tested in partial score test")}
A <- matrix(0, sum(test), dim(dfra)[1])
A[cbind(1:sum(test),which(test))] <- 1
transformed.scoreTest(dfra, msf, fra, A, nullHyp, maxlike)
}
##@summary.sienaTimeTest siena07 summary method for sienaTimeTest objects
summary.sienaTimeTest <- function(object, ...)
{
if (!inherits(object, "sienaTimeTest"))
{
stop("not a legitimate Siena time test object")
}
class(object) <- c("summary.sienaTimeTest", class(object))
object
}
##@print.summary.sienaTimeTest siena07 print method for summary.sienaTimeTest
print.summary.sienaTimeTest <- function(x, ...)
{
if (!inherits(x, "summary.sienaTimeTest"))
{
stop("not a legitimate Siena time test summary object")
}
print.sienaTimeTest(x)
## Additional output to the print will go in here:
## The first condition is merely for compatibility with older versions.
if ((!is.null(x$RequestedExclusions))||(!is.null(x$ExtraExclusions)))
{
if (any(x$RequestedExclusions) || any(x$ExtraExclusions))
{
cat("\nSome effects were excluded from being tested.\n")
}
if (any(x$RequestedExclusions))
{
if (sum(x$RequestedExclusions) <= 1)
{
cat("The effect with number", which(x$RequestedExclusions))
cat(" in", x$sienaFitName, "was excluded.\n")
}
else
{
cat("\nRequested exclusions were the effects numbered\n")
cat(which(x$RequestedExclusions)[1],
paste(", ",
which(x$RequestedExclusions)[2:sum(x$RequestedExclusions)],
sep=""),
" in ", x$sienaFitName, ".\n", sep="")
}
}
if (any(x$ExtraExclusions))
{
if (sum(x$ExtraExclusions) <= 1)
{
cat("The effect with number", which(x$ExtraExclusions),
" in", x$sienaFitName,
" was excluded automatically because of collinearity.\n")
}
else
{
cat("\nThe effects ")
cat(which(x$ExtraExclusions)[1], paste(", ",
which(x$ExtraExclusions)[2:sum(x$ExtraExclusions)],
sep=""),
" in ", x$sienaFitName, ".\n", sep="")
cat("were excluded automatically because of collinearity.\n")
}
}
}
cat("\nIndividual significance tests and one-step estimators:\n")
print(format(round(x$IndividualTest, digits=4)), quote=F)
cat("\nEffect-wise joint significance tests\n")
cat("(i.e. each effect across all dummies):\n")
print(x$EffectTest)
if (x$Waves <=2 && ! x$IndividualTestsOrthogonalized)
{
# cat("\n\nNote that these effect-wise tests have a different form")
# cat("\nthan the individual tests, thus testing with 3 observations")
# cat("\nmay yield different individual and effect-wise values.\n\n")
}
else
{
if (x$IndividualTestsOrthogonalized)
{
cat("\nNote that the 1-df and effect-wise test statistics were",
"\northogonalized with respect to each other (condition=TRUE).")
cat("\n")
}
}
if (!is.null(x$GroupTest))
{
cat(paste("\n",x$GroupsKind,"-wise joint significance tests\n", sep=""))
cat(paste("(i.e. each", tolower(x$GroupsKind),
"across all parameters):\n"))
print(x$GroupTest)
}
tmp <- paste(" (", 1:length(x$BaseRowInD), ") ",
rownames(x$IndividualTest)[x$BaseRowInD], "\n", sep="")
cat("\nUse the following indices for plotting:\n",
tmp[!is.na(x$EffectTestStatistics)])
cat("\nIf you would like to fit time dummies to your model,")
cat("\nuse the includeTimeDummy function.")
cat("\nType \"?sienaTimeTest\" for more information on this output.\n")
invisible(x)
}
##@print.sienaTimeTest siena07 print method for sienaTimeTest objects
print.sienaTimeTest <- function(x, ...)
{
if (!inherits(x, "sienaTimeTest"))
{
stop("not a legitimate Siena time test object")
}
effectNames <- rownames(x$IndividualTest)
dummies <- x$ToTest$toTest
dummyIndex <- paste(" (", format(1:sum(dummies)), ") ",
effectNames[dummies], "\n", sep="")
pvalue <- ifelse((x$JointTest < 0.0001)," < 0.0001",
paste("=",sprintf("%6.4f", x$JointTest)))
cat("Joint significance test of time heterogeneity:\n",
"chi-squared = ",
sprintf("%6.2f", x$JointTestStatistics$testresOverall),
", d.f. = ", sprintf("%1d", length(x$IndividualTestStatistics)),
", p ", pvalue,
", \nwhere H0: The following parameters are zero:\n",
dummyIndex, sep="")
invisible(x)
}
##@plot.sienaTimeTest siena07 plot method for sienaTimeTest objects
plot.sienaTimeTest <- function(x, pairwise=FALSE, effects,
scale=0.2, plevels=c(0.1, 0.05, 0.025), ...)
{
## require(lattice)
tmp <- paste(" (", 1:length(x$BaseRowInD), ") ",
rownames(x$IndividualTest)[x$BaseRowInD], "\n", sep="")
if (missing(effects))
{
effects<- (1:length(tmp))[!is.na(x$EffectTestStatistics)]
}
if (any(!effects %in% 1:length(tmp)))
{
cat("Detected an error with the effects included. For a
parameter-plot, use the following indices:")
stop("\nUse the following indices for plotting the pairwise
moment correlations:\n", tmp)
}
if (pairwise)
{
## On a "pairwise" call, print a pairwise plot of moments
if (length(effects) == 0)
{
x <- x$Moments
}
else
{
x <- x$Moments[, , effects, drop=FALSE]
}
panel.cor <- function(x, y, digits=2, prefix="", cex.cor, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits=digits)[1]
txt <- paste(prefix, txt, sep="")
if(missing(cex.cor)) cex.cor <- 0.8 / strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(apply(x, c(1, 3), sum),
lower.panel=panel.smooth,
upper.panel=panel.cor,
pch=5, ...) ##pch=5 is a diamond
}
else
{
# require(grid)
## Otherwise, make the parameter plots:
## create a vector to split the IndividualTests by
toTest <- x$ToTest
waveNumbers <- x$WaveNumbers
nWaves <- x$Waves
nLevels <- length(plevels)
cis <- lapply(1:nLevels, function(i)
{
cbind(toTest$valsplus - abs(qnorm(plevels[i] / 2,
sd=toTest$dummysd)),
toTest$valsplus + abs(qnorm(plevels[i] / 2,
sd=toTest$dummysd)))
}
)
toTest[, paste(rep(c("lower", "upper"), nLevels),
rep(1:nLevels, each=2), sep="")] <-
do.call(cbind, cis)
## subset the data frame
toTest <- toTest[toTest$baseEffect %in% effects, ]
toTest$effectName <- factor(toTest$effectName,
levels=unique(toTest$effectName))
if (nlevels(toTest$effectName) != length(unique(toTest$baseEffect)))
{
stop("non unique effect names: contact RSiena help")
}
#devAskNewPage(TRUE)
xyplot(valsplus ~ period1 | effectName, data=toTest,
type = "p", as.table=TRUE, bty="n",
xlab="Wave", ylab="Parameter Value", auto.key=TRUE,
scales=list(relation='free', x=list(at=c(0:(nWaves+1)),
labels=c(" ", waveNumbers, " "))), subscripts=TRUE,
par.strip.text=list(lines=2.2), #layout=layout,
xlim=c(0, nWaves + 1), prepanel=function(x, y, ...)
{
ymin <- min(y - scale * abs(y))
ymax <- max(y + scale * abs(y))
list(ylim=c(ymin, ymax))
},
panel=function(x, y, subscripts, ...)
{
sapply(1:length(x), function(i, subscripts)
{ ## pch 45 is a hyphen, pch 20 a filled circle
tmp <- if(toTest$toTest[subscripts][i]) "red" else "gray"
panel.xyplot(c(x[i], x[i]),
c(toTest$lower1[subscripts][i],
toTest$upper1[subscripts][i]),
col=tmp, alpha=.50,
type="l", lend=1, lwd=10)
panel.xyplot(c(x[i], x[i]),
c(toTest$lower1[subscripts][i],
toTest$upper1[subscripts][i]),
col=tmp, alpha=.75,
type="p", pch=45, cex=3)
panel.xyplot(c(x[i], x[i]),
c(toTest$lower2[subscripts][i],
toTest$upper2[subscripts][i]),
col=tmp, alpha=.50,
type="l", lend=1, lwd=10)
panel.xyplot(c(x[i], x[i]),
c(toTest$lower2[subscripts][i],
toTest$upper2[subscripts][i]),
col=tmp, alpha=.75,
type="p", pch=45, cex=3)
panel.xyplot(c(x[i], x[i]),
c(toTest$lower3[subscripts][i],
toTest$upper3[subscripts][i]),
col=tmp, alpha=.25,
type="l", lend=1, lwd=10)
}, subscripts=subscripts)
panel.xyplot(x, y, type="s",
col="black", alpha=.75, pch=2)
panel.xyplot(x, y, type="p", pch=20,
col=1)
panel.abline(a=toTest[subscripts[1], "initial"],
reference=TRUE,
col="black", lwd=2, alpha=.75)
# grid.text(paste("p = ", toTest$effectTest[subscripts[1]]),
# 0.02, 0.1, just="left")
}, ...)
}
}
##@sienaTimeFix siena07 Adds time dummy terms to the effects object
sienaTimeFix <- function(effects, data=NULL, getDocumentation=FALSE)
{
##@addEffect internal sienaTimeFix add one or more effects
addEffect <- function(newEffects, i, newname, effectGroup, shortName,
timeDummy, fix=FALSE, include=TRUE, yName=NULL)
{
tmprows <- createEffects(effectGroup, xName=newname, yName=yName,
name=effects$name[i][1],
groupName=effects$groupName[i][1],
group=effects$group[i][1],
netType=effects$netType[i][1])
tmprows <- tmprows[tmprows$shortName==shortName &
tmprows$type %in% effects$type[i], ]
tmprows$fix <- fix
tmprows$include <- include
tmprows$effectNumber <- max(newEffects$effectNumber) +
(1:nrow(tmprows))
tmprows$timeDummy <- timeDummy
rownames(tmprows) <- paste(newname, effects$type[i], sep=".")
newEffects <- rbind(newEffects, tmprows)
newEffects
}
##@addVarCovar internal sienaTimeFix add a varying covariate
addVarCovar <- function(data, base, nodeSet, newname, dataSub)
{
base <- varCovar(base, nodeSet=nodeSet)
base <- addAttributes.varCovar(base, name=dname)
data[[dataSub]]$vCovars <- c(data[[dataSub]]$vCovars, list(base))
n <- length(data[[dataSub]]$vCovars)
names(data[[dataSub]]$vCovars)[n] <- dname
data
}
if (getDocumentation)
{
return(getInternals())
}
nGroups <- max(unique(effects$group))
# note that period is a character and must be changed to numeric
groupPeriods <- as.numeric(tapply(as.numeric(effects$period), effects$group,
max, na.rm=TRUE)) + 1
intervals <- rep(1:length(groupPeriods), groupPeriods)
localPeriodNos <- unlist(sapply(groupPeriods, function(x)1:x))
nGroups <- length(groupPeriods)
v1 <- do.call(c, lapply(groupPeriods, function(x)1:x))
v2 <- rep(cumsum(c(0, groupPeriods))[1:nGroups], groupPeriods)
v3 <- v1 + v2
periodNos <- v3[!v3 %in% cumsum(groupPeriods)]
if (!is.null(data))
{
atts <- attributes(data)
}
if ((length(periodNos) < 2) && any(effects$timeDummy != ","))
{
warning("Time dummies not relevant with only 2 periods")
effects$timeDummy <- ","
}
structuralRate <- effects$type == "rate" & effects$rateType %in%
"structural"
if (any(effects$timeDummy[structuralRate] != ","))
{
warning("Time dummy effects are not implemented",
" for structural rate effects.")
effects$timeDummy[structuralRate] <- ","
}
# JAL: Implementing these 20-FEB-2011 in RSienaTest
# TODO: Behavioral interactions need to be implemented for these to work.
# Once they are, we can comment lines 575-577 and 717-720 out.
## behaviorNonRateX <- effects$netType =="behavior" & effects$type != "rate"
## if (any(effects$timeDummy[behaviorNonRateX] != ","))
## {
## warning("Time dummy effects are not implemented",
## " for behavior effects of type eval or endow.")
## effects$timeDummy[behaviorNonRateX] <- ","
## }
if (all(effects$timeDummy == ",") )
{
## No time dummy interactions to add, so kick the inputs back.
return(list(effects=effects, data=data))
}
else
{
effects$timeDummy[effects$timeDummy=="all"] <-
paste(periodNos[-1], collapse = ",")
alreadyDummied <- grepl("isDummy", effects$timeDummy)
if (length(alreadyDummied) > 0)
{
## Just remove those effects that have already been dummied so as to
## not messy things up. The assumption is that the user will retain
## all of the previous dummied effects within the column.
effects <- effects[!alreadyDummied, ]
}
timesd <- strsplit(effects$timeDummy, ",| ")
opar <- options(warn=-1)
vals <- sapply(timesd, function(x) any(x=="" | !is.na(as.numeric(x))))
if (any(!vals))
{
stop("Invalid input in timeDummy column", ": ",
paste(effects$timeDummy[!vals], collapse=" and "))
}
options(opar)
timesd <- lapply(timesd, function(x)as.numeric(x[x %in% periodNos]))
dummiedEffects <- sapply(timesd, function(x)length(x) > 0)
baseType <- list(number=which(dummiedEffects),
type=effects$netType[dummiedEffects],
name=effects$name[dummiedEffects])
rateXDummies <- effects$shortName == "RateX" & dummiedEffects
newEffects <- effects
##@getCovar internal sienaTimeFix find the covariate
getCovar <- function(cdvind, vdvind, bdvind, p)
{
if (!is.na(cdvind))
{
tmp <- data[[dataSub]]$cCovars[[cdvind]]
val <- tmp
}
else if (!is.na(vdvind))
{
tmp <- data[[dataSub]]$vCovars[[vdvind]]
val <- tmp[, p]
}
else
{
tmp <- data[[dataSub]]$depvars[[bdvind]]
val <- tmp[, , p]
}
list(covar=tmp, val=val)
}
for (i in which(rateXDummies))
{
effect <- newEffects[i, ]
## Figure out the base values:
if (!is.null(data))
{
cdvind <- match(effect$interaction1, atts$cCovars)
vdvind <- match(effect$interaction1, atts$vCovars)
bdvind <- match(effect$interaction1, atts$netnames)
if (any((is.na(cdvind) & is.na(vdvind) &
(is.na(bdvind))))) #|| atts$types[bdvind] != "behavior"))
{
stop("Having trouble finding the covariate for your rate ",
"effect. Please contact the developers.")
}
}
for (p in timesd[[i]])
{
dname <- paste(effect$interaction1, "Dummy", p, ":",
effect$name, sep="")
dataSub <- intervals[p]
pp <- localPeriodNos[p]
if (!is.null(data))
{
nPer <- attr(data[[dataSub]]$depvars[[effect$name]],
"netdims")[3]
covar <- getCovar(cdvind, vdvind, bdvind, pp)
## add a new varCovar:
nodeSet <- attr(covar$covar, "nodeSet")
base <- matrix(0, nrow=length(covar$val), ncol=nPer - 1)
base[, pp] <- covar$val
data <- addVarCovar(data, base, nodeSet, dname, dataSub)
## also add to other data objects as all must have same set
for (ii in 1:length(data))
{
if (ii != dataSub)
{
nActors <-
length(data[[ii]]$nodeSets[[nodeSet]])
nPer <-
attr(data[[ii]]$depvars[[effect$name]],
"netdims")[3]
base <- matrix(0, nrow=nActors, ncol=nPer - 1)
data <- addVarCovar(data, base, nodeSet, dname, ii)
}
}
}
## add the rate effect row:
newEffects <-
addEffect(newEffects, i, dname, "covarNonSymmetricRate",
'RateX', paste('isDummy', p, i, sep=','))
}
}
## now the non rateX effects
## first add dummies for all requested periods for each
## dependent variable
## then construct the interaction effects with the other effects
## the dummies are fixed unless we need dummies for the density effect
## find which periods we need
use <- sapply(timesd, length) > 0 & !rateXDummies
timeslist <- split(timesd[use], effects$name[use])
timeslist <- lapply(timeslist, function(x)sort(unique(unlist(x))))
timesTypelist <- split(timesd[use],
list(effects$name[use], effects$type[use]))
timesTypelist <- lapply(timesTypelist,
function(x)unique(unlist(x)))
types <- unique(effects$type[use ])
for (depvar in names(timeslist))
{
##if (!is.null(data) && atts$types[[depvar]] == "behavior")
##{
## stop ("Function is not specified for behavior effects")
##}
for (p in timeslist[[depvar]])
{
## create the dummy covariate
pp <- localPeriodNos[p]
dname <- paste("Dummy", p, ":", depvar, sep="")
dataSub <- intervals[p]
if (!is.null(data))
{
dims <- attr(data[[dataSub]]$depvars[[depvar]],
"netdims")
nodeSet <- attr(data[[dataSub]]$depvars[[depvar]],
"nodeSet")[1]
nActors <- dims[1]
nPer <- dims[3]
tmp <- matrix(0, nActors, nPer - 1)
tmp[, pp] <- 1
data <- addVarCovar(data, tmp, nodeSet, dname, dataSub)
## also add to other data objects as all must be the same
for (i in 1:length(data))
{
if ( i != dataSub)
{
dims <- attr(data[[i]]$depvars[[depvar]],
"netdims")
nodeSet <-
attr(data[[i]]$depvars[[depvar]],
"nodeSet")[1]
nActors <- dims[1]
nPer <- dims[3]
base <- matrix(0, nActors, nPer - 1)
data <- addVarCovar(data, base, nodeSet, dname, i)
}
}
}
## This is where the behaviors are treated differently
if (baseType$type[baseType$name==depvar][1] == "behavior")
{
##find the linear effect rows for this depvar
i <- which(effects$name == depvar &
effects$shortName=="linear" &
effects$type %in% types)
if (length(i) == 0)
{
stop("Cannot find 'effect from' effect(s) for ", depvar,
" ", types)
}
## establish whether we want the effFroms
## included, fixed or not
fix <- sapply(timesd[i], function(x)!p %in% x)
typesNames <- paste(depvar, effects$type[i], sep=".")
includeTypes <- !sapply(timesTypelist[typesNames], is.null)
## add one or more effects (depending on types requested)
newEffects <- addEffect(newEffects, i, dname,
"covarBehaviorObjective", "effFrom",
paste('isDummy', p,
effects$effectNumber[i],
sep=','),
fix=fix, include=includeTypes,
yName=dname)
}
else
{
##find the density rows for this depvar
i <- which(effects$name == depvar &
effects$shortName=="density" &
effects$type %in% types)
if (length(i) == 0)
{
stop("Cannot find density effect(s) for ", depvar,
" ", types)
}
## establish whether we want the egoXs included,
## fixed or not
fix <- sapply(timesd[i], function(x)!p %in% x)
typesNames <- paste(depvar, effects$type[i], sep=".")
includeTypes <- !sapply(timesTypelist[typesNames], is.null)
## add one or more effects (depending on types requested)
newEffects <-
addEffect(newEffects, i, dname,
"covarNonSymmetricObjective", "egoX",
paste('isDummy', p,
effects$effectNumber[i], sep=','),
fix=fix, include=includeTypes)
}
}
## now the interactions for any dummied effects for this
## dependent variable
for (j in
which(dummiedEffects & !rateXDummies &
effects$name==depvar &
effects$shortName != "density" &
effects$shortName != "linear"))
{
for (p in timesd[[j]])
{
dname <- paste("Dummy", p, ":", depvar, sep="")
effect <- effects[j, ]
if (baseType$type[baseType$name==depvar][1] == "behavior") {
newEffects <-
includeInteraction(newEffects,
effect$shortName, "effFrom",
character=TRUE,
type=effect$type,
interaction1= c(effect$interaction1,
dname),
interaction2=effect$interaction2,
parameter=effect$parameter,
name=depvar, verbose=FALSE)
} else {
newEffects <-
includeInteraction(newEffects,
effect$shortName, "egoX",
character=TRUE,
type=effect$type,
interaction1= c(effect$interaction1,
dname),
interaction2=effect$interaction2,
parameter=effect$parameter,
name=depvar, verbose=FALSE)
}
## find the row altered
newrow <- newEffects$effect1 == j &
newEffects$effect2 ==
newEffects[paste(dname, effect$type, sep="."),
"effectNumber"]
newEffects$timeDummy[newrow] <-
paste('isDummy', p, j, sep=',')
}
}
}
## reorder the effects so Dummies come before interactions and
## rates where they are expected
e0 <- split(newEffects, newEffects$name)
e0 <- lapply(e0, function(x)
{
e1 <- !grepl("unspecified interaction effect",
x$effectName)
e2 <- grepl("unspecified interaction effect",
x$effectName)
e3 <- x$type == "rate"
rbind(x[e3, ], x[e1 & ! e3, ], x[e2 & ! e3, ])
})
## sort into order of names in original effect
effNames <- unique(effects$name)
order1 <- match(effNames, names(e0))
e0 <- e0[order1]
newEffects <- do.call(rbind, e0)
## now reconstruct the group object to make sure all attributes are
## correct after adding covariates
if (!is.null(data))
{
data <- sienaGroupCreate(data, singleOK=TRUE)
}
list(effects=newEffects, 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
if (any(use))
{
myeff[use, "timeDummy"] <- timeDummy
myeff[use, "include"] <- include
print.data.frame(myeff[use,])
}
else
{
warning("No effects could be found with this specification.")
}
myeff
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.