knitr::opts_chunk$set(error = TRUE, message = FALSE, warning=F)
Data: whamp Package Versions: 3 Model: hiv_caxr_
Preliminaries
Packages
startKnitTime <- Sys.time()
## To update
# if (params$version == 3) {
# source(here("install_master_pkgs.R"))
# } else {
# source(here("install_dev_pkgs.R"))
# }
defaultPaths <- .libPaths()
if(params$version == "3") {
pkglib <- "~/R/GHmaster-library" # ergm 3.11
print("Github Master package versions")
} else if(params$version == "4") {
pkglib <- "~/R/GHdev-library" # ergm 4.0
print("Github Dev package versions")
} else {
stop("You need to specify the R library for the statnet packages")
}
## [1] "Github Master package versions"
.libPaths(c(pkglib, defaultPaths))
.libPaths()
## [1] "C:/Users/Martina Morris/Documents/R/GHmaster-library"
## [2] "C:/Users/Martina Morris/Documents/R/win-library/4.0"
## [3] "C:/Program Files/R/R-4.0.4/library"
devtools::dev_mode(on = T, path = pkglib)
library(ergm.ego, lib.loc = pkglib)
library(EpiModel, lib.loc = pkglib)
sessioninfo::package_info(pkgs = c("network","statnet.common","ergm","tergm",
"ergm.ego","EpiModel"),
dependencies = FALSE)
## package * version date lib
## EpiModel * 2.0.4 2021-02-04 [1]
## ergm * 3.11.0-6010 2021-02-04 [1]
## ergm.ego * 0.6.0-572 2021-02-04 [1]
## network * 1.17.0-666 2021-02-04 [1]
## statnet.common 4.5.0-327 2021-02-04 [1]
## tergm * 3.7.0-2083 2021-02-04 [1]
## source
## Github (statnet/EpiModel@5ca1ef5)
## Github (statnet/ergm@c42b431)
## Github (statnet/ergm.ego@46415e5)
## Github (statnet/network@726c01d)
## Github (statnet/statnet.common@964889d)
## Github (statnet/tergm@4366856)
##
## [1] C:/Users/Martina Morris/Documents/R/GHmaster-library
## [2] C:/Users/Martina Morris/Documents/R/win-library/4.0
## [3] C:/Program Files/R/R-4.0.4/library
# these can be CRAN versions
library(tidyverse)
library(here)
library(latticeExtra)
library(devtools)
# use this for debugging if needed.
#trace(ergm, quote(save(list=ls(), file="ergm_dump.rda")))
if(!dir.exists(here("Data","stergmFits"))){
dir.create(here("Data", "stergmFits"))}
if(!dir.exists(here("Data", "stergmFits","TestingObjects"))){
dir.create(here("Data", "stergmFits","TestingObjects"))}
if(!dir.exists(here("Data", "stergmFits"))){
dir.create(here("Data", "stergmFits"))}
Settings and Data processing
## Set # replicates
reps = 1
## Set ppop size and construction
ppop = 15000
ppopwt = 'round' # Round gives consistent netsize/composition
# Read in data:
# this script is now set to work with WHAMP data -- to use with
# ARTnet, lots of stuff would need to be checked.
# but params$df IS still SET AT THE TOP OF THE DOCUMENT IN THE YAML
# We build epistats first, before filtering for ongoing==1 in
# onglingAlters, so we can calculate the duration estimates for
# inactive ties in epistats.
egodata <- WHAMPData::whampArtnetWA_egodata
source(here("STERGMs","MM","prepData.R")) # keeps inactive alters
## [1] "Before: nodecov(~ age > 66) for all nets"
## nodecov.age>66
## 3.908362
## nodecov.age>66
## 3.559936
## nodecov.age>66
## 2.393845
## [1] "After: nodecov(~ age > 66) for all nets"
## nodecov.age>66
## 0
## nodecov.age>66
## 0
## nodecov.age>66
## 0
## [1] "age.grp tables for egos and all alters"
## ego_age.grp
## 1 2 3 4 5
## 126 252 169 148 137
## mainalter_age.grp
## 1 2 3 4 5
## 48 116 66 33 30
## caslalter_age.grp
## 1 2 3 4 5
## 134 361 260 178 104
## instalter_age.grp
## 1 2 3 4 5
## 196 390 185 100 55
## [1] "main egos 832 ; main alters 293 ; egoWt vector 832"
## [1] "casl egos 832 ; casl alters 1037 ; egoWt vector 832"
## [1] "inst egos 832 ; inst alters 926 ; egoWt vector 832"
epistats <- egonet::build_epistats(egodata = egodata)
## [1] "probability of condom use by main/casl: "
##
## 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
## 1 0.787 0.055 0.043 0.017 0.013 0.021 0.004 0.004 0.013 0.000 0.043
## 2 0.709 0.058 0.049 0.022 0.003 0.058 0.000 0.012 0.014 0.001 0.074
source(here("STERGMs","MM","ongoingAlters.R")) # filters out inactive
## [1] "Alters restricted to ongoing?"
## main
## 1 <NA>
## 244 0
## casl
## 1 <NA>
## 742 0
## Use main and casl objects returned from above
main <- egodata$main
casl <- egodata$casl
The WHAMP survey data sample size is 682 cases. For this analysis we also add 150 WA state resident cases from the ARTnet survey, for a total sample size of 832.
Analysis is restricted to egos who report ever having anal sex.
# Check: Roletype now imputed, so should not be missing
print("Ego roletype: missing if no anal sex partners reported last 12 mos")
## [1] "Ego roletype: missing if no anal sex partners reported last 12 mos"
table(main = main$egos$role.type, useNA = "always")
## main
## 0 1 2 <NA>
## 116 189 527 0
table(casl = casl$egos$role.type, useNA = "always")
## casl
## 0 1 2 <NA>
## 116 189 527 0
# Check: region values for alters outside WA
print("Alter region NA")
## [1] "Alter region NA"
table(main = is.na(main$alters$region))
## main
## FALSE
## 244
table(casl = is.na(casl$alters$region))
## casl
## FALSE
## 742
# Check: HIV imputes for egos & alters
print("HIV imputes: Egos")
## [1] "HIV imputes: Egos"
with(main$egos,
table(status, imputed = ifelse(is.na(imp_prob_hiv), 0, 1), useNA = "al"))
## imputed
## status 0 1 <NA>
## 0 742 1 0
## 1 89 0 0
## <NA> 0 0 0
print("HIV imputes: Main Alters")
## [1] "HIV imputes: Main Alters"
with(main$alters,
table(status = status,
imputed = ifelse(is.na(imp_prob_hiv_alter), 0, 1),
useNA = "al"))
## imputed
## status 0 1 <NA>
## 0 154 58 0
## 1 22 10 0
## <NA> 0 0 0
print("HIV imputes: Casl Alters")
## [1] "HIV imputes: Casl Alters"
with(casl$alters,
table(status = status,
imputed = ifelse(is.na(imp_prob_hiv_alter), 0, 1),
useNA = "al"))
## imputed
## status 0 1 <NA>
## 0 360 285 0
## 1 52 45 0
## <NA> 0 0 0
# Check that num egos is the same for both networks
if (dim(main$egos)[[1]] != dim(casl$egos)[[1]]) {
print("Number of egos in main and casl egodata objects not equal")
}
Look at the mixing matrices (they’re sparse)
writeLines("Mixingmatrix by race: \n\nMain partnerships")
## Mixingmatrix by race:
##
## Main partnerships
round(mixingmatrix(main, "race"),0)
## alter
## ego B H O
## B 3 4 8
## H 0 5 20
## O 11 18 172
writeLines("\nCasual partnerships")
##
## Casual partnerships
round(mixingmatrix(casl, "race"),0)
## alter
## ego B H O
## B 9 11 51
## H 14 17 59
## O 45 73 482
writeLines("\nMixingmatrix by region: \n\nMain partnerships")
##
## Mixingmatrix by region:
##
## Main partnerships
round(mixingmatrix(main, "region"),0)
## alter
## ego EasternWA King WesternWA
## EasternWA 16 1 5
## King 2 138 7
## WesternWA 0 6 67
writeLines("\nCasual partnerships")
##
## Casual partnerships
round(mixingmatrix(casl, "region"),0)
## alter
## ego EasternWA King WesternWA
## EasternWA 51 10 3
## King 7 397 25
## WesternWA 2 26 240
writeLines("\nMixingmatrix by age group: \n\nMain partnerships")
##
## Mixingmatrix by age group:
##
## Main partnerships
round(mixingmatrix(main, "age.grp"),0)
## alter
## ego 1 2 3 4 5
## 1 20 14 4 1 0
## 2 16 44 23 1 2
## 3 1 10 16 11 13
## 4 3 8 9 10 4
## 5 0 5 7 9 10
writeLines("\nCasual partnerships")
##
## Casual partnerships
round(mixingmatrix(casl, "age.grp"),0)
## alter
## ego 1 2 3 4 5
## 1 19 29 7 0 0
## 2 18 78 32 14 6
## 3 10 70 77 42 10
## 4 9 40 43 47 28
## 5 10 39 50 46 38
Model specification (main/casl)
# race*age interaction variables created in prepData.R
# ex: age.B = ifelse(race=="B", age, 0)
model_main <- main ~ edges +
#nodecov("age") + # using race interactions with age
#nodecov(~ age^2) +
nodecov(~ age < 16) + # entering cohort
nodefactor("race", levels = -3) + # White/Other is refcat
nodefactor("region", levels = -2) + # King is refcat
absdiff("sqrt.age") +
nodematch("race", diff=TRUE) +
nodematch("region", diff=TRUE) +
nodecov("age.B") +
nodecov("age.H") +
nodecov("age.O") +
nodecov("age2.B") +
nodecov("age2.H") +
nodecov("age2.O") +
concurrent +
nodefactor(~ deg.casl > 0) +
nodefactor("status") +
nodematch("status") +
offset(nodematch("role.type",
diff = TRUE,
levels = c(1, 2))) +
offset(nodecov(~ age > 66)) # inactive agegrp 6
offset_main <- rep(-Inf, 3)
# race*age^2 interaction variables created in prepData.R
# ex: age2.B = ifelse(race=="B", age*age, 0)
model_casl <- casl ~ edges +
# nodecov("age") + # using race interactions with age
# nodecov(~ age^2) +
nodecov(~ age < 16) + # entering cohort
nodefactor("race", levels = -3) + # White/Other is refcat
nodefactor("region", levels = -2) + # King is refcat
absdiff("sqrt.age") +
nodematch("race", diff=TRUE) +
nodematch("region", diff=TRUE) +
nodecov("age.B") +
nodecov("age.H") +
nodecov("age.O") +
nodecov("age2.B") +
nodecov("age2.H") +
nodecov("age2.O") +
concurrent +
nodefactor(~ deg.main > 0) +
nodefactor("status") +
nodematch("status") +
offset(nodematch("role.type",
diff = TRUE,
levels = c(1, 2))) +
offset(nodecov(~ age > 66)) # inactive agegrp 6
offset_casl <- rep(-Inf, 3)
Target Stats
Note: this uses simple weighted scaling. It should be consistent with the ppop-based network if the ppop.wt = round option is used, but may not be consistent if the ppop.wt = sample option is used.
Main
summary(model_main, scaleto = ppop)
## edges nodecov.age<16
## 2.176323e+03 6.009150e+00
## nodefactor.race.B nodefactor.race.H
## 2.581043e+02 4.722021e+02
## nodefactor.region.EasternWA nodefactor.region.WesternWA
## 3.563686e+02 1.371316e+03
## absdiff.sqrt.age nodematch.race.B
## 1.443308e+03 2.366771e+01
## nodematch.race.H nodematch.race.O
## 4.555126e+01 1.553413e+03
## nodematch.region.EasternWA nodematch.region.King
## 1.421300e+02 1.244442e+03
## nodematch.region.WesternWA nodecov.age.B
## 6.082933e+02 9.015530e+03
## nodecov.age.H nodecov.age.O
## 1.719014e+04 1.352304e+05
## nodecov.age2.B nodecov.age2.H
## 3.358380e+05 6.721190e+05
## nodecov.age2.O concurrent
## 5.633996e+06 2.274761e+02
## nodefactor.deg.casl>0.TRUE nodefactor.status.1
## 1.827319e+03 5.394735e+02
## nodematch.status offset(nodematch.role.type.0)
## 1.881917e+03 0.000000e+00
## offset(nodematch.role.type.1) offset(nodecov.age>66)
## 0.000000e+00 0.000000e+00
Casl
summary(model_casl, scaleto = ppop)
## edges nodecov.age<16
## 6.858195e+03 6.009150e+00
## nodefactor.race.B nodefactor.race.H
## 1.244734e+03 1.724887e+03
## nodefactor.region.EasternWA nodefactor.region.WesternWA
## 1.122984e+03 4.835527e+03
## absdiff.sqrt.age nodematch.race.B
## 5.595518e+03 7.830091e+01
## nodematch.race.H nodematch.race.O
## 1.563387e+02 4.346697e+03
## nodematch.region.EasternWA nodematch.region.King
## 4.638387e+02 3.574675e+03
## nodematch.region.WesternWA nodecov.age.B
## 2.165382e+03 4.960013e+04
## nodecov.age.H nodecov.age.O
## 5.892221e+04 4.569768e+05
## nodecov.age2.B nodecov.age2.H
## 2.085214e+06 2.193525e+06
## nodecov.age2.O concurrent
## 2.098579e+07 3.481179e+03
## nodefactor.deg.main>0.TRUE nodefactor.status.1
## 3.022959e+03 1.866007e+03
## nodematch.status offset(nodematch.role.type.0)
## 5.844100e+03 0.000000e+00
## offset(nodematch.role.type.1) offset(nodecov.age>66)
## 0.000000e+00 0.000000e+00
Control settings
Parallel setup
This is system dependent, so may need tweaking for your setup. To disable parallel set np = 0 below.
# Parallel control parameters:
## system dependent
## don't use all cores
## default parallel.type is PSOCK
## to use/disable parallel, toggle comment the two np lines
np <- max(0, parallel::detectCores() - 4)
# np <- 0
require(parallel)
if(Sys.info()["sysname"] == "Linux") {
cl <- makeCluster(np)
clusterExport(cl, "pkglib")
clusterEvalQ(cl, {
.libPaths(c(pkglib, .libPaths()))
library(ergm.ego)
sessionInfo()
})
} else {
if (np==0) ptype=NULL
}
Estimation controls
if(params$version == "3") {
ergm.controls <- control.ergm(MCMC.interval = 1e5,
parallel = np)
} else {
## proposal tweaking for MCMC and san in 4.0:
## BD for role.type (like sex in het models)
## fmat indicates which mixing types of attr are forbidden
# > network::mixingmatrix(egodata$main, "role.type")
# alter
# ego 0 1 2
# 0 0.000000 25.06249 3.82385
# 1 15.320714 0.00000 27.06031
# 2 5.613572 13.16308 161.69481
fmatrix <- matrix(c(1,0,0,
0,1,0,
0,0,0), 3, 3)
## stratify by race (since B is small group)
mcmc.prop.args <- list(BD_attr = "role.type",
fmat = fmatrix,
Strat_attr = "race",
empirical=TRUE)
san.prop.args <- list(attr="role.type",
fmat = fmatrix)
san.control <- control.san(SAN.prop.weights="BDTNT",
SAN.prop.args = san.prop.args,
SAN.nsteps = 2e7
)
# Add all of these controls to control.ergm
ergm.controls <- control.ergm(MCMLE.effectiveSize=NULL,
parallel = cl,
MCMC.prop.weights="BDStratTNT",
MCMC.prop.args = mcmc.prop.args,
SAN.control = san.control,
init.MPLE.samplesize = 2e6,
#MCMC.interval = 1e5,
#MCMC.burnin = 1e6,
#MCMLE.maxit = 4000
)
}
Fits
If these fits are run with the default control=control.ergm.ego(ppop.wt='round') the nodesets are consistent across the three networks. If run with ppop.wt=‘sample’ you need to ensure consistency manually by modifying the control statement for casl net to pull the network from the main fit. This is now checked and set automatically.
Main
#Note the code is set up to run reps if we need to, from the `ergm.bench` tests.
control_main <- control.ergm.ego(
ppopsize = ppop,
ppop.wt = ppopwt,
ergm.control = ergm.controls
)
time.est.main <- rep(NA, reps)
# for (i in 1:reps) {
startTime <- Sys.time()
fit_main <- ergm.ego(model_main,
offset.coef = offset_main,
control = control_main #,
#verbose = T
)
## Error in ergm.getCluster(control, max(verbose - 1, 0)): Failed to attach package 'ergm' on one or more slave nodes. Make sure it's installed on or accessible from all of them and is in the library path.
time.est.main <- (Sys.time() - startTime)
# time.est.main[i] <- (Sys.time() - startTime)
# }
print("Estimation time:")
## [1] "Estimation time:"
round(time.est.main, 2)
## Time difference of 1.97 mins
summary(fit_main)
## Error in summary(fit_main): object 'fit_main' not found
Casl
# Pull the main network nodeset for consistency, always
mpop <- environment(fit_main$ergm.formula)$popnw
## Error in environment(fit_main$ergm.formula): object 'fit_main' not found
control_casl <- control.ergm.ego(
ppopsize = mpop,
ppop.wt = ppopwt,
ergm.control = ergm.controls
)
## Error in control.ergm.ego(ppopsize = mpop, ppop.wt = ppopwt, ergm.control = ergm.controls): object 'mpop' not found
print(paste("Using ppop.wt =", control_casl$ppop.wt,
"and ppopsize =",
ifelse(class(control_casl$ppopsize)=='numeric',
ppop, "previous network")
))
## Error in paste("Using ppop.wt =", control_casl$ppop.wt, "and ppopsize =", : object 'control_casl' not found
time.est.casl <- rep(NA, reps)
# for (i in 1:reps) {
startTime <- Sys.time()
fit_casl <- ergm.ego(model_casl,
offset.coef = offset_casl,
control = control_casl #,
#verbose = T
)
## Error in get("control", pos = parent.frame()): object 'control_casl' not found
time.est.casl <- (Sys.time() - startTime)
# time.est.casl[i] <- (Sys.time() - startTime)
# }
print("Estimation time:")
## [1] "Estimation time:"
round(time.est.casl, 2)
## Time difference of 0 secs
summary(fit_casl)
## Error in summary(fit_casl): object 'fit_casl' not found
Save the fits
version <- substr(fit_main$ergm_version, 1, 1) # to be really sure
## Error in substr(fit_main$ergm_version, 1, 1): object 'fit_main' not found
if(version != params$version) {stop("Ergm versions inconsistent")}
## Error in eval(expr, envir, enclos): Ergm versions inconsistent
save(fit_main, fit_casl,
file = here("Data", "stergmFits",
paste0("fits_mc_", params$modelout,
params$df, params$version, ".rda")),
compress = TRUE)
## Error in save(fit_main, fit_casl, file = here("Data", "stergmFits", paste0("fits_mc_", : objects 'fit_main', 'fit_casl' not found
MCMC diagnostics
These plots should be used to identify any convergence issues. They are from the last iteration of the MCMC estimation.
Main
startTime <- Sys.time()
mcmc.diagnostics(fit_main, which="plots")
## Error in mcmc.diagnostics(fit_main, which = "plots"): object 'fit_main' not found
time.mcmc.main <- (Sys.time() - startTime)
print("MCMC time:")
## [1] "MCMC time:"
round(time.mcmc.main, 2)
## Time difference of 0 secs
Casual
startTime <- Sys.time()
mcmc.diagnostics(fit_casl, which="plots")
## Error in mcmc.diagnostics(fit_casl, which = "plots"): object 'fit_casl' not found
time.mcmc.casl <- (Sys.time() - startTime)
print("MCMC time:")
## [1] "MCMC time:"
round(time.mcmc.casl, 2)
## Time difference of 0 secs
GOF
We check the model terms here to see if they reproduce the observed target stats in post-estimation simulation.
Main
startTime <- Sys.time()
plot(gof(fit_main),
control.gof.ergm(parallel = np))
## Error in gof(fit_main): object 'fit_main' not found
time.gof.main <- (Sys.time() - startTime)
print("GOF time:")
## [1] "GOF time:"
round(time.gof.main, 2)
## Time difference of 0.01 secs
Casual
startTime <- Sys.time()
plot(gof(fit_casl),
control.gof.ergm(parallel = np))
## Error in gof(fit_casl): object 'fit_casl' not found
time.gof.casl <- (Sys.time() - startTime)
print("GOF time:")
## [1] "GOF time:"
round(time.gof.casl, 2)
## Time difference of 0 secs
Network consistency Checks
Nodeset consistency
orig_egonet_main <- ergm.ego::as.egodata(fit_main$network)$egos
## Error in ergm.ego::as.egodata(fit_main$network): object 'fit_main' not found
orig_egonet_casl <- ergm.ego::as.egodata(fit_casl$network)$egos
## Error in ergm.ego::as.egodata(fit_casl$network): object 'fit_casl' not found
print("nodeset consistency")
## [1] "nodeset consistency"
all.equal(orig_egonet_main, orig_egonet_casl)
## Error in all.equal(orig_egonet_main, orig_egonet_casl): object 'orig_egonet_main' not found
Edge consistency between targetstats and ppop values
Main
index <- ifelse(params$version == "3", 2, 1)
print(rbind("Target edges"= round(fit_main$target.stats[[index]]),
"Edges in starting network"= round(network.edgecount(fit_main$network)),
"Edges in final network" = round(network.edgecount(fit_main$newnetwork))))
## Error in rbind(`Target edges` = round(fit_main$target.stats[[index]]), : object 'fit_main' not found
Casual
print(rbind("Target edges"= round(fit_casl$target.stats[[2]]),
"Edges in starting network"= round(network.edgecount(fit_casl$network)),
"Edges in final network" = round(network.edgecount(fit_casl$newnetwork))))
## Error in rbind(`Target edges` = round(fit_casl$target.stats[[2]]), `Edges in starting network` = round(network.edgecount(fit_casl$network)), : object 'fit_casl' not found
Calculate dissolution coefs and convert to netest
dr <- epistats$init.dr$wtd.mean
print(paste("Using initial death rate of", dr))
## [1] "Using initial death rate of 6.59993247584317e-05"
# Use the wtd medians for the relage estimate -- outliers and a wonky
# upper tail for casl gives it almost the same mean as main partners
main_mean_relage <-
as.numeric(1/(1 - (2^(-1/(epistats$relage.main$wtd.median)))))
casl_mean_relage <-
as.numeric(1/(1 - (2^(-1/(epistats$relage.casl$wtd.median)))))
# Calculate the dissolution coefs
diss_coef_main <-
EpiModel::dissolution_coefs(dissolution = ~offset(edges),
duration = main_mean_relage,
d.rate = epistats$init.dr$wtd.mean)
diss_coef_casl <- EpiModel::dissolution_coefs(dissolution = ~offset(edges),
duration = casl_mean_relage,
d.rate = epistats$init.dr$wtd.mean)
debug(dissolution_coefs)
# Construct EpiModel netest objects using ego.netest
netest_main <- egonet::ego.netest(fit_main, diss_coef_main)
## Error in egonet::ego.netest(fit_main, diss_coef_main): object 'fit_main' not found
netest_casl <- egonet::ego.netest(fit_casl, diss_coef_casl)
## Error in egonet::ego.netest(fit_casl, diss_coef_casl): object 'fit_casl' not found
netest_out <- list("main" = netest_main,
"casl" = netest_casl)
## Error in eval(expr, envir, enclos): object 'netest_main' not found
# Save the epsitats and netest objects
# epistats is version free
saveRDS(epistats,
file = here("Data", "Params",
paste0("epistats_", params$df, ".rds")))
saveRDS(netest_out,
file = here("Data", "stergmFits",
paste0("netest_mc_",
params$modelout,
params$df,
params$version, ".rds")))
## Error in saveRDS(netest_out, file = here("Data", "stergmFits", paste0("netest_mc_", : object 'netest_out' not found
#### Het diss by agegrp -- maybe someday...########################
# main_mean_relage_by_agegrp <-
# as.numeric(1/(1 - (2^(-1/(epistats$relage_by_agegrp.main$wtd.median)))))
# casl_mean_relage_by_agegrp <-
# as.numeric(1/(1 - (2^(-1/(epistats$relage_by_agegrp.casl$wtd.median)))))
# diss_coef_main_by_agegrp <-
# EpiModel::dissolution_coefs(dissolution = ~offset(edges) + offset(nodefactor(age.grp)),
# duration = main_mean_relage_by_agegrp,
# d.rate = epistats$init.dr$wtd.mean)
#
# diss_coef_casl_by_agegrp <-
# EpiModel::dissolution_coefs(dissolution = ~offset(edges) + offset(nodefactor(age.grp)),
# duration = casl_mean_relage_by_agegrp,
# d.rate = epistats$init.dr$wtd.mean)
# netest_main_dagegrp <- egonet::ego.netest(fit_main, diss_coef_main_by_agegrp)
# netest_casl_dagegrp <- egonet::ego.netest(fit_casl, diss_coef_casl_by_agegrp)
# netest_out_dagegrp <- list("main" = netest_main_dagegrp,
# "casl" = netest_casl_dagegrp)
Static Netdx
Static netdx are single draws from the equilibrium distribution of the static graphs. There are no edge dynamics.
Main
startTime <- Sys.time()
dx <- netdx(netest_main, nsims = 50, dynamic = F,
ncores = np)
## Error in netdx(netest_main, nsims = 50, dynamic = F, ncores = np): object 'netest_main' not found
time.netdxs.main <- (Sys.time() - startTime)
print("netdx static simulation time:")
## [1] "netdx static simulation time:"
round(time.netdxs.main, 2)
## Time difference of 0.01 secs
print(dx)
## Error in print(dx): object 'dx' not found
plot(dx, method = "b")
## Error in plot(dx, method = "b"): object 'dx' not found
plot(dx, mean.col = "slateblue", sim.lines = TRUE, plots.joined = FALSE)
## Error in plot(dx, mean.col = "slateblue", sim.lines = TRUE, plots.joined = FALSE): object 'dx' not found
Casl
startTime <- Sys.time()
dx <- netdx(netest_casl, nsims = 50, dynamic = F,
ncores = np)
## Error in netdx(netest_casl, nsims = 50, dynamic = F, ncores = np): object 'netest_casl' not found
time.netdxs.casl <- (Sys.time() - startTime)
print("netdx static simulation time:")
## [1] "netdx static simulation time:"
round(time.netdxs.casl, 2)
## Time difference of 0.01 secs
print(dx)
## Error in print(dx): object 'dx' not found
plot(dx, method = "b")
## Error in plot(dx, method = "b"): object 'dx' not found
plot(dx, mean.col = "slateblue", sim.lines = TRUE, plots.joined = FALSE)
## Error in plot(dx, mean.col = "slateblue", sim.lines = TRUE, plots.joined = FALSE): object 'dx' not found
Dynamic
Note the dynamic netdx are based on a closed population – there are no demographic dynamics. The dx are used to ensure that the targetstats are hit when edge dynamics are included (partnership formation/dissolution).
Main
startTime <- Sys.time()
dx <- netdx(netest_main, nsims = 5, dynamic = T, nsteps = 1000,
set.control.ergm =
control.simulate.ergm(MCMC.burnin = 1e5),
set.control.stergm =
control.simulate.network(MCMC.burnin.min = 1e4),
ncores = np)
## Error in netdx(netest_main, nsims = 5, dynamic = T, nsteps = 1000, set.control.ergm = control.simulate.ergm(MCMC.burnin = 1e+05), : object 'netest_main' not found
time.netdxd.main <- (Sys.time() - startTime)
print("netdx dynamic simulation time:")
## [1] "netdx dynamic simulation time:"
round(time.netdxd.main, 2)
## Time difference of 0 secs
print(dx)
## Error in print(dx): object 'dx' not found
plot(dx, method = "b")
## Error in plot(dx, method = "b"): object 'dx' not found
plot(dx)
## Error in plot(dx): object 'dx' not found
plot(dx, type = "duration")
## Error in plot(dx, type = "duration"): object 'dx' not found
plot(dx, type = "dissolution")
## Error in plot(dx, type = "dissolution"): object 'dx' not found
Casl
startTime <- Sys.time()
dx <- netdx(netest_casl, nsims = 5, dynamic = T, nsteps = 500,
set.control.stergm =
control.simulate.network(MCMC.burnin.min = 1e5),
ncores = np)
## Error in netdx(netest_casl, nsims = 5, dynamic = T, nsteps = 500, set.control.stergm = control.simulate.network(MCMC.burnin.min = 1e+05), : object 'netest_casl' not found
time.netdxd.casl <- (Sys.time() - startTime)
print("netdx dynamic simulation time:")
## [1] "netdx dynamic simulation time:"
round(time.netdxd.casl, 2)
## Time difference of 0 secs
print(dx)
## Error in print(dx): object 'dx' not found
plot(dx, method = "b")
## Error in plot(dx, method = "b"): object 'dx' not found
plot(dx)
## Error in plot(dx): object 'dx' not found
plot(dx, type = "duration")
## Error in plot(dx, type = "duration"): object 'dx' not found
plot(dx, type = "dissolution")
## Error in plot(dx, type = "dissolution"): object 'dx' not found
Computation times by section
knitTime <- Sys.time() - startKnitTime
mainTimes <- c(time.est.main, time.mcmc.main, time.gof.main,
time.netdxs.main, time.netdxd.main) #, time.netdxd.main.agegrp)
caslTimes <- c(time.est.casl, time.mcmc.casl, time.gof.casl,
time.netdxs.casl, time.netdxd.casl) #, time.netdxd.casl.agegrp)
tab <- data.frame(main = round(as.numeric(mainTimes, units="mins"), 2),
casl = round(as.numeric(caslTimes, units="mins"), 2))
rownames(tab) <- c("Estimation", "MCMC", "GOF", "Netdx Static",
"Netdx Dynamic") #, "Netdx Dynamic Agegrp")
tab %>%
DT::datatable(
options = list(
columnDefs = list(
list(className = 'dt-right', targets = 1:2)
)),
caption = paste("Computation Times in Minutes using version",
params$version))
print(paste("Document Knit time: ", round(knitTime, 2)))
## [1] "Document Knit time: 2.15"
sessioninfo::package_info(pkgs = c("network","statnet.common","ergm","tergm",
"ergm.ego","EpiModel"),
dependencies = FALSE)
## package * version date lib
## EpiModel * 2.0.4 2021-02-04 [1]
## ergm * 3.11.0-6010 2021-02-04 [1]
## ergm.ego * 0.6.0-572 2021-02-04 [1]
## network * 1.17.0-666 2021-02-04 [1]
## statnet.common 4.5.0-327 2021-02-04 [1]
## tergm * 3.7.0-2083 2021-02-04 [1]
## source
## Github (statnet/EpiModel@5ca1ef5)
## Github (statnet/ergm@c42b431)
## Github (statnet/ergm.ego@46415e5)
## Github (statnet/network@726c01d)
## Github (statnet/statnet.common@964889d)
## Github (statnet/tergm@4366856)
##
## [1] C:/Users/Martina Morris/Documents/R/GHmaster-library
## [2] C:/Users/Martina Morris/Documents/R/win-library/4.0
## [3] C:/Program Files/R/R-4.0.4/library