glm_epidataCS: Fit an Endemic-Only 'twinstim' as a Poisson-'glm'

In surveillance: Temporal and Spatio-Temporal Modeling and Monitoring of Epidemic Phenomena

Fit an Endemic-Only twinstim as a Poisson-glm

Description

An endemic-only twinstim is equivalent to a Poisson regression model for the aggregated number of events, Y_{[t][\bm{s}],k}, by time-space-type cell. The rate of the corresponding Poisson distribution is e_{[t][\bm{s}]} \cdot \lambda([t],[\bm{s}],k), where e_{[t][\bm{s}]} = |[t]| |[\bm{s}]| is a multiplicative offset. Thus, the glm function can be used to fit an endemic-only twinstim. However, wrapping in glm is usually slower.

Usage

glm_epidataCS(formula, data, ...)

Arguments

formula	an endemic model formula without response, comprising variables of data$stgrid and possibly the variable type for a type-specific model.
data	an object of class "epidataCS".
...	arguments passed to glm. Note that family and offset are fixed internally.

Value

a glm

Author(s)

Sebastian Meyer

Examples

data("imdepi", "imdepifit")

## Fit an endemic-only twinstim() and an equivalent model wrapped in glm()
fit_twinstim <- update(imdepifit, epidemic = ~0, siaf = NULL, subset = NULL,
                       optim.args=list(control=list(trace=0)), verbose=FALSE)
fit_glm <- glm_epidataCS(formula(fit_twinstim)$endemic, data = imdepi)

## Compare the coefficients
cbind(twinstim = coef(fit_twinstim), glm = coef(fit_glm))


### also compare to an equivalent endemic-only hhh4() fit

## first need to aggregate imdepi into an "sts" object
load(system.file("shapes", "districtsD.RData", package="surveillance"))
imdsts <- epidataCS2sts(imdepi, freq = 12, start = c(2002, 1),
                        neighbourhood = NULL, tiles = districtsD,
                        popcol.stgrid = "popdensity")

## determine the correct offset to get an equivalent model
offset <- 2 * rep(with(subset(imdepi$stgrid, !duplicated(BLOCK)),
                  stop - start), ncol(imdsts)) *
          sum(districtsD$POPULATION) * population(imdsts)

## fit the model using hhh4()
fit_hhh4 <- hhh4(imdsts, control = list(
    end = list(
        f = addSeason2formula(~I(start/365-3.5), period=365, timevar="start"),
        offset = offset
    ), family = "Poisson", subset = 1:nrow(imdsts),
    data = list(start=with(subset(imdepi$stgrid, !duplicated(BLOCK)), start))))

summary(fit_hhh4)
stopifnot(all.equal(coef(fit_hhh4), coef(fit_glm), check.attributes=FALSE))

surveillance documentation built on Nov. 4, 2024, 3 a.m.