started to add more general model

R/simulate_serosurvey.R

@@ -223,6 +223,36 @@ probability_seropositive_by_age <- function(
   return(df)
 }
 
+#' Generate probabilities of seropositivity by age based on an age-varying FOI model.
+#'
+#' This function calculates the probabilities of seropositivity by age based on an age-varying FOI model.
+#' It takes into account the FOI and the rate of seroreversion.
+#'
+#' @param foi A dataframe containing the force of infection (FOI) values for different ages.
+#'            It should have two columns: 'age' and 'foi'.
+#' @param seroreversion_rate A non-negative numeric value representing the rate of seroreversion.
+#'
+#' @return A dataframe with columns 'age' and 'seropositivity'.
+probability_seropositive_general_model_by_age <- function(
+    foi,
+    seroreversion_rate) {
+
+  ages <- seq_along(foi$age)
+
+  probabilities <- probability_exact_age_varying(
+    ages = ages,
+    fois = foi$foi,
+    seroreversion_rate = seroreversion_rate
+  )
+
+  df <- data.frame(
+    age = ages,
+    seropositivity = probabilities
+  )
+
+  return(df)
+}
+
 
 #' Add bins based on age intervals.
 #'
@@ -729,3 +759,5 @@ simulate_serosurvey <- function(
 
   return(serosurvey)
 }
+
+

vignettes/simulating_serosurveys.Rmd

