aggregate_moments: Aggregation of stationary or predictive moments

In jbracher/hhh4addon: Extensions to endemic-epidemic timeseries modeling from package surveillance

Aggregation of stationary or predictive moments

Description

Aggregation of stationary or predictive moments as calculated using stationary_moments or predictive_moments.

Usage

aggregate_moments(momentsObj, aggregation_matrix, by_timepoint = FALSE)

Arguments

momentsObj	an object of class moments_hhh4 containing stationary or predictive moments, as returned by stationary_moments or predictive_moments
aggregation_matrix	an aggregation matrix with either momentsObj$n_units columns (for aggregation across units while keeping the temporal structure; set option by_timepoint = TRUE in this case) or length(momentsObj$mu_vector) (for aggregation that does not preserve the temporal structure; set option by_timepoint = FALSE).
by_timepoint	logical: is aggregation only across units while preserving the temporal structure? Note that the new moments_hhh4 object cannot have the condition , mu_matrix, var_matrix and cov_array elements if the temporal structure is given up.

Value

An object of class moments_hhh4 representing the new prediction.

Examples

# load data:
data("noroBL")

########
# fit a bivariate model:
controlBL <- list(end = list(f = addSeason2formula(~ -1 + fe(1, unitSpecific = TRUE))),
                  ar = list(f = ~ -1 + fe(1, unitSpecific = TRUE)),
                  ne = list(f = ~ -1 + fe(1, unitSpecific = TRUE)),
                  family = "NegBinM", subset = 2:260) # not a very parsimonious parametrization, but feasible
fitBL <- hhh4(noroBL, control = controlBL)
pred_mom <- predictive_moments(fitBL, t_condition = 260, lgt = 52, return_Sigma = TRUE)
# Sigma is required in order to aggregate predictions.

#########
# plot predictions for two regions:
par(mfrow = 1:2)
fanplot_prediction(pred_mom, unit = 1, main = "Bremen")
fanplot_prediction(pred_mom, unit = 2, main = "Lower Saxony")

#########
# aggregation 1: combine the two regions
aggr_matr_pool <- matrix(1, ncol = 2)
# specify by_timepoint = TRUE to keep the temporal structure and aggregate only
# counts from the same week:
pred_mom_pooled <- aggregate_moments(pred_mom, aggr_matr_pool, by_timepoint = TRUE)
fanplot_prediction(pred_mom_pooled, unit = 1, ylim = c(0, 500), main = "Aggregation over regions")

#########
# aggregation 2: total burden in the two regions
aggr_matr_total_burden <- matrix(rep(c(1, 0, 0, 1), 52), nrow = 2,
                                 dimnames = list(c("Bremen", "Lower Saxony"),
                                                 NULL))
pred_mom_total_burden <- aggregate_moments(pred_mom, aggr_matr_total_burden)
plot_moments_by_unit(pred_mom_total_burden, main = "Total burdens")

#########
# works also with stationary moments:
stat_mom <- stationary_moments(fitBL, return_Sigma = TRUE)
stat_mom_pooled <- aggregate_moments(stat_mom, aggr_matr_pool, by_timepoint = TRUE)
stat_mom_total_burden <- aggregate_moments(stat_mom, aggr_matr_total_burden, by_timepoint = FALSE)
fanplot_stationary(stat_mom_pooled)
plot_moments_by_unit(stat_mom_total_burden, main = "Total burdens")

jbracher/hhh4addon documentation built on Feb. 16, 2024, 1:45 a.m.