In noamross/age-infects: Package for manuscript and code of age-infection interaction paper
% Comparative Dynamics
% Noam Ross
% 14-08-25 10:10:34
r library(knitr); opts_chunk$set(echo=FALSE, cache=TRUE, message=FALSE, fig.width=12, fig.height=6, dpi=72)
Here I examine the comparative dynamics of age classes in 2-age-class models
of both $SI$, $SIV$, and macroparasite-style dynamic models.
The following parameters are identical between models:
- fecundity $(f)$
- transition rate $(g)$
- carrying capacity $(K)$
The models are also parameterized to have identical apparent mortality rates
at equilibrium. That is the aggregate mortality rates of all infection classes
are the same.
parms_SI = list(
n_stages = 2,
n_parasites = 1,
f = c(1, 1),
g = c(0.1, 0),
d = c(0.01, 0.01),
alpha = c(0.2, 0.2),
lamda = c(0.5, 0.5),
Beta = c(1, 1),
mu = c(0, 0),
xi = c(1, 1),
omega = c(1,1),
K = 1,
progress = 0,
times.min = 0,
times.max = 50,
times.by = 0.1,
inits = 0,
infect_vector = c(0.01, 0.01)
)
parms_SIV = parms_SI
parms_SIV$n_parasites = 2
parms_macro = parms_SI
parms_macro$n_parasites = 10
parms_SI$times.by = 0.01
library(deSolve)
library(plyr)
library(rootSolve)
library(reshape2)
library(ggplot2)
library(data.table)
library(noamtools)
library(numDeriv)
library(gridExtra)
source('../R/process_sodp.R')
source('../R/run_sodp.R')
out_SI = run_sodp(parms_SI)
SIsf = process_sodp(out_SI)
SIsf = sodp_totals(SIsf)
stats_SI = sodp_stats(SIsf)
stats_SI.m = melt(stats_SI, id.vars=c("time", "Species", "SizeClass"))
out_SIV = run_sodp(parms_SIV)
SIVsf = process_sodp(out_SIV)
SIVsf = sodp_totals(SIVsf)
stats_SIV = sodp_stats(SIVsf)
stats_SIV.m = melt(stats_SIV, id.vars=c("time", "Species", "SizeClass"))
out_macro = run_sodp(parms_macro)
macrosf = process_sodp(out_macro)
macrosf = sodp_totals(macrosf)
stats_macro = sodp_stats(macrosf)
stats_macro.m = melt(stats_macro, id.vars=c("time", "Species", "SizeClass"))
times = c(0, 5, 15, 20, 25, 30, 50)
df = data.frame(subset(macrosf, time %in% times))
df = subset(df, SizeClass != "Total" & Population > 1e-5)
df = df[order(df$time),]
df$time = factor(paste("Time =", df$time), levels = paste("Time =", unique(df$time)))
ggplot(df, aes(x=Infected, y=Population, fill=SizeClass)) +
geom_area(position="identity", alpha = 0.5) +
scale_fill_grey(labels = c("Juveniles", "Adults")) +
theme_nr +
xlab("Number of Infections") +
theme(text=element_text(family="Lato")) +
facet_wrap(~time, nrow=1)
df2 = data.frame(subset(SIVsf, time %in% times))
df2 = droplevels(subset(df2, SizeClass != "Total" & Population > 1e-5))
df2 = df2[order(df2$time),]
df2$time = factor(paste("Time =", df2$time), levels = paste("Time =", unique(df2$time)))
ggplot(df2, aes(x=as.factor(Infected), y=Population, fill=SizeClass)) +
geom_bar(stat="identity", position="dodge", alpha = 0.5) +
scale_fill_grey(labels = c("Juveniles", "Adults")) +
scale_x_discrete(labels = c("S", "I", "V")) +
theme_nr +
xlab("Disease Class") +
theme(text=element_text(family="Lato")) +
facet_wrap(~time, nrow=1)
stats = list(SI=stats_SI.m, SIV=stats_SIV.m, macro=stats_macro.m)
statsall = ldply(stats, function(x) x)
mortInf = ldply(stats, function(x) subset(x, variable %in% c("N")))
p1 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 1) +
ylab("Host Population") +
xlab(" ") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
mortInf = ldply(stats, function(x) subset(x, variable %in% c("I")))
p2 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 1) +
ylab("Infected Hosts") +
xlab("Time") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
mortInf = ldply(stats, function(x) subset(x, variable %in% c("MortInfRate")))
p3 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 0.4) +
ylab("Infected Mortality Rate") +
xlab(" ") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
g_legend <- function(p){
tmp <- ggplotGrob(p)
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)}
legend <- g_legend(p1)
lwidth <- sum(legend$width)
plotwidth = (unit(1, "npc") - lwidth)*(1/3)
grid.arrange(p1 + theme(legend.position="none"),
p2 + theme(legend.position="none"),
p3 + theme(legend.position="none"),
legend,
widths=unit.c(plotwidth, plotwidth, plotwidth, lwidth),
main = textGrob("Epidemic Dynamics: 1-infection equivalence", gp=gpar(fontsize=20,fontfamily='Lato')),
nrow=1)
parms_SIVmod = parms_SIV
aa = optim(par = 1, f = function(mod_SIV) {
parms_SIVmod$alpha = parms_SIV$alpha*mod_SIV
parms_SIVmod$lamda = parms_SIV$lamda*mod_SIV
out_SIVmod = run_sodp(parms_SIVmod)
SIVmodsf = process_sodp(out_SIVmod)
SIVmodsf = sodp_totals(SIVmodsf)
stats_SIVmod = sodp_stats(SIVmodsf)
val = subset(stats_SIVmod, time==parms_SIVmod$times.max & SizeClass=="Total")$MortInfRate
if(interactive()) print(val)
return((val - parms_SI$alpha[1])^2)
}, method="Brent", lower=0, upper=2)
mod_SIV = aa$par
parms_SIVmod$alpha = parms_SIV$alpha*mod_SIV
parms_SIVmod$lamda = parms_SIV$lamda*mod_SIV
out_SIVmod = run_sodp(parms_SIVmod)
SIVmodsf = process_sodp(out_SIVmod)
SIVmodsf = sodp_totals(SIVmodsf)
stats_SIVmod = sodp_stats(SIVmodsf)
stats_SIVmod.m = melt(stats_SIVmod, id.vars=c("time", "Species", "SizeClass"))
parms_macromod = parms_macro
bb = optim(par = 0.6118, f = function(mod_macro) {
parms_macromod$alpha = parms_macro$alpha*mod_macro
parms_macromod$lamda = parms_macro$lamda*mod_macro
out_macromod = run_sodp(parms_macromod)
macromodsf = process_sodp(out_macromod)
macromodsf = sodp_totals(macromodsf)
stats_macromod = sodp_stats(macromodsf)
val = subset(stats_macromod, time==parms_macromod$times.max & SizeClass=="Total")$MortInfRate
if(interactive()) print(val)
return((val - parms_SI$alpha[1])^2)
}, method="Brent", lower=0, upper=2)
mod_macro = bb$par
parms_macromod$alpha = parms_macro$alpha*mod_macro
parms_macromod$lamda = parms_macro$lamda*mod_macro
out_macromod = run_sodp(parms_macromod)
macromodsf = process_sodp(out_macromod)
macromodsf = sodp_totals(macromodsf)
stats_macromod = sodp_stats(macromodsf)
stats_macromod.m = melt(stats_macromod, id.vars=c("time", "Species", "SizeClass"))
stats = list(SI=stats_SI.m, SIV=stats_SIVmod.m, macro=stats_macromod.m)
statsall = ldply(stats, function(x) x)
mortInf = ldply(stats, function(x) subset(x, variable %in% c("N")))
p1 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 1) +
ylab("Host Population") +
xlab(" ") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
mortInf = ldply(stats, function(x) subset(x, variable %in% c("I")))
p2 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 1) +
ylab("Infected Hosts") +
xlab("Time") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
mortInf = ldply(stats, function(x) subset(x, variable %in% c("MortInfRate")))
p3 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 0.4) +
ylab("Infected Mortality Rate") +
xlab(" ") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
g_legend <- function(p){
tmp <- ggplotGrob(p)
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)}
legend <- g_legend(p1)
lwidth <- sum(legend$width)
plotwidth = (unit(1, "npc") - lwidth)*(1/3)
grid.arrange(p1 + theme(legend.position="none"),
p2 + theme(legend.position="none"),
p3 + theme(legend.position="none"),
legend,
widths=unit.c(plotwidth, plotwidth, plotwidth, lwidth),
main = textGrob("Epidemic Dynamics: Equilibrium Mortality Equivalence", gp=gpar(fontsize=20,fontfamily='Lato')),
nrow=1)
SI_Final = init_derivs2(parms_SI, SI=TRUE)[,2]
parms_SIVmod2 = parms_SIV
aa = optim(par = c(1,1), fn = function(mod_SIV) {
parms_SIVmod$alpha = parms_SIV$alpha*mod_SIV[1]
parms_SIVmod$lamda = parms_SIV$lamda*mod_SIV[2]
val = init_derivs2(parms_SIVmod, SI=TRUE)[,2]
if(interactive()) print(val)
return(sum(((val - SI_Final)/SI_Final)^2))
})
mod_SIV = aa$par
parms_SIVmod2$alpha = parms_SIV$alpha*mod_SIV[1]
parms_SIVmod2$lamda = parms_SIV$lamda*mod_SIV[2]
out_SIVmod2 = run_sodp(parms_SIVmod2)
SIVmod2sf = process_sodp(out_SIVmod2)
SIVmod2sf = sodp_totals(SIVmod2sf)
stats_SIVmod2 = sodp_stats(SIVmod2sf)
stats_SIVmod2.m = melt(stats_SIVmod2, id.vars=c("time", "Species", "SizeClass"))
parms_macromod2 = parms_macro
bb = optim(par = aa$par, f = function(mod_macro) {
parms_macromod$alpha = parms_macro$alpha*mod_macro[1]
parms_macromod$lamda = parms_macro$lamda*mod_macro[2]
val = init_derivs2(parms_macromod, SI=TRUE)[,2]
if(interactive()) print(val)
return(sum(((val - SI_Final)/SI_Final)^2))
})
mod_macro = bb$par
parms_macromod2$alpha = parms_macro$alpha*mod_macro[1]
parms_macromod2$lamda = parms_macro$lamda*mod_macro[2]
out_macromod2 = run_sodp(parms_macromod2)
macromod2sf = process_sodp(out_macromod2)
macromod2sf = sodp_totals(macromod2sf)
stats_macromod2 = sodp_stats(macromod2sf)
stats_macromod2.m = melt(stats_macromod2, id.vars=c("time", "Species", "SizeClass"))
stats = list(SI=stats_SI.m, SIV=stats_SIVmod2.m, macro=stats_macromod2.m)
statsall = ldply(stats, function(x) x)
mortInf = ldply(stats, function(x) subset(x, variable %in% c("N")))
p1 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 1) +
ylab("Host Population") +
xlab(" ") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
mortInf = ldply(stats, function(x) subset(x, variable %in% c("I")))
p2 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 1) +
ylab("Infected Hosts") +
xlab("Time") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
mortInf = ldply(stats, function(x) subset(x, variable %in% c("MortInfRate")))
p3 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 0.6) +
ylab("Infected Mortality Rate") +
xlab(" ") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
g_legend <- function(p){
tmp <- ggplotGrob(p)
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)}
legend <- g_legend(p1)
lwidth <- sum(legend$width)
plotwidth = (unit(1, "npc") - lwidth)*(1/3)
grid.arrange(p1 + theme(legend.position="none"),
p2 + theme(legend.position="none"),
p3 + theme(legend.position="none"),
legend,
widths=unit.c(plotwidth, plotwidth, plotwidth, lwidth),
main = textGrob("Epidemic Dynamics: Equivalent Initial Derivatives", gp=gpar(fontsize=20,fontfamily='Lato')),
nrow=1)
frac_I = 0.1
sit = subset(stats_SI, SizeClass=="Total", select=c("time", "FracInf"))
sit = sit[order(sit$time),]
timeToFracI_fn = function(sit, frac_I) {
afti = min(which(sit$FracInf >= frac_I), na.rm=TRUE)
if(is.na(afti) || afti==Inf) return(max(sit$time))
befi = min(which(sit$FracInf >= frac_I)) - 1
befx = sit$time[befi]
befy = sit$FracInf[befi]
aftx = sit$time[afti]
afty = sit$FracInf[afti]
slope = (afty-befy)/(aftx-befx)
return(befx + (frac_I - befy)/slope)
}
timeToFracI = timeToFracI_fn(sit, frac_I)
parms_SIVmod3 = parms_SIV
parms_SIVmod3$times.max = 18
parms_SIVmod3$times.by = 0.01
cc = optim(par = 1, fn = function(mod_SIV) {
parms_SIVmod3$alpha = parms_SIV$alpha*mod_SIV
parms_SIVmod3$lamda = parms_SIV$lamda*mod_SIV
out_SIVmod = run_sodp(parms_SIVmod3)
SIVmodsf = process_sodp(out_SIVmod)
SIVmodsf = sodp_totals(SIVmodsf)
stats_SIVmod = sodp_stats(SIVmodsf)
sitSIVmod = subset(stats_SIVmod, SizeClass=="Total", select=c("time", "FracInf"))
sitSIVmod = sitSIVmod[order(sitSIVmod$time),]
timeToFracISIVmod = timeToFracI_fn(sitSIVmod, frac_I)
if(interactive()) print(timeToFracISIVmod)
val = ((timeToFracISIVmod - timeToFracI)/timeToFracI)^2
return(val)
}, method="Brent", lower=0, upper=2)
mod_SIV = cc$par
parms_SIVmod3$alpha = parms_SIV$alpha*mod_SIV
parms_SIVmod3$lamda = parms_SIV$lamda*mod_SIV
parms_SIVmod3$times.max = 50
parms_SIVmod3$times.by = 0.1
out_SIVmod3 = run_sodp(parms_SIVmod3)
SIVmod3sf = process_sodp(out_SIVmod3)
SIVmod3sf = sodp_totals(SIVmod3sf)
stats_SIVmod3 = sodp_stats(SIVmod3sf, dofit=FALSE)
stats_SIVmod3.m = melt(stats_SIVmod3, id.vars=c("time", "Species", "SizeClass"))
parms_macromod3 = parms_macro
parms_macromod3$times.max = 18
parms_macromod3$times.by = 0.01
dd = optim(par = cc$par, f = function(mod_macro) {
parms_macromod3$alpha = parms_macro$alpha*mod_macro
parms_macromod3$lamda = parms_macro$lamda*mod_macro
out_macromod = run_sodp(parms_macromod3)
macromodsf = process_sodp(out_macromod)
macromodsf = sodp_totals(macromodsf)
stats_macromod = sodp_stats(macromodsf)
sitmacromod = subset(stats_macromod, SizeClass=="Total", select=c("time", "FracInf"))
sitmacromod = sitmacromod[order(sitmacromod$time),]
timeToFracImacromod = timeToFracI_fn(sitmacromod, frac_I)
if(interactive()) print(timeToFracImacromod)
val = ((timeToFracImacromod - timeToFracI)/timeToFracI)^2
return(val)
}, method="Brent", lower=0, upper=2)
mod_macro = dd$par
parms_macromod3$alpha = parms_macro$alpha*mod_macro
parms_macromod3$lamda = parms_macro$lamda*mod_macro
parms_macromod3$times.max = 50
parms_macromod3$times.by = 0.1
out_macromod3 = run_sodp(parms_macromod3)
macromod3sf = process_sodp(out_macromod3)
macromod3sf = sodp_totals(macromod3sf)
stats_macromod3 = sodp_stats(macromod3sf, dofit=FALSE)
stats_macromod3.m = melt(stats_macromod3, id.vars=c("time", "Species", "SizeClass"))
stats = list(SI=stats_SI.m, SIV=stats_SIVmod3.m, macro=stats_macromod3.m)
statsall = ldply(stats, function(x) x)
mortInf = ldply(stats, function(x) subset(x, variable %in% c("N")))
p1 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 1) +
ylab("Host Population") +
xlab(" ") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
mortInf = ldply(stats, function(x) subset(x, variable %in% c("I")))
p2 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 1) +
ylab("Infected Hosts") +
xlab("Time") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
mortInf = ldply(stats, function(x) subset(x, variable %in% c("MortInfRate")))
p3 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 0.4) +
ylab("Infected Mortality Rate") +
xlab(" ") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
g_legend <- function(p){
tmp <- ggplotGrob(p)
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)}
legend <- g_legend(p1)
lwidth <- sum(legend$width)
plotwidth = (unit(1, "npc") - lwidth)*(1/3)
grid.arrange(p1 + theme(legend.position="none"),
p2 + theme(legend.position="none"),
p3 + theme(legend.position="none"),
legend,
widths=unit.c(plotwidth, plotwidth, plotwidth, lwidth),
main = textGrob("Epidemic Dynamics: Equivalent Time to 10% infection", gp=gpar(fontsize=20,fontfamily='Lato')),
nrow=1)
frac_I = 0.1
frac_I2 = 0.3
sit = subset(stats_SI, SizeClass=="Total", select=c("time", "FracInf"))
sit = sit[order(sit$time),]
timeToFracI = timeToFracI_fn(sit, frac_I)
timeToFracI[2] = timeToFracI_fn(sit, frac_I2)
parms_SIVmod4 = parms_SIV
parms_SIVmod4$progress = 0
parms_SIVmod4$times.max = 25
parms_SIVmod4$times.by = 0.01
ee = optim(par = c(1, 1), fn = function(mod_SIV) {
parms_SIVmod4$alpha = parms_SIV$alpha*mod_SIV[1]
parms_SIVmod4$lamda = parms_SIV$lamda*mod_SIV[2]
out_SIVmod = run_sodp(parms_SIVmod4)
SIVmodsf = process_sodp(out_SIVmod)
SIVmodsf = sodp_totals(SIVmodsf)
stats_SIVmod = sodp_stats(SIVmodsf)
sitSIVmod = subset(stats_SIVmod, SizeClass=="Total", select=c("time", "FracInf"))
sitSIVmod = sitSIVmod[order(sitSIVmod$time),]
timeToFracISIVmod = timeToFracI_fn(sitSIVmod, frac_I)
timeToFracISIVmod[2] = timeToFracI_fn(sitSIVmod, frac_I2)
if(interactive()) print(timeToFracISIVmod)
val = sum(((timeToFracISIVmod - timeToFracI)/timeToFracI)^2)
return(val)
})
mod_SIV = ee$par
parms_SIVmod4$alpha = parms_SIV$alpha*mod_SIV[1]
parms_SIVmod4$lamda = parms_SIV$lamda*mod_SIV[2]
parms_SIVmod4$times.max = 50
parms_SIVmod4$times.by = 0.1
out_SIVmod4 = run_sodp(parms_SIVmod4)
SIVmod4sf = process_sodp(out_SIVmod4)
SIVmod4sf = sodp_totals(SIVmod4sf)
stats_SIVmod4 = sodp_stats(SIVmod4sf, dofit=FALSE)
stats_SIVmod4.m = melt(stats_SIVmod4, id.vars=c("time", "Species", "SizeClass"))
parms_macromod4 = parms_macro
parms_macromod4$progress = 0
parms_macromod4$times.max = 25
parms_macromod4$times.by = 0.01
ff = optim(par = ee$par, f = function(mod_macro) {
parms_macromod4$alpha = parms_macro$alpha*mod_macro[1]
parms_macromod4$lamda = parms_macro$lamda*mod_macro[2]
out_macromod = run_sodp(parms_macromod4)
macromodsf = process_sodp(out_macromod)
macromodsf = sodp_totals(macromodsf)
stats_macromod = sodp_stats(macromodsf)
sitmacromod = subset(stats_macromod, SizeClass=="Total", select=c("time", "FracInf"))
sitmacromod = sitmacromod[order(sitmacromod$time),]
timeToFracImacromod = timeToFracI_fn(sitmacromod, frac_I)
timeToFracImacromod[2] = timeToFracI_fn(sitmacromod, frac_I2)
if(interactive()) print(timeToFracImacromod)
val = sum(((timeToFracImacromod - timeToFracI)/timeToFracI)^2)
return(val)
})
mod_macro = ff$par
parms_macromod4$alpha = parms_macro$alpha*mod_macro[1]
parms_macromod4$lamda = parms_macro$lamda*mod_macro[2]
parms_macromod4$times.max = 50
parms_macromod4$times.by = 0.1
out_macromod4 = run_sodp(parms_macromod4)
macromod4sf = process_sodp(out_macromod4)
macromod4sf = sodp_totals(macromod4sf)
stats_macromod4 = sodp_stats(macromod4sf, dofit=FALSE)
stats_macromod4.m = melt(stats_macromod4, id.vars=c("time", "Species", "SizeClass"))
stats = list(SI=stats_SI.m, SIV=stats_SIVmod4.m, macro=stats_macromod4.m)
statsall = ldply(stats, function(x) x)
mortInf = ldply(stats, function(x) subset(x, variable %in% c("N")))
p1 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 1) +
ylab("Host Population") +
xlab(" ") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
mortInf = ldply(stats, function(x) subset(x, variable %in% c("I")))
p2 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 1) +
ylab("Infected Hosts") +
xlab("Time") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
mortInf = ldply(stats, function(x) subset(x, variable %in% c("MortInfRate")))
p3 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 0.6) +
ylab("Infected Mortality Rate") +
xlab(" ") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
g_legend <- function(p){
tmp <- ggplotGrob(p)
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)}
legend <- g_legend(p1)
lwidth <- sum(legend$width)
plotwidth = (unit(1, "npc") - lwidth)*(1/3)
grid.arrange(p1 + theme(legend.position="none"),
p2 + theme(legend.position="none"),
p3 + theme(legend.position="none"),
legend,
widths=unit.c(plotwidth, plotwidth, plotwidth, lwidth),
main = textGrob("Epidemic Dynamics: Equivalent time to 10% and 30% infection", gp=gpar(fontsize=20,fontfamily='Lato')),
nrow=1)
frac_I = 0.1
sit = subset(stats_SI, SizeClass=="Total", select=c("time", "FracInf"))
sit = sit[order(sit$time),]
timeToFracI = timeToFracI_fn(sit, frac_I)
timeToFracI[2] = subset(stats_SI, time==parms_SI$times.max & SizeClass=="Total")$FracInf
parms_SIVmod5 = parms_SIV
parms_SIVmod5$progress = 0
parms_SIVmod5$times.by = 0.01
gg = optim(par = c(1, 1), fn = function(mod_SIV) {
parms_SIVmod5$alpha = parms_SIV$alpha*mod_SIV[1]
parms_SIVmod5$lamda = parms_SIV$lamda*mod_SIV[2]
out_SIVmod = run_sodp(parms_SIVmod5)
SIVmodsf = process_sodp(out_SIVmod)
SIVmodsf = sodp_totals(SIVmodsf)
stats_SIVmod = sodp_stats(SIVmodsf)
sitSIVmod = subset(stats_SIVmod, SizeClass=="Total", select=c("time", "FracInf"))
sitSIVmod = sitSIVmod[order(sitSIVmod$time),]
timeToFracISIVmod = timeToFracI_fn(sitSIVmod, frac_I)
timeToFracISIVmod[2]= subset(stats_SIVmod, time==parms_SIVmod5$times.max & SizeClass=="Total")$FracInf
if(interactive()) print(timeToFracISIVmod)
val = sum(((timeToFracISIVmod - timeToFracI)/timeToFracI)^2)
return(val)
})
mod_SIV = gg$par
parms_SIVmod5$alpha = parms_SIV$alpha*mod_SIV[1]
parms_SIVmod5$lamda = parms_SIV$lamda*mod_SIV[2]
parms_SIVmod5$times.by = 0.1
out_SIVmod5 = run_sodp(parms_SIVmod5)
SIVmod5sf = process_sodp(out_SIVmod5)
SIVmod5sf = sodp_totals(SIVmod5sf)
stats_SIVmod5 = sodp_stats(SIVmod5sf, dofit=FALSE)
stats_SIVmod5.m = melt(stats_SIVmod5, id.vars=c("time", "Species", "SizeClass"))
parms_macromod5 = parms_macro
parms_macromod5$progress = 0
parms_macromod5$times.by = 0.01
hh = optim(par = gg$par, f = function(mod_macro) {
parms_macromod5$alpha = parms_macro$alpha*mod_macro[1]
parms_macromod5$lamda = parms_macro$lamda*mod_macro[2]
out_macromod = run_sodp(parms_macromod5)
macromodsf = process_sodp(out_macromod)
macromodsf = sodp_totals(macromodsf)
stats_macromod = sodp_stats(macromodsf)
sitmacromod = subset(stats_macromod, SizeClass=="Total", select=c("time", "FracInf"))
sitmacromod = sitmacromod[order(sitmacromod$time),]
timeToFracImacromod = timeToFracI_fn(sitmacromod, frac_I)
timeToFracImacromod[2] = subset(stats_macromod, time==parms_macromod5$times.max & SizeClass=="Total")$FracInf
if(interactive()) print(timeToFracImacromod)
val = sum(((timeToFracImacromod - timeToFracI)/timeToFracI)^2)
return(val)
})
mod_macro = hh$par
parms_macromod5$alpha = parms_macro$alpha*mod_macro[1]
parms_macromod5$lamda = parms_macro$lamda*mod_macro[2]
parms_macromod5$times.max = 50
parms_macromod5$times.by = 0.1
out_macromod5 = run_sodp(parms_macromod5)
macromod5sf = process_sodp(out_macromod5)
macromod5sf = sodp_totals(macromod5sf)
stats_macromod5 = sodp_stats(macromod5sf, dofit=FALSE)
stats_macromod5.m = melt(stats_macromod5, id.vars=c("time", "Species", "SizeClass"))
stats = list(SI=stats_SI.m, SIV=stats_SIVmod5.m, macro=stats_macromod5.m)
statsall = ldply(stats, function(x) x)
mortInf = ldply(stats, function(x) subset(x, variable %in% c("N")))
p1 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 1) +
ylab("Host Population") +
xlab(" ") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
mortInf = ldply(stats, function(x) subset(x, variable %in% c("I")))
p2 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 1) +
ylab("Infected Hosts") +
xlab("Time") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
mortInf = ldply(stats, function(x) subset(x, variable %in% c("MortInfRate")))
p3 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) +
geom_line() +
ylim(0, 0.6) +
ylab("Infected Mortality Rate") +
xlab(" ") +
scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) +
scale_color_grey() +
scale_size_manual(values=c(1,1, 2), guide='none') +
theme_nr +
theme(text=element_text(family="Lato")) +
guides(linetype = guide_legend(reverse=TRUE))
g_legend <- function(p){
tmp <- ggplotGrob(p)
leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend <- tmp$grobs[[leg]]
return(legend)}
legend <- g_legend(p1)
lwidth <- sum(legend$width)
plotwidth = (unit(1, "npc") - lwidth)*(1/3)
grid.arrange(p1 + theme(legend.position="none"),
p2 + theme(legend.position="none"),
p3 + theme(legend.position="none"),
legend,
widths=unit.c(plotwidth, plotwidth, plotwidth, lwidth),
main = textGrob("Epidemic Dynamics: Equivalent time to 10% infection and equilibrium mortality", gp=gpar(fontsize=20,fontfamily='Lato')),
nrow=1)
Currently matching model by equilibrium mortality. Next, match the models by
initial growth behavior. That is, growth rate at time zero, or R0. How do
models behave comparatively if they have initial growth rate of infected
individuals that is the same? What if we peg both initial rate and final
equilibrium mortality rate, but free both $\alpha$ and $\lambda$ to vary
independently? How will models behave differently?
Need too calculate this rate for the given initial number of infected individuals
Model may need to be modified to only have initial infected individuals be of
lowest disease class, not all disease classes.
At initiation, when all individuals are either susceptible or only have one
infection, the growth in infected individuals is
$$\lambda_{I_1} - d - \alpha$$
Note that the 2nd derivative of $V$ is very small no matter when $I$ is small at $V=H=0$:
$$\frac{d^2V}{dt^2} = - I^{3} \lambda^{2} + I^{2} \left(2 S \lambda^{2} - 4 \alpha \lambda - 3 d \lambda\right)$$
noamross/age-infects documentation built on May 23, 2019, 9:30 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
% Comparative Dynamics % Noam Ross % 14-08-25 10:10:34
r library(knitr); opts_chunk$set(echo=FALSE, cache=TRUE, message=FALSE, fig.width=12, fig.height=6, dpi=72)
Here I examine the comparative dynamics of age classes in 2-age-class models of both $SI$, $SIV$, and macroparasite-style dynamic models.
The following parameters are identical between models:
- fecundity $(f)$
- transition rate $(g)$
- carrying capacity $(K)$
The models are also parameterized to have identical apparent mortality rates at equilibrium. That is the aggregate mortality rates of all infection classes are the same.
parms_SI = list( n_stages = 2, n_parasites = 1, f = c(1, 1), g = c(0.1, 0), d = c(0.01, 0.01), alpha = c(0.2, 0.2), lamda = c(0.5, 0.5), Beta = c(1, 1), mu = c(0, 0), xi = c(1, 1), omega = c(1,1), K = 1, progress = 0, times.min = 0, times.max = 50, times.by = 0.1, inits = 0, infect_vector = c(0.01, 0.01) ) parms_SIV = parms_SI parms_SIV$n_parasites = 2 parms_macro = parms_SI parms_macro$n_parasites = 10 parms_SI$times.by = 0.01
library(deSolve) library(plyr) library(rootSolve) library(reshape2) library(ggplot2) library(data.table) library(noamtools) library(numDeriv) library(gridExtra) source('../R/process_sodp.R') source('../R/run_sodp.R')
out_SI = run_sodp(parms_SI) SIsf = process_sodp(out_SI) SIsf = sodp_totals(SIsf) stats_SI = sodp_stats(SIsf) stats_SI.m = melt(stats_SI, id.vars=c("time", "Species", "SizeClass")) out_SIV = run_sodp(parms_SIV) SIVsf = process_sodp(out_SIV) SIVsf = sodp_totals(SIVsf) stats_SIV = sodp_stats(SIVsf) stats_SIV.m = melt(stats_SIV, id.vars=c("time", "Species", "SizeClass")) out_macro = run_sodp(parms_macro) macrosf = process_sodp(out_macro) macrosf = sodp_totals(macrosf) stats_macro = sodp_stats(macrosf) stats_macro.m = melt(stats_macro, id.vars=c("time", "Species", "SizeClass"))
times = c(0, 5, 15, 20, 25, 30, 50) df = data.frame(subset(macrosf, time %in% times)) df = subset(df, SizeClass != "Total" & Population > 1e-5) df = df[order(df$time),] df$time = factor(paste("Time =", df$time), levels = paste("Time =", unique(df$time))) ggplot(df, aes(x=Infected, y=Population, fill=SizeClass)) + geom_area(position="identity", alpha = 0.5) + scale_fill_grey(labels = c("Juveniles", "Adults")) + theme_nr + xlab("Number of Infections") + theme(text=element_text(family="Lato")) + facet_wrap(~time, nrow=1) df2 = data.frame(subset(SIVsf, time %in% times)) df2 = droplevels(subset(df2, SizeClass != "Total" & Population > 1e-5)) df2 = df2[order(df2$time),] df2$time = factor(paste("Time =", df2$time), levels = paste("Time =", unique(df2$time))) ggplot(df2, aes(x=as.factor(Infected), y=Population, fill=SizeClass)) + geom_bar(stat="identity", position="dodge", alpha = 0.5) + scale_fill_grey(labels = c("Juveniles", "Adults")) + scale_x_discrete(labels = c("S", "I", "V")) + theme_nr + xlab("Disease Class") + theme(text=element_text(family="Lato")) + facet_wrap(~time, nrow=1)
stats = list(SI=stats_SI.m, SIV=stats_SIV.m, macro=stats_macro.m) statsall = ldply(stats, function(x) x) mortInf = ldply(stats, function(x) subset(x, variable %in% c("N"))) p1 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 1) + ylab("Host Population") + xlab(" ") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) mortInf = ldply(stats, function(x) subset(x, variable %in% c("I"))) p2 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 1) + ylab("Infected Hosts") + xlab("Time") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) mortInf = ldply(stats, function(x) subset(x, variable %in% c("MortInfRate"))) p3 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 0.4) + ylab("Infected Mortality Rate") + xlab(" ") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) g_legend <- function(p){ tmp <- ggplotGrob(p) leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box") legend <- tmp$grobs[[leg]] return(legend)} legend <- g_legend(p1) lwidth <- sum(legend$width) plotwidth = (unit(1, "npc") - lwidth)*(1/3) grid.arrange(p1 + theme(legend.position="none"), p2 + theme(legend.position="none"), p3 + theme(legend.position="none"), legend, widths=unit.c(plotwidth, plotwidth, plotwidth, lwidth), main = textGrob("Epidemic Dynamics: 1-infection equivalence", gp=gpar(fontsize=20,fontfamily='Lato')), nrow=1)
parms_SIVmod = parms_SIV aa = optim(par = 1, f = function(mod_SIV) { parms_SIVmod$alpha = parms_SIV$alpha*mod_SIV parms_SIVmod$lamda = parms_SIV$lamda*mod_SIV out_SIVmod = run_sodp(parms_SIVmod) SIVmodsf = process_sodp(out_SIVmod) SIVmodsf = sodp_totals(SIVmodsf) stats_SIVmod = sodp_stats(SIVmodsf) val = subset(stats_SIVmod, time==parms_SIVmod$times.max & SizeClass=="Total")$MortInfRate if(interactive()) print(val) return((val - parms_SI$alpha[1])^2) }, method="Brent", lower=0, upper=2) mod_SIV = aa$par parms_SIVmod$alpha = parms_SIV$alpha*mod_SIV parms_SIVmod$lamda = parms_SIV$lamda*mod_SIV out_SIVmod = run_sodp(parms_SIVmod) SIVmodsf = process_sodp(out_SIVmod) SIVmodsf = sodp_totals(SIVmodsf) stats_SIVmod = sodp_stats(SIVmodsf) stats_SIVmod.m = melt(stats_SIVmod, id.vars=c("time", "Species", "SizeClass")) parms_macromod = parms_macro bb = optim(par = 0.6118, f = function(mod_macro) { parms_macromod$alpha = parms_macro$alpha*mod_macro parms_macromod$lamda = parms_macro$lamda*mod_macro out_macromod = run_sodp(parms_macromod) macromodsf = process_sodp(out_macromod) macromodsf = sodp_totals(macromodsf) stats_macromod = sodp_stats(macromodsf) val = subset(stats_macromod, time==parms_macromod$times.max & SizeClass=="Total")$MortInfRate if(interactive()) print(val) return((val - parms_SI$alpha[1])^2) }, method="Brent", lower=0, upper=2) mod_macro = bb$par parms_macromod$alpha = parms_macro$alpha*mod_macro parms_macromod$lamda = parms_macro$lamda*mod_macro out_macromod = run_sodp(parms_macromod) macromodsf = process_sodp(out_macromod) macromodsf = sodp_totals(macromodsf) stats_macromod = sodp_stats(macromodsf) stats_macromod.m = melt(stats_macromod, id.vars=c("time", "Species", "SizeClass"))
stats = list(SI=stats_SI.m, SIV=stats_SIVmod.m, macro=stats_macromod.m) statsall = ldply(stats, function(x) x) mortInf = ldply(stats, function(x) subset(x, variable %in% c("N"))) p1 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 1) + ylab("Host Population") + xlab(" ") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) mortInf = ldply(stats, function(x) subset(x, variable %in% c("I"))) p2 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 1) + ylab("Infected Hosts") + xlab("Time") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) mortInf = ldply(stats, function(x) subset(x, variable %in% c("MortInfRate"))) p3 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 0.4) + ylab("Infected Mortality Rate") + xlab(" ") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) g_legend <- function(p){ tmp <- ggplotGrob(p) leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box") legend <- tmp$grobs[[leg]] return(legend)} legend <- g_legend(p1) lwidth <- sum(legend$width) plotwidth = (unit(1, "npc") - lwidth)*(1/3) grid.arrange(p1 + theme(legend.position="none"), p2 + theme(legend.position="none"), p3 + theme(legend.position="none"), legend, widths=unit.c(plotwidth, plotwidth, plotwidth, lwidth), main = textGrob("Epidemic Dynamics: Equilibrium Mortality Equivalence", gp=gpar(fontsize=20,fontfamily='Lato')), nrow=1)
SI_Final = init_derivs2(parms_SI, SI=TRUE)[,2] parms_SIVmod2 = parms_SIV aa = optim(par = c(1,1), fn = function(mod_SIV) { parms_SIVmod$alpha = parms_SIV$alpha*mod_SIV[1] parms_SIVmod$lamda = parms_SIV$lamda*mod_SIV[2] val = init_derivs2(parms_SIVmod, SI=TRUE)[,2] if(interactive()) print(val) return(sum(((val - SI_Final)/SI_Final)^2)) }) mod_SIV = aa$par parms_SIVmod2$alpha = parms_SIV$alpha*mod_SIV[1] parms_SIVmod2$lamda = parms_SIV$lamda*mod_SIV[2] out_SIVmod2 = run_sodp(parms_SIVmod2) SIVmod2sf = process_sodp(out_SIVmod2) SIVmod2sf = sodp_totals(SIVmod2sf) stats_SIVmod2 = sodp_stats(SIVmod2sf) stats_SIVmod2.m = melt(stats_SIVmod2, id.vars=c("time", "Species", "SizeClass")) parms_macromod2 = parms_macro bb = optim(par = aa$par, f = function(mod_macro) { parms_macromod$alpha = parms_macro$alpha*mod_macro[1] parms_macromod$lamda = parms_macro$lamda*mod_macro[2] val = init_derivs2(parms_macromod, SI=TRUE)[,2] if(interactive()) print(val) return(sum(((val - SI_Final)/SI_Final)^2)) }) mod_macro = bb$par parms_macromod2$alpha = parms_macro$alpha*mod_macro[1] parms_macromod2$lamda = parms_macro$lamda*mod_macro[2] out_macromod2 = run_sodp(parms_macromod2) macromod2sf = process_sodp(out_macromod2) macromod2sf = sodp_totals(macromod2sf) stats_macromod2 = sodp_stats(macromod2sf) stats_macromod2.m = melt(stats_macromod2, id.vars=c("time", "Species", "SizeClass"))
stats = list(SI=stats_SI.m, SIV=stats_SIVmod2.m, macro=stats_macromod2.m) statsall = ldply(stats, function(x) x) mortInf = ldply(stats, function(x) subset(x, variable %in% c("N"))) p1 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 1) + ylab("Host Population") + xlab(" ") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) mortInf = ldply(stats, function(x) subset(x, variable %in% c("I"))) p2 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 1) + ylab("Infected Hosts") + xlab("Time") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) mortInf = ldply(stats, function(x) subset(x, variable %in% c("MortInfRate"))) p3 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 0.6) + ylab("Infected Mortality Rate") + xlab(" ") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) g_legend <- function(p){ tmp <- ggplotGrob(p) leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box") legend <- tmp$grobs[[leg]] return(legend)} legend <- g_legend(p1) lwidth <- sum(legend$width) plotwidth = (unit(1, "npc") - lwidth)*(1/3) grid.arrange(p1 + theme(legend.position="none"), p2 + theme(legend.position="none"), p3 + theme(legend.position="none"), legend, widths=unit.c(plotwidth, plotwidth, plotwidth, lwidth), main = textGrob("Epidemic Dynamics: Equivalent Initial Derivatives", gp=gpar(fontsize=20,fontfamily='Lato')), nrow=1)
frac_I = 0.1 sit = subset(stats_SI, SizeClass=="Total", select=c("time", "FracInf")) sit = sit[order(sit$time),] timeToFracI_fn = function(sit, frac_I) { afti = min(which(sit$FracInf >= frac_I), na.rm=TRUE) if(is.na(afti) || afti==Inf) return(max(sit$time)) befi = min(which(sit$FracInf >= frac_I)) - 1 befx = sit$time[befi] befy = sit$FracInf[befi] aftx = sit$time[afti] afty = sit$FracInf[afti] slope = (afty-befy)/(aftx-befx) return(befx + (frac_I - befy)/slope) } timeToFracI = timeToFracI_fn(sit, frac_I) parms_SIVmod3 = parms_SIV parms_SIVmod3$times.max = 18 parms_SIVmod3$times.by = 0.01 cc = optim(par = 1, fn = function(mod_SIV) { parms_SIVmod3$alpha = parms_SIV$alpha*mod_SIV parms_SIVmod3$lamda = parms_SIV$lamda*mod_SIV out_SIVmod = run_sodp(parms_SIVmod3) SIVmodsf = process_sodp(out_SIVmod) SIVmodsf = sodp_totals(SIVmodsf) stats_SIVmod = sodp_stats(SIVmodsf) sitSIVmod = subset(stats_SIVmod, SizeClass=="Total", select=c("time", "FracInf")) sitSIVmod = sitSIVmod[order(sitSIVmod$time),] timeToFracISIVmod = timeToFracI_fn(sitSIVmod, frac_I) if(interactive()) print(timeToFracISIVmod) val = ((timeToFracISIVmod - timeToFracI)/timeToFracI)^2 return(val) }, method="Brent", lower=0, upper=2) mod_SIV = cc$par parms_SIVmod3$alpha = parms_SIV$alpha*mod_SIV parms_SIVmod3$lamda = parms_SIV$lamda*mod_SIV parms_SIVmod3$times.max = 50 parms_SIVmod3$times.by = 0.1 out_SIVmod3 = run_sodp(parms_SIVmod3) SIVmod3sf = process_sodp(out_SIVmod3) SIVmod3sf = sodp_totals(SIVmod3sf) stats_SIVmod3 = sodp_stats(SIVmod3sf, dofit=FALSE) stats_SIVmod3.m = melt(stats_SIVmod3, id.vars=c("time", "Species", "SizeClass")) parms_macromod3 = parms_macro parms_macromod3$times.max = 18 parms_macromod3$times.by = 0.01 dd = optim(par = cc$par, f = function(mod_macro) { parms_macromod3$alpha = parms_macro$alpha*mod_macro parms_macromod3$lamda = parms_macro$lamda*mod_macro out_macromod = run_sodp(parms_macromod3) macromodsf = process_sodp(out_macromod) macromodsf = sodp_totals(macromodsf) stats_macromod = sodp_stats(macromodsf) sitmacromod = subset(stats_macromod, SizeClass=="Total", select=c("time", "FracInf")) sitmacromod = sitmacromod[order(sitmacromod$time),] timeToFracImacromod = timeToFracI_fn(sitmacromod, frac_I) if(interactive()) print(timeToFracImacromod) val = ((timeToFracImacromod - timeToFracI)/timeToFracI)^2 return(val) }, method="Brent", lower=0, upper=2) mod_macro = dd$par parms_macromod3$alpha = parms_macro$alpha*mod_macro parms_macromod3$lamda = parms_macro$lamda*mod_macro parms_macromod3$times.max = 50 parms_macromod3$times.by = 0.1 out_macromod3 = run_sodp(parms_macromod3) macromod3sf = process_sodp(out_macromod3) macromod3sf = sodp_totals(macromod3sf) stats_macromod3 = sodp_stats(macromod3sf, dofit=FALSE) stats_macromod3.m = melt(stats_macromod3, id.vars=c("time", "Species", "SizeClass"))
stats = list(SI=stats_SI.m, SIV=stats_SIVmod3.m, macro=stats_macromod3.m) statsall = ldply(stats, function(x) x) mortInf = ldply(stats, function(x) subset(x, variable %in% c("N"))) p1 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 1) + ylab("Host Population") + xlab(" ") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) mortInf = ldply(stats, function(x) subset(x, variable %in% c("I"))) p2 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 1) + ylab("Infected Hosts") + xlab("Time") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) mortInf = ldply(stats, function(x) subset(x, variable %in% c("MortInfRate"))) p3 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 0.4) + ylab("Infected Mortality Rate") + xlab(" ") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) g_legend <- function(p){ tmp <- ggplotGrob(p) leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box") legend <- tmp$grobs[[leg]] return(legend)} legend <- g_legend(p1) lwidth <- sum(legend$width) plotwidth = (unit(1, "npc") - lwidth)*(1/3) grid.arrange(p1 + theme(legend.position="none"), p2 + theme(legend.position="none"), p3 + theme(legend.position="none"), legend, widths=unit.c(plotwidth, plotwidth, plotwidth, lwidth), main = textGrob("Epidemic Dynamics: Equivalent Time to 10% infection", gp=gpar(fontsize=20,fontfamily='Lato')), nrow=1)
frac_I = 0.1 frac_I2 = 0.3 sit = subset(stats_SI, SizeClass=="Total", select=c("time", "FracInf")) sit = sit[order(sit$time),] timeToFracI = timeToFracI_fn(sit, frac_I) timeToFracI[2] = timeToFracI_fn(sit, frac_I2) parms_SIVmod4 = parms_SIV parms_SIVmod4$progress = 0 parms_SIVmod4$times.max = 25 parms_SIVmod4$times.by = 0.01 ee = optim(par = c(1, 1), fn = function(mod_SIV) { parms_SIVmod4$alpha = parms_SIV$alpha*mod_SIV[1] parms_SIVmod4$lamda = parms_SIV$lamda*mod_SIV[2] out_SIVmod = run_sodp(parms_SIVmod4) SIVmodsf = process_sodp(out_SIVmod) SIVmodsf = sodp_totals(SIVmodsf) stats_SIVmod = sodp_stats(SIVmodsf) sitSIVmod = subset(stats_SIVmod, SizeClass=="Total", select=c("time", "FracInf")) sitSIVmod = sitSIVmod[order(sitSIVmod$time),] timeToFracISIVmod = timeToFracI_fn(sitSIVmod, frac_I) timeToFracISIVmod[2] = timeToFracI_fn(sitSIVmod, frac_I2) if(interactive()) print(timeToFracISIVmod) val = sum(((timeToFracISIVmod - timeToFracI)/timeToFracI)^2) return(val) }) mod_SIV = ee$par parms_SIVmod4$alpha = parms_SIV$alpha*mod_SIV[1] parms_SIVmod4$lamda = parms_SIV$lamda*mod_SIV[2] parms_SIVmod4$times.max = 50 parms_SIVmod4$times.by = 0.1 out_SIVmod4 = run_sodp(parms_SIVmod4) SIVmod4sf = process_sodp(out_SIVmod4) SIVmod4sf = sodp_totals(SIVmod4sf) stats_SIVmod4 = sodp_stats(SIVmod4sf, dofit=FALSE) stats_SIVmod4.m = melt(stats_SIVmod4, id.vars=c("time", "Species", "SizeClass")) parms_macromod4 = parms_macro parms_macromod4$progress = 0 parms_macromod4$times.max = 25 parms_macromod4$times.by = 0.01 ff = optim(par = ee$par, f = function(mod_macro) { parms_macromod4$alpha = parms_macro$alpha*mod_macro[1] parms_macromod4$lamda = parms_macro$lamda*mod_macro[2] out_macromod = run_sodp(parms_macromod4) macromodsf = process_sodp(out_macromod) macromodsf = sodp_totals(macromodsf) stats_macromod = sodp_stats(macromodsf) sitmacromod = subset(stats_macromod, SizeClass=="Total", select=c("time", "FracInf")) sitmacromod = sitmacromod[order(sitmacromod$time),] timeToFracImacromod = timeToFracI_fn(sitmacromod, frac_I) timeToFracImacromod[2] = timeToFracI_fn(sitmacromod, frac_I2) if(interactive()) print(timeToFracImacromod) val = sum(((timeToFracImacromod - timeToFracI)/timeToFracI)^2) return(val) }) mod_macro = ff$par parms_macromod4$alpha = parms_macro$alpha*mod_macro[1] parms_macromod4$lamda = parms_macro$lamda*mod_macro[2] parms_macromod4$times.max = 50 parms_macromod4$times.by = 0.1 out_macromod4 = run_sodp(parms_macromod4) macromod4sf = process_sodp(out_macromod4) macromod4sf = sodp_totals(macromod4sf) stats_macromod4 = sodp_stats(macromod4sf, dofit=FALSE) stats_macromod4.m = melt(stats_macromod4, id.vars=c("time", "Species", "SizeClass"))
stats = list(SI=stats_SI.m, SIV=stats_SIVmod4.m, macro=stats_macromod4.m) statsall = ldply(stats, function(x) x) mortInf = ldply(stats, function(x) subset(x, variable %in% c("N"))) p1 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 1) + ylab("Host Population") + xlab(" ") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) mortInf = ldply(stats, function(x) subset(x, variable %in% c("I"))) p2 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 1) + ylab("Infected Hosts") + xlab("Time") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) mortInf = ldply(stats, function(x) subset(x, variable %in% c("MortInfRate"))) p3 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 0.6) + ylab("Infected Mortality Rate") + xlab(" ") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) g_legend <- function(p){ tmp <- ggplotGrob(p) leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box") legend <- tmp$grobs[[leg]] return(legend)} legend <- g_legend(p1) lwidth <- sum(legend$width) plotwidth = (unit(1, "npc") - lwidth)*(1/3) grid.arrange(p1 + theme(legend.position="none"), p2 + theme(legend.position="none"), p3 + theme(legend.position="none"), legend, widths=unit.c(plotwidth, plotwidth, plotwidth, lwidth), main = textGrob("Epidemic Dynamics: Equivalent time to 10% and 30% infection", gp=gpar(fontsize=20,fontfamily='Lato')), nrow=1)
frac_I = 0.1 sit = subset(stats_SI, SizeClass=="Total", select=c("time", "FracInf")) sit = sit[order(sit$time),] timeToFracI = timeToFracI_fn(sit, frac_I) timeToFracI[2] = subset(stats_SI, time==parms_SI$times.max & SizeClass=="Total")$FracInf parms_SIVmod5 = parms_SIV parms_SIVmod5$progress = 0 parms_SIVmod5$times.by = 0.01 gg = optim(par = c(1, 1), fn = function(mod_SIV) { parms_SIVmod5$alpha = parms_SIV$alpha*mod_SIV[1] parms_SIVmod5$lamda = parms_SIV$lamda*mod_SIV[2] out_SIVmod = run_sodp(parms_SIVmod5) SIVmodsf = process_sodp(out_SIVmod) SIVmodsf = sodp_totals(SIVmodsf) stats_SIVmod = sodp_stats(SIVmodsf) sitSIVmod = subset(stats_SIVmod, SizeClass=="Total", select=c("time", "FracInf")) sitSIVmod = sitSIVmod[order(sitSIVmod$time),] timeToFracISIVmod = timeToFracI_fn(sitSIVmod, frac_I) timeToFracISIVmod[2]= subset(stats_SIVmod, time==parms_SIVmod5$times.max & SizeClass=="Total")$FracInf if(interactive()) print(timeToFracISIVmod) val = sum(((timeToFracISIVmod - timeToFracI)/timeToFracI)^2) return(val) }) mod_SIV = gg$par parms_SIVmod5$alpha = parms_SIV$alpha*mod_SIV[1] parms_SIVmod5$lamda = parms_SIV$lamda*mod_SIV[2] parms_SIVmod5$times.by = 0.1 out_SIVmod5 = run_sodp(parms_SIVmod5) SIVmod5sf = process_sodp(out_SIVmod5) SIVmod5sf = sodp_totals(SIVmod5sf) stats_SIVmod5 = sodp_stats(SIVmod5sf, dofit=FALSE) stats_SIVmod5.m = melt(stats_SIVmod5, id.vars=c("time", "Species", "SizeClass")) parms_macromod5 = parms_macro parms_macromod5$progress = 0 parms_macromod5$times.by = 0.01 hh = optim(par = gg$par, f = function(mod_macro) { parms_macromod5$alpha = parms_macro$alpha*mod_macro[1] parms_macromod5$lamda = parms_macro$lamda*mod_macro[2] out_macromod = run_sodp(parms_macromod5) macromodsf = process_sodp(out_macromod) macromodsf = sodp_totals(macromodsf) stats_macromod = sodp_stats(macromodsf) sitmacromod = subset(stats_macromod, SizeClass=="Total", select=c("time", "FracInf")) sitmacromod = sitmacromod[order(sitmacromod$time),] timeToFracImacromod = timeToFracI_fn(sitmacromod, frac_I) timeToFracImacromod[2] = subset(stats_macromod, time==parms_macromod5$times.max & SizeClass=="Total")$FracInf if(interactive()) print(timeToFracImacromod) val = sum(((timeToFracImacromod - timeToFracI)/timeToFracI)^2) return(val) }) mod_macro = hh$par parms_macromod5$alpha = parms_macro$alpha*mod_macro[1] parms_macromod5$lamda = parms_macro$lamda*mod_macro[2] parms_macromod5$times.max = 50 parms_macromod5$times.by = 0.1 out_macromod5 = run_sodp(parms_macromod5) macromod5sf = process_sodp(out_macromod5) macromod5sf = sodp_totals(macromod5sf) stats_macromod5 = sodp_stats(macromod5sf, dofit=FALSE) stats_macromod5.m = melt(stats_macromod5, id.vars=c("time", "Species", "SizeClass"))
stats = list(SI=stats_SI.m, SIV=stats_SIVmod5.m, macro=stats_macromod5.m) statsall = ldply(stats, function(x) x) mortInf = ldply(stats, function(x) subset(x, variable %in% c("N"))) p1 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 1) + ylab("Host Population") + xlab(" ") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) mortInf = ldply(stats, function(x) subset(x, variable %in% c("I"))) p2 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 1) + ylab("Infected Hosts") + xlab("Time") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) mortInf = ldply(stats, function(x) subset(x, variable %in% c("MortInfRate"))) p3 = ggplot(mortInf, aes(x = time, y=value, col=.id, linetype=SizeClass, size=SizeClass)) + geom_line() + ylim(0, 0.6) + ylab("Infected Mortality Rate") + xlab(" ") + scale_linetype_manual(labels=c("Juv", "Adult", "Total"), values=c(2,6,1)) + scale_color_grey() + scale_size_manual(values=c(1,1, 2), guide='none') + theme_nr + theme(text=element_text(family="Lato")) + guides(linetype = guide_legend(reverse=TRUE)) g_legend <- function(p){ tmp <- ggplotGrob(p) leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box") legend <- tmp$grobs[[leg]] return(legend)} legend <- g_legend(p1) lwidth <- sum(legend$width) plotwidth = (unit(1, "npc") - lwidth)*(1/3) grid.arrange(p1 + theme(legend.position="none"), p2 + theme(legend.position="none"), p3 + theme(legend.position="none"), legend, widths=unit.c(plotwidth, plotwidth, plotwidth, lwidth), main = textGrob("Epidemic Dynamics: Equivalent time to 10% infection and equilibrium mortality", gp=gpar(fontsize=20,fontfamily='Lato')), nrow=1)
Currently matching model by equilibrium mortality. Next, match the models by initial growth behavior. That is, growth rate at time zero, or R0. How do models behave comparatively if they have initial growth rate of infected individuals that is the same? What if we peg both initial rate and final equilibrium mortality rate, but free both $\alpha$ and $\lambda$ to vary independently? How will models behave differently?
Need too calculate this rate for the given initial number of infected individuals Model may need to be modified to only have initial infected individuals be of lowest disease class, not all disease classes.
At initiation, when all individuals are either susceptible or only have one infection, the growth in infected individuals is
$$\lambda_{I_1} - d - \alpha$$
Note that the 2nd derivative of $V$ is very small no matter when $I$ is small at $V=H=0$:
$$\frac{d^2V}{dt^2} = - I^{3} \lambda^{2} + I^{2} \left(2 S \lambda^{2} - 4 \alpha \lambda - 3 d \lambda\right)$$
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.