ptpi: Estimate S_0 = S(t_0) using PTPI

Description Usage Arguments Value Mock vital data Missing data References Examples

View source: R/ptpi.R

Description

Using the method of peak-to-peak iteration (PTPI, see References), ptpi() estimates the initial number of susceptibles S0 = S(t0) from time series of incidence (must be roughly periodic), births, and natural mortality, observed at equally spaced time points tk = t0+kΔt (for k = 0,...,n), where Δt denotes the observation interval.

Usage

1
ptpi(df = data.frame(), par_list = list(), a, b, initial_S0_est, iter = 0L)

Arguments

df

A data frame with numeric columns:

Z

Incidence. Z[i] is the number of infections between time t = ti−2 and time t = ti−1. Z must be roughly periodic.

B

Births. B[i] is the number of births between time t = ti−2 and time t = ti−1.

mu

Natural mortality rate. mu[i] is the rate at time t = ti−1 expressed per unit Δt and per capita.

B is optional if hatN0 and nu are defined in par_list, and mu is optional if mu is defined in par_list (see Details).

par_list

A list containing:

hatN0

[ Ñ0 ] Population size at time t = 0 years.

nu

[ νc ] Birth rate expressed per unit Δt and relative to Ñ0 (if modeled as constant).

mu

[ μc ] Natural mortality rate expressed per unit Δt and per capita (if modeled as constant).

hatN0 and nu are optional if B is defined in df, and mu is optional if mu is defined in df (see Details).

a

Integer scalar. Index of first peak in df$Z, possibly obtained via get_peak_times().

b

Integer scalar. Index of last peak in df$Z in phase with first peak, possibly obtained via get_peak_times().

initial_S0_est

Numeric scalar. An initial estimate of S0 = S(t0).

iter

Integer scalar. The number of estimates of S0 = S(t0) to generate before stopping.

Value

A list containing:

S_mat

A numeric matrix with dimensions c(nrow(df), iter+1) containing the susceptible time series generated in each iteration.

S0

A numeric vector listing in order all 1+iter estimates of S0 = S(t0). Equal to S_mat[1, ].

S0_final

The final estimate of S0 = S(t0). Equal to S0[length(S0)].

SA

A numeric vector listing in order all 1+iter estimates of Sa = S(ta), where ta is the time point corresponding to row a in df. Equal to S_mat[a, ].

SA_final

The final estimate of Sa = S(ta). Equal to SA[length(SA)].

A list of the arguments of ptpi() is included as an attribute.

Mock vital data

If df$B is undefined in the function call, then df$B[i] gets the value with(par_list, nu * hatN0 * 1) for all i. If df$mu is undefined in the function call, then df$mu[i] gets the value with(par_list, mu) for all i.

Missing data

Missing values in df[, c("Z", "B", "mu")] are mostly not tolerated. At the moment, ptpi() makes no effort to impute them, so imputation must be done beforehand.

References

deJonge MS, Jagan M, Krylova O, Earn DJD. Fast estimation of time-varying transmission rates for infectious diseases.

Examples

 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
# Simulate 20 years of disease incidence,
# observed weekly
par_list <- make_par_list(dt_weeks = 1)
df <- make_data(
  par_list = par_list,
  n = 20 * 365 / 7, # 20 years is ~1042 weeks
  with_dem_stoch = TRUE,
  seed = 5
)

# Plot incidence time series, and note the
# apparent 1-year period
plot(Z ~ t_years, df,
  type = "l",
  xlab = "Time (years)",
  ylab = "Incidence"
)

# Find peaks in incidence time series
peaks <- get_peak_times(
  x = df$Z,
  period = 365 / 7, # 1 year is ~52 weeks
  bw_mavg = 6,
  bw_peakid = 8
)

# Verify that peaks were identified correctly
abline(v = df$t_years[peaks$all], lty = 2, col = "red")

# Index of first peak
a <- with(peaks, phase[1])

# Index of last peak in phase with first
b <- with(peaks, phase[length(phase)])

# Estimate `S0` from incidence via PTPI,
# starting from an erroneous initial guess
# (mock vital data generated internally)
ptpi_out <- ptpi(
  df = df["Z"],
  par_list = par_list,
  a = a,
  b = b,
  initial_S0_est <- df$S[1] * 4,
  iter = 25
)

# Sequence of estimates
ptpi_out$S0

# Relative error in final estimate
(ptpi_out$S0_final - df$S[1]) / df$S[1]

davidearn/fastbeta documentation built on June 14, 2020, 3:11 p.m.