27-10-2026                       >eR-BioStat and STAEM-3D

---

Modeling infectious diseases using R - the SIR model and serological data

Ziv Shkedy and Rashider Aloni based on Chapter 2-4 & 7-9 in the book Modeling infectious diseases parameters based on serological and social contact data (https://www.simid.be/modeling-infectious-disease-parameters-based-on-serological-and-social-contact-data/)

## load/install libraries
.libPaths(c("./Rpackages",.libPaths()))
library(knitr)
library(tidyverse)
library(deSolve)
library(minpack.lm)
library(ggpubr)
library(readxl)
library(gamlss)
library(data.table)
library(grid)
library(png)
library(nlme)
library(gridExtra)
library(mvtnorm)
library(e1071)
library(lattice)
library(ggplot2)
library(dslabs)
library(NHANES)
library(plyr)
library(dplyr)
library(nasaweather)
library(ggplot2)
library(gganimate)
library(av)
library(gifski)
library(foreach)
library("DAAG")
library(DT)
library(TeachingDemos)
library(gridExtra)

1. The SIR model

The starting point in this tutorial is the SIR model in endemic equilibrium. In particular, we assume that all parameters are constant over time and the focus goes to the age of the individual as the time scale of primary interest.

Formulation of the SIR model

Extending the SIR model in terms of changes over age and time yields the following system of partial differential equations: \[ \begin{eqnarray} \frac{\partial S(a,t)}{\partial a}+\frac{\partial S(a,t)}{\partial t}&=&N \mu(a)-(\lambda(a,t)+\mu(a))S(a,t),\nonumber\\ \frac{\partial I(a,t)}{\partial a}+\frac{\partial I(a,t)}{\partial t}&=&\lambda(a,t)S(a,t)-(\nu+\mu(a))I(a,t),\nonumber\\ \frac{\partial R(a,t)}{\partial a}+\frac{\partial R(a,t)}{\partial t}&=&\nu I(a,t)-\mu(a)R(a,t), \end{eqnarray} \]

Here, \(S(a,t)\) is the number of individual in the population that are susceptible at age \(a\) and time \(t\) and \(I(a,t)\) and \(R(a.t)\) are the number of individual in the population that are infected and immune at age \(a\) and time \(t\), respectivly.

Parameters

• N: population size.
• The force of infection: \(\lambda(a,t)\).
• The average age at infection, A. For a constant force of infection, \[ A=\frac{1}{\lambda}. \]
• Death rate and birth rate: \(\mu(a)\).
• Recovery rate: \(\nu\).
• The average duration in the infected class: \(\frac{1}{\nu}\).

Endemic equilibrium

Under endemic equilibrium, \[ \frac{\partial S(a,t)}{\partial t}=\frac{\partial I(a,t)}{\partial t}=\frac{\partial R(a,t)}{\partial t}=0, \]

the SIR model is simplified

\[ \begin{eqnarray} \frac{d S(a)}{d a}&=&N \mu(a)-(\lambda(a)+\mu(a))S(a),\nonumber\\ \frac{d I(a)}{d a}&=&\lambda(a)S(a)-(\nu+\mu(a))I(a),\nonumber\\ \frac{d R(a)}{d a}&=&\nu I(a)-\mu(a)R(a), \end{eqnarray} \]

2. The SIR model in a closed population

We consider a closed population, i.e., individuals do not enter or exit to/from the population, that is \(\mu(a)=0\). In this case the SIR model is given by

\[ \begin{eqnarray} \frac{d S(a)}{d a}&=&-\lambda(a)S(a),\nonumber\\ \frac{d I(a)}{d a}&=&\lambda(a)S(a)-\nu I(a),\nonumber\\ \frac{d R(a)}{d a}&=&\nu I(a). \end{eqnarray} \]

There are two parameters in the model, the force of infection \(\lambda\) and the recovery rate \(\nu\).

Model assumptions

• No mortality access due to the disease.
• Life long immunity.

Example 1 (\(\lambda=0.2\),D=10 days)

library(deSolve)

SIR model with \(\lambda =0.2\) and \(D=10\) days (\(\nu=36.5\))

For a population in which the average age at infection is 5 years (i.e., the average time that is spent in the susceptible class is 5 years) we have \[ \lambda=\frac{1}{5}=0.2. \]

For a recovery rate of 10 days,

\[ \nu= \left(\frac{10}{365} \right )^{-1}=36.5 \]

In R, the object parameters is used to store the parameter values.

parameters <- c(lambda = 0.2, v=36.5)
parameters

## lambda      v 
##    0.2   36.5

We assume a population size of \(N=5000\) with one infected individual at age=0, \(I(0)=1\). This implies that \(S(0)=N-1=4999\) and \(R(0)=0\). The state vector is given by

state <- c(X=4999,Y=1,Z=0)
state

##    X    Y    Z 
## 4999    1    0

Note that the R objects X, Y and Z correspond to the compartments S, I and R, respectively.

Model formulation

The SIR model

\[ \begin{eqnarray} \frac{d S(a)}{d a}&=&-\lambda(a)S(a),\nonumber\\ \frac{d I(a)}{d a}&=&\lambda(a)S(a)-\nu I(a),\nonumber\\ \frac{d R(a)}{d a}&=&\nu I(a), \end{eqnarray} \]

can be defined in R using theSIR function in the following way

SIR<-function(t,state,parameters)
{
with(as.list(c(state, parameters)),
{
dX <- -lambda*X 
dY <- lambda*X - v*Y
dZ <- v*Y
list(c(dX, dY, dZ))
}) 
}

The object times is the ages that we use to integrate the equation system in order to produce the output. In our example, the age ranges from 0 to 40 years old.

times<-seq(0,40,by=0.01)
times[1:20]

##  [1] 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13 0.14
## [16] 0.15 0.16 0.17 0.18 0.19

Runing the model

We solve the equations system and produce the model’s solution using the function >tt>ode of the R package deSolve. In the panel below, \[ S(0)=4999, I(0)=1,R(0)=0, \] and \[ S(0.05)=4949.259, I(0.05)=22.989501,R(0.05)=27.751380. \]

Numerical solution

The R object out contain the numerical solution of the model. The objects X, Y and Z are the number of susceptible, infected and immune individuals at age a (the R object times).

require(deSolve)
out <- as.data.frame(ode(y=state,times=times,func=SIR,parms=parameters))
head(out)

##   time        X         Y         Z
## 1 0.00 4999.000  1.000000  0.000000
## 2 0.01 4989.012  9.061818  1.926190
## 3 0.02 4979.044 14.641580  6.314481
## 4 0.03 4969.096 18.498345 12.405853
## 5 0.04 4959.168 21.159066 19.673391
## 6 0.05 4949.259 22.989501 27.751380

Graphical solution

Figure 1 shows the numerical solution for the example. We notice that, the fraction infected individuals is relatively small compared to the fraction susceptible and immune. This is due to the different time spent in each compartment. The average duration in the susceptible class is 5 year and after recovery individuals gain life long immunity against reinfection.

par(mfrow=c(1,2), oma=c(0,0,3,0))
plot (times,out$X ,type="l",main="S and R", xlab="age", ylab="-",lwd=2)
lines(times,out$Z,col=3,lwd=2)
legend(20,4000,c("S","R"),lty=c(1,1),col=c(1,3))
plot (times,out$Y ,type="l",main="I", xlab="time", ylab="-",lwd=2)
mtext(outer=TRUE,side=3,"SIR model, D=10 days",cex=1.5)

Figure 1: Solution for a SIR model with lambda=0,2 and D=10 days.

Example 2 (\(\lambda=0.2\),D=60 days)

SIR with \(\lambda =0.2\) D=two months (\(\nu=6.08\))

In this example we change the recovery rate from 10 days to two months (60 days) so individuals stay in the infected class a longer time compare to the first example. In this case \[ \nu= \left(\frac{60}{365} \right )^{-1}=6.08 \]

60/365

## [1] 0.1643836

1/(60/365)

## [1] 6.083333

The new parameters vector is given by

parameters <- c(lambda = 0.2, v=6.083333)
parameters

##   lambda        v 
## 0.200000 6.083333

state <- c(X=4999,Y=1,Z=0)
#state

Runing the model

Note that compared to an average duration of 10 days in the infected class (Example 1), we expect to observed a high proportion of infected individuals at any age which can be clearly seen in Figure 2.

require(deSolve)
out <- as.data.frame(ode(y=state,times=times,func=SIR,parms=parameters))
head(out)

##   time        X        Y         Z
## 1 0.00 4999.000  1.00000 0.0000000
## 2 0.01 4989.012 10.63116 0.3568532
## 3 0.02 4979.044 19.67452 1.2815405
## 4 0.03 4969.096 28.16483 2.7393722
## 5 0.04 4959.168 36.13475 4.6977098
## 6 0.05 4949.259 43.61504 7.1258399

Graphical solution

par(mfrow=c(1,2), oma=c(0,0,3,0))
plot (times,out$X ,type="l",main="S and R", xlab="age", ylab="-",lwd=2)
lines(times,out$Z,col=3,lwd=2)
legend(20,4000,c("S","R"),lty=c(1,1),col=c(1,3))
plot (times,out$Y ,type="l",main="Y", xlab="time", ylab="-",lwd=2)
mtext(outer=TRUE,side=3,"SIR model D=two months", cex=1.5)

Figure 2: Solutions for a SIR model D=60 days.

Example 3: age dependent force of infection

In this example we change the formulation of the force of infection and assume that \(\lambda(a)=\beta \times I(a)\). Here \(\beta\) is the transmission rate and \(I(a)\) is the number of infected individuals at age a.

\(\beta=0.0085\) and \(\lambda=\beta \times I(a)\)

For a recovery rate of 10 days, the new parameters vector is given by

parameters <- c(beta=0.0085, v=36.5)
state <- c(X=4999,Y=1,Z=0)
parameters

##    beta       v 
##  0.0085 36.5000

#state
times<-seq(0,10,by=0.01)
#times[1:20]

Definition of the model

For the SIR model, we use beta Y X instead of lambdaYX

SIR<-function(t,state,parameters)
{
with(as.list(c(state, parameters)),
{
dX <- -beta*Y*X 
dY <- beta*Y*X - v*Y
dZ <- v*Y
list(c(dX, dY, dZ))
}) 
}

Runing the model and graphical output

The solution of the model is presented in Figure 3. Note that the force of infection in this model is not constant but it is proportional for the number of individuals in the infected compartment (Y) presented in Figure 3.

require(deSolve)
out <- as.data.frame(ode(y=state,times=times,func=SIR,parms=parameters))
head(out)

##   time        X        Y         Z
## 1 0.00 4999.000 1.000000 0.0000000
## 2 0.01 4998.562 1.061727 0.3761535
## 3 0.02 4998.097 1.127220 0.7755187
## 4 0.03 4997.604 1.196705 1.1995104
## 5 0.04 4997.080 1.270419 1.6496279
## 6 0.05 4996.524 1.348611 2.1274607

par(mfrow=c(2,2), oma=c(0,0,3,0))
plot (times,out$X ,type="l",main="X", xlab="time", ylab="-")
plot (times,out$Y ,type="l",main="Y", xlab="time", ylab="-")
plot (times,out$Z ,type="l",main="Y", xlab="time", ylab="-")
mtext(outer=TRUE,side=3,"SIR model, lambda=beta*I",cex=1.5)

Figure 3: Solution for a SIR model with an age dependent force of infection.

3. Modeling serological data

Current status data

Consider an age-specific cross-sectional prevalence sample of size \(N\) and let \(a_{i}\) be the age and \(Z_i\) the antibody activity level (in U/ml) of the \(i\)th subject. The antibody activity levels \(Z_i\) are quantitative results calculated from the optical densities obtained from commercial IgG enzyme immune-assays (EIAs).

We define the dichotomized version of \(Z_i\) as the binary variable \(Y_{i}\) \[\begin{equation} \label{binom} Y_{i}= \left\{ \begin{array}{ll} 1 & \mbox{if $Z_i > \tau_u$}, \\ 0 & \mbox{if $Z_i < \tau_\ell$}, \\ \end{array} \right. \end{equation}\] and missing if \(\tau_\ell < Z_i < \tau_u\) where \(\tau_\ell\) and \(\tau_u\) are the lower and upper threshold values (20 resp. 24 U/ml for Parvo B19). So, for each individual \(i\) we observe whether s/he has experienced infection before age \(a_i\). We do not observe the age at infection, but rather a censored value: right censored if \(Y_i = 0\) and left censored if \(Y_i=1\). This censored data point of view and this link to survival analysis will regularly show up.

The prevalence

Denote \(\pi(a_{i})\) the probability to be infected before age \(a_{i}\). This probability will be referred to as the so-called seroprevalence. It is important to note that, due to misclassification, the seroprevalence is not exactly equal to the disease prevalence.

Methods to model \(\pi(a_{i})\) are based on maximizing the log likelihood as a function of the model parameters, based on the data \(\{(a_i,Y_i)\}_{i=1}^N\), where \(N\) denotes the sample size (excluding the equivocal cases) and is given by

\[ \begin{equation} L=\sum_{i=1}^{N}Y_{i}\log \left \{\pi(a_{i}) \right \} +(1-Y_{i}) \log \left \{1-\pi(a_{i}) \right \}. \end{equation} \]

GLM for the prevalence

In the terminology of generalized linear models, the age-dependent probability \(\pi(a)\) is modeled as

\[ \begin{equation} \pi(a) = g^{-1}(\eta(a)), \end{equation} \] where \(\eta(a)\) is the so-called linear predictor and \(g\) is the so-called link function. For binary responses, is often taken to be a logit link function, \(\mbox{log}(\pi/(1-\pi))\), but other link functions such as the complementary log-log link, \(\log(-\log(1-\pi)),\) and log link, \(-\log(1-\pi)\), can be used as well. In the next section the age-dependent force of infection \(\lambda(a)\) is derived from the seroprevalence \(\pi(a)\), for different link functions \(g\).

Age-Dependent Force of Infection

The force of infection can be easily derived from the seroprevalence under assumptions including time homogeneity and life-long immunity. In terminology of survival analysis, the FOI is nothing else than the hazard function. Using a so-called catalytic model with log link (i.e. \(\pi(a)=1- e^{-\eta \left( a \right)}\)) leads to a simple interpretation of the first derivative of the linear predictor. Indeed, \(\eta(a)\) is the nothing else than the cumulative hazard and hence the force of infection is simply the first derivative of the linear predictor.

\[ \begin{equation} \lambda(a)=\frac{\pi'(a)}{1-\pi(a)}= \frac{ \eta'(a) e^{-\eta(a)}}{e^{-\eta(a)}}=\eta'(a). \end{equation} \] In general, when the link function is not restricted to be the log link, the force of infection can still be derived analogously. The force of infection can be expressed as a product of two functions: \[\begin{equation} \lambda(a)=\eta'(a)\delta(\eta(a)). \end{equation}\] Here, \(\delta(\cdot)\) is a known function for which the form is determined by the link function \(g\).

4. Serological samples

For illustration of all models we discussed below we use a serological data from UK (Rubella) and Bulgaria (Hepatitis A).

Rubela in UK

robula<- list(
AGE = c(1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 
    10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5,
    22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5,
    34.5, 35.5, 36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5),
POS = c(31, 30, 34, 57, 95, 104, 90, 96, 134, 110, 111, 147, 138, 141, 53,
    49, 73, 69, 97, 65, 74, 84, 82, 79, 90, 84, 81, 72, 71, 51, 45, 45, 35,
    39, 36, 37, 37, 37, 28, 26, 25, 21, 18, 18), 
NEG = c(175, 116, 134, 
    132, 124, 91, 74, 49, 46, 50, 37, 31, 39, 24, 14, 9, 8, 10, 14, 11, 8, 
    17, 6, 6, 4, 7, 8, 4, 8, 5, 7, 3, 2, 2, 4, 1, 2, 4, 2, 1, 0, 1, 1, 0))
Rub1<-data.frame(robula$AGE,robula$POS,robula$NEG,robula$POS+robula$NEG)
names(Rub1)<-c("AGE","POS","NEG","NTOT")
head(Rub1)

##   AGE POS NEG NTOT
## 1 1.5  31 175  206
## 2 2.5  30 116  146
## 3 3.5  34 134  168
## 4 4.5  57 132  189
## 5 5.5  95 124  219
## 6 6.5 104  91  195

plot(Rub1$AGE,Rub1$POS/Rub1$NTOT)

Figure 4: Rubela in UK.

Hepatatis A in Bulgria

kei1<-list(
age=c(
   1, 2, 3, 4, 5, 6,7, 8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22,
  23,24,25,26,27,28,29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44,
  45,46,47,48,49,50,51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66,
  67,68,69,70,72,73,74, 75 ,76, 77, 78, 79, 80, 81, 82, 83, 86)+0.5,
posi=c(
  3, 3, 3, 4, 7, 4, 3, 4, 7, 8, 2, 3, 2, 0, 5,13, 1, 3,15,22,15, 7,
  8, 7,12, 5,10,15, 9, 9, 9, 8, 9, 8, 9,13, 6,15,11, 6, 8,13, 7, 5,
  7, 9, 9,22, 6,10, 6,13, 8, 7,13,11, 8, 8, 9,13, 5, 5, 5, 5,10, 8,
  4, 5, 4, 8, 9, 1, 4, 7, 6, 2, 3, 2, 4, 1, 1, 2, 1),
ni=c(
  16,15,16,13,12,15,12,11,10,15, 7, 7,11, 1,16,41, 2, 6,32,37,24,10,
  10,11,15,10,13,19,12, 9,14,10,11, 9,14,14, 7,16,13, 8, 8,14,10, 5,
   7, 9, 9,22, 7,10, 6,14, 8, 7,13,11, 8, 8,10,16, 5, 6, 5, 5,10, 8,
   4, 5, 5, 8, 9, 1, 4, 7, 6, 2, 3, 2, 4, 1, 1, 2, 1))
Kei1<-data.frame(kei1$age,kei1$posi,kei1$ni-kei1$posi,kei1$ni)
names(Kei1)<-c("AGE","POS","NEG","NTOT")
head(Kei1)

##   AGE POS NEG NTOT
## 1 1.5   3  13   16
## 2 2.5   3  12   15
## 3 3.5   3  13   16
## 4 4.5   4   9   13
## 5 5.5   7   5   12
## 6 6.5   4  11   15

plot(Kei1$AGE,Kei1$POS/Kei1$NTOT)

Figure 5: Hepatitis A in Bulgaria.

5. Modeling serological data: early models

Munch’s model

The first references on the estimation of the force of infection goes back to the early work of Muench, followed by Griffiths. Muench suggested to model the infection process with a so-called catalytic model, in which the distribution of the time spent in the susceptible class in the SIR model is exponential with rate \(\beta\). More precisely,Muench proposed to model the prevalence by \[ \pi(a)=k\left \{ 1-\exp(-\beta a)\right \}, \] where \(1-k\) is the proportion of the population staying uninfected for lifetime. Under this catalytic model and assuming that \(k=1\), it follows that \(\lambda(a) = \beta\). So in this case the force of infection is age independent.

Application to Hepatitis A in Bulgria

in R, Munch’s model can be fitted using the glm() function as a GLM with \(log\) link function and linear predictor \(\eta(a)=a\) in the following way glm(sero negative/ total ~ -1+a, family = binomial(link=log)). The parameter estimate for the (constant) force of infection is equal to \(\lambda=0.05185577\).

fit.munch<-glm(NEG/NTOT ~ -1+AGE, family = binomial(link=log), data = Kei1)
summary(fit.munch)

## 
## Call:
## glm(formula = NEG/NTOT ~ -1 + AGE, family = binomial(link = log), 
##     data = Kei1)
## 
## Coefficients:
##      Estimate Std. Error z value Pr(>|z|)    
## AGE -0.051856   0.008321  -6.232  4.6e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance:    Inf  on 83  degrees of freedom
## Residual deviance: 10.104  on 82  degrees of freedom
## AIC: 39.321
## 
## Number of Fisher Scoring iterations: 6

Figure 6 shows the data and the estimated prevalence and force of infection.

par(mfrow=c(1,2))
plot(Kei1$AGE,1-Kei1$NEG/Kei1$NTOT,xlab="AGE",ylab="Prevalence")
lines(Kei1$AGE,1-fit.munch$fit)
plot(Kei1$AGE,rep(-fit.munch$coeff,length(Kei1$AGE)),type="l",xlab="AGE",ylab="force of infection")

Figure 6: Hepatitis in Bulgria. Right panel: prevalence. Left Panel: the force of infection

Griffiths’ model

Griffiths proposed a model for measles in which the force of infection increases linearly in the age range \((\tau,10)\), \(\tau \ge 0\) \[\lambda(a) = \left\{ \begin{array}{ll} \beta_ {1}(a+\beta_{0}) & a > \tau, \\ 0 & a \le \tau. \end{array} \right. \]

He focused on the interval \((\tau,10)\) as the interval \((0,\tau)\) represents the period of inherited immunity and since more than 95% of the measles cases occur before the age of 10. Note that, since specified \(\tau\) as a parameter in the model, Griffiths’ model should be interpreted as a change-point model. Interestingly, justified his choice of a linear force of infection by using a non-parametric estimate for the force of infection, \(\lambda(a)=\Delta \pi/(1-\pi(a))\), which he plotted against age and showed the linear trend of the force of infection. In coming chapters, we will discuss more advanced non-parametric models for the force of infections. himself mentioned that his model for the prevalence corresponds to a model in which the linear predictor is a quadratic function of age.

Griffith’s model with linear force of infection \(\lambda(a)=\beta_{1}+2 \beta_{2} a\) can be estimated using a GLM for which the linear predictor is a quadratic function of the age, i.e. \(\eta(a)=\beta_{1} a+\beta_{2}a^{2}\). In the code below, the variables and correspond to \(a\) and \(a^{2}\), respectively.

Application to Hepatitis A in bulgria

Grifiths’ model can be ftted suing the glm() function as a GLM with log link function in the following way: glm(sero negative/ total ~-1+a+a^2, family = binomial(link=log)). Note that the linear predictor in this case is \(\eta(a)=a+a^2\). The panel below shows the parameter estimates for the model

age1<-Kei1$AGE
age2<-Kei1$AGE^2
age3<-Kei1$AGE^3
fit.griffiths<-glm(Kei1$NEG/Kei1$NTOT~-1+age1+age2, family = binomial(link=log))
summary(fit.griffiths)

## 
## Call:
## glm(formula = Kei1$NEG/Kei1$NTOT ~ -1 + age1 + age2, family = binomial(link = log))
## 
## Coefficients:
##        Estimate Std. Error z value Pr(>|z|)  
## age1 -0.0415572  0.0183284  -2.267   0.0234 *
## age2 -0.0002707  0.0004581  -0.591   0.5545  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance:   Inf  on 83  degrees of freedom
## Residual deviance: 9.704  on 81  degrees of freedom
## AIC: 40.411
## 
## Number of Fisher Scoring iterations: 7

Figure 7 shows the data with estimated prevalence and the estimated force of infection. The linear force of infection, \(\lambda(a)=\beta_{1}+2 \beta_{2} a=0.0415572+2 \times 0.0002707\) as shown in the right panel.

foi.linear<- -fit.griffiths$coeff[1]+ -fit.griffiths$coeff[2]*Kei1$AGE
par(mfrow=c(1,2))
plot(Kei1$AGE,1-Kei1$NEG/Kei1$NTOT,xlab="AGE",ylab="Prevalence")
lines(Kei1$AGE,1-fit.griffiths$fit)
plot(Kei1$AGE,foi.linear,type="l",xlab="AGE",ylab="force of infection")

Figure 7: Hepatitis in Bulgria. Right panel: prevalence. Left Panel: the force of infection

Grenfell and Anderson’s model

Grenfell extended the models of Muench and Griffiths further and used polynomial functions to model the force of infection. The advantage of higher order polynomials is their flexible curve shapes. Grenfell did not restrict a priori the force of infection to be constant or linear but let the data determine the order of the polynomial. Their model assumes that \(\pi(a)=1-e^{-\Sigma \beta_{i} a^{i}}\) which implies that, again using a GLM with log link., the force of infection equals \(\lambda(a)=\sum \beta_{i} i a^{i-1}\). R code for a possible model with quadratic force of infection such as \(\lambda(a)=\eta`(a)=\beta_{1}+2 \beta_{2} a + 3 \beta_{3} a^{2}\) is given by

Application to Hepatitis A in Bulgaria

We use glm(sero negative/total ~ a + a^2 + a^3, family = binomial(link=log)) to fit Grenfel and Anderson model as third order polynomial GLM with log link function and \(\eta(a)=\beta_{0}+\beta_{1}a+ \beta_{2} a^2 + \beta_{3} a^{3}\)

fit.grenfell<-glm(Kei1$NEG/Kei1$NTOT ~ age1+age2+age3, family = binomial(link=log))
summary(fit.grenfell)

## 
## Call:
## glm(formula = Kei1$NEG/Kei1$NTOT ~ age1 + age2 + age3, family = binomial(link = log))
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.314e-01  4.063e-01  -0.816    0.415
## age1         9.012e-03  6.692e-02   0.135    0.893
## age2        -1.903e-03  2.680e-03  -0.710    0.478
## age3         1.364e-05  2.762e-05   0.494    0.622
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38.9554  on 82  degrees of freedom
## Residual deviance:  8.7659  on 79  degrees of freedom
## AIC: 45.755
## 
## Number of Fisher Scoring iterations: 8

The age dependent force of infection, \(\lambda(a)=\eta`(a)=\beta_{1}+2 \beta_{2} a + 3 \beta_{3} a^{2}\), is shown in Figure 8 (right panel).

foi.flex<- -fit.grenfell$coeff[1]+ -fit.grenfell$coeff[2]*Kei1$AGE+ -fit.grenfell$coeff[3]*Kei1$AGE^2
par(mfrow=c(1,2))
plot(Kei1$AGE,1-Kei1$NEG/Kei1$NTOT,xlab="AGE",ylab="Prevalence")
lines(Kei1$AGE,1-fit.grenfell$fit)
plot(Kei1$AGE,foi.flex,type="l",xlab="AGE",ylab="force of infection")

Figure 8: Hepatitis in Bulgria. Right panel: prevalence. Left Panel: the force of infection

6. Modeling the focrce of infection using GLMs

In Section 5, all fitted models were formulated as GLMs with \(log\) link functions. In this section we illustrate how to estimate the force of infection with different link functions.

Age dependent focrce of infection: Weibul models

Other parametric models, fitted within the framework of generalized linear models (GLM) with binomial error, were discussed by Becker1989, Diamond and Keiding, who used models with link function in order to parameterize the prevalence and the force of infection as a Weibull model. For a Weibull model, the prevalence is given by \[ \pi(a)=1-e^{-\beta_{0} a^{\beta_{1}}}. \] The linear predictor for a model with complementary log-log link function is given by \[ \eta(a)=\log(\beta_{0})+\beta_{1} \log(a). \] The Weibull model implies that the force of infection is a monotone function of the exposure time, \[ \lambda(a)=\beta_{0}\beta_{1} a ^{\beta_{1}-1}. \] In case that \(\beta_{1} > 1\), the force of infection is monotone increasing, while \(\beta_{1}<1\) implies that the force of infection decreases with the exposure time.

Rubela in UK

fit.weibul<-glm(POS/NTOT ~ log(AGE), family = binomial(link=cloglog),data = Rub1)
summary(fit.weibul)

## 
## Call:
## glm(formula = POS/NTOT ~ log(AGE), family = binomial(link = cloglog), 
##     data = Rub1)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -1.8924     0.9236  -2.049  0.04046 * 
## log(AGE)      0.8498     0.3082   2.758  0.00582 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 12.0671  on 43  degrees of freedom
## Residual deviance:  1.0133  on 42  degrees of freedom
## AIC: 21.276
## 
## Number of Fisher Scoring iterations: 6

fit.weibul$coeff

## (Intercept)    log(AGE) 
##  -1.8923888   0.8497694

alpha<-exp(fit.weibul$coeff[1])
beta<-fit.weibul$coeff[2]
lambda2<-alpha*beta*Rub1$AGE^(beta-1)
par(mfrow=c(1,2))
plot(Rub1$AGE,Rub1$POS/Rub1$NTOT,xlab="AGE",ylab="Prevalence")
lines(Rub1$AGE,fit.weibul$fit,xlab="AGE",ylab="force of infection")
plot(Rub1$AGE,lambda2,type="l")

Figure 9: Rubela in UK. Right panel: prevalence. Left Panel: the force of infection

Constant force of infection

A model with constant force of infection can be estimated as a GLM with \(log\) link (as in Section 5) or as a GLM with \(c-log-log\) link for which \(\beta_{1}=1\),

\[ \eta(a)=\log(\beta_{0})+ \log(a). \] The Weibull model with \(\beta_{1}=1\) implies that the force of infection is constant and equal to, \[ \lambda(a)=\beta_{0}\beta_{1} a ^{\beta_{1}-1}=\beta_{0}. \]

`

Rubela in UK

In %, we can fit a GLM with \(c-log-log\) link using the flowing code for the glm function: glm(sero positive/ total ~ offset(log(a)), family = binomial(link=cloglog)). Note that we use log(age) as an offset predictor (so \(\beta_{1}=1\)). For this model \(\hat{\beta_{0}=-2.336\) so the estimated force of infection is equal to

exp(-2.336)

## [1] 0.09671372

The panel below presents the estimated Weibul model with \(\beta_{1}=1\).

fit.exp<-glm(POS/NTOT ~ offset(log(AGE)), family = binomial(link=cloglog), data = Rub1)
summary(fit.exp)

## 
## Call:
## glm(formula = POS/NTOT ~ offset(log(AGE)), family = binomial(link = cloglog), 
##     data = Rub1)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -2.336      0.219  -10.67   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1.2493  on 43  degrees of freedom
## Residual deviance: 1.2493  on 43  degrees of freedom
## AIC: 18.914
## 
## Number of Fisher Scoring iterations: 5

For comparison, exponential model with constant force of infection is shown below and in Figure 10. Note that the two models are equivalent.

fit.exp<-glm(NEG/NTOT ~ -1+AGE, family = binomial(link=log), data = Rub1)
summary(fit.exp)

## 
## Call:
## glm(formula = NEG/NTOT ~ -1 + AGE, family = binomial(link = log), 
##     data = Rub1)
## 
## Coefficients:
##     Estimate Std. Error z value Pr(>|z|)    
## AGE -0.09668    0.02117  -4.566 4.97e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance:    Inf  on 44  degrees of freedom
## Residual deviance: 1.2493  on 43  degrees of freedom
## AIC: 18.914
## 
## Number of Fisher Scoring iterations: 6

lambda.exp<-fit.exp$coeff
par(mfrow=c(1,2))
plot(Rub1$AGE,1-Rub1$NEG/Rub1$NTOT,xlab="AGE",ylab="Prevalence")
lines(Rub1$AGE,1-fit.exp$fit)
fit.exp$coeff

##         AGE 
## -0.09668382

plot(Rub1$AGE,rep(-fit.exp$coeff,length(Rub1$AGE)),type="l",xlab="AGE",ylab="force of infection")

Figure 10: Rubela in UK. Right panel: prevalence. Left Panel: the force of infection

Age dependent force of infection: log-logistic models

The log-logistic model can be fitted as a GLM with \(logit\) link with log(age) as predictor, that is \[ \eta(a)=\beta_{0}+\beta_{1}log(a). \]

For this model the force of infection is given by

\[ \lambda(a)=\frac{\beta_{0}a^{\beta_{0}-1}\beta_{1}}{1+\beta_{1}a^{\beta_{0}}}. \]

Rubela in UK

in R we use glm(sero positive/ total ~ log(a), family = binomial(link=logit)).

fit.loglogistic<-glm(POS/NTOT ~ log(AGE), family = binomial(link=logit),data = Rub1)
summary(fit.loglogistic)

## 
## Call:
## glm(formula = POS/NTOT ~ log(AGE), family = binomial(link = logit), 
##     data = Rub1)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -2.9380     1.5088  -1.947  0.05151 . 
## log(AGE)      1.6284     0.5725   2.844  0.00445 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 12.06714  on 43  degrees of freedom
## Residual deviance:  0.67183  on 42  degrees of freedom
## AIC: 19.459
## 
## Number of Fisher Scoring iterations: 5

Figure 11 shows the prevalence and force of infection.

beta<-fit.loglogistic$coeff[1]
alpha<-fit.loglogistic$coeff[2]
aa<-alpha*(Rub1$AGE^(alpha-1))*beta
bb<-1+beta*(Rub1$AGE^(alpha))
lambda3<-aa/bb
par(mfrow=c(1,2))
plot(Rub1$AGE,Rub1$POS/Rub1$NTOT,xlab="AGE",ylab="Prevalence")
lines(Rub1$AGE,fit.loglogistic$fit)
plot(Rub1$AGE,lambda3,type="l",xlab="AGE",ylab="force of infection")

Figure 11: Rubela in UK. Right panel: prevalence. Left Panel: the force of infection

Constant and age-dependent focrce of infection

In this tutorials, we discuss different models to estimated both the prevalence and force of infection. All models were fitted within the GLM modeling framework and were used, as shown in Figure 12, to estimated constant, monotone and flexible models for the force of infection.

lambda<-0.05
agei<-seq(from=0,to=40,length=1000)
lambda2<-0.5*1*agei^(-0.5)
par(mfrow=c(2,2))
pro1<-1-exp(-lambda*agei)
plot(agei,rep(lambda,length(agei)),xlab="age",ylab="foi",type="l")
lines(agei,rep(lambda,length(agei)),lty=1,col=2,lwd=3)
title("Exponential")
plot(agei,lambda2,xlab="age",ylab="foi",type="l")
lines(agei,lambda2,lty=1,col=2,lwd=3)
title("Weibull")
alpha<- 1.5
beta<- 0.25
aa<-alpha*(agei^(alpha-1))*beta
bb<-1+beta*(agei^(alpha))
lambda3<-aa/bb
plot(agei,lambda3,xlab="age",ylab="foi",type="l")
lines(agei,lambda3,lty=1,col=2,lwd=3)
title("Log-logistics")

Figure 12: Rubela in UK. Right panel: prevalence. Left Panel: the force of infection

7. The SIR model and prevalence estimation

In this section we explore the connection between a solution obtained for a SIR model the estimated prevalence obtained for a model with a constant force of infection. We consider a SIR model for a closed population of the form

\[ \begin{eqnarray} \frac{d S(a)}{d a}&=&-\lambda(a)S(a),\nonumber\\ \frac{d I(a)}{d a}&=&\lambda(a)S(a)-\nu I(a),\nonumber\\ \frac{d R(a)}{d a}&=&\nu I(a). \end{eqnarray} \]

with a constant force of infection, \(\lambda(a)=0.05\) (i.e., \(A=20\) years) and a recovery rate of 60 days.

N<-5000
parameters <- c(lambda = 0.05, v=6.083333)
parameters

##   lambda        v 
## 0.050000 6.083333

state <- c(X=4999,Y=1,Z=0)
#state
times<-seq(0,85,by=0.01)
#times

SIR<-function(t,state,parameters)
{
with(as.list(c(state, parameters)),
{
dX <- -lambda*X 
dY <- lambda*X - v*Y
dZ <- v*Y
list(c(dX, dY, dZ))
}) 
}

We run the model for 85 years.

require(deSolve)
out <- as.data.frame(ode(y=state,times=times,func=SIR,parms=parameters))
head(out)

##   time        X         Y         Z
## 1 0.00 4999.000  1.000000 0.0000000
## 2 0.01 4996.501  3.365359 0.1335165
## 3 0.02 4994.003  5.589903 0.4065979
## 4 0.03 4991.507  7.681944 0.8109345
## 5 0.04 4989.012  9.649302 1.3387062
## 6 0.05 4986.518 11.499337 1.9825540

Figure 13 shows the model’s solution (in proportions).

par(mfrow=c(1,2), oma=c(0,0,3,0))
plot (times,out$X/N ,type="l",main="S and R", xlab="age", ylab="-",lwd=2
      ,ylim=c(0,1))
lines(times,out$Z/N,col=3,lwd=2)
plot (times,out$Y/N ,type="l",main="I", xlab="time", ylab="-",lwd=2)
mtext(outer=TRUE,side=3,"SIR model D=two months", cex=1.5)

Figure 13: Solutions for the SIR model

Figure 14, shows the proportion of immune in the serological sample of Hepatitis A in Bulgaria, a prevalence estimate for a model with a constant force of infection (red line) and the SIR model solution for the proportion immune (green line). Clearly, the two models predict the same pattern for the prevalence. However, there is a fundamental difference between the two models. For the GLM estimate the parameters were estimated from the data while the solution of a SIR model was obtained from a model in which the parameters were assumed to be known. Furthermore, the SIR model assume a specific transmission process for the diseases while the GLM assume an underlying distribution for the time that susceptible spent in the susceptible compartment.

par(mfrow=c(1,1))
plot(Kei1$AGE,1-Kei1$NEG/Kei1$NTOT,ylim=c(0,1),xlab="age",ylab="prevalence")
lines(times,out$Z/N,lty=1,lwd=2,col=3)
lines(Kei1$AGE,1-fit.munch$fit,lwd=2,col=2)
#lines(times,out$Z/N,col=3,lwd=2)
#plot (times,out$Y/N ,type="l",main="Y", xlab="time", ylab="-",lwd=2)
#mtext(outer=TRUE,side=3,"SIR model D=two months", cex=1.5)
legend(60,0.2,c("GLM","SIR"),lty=c(1,1),col=c(2,3))

Figure 14: Prevalance estimates for hepatitis A in Bulgaria.

8. Non parametric models

The PAVA algorethm

pava <- function (x, wt=rep(1,length(x)))

#  Compute the isotonic regression of numeric vector 'x', with
#  weights 'wt', with respect to simple order.  The pool-adjacent-
#  violators algorithm is used.  Returns a vector of the same length
#  as 'x' containing the regression.
#  02 Sep 1994 / R.F. Raubertas

{

   n <- length(x)

   if (n <= 1) return (x)

   if (any(is.na(x)) || any(is.na(wt))) {

      stop ("Missing values in 'x' or 'wt' not allowed")

   }

   lvlsets <- (1:n)

   repeat {

      viol <- (as.vector(diff(x)) < 0)  # Find adjacent violators

      if (!(any(viol))) break

 

      i <- min( (1:(n-1))[viol])        # Pool first pair of violators

      lvl1 <- lvlsets[i]

      lvl2 <- lvlsets[i+1]

      ilvl <- (lvlsets == lvl1 | lvlsets == lvl2)

      x[ilvl] <- sum(x[ilvl]*wt[ilvl]) / sum(wt[ilvl])

      lvlsets[ilvl] <- lvl1

   }

  return(x)

}

Application to Rubela in UK

Estimating the prevalence

xi<-Rub1$AGE
wi<-Rub1$NTOT
yi<-Rub1$POS/Rub1$NTOT
par(mfrow=c(1,1))
plot(xi,yi,pch="*",xlab="age",ylab="prevalence")
lines(xi,yi,lty=2)
iso1<-pava(yi,wt=wi)
#iso1
#length(unique(iso1))
lines(xi,iso1)

Figure 15: Rubela in UK. Right panel: prevalence. Left Panel: the force of infection

Estimation the force of infection

d.iso1<-diff(iso1)
foi.iso<-d.iso1/(1-iso1[-c(1)])
plot(xi[-c(1)],foi.iso,xlab="age",ylab="force of infection")
foi.iso1<-foi.iso[-c(43)]
fit.foi<-loess(foi.iso1~xi[-c(1,44)],span=0.5)
lines(xi[-c(1,44)],fit.foi$fit)

Figure 16: Rubela in UK. Right panel: prevalence. Left Panel: the force of infection

Hepatitis A in Bulgaria

Estimating the prevalence

xi<-Kei1$AGE
wi<-Kei1$NTOT
yi<-Kei1$POS/Kei1$NTOT
par(mfrow=c(1,1))
plot(xi,yi,pch="*",xlab="age",ylab="prevalence")
lines(xi,yi,lty=2)
iso1<-pava(yi,wt=wi)
#iso1
#length(unique(iso1))
lines(xi,iso1)

Figure 17: Rubela in UK. Right panel: prevalence. Left Panel: the force of infection

Estimating the force of infection

d.iso1<-diff(iso1)
#length(d.iso1)
#length(1-iso1[-c(1)])
foi.iso<-d.iso1/(1-iso1[-c(1)])
plot(xi[-c(1)],foi.iso)

Figure 18: Rubela in UK. Right panel: prevalence. Left Panel: the force of infection