data("measlesWeserEms")
## data contains adjaceny orders as required for parametric weights
plot(measlesWeserEms, type = observed ~ unit, labels = TRUE)
neighbourhood(measlesWeserEms)[1:6,1:6]
max(neighbourhood(measlesWeserEms)) # max order is 5
## fit a power-law decay of spatial interaction
## in a hhh4 model with seasonality and random intercepts in the endemic part
measlesModel <- list(
ar = list(f = ~ 1),
ne = list(f = ~ 1, weights = W_powerlaw(maxlag=5)),
end = list(f = addSeason2formula(~-1 + ri(), S=1, period=52)),
family = "NegBin1")
## fit the model
set.seed(1) # random intercepts are initialized randomly
measlesFit <- hhh4(measlesWeserEms, measlesModel)
summary(measlesFit) # "neweights.d" is the decay parameter d
coefW(measlesFit)
## plot the spatio-temporal weights o_ji^-d / sum_k o_jk^-d
## as a function of adjacency order
plot(measlesFit, type = "neweights", xlab = "adjacency order")
## normalization => same distance does not necessarily mean same weight.
## to extract the whole weight matrix W: getNEweights(measlesFit)
## visualize contributions of the three model components
## to the overall number of infections (aggregated over all districts)
plot(measlesFit, total = TRUE)
## little contribution from neighbouring districts
if (surveillance.options("allExamples")) {
## simpler model with autoregressive effects captured by the ne component
measlesModel2 <- list(
ne = list(f = ~ 1, weights = W_powerlaw(maxlag=5, from0=TRUE)),
end = list(f = addSeason2formula(~-1 + ri(), S=1, period=52)),
family = "NegBin1")
measlesFit2 <- hhh4(measlesWeserEms, measlesModel2)
## omitting the separate AR component simplifies model extensions/selection
## and interpretation of covariate effects (only two predictors left)
plot(measlesFit2, type = "neweights", exclude = NULL, xlab = "adjacency order")
## strong decay, again mostly within-district transmission
## (one could also try a purely autoregressive model)
plot(measlesFit2, total = TRUE,
legend.args = list(legend = c("epidemic", "endemic")))
## almost the same RMSE as with separate AR and NE effects
c(rmse1 = sqrt(mean(residuals(measlesFit, "response")^2)),
rmse2 = sqrt(mean(residuals(measlesFit2, "response")^2)))
}
Run the code above in your browser using DataLab