In collinmmccabe/enss: Effective Network Size Simulation

library(scatterplot3d)
library(stargazer)
library(lmodel2)

Correlations between Network Size, Spread Time, and Prevalence

We tested for associations between the size of a network and both a.) time until extinction, saturation, or equilibrium of a disease as well as b.) maximum prevalence of the disease. We tested these using ordinary least squares and model 2 regression models. The relationships were not quite as tight as hoped, but do still show general trends toward predictive relationships. The results are broken up again by transmission mode.

SIS results

Graph

Effective_SIS_max <- readRDS("Effective_SIS_max")
SIS_max_df <- data.frame(days=Effective_SIS_max[[6]][1,], S_final=Effective_SIS_max[[6]][2,], I_final=Effective_SIS_max[[6]][3,], n=6)
for (i in 7:max_n) {
  SIS_max_df <- rbind(SIS_max_df, cbind(days=Effective_SIS_max[[i]][1,], S_final=Effective_SIS_max[[i]][2,], I_final=Effective_SIS_max[[i]][3,], n=i))
}
SIS_max_df <- cbind(SIS_max_df, equilibrium=(SIS_max_df$I_final / SIS_max_df$n))
scatterplot3d(data.frame(days=SIS_max_df$days, eq=SIS_max_df$equilibrium, n=SIS_max_df$n),
              main="SIS Model Outcomes\n(Results of 1000 Iterations)",
              xlab = "Outbreak Time (days)",
              ylab = "Prevalence of Infection at Equilibrium",
              zlab = "Network size (n)")

Linear Models

modSIS.1=lm(n ~ days + equilibrium, data=SIS_max_df)
modSIS.2=lm(n ~ days, data=SIS_max_df)
modSIS.3=lm(n ~ equilibrium, data=SIS_max_df)
stargazer(modSIS.1, modSIS.2, modSIS.3, header=FALSE, title = "SIS Linear Models", omit.stat="f", intercept.bottom=FALSE,  covariate.labels=c("Intercept", "Outbreak Duration (days)", "Equilibrium Number of Infected (n)"), dep.var.labels="Network Size (n)")

lmodel2(n ~ days + equilibrium, data=SIS_max_df)

These results suggest that maybe we should separate instances of disease going extinct vs reaching equilibrium, as these show distinct splits in the graph

SI results

Graph

Effective_SI_max <- readRDS("./sims_maximal/max/Effective_SI_max")
SI_max_df <- data.frame(days=Effective_SI_max[[6]][1,], S_final=Effective_SI_max[[6]][2,], I_final=Effective_SI_max[[6]][3,], n=6)
for (i in 7:max_n) {
  SI_max_df <- rbind(SI_max_df, cbind(days=Effective_SI_max[[i]][1,], S_final=Effective_SI_max[[i]][2,], I_final=Effective_SI_max[[i]][3,], n=i))
}
scatterplot3d(data.frame(days=SI_max_df$days, eq=SI_max_df$I_final, n=SI_max_df$n),
              main="SI Model Outcomes\n(Results of 1000 Iterations)",
              xlab = "Outbreak Time (days)",
              ylab = "Final Prevalence of Infection",
              zlab = "Network size (n)")

Linear Models

modSI.1=lm(n ~ days + I_final, data=SI_max_df)
modSI.2=lm(n ~ days, data=SI_max_df)
modSI.3=lm(n ~ I_final, data=SI_max_df)
stargazer(modSI.1, modSI.2, modSI.3, header=FALSE, title = "SI Linear Models", omit.stat="f", intercept.bottom=FALSE,  covariate.labels=c("Intercept", "Outbreak Duration (days)", "Final Number of Infected (n)"), dep.var.labels="Network Size (n)")

lmodel2(n ~ days + I_final, data=SI_max_df)

In the SI models, I think beta might be too high to get good estimates of spread time, since everything appears to go to saturation so quickly.

STD results

Graph

Effective_STD_max <- readRDS("./sims_maximal/max/Effective_STD_max")
STD_max_df <- data.frame(days=Effective_STD_max[[6]][1,], S_final=Effective_STD_max[[6]][2,], I_final=Effective_STD_max[[6]][3,], n=6)
for (i in 7:max_n) {
  STD_max_df <- rbind(STD_max_df, cbind(days=Effective_STD_max[[i]][1,], S_final=Effective_STD_max[[i]][2,], I_final=Effective_STD_max[[i]][3,], n=i))
}
scatterplot3d(data.frame(days=STD_max_df$days, eq=STD_max_df$I_final, n=STD_max_df$n),
              main="STD Model Outcomes\n(Results of 1000 Iterations)",
              xlab = "Outbreak Time (days)",
              ylab = "Final Prevalence of Infection",
              zlab = "Network size (n)")

Linear Models

modSTD.1=lm(n ~ days + I_final, data=STD_max_df)
modSTD.2=lm(n ~ days, data=STD_max_df)
modSTD.3=lm(n ~ I_final, data=STD_max_df)
stargazer(modSTD.1, modSTD.2, modSTD.3, header=FALSE, title = "STD Linear Models", omit.stat="f", intercept.bottom=FALSE,  covariate.labels=c("Intercept", "Outbreak Duration (days)", "Final Number of Infected (n)"), dep.var.labels="Network Size (n)")

lmodel2(n ~ days + I_final, data=STD_max_df)

...again, the STD results might indicate an issue of beta being too high.

SIR results

Graph

Effective_SIR_max <- readRDS("./sims_maximal/max/Effective_SIR_max")
SIR_max_df <- data.frame(days=Effective_SIR_max[[6]][1,], S_final=Effective_SIR_max[[6]][2,], I_final=Effective_SIR_max[[6]][3,], R_final=Effective_SIR_max[[6]][4,], I_max=Effective_SIR_max[[6]][5,], n=6)
for (i in 7:max_n) {
  SIR_max_df <- rbind(SIR_max_df, cbind(days=Effective_SIR_max[[i]][1,], S_final=Effective_SIR_max[[i]][2,], I_final=Effective_SIR_max[[i]][3,], R_final=Effective_SIR_max[[i]][4,], I_max=Effective_SIR_max[[i]][5,], n=i))
}
scatterplot3d(data.frame(days=SIR_max_df$days, eq=SIR_max_df$I_max, n=SIR_max_df$n),
              main="SIR Model Outcomes\n(Results of 1000 Iterations)",
              xlab = "Outbreak Time (days)",
              ylab = "Maximum Prevalence of Infection",
              zlab = "Network size (n)")

Linear Models

Should we ignore simulations with early extinctions, because again, there is grouping of those results on the graph? Or include proportion of iterations going extinct as another variable for prediction?

modSIR.1=lm(n ~ days + I_max, data=SIR_max_df)
modSIR.2=lm(n ~ days, data=SIR_max_df)
modSIR.3=lm(n ~ I_max, data=SIR_max_df)
stargazer(modSIR.1, modSIR.2, modSIR.3, header=FALSE, title="SIR Linear Models", omit.stat="f", intercept.bottom=FALSE,  covariate.labels=c("Intercept", "Outbreak Duration (days)", "Maximum Infected at once (n)"), dep.var.labels="Network Size (n)")

lmodel2(n ~ days + I_max, data=SIR_max_df)

collinmmccabe/enss documentation built on May 5, 2024, 6:23 a.m.