[Rsiena-commits] r93 - in pkg/RSienaTest: R doc inst/examples man src/model src/model/ml tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jun 4 14:25:08 CEST 2010
Author: ripleyrm
Date: 2010-06-04 14:25:08 +0200 (Fri, 04 Jun 2010)
New Revision: 93
Added:
pkg/RSienaTest/doc/Siena_algorithms4.tex
pkg/RSienaTest/inst/examples/algorithms.r
pkg/RSienaTest/inst/examples/runalg.r
pkg/RSienaTest/man/includeTimeDummy.Rd
pkg/RSienaTest/man/plot.sienaTimeTest.Rd
Modified:
pkg/RSienaTest/R/effectsMethods.r
pkg/RSienaTest/R/globals.r
pkg/RSienaTest/R/maxlikec.r
pkg/RSienaTest/R/phase1.r
pkg/RSienaTest/R/phase2.r
pkg/RSienaTest/R/phase3.r
pkg/RSienaTest/R/sienaDataCreate.r
pkg/RSienaTest/R/sienaTimeTest.r
pkg/RSienaTest/R/sienaeffects.r
pkg/RSienaTest/R/sienaprint.r
pkg/RSienaTest/R/simstatsc.r
pkg/RSienaTest/doc/RSienaDeveloper.tex
pkg/RSienaTest/man/includeEffects.Rd
pkg/RSienaTest/man/includeInteraction.Rd
pkg/RSienaTest/man/setEffect.Rd
pkg/RSienaTest/man/sienaTimeTest.Rd
pkg/RSienaTest/src/model/EpochSimulation.cpp
pkg/RSienaTest/src/model/EpochSimulation.h
pkg/RSienaTest/src/model/Model.cpp
pkg/RSienaTest/src/model/Model.h
pkg/RSienaTest/src/model/ml/BehaviorChange.cpp
pkg/RSienaTest/src/model/ml/BehaviorChange.h
pkg/RSienaTest/src/model/ml/Chain.cpp
pkg/RSienaTest/src/model/ml/Chain.h
pkg/RSienaTest/src/model/ml/MLSimulation.cpp
pkg/RSienaTest/src/model/ml/MiniStep.cpp
pkg/RSienaTest/src/model/ml/MiniStep.h
pkg/RSienaTest/src/model/ml/NetworkChange.cpp
pkg/RSienaTest/src/model/ml/NetworkChange.h
pkg/RSienaTest/tests/parallel.R
pkg/RSienaTest/tests/parallel.Rout.save
Log:
Algorithms, maximum likelihood, help pages
Modified: pkg/RSienaTest/R/effectsMethods.r
===================================================================
--- pkg/RSienaTest/R/effectsMethods.r 2010-06-04 11:51:56 UTC (rev 92)
+++ pkg/RSienaTest/R/effectsMethods.r 2010-06-04 12:25:08 UTC (rev 93)
@@ -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/RSienaTest/R/globals.r
===================================================================
--- pkg/RSienaTest/R/globals.r 2010-06-04 11:51:56 UTC (rev 92)
+++ pkg/RSienaTest/R/globals.r 2010-06-04 12:25:08 UTC (rev 93)
@@ -23,6 +23,7 @@
x <- x
beverbose <- verbose
besilent <- silent
+ noReportFile <- FALSE
function(txt, dest, fill=FALSE, sep=" ", hdest, openfiles=FALSE,
closefiles=FALSE, type=c("a", "w", "n"), projname="Siena" ,
verbose=FALSE, silent=FALSE)
@@ -40,6 +41,10 @@
{
x$outf <<- file(paste(projname, ".out", sep=""), open="a")
}
+ else if (type == "n")
+ {
+ noReportFile <<- TRUE
+ }
}
else if (closefiles)
@@ -69,7 +74,10 @@
}
else
{
- cat(txt, file = x[[hdest]], fill = fill, sep = sep)
+ if (!noReportFile)
+ {
+ cat(txt, file = x[[hdest]], fill = fill, sep = sep)
+ }
}
}
else
@@ -85,12 +93,18 @@
{
if (is.null(x[[deparse(substitute(dest))]]))
{
- cat(txt, fill=fill, sep=sep)
+ if (!besilent)
+ {
+ cat(txt, fill=fill, sep=sep)
+ }
}
else
{
- cat(txt, file=x[[deparse(substitute(dest))]],
- fill=fill, sep=sep)
+ if (!noReportFile)
+ {
+ cat(txt, file=x[[deparse(substitute(dest))]],
+ fill=fill, sep=sep)
+ }
}
}
}
Modified: pkg/RSienaTest/R/maxlikec.r
===================================================================
--- pkg/RSienaTest/R/maxlikec.r 2010-06-04 11:51:56 UTC (rev 92)
+++ pkg/RSienaTest/R/maxlikec.r 2010-06-04 12:25:08 UTC (rev 93)
@@ -33,7 +33,10 @@
f$pModel <- NULL
f$pData <- NULL
FRANstore(NULL) ## clear the stored object
- PrintReport(z, x)
+ if (is.null(z$print) || z$print)
+ {
+ PrintReport(z, x)
+ }
if (sum(z$test))
{
z$fra <- colMeans(z$sf, na.rm=TRUE)
@@ -41,7 +44,10 @@
z <- c(z, ans)
TestOutput(z, x)
}
- dimnames(z$dfra)[[1]] <- as.list(z$requestedEffects$shortName)
+ if (!is.null(z$dfra))
+ {
+ dimnames(z$dfra)[[1]] <- as.list(z$requestedEffects$shortName)
+ }
return(z)
}
######################################################################
@@ -50,14 +56,14 @@
## retrieve stored information
f <- FRANstore()
## browser()
- if (z$Phase == 2)
- {
- returnDeps <- FALSE
- }
- else
- {
- returnDeps <- f$returnDeps
- }
+ ##if (z$Phase == 2)
+ ##{
+ ## returnDeps <- FALSE
+ ##}
+ ##else
+ ##{
+ returnDeps <- z$returnDeps
+ ##}
## create a grid of periods with group names in case want to parallelize
## using this
groupPeriods <- attr(f, "groupPeriods")
@@ -68,10 +74,10 @@
## we are not
if (z$int2==1 || nrow(callGrid) == 1)
{
- ## for now!
- ans <- .Call('mlPeriod', PACKAGE=pkgname, z$Deriv, f$pData,
- f$pModel, f$myeffects, f$pMLSimulation, z$theta,
- returnDeps, 1, 1)
+ ## for now!
+ ans <- .Call('mlPeriod', PACKAGE=pkgname, z$Deriv, f$pData,
+ f$pModel, f$myeffects, f$pMLSimulation, z$theta,
+ returnDeps, 1, 1, z$nrunMH)
}
else
{
@@ -131,9 +137,9 @@
dffraw <- ans[[7]]
dff <- matrix(0, nrow=z$pp, ncol=z$pp)
start <- 1
+ rawsub <- 1
for (i in 1:length(f$myeffects))
{
- rawsub <- 1
rows <- nrow(f$myeffects[[i]])
nRates <- sum(f$myeffects[[i]]$type == 'rate')
nonRates <- rows - nRates
@@ -145,8 +151,10 @@
start : (start + nonRates - 1)] <-
dffraw[rawsub:(rawsub + nonRates * nonRates - 1)]
start <- start + nonRates
+ rawsub <- rawsub + nonRates * nonRates
}
- dff[upper.tri(dff)] <- dff[lower.tri(dff)]
+ dff <- dff + t(dff)
+ diag(dff) <- diag(dff) / 2
dff <- -dff
}
else
@@ -156,7 +164,18 @@
# browser()
list(fra = fra, ntim0 = NULL, feasible = TRUE, OK = TRUE,
- sims=sims, dff = dff, chain = ans[[6]], accepts=ans[[8]],
+ sims=sims, dff = dff, chain = list(ans[[6]]), accepts=ans[[8]],
rejects= ans[[9]])
}
+dist2full<-function(dis) {
+
+ n<-attr(dis,"Size")
+
+ full<-matrix(0,n,n)
+
+ full[lower.tri(full)]<-dis
+
+ full+t(full)
+
+}
Modified: pkg/RSienaTest/R/phase1.r
===================================================================
--- pkg/RSienaTest/R/phase1.r 2010-06-04 11:51:56 UTC (rev 92)
+++ pkg/RSienaTest/R/phase1.r 2010-06-04 12:25:08 UTC (rev 93)
@@ -568,5 +568,7 @@
zsmall$cconditional <- z$cconditional
zsmall$condvar <- z$condvar
zsmall$pp <- z$pp
+ zsmall$nrunMH <- z$nrunMH
+ zsmall$returnDeps <- z$returnDeps
zsmall
}
Modified: pkg/RSienaTest/R/phase2.r
===================================================================
--- pkg/RSienaTest/R/phase2.r 2010-06-04 11:51:56 UTC (rev 92)
+++ pkg/RSienaTest/R/phase2.r 2010-06-04 12:25:08 UTC (rev 93)
@@ -63,7 +63,7 @@
z
}
##@proc2subphase siena07 Do one subphase of phase 2
-proc2subphase<- function(z, x, subphase, ...)
+proc2subphase <- function(z, x, subphase, useAverage=TRUE, ...)
{
## init subphase of phase 2
z <- AnnouncePhase(z, x, subphase)
@@ -87,14 +87,15 @@
z$time1 <- proc.time()[3]
z$thav <- z$theta
z$thavn <- 1
- ## cat(z$thav, z$theta, '\n')
+ ## cat(z$thav, z$theta, '\n')
z$prod0 <- rep(0, z$pp)
z$prod1 <- rep(0, z$pp)
## ###############################################
## do the iterations for this repeat of this subphase
## ##############################################
+ ##z <- doIterationsCopy(z, x, subphase, ...)
z <- doIterations(z, x, subphase, ...)
- ## if (z$nit == 50) browser()
+ ## if (z$nit == 50) browser()
if (!z$OK || UserInterruptFlag() || UserRestartFlag() ||
EarlyEndPhase2Flag())
{
@@ -155,7 +156,20 @@
' = ', format(subphaseTime, nsmall=4, digits=4),
'\n', sep=''), lf)
}
- z$theta <- z$thav / z$thavn #(z$nit + 1)
+ if (useAverage)
+ {
+ z$theta <- z$thav / z$thavn #(z$nit + 1)
+ }
+ else
+ {
+ ## use regession
+ cat(z$thav / z$thavn, '\n') #(z$nit + 1)
+ mylm <- lm(z$sf[1:z$nit, ] ~ z$thetaStore[1:z$nit, ])
+ coefs <- coef(mylm)[-1, ]
+ newvals <- solve(t(coefs), - mylm$coef[1, ])
+ z$theta <- newvals
+ cat(z$theta, '\n')
+ }
DisplayThetaAutocor(z)
## cat('it',z$nit,'\n')
##recalculate autocor using -1 instead of -2 as error
@@ -192,6 +206,7 @@
ac <- 0
xsmall <- NULL
zsmall <- makeZsmall(z)
+ z$returnDeps <- FALSE
repeat
{
z$n <- z$n+1
@@ -340,3 +355,446 @@
}
z
}
+##@doIterationsCopy siena07 Do all iterations for 1 repeat of 1 subphase of phase 2
+doIterationsCopy <- function(z, x, subphase, numberIterations=10,
+ ...)
+{
+ ## designed to try things out!
+ if (z$cconditional)
+ {
+ stop("cannot use this routine with conditional estimation")
+ }
+ z$nit <- 0
+ ac <- 0
+ xsmall <- NULL
+ zsmall <- makeZsmall(z)
+ z <- setUpPhase2Storage(z, numberIterations)
+ repeat
+ {
+ z$n <- z$n+1
+ z$nit <- z$nit + 1
+ if (subphase == 1 && z$nit > 1)
+ z$time1 <- proc.time()[[3]]
+ if (subphase == 1 && z$nit > 10)
+ {
+ time1 <- proc.time()[[3]] - z$time1
+ if (time1 > 1e-5)
+ {
+ z$writefreq <- max(1, round(20.0 / time1))
+ }
+ else
+ {
+ z$writefreq <- 20
+ }
+ ## if (is.batch())
+ ## {
+ ## z$writefreq <- z$writefreq * 2 ##compensation for it
+ ## ## running faster with no tcl/tk
+ ## }
+ z$writefreq <- roundfreq(z$writefreq)
+ z$writefreq <- 1
+ }
+ if ((z$nit <= 10) || (z$nit %% z$writefreq ==0))
+ {
+ DisplayIteration(z)
+ if (is.batch())
+ {
+ val <- getProgressBar(z$pb)
+ increment <- ifelse(z$nit <= 10, 1, z$writefreq)
+ Report(paste('Phase ', z$Phase, ' Subphase ', subphase,
+ ' Iteration ', z$nit,' Progress: ',
+ round((increment + val) /
+ z$pb$pbmax * 100),
+ '%\n', sep = ''))
+ z$pb <- setProgressBar(z$pb, val + increment)
+ }
+ else
+ {
+ if (z$nit > 1)
+ {
+ DisplayDeviations(z, fra)
+ }
+ if (z$nit %% z$writefreq == 0)
+ {
+ val <- getProgressBar(z$pb)
+ z$pb <-setProgressBar(z$pb, val + z$writefreq)
+ }
+ }
+ }
+ zsmall$nit <- z$nit
+
+ zsmall$addChainToStore <- FALSE
+ zsmall$needChangeContributions <- FALSE
+
+ if (z$int == 1) ## not using a cluster
+ {
+ fra <- z$targets - z$targets
+ for (i in 1:numberIterations)
+ {
+ zz <- x$FRAN(zsmall, xsmall, returnDeps=TRUE)
+ ## browser()
+ fra0 <- colSums(zz$fra) - z$targets
+ fra <- fra + fra0
+ if (!zz$OK)
+ {
+ z$OK <- zz$OK
+ break
+ }
+ z <- storeChainsAndFra(z, zz, i, fra0)
+ }
+ if (!z$OK)
+ {
+ break
+ }
+ fra <- fra / numberIterations
+ z <- calculateLikelihoods(z)
+ }
+ else
+ {
+ zz <- clusterCall(z$cl, usesim, zsmall, xsmall)
+ fra <- sapply(zz, function(x) colSums(x$fra) - z$targets)
+ dim(fra) <- c(z$pp, z$int)
+ fra <- rowMeans(fra)
+ zz$OK <- sapply(zz, function(x) x$OK)
+ if (!all(zz$OK))
+ {
+ z$OK <- FALSE
+ break
+ }
+ }
+ if (x$maxlike)
+ {
+ z$phase2fras[subphase, ,z$nit] <- fra
+ ## z$rejectprops[subphase, , z$nit] <- zz$rejectprop
+ }
+ if (z$nit %% 2 == 1)
+ {
+ prev.fra <- fra
+ }
+ else
+ {
+ z$prod0 <- z$prod0 + fra * fra
+ z$prod1 <- z$prod1 + fra * prev.fra
+ ac <- ifelse (z$prod0 > 1e-12, z$prod1 / z$prod0, -2)
+ z$maxacor <- max(-99, ac[!z$fixed]) ##note -2 > -99
+ z$minacor <- min(1, ac[(!z$fixed) & ac > -1.0])
+ z$ac <- ac
+ if (z$nit %% z$writefreq == 0)
+ {
+ DisplayThetaAutocor(z)
+ }
+ }
+ z <- doChangeStep(z, x, fra)
+ zsmall$theta <- z$theta
+
+ if (x$maxlike && !is.null(x$moreUpdates) && x$moreUpdates > 0)
+ {
+ z <- doMoreUpdates(z, x, x$moreUpdates * subphase)
+ zsmall$theta <- z$theta
+ }
+ ## importance sampling steps
+ importSub <- 1
+ repeat
+ {
+ ## browser()
+ z <- predictOutcomes(z)
+ #cat(z$varWeights,'\n')
+ if (is.na(z$varWeights) ||
+ z$varWeights > var(c(rep(0, numberIterations-1), 1))/10)
+ {
+ break
+ }
+ z <- doChangeStep(z, x, z$predictedStats)
+ zsmall$theta <- z$theta
+ importSub <- importSub + 1
+ if (importSub > 25)
+ {
+ break
+ }
+ }
+ if (zsmall$addChainToStore)
+ {
+ clearStoredChains()
+ }
+ ##check for user interrupt
+ ## browser()
+ CheckBreaks()
+ if (UserInterruptFlag() || UserRestartFlag() || EarlyEndPhase2Flag())
+ {
+ break
+ }
+ ## do we stop?
+ if ( (z$nit >= z$n2min && z$maxacor < 1e-10)
+ || (z$nit >= z$n2max) || (z$nit >= 50 && z$minacor < -0.8 &&
+ z$repeatsubphase < z$maxrepeatsubphase))
+ {
+ break
+ }
+ }
+ z
+}
+
+##@setUpPhase2Storage siena07 create stores and values in nonstandard Phase2
+setUpPhase2Storage <- function(z, niter)
+{
+ nGroup <- z$f$nGroup
+ z$chains <- lapply(1:nGroup, function(x)
+ lapply(1:z$groupPeriods[x], function(y)
+ vector("list", niter)))
+ z$lik0 <- rep(0, niter * sum(z$groupPeriods - 1))
+ z$iterFra <- matrix(0, ncol=z$pp, nrow=niter)
+ z$thetaStore <- matrix(0, ncol=z$pp, nrow=z$n2max)
+ z$sf <- matrix(0, ncol=z$pp, nrow=z$n2max)
+ z
+}
+
+##@storeChainaAndFra siena07 update step for storage in non-standard phase 2
+storeChainsAndFra <- function(z, zz, i, fra)
+{
+ for (ii in 1:z$nGroup)
+ {
+ for (jj in 1:z$groupPeriods[ii])
+ {
+ z$chains[[ii]][[jj]][[i]] <- zz$chain[[ii]][[jj]]
+ }
+ }
+ z$iterFra[i, ] <- fra
+ z$thetaStore[z$nit, ] <- z$theta
+ z$sf[z$nit, ] <- fra
+ z
+}
+
+##@calculateLikelihoods siena07 for use in non-standard phase 2
+calculateLikelihoods <- function(z)
+{
+ storeSub <- 1
+ for (ii in 1:z$nGroup)
+ {
+ for (jj in 1:z$groupPeriods[ii])
+ {
+ chains <- z$chains[[ii]][[jj]]
+ for (chain in chains)
+ {
+ if (length(chain) == 0)
+ {
+ stop("no events with this theta")
+ }
+ z$lik0[storeSub] <-
+ getLikelihoodPhase2(chain, z$nactors[[ii]],
+ z$theta[z$rateParameterPosition[[ii]][[jj]]])
+ storeSub <- storeSub + 1
+ }
+ }
+ }
+ z
+}
+
+##@getLikelihoodPhase2 siena07 for use in non-standard phase 2
+getLikelihoodPhase2 <- function(chain, nactors, lambda)
+{
+ loglik <- 0
+ ncvals <- sapply(chain, function(x)x[[3]])
+ nc <- nactors
+ nc[] <- 0
+ ncvals <- table(ncvals)
+ nc[names(ncvals)] <- ncvals
+ logChoiceProb <- sapply(chain, function(x)x[[9]])
+ logOptionSetProb <- sapply(chain, function(x)x[[8]])
+ loglik <- sum(logChoiceProb) # + sum(logOptionSetProb)
+ #print(sum(logOptionSetProb))
+ loglik <- loglik - sum(nactors * lambda) + sum(nc * log(lambda))
+ loglik
+}
+
+##@predictOutcomes siena07 for use in non-standard phase 2
+predictOutcomes <- function(z)
+{
+ ## now the likelihood for the new theta
+ ps <- getProbabilitiesPhase2(z$chain, z$theta, z$nactors,
+ z$rateParameterPosition, z$cl)$lik
+ ps <- exp(unlist(ps) - z$lik0)
+ ps <- ps / sum(ps)
+ z$predictedStats <- apply(z$iterFra, 2, function(x) sum(x * ps))
+ z$varWeights <- var(ps)
+ z
+}
+##@doChangeStep siena07 for use in non-standard phase 2, or outside siena07
+doChangeStep <- function(z, x, fra)
+{
+ ## limit change. Reporting is delayed to
+ ## end of phase.
+ ## browser()
+ if (x$diag)## !maxlike at present
+ {
+ maxrat<- max(ifelse(z$sd, abs(fra)/ z$sd, 1.0))#### check this
+ if (maxrat > x$maxmaxrat)
+ {
+ maxrat <- x$maxmaxrat / maxrat
+ z$truncated[z$nit] <- TRUE
+ }
+ else
+ maxrat <- 1.0
+ fchange<- z$gain * fra * maxrat / diag(z$dfra)
+ }
+ else
+ {
+ fchange <- as.vector(z$gain * fra %*% z$dinv)
+ }
+ ## browser()
+ fchange[z$fixed] <- 0.0
+ ## check positivity restriction
+ if (!is.null(z$nit))
+ {
+ z$positivized[z$nit, ] <- z$posj & (fchange > z$theta)
+ }
+ fchange <- ifelse(z$posj & (fchange > z$theta), z$theta * 0.5, fchange)
+
+ z$theta <- z$theta - fchange
+ if (!is.null(z$thav))
+ {
+ z$thav <- z$thav + z$theta
+ z$thavn <- z$thavn + 1
+ }
+ z
+}
+
+##@clearStoredChains algorithms Clears storage used for chains in C.
+clearStoredChains <- function()
+{
+ f <- FRANstore()
+ .Call("clearStoredChains", PACKAGE=pkgname, f$pModel)
+}
+
+##@getProbabilitiesPhase2 algorithms Recalculates change contributions
+getProbabilitiesPhase2 <- function(chain, theta, nactors, rateParameterPosition,
+ cl=NULL)
+{
+ f <- RSiena:::FRANstore()#
+ getScores <- FALSE
+ ##print(theta)
+ for (i in 1:length(chain)) # group
+ {
+ for (j in 1:length(chain[[i]])) # period
+ {
+ if (!is.null(cl) )
+ {
+ tmp <- parLapply(cl, chain[[i]][[j]],
+ function(x, i, j, theta, getScores, k, n)
+ {
+ f <- RSiena:::FRANstore()
+ resp <- .Call("getChainProbabilitiesList",
+ PACKAGE = pkgname, x,
+ f$pData, f$pModel, as.integer(i),
+ as.integer(j), f$myeffects,
+ theta, getScores)
+ lik <-
+ getLikelihoodPhase2(resp[[1]], nactors[[i]],
+ theta[k])
+ if (getScores)
+ {
+ sc <- resp[[2]]
+ }
+ else
+ {
+ sc <- NULL
+ }
+ list(lik=lik, sc=sc)
+ }, i=i, j=j, theta=theta, getScores=getScores,
+ k=rateParameterPosition[[i]][[j]],
+ n=nactors[[i]]
+ )
+ lik <- sapply(tmp, function(x)x[[1]])
+ if (getScores)
+ {
+ sc <- t(sapply(tmp, function(x)x[[2]]))
+ }
+ else
+ {
+ sc <- NULL
+ }
+ }
+ else
+ {
+ lik <- rep(0, length(chain[[i]][[j]]))
+ sc <- matrix(0, nrow=length(chain[[i]][[j]]),
+ ncol=length(theta))
+ for (k in 1:length(chain[[i]][[j]])) # one chain
+ {
+ resp <- .Call("getChainProbabilitiesList",
+ PACKAGE = pkgname,
+ chain[[i]][[j]][[k]],
+ f$pData, f$pModel, as.integer(i),
+ as.integer(j), f$myeffects, theta,
+ getScores)
+
+ lik[k] <-
+ getLikelihoodPhase2(resp[[1]], nactors[[i]],
+ theta[rateParameterPosition[[i]][[j]]])
+ if (getScores)
+ {
+ sc[k,] <- resp[[2]]
+ }
+ }
+ }
+ }
+ }
+ list(lik=lik, sc=sc)
+}
+##@initForAlgorithms siena07 stores values for use in nonstandard Phase2
+initForAlgorithms <- function(z)
+{
+ if (z$cconditional)
+ {
+ return(z)
+ }
+ nGroup <- z$f$nGroup
+ z$nDependentVariables <- length(z$f$depNames)
+ atts <- attributes(z$f)
+ z$groupPeriods <- atts$groupPeriods - 1
+ netnames <- names(z$f[[1]]$depvars)
+ z$nactors <- lapply(1:nGroup, function(i, periods, data)
+ {
+ tmp <- sapply(data[[i]]$depvars, function(x)
+ dim(x)[1])
+ tmp <- tmp[match(netnames, names(data[[i]]$depvars))]
+ tmp
+ }, periods=z$groupPeriods, data=z$f[1:nGroup]
+ )
+ z$rateParameterPosition <-
+ lapply(1:nGroup, function(i, periods, data)
+ {
+ lapply(1:periods[i], function(j)
+ {
+ rateEffects <-
+ z$effects[z$effects$basicRate &
+ z$effects$period == j &
+ z$effects$group == i,]
+ rateEffects <-
+ rateEffects[match(netnames,
+ rateEffects$name), ]
+ tmp <- as.numeric(row.names(rateEffects))
+ names(tmp) <- netnames
+ tmp
+ }
+ )
+ }, periods=z$groupPeriods, data=z$f[1:nGroup]
+ )
+ z$evalParameterPosition <-
+ lapply(netnames, function(x)
+ {
+ as.numeric(row.names(z$effects[z$effects$name == x &
+ z$effect$type =="eval", ]))
+ }
+ )
+ names(z$evalParameterPosition) <- netnames
+ z$endowParameterPosition <-
+ lapply(netnames, function(x)
+ {
+ as.numeric(row.names(z$effects[z$effects$name == x &
+ z$effect$type =="endow", ]))
+ }
+ )
+ z$basicRate <- z$effects$basicRate
+ z$nGroup <- nGroup
+ z
+}
Modified: pkg/RSienaTest/R/phase3.r
===================================================================
--- pkg/RSienaTest/R/phase3.r 2010-06-04 11:51:56 UTC (rev 92)
+++ pkg/RSienaTest/R/phase3.r 2010-06-04 12:25:08 UTC (rev 93)
@@ -19,6 +19,7 @@
DisplayTheta(z)
z$Phase <- 3
int <- z$int
+ z$returnDeps <- z$returnDepsStored
if (x$checktime) z$ctime <- proc.time()[3]
## fix up iteration numbers if using multiple processors
Modified: pkg/RSienaTest/R/sienaDataCreate.r
===================================================================
--- pkg/RSienaTest/R/sienaDataCreate.r 2010-06-04 11:51:56 UTC (rev 92)
+++ pkg/RSienaTest/R/sienaDataCreate.r 2010-06-04 12:25:08 UTC (rev 93)
@@ -158,39 +158,8 @@
{
nm[fixup] <- dep
}
- ## Josh inserting this code: #####################################
- ## Adding this list flattening routine so that sienaDataCreate can
- ## accept a list of lists. Useful for expanding time dummies, etc.
- ## Also attaches meaningful names to the flattened items.
- listContainsList <- function (x) {
- which(sapply(seq(along=x), function(i) class(x[[i]])=='list'))
- }
- args <- list(...)
- toFlatten <- listContainsList(args)
- if (length(toFlatten)>0) {
- numberOfOldArgs <- length(args[-toFlatten])
- numberOfNewArgs <- sum(sapply(toFlatten, function(i) length(args[[i]])))
- flattened <- vector("list", numberOfNewArgs + numberOfOldArgs)
- newNames <- rep("", numberOfNewArgs + numberOfOldArgs)
- newNames[1:numberOfOldArgs] <- nm[-toFlatten]
- flattened[seq(numberOfOldArgs)] <- args[-toFlatten]
- count <- numberOfOldArgs
- for (i in toFlatten) {
- for (j in seq(along=args[[i]])) {
- count = count + 1
- flattened[[count]] <- args[[i]][[j]]
- newNames[count] <- names(args[[i]][j])
- }
- }
- names(flattened) <- newNames
- narg <- length(flattened)
- nm <- newNames
- dots <- flattened
- } else {
- dots <- list(...)
- }
- ################################################################
-
+ dots <- list(...)
+ names(dots) <- nm
if (any(duplicated(nm)))
{
stop('names must be unique')
@@ -791,6 +760,7 @@
z$dyvCovars <- dyvCovars
z$compositionChange <- compositionChange
z <- checkConstraints(z)
+ z <- covarDist2(z)
class(z) <- 'siena'
z
}
@@ -1803,3 +1773,123 @@
}
ranges
}
+
+##@covarDist2 DataCreate calculate average alter values for covariate wrt each net
+covarDist2 <- function(z)
+{
+ netNames <- names(z$depvars)
+ netTypes <- sapply(z$depvars, function(x)attr(x, "type"))
+ netActorSet <- sapply(z$depvars, function(x)
+ if (attr(x, "type") == "bipartite")
+ {
+ attr(x, "nodeSet")[1]
+ }
+ else
+ {
+ attr(x, "nodeSet")
+ }
+ )
+ ## find the constant covariates
+ for (i in seq(along=z$cCovars))
+ {
+ nodeSet <- attr(z$cCovars[[i]], "nodeSet")
+ use <- (netTypes != "behavior" & netActorSet == nodeSet)
+ simMeans <- namedVector(NA, netNames[use])
+ for (j in seq(along=z$depvars[use]))
+ {
+ simMeans[netNames[use][j]] <-
+ calcCovarDist2(matrix(z$cCovars[[i]], ncol=1),
+ z$depvars[[j]])
+ }
+ attr(z$cCovars[[i]], "simMeans") <- simMeans
+ }
+ for (i in seq(along=z$vCovars))
+ {
+ nodeSet <- attr(z$vCovars[[i]], "nodeSet")
+ use <- (netTypes != "behavior" & netActorSet == nodeSet)
+ simMeans <- namedVector(NA, netNames[use])
+ for (j in seq(along=z$depvars[use]))
+ {
+ simMeans[netNames[use][j]] <-
+ calcCovarDist2(z$vCovars[[i]], z$depvars[[j]])
+ }
+ attr(z$vCovars[[i]], "simMeans") <- simMeans
+
+ }
+ for (i in seq(along=z$depvars))
+ {
+ if (netTypes[i] == "behavior")
+ {
+ beh <- z$depvars[[i]][, 1, ]
+ ## take off the mean NB no structurals yet!
+ beh <- beh - mean(beh, na.rm=TRUE)
+ nodeSet <- netActorSet[i]
+ use <- (netTypes != "behavior" & netActorSet == nodeSet)
+ simMeans <- namedVector(NA, netNames[use])
+ for (j in seq(along=z$depvars[use]))
+ {
+ simMeans[netNames[use][j]] <-
+ calcCovarDist2(beh[, -ncol(beh), drop=FALSE],
+ z$depvars[[j]], rval=range(beh))
+ }
+ attr(z$depvars[[i]], "simMeans") <- simMeans
+ }
+ }
+ z
+}
+##@ calcCovarDist2 DataCreate similarity mean for alter of this depvar
+calcCovarDist2 <- function(covar, depvar, rval=NULL)
+{
+ ## remove final obs from depvars
+ observations <- dim(depvar)[3] - 1
+
+ simTotal <- rep(NA, observations)
+ simCnt <- rep(NA, observations)
+
+ for (i in 1:observations)
+ {
+ if (ncol(covar) == 1) # maybe constant covariate used for each obs
+ {
+ xx <- covar[, 1]
+ }
+ else
+ {
+ xx <- covar[, i]
+ }
+ if (attr(depvar, "sparse"))
+ {
+ dep <- depvar[[i]]
+ }
+ else
+ {
+ dep <- depvar[, , i]
+ }
+ dep[dep %in% c(10, 11)] <- dep[dep %in% c(10, 11)] - 10
+ vi <- apply(dep, 1, function(x)
+ {
+ if (sum(x, na.rm=TRUE) > 0)
+ {
+ nonMissing <- sum(!is.na(xx[!is.na(x) & x > 0]))
+ if (nonMissing > 0)
+ {
+ sum(x * xx, na.rm=TRUE)/ sum(x, na.rm=TRUE)
+ }
+ else
+ {
+ NA
+ }
+ }
+ else
+ {
+ 0
+ }
+
+ }
+ )
+
+ tmp <- rangeAndSimilarity(vi, rval)
+ simTotal[i] <- tmp$simTotal
+ simCnt[i] <- tmp$simCnt
+ }
+ sum(simTotal)/sum(simCnt)
+}
Modified: pkg/RSienaTest/R/sienaTimeTest.r
===================================================================
--- pkg/RSienaTest/R/sienaTimeTest.r 2010-06-04 11:51:56 UTC (rev 92)
+++ pkg/RSienaTest/R/sienaTimeTest.r 2010-06-04 12:25:08 UTC (rev 93)
@@ -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
## *
@@ -51,7 +51,7 @@
# 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",
+ & length(grep("Dummy",
sienaFit$effects[-escreen,]$effectName)) > 0)
if (length(dscreen)==0)
{
@@ -90,7 +90,7 @@
## this information in dummyByEffect -- this is used
## extensively in plot.sienaTimeTest
dummyByEffect[rownames(toTest)==as.numeric(tmp[3]),
- colnames(toTest)==as.numeric(tmp[2])] <-
+ colnames(toTest)==as.numeric(tmp[2])] <-
which(sienaFit$effects[-escreen,]$
effectNumber[-c(rscreen,dscreen)]==i)
}
@@ -638,7 +638,7 @@
"Dummy",p,sep="")
base <- matrix(0,nact,nper-1)
## Figure out the base values:
- dvind <- which(names(data$cCovars) ==
+ dvind <- which(names(data$cCovars) ==
effects$interaction1[effects$effectNumber==i])
## Stick them into the right time spot
base[,p] <- data$cCovars[[dvind]]
@@ -752,10 +752,10 @@
}
##@includeTimeDummy DataCreate
includeTimeDummy <- function(myeff, ..., timeDummy="all", name=myeff$name[1],
- type="eval", interaction1="", interaction2="",
+ type="eval", interaction1="", interaction2="", include=TRUE,
character=FALSE)
{
-
+
if (character)
{
dots <- sapply(list(...), function(x)x)
@@ -782,5 +782,7 @@
myeff$interaction1 == interaction1 &
myeff$interaction2 == interaction2
myeff[use, "timeDummy"] <- timeDummy
+ myeff[use, "include"] <- include
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/rsiena -r 93
More information about the Rsiena-commits
mailing list