### LOAD REQUIRED PACKAGES AND AUXILIARY FUNCTIONS
library(RSiena)
source('00_functions.R')
# estimate models for each village separately
for (i in 1:length(village.codes)) {
# define data objects
saom.talk <- eval(parse(text=paste0('v', i, '.adjmatrix')))
saom.covars <- eval(parse(text=paste0('v', i, '.attributes')))
# define two-mode network: use of preventive measures
saom.use <- as.matrix(saom.covars[,vars.dep])
# identify and healer
who.asha <- which(saom.covars[,'ASHA']==1)
who.healer <- which(saom.covars[,'Traditional Healer']==1)
# identify measures used by asha and healer
saom.use.asha <- saom.use[who.asha,]
if (length(who.healer)>0) {
saom.use.healer <- saom.use[who.healer,]
} else {
saom.use.healer <- rep(0, 8)
}
# remove asha and healer from networks and individual vars objects
saom.covars <- saom.covars[-c(who.asha, who.healer),]
saom.talk <- saom.talk[-c(who.asha, who.healer),-c(who.asha, who.healer)]
saom.use <- saom.use[-c(who.asha, who.healer),]
# define household network
saom.hh <- outer(saom.covars[,'Household'], saom.covars[,'Household'], '==')*1
# recode household to numeric
saom.covars[,'Household'] <- as.numeric(as.factor(saom.covars[,'Household']))
# modified talk network for "first wave"
saom.talk.mod <- saom.talk
y <- sample(which(saom.talk.mod>0), 1)
saom.talk.mod[y] <- saom.talk.mod[y] - 1
# modified use network for "first wave"
saom.use.mod <- saom.use
y <- sample(which(saom.use.mod>0), 1)
saom.use.mod[y] <- saom.use.mod[y] - 1
rm(y)
# set up Siena objects
s.nodes1 <- sienaNodeSet(nrow(saom.talk), nodeSetName='Actors')
s.nodes2 <- sienaNodeSet(ncol(saom.use), nodeSetName='Measures')
s.talk <- sienaDependent(array(c(saom.talk.mod, saom.talk),
dim=c(nrow(saom.talk), ncol(saom.talk), 2)),
allowOnly=FALSE, nodeSet=c('Actors'))
s.use <- sienaDependent(array(c(saom.use.mod, saom.use),
dim=c(nrow(saom.use), ncol(saom.use), 2)),
allowOnly=FALSE, type='bipartite', nodeSet=c('Actors', 'Measures'))
s.hh.dyad <- coDyadCovar(saom.hh, centered=F)
s.hh.indiv <- coCovar(saom.covars[,'Household'], centered=F)
s.female <- coCovar(saom.covars[,'Female'], centered=F)
s.headhouse <- coCovar(saom.covars[,'Head of household'], centered=F)
s.caresick <- coCovar(saom.covars[,'Carer for a sick person'], centered=F)
s.workfields <- coCovar(saom.covars[,'Works in fields'], centered=F)
s.age <- coCovar(saom.covars[,'Age'], centered=F)
s.educ <- coCovar(saom.covars[,'Education'], centered=F)
s.nonresptalk <- coCovar(saom.covars[,'Network size among non-respondents'], centered=F)
s.outvilltalk <- coCovar(saom.covars[,'Network size out of village'], centered=F)
s.ashatalk <- coCovar(saom.covars[,'Talks to ASHA'], centered=F)
s.ashause <- coCovar(saom.use.asha, centered=F, nodeSet='Measures')
s.healertalk <- coCovar(saom.covars[,'Talks to Healer'], centered=F)
s.healeruse <- coCovar(saom.use.healer, centered=F, nodeSet='Measures')
s.data <- sienaDataCreate(talk=s.talk, use=s.use,
household_dyadic=s.hh.dyad, household_indiv=s.hh.indiv,
female=s.female, headhouse=s.headhouse, caresick=s.caresick,
workfields=s.workfields, age=s.age, educ=s.educ,
nonresptalk=s.nonresptalk, outvilltalk=s.outvilltalk,
ashatalk=s.ashatalk, ashause=s.ashause,
healertalk=s.healertalk, healeruse=s.healeruse,
nodeSets=list(s.nodes1, s.nodes2))
# specify model
s.effects <- getEffects(s.data)
# fix rates
s.effects <- setEffect(s.effects, Rate, name='use', type='rate', initialValue=rates, fix=TRUE)
s.effects <- setEffect(s.effects, Rate, name='talk', type='rate', initialValue=rates, fix=TRUE)
# dep: talk
s.effects <- includeEffects(s.effects, recip, transTrip, transRecTrip, outAct, inPop, name='talk')
s.effects <- includeEffects(s.effects, sameX, name='talk', interaction1='household_indiv')
s.effects <- includeEffects(s.effects, egoX, altX, sameX, name='talk', interaction1='ashatalk')
s.effects <- includeEffects(s.effects, egoX, altX, sameX, name='talk', interaction1='healertalk')
if (length(who.healer)==0) {
s.effects <- setEffect(s.effects, egoX, name='talk', interaction1='healertalk', initialValue=0, fix=TRUE)
s.effects <- setEffect(s.effects, altX, name='talk', interaction1='healertalk', initialValue=0, fix=TRUE)
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='healertalk', initialValue=0, fix=TRUE)
}
s.effects <- includeEffects(s.effects, egoX, altX, sameX, name='talk', interaction1='female')
s.effects <- includeEffects(s.effects, egoX, altX, sameX, name='talk', interaction1='headhouse')
s.effects <- includeEffects(s.effects, egoX, altX, sameX, name='talk', interaction1='caresick')
s.effects <- includeEffects(s.effects, egoX, altX, sameX, name='talk', interaction1='workfields')
s.effects <- includeEffects(s.effects, egoX, altX, simX, name='talk', interaction1='age')
s.effects <- includeEffects(s.effects, egoX, altX, simX, name='talk', interaction1='educ')
s.effects <- includeEffects(s.effects, egoX, altX, simX, name='talk', interaction1='nonresptalk')
s.effects <- includeEffects(s.effects, egoX, altX, simX, name='talk', interaction1='outvilltalk')
s.effects <- includeEffects(s.effects, from, outActIntn, outPopIntn, name='talk', interaction1='use')
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use', par=1)
s.effects <- setEffect(s.effects, outPopIntn, name='talk', interaction1='use', par=1)
# interaction of hh and talk
s.effects <- includeEffects(s.effects, covNetNet, name='talk', interaction1='use', interaction2='household_indiv')
# dep: use
s.effects <- includeEffects(s.effects, cycle4, outAct, inPop, name='use')
s.effects <- includeEffects(s.effects, sameXCycle4, name='use', interaction1='household_indiv')
s.effects <- includeEffects(s.effects, egoX, name='use', interaction1='ashatalk')
s.effects <- includeEffects(s.effects, altX, name='use', interaction1='ashause')
s.effects <- includeInteraction(s.effects, egoX, altX, name='use', interaction1=c('ashatalk','ashause'))
s.effects <- includeEffects(s.effects, egoX, name='use', interaction1='healertalk')
s.effects <- includeEffects(s.effects, altX, name='use', interaction1='healeruse')
s.effects <- includeInteraction(s.effects, egoX, altX, name='use', interaction1=c('healertalk','healeruse'))
if (length(who.healer)==0) {
s.effects <- setEffect(s.effects, egoX, name='use', interaction1='healertalk', initialValue=0, fix=TRUE)
s.effects <- setEffect(s.effects, altX, name='use', interaction1='healeruse', initialValue=0, fix=TRUE)
s.effects <- includeInteraction(s.effects, egoX, altX, name='use', interaction1=c('healertalk','healeruse'), fix=TRUE)
}
s.effects <- includeEffects(s.effects, egoX, name='use', interaction1='female')
s.effects <- includeEffects(s.effects, egoX, name='use', interaction1='headhouse')
s.effects <- includeEffects(s.effects, egoX, name='use', interaction1='caresick')
s.effects <- includeEffects(s.effects, egoX, name='use', interaction1='workfields')
s.effects <- includeEffects(s.effects, egoX, name='use', interaction1='age')
s.effects <- includeEffects(s.effects, egoX, name='use', interaction1='educ')
s.effects <- includeEffects(s.effects, egoX, name='use', interaction1='nonresptalk')
s.effects <- includeEffects(s.effects, egoX, name='use', interaction1='outvilltalk')
s.effects <- includeEffects(s.effects, to, outActIntn, inActIntn, name='use', interaction1='talk')
s.effects <- setEffect(s.effects, outActIntn, name='use', interaction1='talk', par=1)
s.effects <- setEffect(s.effects, inActIntn, name='use', interaction1='talk', par=1)
# interaction of hh and talk
s.effects <- includeEffects(s.effects, sameWXClosure, name='use', interaction1='talk', interaction2='household_indiv')
# village-specific effect fixing
# village 1
if (i==1) {
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='healertalk',
initialValue=0, fix=TRUE, test=T)
}
# village 4
if (i==4) {
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='healertalk',
initialValue=-0.1, fix=T, test=T)
}
# village 5
if (i==5) {
s.effects <- setEffect(s.effects, inActIntn, name='use', interaction1='talk',
initialValue=0.1, fix=T, test=T)
}
# village 7
if (i==7) {
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, outPopIntn, name='talk', interaction1='use',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='caresick',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, egoX, name='talk', interaction1='caresick',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, altX, name='talk', interaction1='caresick',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, egoX, name='use', interaction1='caresick',
initialValue=1, fix=T, test=T)
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='workfields',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, egoX, name='talk', interaction1='workfields',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, altX, name='talk', interaction1='workfields',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, egoX, name='use', interaction1='ashatalk',
initialValue=0.5, fix=T, test=T)
}
# village 8
if (i==8) {
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=2, fix=T, test=T)
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='caresick',
initialValue=0, fix=T, test=T)
s.effects <- includeInteraction(s.effects, egoX, altX, name='use',
interaction1=c('healertalk','healeruse'),
fix=T, test=T)
s.effects <- setEffect(s.effects, altX, name='use', interaction1='healeruse',
initialValue=1, fix=T, test=T)
s.effects <- setEffect(s.effects, outActIntn, name='use', interaction1='talk',
initialValue=0, fix=T, test=T)
}
# village 10
if (i==10) {
s.effects <- setEffect(s.effects, egoX, name='talk', interaction1='caresick',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, altX, name='talk', interaction1='caresick',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='caresick',
initialValue=0, fix=T, test=T)
}
# fix mixed triad closure for GoF check
if (fix.netexp) {
s.effects <- setEffect(s.effects, from, name='talk', interaction1='use',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, to, name='use', interaction1='talk',
initialValue=0, fix=T)
}
# fix other crossnet effects for GoF check
if (fix.netsize) {
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, outPopIntn, name='talk', interaction1='use',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, outActIntn, name='use', interaction1='talk',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, inActIntn, name='use', interaction1='talk',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, egoX, name='use', interaction1='outvilltalk',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, egoX, name='use', interaction1='nonresptalk',
initialValue=0, fix=T)
}
# fix household effects for GoF check
if (fix.hhexp) {
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='household_indiv',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, sameXCycle4, name='use', interaction1='household_indiv',
initialValue=0, fix=T)
}
# fix asha and healer effects for GoF check
if (fix.oplead) {
s.effects <- setEffect(s.effects, egoX, name='talk', interaction1='ashatalk',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, altX, name='talk', interaction1='ashatalk',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='ashatalk',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, egoX, name='talk', interaction1='healertalk',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, altX, name='talk', interaction1='healertalk',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='healertalk',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, egoX, name='use', interaction1='ashatalk',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, altX, name='use', interaction1='ashause',
initialValue=0, fix=T)
s.effects <- includeInteraction(s.effects, egoX, altX, name='use', interaction1=c('ashatalk','ashause'),
fix=T)
s.effects[s.effects$effect1==2663 & s.effects$effect2==2708, 'initialValue'] <- 0
s.effects <- setEffect(s.effects, egoX, name='use', interaction1='healertalk',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, altX, name='use', interaction1='healeruse',
initialValue=0, fix=T)
s.effects <- includeInteraction(s.effects, egoX, altX, name='use', interaction1=c('healertalk','healeruse'),
fix=T)
s.effects[s.effects$effect1==2729 & s.effects$effect2==2774, 'initialValue'] <- 0
}
# fix individual effects for GoF check
if (fix.individual) {
s.effects <- setEffect(s.effects, egoX, name='talk', interaction1='female',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, altX, name='talk', interaction1='female',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='female',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, egoX, name='talk', interaction1='headhouse',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, altX, name='talk', interaction1='headhouse',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='headhouse',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, egoX, name='talk', interaction1='caresick',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, altX, name='talk', interaction1='caresick',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='caresick',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, egoX, name='talk', interaction1='workfields',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, altX, name='talk', interaction1='workfields',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='workfields',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, egoX, name='talk', interaction1='age',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, altX, name='talk', interaction1='age',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, simX, name='talk', interaction1='age',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, egoX, name='talk', interaction1='educ',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, altX, name='talk', interaction1='educ',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, simX, name='talk', interaction1='educ',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, egoX, name='use', interaction1='female',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, egoX, name='use', interaction1='headhouse',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, egoX, name='use', interaction1='caresick',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, egoX, name='use', interaction1='workfields',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, egoX, name='use', interaction1='age',
initialValue=0, fix=T)
s.effects <- setEffect(s.effects, egoX, name='use', interaction1='educ',
initialValue=0, fix=T)
}
# fine-tune fixing for some subset specifications
if (fix.netexp & !fix.hhexp & !fix.netsize & !fix.oplead &
!fix.individual) {
if (i==7) {
s.effects <- setEffect(s.effects, outPopIntn, name='talk', interaction1='use',
initialValue=0.25, fix=T, test=T)
}
if (i==10) {
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, egoX, name='talk', interaction1='caresick',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, altX, name='talk', interaction1='caresick',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='caresick',
initialValue=0, fix=T, test=T)
}
}
if (!fix.netexp & fix.hhexp & !fix.netsize & !fix.oplead &
!fix.individual) {
if (i==8) {
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=0.6, fix=T, test=T)
}
}
if (fix.netexp & !fix.hhexp & !fix.netsize & !fix.oplead &
!fix.individual) {
if (i==8) {
s.effects <- setEffect(s.effects, outActIntn, name='use', interaction1='talk',
initialValue=0.5, fix=T, test=T)
}
}
if (fix.netexp & fix.hhexp & !fix.netsize & !fix.oplead &
!fix.individual) {
if (i==8) {
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=0.6, fix=T, test=T)
}
}
# further fixing for rate parameter checks
if (rates==6) {
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=-0, fix=T, test=T)
s.effects <- setEffect(s.effects, outPopIntn, name='talk', interaction1='use',
initialValue=-0, fix=T, test=T)
s.effects <- setEffect(s.effects, outActIntn, name='use', interaction1='talk',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, inActIntn, name='use', interaction1='talk',
initialValue=0, fix=T, test=T)
}
if (rates==5) {
if (i==2) {
s.effects <- setEffect(s.effects, outActIntn, name='use', interaction1='talk',
initialValue=-1, fix=T, test=T)
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=-1, fix=T, test=T)
}
if (i==6) {
s.effects <- setEffect(s.effects, outActIntn, name='use', interaction1='talk',
initialValue=-1, fix=T, test=T)
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=-1, fix=T, test=T)
}
if (i==7) {
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=0.5, fix=T, test=T)
s.effects <- setEffect(s.effects, outPopIntn, name='talk', interaction1='use',
initialValue=0.5, fix=T, test=T)
s.effects <- setEffect(s.effects, outActIntn, name='use', interaction1='talk',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='household_indiv',
initialValue=2, fix=T, test=T)
}
if (i==8) {
s.effects <- setEffect(s.effects, outActIntn, name='use', interaction1='talk',
initialValue=1, fix=T, test=T)
}
if (i==9) {
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=1, fix=T, test=T)
s.effects <- setEffect(s.effects, outPopIntn, name='talk', interaction1='use',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='household_indiv',
initialValue=4, fix=T, test=T)
}
if (i==10) {
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=-1, fix=T, test=T)
s.effects <- setEffect(s.effects, egoX, name='use', interaction1='caresick',
initialValue=0.5, fix=T, test=T)
}
}
if (rates==8) {
if (i==1) {
s.effects <- setEffect(s.effects, sameX, name='talk', interaction1='household_indiv',
initialValue=2, fix=T, test=T)
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=-1, fix=T, test=T)
}
if (i==2) {
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=-1, fix=T, test=T)
s.effects <- setEffect(s.effects, outPopIntn, name='talk', interaction1='use',
initialValue=-1, fix=T, test=T)
}
if (i==3) {
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, outPopIntn, name='talk', interaction1='use',
initialValue=-0.5, fix=T, test=T)
}
if (i==4) {
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=0, fix=T, test=T)
}
if (i==10) {
s.effects <- setEffect(s.effects, outActIntn, name='talk', interaction1='use',
initialValue=0, fix=T, test=T)
s.effects <- setEffect(s.effects, outPopIntn, name='talk', interaction1='use',
initialValue=0, fix=T, test=T)
}
}
# first round of estimation
s.algo <- sienaAlgorithmCreate(projname='saom-stationary', n3=5000, maxlike=F,
firstg=0.01, nsub=6)
s.ans <- siena07(s.algo, data=s.data, effects=s.effects, returnDeps=T,
batch=T, verbose=F, useCluster=T, nbrNodes=4)
# continued estimation
s.algo <- sienaAlgorithmCreate(projname='saom-stationary', n3=5000, maxlike=F,
nsub=1, n2start=500)
s.result <- siena07ToConvergence(alg=s.algo, dat=s.data, eff=s.effects, returnDeps=T,
batch=T, verbose=F, useCluster=T, nbrNodes=4, ans0=s.ans)
# save result
if (!fix.individual & !fix.oplead & !fix.netsize & !fix.netexp & !fix.hhexp)
objname <- paste0('saom.v', i, '.rate', rates, '.', saom.models[7])
if (!fix.individual & !fix.oplead & !fix.netsize & fix.netexp & !fix.hhexp)
objname <- paste0('saom.v', i, '.rate', rates, '.', saom.models[6])
if (!fix.individual & !fix.oplead & !fix.netsize & !fix.netexp & fix.hhexp)
objname <- paste0('saom.v', i, '.rate', rates, '.', saom.models[5])
if (!fix.individual & !fix.oplead & !fix.netsize & fix.netexp & fix.hhexp)
objname <- paste0('saom.v', i, '.rate', rates, '.', saom.models[4])
if (!fix.individual & !fix.oplead & fix.netsize & fix.netexp & fix.hhexp)
objname <- paste0('saom.v', i, '.rate', rates, '.', saom.models[3])
if (!fix.individual & fix.oplead & fix.netsize & fix.netexp & fix.hhexp)
objname <- paste0('saom.v', i, '.rate', rates, '.', saom.models[2])
if (fix.individual & fix.oplead & fix.netsize & fix.netexp & fix.hhexp)
objname <- paste0('saom.v', i, '.rate', rates, '.', saom.models[1])
assign(objname, eval(s.result))
# print results to files if requested
if (print.results) siena.table(s.result, file=paste0(objname, '.html'),
type='html', sig=T, d=2)
}
rm(i, saom.talk, saom.talk.mod, saom.covars, saom.use, saom.use.mod, saom.use.asha,
saom.use.healer, saom.hh, s.nodes1, s.nodes2,
s.talk, s.use, s.hh.dyad, s.hh.indiv, s.female, s.headhouse, s.caresick, s.workfields,
s.age, s.educ, s.nonresptalk, s.outvilltalk, s.ashatalk, s.ashause, s.healertalk, s.healeruse,
s.result, s.algo, s.data, s.effects, s.ans,
who.asha, who.healer, objname)
rm(siena07ToConvergence, household.aggregate, predictive.acc, reg.backward,
GeodesicDistribution, CliqueCensus, mixedTriadCensus.agree)