ptpi: Estimate S_0 = S(t_0) using PTPI In davidearn/fastbeta: Fast Estimation of Time-Varying Transmission Rates

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.