Skip to contents

This function implements an EM algorithm to estimate the relative case fatality ratio between two groups when reporting rates are time-varying and deaths are lagged because of survival time.

Usage

EMforCFR(assumed.nu, alpha.start.values, full.data, max.iter = 50,
  verb = FALSE, tol = 1e-10, SEM.var = TRUE)

Arguments

assumed.nu

a vector of probabilities corresponding to the survival distribution, i.e. nu[i]=Pr(surviving i days | fatal case)

alpha.start.values

a vector starting values for the reporting rate parameter of the GLM model. This must have length which corresponds to one less than the number of unique integer values of full.dat[,"new.times"].

full.data

A matrix of observed data. See description below.

max.iter

The maximum number of iterations for the EM algorithm and the accompanying SEM algorithm (if used).

verb

An indicator for whether the function should print results as it runs.

tol

A tolerance to use to test for convergence of the EM algorithm.

SEM.var

If TRUE, the SEM algorithm will be run in addition to the EM algorithm to calculate the variance of the parameter estimates.

Value

A list with the following elements

naive.rel.cfr

the naive estimate of the relative case fatality ratio

glm.rel.cfr

the reporting-rate-adjusted estimate of the relative case fatality ratio

EM.rel.cfr

the lag-adjusted estimate of the relative case fatality ratio

EM.re.cfr.var

the variance for the log-scale lag-adjusted estimator taken from the final M-step

EM.rel.cfr.var.SEM

the Supplemented EM algorithm variance for the log-scale lag-adjusted estimator

EM.rel.cfr.chain

a vector of the EM algorithm iterates of the lag-adjusted relative CFR estimates

EMiter

the number of iterations needed for the EM algorithm to converge

EMconv

indicator for convergence of the EM algorithm. 0 indicates all parameters converged within max.iter iterations. 1 indicates that the estimate of the relative case fatality ratio converged but other did not. 2 indicates that the relative case fatality ratio did not converge.

SEMconv

indicator for convergence of SEM algorithm. Same scheme as EMconv.

ests

the coefficient estimates for the model

ests.chain

a matrix with all of the coefficient estimates, at each EM iteration

DM

the DM matrix from the SEM algorithm

DMiter

a vector showing how many iterations it took for the variance component to converge in the SEM algorithm

Details

The data matrix full.data must have the following columns:

grp

a 1 or a 2 indicating which of the two groups, j, the observation is for.

new.times

an integer value representing the time, t, of observation.

R

the count of recovered cases with onset at time t in group j.

D

the count of deaths which occurred at time t in groupo j (note that these deaths did not have disease onset at time t but rather died at time t).

N

the total cases at t, j, or the sum of R and D columns.

Examples

    ## This is code from the CFR vignette provided in the documentation.

data(simulated.outbreak.deaths)
min.cases <- 10
N.1 <- simulated.outbreak.deaths[1:60, "N"]
N.2 <- simulated.outbreak.deaths[61:120, "N"]
first.t <- min(which(N.1 > min.cases & N.2 > min.cases))
last.t <- max(which(N.1 > min.cases & N.2 > min.cases))
idx.for.Estep <- first.t:last.t
new.times <- 1:length(idx.for.Estep)
simulated.outbreak.deaths <- cbind(simulated.outbreak.deaths, new.times = NA)
simulated.outbreak.deaths[c(idx.for.Estep, idx.for.Estep + 60), "new.times"] <- rep(new.times, + 2)
assumed.nu = c(0, 0.3, 0.4, 0.3)
alpha.start <- rep(0, 22)

## caution! this next line may take several minutes (5-10, depanding on
##    the speed of your machine) to run.
if (FALSE) cfr.ests <- EMforCFR(assumed.nu = assumed.nu,
                              alpha.start.values = alpha.start,
                              full.data = simulated.outbreak.deaths,
                              verb = FALSE,
                              SEM.var = TRUE,
                              max.iter = 500,
                              tol = 1e-05)