In bbolker/McMasterPandemic: Pandemic Model
knitr::opts_chunk$set(echo = TRUE)
Notes on other references: (1) testing_flow.md(I am using the notation in this doc)
On Math and Weighting
Model and Definitions
- Assume there are 3 groups of individuals: $S$, $I$ and $R$.
- Each group splits into 4 subgroups to incorporate the testing mechanism, so overall the model framework has 12 compartments as follows. $S_u,$ $S_n$, $S_p$, $S_t$, $I_u$, $I_n$, $I_p$, $I_t$, $R_u$, $R_n$, $R_p$, $R_t$.
- Further, assume that (1) a perfectly specific test, i.e. no false positives, then $S_p$ and $S_t$ are always zero (2) there are two 'integrator' or 'accumulator' compartments, $N$ and $P$, to collect cumulative reported negative or positive tests.
Per Capita Rate of Flow from Untested to Positive or Negative Compartments
- $t$: testing intensity with dimension of $1/time$.
-
$\omega$: the rate of onward flow from the awaiting positive compartment, $p$, to reported/tested compartment, $t$, or from awaiting negative compartment, $n$, back to $u$. It has units of $1/time$.
-
Note: The complicated part is the flow rate from $Z_u$ to $Z_n$, $Z_p$. The (total, not per capita) flow from $Z_u$ to $Z_n$ is $t (1-P_Z) F(W_Z) Z_u$, from $Z_u$ to $Z_p$ is $t P_Z F(W_Z) Z_u$. Where, $P_Z$ is the probability of positive test.
-
At the simplest $F(W_Z)=W_Z$. The flows from $Z_n$ to $Z_u$ and $Z_p$ to $Z_t$ are $\omega Z_n$ and $\omega Z_p$ respectively. We are using a dimensionless scaling where $F(W_Z) = (W_Z \cdot N_0)/\sum_{Z'} \left( W_{Z'} {Z'}_{u} \right)$, i.e. the denominator is the sum of the product of the weights and the current occupancy of the untested compartments.
Model
Let the force of infection $\Lambda=\beta \frac{(I_u+\eta_w I_n+\eta_w I_p+ \eta_t I_t)}{N_0}$, where $\eta$ is the isolation parameter with the following assumptions:
(1) All susceptible individuals are exposed to all infectious individuals regardless of the group factor.
(2) $\eta_t<\eta_w$, i.e., the awaiting individuals for test results have higher transmission probability than the reported individuals.
(3) The scaling parameter is defined as $sc = t N_0/(W_s S_u + W_i I_u + W_r R_u)$ and $F_s = sc \times W_s$, $F_i=sc \times W_i$ and $F_r = sc\times W_r$.
The model is
\begin{align}
\label{mod:1}
1)\ d S_u/dt &= -\Lambda S_u - F_s S_u + \omega S_n \
2)\ d S_n/dt &= -\Lambda S_n + F_s S_u - \omega S_n \
3)\ d I_u/dt &= \Lambda S_u - F_i I_u + \omega I_n - \gamma I_u \
4)\ d I_n/dt &= \Lambda S_n + (1-P_i) F_i I_u - \omega I_n -\gamma I_n \
5)\ d I_p/dt &= P_i F_i I_u - \omega I_p -\gamma I_p \
6)\ d I_t/dt &= \omega I_p - \gamma I_t \
7)\ d R_u/dt &= \gamma I_u - F_r R_u + \omega R_n \
8)\ d R_n/dt &= \gamma I_n + (1-P_r) F_r R_u - \omega R_n \
9)\ d R_p/dt &= \gamma I_p + P_r F_r R_u - \omega R_p \
10)\ d R_t/dt&= \gamma I_t + \omega R_p \
11)\ dN/dt &= \omega (S_n + I_n + R_n) \
12)\ dP/dt &= \omega(I_p + R_p).
\end{align}
library(deSolve)
library(ggplot2)
library(tidyr)
library(McMasterPandemic)
library(directlabels)
unpack <- McMasterPandemic::unpack
sir.model <- function(time,state,params){
unpack(as.list(c(state,params)))
## Force of Infection
Lambda <- beta * (I_u + eta_w*(I_n+I_p) + eta_t*I_t)/N0
## scaling the weights
sc <- t*N0/(W_s*S_u+W_i*I_u+W_r*R_u)
# sc <- t*N0
F_s <- sc*W_s
F_i <- sc*W_i
F_r <- sc*W_r
dS_u.dt <- -Lambda*S_u - (1-P_s) * F_s * S_u + omega * S_n
dS_n.dt <- -Lambda*S_n + (1-P_s) * F_s * S_u - omega * S_n
dI_u.dt <- Lambda*S_u + omega * I_n - F_i * I_u - gamma * I_u
dI_n.dt <- Lambda*S_n + (1-P_i) * F_i * I_u - omega * I_n - gamma * I_n
dI_p.dt <- P_i * F_i * I_u - omega * I_p - gamma * I_p
dI_t.dt <- omega * I_p - gamma * I_t
dR_u.dt <- gamma * I_u + omega * R_n - F_r * R_u
dR_n.dt <- gamma * I_n + (1-P_r) * F_r * R_u - omega * R_n
dR_p.dt <- gamma * I_p + P_r * F_r * R_u - omega * R_p
dR_t.dt <- gamma * I_t + omega * R_p
dN.dt <- omega * (S_n + I_n + R_n)
dP.dt <- omega *(I_p + R_p)
# return the rate of change
dxdt <- c(dS_u.dt,dS_n.dt,dI_u.dt,dI_n.dt,dI_p.dt,dI_t.dt,dR_u.dt,dR_n.dt,dR_p.dt,dR_t.dt,dN.dt,dP.dt)
## }
return(list(dxdt))
}
params <- c(N0=1000000, beta=1,gamma=1/3, omega=0.25, t=0.01,
W_s=0.01, W_i=1, W_r=0.01, #weights
P_s=0, P_i=0.5, P_r=0.5, #prob of waiting for being tested positive
eta_w=0.02, eta_t=0.01 #isolation parameter
)
class(params) <- "params_pansim"
state_init <- c(S_u=params[["N0"]], S_n=0,
I_u=1,I_n=0,I_p=0,I_t=0,
R_u=0,R_n=0,R_p=0,R_t=0,
N=0,P=0)
# Disease free equilibrium
Sn_dfe <-function(params){
unpack(as.list(params))
return((t*(1-P_s)*N0)/omega)
}
Su_dfe<- params[["N0"]]-Sn_dfe(params)
state_dfe <- c(S_u=Su_dfe, S_n=Sn_dfe(params),
I_u=0,I_n=0,I_p=0,I_t=0,
R_u=0,R_n=0,R_p=0,R_t=0,
N=0,P=0)
d <- seq(0,200,by=0.1)
# Conditional stop for desolver
rootfun <- function (time, state,params) {
unpack(as.list(c(state,params)))
return(I_u - 1) }
Test gradient
g1 <- sir.model(time=0,state=state_init, params=params)
all.equal(sum(g1[[1]]),0)
## debug(sir.model)
## g1 <- sir.model(time=0,state=state_init, params=params)
out <- as.data.frame(
ode(
func=sir.model,
y=state_init,
times= d,
parms=update(params,beta=2),
# atol = 1e+1, rtol = 1e+1
# rootfun = rootfun,
# method="lsodar"
)
)
out2 <- out %>%
pivot_longer(c(S_u,S_n,I_u,I_n,I_p,I_t, R_u,R_n,R_p,R_t), names_to = "compartment", values_to = "value")
gg1 <- (ggplot(data=out2, aes(x=time,y=value,col=compartment))
+ geom_line()
+ scale_y_log10(limits=c(0.1,params[["N0"]]))
## + scale_x_continuous(limits=c(0,240))
)
direct.label(gg1,"last.bumpup")
d_long <- seq(0,300,by=0.1)
## capturing badly behaved integrator warning messages
cc <- capture.output(
out_long <- as.data.frame(
ode(
func=sir.model,
y=state_init,
times= d_long,
parms=update(params,beta=2)
# rootfun = rootfun,
# method="lsodar"
)
))
out_long2 <- out_long %>%
pivot_longer(c(S_u,S_n,I_u,I_n,I_p,I_t, R_u,R_n,R_p,R_t),
names_to = "compartment", values_to = "value")
direct.label(gg1 %+% out_long2,"last.bumpup")
Steady States
The Disease-Free Equilibrium, namely $E_1$, for the above SIR model is
\label{exp:dfe} $S_n^= \frac{t (1-P_s)}{\omega} N_0, S_u^= N_0-S_n^*$, and $I_j=R_j=0 \ \forall j$.
Basic Reproduction Number
The next generation matrix, $G = F V^{-1}$, is in the lower triangular form with $F$ and $G$ as follows.
$$R_0= (A \times S_u^ + B \times S_n^) \times C,$$ where
$A=\gamma(\omega+\gamma) + (\gamma \eta_w + \omega \eta_t P_i) F_i$,
$B=\big(\omega+(F_i+\gamma)\eta_w\big) \gamma+\frac{(\eta_w \gamma+ \eta_t\omega) \omega P_i F_i }{\omega+\gamma}$ and
$C=\frac{\beta/\gamma}{N_0 (\gamma(\omega+\gamma)+F_i(\gamma+\omega P_i))}$.
\begin{equation}\label{FV}
F= \beta/N_0 \left[ \begin {array}{cccc}
S_u&\eta_w\,S_u&\eta_w\,S_u&\eta_t\,S_u\
S_n&\eta_w\,S_n&\eta_w\,S_n&\eta_t\,S_n\
0&0&0&0\
0&0&0&0
\end {array} \right],
V=
\left[ \begin {array}{cccc}
F_i+\gamma&-\omega&0&0\
-(1-P_i)F_i&\omega+\gamma&0&0\
-P_i F_i&0&\omega+\gamma&0\
0&0&-\omega&\gamma
\end {array} \right],
\end{equation}
\begin{equation}
G = \left[ \begin {array}{cc}
G_{11}&G_{12}\
0&0
\end {array} \right], \text{ where } \
G_{11} =C
\left[\begin {array}{cc}
A\,S_u & B\,S_u\
A\,S_n & B\,S_n
\end {array}\right],
\end{equation}
The Eigenvalues of $G_{11}$, are 0 and $R_0$.
F_i <- function(state,params){
unpack(as.list(c(state,params)))
return(t*N0*W_i/(W_s*S_u+W_i*I_u+W_r*R_u))
}
# just a test
F_i(state=state_dfe,params=update(params,t=0))
R0<-function(params){
unpack(as.list(params))
Sn_dfe<- Sn_dfe(params)
Su_dfe<- N0-Sn_dfe
F_i <-t*N0*W_i/(W_s*Su_dfe) #at DFE
# F_i <-t*N0*W_i
A <-gamma*(omega+gamma+eta_w*F_i) + omega*eta_t*P_i*F_i
B <-(omega+eta_w*(F_i+gamma))*gamma+(eta_w*gamma+eta_t*omega)*(omega*P_i*F_i)/(omega+gamma)
C <- beta/(N0*gamma*(gamma*(omega+gamma)+F_i*(gamma+omega*P_i)))
return((A*Su_dfe+B*Sn_dfe)*C)
}
# R0<1:
R0(params=params)
# R0>1:
R0(params=update(params,beta=2))
Can we recover $R_0=\beta/\gamma$ of the simple SIR? -yes
# checking if simple R_0 for SIR can be recovered?
# Case #1: t=0
all.equal(R0(params=update(params,t=0)), params[["beta"]]/params[["gamma"]])
# Case #2:
all.equal(R0(params=update(params,eta_t=1,eta_w=1)), params[["beta"]]/params[["gamma"]])
On Behaviour of Testing at DFE Boundry
(1) Note that once the only untested people are susceptibles, the FoI will be $\Lambda=0$, testing rate $F_s=t N_0/S_u$. Thus, eq(1) of the model will be $d S_u/dt = - t N_0 + \omega S_n$ which is a linear rate of leaving the $S_u$ compartment. That means that the testing rate will be no longer dependent on $S_u$.
Contour plot of $R_0=1$ and interpretations
# specify the ranges
beta <- 0.8 ## set beta high enough to allow R0>1 in worst case
eta_w <- seq(0,1, length.out=5)
eta_t <- seq(0,1, length.out=5)
t <- seq(0,0.01, length.out=20)
omega <- seq(0.1,1 , length.out=20) #note omega must be > t
df1 <- expand.grid(N0=params[["N0"]],beta=beta,gamma=params[["gamma"]],omega=omega,t=t,
W_s=params[["W_s"]],W_i=params[["W_i"]],W_r=params[["W_r"]],
P_s=params[["P_s"]],P_i=params[["P_i"]],P_r=params[["P_r"]],
eta_w=eta_w,eta_t=eta_t)
df2<- data.frame(df1,R0=apply(df1,1,function(params_in)R0(params=params_in)),
Fi= apply(df1,1,function(params_in)F_i(state=state_dfe,params=params_in)))
which(df2[]<0)
## hacked (not very robustly) to chain label_both() and label_parsed() ...
label_special <- function (labels, multi_line=FALSE, sep = "== ") {
value <- ggplot2:::label_value(labels, multi_line = multi_line)
variable <- ggplot2:::label_variable(labels, multi_line = multi_line)
variable[[1]] <- gsub("_(.)$","[\\1]",variable[[1]])
## not using multiple faceting variables on each margin, so forget paste
out <- Map(paste, variable, value, sep = sep)
out[[1]] <- lapply(out[[1]], function(expr) c(parse(text=expr)))
out
}
p1 <- ggplot(df2,aes(x=omega,y=t,z=R0))+ theme_bw() +
theme(panel.spacing=grid::unit(0,"lines"))
(p1
+ geom_contour_filled(breaks=pretty(unique(df2$R0),5),
alpha=0.5) ## breaks=0.2)
+ geom_contour(breaks=1,alpha=0.5,colour="black")
+ facet_grid(eta_w~eta_t, labeller=label_special)
+ scale_x_continuous(expand=expansion(c(0,0)), n.breaks=3)
+ scale_y_continuous(expand=expansion(c(0,0)), n.breaks=3)
+ scale_fill_viridis_d(name=parse(text="R[0]"))
)
Let call the $F_i$ the critical testing rate in the above plot.
Propositions
- $\partial{R_0}/\partial{\eta_t} \geq 0$ and $\partial{R_0}/\partial{\eta_w} \geq 0$.
- $\partial{R_0}/\partial{t} \leq 0$; this is concluded from the Taylor expansion of $R_0$ around $t=0$ which is
$R_0 \approx \beta/\gamma + \frac{\beta t}{\omega (\omega+\gamma) \gamma^2 W_s} \Big(\gamma(\eta_w-1)(\gamma W_s+\omega W_i) + (\eta_t -1)P_iW_i \omega^2 \Big) + \mathcal{O}(t^2)$.
- Based on the above Taylor expansion around $t=0$,
$$\partial{R_0}/\partial{\omega}= \frac{-\beta t}{\gamma W_s\omega^2 (\gamma+\omega)^2} (a \omega^2 + b \omega + c), $$
where $a=(\eta_w-1)W_i-(\eta_t-1)P_iW_i$, $b=2(\eta_w-1)\gamma W_s$ and $c=(\eta_w-1)\gamma^2 W_s$.
Given that $\eta_t\leq \eta_w$, it is inferred that $a\geq 0$, $b\leq 0$ and $c \leq 0$.
Let $\omega_1<0$ and $\omega_2>0$ be the roots of the quadratic expression in $\partial{R_0}/\partial{\omega}$. Thus, $\partial{R_0}/\partial{\omega}<0$ for $0<\omega<\omega_2$ and $\partial{R_0}/\partial{\omega}>0$ for $\omega>\omega_2$.
Interpretation of the R0 contour plot
Note that $\eta_t\leq \eta_w$, i.e., reported people have a lower transmission probability than awaiting individuals. The perfect isolation/quarantine occurs when $\eta_w=\eta_t=0$ (see the upper left corner panel). Thus, we only consider the lower-triangle panels in the plot [not sure how to extract them in ggplot?].
(1) Looking at the first column of panels, there seems to be a compensatory relation between $t$ and $\omega$. i.e., higher $t$ requires smaller $\omega$ to keep R0 lower than 1. Also, for a fixed $\omega$, increasing $\eta_w$ requires higher $t$ to achive the critical R0.
Literature Review
\cite{bergstrom2020frequency}
(1) Model, assumptions: They developed a function, namely expected expoisour $E(C,\tau)$, to approximate trade-offs between the frequency of testing, n, the sensitivity of testing, q, and the delay between
testing and results, d. This function is explicitly derived and was connected the effective reproduction number $R=R_0 S$, where $S$ is the proportion of population susciptable.
assumption that transmission rates are a step function: individuals who
have COVID go from non-infectious to fully infectious instantaneously,
and remain fully infectious until they are no longer able to transmit disease. Test sensitivity takes the same form over the course of infection.
More sophisticated models could allow varying infectiousness and varying
sensitivity over time, as in \citep{larremore2020test}.
bbolker/McMasterPandemic documentation built on Aug. 25, 2024, 6:35 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.
knitr::opts_chunk$set(echo = TRUE)
Notes on other references: (1) testing_flow.md(I am using the notation in this doc)
On Math and Weighting
Model and Definitions
- Assume there are 3 groups of individuals: $S$, $I$ and $R$.
- Each group splits into 4 subgroups to incorporate the testing mechanism, so overall the model framework has 12 compartments as follows. $S_u,$ $S_n$, $S_p$, $S_t$, $I_u$, $I_n$, $I_p$, $I_t$, $R_u$, $R_n$, $R_p$, $R_t$.
- Further, assume that (1) a perfectly specific test, i.e. no false positives, then $S_p$ and $S_t$ are always zero (2) there are two 'integrator' or 'accumulator' compartments, $N$ and $P$, to collect cumulative reported negative or positive tests.
Per Capita Rate of Flow from Untested to Positive or Negative Compartments
- $t$: testing intensity with dimension of $1/time$.
-
$\omega$: the rate of onward flow from the awaiting positive compartment, $p$, to reported/tested compartment, $t$, or from awaiting negative compartment, $n$, back to $u$. It has units of $1/time$.
-
Note: The complicated part is the flow rate from $Z_u$ to $Z_n$, $Z_p$. The (total, not per capita) flow from $Z_u$ to $Z_n$ is $t (1-P_Z) F(W_Z) Z_u$, from $Z_u$ to $Z_p$ is $t P_Z F(W_Z) Z_u$. Where, $P_Z$ is the probability of positive test.
-
At the simplest $F(W_Z)=W_Z$. The flows from $Z_n$ to $Z_u$ and $Z_p$ to $Z_t$ are $\omega Z_n$ and $\omega Z_p$ respectively. We are using a dimensionless scaling where $F(W_Z) = (W_Z \cdot N_0)/\sum_{Z'} \left( W_{Z'} {Z'}_{u} \right)$, i.e. the denominator is the sum of the product of the weights and the current occupancy of the untested compartments.
Model
Let the force of infection $\Lambda=\beta \frac{(I_u+\eta_w I_n+\eta_w I_p+ \eta_t I_t)}{N_0}$, where $\eta$ is the isolation parameter with the following assumptions:
(1) All susceptible individuals are exposed to all infectious individuals regardless of the group factor. (2) $\eta_t<\eta_w$, i.e., the awaiting individuals for test results have higher transmission probability than the reported individuals. (3) The scaling parameter is defined as $sc = t N_0/(W_s S_u + W_i I_u + W_r R_u)$ and $F_s = sc \times W_s$, $F_i=sc \times W_i$ and $F_r = sc\times W_r$.
The model is \begin{align} \label{mod:1} 1)\ d S_u/dt &= -\Lambda S_u - F_s S_u + \omega S_n \ 2)\ d S_n/dt &= -\Lambda S_n + F_s S_u - \omega S_n \ 3)\ d I_u/dt &= \Lambda S_u - F_i I_u + \omega I_n - \gamma I_u \ 4)\ d I_n/dt &= \Lambda S_n + (1-P_i) F_i I_u - \omega I_n -\gamma I_n \ 5)\ d I_p/dt &= P_i F_i I_u - \omega I_p -\gamma I_p \ 6)\ d I_t/dt &= \omega I_p - \gamma I_t \ 7)\ d R_u/dt &= \gamma I_u - F_r R_u + \omega R_n \ 8)\ d R_n/dt &= \gamma I_n + (1-P_r) F_r R_u - \omega R_n \ 9)\ d R_p/dt &= \gamma I_p + P_r F_r R_u - \omega R_p \ 10)\ d R_t/dt&= \gamma I_t + \omega R_p \ 11)\ dN/dt &= \omega (S_n + I_n + R_n) \ 12)\ dP/dt &= \omega(I_p + R_p). \end{align}
library(deSolve) library(ggplot2) library(tidyr) library(McMasterPandemic) library(directlabels) unpack <- McMasterPandemic::unpack
sir.model <- function(time,state,params){ unpack(as.list(c(state,params))) ## Force of Infection Lambda <- beta * (I_u + eta_w*(I_n+I_p) + eta_t*I_t)/N0 ## scaling the weights sc <- t*N0/(W_s*S_u+W_i*I_u+W_r*R_u) # sc <- t*N0 F_s <- sc*W_s F_i <- sc*W_i F_r <- sc*W_r dS_u.dt <- -Lambda*S_u - (1-P_s) * F_s * S_u + omega * S_n dS_n.dt <- -Lambda*S_n + (1-P_s) * F_s * S_u - omega * S_n dI_u.dt <- Lambda*S_u + omega * I_n - F_i * I_u - gamma * I_u dI_n.dt <- Lambda*S_n + (1-P_i) * F_i * I_u - omega * I_n - gamma * I_n dI_p.dt <- P_i * F_i * I_u - omega * I_p - gamma * I_p dI_t.dt <- omega * I_p - gamma * I_t dR_u.dt <- gamma * I_u + omega * R_n - F_r * R_u dR_n.dt <- gamma * I_n + (1-P_r) * F_r * R_u - omega * R_n dR_p.dt <- gamma * I_p + P_r * F_r * R_u - omega * R_p dR_t.dt <- gamma * I_t + omega * R_p dN.dt <- omega * (S_n + I_n + R_n) dP.dt <- omega *(I_p + R_p) # return the rate of change dxdt <- c(dS_u.dt,dS_n.dt,dI_u.dt,dI_n.dt,dI_p.dt,dI_t.dt,dR_u.dt,dR_n.dt,dR_p.dt,dR_t.dt,dN.dt,dP.dt) ## } return(list(dxdt)) }
params <- c(N0=1000000, beta=1,gamma=1/3, omega=0.25, t=0.01, W_s=0.01, W_i=1, W_r=0.01, #weights P_s=0, P_i=0.5, P_r=0.5, #prob of waiting for being tested positive eta_w=0.02, eta_t=0.01 #isolation parameter ) class(params) <- "params_pansim" state_init <- c(S_u=params[["N0"]], S_n=0, I_u=1,I_n=0,I_p=0,I_t=0, R_u=0,R_n=0,R_p=0,R_t=0, N=0,P=0) # Disease free equilibrium Sn_dfe <-function(params){ unpack(as.list(params)) return((t*(1-P_s)*N0)/omega) } Su_dfe<- params[["N0"]]-Sn_dfe(params) state_dfe <- c(S_u=Su_dfe, S_n=Sn_dfe(params), I_u=0,I_n=0,I_p=0,I_t=0, R_u=0,R_n=0,R_p=0,R_t=0, N=0,P=0) d <- seq(0,200,by=0.1) # Conditional stop for desolver rootfun <- function (time, state,params) { unpack(as.list(c(state,params))) return(I_u - 1) }
Test gradient
g1 <- sir.model(time=0,state=state_init, params=params) all.equal(sum(g1[[1]]),0) ## debug(sir.model) ## g1 <- sir.model(time=0,state=state_init, params=params)
out <- as.data.frame( ode( func=sir.model, y=state_init, times= d, parms=update(params,beta=2), # atol = 1e+1, rtol = 1e+1 # rootfun = rootfun, # method="lsodar" ) ) out2 <- out %>% pivot_longer(c(S_u,S_n,I_u,I_n,I_p,I_t, R_u,R_n,R_p,R_t), names_to = "compartment", values_to = "value")
gg1 <- (ggplot(data=out2, aes(x=time,y=value,col=compartment)) + geom_line() + scale_y_log10(limits=c(0.1,params[["N0"]])) ## + scale_x_continuous(limits=c(0,240)) ) direct.label(gg1,"last.bumpup")
d_long <- seq(0,300,by=0.1) ## capturing badly behaved integrator warning messages cc <- capture.output( out_long <- as.data.frame( ode( func=sir.model, y=state_init, times= d_long, parms=update(params,beta=2) # rootfun = rootfun, # method="lsodar" ) )) out_long2 <- out_long %>% pivot_longer(c(S_u,S_n,I_u,I_n,I_p,I_t, R_u,R_n,R_p,R_t), names_to = "compartment", values_to = "value")
direct.label(gg1 %+% out_long2,"last.bumpup")
Steady States
The Disease-Free Equilibrium, namely $E_1$, for the above SIR model is \label{exp:dfe} $S_n^= \frac{t (1-P_s)}{\omega} N_0, S_u^= N_0-S_n^*$, and $I_j=R_j=0 \ \forall j$.
Basic Reproduction Number
The next generation matrix, $G = F V^{-1}$, is in the lower triangular form with $F$ and $G$ as follows.
$$R_0= (A \times S_u^ + B \times S_n^) \times C,$$ where $A=\gamma(\omega+\gamma) + (\gamma \eta_w + \omega \eta_t P_i) F_i$, $B=\big(\omega+(F_i+\gamma)\eta_w\big) \gamma+\frac{(\eta_w \gamma+ \eta_t\omega) \omega P_i F_i }{\omega+\gamma}$ and $C=\frac{\beta/\gamma}{N_0 (\gamma(\omega+\gamma)+F_i(\gamma+\omega P_i))}$.
\begin{equation}\label{FV}
F= \beta/N_0 \left[ \begin {array}{cccc}
S_u&\eta_w\,S_u&\eta_w\,S_u&\eta_t\,S_u\
S_n&\eta_w\,S_n&\eta_w\,S_n&\eta_t\,S_n\
0&0&0&0\
0&0&0&0
\end {array} \right],
V=
\left[ \begin {array}{cccc}
F_i+\gamma&-\omega&0&0\
-(1-P_i)F_i&\omega+\gamma&0&0\
-P_i F_i&0&\omega+\gamma&0\
0&0&-\omega&\gamma
\end {array} \right],
\end{equation}
\begin{equation} G = \left[ \begin {array}{cc} G_{11}&G_{12}\ 0&0 \end {array} \right], \text{ where } \ G_{11} =C \left[\begin {array}{cc} A\,S_u & B\,S_u\ A\,S_n & B\,S_n \end {array}\right], \end{equation}
The Eigenvalues of $G_{11}$, are 0 and $R_0$.
F_i <- function(state,params){ unpack(as.list(c(state,params))) return(t*N0*W_i/(W_s*S_u+W_i*I_u+W_r*R_u)) } # just a test F_i(state=state_dfe,params=update(params,t=0)) R0<-function(params){ unpack(as.list(params)) Sn_dfe<- Sn_dfe(params) Su_dfe<- N0-Sn_dfe F_i <-t*N0*W_i/(W_s*Su_dfe) #at DFE # F_i <-t*N0*W_i A <-gamma*(omega+gamma+eta_w*F_i) + omega*eta_t*P_i*F_i B <-(omega+eta_w*(F_i+gamma))*gamma+(eta_w*gamma+eta_t*omega)*(omega*P_i*F_i)/(omega+gamma) C <- beta/(N0*gamma*(gamma*(omega+gamma)+F_i*(gamma+omega*P_i))) return((A*Su_dfe+B*Sn_dfe)*C) } # R0<1: R0(params=params) # R0>1: R0(params=update(params,beta=2))
Can we recover $R_0=\beta/\gamma$ of the simple SIR? -yes
# checking if simple R_0 for SIR can be recovered? # Case #1: t=0 all.equal(R0(params=update(params,t=0)), params[["beta"]]/params[["gamma"]]) # Case #2: all.equal(R0(params=update(params,eta_t=1,eta_w=1)), params[["beta"]]/params[["gamma"]])
On Behaviour of Testing at DFE Boundry
(1) Note that once the only untested people are susceptibles, the FoI will be $\Lambda=0$, testing rate $F_s=t N_0/S_u$. Thus, eq(1) of the model will be $d S_u/dt = - t N_0 + \omega S_n$ which is a linear rate of leaving the $S_u$ compartment. That means that the testing rate will be no longer dependent on $S_u$.
Contour plot of $R_0=1$ and interpretations
# specify the ranges beta <- 0.8 ## set beta high enough to allow R0>1 in worst case eta_w <- seq(0,1, length.out=5) eta_t <- seq(0,1, length.out=5) t <- seq(0,0.01, length.out=20) omega <- seq(0.1,1 , length.out=20) #note omega must be > t df1 <- expand.grid(N0=params[["N0"]],beta=beta,gamma=params[["gamma"]],omega=omega,t=t, W_s=params[["W_s"]],W_i=params[["W_i"]],W_r=params[["W_r"]], P_s=params[["P_s"]],P_i=params[["P_i"]],P_r=params[["P_r"]], eta_w=eta_w,eta_t=eta_t) df2<- data.frame(df1,R0=apply(df1,1,function(params_in)R0(params=params_in)), Fi= apply(df1,1,function(params_in)F_i(state=state_dfe,params=params_in))) which(df2[]<0)
## hacked (not very robustly) to chain label_both() and label_parsed() ... label_special <- function (labels, multi_line=FALSE, sep = "== ") { value <- ggplot2:::label_value(labels, multi_line = multi_line) variable <- ggplot2:::label_variable(labels, multi_line = multi_line) variable[[1]] <- gsub("_(.)$","[\\1]",variable[[1]]) ## not using multiple faceting variables on each margin, so forget paste out <- Map(paste, variable, value, sep = sep) out[[1]] <- lapply(out[[1]], function(expr) c(parse(text=expr))) out }
p1 <- ggplot(df2,aes(x=omega,y=t,z=R0))+ theme_bw() + theme(panel.spacing=grid::unit(0,"lines")) (p1 + geom_contour_filled(breaks=pretty(unique(df2$R0),5), alpha=0.5) ## breaks=0.2) + geom_contour(breaks=1,alpha=0.5,colour="black") + facet_grid(eta_w~eta_t, labeller=label_special) + scale_x_continuous(expand=expansion(c(0,0)), n.breaks=3) + scale_y_continuous(expand=expansion(c(0,0)), n.breaks=3) + scale_fill_viridis_d(name=parse(text="R[0]")) )
Let call the $F_i$ the critical testing rate in the above plot.
Propositions
- $\partial{R_0}/\partial{\eta_t} \geq 0$ and $\partial{R_0}/\partial{\eta_w} \geq 0$.
- $\partial{R_0}/\partial{t} \leq 0$; this is concluded from the Taylor expansion of $R_0$ around $t=0$ which is $R_0 \approx \beta/\gamma + \frac{\beta t}{\omega (\omega+\gamma) \gamma^2 W_s} \Big(\gamma(\eta_w-1)(\gamma W_s+\omega W_i) + (\eta_t -1)P_iW_i \omega^2 \Big) + \mathcal{O}(t^2)$.
- Based on the above Taylor expansion around $t=0$, $$\partial{R_0}/\partial{\omega}= \frac{-\beta t}{\gamma W_s\omega^2 (\gamma+\omega)^2} (a \omega^2 + b \omega + c), $$ where $a=(\eta_w-1)W_i-(\eta_t-1)P_iW_i$, $b=2(\eta_w-1)\gamma W_s$ and $c=(\eta_w-1)\gamma^2 W_s$. Given that $\eta_t\leq \eta_w$, it is inferred that $a\geq 0$, $b\leq 0$ and $c \leq 0$. Let $\omega_1<0$ and $\omega_2>0$ be the roots of the quadratic expression in $\partial{R_0}/\partial{\omega}$. Thus, $\partial{R_0}/\partial{\omega}<0$ for $0<\omega<\omega_2$ and $\partial{R_0}/\partial{\omega}>0$ for $\omega>\omega_2$.
Interpretation of the R0 contour plot
Note that $\eta_t\leq \eta_w$, i.e., reported people have a lower transmission probability than awaiting individuals. The perfect isolation/quarantine occurs when $\eta_w=\eta_t=0$ (see the upper left corner panel). Thus, we only consider the lower-triangle panels in the plot [not sure how to extract them in ggplot?]. (1) Looking at the first column of panels, there seems to be a compensatory relation between $t$ and $\omega$. i.e., higher $t$ requires smaller $\omega$ to keep R0 lower than 1. Also, for a fixed $\omega$, increasing $\eta_w$ requires higher $t$ to achive the critical R0.
Literature Review
\cite{bergstrom2020frequency} (1) Model, assumptions: They developed a function, namely expected expoisour $E(C,\tau)$, to approximate trade-offs between the frequency of testing, n, the sensitivity of testing, q, and the delay between testing and results, d. This function is explicitly derived and was connected the effective reproduction number $R=R_0 S$, where $S$ is the proportion of population susciptable. assumption that transmission rates are a step function: individuals who have COVID go from non-infectious to fully infectious instantaneously, and remain fully infectious until they are no longer able to transmit disease. Test sensitivity takes the same form over the course of infection. More sophisticated models could allow varying infectiousness and varying sensitivity over time, as in \citep{larremore2020test}.
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.