Fit an Endemic-Only twinstim as a Poisson-glm
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 λ([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.
glm_epidataCS(formula, data, ...)
a glm
Sebastian Meyer
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))
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.