@@ -320,3 +320,237 @@ serosurvey_serorevert %>%
   facet_wrap(~type)
 ```
 
+# Simulating from a general serological model
+We now illustrate an example age- and time-dependent FOI model for representing the prevalence of HIV (in the absence of treatment). In the absence of treatment, HIV invariably ends in death, and death usually occurs after the onset of AIDS which can occur typically many years after death. This means that the time from infection to death has a characteristic waiting time, which year we model as being around 10 years long.
+
+Our model considers a single susceptible group, $S^{\tau}(t)$ indexed by their birth year, $\tau$, and calendar time, $t$. There are then a series of consecutive seropositive states, $X_i^{\tau}(t)$, for $i\in[1,10]$; then finally a dead state, $D^{\tau}(t)$. By distributing the seropositive population across a range of states, this model effectively institutes a delay from infection until death. Our model takes the following form:
+
+\begin{equation}\label{eq:serocatalytic_hiv}
+    \begin{aligned}
+        \frac{dS^{\tau}(t)}{dt} &= -u(t-\tau)v(t) S^{\tau}(t)\\
+        \frac{dX_1^{\tau}(t)}{dt} &= u(t-\tau)v(t) S^{\tau}(t) - X_1^{\tau}(t)\\
+        \frac{dX_i^{\tau}(t)}{dt} &= X_{i-1}^{\tau}(t) - X_{i}^{\tau}(t), \quad \forall i \in [2,10],\\
+        \frac{dD^{\tau}(t)}{dt} &= X_{10},
+    \end{aligned}
+\end{equation}
+
+where $S^{\tau}(0)=1$ and all other compartments have initial conditions set to zero. We assume that $u(.)$ and $v(.)$ are piecewise-constant with one-year pieces, with their values the same as those from the sexually transmitted disease example above.
+
+This model does not fit neatly into our other classes of serocatalytic models and, more generally, we cannot produce specific methods for any foreseeable model structure. But it is possible to produce a general strategy for solving such models by rewriting them as linear systems of ordinary differential equations.
+
+Since the system is linear, we can write it as a vector differential equation by defining $\boldsymbol{Y}^{\tau}(t):=[S^{\tau}(t),X_1^{\tau}(t),...,X_{10}^{\tau}(t),D^{\tau}(t)]'$ as a vector of states and $\bar{\boldsymbol{A}}_{t,\tau}$ as a matrix of constants that are fixed for a given year and year of birth:
+
+\begin{equation}
+    \bar{\boldsymbol{A}}_{t,\tau} := \begin{bmatrix}
+-\bar u(t-\tau) \bar v(t) & 0 & \cdots \\
+\bar u(t-\tau) \bar v(t) & -1 & 0 & \cdots\\
+0 & 1 & -1 & 0 & \cdots\\
+0 & 0 & 1 & -1 & 0 & \cdots\\
+\vdots & \vdots & \ddots & \ddots & \ddots & \ddots\\
+0 & \cdots & & 0 & 1 & 0
+\end{bmatrix}.
+\end{equation}
+
+Using this matrix, we can write and solve the system of equations:
+
+\begin{equation}
+    \frac{d\boldsymbol{Y}^{\tau}(t)}{dt} = \bar{\boldsymbol{A}}_{t,\tau} \boldsymbol{Y}^{\tau}(t) \implies \boldsymbol{Y}^{\tau}(t+1) = \exp(\bar{\boldsymbol{A}}_{t,\tau}) \boldsymbol{Y}^{\tau}(t).
+\end{equation}
+
+where we have used matrix exponentials to find the solution because the individual algebraic expressions are unwieldy with this many compartments. This matrix form means that a general solution for the system can be written down as follows:
+
+\begin{equation}
+    \boldsymbol{Y}^{\tau}(t) = \exp\left(\sum_{t'=1}^{t}\bar{\boldsymbol{A}}_{t',\tau}\right) \boldsymbol{Y}^{\tau}(0).
+\end{equation}
+
+We now illustrate how such a system can be solved and used to simulate a serosurvey.
+
+```{r}
+u <- foi_df_age$foi
+v <- foi_df_time$foi
+
+construct_A <- function(t, tau, u, v) {
+  u_bar <- u[t - tau]
+  v_bar <- v[t]
+  
+  A <- diag(-1, ncol = 12, nrow = 12)
+  A[row(A) == (col(A) + 1)] <- 1
+  A[1, 1] <- -u_bar * v_bar
+  A[2, 1] <- u_bar * v_bar
+  A[12, 12] <- 0
+  
+  A
+}
+
+initial_conditions <- rep(0, 12)
+initial_conditions[1] <- 1
+
+sum_of_A <- function(t, tau, construct_A_fn, ...) {
+  k <- 1
+  for(i in (tau + 1):t) {
+    if(k == 1)
+      A <- construct_A_fn(i, tau, ...)
+    else {
+      tmp <- construct_A_fn(i, tau, ...)
+      A <- A + tmp
+    }
+    k <- k + 1
+  }
+  A
+}
+
+library(Matrix)
+# solving for all groups born between tau=1 and tau=80
+prop_seropositive_of_alive <- vector(length = 80)
+for(i in seq_along(prop_seropositive_of_alive)) {
+  A_sum <- sum_of_A(80, 80 - i, construct_A, u, v)
+  Y <- Matrix::expm(A_sum) %>% as.matrix() %*% initial_conditions
+  prop_seropositive_of_alive[i] <- sum(Y[2:11]) / (1 - Y[12])
+}
+
+plot(1:80, prop_seropositive_of_alive)
+```
+
+Compare with time-varying model
+```{r}
+construct_A <- function(t, tau, lambda) {
+  A <- matrix(0, ncol = 2, nrow = 2)
+  A[1, 1] <- -lambda[t]
+  A[2, 1] <- lambda[t]
+  A
+}
+
+lambda <- runif(80, 0, 0.01)
+initial_conditions <- c(1, 0)
+prop_seropositive_of_alive <- vector(length = 80)
+for(i in seq_along(prop_seropositive_of_alive)) {
+  A_sum <- sum_of_A(80, 80 - i, construct_A, lambda)
+  Y <- Matrix::expm(A_sum) %>% as.matrix() %*% initial_conditions
+  prop_seropositive_of_alive[i] <- Y[2]
+}
+
+plot(1:80, prop_seropositive_of_alive)
+seroreversion_system <- data.frame(
+  age=1:80,
+  seropositivity=prop_seropositive_of_alive
+) %>% 
+  mutate(type="system")
+
+
+foi_df <- data.frame(
+  year=seq(1, 80, 1),
+  foi=lambda
+)
+
+seropositive_true <- probability_seropositive_by_age(
+  model = "time",
+  foi = foi_df,
+  seroreversion_rate = 0
+) %>% 
+  mutate(type="official")
+
+seroreversion_system %>% 
+  bind_rows(seropositive_true) %>% 
+  ggplot(aes(x=age, y=seropositivity)) +
+  geom_line(aes(colour=type))
+
+```
+
+Compare with an age-varying model
+```{r}
+construct_A <- function(t, tau, lambda) {
+  A <- matrix(0, ncol = 2, nrow = 2)
+  A[1, 1] <- -lambda[t - tau]
+  A[2, 1] <- lambda[t - tau]
+  A
+}
+
+lambda <- runif(80, 0, 0.01)
+initial_conditions <- c(1, 0)
+prop_seropositive_of_alive <- vector(length = 80)
+for(i in seq_along(prop_seropositive_of_alive)) {
+  A_sum <- sum_of_A(80, 80 - i, construct_A, lambda)
+  Y <- Matrix::expm(A_sum) %>% as.matrix() %*% initial_conditions
+  prop_seropositive_of_alive[i] <- Y[2]
+}
+
+plot(1:80, prop_seropositive_of_alive)
+seroreversion_system <- data.frame(
+  age=1:80,
+  seropositivity=prop_seropositive_of_alive
+) %>% 
+  mutate(type="system")
+
+
+foi_df <- data.frame(
+  age=seq(1, 80, 1),
+  foi=lambda
+)
+
+seropositive_true <- probability_seropositive_by_age(
+  model = "age",
+  foi = foi_df,
+  seroreversion_rate = 0
+) %>% 
+  mutate(type="official")
+
+seroreversion_system %>% 
+  bind_rows(seropositive_true) %>% 
+  ggplot(aes(x=age, y=seropositivity)) +
+  geom_line(aes(colour=type))
+```
+
+Check age and time
+```{r}
+construct_A <- function(t, tau, u, v) {
+  A <- matrix(0, ncol = 2, nrow = 2)
+  u_bar <- u[t - tau]
+  v_bar <- v[t]
+  
+  A[1, 1] <- -u_bar * v_bar
+  A[2, 1] <- u_bar * v_bar
+  A
+}
+
+u <- foi_df_age$foi
+v <- foi_df_time$foi
+initial_conditions <- c(1, 0)
+prop_seropositive_of_alive <- vector(length = 80)
+for(i in seq_along(prop_seropositive_of_alive)) {
+  A_sum <- sum_of_A(80, 80 - i, construct_A, u, v)
+  Y <- Matrix::expm(A_sum) %>% as.matrix() %*% initial_conditions
+  prop_seropositive_of_alive[i] <- Y[2]
+}
+
+seroreversion_system <- data.frame(
+  age=1:80,
+  seropositivity=prop_seropositive_of_alive
+) %>% 
+  mutate(type="system")
+
+foi_df <- expand.grid(
+  year=foi_df_time$year,
+  age=foi_df_age$age
+) %>% 
+  left_join(foi_df_age, by="age") %>% 
+  rename(foi_age=foi) %>% 
+  left_join(foi_df_time, by="year") %>% 
+  rename(foi_time=foi) %>% 
+  mutate(foi = foi_age * foi_time) %>% 
+  select(-c("foi_age", "foi_time")) %>% 
+  mutate(birth_year = year - age) %>% 
+  filter(birth_year >= 1945) %>% 
+  arrange(birth_year, age)
+
+serosurvey_true <- probability_seropositive_by_age(
+  model = "age-time",
+  foi = foi_df %>% select(-birth_year)
+) %>% 
+  mutate(type="official")
+
+seroreversion_system %>% 
+  bind_rows(serosurvey_true) %>% 
+  ggplot(aes(x=age, y=seropositivity)) +
+  geom_line(aes(colour=type))
+```
+