Parametric models

In serosv: Model Infectious Disease Parameters from Serosurveys

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(serosv)

Polynomial models

Refer to Chapter 6.1.1

Use polynomial_model() to fit a polynomial model.

We will use the Hepatitis A data from Belgium 1993--1994 for this example.

a <- hav_bg_1964
neg <- a$tot -a$pos
pos <- a$pos
age <- a$age
tot <- a$tot

Muench model

Proposed model

[@muench_derivation_1934] suggested to model the infection process with so-called "catalytic model", in which the distribution of the time spent in the susceptible class in SIR model is exponential with rate $\beta$

$$ \pi(a) = k(1 - e^{-\beta a} ) $$

Where:

• $\pi$ is the seroprevalence at age $a$

• $1 - k$ is the proportion of population that stay uninfected for a lifetime

• $a$ is the variable age

Under this catalytic model and assuming that $k = 1$, force infection would be $\lambda(a) = \beta$

Fitting data

Muench's model can be estimated by either defining k = 1 (a degree one linear predictor, note that it is irrelevant to the k in the proposed model) or setting the type = "Muench".

muench1 <- polynomial_model(age, pos = pos, tot = tot, k = 1)
summary(muench1$info)

muench2 <- polynomial_model(age, pos = pos, tot = tot, type = "Muench")
summary(muench2$info)

We can plot any model with the plot() function.

plot(muench2) 

Griffith model

Proposed model

Griffith proposed a model for force of infection as followed

$$ \lambda(a) = \beta_1 + 2\beta_2a
$$

Which can be estimated using a GLM where the for which the linear predictor was $\eta(a) = \beta_1 + \beta_2a^{2}$

Fitting data

Similarly, we can estimate Griffith's model either by defining k = 2, or setting the type = "Griffith"

gf_model <- polynomial_model(age, pos = pos, tot = tot, type = "Griffith")
plot(gf_model)

Grenfell and Anderson model

Proposed model

[@grenfell_estimation_1985] extended the models of Muench and Griffiths further suggest the use of higher order polynomial functions to model the force of infection which assumes prevalence model as followed

$$ \pi(a) = 1 - e^{-\Sigma_i \beta_i a^i} $$

Which implies that force of infection equals $\lambda(a) = \Sigma \beta_i i a^{i-1}$

Fitting data

And Grenfell and Anderson's model.

grf_model <- polynomial_model(age, pos = pos, tot = tot, type = "Grenfell")
plot(grf_model)

---

Nonlinear models

Refer to Chapter 6.1.2

Farrington model

Proposed model

For Farrington's model, the force of infection was defined non-negative for all a $\lambda(a) \geq 0$ and increases to a peak in a linear fashion followed by an exponential decrease

$$ \lambda(a) = (\alpha a - \gamma)e^{-\beta a} + \gamma $$

Where $\gamma$ is called the long term residual for FOI, as $a \rightarrow \infty$ , $\lambda (a) \rightarrow \gamma$

Integrating $\lambda(a)$ would results in the following non-linear model for prevalence

$$ \pi (a) = 1 - e^{-\int_0^a \lambda(s) ds} \ = 1 - exp{ \frac{\alpha}{\beta}ae^{-\beta a} + \frac{1}{\beta}(\frac{\alpha}{\beta} - \gamma)(e^{-\beta a} - 1) -\gamma a } $$

Fitting data

Use farrington_model() to fit a Farrington's model.

rubella <- rubella_uk_1986_1987
rubella$neg <- rubella$tot - rubella$pos

farrington_md <- farrington_model(
   rubella$age, pos = rubella$pos, tot = rubella$tot,
   start=list(alpha=0.07,beta=0.1,gamma=0.03)
   )
plot(farrington_md)

Weibull model

Proposed model

For a Weibull model, the prevalence is given by

$$ \pi (d) = 1 - e^{ - \beta_0 d ^ {\beta_1}} $$

Where $d$ is exposure time (difference between age of injection and age at test)

The model was reformulated as a GLM model with log - log link and linear predictor using log(d)

$$\eta(d) = log(\beta_0) + \beta_1 log(d)$$

Thus implies that the force of infection is a monotone function of the exposure time as followed

$$ \lambda(d) = \beta_0 \beta_1 d^{\beta_1 - 1} $$

Fitting data

Use weibull_model() to fit a Weibull model.

hcv <- hcv_be_2006[order(hcv_be_2006$dur), ]
dur <- hcv$dur
infected <- hcv$seropositive

wb_md <- weibull_model(
   t = dur,
   status = infected
   )
plot(wb_md) 

---

Fractional polynomial model

Refer to Chapter 6.2

Proposed model

Fractional polynomial model generalize conventional polynomial class of functions. In the context of binary responses, a fractional polynomial of degree $m$ for the linear predictor is defined as followed

$$ \eta_m(a, \beta, p_1, p_2, ...,p_m) = \Sigma^m_{i=0} \beta_i H_i(a) $$

Where $m$ is an integer, $p_1 \le p_2 \le... \le p_m$ is a sequence of powers, and $H_i(a)$ is a transformation given by

$$ H_i = \begin{cases} a^{p_i} & \text{ if } p_i \neq p_{i-1}, \ H_{i-1}(a) \times log(a) & \text{ if } p_i = p_{i-1}, \end{cases} $$

Best power selection

Use find_best_fp_powers() to find the powers which gives the lowest deviance score

hav <- hav_be_1993_1994
best_p <- find_best_fp_powers(
  hav$age, pos = hav$pos,tot = hav$tot,
  p=seq(-2,3,0.1), mc=FALSE, degree=2, link="cloglog"
)
best_p

Fitting data

Use fp_model() to fit a fractional polynomial model

model <- fp_model(
  hav$age, pos = hav$pos, tot = hav$tot,
  p=c(1.5, 1.6), link="cloglog")
compute_ci.fp_model(model)
plot(model)

serosv documentation built on Oct. 18, 2024, 5:07 p.m.