History of Exposure

PR time series to EIR time series

This vignette relies on software developed for RAMP Simulation-Based Analytics (see SimBA)

Pseudo-Data

Suppose we have data describing PfPR and we have case data from facilities. From these two signals, we think the PfPR over the last six years looks something like this,

set.seed(23)
t = c(0:6)*365 
tt = seq(0, 6*365, by = 30)
y = c(0.27, .28, .32, .15, .13, .18, 0.23)
seas = sin(2*pi*(tt-120)/365)*.03
trend = spline(t, y, xout=tt)$y
noise = rnorm(length(tt), 0, .3*y)^2


signal = trend+seas+noise
plot(tt, signal, type = "l", ylim = c(0, 0.45), xlab = "time", ylab = "trend", lwd=2)
lines(tt, trend, col = "darkred")
points(t,y)

#segments(950, 0, 950, 0.4, col = "sienna", lwd=2)
#segments(1315, 0, 1315, 0.4, col = "sienna", lwd=2)

Setting up a model

We start by making assumptions about the model for infection and immunity, about some aspects of seasonality, and etc.

We would like to start with a simple model. We use ramp.xds and we want to use the SIP module from ramp.library. We will also need some functions from ramp.work

library(ramp.xds)
library(ramp.library)
library(ramp.work)

Infection Dynamics

Each model asks us to make some assumptions about the parameters affecting malaria infection and immunity, including drug taking. The SIP model is explained in the documentation:

Xo = list(xi=2/365, rho=0.1)

Seasonality

We need to make some kind of assumption about the form for seasonal exposure to the EIR. We’ll adjust the phase in a bit, but for now, we’re going to specify the shape. This is explained in the documentation for make_function.sin in ramp.xds (offsite).

p0 <- makepar_F_sin(floor=2)
F_s0 = make_function(p0) 
plot(tt, F_s0(tt), type = "l", ylim = c(c(0, 2)), ylab = "Seasonal Effect", xlab = "Time")

Travel Malaria

If you want to think about travel malaria, it’s pretty easy. I assume this is not of great interest.

Build Model

We’re going to be using xds_setup_cohort for this.

mod0 <- xds_setup_cohort(Xname = "SIP", Xopts=Xo, F_season=F_s0)

Progressive Calibration

We want to develop a model that matches the features of some data set. To do so, we have described a protocol for doing this in an orderly way:

Mean EIR
Phase matching to align the peaks

Scaling

mod0 <- xde_scaling_eir(mod0)

eirpr0 <- xde_pr2eir(mean(trend), mod0)
eirpr0

$pr
[1] 0.2168994

$eir
[1] 0.007474159

$errors
numeric(0)

The following shows the scaling relationship and the fitted point.

plot_eirpr(mod0)
with(eirpr0, points(365*eir, pr))

This gives us an annual EIR value:

eirpr0$eir-> eir0
eir0*365

[1] 2.728068

Now, we have a model that outputs a time series that has the right mean.

mod0 <- set_eir(eir0, mod0)  
mod0 <- xds_solve_cohort(mod0, A=3)
mod0 <- last_to_inits(mod0) 
mod0 <- xds_solve_cohort(mod0, A=7)
XH <- get_XH(mod0)

with(XH, plot(time, true_pr, type = "l", ylim = c(0,0.45), col = "darkred"))
lines(tt, signal)

Phase Matching

An easy way to get the phase right is to compare the phase in data to the phase of the simulation.

ramp.work has a function called mean_phase_peak that computes the average peak. In the pseudo-data I generated, malaria peaks around the 210^{th} day of the year:

mean_phase_peak(tt, signal)$mean_peak -> sig_pk 
mean_phase_peak(tt, smooth(signal))$mean_peak -> sig_pk2
c(sig_pk, sig_pk2)

[1] 212.7703 198.9865

mod_pk <- with(XH,mean_phase_peak(time, true_pr)$mean_peak) 
mod_pk

[1] 142.6953

The difference is about 67 days:

sig_pk - mod_pk -> phase_shift
phase_shift

[1] 70.07496

The new seasonal signal looks a lot like the old one, but it is shifted by around 67 days.

p1 <- makepar_F_sin(floor=2, phase=phase_shift)
F_s1 = make_function(p1) 
plot(tt, F_s0(tt), type = "l", ylim = c(c(0, 2)), ylab = "Seasonal Effect", xlab = "Time")
lines(tt, F_s1(tt), col = "darkred")

We’re trying to get close to 210

mod1 <- xds_setup_cohort(Xname = "SIP", Xopts=Xo, F_season=F_s1, eir=eir0)
mod1 <- xds_solve_cohort(mod1, A=3)
mod1 <- last_to_inits(mod1) 
mod1 <- xds_solve_cohort(mod1, A=7)
XH <- get_XH(mod1)
with(XH, mean_phase_peak(time, true_pr)$mean_peak)

[1] 211.2891

with(XH, plot(time, true_pr, type = "l", ylim = c(0,0.45), col = "darkred"))
lines(tt, signal)

Trend

The next step is to configure a function to correct the trend. Spline functions work pretty well, but any function of time would do.

F_spline = function(tt){
  t = c(0, 1.8, 3, 4, 6)*365
  y = c(1.1, 1.6, .6, .7, .8)
  spline(t, y, xout=(tt))$y
}
plot(tt, F_spline(tt), type = "l")

It’s not perfect, but it gets a little closer.

make_F_yspline = function(y, tt, p0=-90){
  t = seq(p0, max(tt), length.out=length(y))
  return(function(tt){spline(t, y, xout=(tt))$y})
} 
yy = c(1.1, 1.6, .6, .7, .8)
F_yspline = make_F_yspline(yy, tt)
plot(tt, F_yspline(tt), type = "l")

mod2 <- last_to_inits(mod1) 
mod2 <- xds_solve_cohort(mod2, A=2)
mod2$EIRpar$F_trend  <- F_spline 
mod2 <- xds_solve_cohort(mod2, A=7)
XH <- get_XH(mod2)
plot(tt, signal, type = "l", lwd=2,ylim = c(0, 0.45))
with(XH, lines(time, true_pr, col = "darkred"))

fit_it = function(n, data, mod){
  gof = function(y, signal, mod){
    tt = seq(0, 6*365, by = 30)
    mod$EIRpar$F_trend = make_F_yspline(y, tt)
    mod <- xds_solve_cohort(mod, A=6, da=30)
    pr <- get_XH(mod)$true_pr
    sum((signal-pr)^2) 
  }
  nlm(gof, rep(1, n), signal=data, mod=mod)$estimate
}

mod3 <- mod2
fit_it(8, signal, mod2) -> yy
newF <- make_F_yspline(yy, tt)

plot(tt, newF(tt), type = "l", ylab = "Fitted Trend", xlab = "Time")
segments(0, 1, max(tt), 1, col = grey(0.5))

mod3$EIRpar$F_trend = newF
mod3 <- xds_solve_cohort(mod3, A=6, da=30)
XH <- get_XH(mod3)
plot(tt, signal, type = "l", lwd=2,ylim = c(0, 0.45))
with(XH, lines(time, true_pr, col = "darkred"))

hist(signal-XH$true_pr)