Nothing
#/******************************************************************************
# * SIENA: Simulation Investigation for Empirical Network Analysis
# *
# * Web: http://www.stats.ox.ac.uk/~snijders/siena
# *
# * File: initializeFRAN.r
# *
# * Description: This module contains the code for setting up the data in C++
# * and for ML making the initial chain.
# *****************************************************************************/
##@initializeFRAN siena07 reformat data and send to C. get targets.
initializeFRAN <- function(z, x, data, effects, prevAns=NULL, initC,
profileData=FALSE, returnDeps=FALSE,
returnChains=FALSE, byGroup=FALSE,
returnDataFrame=FALSE, byWave=FALSE,
returnLoglik=FALSE, onlyLoglik=FALSE)
{
##@checkNames used in initializeFRAN for siena07
checkNames <- function(xx, types) {
# checks whether names of named vectors are in names of dependent variables.
# returns xx; if is.null(xx), returns named vector
# with variable names and all values 1.
if (('oneMode' %in% types) || ('bipartite' %in% types))
{
thetype <- 'network'
}
else
{
thetype <- 'behavior'
}
if (inherits(data, "sienaGroup"))
{
theVars <- sapply(data[[1]]$depvars, function(x){attr(x,'type') %in% types})
theVarNames <- names(data[[1]]$depvars)[theVars]
}
else
{
theVars <- sapply(data$depvars, function(x){attr(x,'type') %in% types})
theVarNames <- names(data$depvars)[theVars]
}
if ((length(xx) >= 1) && (!(length(xx) == sum(theVars))))
{
cat(deparse(substitute(xx)), ' = (',
paste(names(xx), xx, sep=" = ", collapse="; ") ,') ,\n')
if (sum(theVars)==1)
{
cat('but there is only one dependent variable.\n')
}
else
{
cat('but there are ', sum(theVars), ' dependent variables.\n')
}
cat('Length of',deparse(substitute(xx)),
'should be equal to number of', thetype,'variables.\n')
stop('Invalid algorithm-data combination.')
}
# For this check:
# if there are no names, then (names(xx) %in% names(data$depvars))
# will be logical(0), and all(logical(0)) is TRUE.
if (is.null(xx))
{
wrongName <- FALSE
xx <- rep(1, sum(theVars))
names(xx) <- theVarNames
}
else
{
wrongName <- (!all(names(xx) == theVarNames))
if (is.null(names(xx)))
{
wrongName <- TRUE
}
else
{
if (any(is.na(names(xx))))
{
wrongName <- TRUE
}
}
if (wrongName)
{
cat(deparse(substitute(xx)),
'in the algorithm object x should be a named vector with\n')
cat(' the names of dependent', thetype, 'variables')
cat(' in the data set, in the same order.\n')
cat(' Names in algorithm object: ', names(xx), '\n')
cat(' Names in data set: ', theVarNames, '\n')
stop('Invalid algorithm-data combination.')
}
}
xx
}
# start of initializeFRAN proper
z$effectsName <- deparse(substitute(effects))
## fix up the interface so can call from outside robmon framework
if (is.null(z$FinDiff.method))
{
z$FinDiff.method <- FALSE
}
if (is.null(z$int))
{
z$int <- 1
}
if (is.null(z$int2))
{
z$int2 <- 1
}
if (!initC) ## ie first time round
{
if (!inherits(data,"siena"))
{
stop("not valid siena data object")
}
## check the effects object
defaultEffects <- getEffects(data)
if (is.null(effects))
{
effects <- defaultEffects
}
else
{
## todo check that the effects match the data dependent variables
userlist <- apply(effects[effects$include,], 1, function(x)
paste(x[c("name", "effectName",
"type", "groupName")],
collapse="|"))
deflist <- apply(defaultEffects, 1, function(x)
paste(x[c("name", "effectName",
"type", "groupName")],
collapse="|"))
if (!all(userlist %in% deflist))
{
bad <- which(!(userlist %in% deflist))
print(userlist[bad])
cat("invalid effect requested: see above; \n")
cat("there seems to be a mismatch between data set and effects object.\n")
stop("Perhaps the effects object must be created from scratch.")
}
}
if (!inherits(effects, "data.frame"))
{
stop("effects is not a data.frame")
}
if (x$useStdInits)
{
if (any(effects$effectName != defaultEffects$effectName))
{
stop("Cannot use standard initialisation with a ",
"different effect list")
}
effects$initialValue <- defaultEffects$initialValue
}
# Check that the following attributes have correct names
# and change NULL to default values for x$modelType and x$behModelType
# For symmetric networks, default is changed below from 1 to 2.
checkNames(x$MaxDegree, c('oneMode','bipartite'))
checkNames(x$UniversalOffset, c('oneMode','bipartite'))
x$modelType <- checkNames(x$modelType, c('oneMode','bipartite'))
x$behModelType <- checkNames(x$behModelType, 'behavior')
# The following error will occur if ML estimation is requested
# and there are any impossible changes from structural values
# to different observed values.
if (x$maxlike)
{
if (inherits(data, "sienaGroup"))
{
impossible <- FALSE
zerochange <- FALSE
imp.change <- lapply(data, checkImpossibleChanges)
if (sum(unlist(imp.change)) > 0)
{
cat('For some groups, there are impossible changes;\n')
cat('This occurs for groups\n', which(sapply(imp.change,
function(ic){(sum(ic) > 0)})),'\n')
impossible <- TRUE
}
z.changes <- sapply(data, checkZeroChanges)
if (sum(z.changes >= 1))
{
cat(' For some groups and some period, there is no change.\n')
cat('This occurs for groups ', which(sapply(z.changes,
function(zc){(zc > 0)})),'\n')
zerochange <- TRUE
}
}
else
{
imp.change <- checkImpossibleChanges(data)
if (imp.change %in% c(1,3))
{
cat('There are some changes from structural zero to observed 1.\n')
}
if (imp.change %in% c(2,3))
{
cat('There are some changes from structural one to observed 0.\n')
}
if (imp.change >= 4)
{
cat('There are some impossible changes between structural and observed values,\n')
cat('some through intermediary NA.\n')
}
impossible <- (imp.change > 0)
zerochange <- FALSE
z.changes <- checkZeroChanges(data)
if (sum(z.changes >= 1))
{
cat(' For some period, there is no change.\n')
zerochange <- TRUE
}
}
if (impossible || zerochange)
{
cat('This is not allowed for likelihood-based estimation.\n')
if (impossible)
{
cat('Possible remedies are: represent periods by multiple groups,\n')
cat('or change the offending values.\n')
}
stop('Impossible changes')
}
}
## get data object into group format to save coping with two
## different formats
if (inherits(data, "sienaGroup"))
{
nGroup <- length(data)
}
else
{
nGroup <- 1
data <- sienaGroupCreate(list(data), singleOK=TRUE)
}
## add any effects needed for time dummies
tmp <- sienaTimeFix(effects, data)
data <- tmp$data
effects <- tmp$effects
if (!x$useStdInits)
{
if (!is.null(prevAns) && inherits(prevAns, "sienaFit"))
{
effects <- updateTheta(effects, prevAns)
}
}
## add any effects needed for settings model
# this now is replaced by adding the settings in getEffects,
# which is the more logical place.
# If all works, this can be deleted,
# and also the function addSettingsEffects can be deleted.
# I used this function as a template for the change to getEffects.
# I wonder why the next 8 lines cannot be dropped;
# gives error message "cannot find setting col".
# if (!is.null(x$settings))
# {
# effects <- addSettingsEffects(effects, x)
# }
# else
# {
# effects$setting <- rep("", nrow(effects))
# }
## find any effects not included which are needed for interactions
tmpEffects <- effects[effects$include, ]
interactionNos <- unique(c(tmpEffects$effect1, tmpEffects$effect2,
tmpEffects$effect3))
interactionNos <- interactionNos[interactionNos > 0]
interactions <- effects$effectNumber %in%
interactionNos
effects$requested <- effects$include
requestedEffects <- effects[effects$include, ]
effects$include[interactions] <- TRUE
effects <- effects[effects$include, ]
## split and rejoin both versions before continuing
depvarnames <- names(data[[1]]$depvars)
effects1 <- split(requestedEffects, requestedEffects$name)
effects1order <- match(c(depvarnames, "sde"), names(effects1))
requestedEffects <- do.call(rbind, effects1[effects1order])
row.names(requestedEffects) <- 1:nrow(requestedEffects)
effects1 <- split(effects, effects$name)
effects1order <- match(c(depvarnames, "sde"), names(effects1))
effects <- do.call(rbind, effects1[effects1order])
row.names(effects) <- 1:nrow(effects)
## now set up z provisionally
z$theta <- requestedEffects$initialValue
z$fixed <- requestedEffects$fix
z$test <- requestedEffects$test
z$pp <- length(z$test)
z$posj <- rep(FALSE, z$pp)
z$posj[requestedEffects$basicRate] <- TRUE
z$BasicRateFunction <- z$posj
## sort out names of user specified interaction effects
effects <- fixUpEffectNames(effects)
## copy interaction names to the requested effects
requestedEffects$effectName <- effects[effects$requested,
"effectName"]
requestedEffects$functionName <- effects[effects$requested,
"functionName"]
if (attr(data, "compositionChange"))
{
if (x$maxlike)
{
stop("Not possible to use maximum likelihood estimation ",
"with composition change")
}
}
## if not specified whether conditional or nor, set to conditional
## iff there is only one dependent variable (therefore number 1)
## and not maxlike
if (is.na(x$cconditional))
{
x$cconditional <- !x$maxlike && (length(depvarnames) == 1)
if (x$cconditional)
{
x$condvarno <- 1
}
}
if (any(hasSettings(data)))
{
if (x$cconditional)
{
stop("Estimation for a settings model requires ",
"conditional estimation")
}
}
types <- sapply(data[[1]]$depvars, function(x) attr(x, "type"))
## now check if conditional estimation is OK and copy to z if so
z$cconditional <- FALSE
if ((x$cconditional) & (!(attr(data, "compositionChange"))))
{
if (x$maxlike)
{
stop("Conditional estimation is not possible with",
"maximum likelihood method")
}
## if (nets == 1) not sure if this is necessary
## {
z$cconditional <- TRUE
## find the conditioning variable
observations <- attr(data, "observations")
## this is actual number of waves to process
if (x$condname != "")
{
z$condvarno <- match(x$condname, attr(data, "netnames"))
if (is.na(z$condvarno))
{
cat("\nNo match between variable name in algorithm object\n")
cat("and those in data set.\n")
cat("Algorithm object: ", x$condname, "\n")
cat("Data set: ", attr(data, "netnames"), "\n")
stop("Incorrect variable name.\n")
}
z$condname <- x$condname
}
else
{
z$condvarno <- x$condvarno
z$condname <- attr(data, "netnames")[x$condvarno]
}
z$condtype <- attr(data, "types")[z$condvarno]
if (z$condtype == "oneMode")
z$symmetric <- attr(data, "symmetric")[[z$condvarno]]
else
z$symmetric <- FALSE
## find the positions of basic rate effects for this network
z$condvar <-
(1:nrow(requestedEffects))[requestedEffects$name==
z$condname][1:observations]
z$theta<- z$theta[-z$condvar]
z$fixed<- z$fixed[-z$condvar]
z$test<- z$test[-z$condvar]
z$pp<- z$pp - length(z$condvar)
z$scale<- z$scale[-z$condvar]
z$BasicRateFunction <- z$posj[-z$condvar]
z$posj <- z$posj[-z$condvar]
z$theta[z$posj] <-
z$theta[z$posj] / requestedEffects$initialValue[z$condvar]
}
z$qq <- z$pp
## unpack data and put onto f anything we may need next time round.
f <- lapply(data, function(xx, x) unpackData(xx, x), x=x)
attr(f, "netnames") <- attr(data, "netnames")
attr(f, "symmetric") <- attr(data, "symmetric")
attr(f, "allUpOnly") <- attr(data, "allUpOnly")
attr(f, "allDownOnly") <- attr(data, "allDownOnly")
attr(f, "allHigher") <- attr(data, "allHigher")
attr(f, "allDisjoint") <- attr(data, "allDisjoint")
attr(f, "allAtLeastOne") <- attr(data, "allAtLeastOne")
attr(f, "anyUpOnly") <- attr(data, "anyUpOnly")
attr(f, "anyDownOnly") <- attr(data, "anyDownOnly")
attr(f, "anyHigher") <- attr(data, "anyHigher")
attr(f, "anyDisjoint") <- attr(data, "anyDisjoint")
attr(f, "anyAtLeastOne") <- attr(data, "anyAtLeastOne")
attr(f, "types") <- attr(data, "types")
attr(f, "observations") <- attr(data, "observations")
attr(f, "compositionChange") <- attr(data, "compositionChange")
attr(f, "exooptions") <- attr(data, "exooptions")
attr(f, "groupPeriods") <- attr(data, "groupPeriods")
attr(f, "periodNos") <- attr(data, "periodNos")
attr(f, "numberNonMissingNetwork") <-
attr(data, "numberNonMissingNetwork")
attr(f, "numberMissingNetwork") <- attr(data, "numberMissingNetwork")
attr(f, "numberNonMissingBehavior") <-
attr(data, "numberNonMissingBehavior")
attr(f, "numberMissingBehavior") <- attr(data, "numberMissingBehavior")
## attr(f, "totalMissings") <- attr(data, "totalMissings")
if (x$maxlike && x$FinDiff.method)
{
stop("Finite difference method for derivatives not available",
"with Maximum likelihood method")
}
## maxlike not available for symmetric networks; or is it?.
## check model type: default for symmetric is type 2 (forcing model).
## maxlike only for some model types? ABCD
syms <- attr(data,"symmetric")[ attr(data,"types") %in% c("oneMode","bipartite")]
z$FinDiffBecauseSymmetric <- FALSE
z$modelType <- x$modelType
z$behModelType <- x$behModelType
if (any(!is.na(syms) & syms))
{
## z$FinDiff.method <- TRUE
## z$FinDiffBecauseSymmetric <- TRUE
if (x$maxlike)
{
# stop("Maximum likelihood method not implemented",
# "for symmetric networks")
}
symms <- syms
symms[is.na(symms)] <- FALSE
z$modelType[(z$modelType == 1) & symms] <- 2
}
if (z$cconditional)
{
attr(f, "change") <-
sapply(f, function(xx)as.integer(attr(xx$depvars[[z$condname]],
"distance")))
attr(f,"condEffects") <- requestedEffects[z$condvar,]
effcondvar <-
(1:nrow(effects))[effects$name==
z$condname][1:observations]
effects <- effects[-effcondvar, ]
requestedEffects <- requestedEffects[-z$condvar,]
}
## use previous dfra only if everything matches including method
if (!is.null(prevAns) && inherits(prevAns, "sienaFit") &&
prevAns$maxlike == z$maxlike)
{
if (!is.null(prevAns$dfra) &&
nrow(prevAns$dfra) == nrow(requestedEffects) &&
all(rownames(prevAns$dfra) == requestedEffects$shortName) &&
all(prevAns$fix == requestedEffects$fix) &&
!is.null(prevAns$sf))
{
z$haveDfra <- TRUE
z$dfra <- prevAns$dfra
z$dinv <- prevAns$dinv
# z$dinvv must not be taken from prevAns,
# because the value of diagonalize
# is defined in x and may have changed.
# Therefore here we copy the corresponding lines
# from phase1.r.
# Partial diagonalization of derivative matrix
# for use if 0 < x$diagonalize < 1.
temp <- (1-x$diagonalize)*z$dfra +
x$diagonalize*diag(diag(z$dfra), nrow=dim(z$dfra)[1])
temp[z$fixed, ] <- 0.0
temp[, z$fixed] <- 0.0
diag(temp)[z$fixed] <- 1.0
# Invert this matrix
z$dinvv <- solve(temp)
z$sf <- prevAns$sf
# check for backward compatibility with pre-1.1-220 versions:
if (is.null(prevAns$regrCoef))
{
z$regrCoef <- rep(0, z$pp)
z$regrCor <- rep(0, z$pp)
}
else
{
z$regrCoef <- prevAns$regrCoef
z$regrCor <- prevAns$regrCor
}
}
}
z$effects <- effects
z$requestedEffects <- requestedEffects
}
else ## initC, i.e just send already set up data into new processes
{
f <- FRANstore()
## Would like f to be just the data objects plus the attributes
## but need the effects later. Also a few other things,
## which probably could be attributes but are not!
## They will be automatically removed: if needed they must be readded
ff <- f
nGroup <- f$nGroup
f[(nGroup + 1): length(f)] <- NULL
}
pData <- .Call(C_setupData, PACKAGE=pkgname,
lapply(f, function(x)(as.integer(x$observations))),
lapply(f, function(x)(x$nodeSets)))
ans <- .Call(C_OneMode, PACKAGE=pkgname,
pData, lapply(f, function(x)x$nets))
ans <- .Call(C_Bipartite, PACKAGE=pkgname,
pData, lapply(f, function(x)x$bipartites))
ans <- .Call(C_Behavior, PACKAGE=pkgname,
pData, lapply(f, function(x)x$behavs))
ans <- .Call(C_Continuous, PACKAGE=pkgname,
pData, lapply(f, function(x)x$contbehavs))
ans <-.Call(C_ConstantCovariates, PACKAGE=pkgname,
pData, lapply(f, function(x)x$cCovars))
ans <-.Call(C_ChangingCovariates, PACKAGE=pkgname,
pData, lapply(f, function(x)x$vCovars))
ans <-.Call(C_DyadicCovariates, PACKAGE=pkgname,
pData, lapply(f, function(x)x$dycCovars))
ans <-.Call(C_ChangingDyadicCovariates, PACKAGE=pkgname,
pData, lapply(f, function(x)x$dyvCovars))
ans <-.Call(C_ExogEvent, PACKAGE=pkgname,
pData, lapply(f, function(x)x$exog))
## split the names of the constraints
higher <- attr(f, "allHigher")
disjoint <- attr(f, "allDisjoint")
atLeastOne <- attr(f, "allAtLeastOne")
froms <- sapply(strsplit(names(higher), ","), function(x)x[1])
tos <- sapply(strsplit(names(higher), ","), function(x)x[2])
ans <- .Call(C_Constraints, PACKAGE=pkgname,
pData, froms[higher], tos[higher],
froms[disjoint], tos[disjoint],
froms[atLeastOne], tos[atLeastOne])
##store the address
f$pData <- pData
## register a finalizer
ans <- reg.finalizer(f$pData, clearData, onexit = FALSE)
if (!initC)
{
storage.mode(effects$parm) <- "integer"
storage.mode(effects$group) <- "integer"
storage.mode(effects$period) <- "integer"
effects$effectPtr <- rep(NA, nrow(effects))
if (any(attr(f,"types") == "continuous"))
{
splitFactor <- factor(effects$name, levels=c(attr(f, "netnames"), "sde"))
}
else
{
splitFactor <- factor(effects$name, levels=attr(f, "netnames"))
}
if (!all(attr(f,"netnames") %in% effects$name))
{
stop("Must have at least one effect for each dependent variable")
}
myeffects <- split(effects, splitFactor)
myCompleteEffects <- myeffects
## remove interaction effects and save till later
basicEffects <-
lapply(myeffects, function(x)
{
x[!x$shortName %in% c("unspInt", "behUnspInt"), ]
}
)
basicEffectsl <-
lapply(myeffects, function(x)
{
!x$shortName %in% c("unspInt", "behUnspInt")
}
)
interactionEffects <-
lapply(myeffects, function(x)
{
x[x$shortName %in% c("unspInt", "behUnspInt"), ]
}
)
interactionEffectsl <-
lapply(myeffects, function(x)
{
x$shortName %in% c("unspInt", "behUnspInt")
}
)
## store effects objects as we may need to recreate them
f$interactionEffects <- interactionEffects
f$basicEffects <- basicEffects
f$interactionEffectsl <- interactionEffectsl
f$basicEffectsl <- basicEffectsl
}
else
{
myCompleteEffects <- ff$myCompleteEffects
basicEffects <- ff$basicEffects
interactionEffects <- ff$interactionEffects
basicEffectsl <- ff$basicEffectsl
interactionEffectsl <- ff$interactionEffectsl
types <- ff$types
}
ans <- .Call(C_effects, PACKAGE=pkgname, pData, basicEffects)
pModel <- ans[[1]][[1]]
for (i in seq(along=(ans[[2]]))) ## ans[[2]] is a list of lists of
## pointers to effects. Each list, here i, corresponds to one
## dependent variable
## The first are NULL, for rate effects, which are not used for interactions.
{
effectPtr <- ans[[2]][[i]]
basicEffects[[i]]$effectPtr <- effectPtr
interactionEffects[[i]]$effect1 <-
basicEffects[[i]]$effectPtr[match(interactionEffects[[i]]$effect1,
basicEffects[[i]]$effectNumber)]
interactionEffects[[i]]$effect2 <-
basicEffects[[i]]$effectPtr[match(interactionEffects[[i]]$effect2,
basicEffects[[i]]$effectNumber)]
interactionEffects[[i]]$effect3 <-
basicEffects[[i]]$effectPtr[match(interactionEffects[[i]]$effect3,
basicEffects[[i]]$effectNumber)]
}
ans <- .Call(C_interactionEffects, PACKAGE=pkgname,
pModel, interactionEffects)
## copy these pointers to the interaction effects and then insert in
## effects object in the same rows for later use
for (i in 1:length(ans[[1]])) ## ans is a list of lists of
## pointers to effects. Each list corresponds to one
## dependent variable
{
if (nrow(interactionEffects[[i]]) > 0)
{
effectPtr <- ans[[1]][[i]]
interactionEffects[[i]]$effectPtr <- effectPtr
}
myCompleteEffects[[i]][basicEffectsl[[i]], ] <- basicEffects[[i]]
myCompleteEffects[[i]][interactionEffectsl[[i]],] <-
interactionEffects[[i]]
##myeffects[[i]] <- myeffects[[i]][order(myeffects[[i]]$effectNumber),]
}
## remove the effects only created as underlying effects
## for interaction effects. first store the original for use next time
myeffects <- lapply(myCompleteEffects, function(x)
{
x[x$requested, ]
}
)
##store address of model
f$pModel <- pModel
ans <- reg.finalizer(f$pModel, clearModel, onexit = FALSE)
if (x$MaxDegree == 0 || is.null(x$MaxDegree))
{
MAXDEGREE <- NULL
}
else
{
MAXDEGREE <- x$MaxDegree
storage.mode(MAXDEGREE) <- "integer"
}
if (x$UniversalOffset == 0 || is.null(x$UniversalOffset))
{
UNIVERSALOFFSET <- NULL
}
else
{
UNIVERSALOFFSET <- x$UniversalOffset
storage.mode(UNIVERSALOFFSET) <- "double"
}
if ((length(x$modelType) == 0)||all(x$modelType == 0) || is.null(x$modelType))
{
MODELTYPE <- NULL
}
else
{
MODELTYPE <- x$modelType
storage.mode(MODELTYPE) <- "integer"
}
if ((length(x$behModelType) == 0)||(x$behModelType == 0) || is.null(x$behModelType))
{
BEHMODELTYPE <- NULL
}
else
{
BEHMODELTYPE <- x$behModelType
storage.mode(BEHMODELTYPE) <- "integer"
}
if (z$cconditional)
{
CONDVAR <- z$condname
CONDTARGET <- attr(f, "change")
## cat(CONDTARGET, "\n")
}
else
{
CONDVAR <- NULL
CONDTARGET <- NULL
}
simpleRates <- TRUE
if (any(!z$effects$basicRate & z$effects$type =="rate"))
{
simpleRates <- FALSE
}
z$simpleRates <- simpleRates
ans <- .Call(C_setupModelOptions, PACKAGE=pkgname,
pData, pModel, MAXDEGREE, UNIVERSALOFFSET, CONDVAR, CONDTARGET,
profileData, z$parallelTesting, MODELTYPE, BEHMODELTYPE,
z$simpleRates, x$normSetRates)
if (!initC)
{
ans <- .Call(C_getTargets, PACKAGE=pkgname, pData, pModel, myeffects,
z$parallelTesting, returnActorStatistics=FALSE,
returnStaticChangeContributions=FALSE)
##stop("done")
## create a grid of periods with group names in case want to
## parallelize using this or to access chains easily
groupPeriods <- attr(f, "groupPeriods")
z$callGrid <- cbind(rep(1:nGroup, groupPeriods - 1),
as.vector(unlist(sapply(groupPeriods - 1,
function(x) 1:x))))
if (!x$maxlike)
{
z$targets <- rowSums(ans)
z$targets2 <- ans
# For the moment, the following is an undocumented and hidden option.
# This replaces the targets calculated from the data
# by user-defined targets.
if ((!is.null(x$targets[1])) & (nGroup == 1) & (groupPeriods[1] == 2))
{
if (length(z$targets) == length(x$targets))
{
z$targets <- x$targets
z$targets2 <- matrix(x$targets, length(x$targets), 1)
message('Note: targets taken from algorithm object.')
cat('\n')
print(z$targets)
cat('\n')
}
}
}
else
{
z$targets <- rep(0, z$pp)
z$targets2 <- ans
z$targets2[] <- 0
z$maxlikeTargets <- rowSums(ans)
z$maxlikeTargets2 <- ans
z$mult <- x$mult
length.nrunMH <- length(colSums(z$maxlikeTargets2[z$requestedEffects$basicRate,
, drop=FALSE ]))
if ((length(z$mult) >= 2) &&
(length(z$mult) != length.nrunMH))
{
stop(paste('Length of parameter mult in the algorithm object is incorrect',
'(should be 1 or', length.nrunMH, ').'))
}
z$nrunMH <- floor(z$mult *
colSums(z$maxlikeTargets2[z$requestedEffects$basicRate,
, drop=FALSE ]))
## make the number pretty
z$nrunMH <- ifelse (z$nrunMH > 100,
round(z$nrunMH / 100) * 100, z$nrunMH)
##thetaMat is to allow different thetas for each group in Bayes
z$thetaMat <- matrix(z$theta, nrow=nGroup, ncol=z$pp, byrow=TRUE)
}
}
# Here came an error
# Error: INTEGER() can only be applied to a 'integer', not a 'double'
# This was because storage.mode had not been set properly for some variable
if (x$maxlike)
{
if (!initC)
{
## set up chains and do initial steps
types <- attr(f, "types")
nbrNonMissNet <- attr(f, "numberNonMissingNetwork")
nbrMissNet <- attr(f, "numberMissingNetwork")
nbrNonMissBeh <- attr(f, "numberNonMissingBehavior")
nbrMissBeh <- attr(f, "numberMissingBehavior")
if (sum(nbrMissNet + nbrNonMissNet) > 0)
{
z$prmin <- nbrMissNet/ (nbrMissNet + nbrNonMissNet)
}
else
{
z$prmin <- rep(0, length(nbrMissNet))
}
if (sum(nbrMissBeh + nbrNonMissBeh) > 0)
{
z$prmib <- nbrMissBeh/ (nbrMissBeh + nbrNonMissBeh)
}
else
{
z$prmib <- rep(0, length(nbrMissBeh))
}
## localML
if (is.null(x$localML))
{
z$localML <- FALSE
} else {
z$localML <- x$localML
}
local <- ifelse(is.na(effects$local[effects$include]),
FALSE, effects$local[effects$include])
if (z$localML & any(!local))
{
stop("Non-local effect chosen.")
}
z$probs <- c(x$pridg, x$prcdg, x$prper, x$pripr, x$prdpr, x$prirms,
x$prdrms)
ans <- .Call(C_mlMakeChains, PACKAGE=pkgname, pData, pModel,
z$probs, z$prmin, z$prmib,
x$minimumPermutationLength,
x$maximumPermutationLength,
x$initialPermutationLength,
z$localML)
f$minimalChain <- ans[[1]]
f$chain <- ans[[2]]
}
else ## set up the initial chains in the sub processes
{
ans <- .Call(C_mlInitializeSubProcesses,
PACKAGE=pkgname, pData, pModel,
z$probs, z$prmin, z$prmib,
x$minimumPermutationLength,
x$maximumPermutationLength,
x$initialPermutationLength, ff$chain,
z$localML)
f$chain <- ff$chain
}
}
f$simpleRates <- simpleRates
f$myeffects <- myeffects
f$myCompleteEffects <- myCompleteEffects
if (!initC)
{
if (is.null(z$print) || z$print)
{
DataReport(z, x, f)
}
f$randomseed2 <- z$randomseed2
}
else
{
f$randomseed2 <- ff$randomseed2
}
f$observations <- attr(f, "observations") + 1
f$depNames <- names(f[[1]]$depvars)
f$groupNames <- names(f)[1:nGroup]
f$nGroup <- nGroup
f$types <- types
if (!initC)
{
z$f <- f
##z <- initForAlgorithms(z)
z$periodNos <- attr(data, "periodNos")
z$f$myeffects <- NULL
z$f$myCompleteEffects <- NULL
if (!returnDeps)
{
z$f[1:nGroup] <- NULL
}
}
if (initC || (z$int == 1 && z$int2 == 1 &&
(is.null(z$nbrNodes) || z$nbrNodes == 1)))
{
f[1:nGroup] <- NULL
}
FRANstore(f) ## store f in FRANstore
z$returnDeps <- returnDeps
z$returnDepsStored <- returnDeps
z$observations <- f$observations
z$returnChains <- returnChains
z$byGroup <- byGroup
z$byWave <- byWave
z$returnDataFrame <- returnDataFrame
z$nDependentVariables <- length(z$f$depNames)
z$x <- x # some elements -- named vectors -- may have changed
if (initC)
{
NULL
}
else
{
z
}
}
##@createEdgeLists siena07 Reformat data for C++
createEdgeLists<- function(mat, matorig, bipartite)
{
## mat1 is basic values, with missings and structurals replaced
tmp <- lapply(1 : nrow(mat), function(x, y)
{
mymat <- matrix(0, nrow = sum(y[x, ] > 0), ncol = 3)
mymat[, 1] <- x
mymat[, 2] <- which(y[x, ] != 0)
mymat[, 3] <- y[x, mymat[, 2]]
mymat
}, y = mat)
mat1 <- do.call(rbind, tmp)
## mat2 reverts to matorig to get the missing values
tmp <- lapply(1 : nrow(matorig), function(x, y)
{
mymat <- matrix(0, nrow = sum(is.na(y[x, ])), ncol = 3)
mymat[, 1] <- x
mymat[, 2] <- which(is.na(y[x, ]))
mymat[, 3] <- 1
mymat
}, y = matorig)
mat2 <- do.call(rbind, tmp)
## remove the diagonal if not bipartite
if (!bipartite)
{
mat2 <- mat2[mat2[, 1] != mat2[, 2], , drop=FALSE]
}
## mat3 structurals
struct <- mat1[,3] %in% c(10, 11)
mat1[struct, 3] <- mat1[struct,3] - 10
mat3 <- mat1[struct, , drop=FALSE]
mat3[, 3] <- 1
mat1 <- mat1[!mat1[,3] == 0, , drop=FALSE] ##remove any zeros just created
##fix up storage mode to be integer
storage.mode(mat1) <- "integer"
storage.mode(mat2) <- "integer"
storage.mode(mat3) <- "integer"
## add attribute of size
if (bipartite)
{
attr(mat1, "nActors") <- c(nrow(mat), ncol(mat))
attr(mat2, "nActors") <- c(nrow(mat), ncol(mat))
attr(mat3, "nActors") <- c(nrow(mat), ncol(mat))
}
else
{
attr(mat1, "nActors") <- nrow(mat)
attr(mat2, "nActors") <- nrow(mat)
attr(mat3, "nActors") <- nrow(mat)
}
list(mat1 = t(mat1), mat2 = t(mat2), mat3 = t(mat3))
}
##@createCovarEdgeLists siena07 Reformat data for C++
createCovarEdgeList<- function(mat, matorig)
{
tmp <- lapply(1 : nrow(mat), function(x, y)
{
mymat <- matrix(0, nrow = sum(y[x, ] != 0), ncol = 3)
mymat[, 1] <- x
mymat[, 2] <- which(y[x, ] != 0)
mymat[, 3] <- y[x, mymat[, 2]]
mymat
}, y = mat)
mat1 <- do.call(rbind, tmp)
##mat2 reverts to matorig to get the missing values
tmp <- lapply(1 : nrow(matorig), function(x, y)
{
mymat <- matrix(0, nrow = sum(is.na(y[x, ])), ncol = 3)
mymat[, 1] <- x
mymat[, 2] <- which(is.na(y[x, ]))
mymat[, 3] <- 1
mymat
}, y = matorig)
mat2 <- do.call(rbind, tmp)
## add attribute of size
attr(mat1, "nActors1") <- nrow(mat)
attr(mat1, "nActors2") <- ncol(mat)
list(mat1=t(mat1), mat2=t(mat2))
}
##@unpackOneMode siena07 Reformat data for C++
unpackOneMode <- function(depvar, observations, compositionChange)
{
edgeLists <- vector("list", observations)
networks <- vector("list", observations)
actorSet <- attr(depvar, "nodeSet")
compActorSets <- sapply(compositionChange, function(x)attr(x, "nodeSet"))
thisComp <- match(actorSet, compActorSets)
compChange <- !is.na(thisComp)
if (compChange)
{
action <- attr(compositionChange[[thisComp]], "action")
ccOption <- attr(compositionChange[[thisComp]], "ccOption")
}
else
{
ccOption <- 0
action <- matrix(0, nrow=attr(depvar, "netdims")[1], ncol=observations)
}
## sort out composition change
## convertToStructuralZeros()?
sparse <- attr(depvar, "sparse")
allowOnly <- attr(depvar, "allowOnly")
if (sparse)
{
## require(Matrix)
## have a list of sparse matrices in triplet format
## with missings and structurals embedded and 0 based indices!
netmiss <- vector("list", observations)
for (i in 1:observations)
{
## extract this matrix
networks[[i]] <- depvar[[i]]
nActors <- nrow(depvar[[i]])
## stop if any duplicates
netmat <- cbind(networks[[i]]@i+1, networks[[i]]@j+1,
networks[[i]]@x)
if (any(duplicated(netmat[, 1:2])))
{
stop("duplicate entries in sparse matrix")
}
## 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
if (i == 1)
{
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)
{
mat1 <- networks[[i]]
## drop the diagonal, if present
diag(mat1) <- 0
mat1 <- cbind(mat1@i + 1, mat1@j + 1, mat1@x)
##missing edgelist
mat2 <- netmiss[[i]]
mat2[, 3] <- 1
## rows of mat1 with structural values
struct <- mat1[, 3] %in% c(10, 11)
## reset real data
mat1[struct, 3] <- mat1[struct, 3] - 10
## copy reset data to structural edgelist
mat3 <- mat1[struct, , drop = FALSE]
mat3[, 3] <- 1
## now remove the zeros from reset data
mat1 <- mat1[!mat1[, 3] == 0, ]
## do comp change
if (compChange)
{
## revert to sparse matrices temporarily
mat1 <- spMatrix(nrow=nActors, ncol=nActors, i = mat1[, 1],
j=mat1[, 2], x=mat1[, 3])
mat2 <- spMatrix(nrow=nActors, ncol=nActors, i = mat2[, 1],
j=mat2[, 2], x=mat2[, 3])
mat3 <- spMatrix(nrow=nActors, ncol=nActors, i = mat3[, 1],
j=mat3[, 2], x=mat3[, 3])
ones <- which(action[, i] == 1)
twos <- which(action[, i] == 2)
threes <- which(action[, i] == 3)
for (j in ones) ## False data is not preceded by anything real
{
if (ccOption %in% c(1, 2))
{
## find missing values for this actor
use <- mat2[j, ] > 0
## remove from real data (i.e. zero)
mat1[j, use] <- 0
mat1[use, j] <- 0
## remove from missing data
mat2[j, use] <- 0
mat2[use, j] <- 0
## remove from raw data for distances later
depvar[[i]][j, use] <- 0 ## zero
depvar[[i]][use, j] <- 0
depvar[[i]][j, j] <- NA
}
else if (ccOption == 3)
{
## add the row and column to the missing data
mat2[j, ] <- 1
mat2[, j] <- 1
mat2[j, j] <- 0
## set to missing in raw data for distances later
depvar[[i]][j, ] <- NA
depvar[[i]][, j] <- NA
}
}
for (j in threes) ## False data is preceded and followed by real
{
if (ccOption %in% c(1, 2))
{
## find missing values for this actor
use <- mat2[j, ] > 0
## remove these from mat2, the missing data
mat2[j, use] <- 0
mat2[use, j] <- 0
## carry forward
if (i == 1)
{
## 0 any matches from mat1, the real data
mat1[j, use] <- 0
mat1[use, j] <- 0
}
else
{
mat1[j, use] <- networks[[i-1]][j, use]
mat1[use, j] <- networks[[i-1]][use, j]
}
depvar[[i]][j, use] <- 0 ## not missing
depvar[[i]][use, j] <- 0
depvar[[i]][j, j] <- NA
}
else if (ccOption == 3)
{
## add the row and column to the missing data
mat2[j, ] <- 1
mat2[, j] <- 1
mat2[j, j] <- 0
depvar[[i]][j, ] <- NA
depvar[[i]][, j] <- NA
}
}
for (j in twos) ## False data is not followed by anything real
{
if (ccOption == 1)
{
## find missing values for this actor
use <- mat2[j, ] > 0
## remove these from mat2, the missing data
mat2[j, use] <- 0
mat2[use, j] <- 0
depvar[[i]][j, use] <- 0 ## not missing
depvar[[i]][use, j] <- 0
depvar[[i]][j, j] <- NA
## carry forward
if (i == 1)
{
## 0 any matches from mat1, the real data
mat1[j, use] <- 0
mat1[use, j] <- 0
}
else
{
mat1[j, use] <- networks[[i-1]][j , use]
mat1[use, j] <- networks[[i-1]][use, j]
}
}
else if (ccOption %in% c(2, 3))
{
## add the row and column to the missing data
mat2[j, ] <- 1
mat2[, j] <- 1
mat2[j, j] <- 0
depvar[[i]][j, ] <- NA
depvar[[i]][, j] <- NA
}
}
## now revert to triplet matrices, after updating networks
networks[[i]] <- mat1
mat1 <- cbind(mat1@i + 1, mat1@j + 1, mat1@x)
mat2 <- cbind(mat2@i + 1, mat2@j + 1, mat2@x)
mat3 <- cbind(mat3@i + 1, mat3@j + 1, mat3@x)
if (any (mat1[, 3] == 0) || any (mat2[, 3] == 0) ||
any (mat3[, 3] == 0))
{
stop("zero values in sparse matrices")
}
if (any (duplicated(mat1[, -3])) ||
any (duplicated(mat2[, -3])) ||
any (duplicated(mat3[, -3])))
{
stop("duplicate values in sparse matrices")
}
if (any (mat1[, 1] == mat1[, 2]) ||
any (mat2[, 1] == mat2[, 2]) ||
any (mat3[, 1] == mat3[, 2]))
{
stop("loop values in sparse matrices")
}
}
##fix up storage mode to be integer
storage.mode(mat1) <- "integer"
storage.mode(mat2) <- "integer"
storage.mode(mat3) <- "integer"
## add attribute of size
attr(mat1,"nActors") <- nActors
attr(mat2,"nActors") <- nActors
attr(mat3,"nActors") <- nActors
if (i < observations)
{
## recreate the distance etc
mymat1 <- depvar[[i]]
mymat2 <- depvar[[i + 1]]
##remove structural values
x1 <- mymat1@x
x2 <- mymat2@x
x1[x1 %in% c(10, 11)] <- NA
x2[x2 %in% c(10, 11)] <- NA
mymat1@x <- x1
mymat2@x <- x2
diag(mymat1) <- 0
diag(mymat2) <- 0
mydiff <- mymat2 - mymat1
attr(depvar, "distance")[i] <- sum(mydiff != 0,
na.rm = TRUE)
if (allowOnly)
{
if (all(mydiff@x >= 0, na.rm=TRUE))
{
attr(depvar, "uponly")[i] <- TRUE
}
if (all(mydiff@x <= 0, na.rm=TRUE))
{
attr(depvar, "downonly")[i] <- TRUE
}
}
}
edgeLists[[i]] <- list(mat1 = t(mat1), mat2 = t(mat2),
mat3 = t(mat3))
}
}
else
{
for (i in 1:observations) ## carry missings forward if exist
{
networks[[i]] <- depvar[, , i]
if (i == 1)
networks[[i]][is.na(depvar[, , i])] <-0
else ##carry missing forward!
networks[[i]][is.na(depvar[, , i])] <-
networks[[i-1]][is.na(depvar[, , i])]
}
for (i in 1:observations)
{
ones <- which(action[, i] == 1)
twos <- which(action[, i] == 2)
threes <- which(action[, i] == 3)
for (j in ones) ## False data is not preceded by anything real
{
if (ccOption %in% c(1, 2))
{
use <- is.na(depvar[j, , i])
depvar[j, use, i] <- 0 ## not missing
depvar[use, j, i] <- 0
depvar[j, j, i] <- NA
networks[[i]][j, use] <- 0 ## zero
networks[[i]][use, j] <- 0
}
else if (ccOption == 3)
{
depvar[j, , i] <- NA ## missing
depvar[, j, i] <- NA
}
}
for (j in threes) ## False data is preceded and followed by real
{
if (ccOption %in% c(1, 2))
{
use <- is.na(depvar[j, , i])
depvar[j, use, i] <- 0 ## not missing
depvar[use, j, i] <- 0
depvar[j, j, i] <- NA
## carry forward already done
if (i == 1)
{
networks[[i]][j, use] <- 0
networks[[i]][use, j] <- 0
}
else
{
networks[[i]][j, use] <- networks[[i-1]][j, use]
networks[[i]][use, j] <- networks[[i-1]][use, j]
}
}
else if (ccOption == 3)
{
depvar[j, , i] <- NA ## missing
depvar[, j, i] <- NA
}
}
for (j in twos) ## False data is not followed by anything real
{
if (ccOption == 1)
{
use <- is.na(depvar[j, , i])
depvar[j, use, i] <- 0 ## not missing
depvar[use, j, i] <- 0
depvar[j, j, i] <- NA
## carry forward already done
if (i == 1)
{
networks[[i]][j, use] <- 0
networks[[i]][use, j] <- 0
}
else
{
networks[[i]][j, use] <- networks[[i-1]][j, use]
networks[[i]][use, j] <- networks[[i-1]][use, j]
}
}
else if (ccOption %in% c(2, 3))
{
depvar[j, , i] <- NA ## missing
depvar[, j, i] <- NA
}
}
}
for (i in 1:observations)
{
if (i < observations)
{
## recreate distances, as we have none in c++. (no longer true)
mymat1 <- depvar[, ,i, drop=FALSE]
mymat2 <- depvar[, ,i + 1, drop=FALSE]
##remove structural values
mymat1[mymat1 %in% c(10, 11)] <- NA
mymat2[mymat2 %in% c(10, 11)] <- NA
## and the diagonal
diag(mymat1[, ,1]) <- 0
diag(mymat2[, ,1]) <- 0
mydiff <- mymat2 - mymat1
attr(depvar, "distance")[i] <- sum(mydiff != 0,
na.rm = TRUE)
if (allowOnly)
{
if (all(mydiff >= 0, na.rm=TRUE))
{
attr(depvar, "uponly")[i] <- TRUE
}
if (all(mydiff <= 0, na.rm=TRUE))
{
attr(depvar, "downonly")[i] <- TRUE
}
}
}
diag(networks[[i]]) <- 0
edgeLists[[i]] <- createEdgeLists(networks[[i]], depvar[, , i],
FALSE)
}
}
## add attribute of nodeset
attr(edgeLists, "nodeSet") <- attr(depvar, "nodeSet")
## add attribute of name
attr(edgeLists, "name") <- attr(depvar, "name")
## add attribute of distance
attr(edgeLists, "distance") <- attr(depvar, "distance")
## attr uponly and downonly
attr(edgeLists, "uponly") <- attr(depvar, "uponly")
attr(edgeLists, "downonly") <- attr(depvar, "downonly")
## attr symmetric
attr(edgeLists, "symmetric") <- attr(depvar, "symmetric")
## attr balmean
attr(edgeLists, "balmean") <- attr(depvar, "balmean")
attr(edgeLists, "structmean") <- attr(depvar, "structmean")
attr(edgeLists, "averageInDegree") <- attr(depvar, "averageInDegree")
attr(edgeLists, "averageOutDegree") <- attr(depvar, "averageOutDegree")
attr(edgeLists, "settingsinfo") <- attr(depvar, "settingsinfo")
return(edgeLists = edgeLists)
}
##@unpackBipartite siena07 Reformat data for C++
unpackBipartite <- function(depvar, observations, compositionChange)
{
edgeLists <- vector("list", observations)
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")
allowOnly <- attr(depvar, "allowOnly")
if (sparse)
{
## require(Matrix)
## have a list of sparse matrices in triplet format
## with missings and structurals embedded and 0 based indices!
netmiss <- vector("list", observations)
for (i in 1:observations)
{
## extract this matrix
networks[[i]] <- depvar[[i]]
nActors <- nrow(depvar[[i]])
nReceivers <- ncol(depvar[[i]])
## stop if any duplicates
netmat <- cbind(networks[[i]]@i+1, networks[[i]]@j+1,
networks[[i]]@x)
if (any(duplicated(netmat[, 1:2])))
{
stop("duplicate entries in sparse matrix")
}
## extract missing entries
netmiss[[i]] <- netmat[is.na(netmat[, 3]), , drop = FALSE]
## carry forward missing values if any
if (i == 1) # set missings to zero
{
netmat <- netmat[!is.na(netmat[,3]), ]
networks[[i]] <- spMatrix(nActors, nReceivers, 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)
{
mat1 <- networks[[i]]
mat1 <- cbind(mat1@i + 1, mat1@j + 1, mat1@x)
##missing edgelist
mat2 <- netmiss[[i]]
mat2[, 3] <- 1
## rows of mat1 with structural values
struct <- mat1[, 3] %in% c(10, 11)
## reset real data
mat1[struct, 3] <- mat1[struct, 3] - 10
## copy reset data to structural edgelist
mat3 <- mat1[struct, , drop = FALSE]
## now remove the zeros from reset data
mat1 <- mat1[!mat1[, 3] == 0, ]
## do comp change
if (compChange)
{
## revert to sparse matrices temporarily
mat1 <- spMatrix(nrow=nActors, ncol=nReceivers, i = mat1[, 1],
j=mat1[, 2], x=mat1[, 3])
mat2 <- spMatrix(nrow=nActors, ncol=nReceivers, i = mat2[, 1],
j=mat2[, 2], x=mat2[, 3])
mat3 <- spMatrix(nrow=nActors, ncol=nReceivers, i = mat3[, 1],
j=mat3[, 2], x=mat3[, 3])
ones <- which(action[, i] == 1)
twos <- which(action[, i] == 2)
threes <- which(action[, i] == 3)
for (j in ones) ## False data is not preceded by anything real
{
if (ccOption %in% c(1, 2))
{
## find missing values for this actor
use <- mat2[j, ] > 0
## remove from real data (i.e. zero)
mat1[j, use] <- 0
## remove from missing data
mat2[j, use] <- 0
## remove from raw data for distances later
depvar[[i]][j, use] <- 0 ## zero
}
else if (ccOption == 3)
{
## add the row to the missing data
mat2[j, ] <- 1
## set to missing in raw data for distances later
depvar[[i]][j, ] <- NA
}
}
for (j in threes) ## False data is preceded and followed by real
{
if (ccOption %in% c(1, 2))
{
## find missing values for this actor
use <- mat2[j, ] > 0
## remove these from mat2, the missing data
mat2[j, use] <- 0
## carry forward
if (i == 1)
{
## 0 any matches from mat1, the real data
mat1[j, use] <- 0
}
else
{
mat1[j, use] <- networks[[i-1]][j, use]
}
depvar[[i]][j, use] <- 0 ## not missing
}
else if (ccOption == 3)
{
## add the row to the missing data
mat2[j, ] <- 1
depvar[[i]][j, ] <- NA
}
}
for (j in twos) ## False data is not followed by anything real
{
if (ccOption == 1)
{
## find missing values for this actor
use <- mat2[j, ] > 0
## remove these from mat2, the missing data
mat2[j, use] <- 0
depvar[[i]][j, use] <- 0 ## not missing
## carry forward
if (i == 1)
{
## 0 any matches from mat1, the real data
mat1[j, use] <- 0
}
else
{
mat1[j, use] <- networks[[i-1]][j , use]
}
}
else if (ccOption %in% c(2, 3))
{
## add the row to the missing data
mat2[j, ] <- 1
depvar[[i]][j, ] <- NA
}
}
## now revert to triplet matrices, after updating networks
networks[[i]] <- mat1
mat1 <- cbind(mat1@i + 1, mat1@j + 1, mat1@x)
mat2 <- cbind(mat2@i + 1, mat2@j + 1, mat2@x)
mat3 <- cbind(mat3@i + 1, mat3@j + 1, mat3@x)
if (any (mat1[, 3] == 0) || any (mat2[, 3] == 0) ||
any (mat3[, 3] == 0))
{
stop("zero values in sparse matrices")
}
if (any (duplicated(mat1[, -3])) ||
any (duplicated(mat2[, -3])) ||
any (duplicated(mat3[, -3])))
{
stop("duplicate values in sparse matrices")
}
}
##fix up storage mode to be integer
storage.mode(mat1) <- "integer"
storage.mode(mat2) <- "integer"
storage.mode(mat3) <- "integer"
## add attribute of size
attr(mat1,"nActors") <- c(nActors, nReceivers)
attr(mat2,"nActors") <- c(nActors, nReceivers)
attr(mat3,"nActors") <- c(nActors, nReceivers)
if (i < observations)
{
## recreate the distance etc
mymat1 <- depvar[[i]]
mymat2 <- depvar[[i + 1]]
##remove structural values
x1 <- mymat1@x
x2 <- mymat2@x
x1[x1 %in% c(10, 11)] <- NA
x2[x2 %in% c(10, 11)] <- NA
mymat1@x <- x1
mymat2@x <- x2
mydiff <- mymat2 - mymat1
attr(depvar, "distance")[i] <- sum(mydiff != 0,
na.rm = TRUE)
if (allowOnly)
{
if (all(mydiff@x >= 0, na.rm=TRUE))
{
attr(depvar, "uponly")[i] <- TRUE
}
if (all(mydiff@x <= 0, na.rm=TRUE))
{
attr(depvar, "downonly")[i] <- TRUE
}
}
}
edgeLists[[i]] <- list(mat1 = t(mat1), mat2 = t(mat2),
mat3 = t(mat3))
}
}
else
{
for (i in 1:observations) ## carry missings forward if exist
{
networks[[i]] <- depvar[, , i]
if (i == 1)
networks[[i]][is.na(depvar[, , i])] <-0
else ##carry missing forward!
networks[[i]][is.na(depvar[, , i])] <-
networks[[i-1]][is.na(depvar[, , i])]
}
for (i in 1:observations)
{
ones <- which(action[, i] == 1)
twos <- which(action[, i] == 2)
threes <- which(action[, i] == 3)
for (j in ones) ## False data is not preceded by anything real
{
if (ccOption %in% c(1, 2))
{
use <- is.na(depvar[j, , i])
depvar[j, use, i] <- 0 ## not missing
networks[[i]][j, use] <- 0 ## zero
}
else if (ccOption == 3)
{
depvar[j, , i] <- NA ## missing
}
}
for (j in threes) ## False data is preceded and followed by real
{
if (ccOption %in% c(1, 2))
{
use <- is.na(depvar[j, , i])
depvar[j, use, i] <- 0 ## not missing
## carry forward already done
if (i == 1)
{
networks[[i]][j, use] <- 0
}
else
{
networks[[i]][j, use] <- networks[[i-1]][j, use]
}
}
else if (ccOption == 3)
{
depvar[j, , i] <- NA ## missing
}
}
for (j in twos) ## False data is not followed by anything real
{
if (ccOption == 1)
{
use <- is.na(depvar[j, , i])
depvar[j, use, i] <- 0 ## not missing
## carry forward already done
if (i == 1)
{
networks[[i]][j, use] <- 0
}
else
{
networks[[i]][j, use] <- networks[[i-1]][j, use]
}
}
else if (ccOption %in% c(2, 3))
{
depvar[j, , i] <- NA ## missing
}
}
}
for (i in 1:observations)
{
if (i < observations)
{
## recreate distances, as we have none in c++. (no longer true)
mymat1 <- depvar[,,i, drop=FALSE]
mymat2 <- depvar[,,i + 1,drop=FALSE]
##remove structural values
mymat1[mymat1 %in% c(10,11)] <- NA
mymat2[mymat2 %in% c(10,11)] <- NA
mydiff <- mymat2 - mymat1
attr(depvar, "distance")[i] <- sum(mydiff != 0,
na.rm = TRUE)
if (allowOnly)
{
if (all(mydiff >= 0, na.rm=TRUE))
{
attr(depvar, "uponly")[i] <- TRUE
}
if (all(mydiff <= 0, na.rm=TRUE))
{
attr(depvar, "downonly")[i] <- TRUE
}
}
}
edgeLists[[i]] <- createEdgeLists(networks[[i]], depvar[, , i], TRUE)
}
}
## add attribute of nodeset
attr(edgeLists, "nodeSet") <- attr(depvar, "nodeSet")
## add attribute of name
attr(edgeLists, "name") <- attr(depvar, "name")
## add attribute of distance
attr(edgeLists, "distance") <- attr(depvar, "distance")
## attr uponly and downonly
attr(edgeLists, "uponly") <- attr(depvar, "uponly")
attr(edgeLists, "downonly") <- attr(depvar, "downonly")
## attr symmetric
attr(edgeLists, "symmetric") <- attr(depvar, "symmetric")
## attr balmean
attr(edgeLists, "balmean") <- attr(depvar, "balmean")
## attr structmean
attr(edgeLists, "structmean") <- attr(depvar, "structmean")
attr(edgeLists, "averageOutDegree") <- attr(depvar, "averageOutDegree")
return(edgeLists = edgeLists)
}
##@unpackBehavior siena07 Reformat data for C++
unpackBehavior<- function(depvar, observations)
{
beh <- depvar[, 1, ]
behmiss <- is.na(beh)
allna <- apply(beh, 1, function(x)all(is.na(x)))
modes <- attr(depvar, "modes")
## in case of missing data, use imputation values, provided by user
## else carry forward missings (otherwise use the mode / median)
## allNAs: use modes, in case of continuous behavior depvar$modes
## contains medians
if (!is.null(attr(depvar, "imputationValues")) )
{
beh[behmiss] <- attr(depvar, "imputationValues")[behmiss]
}
else
{
beh[allna, ] <- rep(modes, each=sum(allna))
for (i in 2:observations)
{
beh[is.na(beh[, i]), i] <- beh[is.na(beh[, i]), i - 1]
}
for (i in (observations-1):1)
{
beh[is.na(beh[, i]), i] <- beh[is.na(beh[, i]), i +1]
}
}
## need a better definition of structural for behavior variables
## struct <- beh %in% c(10, 11)
## beh[struct] <- beh[struct] - 10
## behstruct <- beh
## behstruct[!struct] <- 0
## add attribute of nodeset
attr(beh, "nodeSet") <- attr(depvar, "nodeSet")
## add attribute of name
attr(beh, "name") <- attr(depvar, "name")
## attr uponly and downonly
attr(beh, "uponly") <- attr(depvar, "uponly")
attr(beh, "downonly") <- attr(depvar, "downonly")
## attr symmetric
attr(beh, "symmetric") <- attr(depvar, "symmetric")
## attr distance
attr(beh, "distance") <- attr(depvar, "distance")
## attr simMean
attr(beh, "simMean") <- attr(depvar, "simMean")
## attr simMeans
attr(beh, "simMeans") <- attr(depvar, "simMeans")
if (attr(depvar, "type") == "behavior")
{
beh <- round(beh)
storage.mode(beh) <- "integer"
}
else # type == "continuous"
{
storage.mode(beh) <- "double"
}
list(beh=beh, behmiss=behmiss)
}
##@convertToStructuralZeros Miscellaneous To be implemented
convertToStructuralZeros <- function()
{
}
##@unpackCDyad siena07 Reformat data for C++
unpackCDyad<- function(dycCovar)
{
sparse <- attr(dycCovar, "sparse")
bipartite <- attr(dycCovar, "type") == "bipartite"
if (sparse)
{
## have a list containing 1 sparse matrix in triplet format
## with missings embedded
## with 0 based indices!
varmat <- cbind(dycCovar[[1]]@i+1, dycCovar[[1]]@j+1, dycCovar[[1]]@x)
if (any(duplicated(varmat[, 1:2])))
{
stop("duplicate entries in sparse matrix dyadic covariate")
}
##drop the diagonal, if present - not for bipartite
if (!bipartite)
{
varmat <- varmat[varmat[,1] != varmat[, 2],]
}
mat1 <- varmat
mat1[is.na(varmat[, 3]), 3] <- attr(dycCovar, "mean")
mat1 <- mat1[!mat1[, 3] == 0, ]
## add attribute of dim
attr(mat1, "nActors1") <- nrow(dycCovar[[1]])
attr(mat1, "nActors2") <- ncol(dycCovar[[1]])
mat2 <- varmat[is.na(varmat[, 3]), , drop=FALSE]
mat2[, 3] <- 1
## add attribute of dim
attr(mat2,"nActors1") <- nrow(dycCovar[[1]])
attr(mat2,"nActors2") <- ncol(dycCovar[[1]])
edgeLists <- list(t(mat1), t(mat2))
}
else
{
if (!bipartite)
{
diag(dycCovar) <- 0
}
dycCovar1 <- dycCovar
dycCovar1[is.na(dycCovar1)] <- attr(dycCovar, "mean")
edgeLists <- createCovarEdgeList(dycCovar1, dycCovar)
}
## add attribute of nodesets
attr(edgeLists, "nodeSet") <- attr(dycCovar, "nodeSet")
## add attribute of type
attr(edgeLists, "type") <- attr(dycCovar, "type")
## add attribute of name
attr(edgeLists, "name") <- attr(dycCovar, "name")
## add attribute of mean
attr(edgeLists, "mean") <- attr(dycCovar, "mean")
return(edgeLists = edgeLists)
}
##@unpackVDyad siena07 Reformat data for C++
unpackVDyad<- function(dyvCovar, observations)
{
edgeLists <- vector("list", observations)
sparse <- attr(dyvCovar, "sparse")
means <- attr(dyvCovar, "meanp")
bipartite <- attr(dyvCovar, "type") == "bipartite"
if (sparse)
{
## have a list of sparse matrices in triplet format
## with 0 based indices!
for (i in 1:(observations - 1))
{
thisvar <- dyvCovar[[i]]
varmat <- cbind(thisvar@i+1, thisvar@j+1, thisvar@x)
## drop the diagonal, if present no - bipartite?
if (!bipartite)
{
varmat <- varmat[varmat[,1] != varmat[, 2],]
}
mat1 <- varmat
mat1[is.na(varmat[, 3]), 3] <- means[i]
mat1 <- mat1[!mat1[, 3] == 0, ]
mat2 <- varmat[is.na(varmat[, 3]),, drop=FALSE ]
mat2[, 3] <- 1
## add attribute of size
attr(mat1, "nActors1") <- nrow(dyvCovar[[i]])
attr(mat1, "nActors2") <- ncol(dyvCovar[[i]])
attr(mat2, "nActors1") <- nrow(dyvCovar[[i]])
attr(mat2, "nActors2") <- ncol(dyvCovar[[i]])
edgeLists[[i]] <- list(t(mat1), t(mat2))
}
}
else
{
for (i in 1:(observations - 1))
{
if (!bipartite)
{
diag(dyvCovar[, , i]) <- 0
}
thisvar <- dyvCovar[, , i]
thisvar[is.na(thisvar) ] <- means[i]
edgeLists[[i]] <- createCovarEdgeList(thisvar, dyvCovar[, , i])
}
}
## add attribute of nodeset
attr(edgeLists, "nodeSet") <- attr(dyvCovar, "nodeSet")
## add attribute of type
attr(edgeLists, "type") <- attr(dyvCovar, "type")
## add attribute of name
attr(edgeLists, "name") <- attr(dyvCovar, "name")
## add attribute of mean
attr(edgeLists, "mean") <- attr(dyvCovar, "mean")
return(edgeLists = edgeLists)
}
##@unpackData siena07 Reformat data for C++
unpackData <- function(data, x)
{
f <- NULL
observations<- data$observations
types <- sapply(data$depvars, function(x) attr(x, "type"))
f$nDepvars <- length(data$depvars)
oneModes <- data$depvars[types == "oneMode"]
Behaviors <- data$depvars[types == "behavior"]
continuousBehaviors <- data$depvars[types == "continuous"]
bipartites <- data$depvars[types == "bipartite"]
## add the settings
# oneModes <- lapply(oneModes, function(depvar) {
# name <- attr(depvar, "name")
# if (name %in% names(x$settings)) {
# # attr(depvar, "settings") <- c("universal", "primary", x$settings[[name]])
# attr(depvar, "settings") <- c(x$settings[[name]])
# }
# depvar
# })
f$nets <- lapply(oneModes, function(x, n, comp)
unpackOneMode(x, n, comp),
n = observations, comp=data$compositionChange)
names(f$nets) <- names(oneModes)
f$bipartites <- lapply(bipartites, function(x, n, comp)
unpackBipartite(x, n, comp),
n = observations, comp=data$compositionChange)
names(f$bipartites) <- names(bipartites)
f$behavs <- lapply(Behaviors, function(x, n) unpackBehavior(x, n),
n = observations)
names(f$behavs) <- names(Behaviors)
f$contbehavs <- lapply(continuousBehaviors, function(x, n)
unpackBehavior(x, n), n = observations)
names(f$contbehavs) <- names(continuousBehaviors)
f$observations <- observations
f$seed<- vector("list", observations - 1)
f$depvars <- data$depvars
f$nodeSets <- data$nodeSets
f$oneModes <- oneModes
f$Behaviors <- Behaviors
f$continuousBehaviors <- continuousBehaviors
f$oneModeUpOnly <- sapply(oneModes, function(x) attr(x, "uponly"))
f$oneModeDownOnly <- sapply(oneModes, function(x) attr(x, "downonly"))
f$behaviorUpOnly <- sapply(Behaviors, function(x) attr(x, "uponly"))
f$behaviorDownOnly <- sapply(Behaviors, function(x) attr(x,
"downonly"))
f$distances <- sapply(data$depvars, function(x) attr(x, "distance"))
f$cCovars <- data$cCovars
f$vCovars <- data$vCovars
## dyadic covars need to be edgelists
f$dycCovars <- lapply(data$dycCovars, function(x) unpackCDyad(x))
f$dyvCovars <- lapply(data$dyvCovars, function(x,n) unpackVDyad(x,n),
n=observations)
## create the composition change event lists
f$exog <- lapply(data$compositionChange, function(x)
unpackCompositionChange(x))
f
}
##@unpackCompositionChange siena07 Reformat data for C++
unpackCompositionChange <- function(compositionChange)
{
atts <- attributes(compositionChange)
events <- atts$events
activeStart <- atts$activeStart
observations <- ncol(activeStart)
## check that there is someone there always
for (i in 1:(observations - 1))
{
activeAll <- sum(activeStart[, i] & activeStart[, i + 1])
if (activeAll < 2)
{
active <- sum(activeStart[, i])
if (active == 1)
stop("Only one active actor at start of period", i)
else if (active == 0)
stop("No active actors at start of period", i)
perEvents <- events[events$period == i,]
perEvents <- perEvents[order(perEvents$time),]
changes <- c(1, -1)[as.numeric(as.character(perEvents$event))]
active <- active + cumsum(changes)
if (any(active < 2))
{
stop("No/only one active actor(s) left.")
}
}
}
events <- events[events$time > 1e-10,]
exog <- list(events=events, activeStart=activeStart)
attr(exog, "nodeSet") <- attr(compositionChange, "nodeSet")
exog
}
##@fixUpEffectNames siena07 Replace # and construct interaction names
fixUpEffectNames <- function(effects)
{
## replace # by the parm value in function and effect names
effects$effectName <-
sapply(1:nrow(effects), function(x, y)
{
y <- y[x, ]
gsub("#", y$parm, y$effectName)
}, effects)
effects$functionName <-
sapply(1:nrow(effects), function(x, y)
{
y <- y[x, ]
gsub("#", y$parm, y$functionName)
}, y=effects)
##validate user-specified network interactions
interactions <- effects[effects$shortName == "unspInt" & effects$include &
effects$effect1 > 0, ]
if (nrow(interactions) > 0)
{
unspIntNames <-
sapply(1:nrow(interactions), function(x, y, z)
{
y <- y[x, ] ## get the interaction effect
twoway <- y$effect3 == 0
## now get the rows which are to interact
inter1 <- z[z$effectNumber == y$effect1, ]
if (nrow(inter1) != 1 )
{
stop("invalid network interaction specification: ",
"effect number 1")
}
inter2 <- z[z$effectNumber == y$effect2, ]
if (nrow(inter2) != 1 )
{
stop("invalid network interaction specification: ",
"effect number 2")
}
if (!twoway)
{
inter3 <- z[z$effectNumber == y$effect3, ]
if (nrow(inter3) != 1)
{
stop("invalid network interaction specification: ",
"effect number 3")
}
}
else
{
inter3 <- z[is.na(z$effectNumber), ]
## should be empty row
}
if (twoway)
{
if (inter1$name != inter2$name)
{
stop("invalid network interaction specification: ",
"must all be same network")
}
if (inter1$type != inter2$type)
{
stop("invalid network interaction specification: ",
"must all be same type: ",
"evaluation, endowment or creation")
}
}
else
{
if (inter1$name != inter2$name ||
inter1$name != inter3$name)
{
stop("invalid network interaction specification: ",
"must all be same network")
}
if (inter1$type != inter2$type ||
inter1$type != inter3$type)
{
stop("invalid network interaction specification: ",
"must all be ",
"same type: evaluation, endowment or creation ")
}
}
## check types
inters <- rbind(inter1, inter2, inter3)
egos <- which(inters$interactionType == "ego")
egoCount <- length(egos)
dyads <- which(inters$interactionType == "dyadic")
dyadCount <- length(dyads)
if (twoway)
{
if (egoCount < 1 && dyadCount != 2)
{
stop("invalid network interaction specification: ",
"must be at least one ego or both dyadic ",
"effects")
}
}
else
{
if (egoCount < 2 && (egoCount + dyadCount < 3))
{
stop("invalid network 3-way interaction specification: ",
"must be at least two ego effects ",
"or all ego or dyadic effects")
}
}
## construct a name
## make sure the egos are at the front of inters
if (egoCount > 0)
{
inters <- rbind(inters[egos, ], inters[-egos, ])
}
tmpnames <- inters$effectName
tmpnames[-1] <- sub(paste(inters$name[1], ": ",
sep=""), "", tmpnames[-1])
tmpname <- paste(tmpnames, collapse = " x ")
if (twoway && nchar(tmpname) < 38)
{
tmpname <- paste("int. ", tmpname)
}
if (!twoway)
{
tmpname <- paste("i3.", tmpname)
}
tmpname
}, y=interactions, z=effects)
effects[effects$shortName == "unspInt" & effects$include &
!is.na(effects$effect1), c("effectName", "functionName")] <-
unspIntNames
}
##validate user-specified behavior interactions
interactions <- effects[effects$shortName == "behUnspInt" &
effects$include &
effects$effect1 > 0, ]
if (nrow(interactions) > 0)
{
unspIntNames <-
sapply(1:nrow(interactions), function(x, y, z)
{
y <- y[x, ] ## get the interaction effect
twoway <- y$effect3 == 0
## now get the rows which are to interact
inter1 <- z[z$effectNumber == y$effect1, ]
if (nrow(inter1) != 1 )
{
stop("invalid behavior interaction specification: ",
"effect number 1")
}
inter2 <- z[z$effectNumber == y$effect2, ]
if (nrow(inter2) != 1 )
{
stop("invalid behavior interaction specification: ",
"effect number 2")
}
if (!twoway)
{
inter3 <- z[z$effectNumber == y$effect3, ]
if (nrow(inter3) != 1)
{
stop("invalid behavior interaction specification: ",
"effect number 3")
}
}
else
{
inter3 <- z[is.na(z$effectNumber), ]
## should be empty row
}
if (twoway)
{
if (inter1$name != inter2$name)
{
stop("invalid behavior interaction specification: ",
"must all be same behavior variable")
}
if (inter1$type != inter2$type)
{
stop("invalid behavior interaction specification: ",
"must be same type: evaluation, endowment ",
"or creation")
}
}
else
{
if (inter1$name != inter2$name ||
inter1$name != inter3$name)
{
stop("invalid behavior interaction specification: ",
"must all be same behavior variable")
}
if (inter1$type != inter2$type ||
inter1$type != inter3$type)
{
stop("invalid behavior interaction specification: ",
"must all be ",
"same type: evaluation, endowment or creation")
}
}
## check types - at most one should be not OK here
inters <- rbind(inter1, inter2, inter3)
if (length(which(inters$interactionType != "OK")) > 1)
{
stop("invalid behavior interaction specification: ",
"at most one effect with interactionType ",
"not OK is allowed")
}
##if (any(inters$interactionType != "OK"))
##{
## stop("invalid behavior interaction specification: ",
## "only effects with interactionType OK are allowed")
##}
## construct a name
tmpnames <- inters$effectName
tmpnames[-1] <- sub(paste("behavior ", inters$name[1], " ",
sep=""), "", tmpnames[-1])
tmpname <- paste(tmpnames, collapse = " x ")
if (twoway && nchar(tmpname) < 38)
{
tmpname <- paste("int. ", tmpname)
}
if (!twoway)
{
tmpname <- paste("i3.", tmpname)
}
tmpname
}, y=interactions, z=effects)
effects[effects$shortName == "behUnspInt" & effects$include &
!is.na(effects$effect1), c("effectName", "functionName")] <-
unspIntNames
}
##validate user-specified continuous interactions
interactions <- effects[effects$shortName == "contUnspInt" &
effects$include &
effects$effect1 > 0, ]
if (nrow(interactions) > 0)
{
unspIntNames <-
sapply(1:nrow(interactions), function(x, y, z)
{
#browser()
y <- y[x, ] ## get the interaction effect
twoway <- y$effect3 == 0
## now get the rows which are to interact
inter1 <- z[z$effectNumber == y$effect1, ]
if (nrow(inter1) != 1 )
{
stop("invalid behavior interaction specification: ",
"effect number 1")
}
inter2 <- z[z$effectNumber == y$effect2, ]
if (nrow(inter2) != 1 )
{
stop("invalid behavior interaction specification: ",
"effect number 2")
}
if (!twoway)
{
inter3 <- z[z$effectNumber == y$effect3, ]
if (nrow(inter3) != 1)
{
stop("invalid behavior interaction specification: ",
"effect number 3")
}
}
else
{
inter3 <- z[is.na(z$effectNumber), ]
## should be empty row
}
if (twoway)
{
if (inter1$name != inter2$name)
{
stop("invalid behavior interaction specification: ",
"must all be same behavior variable")
}
if (inter1$type != inter2$type)
{
stop("invalid behavior interaction specification: ",
"must be same type: evaluation")
}
}
else
{
if (inter1$name != inter2$name ||
inter1$name != inter3$name)
{
stop("invalid behavior interaction specification: ",
"must all be same behavior variable")
}
if (inter1$type != inter2$type ||
inter1$type != inter3$type)
{
stop("invalid behavior interaction specification: ",
"must all be same type: evaluation")
}
}
## check types - at most one should be not OK here
inters <- rbind(inter1, inter2, inter3)
if (length(which(inters$interactionType != "OK")) > 1)
{
stop("invalid behavior interaction specification: ",
"at most one effect with interactionType ",
"not OK is allowed")
}
##if (any(inters$interactionType != "OK"))
##{
## stop("invalid behavior interaction specification: ",
## "only effects with interactionType OK are allowed")
##}
## construct a name
tmpnames <- inters$effectName
tmpnames[-1] <- sub(paste("behavior ", inters$name[1], " ",
sep=""), "", tmpnames[-1])
tmpname <- paste(tmpnames, collapse = " x ")
if (twoway && nchar(tmpname) < 38)
{
tmpname <- paste("int. ", tmpname)
}
if (!twoway)
{
tmpname <- paste("i3.", tmpname)
}
tmpname
}, y=interactions, z=effects)
effects[effects$shortName == "contUnspInt" & effects$include &
!is.na(effects$effect1), c("effectName", "functionName")] <-
unspIntNames
}
effects
}
##@updateTheta siena07 Copy theta values from previous fit
updateTheta <- function(effects, prevAns, varName=NULL)
{
if (!inherits(effects, "data.frame"))
{
stop("effects is not a data.frame")
}
if (!inherits(prevAns, "sienaFit"))
{
stop("prevAns is not an RSiena fit object")
}
prevEffects <- prevAns$requestedEffects[which(prevAns$requestedEffects$type != 'gmm'),]
prevEffects$initialValue <- prevAns$theta
if (prevAns$cconditional)
{
condEffects <- attr(prevAns$f, "condEffects")
condEffects$initialValue <- prevAns$rate
prevEffects <- rbind(prevEffects, condEffects)
}
if (!is.null(varName))
{
prevEffects$name[!( prevEffects$name %in% varName)] <- ' '
}
oldlist <- apply(prevEffects, 1, function(x)
paste(x[c("name", "shortName",
"type", "groupName",
"interaction1", "interaction2",
"period", "effect1", "effect2", "effect3")],
collapse="|"))
efflist <- apply(effects, 1, function(x)
paste(x[c("name", "shortName",
"type", "groupName",
"interaction1", "interaction2",
"period", "effect1", "effect2", "effect3")],
collapse="|"))
use <- efflist %in% oldlist
effects$initialValue[use] <-
prevEffects$initialValue[match(efflist, oldlist)][use]
effects
}
##@addSettingseffects siena07 add extra rate effects for settings model
# addSettingsEffects <- function(effects, x)
# {
# ## need a list of settings (constant dyadic covariate name) for some or all
# ## dependent networks. If symmetric this is equivalent to a forcing model.
# ## The covariate actor sets should match the network actor sets.
# ## ?? would it ever make sense for bipartites? Allow it here for now and see!
# ## maybe better to add the settings pars to the data object but for now
# ## they are on the model with maxdegree. TODO write some validation
# varNames <- names(x$settings)
# nbrSettings <- sapply(x$settings, length)
# tmp <-
# lapply(1:length(x$settings), function(i)
# {
# ## find effects for this variable
# varEffects <- effects[effects$name == varNames[i], ]
# ## find relevant rate effects
# dupl <- varEffects[varEffects$basicRate, ]
# ## make extra copies
# newEffects <- dupl[rep(1:nrow(dupl),
# each = nbrSettings[i] + 2), ]
# newEffects <- split(newEffects, list(newEffects$group,
# newEffects$period))
# newEffects <- lapply(newEffects, function(dd)
# {
# dd$setting <- c("universal", "primary",
# x$settings[[i]])
# i1 <- regexpr("rate", dd$effectName)
# dd$effectName <-
# paste(substr(dd$effectName, 1, i1 - 2),
# dd$setting,
# substring(dd$effectName, i1))
# dd
# }
# )
# newEffects <- do.call(rbind, newEffects)
# ## add extra column to non basic rate effects
# varEffects$setting <- rep("", nrow(varEffects))
# ## combine them
# rbind(newEffects, varEffects[!varEffects$basicRate, ])
# })
# ## join them together again
# do.call(rbind, tmp)
# }
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.