## For other simulations, just use a different file here.
## However, do keep the name sim.
suppressPackageStartupMessages(library(tidyverse))
options(dplyr.summarise.inform = FALSE)
## Set path to model run you want to get diagnostics for:
# sim <- readRDS("../../../Data/EpiModelSims/ergm3_ca_nohiv_sim_with15yo_boost.rds")
if (params$sim == "None") {
sim <- readRDS("../../../EpiModel/AE/sim_epimodel3/sim_at_2033_jan12.rds")
}else{
sim <- params$sim
}
## Registered S3 method overwritten by 'tergm':
## method from
## simulate_formula.network ergm
## Creating the netest object (make sure it is the same as what was
## used to run the simulation).
## No HIV
# netest_mc <- readRDS("../../../Data/stergmFits/netest_mc_nohiv_ca_whamp3.rds")
# load("../../../Data/stergmFits/netest_inst_nohiv_whamp3.rda")
## HIV
if (params$netest == "None") {
netest_mc <- readRDS("../../../Data/stergmFits/netest_mc_hiv_caxr_whamp3.rds")
load("../../../Data/stergmFits/netest_inst_hiv_caxr_whamp3.rda")
netest_obj <- list("main" = netest_mc$main,
"casl" = netest_mc$casl,
"inst" = netest_inst)
}else{
netest_obj <- params$netest
}
make_mean_deg_plot <- function(sim, var_type = "race",
nest_obj = NULL, simnum = 1){
edge_counts <- sim$epi[[simnum]]$demog_edge_counts
edge_counts$at[edge_counts$at == 1] <-
min(edge_counts$at[edge_counts$at != 1]) - 1
edge_counts[, "col_var"] <- edge_counts[, var_type]
f_edg_cnts <- edge_counts %>% filter(ptype == 0) %>%
mutate(rep = 3) %>% uncount(rep, .id = "ptype")
edge_counts <- bind_rows(
edge_counts %>% filter(ptype != 0),
f_edg_cnts
)
cmpr_edges <- pivot_wider(
edge_counts,
values_from = "count", names_from = "type",
id_cols = colnames(edge_counts)
) %>% filter(age.grp != 6) %>%
group_by(ptype, col_var, at) %>%
summarise(partner = sum(partner, na.rm = TRUE),
total = sum(total, na.rm = TRUE)) %>%
mutate(
"mean_deg" = partner / total,
"year" = round(
(at - sim$control$start) / 52 + sim$control$year_start, 4
)
)
cmpr_edges$ptype <- factor(cmpr_edges$ptype, levels = 1:3,
labels = c("Main", "Casl", "Inst"))
if (is.null(nest_obj)) {
cmpr_edges %>% ungroup() %>%
ggplot(aes(x = year, y = mean_deg,
color = as.character(col_var))) +
geom_point(alpha = 0.5) + geom_smooth() +
facet_wrap(~ ptype, scales = "free_y" ) +
geom_vline(xintercept = 2019, color = "black") +
scale_colour_discrete(var_type) + ylab("Mean Degree") +
xlab("Year") + theme(legend.position = "bottom")
}else{
degs <- lapply(nest_obj, "[[", "fit") %>% lapply("[[", "newnetwork") %>%
sapply(sna::degree)
ego_dat <- ergm.ego::as.egodata(nest_obj$main$fit$newnetwork)$egos
ego_dat <- data.frame(ego_dat, degs) %>%
filter(age.grp != 6) %>%
select(age.grp, race, region,
"Main" = main, "Casl" = casl, "Inst" = inst)
ego_dat[, "col_var"] <- ego_dat[, var_type]
long_dat <- pivot_longer(
ego_dat, names_to = "ptype", values_to = "deg",
# id_cols = c("age.grp", "race", "region", "col_var"),
cols = c("Main", "Casl", "Inst")
)
mean_degs <- long_dat %>% group_by(ptype, col_var) %>%
summarise(Mean_deg = mean(deg) / 2)
mean_degs$ptype <- factor(mean_degs$ptype,
levels = c("Main", "Casl", "Inst"),
labels = c("Main", "Casl", "Inst"))
cmpr_edges %>% ungroup() %>%
ggplot(aes(x = year, y = mean_deg,
color = as.character(col_var))) +
geom_point(alpha = 0.5) + geom_smooth() +
geom_hline(data = mean_degs, aes(yintercept = Mean_deg,
color = as.factor(col_var))) +
facet_wrap(~ ptype, scales = "free_y" ) +
geom_vline(xintercept = 2019, color = "black") +
scale_colour_discrete(var_type) + ylab("Mean Degree") +
xlab("Year") + theme(legend.position = "bottom")
}
}
make_mean_agextime <- function(sim, var_type = "race",
yrs = NULL, targ_year = 2019,
nest_obj = NULL,
plot_dif = FALSE){
if (length(targ_year) != 1) {stop("Must specify one target year")}
edge_counts <- sim$epi[[1]]$demog_edge_counts
if (var_type == "none") {
edge_counts[, "facet_var"] <- "pooled"
}else{
edge_counts[, "facet_var"] <- edge_counts[, var_type]
}
f_edg_cnts <- edge_counts %>% filter(ptype == 0) %>%
mutate(rep = 3) %>% uncount(rep, .id = "ptype")
edge_counts <- bind_rows(
edge_counts %>% filter(ptype != 0),
f_edg_cnts
)
cmpr_edges <- pivot_wider(
edge_counts,
values_from = "count", names_from = "type",
id_cols = colnames(edge_counts)
) %>% group_by(ptype, age, at, facet_var) %>%
summarise(partner = sum(partner, na.rm = TRUE),
total = sum(total, na.rm = TRUE)) %>%
mutate(
"mean_deg" = partner / total,
"year" = round(
(at - sim$control$start) / 52 + sim$control$year_start
)
)
cmpr_edges$ptype <- factor(cmpr_edges$ptype, levels = 1:3,
labels = c("Main", "Casl", "Inst"))
all_yrs <- sort(unique(cmpr_edges$year), decreasing = FALSE)[-1]
if (is.null(yrs)) {
keep_yrs <- c(min(all_yrs), all_yrs[all_yrs %% 10 == 0])
}else{
keep_yrs <- yrs
}
keep_yrs <- unique(c(yrs, targ_year))
cmpr_edges <- cmpr_edges %>% filter(year %in% keep_yrs,
year < 2031)
cmpr_edges$year <- factor(
cmpr_edges$year,
levels = sort(cmpr_edges$year),
labels = sort(cmpr_edges$year)
)
if (is.null(nest_obj)) {
cmpr_edges %>% ungroup() %>%
ggplot(aes(x = age, y = mean_deg, col = facet_var)) +
# geom_line(alpha = 0.5) +
geom_smooth(se = FALSE, span = 4) + #, size = 2) +
facet_grid(year~ptype) +
geom_vline(xintercept = 2019, color = "black") +
scale_colour_discrete("Partnership Type") +
ylab("Mean Degree") +
xlab("Age") + theme(legend.position = "bottom")
}else{
degs <- lapply(nest_obj, "[[", "fit") %>% lapply("[[", "newnetwork") %>%
sapply(sna::degree)
ego_dat <- ergm.ego::as.egodata(nest_obj$main$fit$newnetwork)$egos
ego_dat <- data.frame(ego_dat, degs) %>%
select(age, race, region,
"Main" = main, "Casl" = casl, "Inst" = inst) %>%
mutate(age = floor(age))
if (var_type == "none") {
ego_dat[, "facet_var"] <- "pooled"
}else{
ego_dat[, "facet_var"] <- ego_dat[, var_type]
}
long_dat <- pivot_longer(
ego_dat, names_to = "ptype", values_to = "deg",
cols = c("Main", "Casl", "Inst")
)
mean_degs <- long_dat %>% group_by(ptype, age, facet_var) %>%
summarise(mean_deg = mean(deg) / 2)
mean_degs$ptype <- factor(mean_degs$ptype,
levels = c("Main", "Casl", "Inst"),
labels = c("Main", "Casl", "Inst"))
r_mean_degs <- mean_degs %>% mutate(rep = length(keep_yrs)) %>%
uncount(rep, .id = "year") %>% mutate(year = keep_yrs[year]) %>%
filter(year < 2031)
r_mean_degs$year <- factor(
r_mean_degs$year,
levels = sort(r_mean_degs$year),
labels = sort(r_mean_degs$year)
)
if (var_type == "none") {alph <- 0.3}else{alph = 0}
all_md <- bind_rows(
bind_cols(cmpr_edges, "type" = "Simulated"),
bind_cols(r_mean_degs, "type" = "Start/Target (from netest$newnetwork)")
)
if (plot_dif) {
wd <- all_md %>% ungroup() %>% select(type, facet_var,
mean_deg, ptype, year, age, highlight) %>%
pivot_wider(names_from = type,
values_from = mean_deg) %>%
mutate(dif = Simulated - `Start/Target (from netest$newnetwork)` )
wd %>%
ggplot(aes(x = age, y = dif,
col = facet_var)) +
geom_hline(yintercept = 0, color = "black") +
geom_vline(xintercept = 65, color = "gray") +
geom_smooth(se = FALSE) + #, size = 2, span = 2.5) +
facet_wrap(~ptype + year, scales = "free_y", dir = "v",
ncol = 3) +
scale_colour_discrete(var_type) +
coord_cartesian(ylim = c(-0.6, 0.4)) +
ylab("Mean Degree (Observed - Target)") +
xlab("Age") + theme(legend.position = "bottom")
}else{
fac_df <- expand.grid(year = unique(all_md$year),
ptype = unique(all_md$ptype))
fac_df$highl <- as.numeric(fac_df$year == targ_year)
all_md %>% ungroup() %>%
ggplot() +
# geom_point(alpha = alph) +
geom_rect(data = fac_df %>% filter(highl == 1),
aes(fill = as.factor(highl)),
xmin = -Inf, xmax = Inf,
ymin = -Inf, ymax = Inf, alpha = 0.3) +
geom_vline(xintercept = 65, color = "gray") +
geom_smooth(aes(x = age, y = mean_deg,
color = facet_var,
linetype = type),
se = FALSE) + #, size = 2, span = 2.5) +
facet_grid(ptype ~ year, scales = "free_y") +
scale_colour_discrete(var_type) +
scale_fill_manual(values = "black") +
ylab("Mean Degree") + guides(fill = FALSE) +
xlab("Age") + theme(legend.position = "bottom") +
coord_cartesian(ylim = c(0, NA))
}
}
}
Mean Degree plots
Within race group
make_mean_deg_plot(sim, var_type = "race", nest_obj = netest_obj)
Within Region
make_mean_deg_plot(sim, var_type = "region", netest_obj)
Mean Degree across Age
Within race group
make_mean_agextime(sim, "race", yrs = c(2005, 2015, 2019, 2025, 2030),
nest_obj = netest_obj)
Within Region
make_mean_agextime(sim, "region", yrs = c(2005, 2015, 2019, 2025, 2030),
nest_obj = netest_obj)
Session Info
sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] forcats_0.5.0 stringr_1.4.0 dplyr_1.0.2 purrr_0.3.4
## [5] readr_1.3.1 tidyr_1.1.2 tibble_3.0.4 ggplot2_3.3.2
## [9] tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-148 fs_1.4.2 lubridate_1.7.9
## [4] doParallel_1.0.15 RColorBrewer_1.1-2 httr_1.4.2
## [7] rprojroot_1.3-2 tools_4.0.2 backports_1.1.10
## [10] R6_2.5.0 mgcv_1.8-31 DBI_1.1.0
## [13] lazyeval_0.2.2 colorspace_1.4-1 withr_2.3.0
## [16] tidyselect_1.1.0 compiler_4.0.2 ergm.multi_4.0-3965
## [19] cli_2.2.0 rvest_0.3.6 xml2_1.3.2
## [22] network_1.17.0-585 labeling_0.3 scales_1.1.1
## [25] DEoptimR_1.0-8 robustbase_0.93-6 digest_0.6.27
## [28] ergm.ego_0.5-471 rmarkdown_2.4 networkDynamic_0.10.2
## [31] pkgconfig_2.0.3 htmltools_0.5.0 egonet_0.0.1.0
## [34] dbplyr_1.4.4 rlang_0.4.9 readxl_1.3.1
## [37] rstudioapi_0.11 farver_2.0.3 generics_0.1.0
## [40] jsonlite_1.7.1 statnet.common_4.4.0-300 magrittr_2.0.1
## [43] Matrix_1.2-18 Rcpp_1.0.5 munsell_0.5.0
## [46] fansi_0.4.1 ape_5.4-1 lifecycle_0.2.0
## [49] stringi_1.5.3 yaml_2.2.1 MASS_7.3-51.6
## [52] grid_4.0.2 blob_1.2.1 parallel_4.0.2
## [55] crayon_1.3.4 lattice_0.20-41 haven_2.3.1
## [58] splines_4.0.2 hms_0.5.3 sna_2.5
## [61] knitr_1.30 pillar_1.4.7 EpiModel_2.0.3
## [64] codetools_0.2-16 lpSolve_5.6.15 rle_0.9.2-223
## [67] srvyr_0.3.10 reprex_0.3.0 glue_1.4.2
## [70] evaluate_0.14 trust_0.1-8 mitools_2.4
## [73] modelr_0.1.8 deSolve_1.28 vctrs_0.3.6
## [76] foreach_1.5.0 cellranger_1.1.0 gtable_0.3.0
## [79] assertthat_0.2.1 tergmLite_2.2.0 tergm_4.0.0-2013
## [82] xfun_0.18 broom_0.7.1 survey_4.0
## [85] ergm_4.0-5824 coda_0.19-4 survival_3.2-7
## [88] iterators_1.0.12 ellipsis_0.3.1 here_0.1