In parksw3/fitode: Tools for Ordinary Differential Equations Model Fitting
[
\newcommand{\rzero}{{\cal R}_0}
]
I'm going to try to understand why the Bombay data are so hard to fit/find reasonable confidence intervals for the parameters of an SIR model. We already know from Bacaer's paper that the closed-population, constant-transmission SIR is not a good biological model for the Bombay plague epidemic, but it is "SIR-shaped" and, for historical comparisons (i.e. to compare with Kermack and McKendrick's paper) we would like to be able to fit it with an SIR model even if that's not biologically sensible.
My thoughts on identifiability:
- It makes sense that $r = \beta-\gamma$ and $I(0)$ should be relatively easily
identifiable (from the initial epidemic curve). If we assume we're
starting from zero recovered individuals, then $S(0) \approx N$. Because $\rzero$ is monotonically related to final size (as a proportion of the
population size), so $\rzero$ and $N$ (effective population size) should be
jointly unidentifiable.
- If we reparameterized as ${r, \rzero, N, I_0}$ then we should be able to estimate ${r, \rzero, I(0)}$ if we fix $N$.
- On the other hand, fixing $\gamma$ should work in principle (once $\gamma$
is fixed the initial slope should determine beta, which will determine
$\rzero$, so we can estimate $S(0)$ ($\approx N$) independently from final size ...), but ...
- estimating dispersion (negative binomial $k$ or Gaussian residual variance) should be independent of these \ldots
- We can also reason from the fact that Kermack and McKendrick's formula works (fitting an SIR curve with only three parameters [which are not the same as, but correspond to, epidemic peak height, peak time, and half-initial growth rate]), that we only have three free parameters to work with.
I'm going to work on this by brute force (optim() + deSolve() only) to avoid the possibilities that something weird is going on in bbmle and/or fitode.
To give ourselves the best chance of having everything work, we're going to parameterize this robustly: ${ \log(r), \log(\rzero-1), \log(I(0)) }$ (thought experiment: given the choices ${\beta, \gamma, r, \rzero, G}$ (where $G$ is generation time), which are the closest to independent conditional on an observed epidemic curve? This relates back to some of JD's stuff \ldots
We'll stick with the assumption of 100\% infection fatality rate for now.
library(deSolve)
data("bombay", package = "fitode")
## incidence-based version
bombay2 <- rbind(
c(times=bombay$week[1] -
diff(bombay$week)[1], mort=NA),
bombay
)
## set up sir.grad so we can switch parameterizations later on
options(sir.grad_param = "R0m1_r")
sir.grad <- function (t, x, params) {
with(as.list(c(x, params)), {
parameterization <- getOption("sir.grad_param", "R0m1_r")
if (parameterization == "R0m1_r") {
R0m1 <- exp(logR0m1)
r <- exp(logr)
gamma <- r/R0m1
beta <- (r + gamma)/N
} else if (parameterization == "R0m1_loggamma") {
R0m1 <- exp(logR0m1)
gamma <- exp(loggamma)
beta <- (R0m1+1)*gamma/N
} else if (parameterization == "logbeta_loggamma") {
gamma <- exp(loggamma)
beta <- exp(logbeta)
}
dS.dt <- -beta*S*I
dI.dt <- beta*S*I-gamma*I
dR.dt <- gamma*I
dxdt <- c(dS.dt,dI.dt,dR.dt)
list(dxdt)
})
}
sir.pred <- function(p) {
times <- bombay2$week
I0 <- exp(p[["logI0"]])
S0 <- exp(p[["logS0"]])
ode.params <- c(xfun(p, c("logI0", "logS0", "logk")), N = S0)
xstart <- c(S=S0, I=I0, R=0)
out <- deSolve::ode(func=sir.grad,
y=xstart,
times=times,
parms=ode.params
) |> as.data.frame()
diff(out$R)
}
xfun <- function(x, n) x[setdiff(names(x), n)]
sir.nll <- function(p, N = N) {
## set S0 to N, ignore I0
sirpred <- sir.pred(c(xfun(p, "logk"), logS0 = log(N)))
-sum(dnbinom(x = bombay2$mort[-1],
mu = sirpred,
size = exp(p["logk"]), log = TRUE))
}
start0 <- c(logR0m1 = log(4), logr = log(log(2)/1.5), logI0 = log(0.01),
logk = log(10))
p0 <- sir.pred(c(xfun(start0, "logk"), logS0 = 0))
plot(p0, log = "y")
I'm going to use 820,000 for $N$ (from Wikipedia this was the value from the 1891 census, five years earlier - obviously the effective population size could be different from that but ...)
bombay_pop <- 820e3
Eyeballing to adjust ($I(0)$ needs to be larger than we would have thought, to account for the delay from infection to mortality)
start1 <- start0
start1[["logI0"]] <- log(3)
start1[["logR0m1"]] <- log(0.1)
plot(mort ~ week, data = bombay2, log = "y", ylim = c(1,10000))
p1 <- sir.pred(c(xfun(start1, "logk"), logS0 = log(bombay_pop)))
lines(p1)
This is actually pretty terrible, but should be close enough for starting (should spend more time thinking about how
sir.nll(start1, N = bombay_pop)
opt1 <- optim(start1, fn = sir.nll, N = bombay_pop, control = list(maxit = 1e6),
hessian = TRUE)
plot(mort ~ week, data = bombay, log = "y")
p2 <- sir.pred(c(xfun(opt1$par, "logk"), logS0 = log(bombay_pop)))
lines(p2)
v <- solve(opt1$hessian)
sds <- sqrt(diag(v))
raw_pars <- cbind(est = opt1$par, lwr = opt1$par-1.96*sds, upr = opt1$par+1.96*sds)
human_pars <- rbind(R0 = exp(raw_pars["logR0m1",])+1,
r = exp(raw_pars["logr",]),
I0 = exp(raw_pars["logI0",]),
k = exp(raw_pars["logk",]))
print(human_pars, digits = 2)
Parameters are correlated, but not pathologically so:
print(cov2cor(v), digits = 2)
Simulate outcomes (confidence intervals and prediction intervals by "population prediction interval", i.e. MVN simulation + nbinom observation error)
nsim <- 1000
nt <- length(p1)
confmat <- predmat <- matrix(NA_real_, nrow = nt, ncol = nsim)
set.seed(101)
parmat <- MASS::mvrnorm(nsim, mu = opt1$par, Sigma = v)
for (i in 1:nsim) {
confmat[,i] <- sir.pred(c(xfun(parmat[i,], "logk"), logS0 = log(bombay_pop)))
predmat[,i] <- rnbinom(nt, mu = confmat[,i], size = exp(opt1$par[["logk"]]))
}
pconf <- t(apply(confmat, 1, quantile, c(0.025, 0.975)))
ppred <- t(apply(predmat, 1, quantile, c(0.025, 0.975)))
par(las=1, bty = "l")
plot(mort ~ week, data = bombay, log = "y", ylim = c(3, 1500))
lines(p2)
matlines(pconf, col = 4, lty = 2)
matlines(ppred, col = 6, lty = 3)
legend("topleft", lty = c(1, 2, 3),
legend = c("best fit", "95% conf int", "95% pred int"),
col = c(1, 4, 6))
What about translating back from parameters all the way to $\beta/\gamma$?
trfun <- function(logr, logR0m1, logI0, logk, logS0, N = log(bombay_pop)) {
I0 <- exp(logI0)
if (!missing(logS0)) {
##
S0 <- exp(logS0)
N <- S0 + I0
} else {
S0 <- N - I0
}
c(beta = exp(logr)*(1+1/exp(logR0m1))/N,
gamma = exp(logr - logR0m1),
S0 = S0,
I0 = I0,
k = exp(logk))
}
trmat <- apply(parmat, 1, \(x) do.call(trfun, as.list(x)))
print(cbind(est = do.call(trfun, as.list(opt1$par)),
lwr = apply(trmat, 1, quantile, 0.025),
upr = apply(trmat, 1, quantile, 0.975)),
digits = 3)
The $\gamma$ estimate is obviously silly (95% CI approx. 70-80).
What if we fix $I(0)=1$ as done in the paper and allow $S(0)$ (effective population size) to vary instead?
sir.nll2 <- function(p, logI0 = log(1)) {
sirpred <- sir.pred(c(xfun(p, "logk"), logI0 = logI0))
-sum(dnbinom(x = bombay2$mort[-1],
mu = sirpred,
size = exp(p["logk"]), log = TRUE))
}
start2 <- c(xfun(start1, "logI0"), logS0 = log(bombay_pop))
opt2 <- optim(start2, fn = sir.nll2, logI0 = log(1),
control = list(maxit = 1e6),
hessian = TRUE)
p2B <- sir.pred(c(xfun(opt2$par, "logk"), logI0 = log(1)))
v2 <- solve(opt2$hessian)
sds2 <- sqrt(diag(v))
raw_pars2 <- cbind(est = opt2$par,
lwr = opt2$par-1.96*sds2, upr = opt2$par+1.96*sds2)
human_pars2 <- rbind(R0 = exp(raw_pars2["logR0m1",])+1,
r = exp(raw_pars2["logr",]),
S0 = exp(raw_pars2["logS0",]),
k = exp(raw_pars2["logk",]))
print(human_pars2, digits = 3)
Stronger correlations (although still not insane)
print(cov2cor(v2), digits = 2)
nsim <- 1000
nt <- length(p1)
confmat2 <- predmat2 <- matrix(NA_real_, nrow = nt, ncol = nsim)
set.seed(101)
parmat2 <- MASS::mvrnorm(nsim, mu = opt2$par, Sigma = v)
for (i in 1:nsim) {
confmat2[,i] <- sir.pred(c(xfun(parmat2[i,], "logk"), logI0 = log(1)))
predmat2[,i] <- rnbinom(nt, mu = confmat2[,i], size = exp(opt2$par[["logk"]]))
}
## DRY ...
pconf2 <- t(apply(confmat2, 1, quantile, c(0.025, 0.975)))
ppred2 <- t(apply(predmat2, 1, quantile, c(0.025, 0.975)))
par(las=1, bty = "l")
plot(mort ~ week, data = bombay, log = "y", ylim = c(3, 1500))
lines(p2B)
lines(p2, col = 2)
matlines(pconf2, col = 4, lty = 2)
matlines(ppred2, col = 6, lty = 3)
legend("bottom", lty = c(1, 2, 3),
legend = c("best fit (fix I0/gamma)",
"best fit (fix N/?)",
"95% conf int", "95% pred int"),
col = c(1, 2, 4, 6))
trfun2 <- function(logr, logR0m1, logS0, logk) {
c(beta = exp(logr)*(1+1/exp(logR0m1))/exp(logS0),
gamma = exp(logr - logR0m1),
S0 = exp(logS0),
I0 = 1,
k = exp(logk))
}
trmat2 <- apply(parmat2, 1, \(x) do.call(trfun2, as.list(x)))
print(cbind(est = do.call(trfun2, as.list(opt2$par)),
lwr = apply(trmat2, 1, quantile, 0.025),
upr = apply(trmat2, 1, quantile, 0.975)),
digits = 2)
Correlations among beta/gamma/S0/k aren't too bad ...
cor(t(trmat2[-4,]))
The value of $\gamma$ can compensate for a fixed value of $I(0)$ because we are measuring incidence (i.e. the initial observed mortality rate is approximately $\gamma I(0)$).
Do we need to worry about incidence vs prevalence when translating K-M parameters to biological parameters ... ?
Gamma vs. N
options(sir.grad_param = "R0m1_loggamma")
start3 <- c(xfun(start1, c("logI0", "logr")),
loggamma = log(4.77), logS0 = log(bombay_pop))
sir.nll2(start3, logI0 = log(1))
opt3 <- optim(start3, fn = sir.nll2, logI0 = log(1),
control = list(maxit = 1e6),
hessian = TRUE)
exp(opt3$par)
## NLL with fixed value of gamma
sir.nll3 <- function(p, logI0 = log(1), loggamma = log(4.77)) {
sirpred <- sir.pred(c(xfun(p, "logk"), logI0 = logI0, loggamma = loggamma))
-sum(dnbinom(x = bombay2$mort[-1],
mu = sirpred,
size = exp(p["logk"]), log = TRUE))
}
opt4 <- optim(xfun(start3, "loggamma"),
fn = sir.nll3,
logI0 = log(1), loggamma = log(4.77),
control = list(maxit = 1e6),
hessian = FALSE)
lgammavec <- seq(log(1), log(100), length = 51)
profmat <- matrix(NA_real_, nrow = length(lgammavec), ncol = 4,
dimnames = list(NULL, c(names(opt4$par), "nll")))
for (i in seq_along(lgammavec)) {
pstart <- start3
## cat(i, lgammavec[i], "\n")
oo <- try(optim(
xfun(start3, "loggamma"),
fn = sir.nll3,
logI0 = log(1), loggamma = lgammavec[i],
control = list(maxit = 1e6),
hessian = FALSE),
silent = TRUE
)
if (!is(oo, "try-error")) {
pstart[names(oo$par)] <- oo$par
profmat[i,] <- c(oo$par, oo$value)
}
}
par(mfrow = c(2,2), las = 1, bty = "l", mar = c(3,4,1,1),
mgp = c(2, 1, 0))
for (i in 1:4) {
plot(lgammavec, profmat[,i], ylab = colnames(profmat)[i])
}
This fails above a certain fixed gamma value, presumably because we aren't allowing $I(0)$ to decrease ...
fixing gamma
tl;dr with this set of constraints we ...
SIR_start0 <- c(beta = 7.89e-05, gamma = 4.47, I0 = 1, S0 = 61500, k = 50)
start5 <- c(xfun(opt2$par, "logr"), logI0 = log(1))
options(sir.grad_param = "R0m1_loggamma")
sir.nll4 <- function(p, loggamma = log(SIR_start0[["gamma"]])) {
sirpred <- sir.pred(c(xfun(p, c("logk", "loggamma")), loggamma = loggamma))
-sum(dnbinom(x = bombay2$mort[-1],
mu = sirpred,
size = exp(p["logk"]), log = TRUE))
}
sir.nll4(start5)
opt5 <- optim(start5,
fn = sir.nll4,
loggamma = log(4.77),
control = list(maxit = 1e6),
hessian = TRUE)
v5 <- solve(opt5$hessian)
sqrt(diag(v5))
Still not identical to what we were doing with fitode since we have
substituted logR0m1 for beta in the parameterization ...
no constraints
JD thinks that all four parameters should be identifiable ... ??
start6 <- c(opt5$par, loggamma = log(4.77))
options(sir.grad_param = "R0m1_loggamma")
sir.nll5 <- function(p) {
sirpred <- sir.pred(xfun(p, c("logk")))
-sum(dnbinom(x = bombay2$mort[-1],
mu = sirpred,
size = exp(p["logk"]), log = TRUE))
}
sir.nll5(start6)
opt6 <- optim(start6,
fn = sir.nll5,
control = list(maxit = 1e6),
hessian = TRUE)
data.frame(est = exp(opt6$par), se = sqrt(diag(solve(opt6$hessian)))*
exp(opt6$par))
with fitode
Not done ...
library(fitode)
SIR_model <- odemodel(
name="SIR model",
model=list(
S ~ - beta * S * I,
I ~ beta * S * I - gamma * I,
R ~ gamma * I
),
observation = list(
mort ~ dnbinom(mu = R, size = k)
),
diffnames="R",
initial=list(
S ~ N-I0,
I ~ I0,
R ~ 0
),
par=c("beta", "gamma", "S0", "I0", "k")
)
SIR_start2 <- do.call(trfun, as.list(c(opt2$par, logI0 = log(1))))
SIR_start0 <- c(beta = 7.89e-05, gamma = 4.47, I0 = 1, S0 = 61500, k = 50)
print(data.frame(s0 = SIR_start0, s2 = SIR_start2[names(SIR_start0)]), digits = 3)
Results from opt2 are similar to our original fitode-chosen starting values ...
SIR_fit1 <- fitode(
model = SIR_model,
data = bombay2,
fixed = list(gamma = SIR_start0[["gamma"]]),
start = SIR_start0,
tcol = "week")
coef(SIR_fit1)
To do
- work back toward doing this with
fitode
- fix $\gamma$ instead of N?
- How do estimates match with Bacaer, K&M, etc. ? (Is this where we came in?)
- explore multivariate sensitivity (e.g. 2D profiles for R0, r? should re-optimize I0 in each case because of correlation between I0, R)
- say something about biologically plausible range of $\gamma$? (go back to louseplague paper - how do we translate rat infectious periods to SIR infectious periods ... ? Louse plague prior fig gives $1/\gamma_r$ info; estimates from Tollenaere et al 2010 (4-10), Keeling and Gilligan 2000 (outlier, 18), Engelthaler et al 2000 (2-5)). Various other 'gamma' parameters are there, could dig through them ...)
parksw3/fitode documentation built on April 3, 2024, 7:45 a.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.
[ \newcommand{\rzero}{{\cal R}_0} ]
I'm going to try to understand why the Bombay data are so hard to fit/find reasonable confidence intervals for the parameters of an SIR model. We already know from Bacaer's paper that the closed-population, constant-transmission SIR is not a good biological model for the Bombay plague epidemic, but it is "SIR-shaped" and, for historical comparisons (i.e. to compare with Kermack and McKendrick's paper) we would like to be able to fit it with an SIR model even if that's not biologically sensible.
My thoughts on identifiability:
- It makes sense that $r = \beta-\gamma$ and $I(0)$ should be relatively easily identifiable (from the initial epidemic curve). If we assume we're starting from zero recovered individuals, then $S(0) \approx N$. Because $\rzero$ is monotonically related to final size (as a proportion of the population size), so $\rzero$ and $N$ (effective population size) should be jointly unidentifiable.
- If we reparameterized as ${r, \rzero, N, I_0}$ then we should be able to estimate ${r, \rzero, I(0)}$ if we fix $N$.
- On the other hand, fixing $\gamma$ should work in principle (once $\gamma$ is fixed the initial slope should determine beta, which will determine $\rzero$, so we can estimate $S(0)$ ($\approx N$) independently from final size ...), but ...
- estimating dispersion (negative binomial $k$ or Gaussian residual variance) should be independent of these \ldots
- We can also reason from the fact that Kermack and McKendrick's formula works (fitting an SIR curve with only three parameters [which are not the same as, but correspond to, epidemic peak height, peak time, and half-initial growth rate]), that we only have three free parameters to work with.
I'm going to work on this by brute force (optim() + deSolve() only) to avoid the possibilities that something weird is going on in bbmle and/or fitode.
To give ourselves the best chance of having everything work, we're going to parameterize this robustly: ${ \log(r), \log(\rzero-1), \log(I(0)) }$ (thought experiment: given the choices ${\beta, \gamma, r, \rzero, G}$ (where $G$ is generation time), which are the closest to independent conditional on an observed epidemic curve? This relates back to some of JD's stuff \ldots
We'll stick with the assumption of 100\% infection fatality rate for now.
library(deSolve) data("bombay", package = "fitode") ## incidence-based version bombay2 <- rbind( c(times=bombay$week[1] - diff(bombay$week)[1], mort=NA), bombay )
## set up sir.grad so we can switch parameterizations later on options(sir.grad_param = "R0m1_r") sir.grad <- function (t, x, params) { with(as.list(c(x, params)), { parameterization <- getOption("sir.grad_param", "R0m1_r") if (parameterization == "R0m1_r") { R0m1 <- exp(logR0m1) r <- exp(logr) gamma <- r/R0m1 beta <- (r + gamma)/N } else if (parameterization == "R0m1_loggamma") { R0m1 <- exp(logR0m1) gamma <- exp(loggamma) beta <- (R0m1+1)*gamma/N } else if (parameterization == "logbeta_loggamma") { gamma <- exp(loggamma) beta <- exp(logbeta) } dS.dt <- -beta*S*I dI.dt <- beta*S*I-gamma*I dR.dt <- gamma*I dxdt <- c(dS.dt,dI.dt,dR.dt) list(dxdt) }) } sir.pred <- function(p) { times <- bombay2$week I0 <- exp(p[["logI0"]]) S0 <- exp(p[["logS0"]]) ode.params <- c(xfun(p, c("logI0", "logS0", "logk")), N = S0) xstart <- c(S=S0, I=I0, R=0) out <- deSolve::ode(func=sir.grad, y=xstart, times=times, parms=ode.params ) |> as.data.frame() diff(out$R) } xfun <- function(x, n) x[setdiff(names(x), n)] sir.nll <- function(p, N = N) { ## set S0 to N, ignore I0 sirpred <- sir.pred(c(xfun(p, "logk"), logS0 = log(N))) -sum(dnbinom(x = bombay2$mort[-1], mu = sirpred, size = exp(p["logk"]), log = TRUE)) }
start0 <- c(logR0m1 = log(4), logr = log(log(2)/1.5), logI0 = log(0.01), logk = log(10)) p0 <- sir.pred(c(xfun(start0, "logk"), logS0 = 0)) plot(p0, log = "y")
I'm going to use 820,000 for $N$ (from Wikipedia this was the value from the 1891 census, five years earlier - obviously the effective population size could be different from that but ...)
bombay_pop <- 820e3
Eyeballing to adjust ($I(0)$ needs to be larger than we would have thought, to account for the delay from infection to mortality)
start1 <- start0 start1[["logI0"]] <- log(3) start1[["logR0m1"]] <- log(0.1) plot(mort ~ week, data = bombay2, log = "y", ylim = c(1,10000)) p1 <- sir.pred(c(xfun(start1, "logk"), logS0 = log(bombay_pop))) lines(p1)
This is actually pretty terrible, but should be close enough for starting (should spend more time thinking about how
sir.nll(start1, N = bombay_pop) opt1 <- optim(start1, fn = sir.nll, N = bombay_pop, control = list(maxit = 1e6), hessian = TRUE) plot(mort ~ week, data = bombay, log = "y") p2 <- sir.pred(c(xfun(opt1$par, "logk"), logS0 = log(bombay_pop))) lines(p2) v <- solve(opt1$hessian) sds <- sqrt(diag(v)) raw_pars <- cbind(est = opt1$par, lwr = opt1$par-1.96*sds, upr = opt1$par+1.96*sds) human_pars <- rbind(R0 = exp(raw_pars["logR0m1",])+1, r = exp(raw_pars["logr",]), I0 = exp(raw_pars["logI0",]), k = exp(raw_pars["logk",])) print(human_pars, digits = 2)
Parameters are correlated, but not pathologically so:
print(cov2cor(v), digits = 2)
Simulate outcomes (confidence intervals and prediction intervals by "population prediction interval", i.e. MVN simulation + nbinom observation error)
nsim <- 1000 nt <- length(p1) confmat <- predmat <- matrix(NA_real_, nrow = nt, ncol = nsim) set.seed(101) parmat <- MASS::mvrnorm(nsim, mu = opt1$par, Sigma = v) for (i in 1:nsim) { confmat[,i] <- sir.pred(c(xfun(parmat[i,], "logk"), logS0 = log(bombay_pop))) predmat[,i] <- rnbinom(nt, mu = confmat[,i], size = exp(opt1$par[["logk"]])) } pconf <- t(apply(confmat, 1, quantile, c(0.025, 0.975))) ppred <- t(apply(predmat, 1, quantile, c(0.025, 0.975)))
par(las=1, bty = "l") plot(mort ~ week, data = bombay, log = "y", ylim = c(3, 1500)) lines(p2) matlines(pconf, col = 4, lty = 2) matlines(ppred, col = 6, lty = 3) legend("topleft", lty = c(1, 2, 3), legend = c("best fit", "95% conf int", "95% pred int"), col = c(1, 4, 6))
What about translating back from parameters all the way to $\beta/\gamma$?
trfun <- function(logr, logR0m1, logI0, logk, logS0, N = log(bombay_pop)) { I0 <- exp(logI0) if (!missing(logS0)) { ## S0 <- exp(logS0) N <- S0 + I0 } else { S0 <- N - I0 } c(beta = exp(logr)*(1+1/exp(logR0m1))/N, gamma = exp(logr - logR0m1), S0 = S0, I0 = I0, k = exp(logk)) } trmat <- apply(parmat, 1, \(x) do.call(trfun, as.list(x))) print(cbind(est = do.call(trfun, as.list(opt1$par)), lwr = apply(trmat, 1, quantile, 0.025), upr = apply(trmat, 1, quantile, 0.975)), digits = 3)
The $\gamma$ estimate is obviously silly (95% CI approx. 70-80).
What if we fix $I(0)=1$ as done in the paper and allow $S(0)$ (effective population size) to vary instead?
sir.nll2 <- function(p, logI0 = log(1)) { sirpred <- sir.pred(c(xfun(p, "logk"), logI0 = logI0)) -sum(dnbinom(x = bombay2$mort[-1], mu = sirpred, size = exp(p["logk"]), log = TRUE)) } start2 <- c(xfun(start1, "logI0"), logS0 = log(bombay_pop)) opt2 <- optim(start2, fn = sir.nll2, logI0 = log(1), control = list(maxit = 1e6), hessian = TRUE) p2B <- sir.pred(c(xfun(opt2$par, "logk"), logI0 = log(1))) v2 <- solve(opt2$hessian) sds2 <- sqrt(diag(v)) raw_pars2 <- cbind(est = opt2$par, lwr = opt2$par-1.96*sds2, upr = opt2$par+1.96*sds2) human_pars2 <- rbind(R0 = exp(raw_pars2["logR0m1",])+1, r = exp(raw_pars2["logr",]), S0 = exp(raw_pars2["logS0",]), k = exp(raw_pars2["logk",])) print(human_pars2, digits = 3)
Stronger correlations (although still not insane)
print(cov2cor(v2), digits = 2)
nsim <- 1000 nt <- length(p1) confmat2 <- predmat2 <- matrix(NA_real_, nrow = nt, ncol = nsim) set.seed(101) parmat2 <- MASS::mvrnorm(nsim, mu = opt2$par, Sigma = v) for (i in 1:nsim) { confmat2[,i] <- sir.pred(c(xfun(parmat2[i,], "logk"), logI0 = log(1))) predmat2[,i] <- rnbinom(nt, mu = confmat2[,i], size = exp(opt2$par[["logk"]])) } ## DRY ... pconf2 <- t(apply(confmat2, 1, quantile, c(0.025, 0.975))) ppred2 <- t(apply(predmat2, 1, quantile, c(0.025, 0.975)))
par(las=1, bty = "l") plot(mort ~ week, data = bombay, log = "y", ylim = c(3, 1500)) lines(p2B) lines(p2, col = 2) matlines(pconf2, col = 4, lty = 2) matlines(ppred2, col = 6, lty = 3) legend("bottom", lty = c(1, 2, 3), legend = c("best fit (fix I0/gamma)", "best fit (fix N/?)", "95% conf int", "95% pred int"), col = c(1, 2, 4, 6))
trfun2 <- function(logr, logR0m1, logS0, logk) { c(beta = exp(logr)*(1+1/exp(logR0m1))/exp(logS0), gamma = exp(logr - logR0m1), S0 = exp(logS0), I0 = 1, k = exp(logk)) } trmat2 <- apply(parmat2, 1, \(x) do.call(trfun2, as.list(x))) print(cbind(est = do.call(trfun2, as.list(opt2$par)), lwr = apply(trmat2, 1, quantile, 0.025), upr = apply(trmat2, 1, quantile, 0.975)), digits = 2)
Correlations among beta/gamma/S0/k aren't too bad ...
cor(t(trmat2[-4,]))
The value of $\gamma$ can compensate for a fixed value of $I(0)$ because we are measuring incidence (i.e. the initial observed mortality rate is approximately $\gamma I(0)$).
Do we need to worry about incidence vs prevalence when translating K-M parameters to biological parameters ... ?
Gamma vs. N
options(sir.grad_param = "R0m1_loggamma") start3 <- c(xfun(start1, c("logI0", "logr")), loggamma = log(4.77), logS0 = log(bombay_pop)) sir.nll2(start3, logI0 = log(1)) opt3 <- optim(start3, fn = sir.nll2, logI0 = log(1), control = list(maxit = 1e6), hessian = TRUE) exp(opt3$par) ## NLL with fixed value of gamma sir.nll3 <- function(p, logI0 = log(1), loggamma = log(4.77)) { sirpred <- sir.pred(c(xfun(p, "logk"), logI0 = logI0, loggamma = loggamma)) -sum(dnbinom(x = bombay2$mort[-1], mu = sirpred, size = exp(p["logk"]), log = TRUE)) } opt4 <- optim(xfun(start3, "loggamma"), fn = sir.nll3, logI0 = log(1), loggamma = log(4.77), control = list(maxit = 1e6), hessian = FALSE) lgammavec <- seq(log(1), log(100), length = 51) profmat <- matrix(NA_real_, nrow = length(lgammavec), ncol = 4, dimnames = list(NULL, c(names(opt4$par), "nll"))) for (i in seq_along(lgammavec)) { pstart <- start3 ## cat(i, lgammavec[i], "\n") oo <- try(optim( xfun(start3, "loggamma"), fn = sir.nll3, logI0 = log(1), loggamma = lgammavec[i], control = list(maxit = 1e6), hessian = FALSE), silent = TRUE ) if (!is(oo, "try-error")) { pstart[names(oo$par)] <- oo$par profmat[i,] <- c(oo$par, oo$value) } } par(mfrow = c(2,2), las = 1, bty = "l", mar = c(3,4,1,1), mgp = c(2, 1, 0)) for (i in 1:4) { plot(lgammavec, profmat[,i], ylab = colnames(profmat)[i]) }
This fails above a certain fixed gamma value, presumably because we aren't allowing $I(0)$ to decrease ...
fixing gamma
tl;dr with this set of constraints we ...
SIR_start0 <- c(beta = 7.89e-05, gamma = 4.47, I0 = 1, S0 = 61500, k = 50) start5 <- c(xfun(opt2$par, "logr"), logI0 = log(1)) options(sir.grad_param = "R0m1_loggamma") sir.nll4 <- function(p, loggamma = log(SIR_start0[["gamma"]])) { sirpred <- sir.pred(c(xfun(p, c("logk", "loggamma")), loggamma = loggamma)) -sum(dnbinom(x = bombay2$mort[-1], mu = sirpred, size = exp(p["logk"]), log = TRUE)) } sir.nll4(start5) opt5 <- optim(start5, fn = sir.nll4, loggamma = log(4.77), control = list(maxit = 1e6), hessian = TRUE) v5 <- solve(opt5$hessian) sqrt(diag(v5))
Still not identical to what we were doing with fitode since we have
substituted logR0m1 for beta in the parameterization ...
no constraints
JD thinks that all four parameters should be identifiable ... ??
start6 <- c(opt5$par, loggamma = log(4.77)) options(sir.grad_param = "R0m1_loggamma") sir.nll5 <- function(p) { sirpred <- sir.pred(xfun(p, c("logk"))) -sum(dnbinom(x = bombay2$mort[-1], mu = sirpred, size = exp(p["logk"]), log = TRUE)) } sir.nll5(start6) opt6 <- optim(start6, fn = sir.nll5, control = list(maxit = 1e6), hessian = TRUE) data.frame(est = exp(opt6$par), se = sqrt(diag(solve(opt6$hessian)))* exp(opt6$par))
with fitode
Not done ...
library(fitode) SIR_model <- odemodel( name="SIR model", model=list( S ~ - beta * S * I, I ~ beta * S * I - gamma * I, R ~ gamma * I ), observation = list( mort ~ dnbinom(mu = R, size = k) ), diffnames="R", initial=list( S ~ N-I0, I ~ I0, R ~ 0 ), par=c("beta", "gamma", "S0", "I0", "k") ) SIR_start2 <- do.call(trfun, as.list(c(opt2$par, logI0 = log(1)))) SIR_start0 <- c(beta = 7.89e-05, gamma = 4.47, I0 = 1, S0 = 61500, k = 50) print(data.frame(s0 = SIR_start0, s2 = SIR_start2[names(SIR_start0)]), digits = 3)
Results from opt2 are similar to our original fitode-chosen starting values ...
SIR_fit1 <- fitode( model = SIR_model, data = bombay2, fixed = list(gamma = SIR_start0[["gamma"]]), start = SIR_start0, tcol = "week") coef(SIR_fit1)
To do
- work back toward doing this with
fitode - fix $\gamma$ instead of N?
- How do estimates match with Bacaer, K&M, etc. ? (Is this where we came in?)
- explore multivariate sensitivity (e.g. 2D profiles for R0, r? should re-optimize I0 in each case because of correlation between I0, R)
- say something about biologically plausible range of $\gamma$? (go back to louseplague paper - how do we translate rat infectious periods to SIR infectious periods ... ? Louse plague prior fig gives $1/\gamma_r$ info; estimates from Tollenaere et al 2010 (4-10), Keeling and Gilligan 2000 (outlier, 18), Engelthaler et al 2000 (2-5)). Various other 'gamma' parameters are there, could dig through them ...)
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.