Supplementary material

Area-level excess mortality in times of COVID-19 in Switzerland: geographical, socioeconomic and political determinants

Data

summary_table(exp_deaths_2020_year) %>% 
  flextable::flextable()

Age group	Sex	Observed	Expected (median)	Expected (lower bound)	Expected (upper bound)	Relative excess (median)	Relative excess (lower bound)	Relative excess (upper bound)
40-59	Female	1,713	1,769	1,642	1,895	0.97	0.90	1.04
40-59	Male	2,966	2,230	2,074	2,396	1.33	1.24	1.43
60-69	Female	2,611	2,592	2,421	2,753	1.01	0.95	1.08
60-69	Male	4,478	3,201	3,013	3,449	1.40	1.30	1.49
70-79	Female	6,203	5,028	4,708	5,376	1.23	1.15	1.32
70-79	Male	8,972	5,953	5,615	6,310	1.51	1.42	1.60
80+	Female	27,541	17,205	16,453	18,284	1.60	1.51	1.67
80+	Male	20,292	17,677	16,791	18,900	1.15	1.07	1.21
Total	Total	74,776	55,676	53,865	57,821	1.34	1.29	1.39

correlogram(exp_deaths_2020_year)

Models of observed and expected deaths by municipality

Step 1: iterative model development

To facilitate model development we start by only using the median excess mortality by municipality, age group and sex in 2020.

data1 = exp_deaths_2020_year %>% 
  group_by(canton, GMDNR, GMDNAME, age_group, id_space, sex, munici_observed, munici_pop,
           density_high, density_low, across(starts_with("sep")), border, lang_fr, lang_it,
           across(starts_with("vote")),type_urban,type_rural,type_periurban) %>% 
  summarise(munici_exp_deaths=median(munici_exp_deaths),
            munici_excess=median(munici_excess)) %>% 
  mutate(E=log(ifelse(munici_exp_deaths==0,1e-4,munici_exp_deaths))) %>% 
  ungroup()
data2 = exp_deaths_2020_year %>% 
  group_by(canton, GMDNR, GMDNAME, age_group, id_space, sex, munici_observed, munici_pop,
           density_high, density_low, across(starts_with("sep")), border, lang_fr, lang_it,
           across(starts_with("vote")),type_urban,type_rural,type_periurban) %>% 
  summarise(munici_exp_deaths=mean(munici_exp_deaths),
            munici_excess=mean(munici_excess)) %>% 
  mutate(E=log(ifelse(munici_exp_deaths==0,1e-4,munici_exp_deaths))) %>% 
  ungroup()
rm(exp_deaths_2020_year)

# gc()

hyper.iid = list(theta = list(prior = "pc.prec", param = c(1, 0.01)))
hyper.bym2 = list(theta1 = list("PCprior", c(1, 0.01)), 
                  theta2 = list("PCprior", c(0.5, 0.5)))
threads = parallel::detectCores()

Model 1.0: no covariates

We use a model structure similar to Poisson regression, where \(O_{t,i,j,k}\), the number of observed deaths during week \(t\) in municipality \(i\), age group \(j\) and sex group \(k\), depends on the number of expected deaths \(E_{t,i,j,k}\) based on historical data and a log-linear predictor \(\log \lambda = \alpha + \beta X\).

\[ O_i \sim \text{Poisson}(\lambda E_i ) \\ \] At start, \(\lambda\) only includes one intercept parameter \(\alpha\), so that the estimate of \(\exp(\alpha)\) can be interpreted as an average relative excess mortality (that is, the ratio of observed on expected) for 2020. By adding covariates to \(\lambda\), we aim to disentangle the various factors that are associated with excess mortality at the local level.

We implement this model in R-INLA, a Bayesian inference package that is especially adapted to spatial data. This is achieved in practice by including \(\log (E_{i,j,k})\) as an offset (although an alternative formulation based on the E argument exists). During model development, we compare different model versions based on the WAIC (lower values imply a better fit).

model1.0 = INLA::inla(munici_observed ~ 1 + offset(E),
                      data = data1,
                      family = "Poisson",
                      control.compute = list(config = TRUE, waic = TRUE),
                      quantiles = c(0.025, 0.5, 0.975),
                      num.threads = threads,
                      safe = TRUE)
summary(model1.0)


Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, 
   quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, 
   Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = 
   link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, 
   control.compute = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, control.fixed = 
   control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " 
   control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = 
   control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg 
   = inla.arg, num.threads = num.threads, ", " blas.num.threads = blas.num.threads, keep = 
   keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, 
   safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
Time used:
    Pre = 0.795, Running = 0.549, Post = 0.11, Total = 1.45 
Fixed effects:
             mean    sd 0.025quant 0.5quant 0.975quant  mode kld
(Intercept) 0.325 0.004      0.317    0.325      0.332 0.325   0

Watanabe-Akaike information criterion (WAIC) ...: 49278.96
Effective number of parameters .................: 3.53

Marginal log-Likelihood:  -24648.51 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.0$summary.fixed)[c(1,3,5)]

                mean 0.025quant 0.975quant
(Intercept) 1.383534   1.373679    1.39346

sum(data1$munici_observed)/sum(data1$munici_exp_deaths)

[1] 1.383562

As a sanity check, we find a relative excess mortality of 38% for 2020, that is coherent with a simple calculation (74,776 observed / 54,046 expected = 1.38). Remember that we excluded the age group 0-40, which explains why this is higher than numbers reported for Switzerland, generally around 10% for 2020. We can also look at the model fit and at the residuals. Obviously the model fit is not good here, as this basic model assumes a unique relative excess mortality for all areas, sexes and age groups.

Model 1.1: age and sex

We hypothesize that excess mortality affected different age and sex groups differently. We thus add the age group, the sex and the interaction of the two as covariates.

model1.1 = INLA::inla(munici_observed ~ - 1 + offset(E) +
                        sex:age_group,
                      data = data1,
                      family = "Poisson",
                      control.compute = list(config = TRUE, waic = TRUE),
                      quantiles = c(0.025, 0.5, 0.975),
                      num.threads = threads,
                      safe = TRUE)
summary(model1.1)


Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, 
   quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, 
   Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = 
   link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, 
   control.compute = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, control.fixed = 
   control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " 
   control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = 
   control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg 
   = inla.arg, num.threads = num.threads, ", " blas.num.threads = blas.num.threads, keep = 
   keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, 
   safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
Time used:
    Pre = 0.515, Running = 0.842, Post = 0.0968, Total = 1.45 
Fixed effects:
                          mean    sd 0.025quant 0.5quant 0.975quant  mode kld
sexFemale:age_group40-59 0.008 0.024     -0.040    0.008      0.055 0.008   0
sexMale:age_group40-59   0.354 0.018      0.318    0.354      0.390 0.354   0
sexFemale:age_group60-69 0.035 0.020     -0.003    0.035      0.074 0.035   0
sexMale:age_group60-69   0.417 0.015      0.388    0.417      0.447 0.417   0
sexFemale:age_group70-79 0.235 0.013      0.210    0.235      0.260 0.235   0
sexMale:age_group70-79   0.473 0.011      0.453    0.473      0.494 0.473   0
sexFemale:age_group80+   0.493 0.006      0.481    0.493      0.504 0.493   0
sexMale:age_group80+     0.149 0.007      0.136    0.149      0.163 0.149   0

Watanabe-Akaike information criterion (WAIC) ...: 47190.67
Effective number of parameters .................: 4.17

Marginal log-Likelihood:  -23651.68 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.1$summary.fixed)[c(1,3,5)]

                             mean 0.025quant 0.975quant
sexFemale:age_group40-59 1.007864  0.9612512   1.056738
sexMale:age_group40-59   1.424950  1.3745805   1.477166
sexFemale:age_group60-69 1.036069  0.9970809   1.076581
sexMale:age_group60-69   1.517741  1.4739324   1.562852
sexFemale:age_group70-79 1.264505  1.2334258   1.296368
sexMale:age_group70-79   1.605326  1.5724496   1.638889
sexFemale:age_group80+   1.636631  1.6174162   1.656075
sexMale:age_group80+     1.161138  1.1452712   1.177224

model1.1$waic$waic - model1.0$waic$waic

[1] -2088.29

As expected, the relative excess mortality varies a lot across age and sex groups. It’s very small in females aged 40-59 and 60-69 (in fact the data is compatible with no excess in both cases). It increases in females aged 70-79, and even more so aged 80+. It’s comparatively higher in males below 80, but somewhat surprisingly lower in males in age group 80+. This still corresponds to basic sanity checks with the data.

data1 %>% 
  group_by(sex,age_group) %>% 
  summarise(munici_observed=sum(munici_observed),
            munici_exp_deaths=sum(munici_exp_deaths)) %>% 
  mutate(ratio=munici_observed/munici_exp_deaths)

# A tibble: 8 × 5
# Groups:   sex [2]
  sex    age_group munici_observed munici_exp_deaths ratio
  <chr>  <chr>               <int>             <dbl> <dbl>
1 Female 40-59                1713             1699   1.01
2 Female 60-69                2611             2520.  1.04
3 Female 70-79                6203             4905   1.26
4 Female 80+                 27541            16828.  1.64
5 Male   40-59                2966             2081   1.43
6 Male   60-69                4478             2950   1.52
7 Male   70-79                8972             5588.  1.61
8 Male   80+                 20292            17476.  1.16

We observe an improvement of the model fit, not easy to spot on the plot because of the large number of points, but made clear by the large decrease in WAIC.

Model 1.2: spatial variability

We now account for spatial variability, first in a simple way using an i.i.d. random effect, so that all municipalities can vary independently from each other around a global average. Note that this “municipality effect” applies the same to all age and sex groups.

model1.2 = INLA::inla(munici_observed ~ - 1 + offset(E) +
                        sex:age_group +
                        f(id_space, model = "iid"),
                      data = data1,
                      family = "Poisson",
                      control.compute = list(config = TRUE, waic = TRUE),
                      quantiles = c(0.025, 0.5, 0.975),
                      num.threads = threads,
                      safe = TRUE)
summary(model1.2)


Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, 
   quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, 
   Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = 
   link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, 
   control.compute = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, control.fixed = 
   control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " 
   control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = 
   control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg 
   = inla.arg, num.threads = num.threads, ", " blas.num.threads = blas.num.threads, keep = 
   keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, 
   safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
Time used:
    Pre = 0.53, Running = 1.66, Post = 0.192, Total = 2.39 
Fixed effects:
                          mean    sd 0.025quant 0.5quant 0.975quant  mode kld
sexFemale:age_group40-59 0.011 0.024     -0.036    0.011      0.058 0.011   0
sexMale:age_group40-59   0.358 0.018      0.322    0.358      0.394 0.358   0
sexFemale:age_group60-69 0.039 0.020      0.001    0.039      0.078 0.039   0
sexMale:age_group60-69   0.420 0.015      0.391    0.420      0.450 0.420   0
sexFemale:age_group70-79 0.238 0.013      0.213    0.238      0.263 0.238   0
sexMale:age_group70-79   0.477 0.011      0.456    0.477      0.498 0.477   0
sexFemale:age_group80+   0.497 0.006      0.485    0.497      0.509 0.497   0
sexMale:age_group80+     0.153 0.007      0.139    0.153      0.167 0.153   0

Random effects:
  Name    Model
    id_space IID model

Model hyperparameters:
                          mean      sd 0.025quant 0.5quant 0.975quant    mode
Precision for id_space 3453.45 2510.49    1339.89  2741.78    9707.17 2116.51

Watanabe-Akaike information criterion (WAIC) ...: 47159.19
Effective number of parameters .................: 14.49

Marginal log-Likelihood:  -23646.50 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.2$summary.fixed)[c(1,3,5)]

                             mean 0.025quant 0.975quant
sexFemale:age_group40-59 1.011058  0.9641761   1.060222
sexMale:age_group40-59   1.430421  1.3796297   1.483087
sexFemale:age_group60-69 1.039892  1.0006107   1.080718
sexMale:age_group60-69   1.522570  1.4783840   1.568084
sexFemale:age_group70-79 1.268978  1.2375430   1.301220
sexMale:age_group70-79   1.610709  1.5773914   1.644743
sexFemale:age_group80+   1.643565  1.6234816   1.663948
sexMale:age_group80+     1.165188  1.1488960   1.181732

model1.2$waic$waic - model1.1$waic$waic

[1] -31.47579

The age and sex effect remains similar, but the model fit as measured by the WAIC is improved now that we account for local differences.

We find noisy estimates in some places, suggesting issues related to small area estimation. One solution is to partially pool information between municipalities that are geographically linked based on a spatial structure.

Model 1.3: structured spatial variability

We still focus on spatial variability, but now the municipalities are no longer independent: we account for the correlation between neighboring municipalities with a BYM model. Neighboring municipalities are defined as municipalities sharing a border. This will allow us to differentiate between what can be attributed to a municipality in particular, and what can be attributed to regional effects (like a COVID wave).

model1.3 = INLA::inla(munici_observed ~ - 1 + offset(E) +
                        sex:age_group +
                        f(id_space, model = "bym2", graph = "data/nb/gg_wm_q.adj", scale.model = TRUE, 
                          hyper = hyper.bym2, constr=TRUE),
                      data = data1,
                      family = "Poisson",
                      control.compute = list(config = TRUE, waic = TRUE),
                      quantiles = c(0.025, 0.5, 0.975),
                      num.threads = threads,
                      safe = TRUE)
summary(model1.3)


Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, 
   quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, 
   Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = 
   link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, 
   control.compute = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, control.fixed = 
   control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " 
   control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = 
   control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg 
   = inla.arg, num.threads = num.threads, ", " blas.num.threads = blas.num.threads, keep = 
   keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, 
   safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
Time used:
    Pre = 8.84, Running = 6.17, Post = 0.451, Total = 15.5 
Fixed effects:
                          mean    sd 0.025quant 0.5quant 0.975quant  mode kld
sexFemale:age_group40-59 0.015 0.024     -0.032    0.015      0.063 0.015   0
sexMale:age_group40-59   0.362 0.018      0.326    0.362      0.398 0.362   0
sexFemale:age_group60-69 0.044 0.020      0.005    0.044      0.083 0.044   0
sexMale:age_group60-69   0.427 0.015      0.398    0.427      0.457 0.427   0
sexFemale:age_group70-79 0.243 0.013      0.218    0.243      0.269 0.243   0
sexMale:age_group70-79   0.482 0.011      0.461    0.482      0.503 0.482   0
sexFemale:age_group80+   0.502 0.006      0.489    0.502      0.514 0.502   0
sexMale:age_group80+     0.158 0.007      0.144    0.158      0.172 0.158   0

Random effects:
  Name    Model
    id_space BYM2 model

Model hyperparameters:
                          mean      sd 0.025quant 0.5quant 0.975quant    mode
Precision for id_space 708.735 153.664    454.100  692.633   1056.562 661.505
Phi for id_space         0.976   0.033      0.885    0.987      0.999   0.997

Watanabe-Akaike information criterion (WAIC) ...: 46981.09
Effective number of parameters .................: 13.86

Marginal log-Likelihood:  -22734.63 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.3$summary.fixed)[c(1,3,5)]

                             mean 0.025quant 0.975quant
sexFemale:age_group40-59 1.015551  0.9683642   1.065039
sexMale:age_group40-59   1.436459  1.3853036   1.489507
sexFemale:age_group60-69 1.045067  1.0055022   1.086191
sexMale:age_group60-69   1.533087  1.4883961   1.579125
sexFemale:age_group70-79 1.275489  1.2437187   1.308077
sexMale:age_group70-79   1.619873  1.5860552   1.654422
sexFemale:age_group80+   1.651260  1.6306267   1.672193
sexMale:age_group80+     1.171041  1.1543556   1.187984

model1.3$waic$waic - model1.2$waic$waic

[1] -178.0986

We see that the structure accounts for a large part of the spatial variability (Phi estimated to 98%). This addition also improves the model fit as measured by the WAIC. The following interactive map allows to look at specific municipality effects.

We observe that many of the municipalities with higher relative excess mortality are in the western and southern parts, the ones that were hit first by COVID-19 in 2020. We also observe areas with higher excess in the North and Northeastern parts. These largely correspond to areas that were hit the most during the first and the second COVID-19 waves of spring and fall 2020 (Konstantinoudis et al. 2022).

Model 1.4: local characteristics

Having accounted for regional variability (arguably caused by COVID-19 waves of different timings and scales), we move on to explore the effect of local characteristics at the municipality level.

Rural/urban

The Federal Statistical Office classifies Swiss municipalities in 3 classes: urban, peri-urban or rural (https://www.bfs.admin.ch/bfs/en/home/statistics/territory-environment/nomenclatures/gemtyp.html). We add this covariate to the model taking the “rural” category as the reference.

model1.4a = INLA::inla(munici_observed ~ - 1 + offset(E) +
                         sex:age_group +
                         f(id_space, model = "bym2", graph = "data/nb/gg_wm_q.adj", scale.model = TRUE, 
                           hyper = hyper.bym2, constr=TRUE) +
                         type_periurban + type_urban,
                       data = data1,
                       family = "Poisson",
                       control.compute = list(config = TRUE, waic = TRUE),
                       quantiles = c(0.025, 0.5, 0.975),
                       num.threads = threads,
                       safe = TRUE)
summary(model1.4a)


Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, 
   quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, 
   Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = 
   link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, 
   control.compute = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, control.fixed = 
   control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " 
   control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = 
   control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg 
   = inla.arg, num.threads = num.threads, ", " blas.num.threads = blas.num.threads, keep = 
   keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, 
   safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
Time used:
    Pre = 9.14, Running = 3.19, Post = 0.481, Total = 12.8 
Fixed effects:
                           mean    sd 0.025quant 0.5quant 0.975quant   mode kld
type_periurban           -0.032 0.013     -0.057   -0.032     -0.007 -0.032   0
type_urban               -0.059 0.011     -0.080   -0.059     -0.037 -0.059   0
sexFemale:age_group40-59  0.055 0.026      0.005    0.055      0.106  0.055   0
sexMale:age_group40-59    0.402 0.020      0.363    0.402      0.442  0.402   0
sexFemale:age_group60-69  0.085 0.021      0.043    0.085      0.127  0.085   0
sexMale:age_group60-69    0.467 0.017      0.433    0.467      0.501  0.467   0
sexFemale:age_group70-79  0.284 0.015      0.254    0.284      0.314  0.284   0
sexMale:age_group70-79    0.523 0.014      0.496    0.523      0.550  0.523   0
sexFemale:age_group80+    0.543 0.011      0.522    0.543      0.564  0.543   0
sexMale:age_group80+      0.198 0.011      0.177    0.198      0.220  0.198   0

Random effects:
  Name    Model
    id_space BYM2 model

Model hyperparameters:
                         mean     sd 0.025quant 0.5quant 0.975quant   mode
Precision for id_space 772.16 170.72     493.03   752.89    1162.15 714.88
Phi for id_space         1.00   0.00       1.00     1.00       1.00   1.00

Watanabe-Akaike information criterion (WAIC) ...: 46955.41
Effective number of parameters .................: 12.93

Marginal log-Likelihood:  -22741.70 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.4a$summary.fixed)[c(1,3,5)]

                              mean 0.025quant 0.975quant
type_periurban           0.9684819  0.9446374  0.9929292
type_urban               0.9430452  0.9228926  0.9636496
sexFemale:age_group40-59 1.0568663  1.0049748  1.1114377
sexMale:age_group40-59   1.4953383  1.4369642  1.5560847
sexFemale:age_group60-69 1.0882576  1.0434695  1.1349688
sexMale:age_group60-69   1.5957264  1.5426347  1.6506476
sexFemale:age_group70-79 1.3285483  1.2890621  1.3692458
sexMale:age_group70-79   1.6872512  1.6424280  1.7333007
sexFemale:age_group80+   1.7214371  1.6856810  1.7579569
sexMale:age_group80+     1.2194751  1.1930976  1.2464389

model1.4a$waic$waic - model1.3$waic$waic

[1] -25.68285

drivers_plot(model1.4a,data1)

On average, peri-urban and urban municipalities appear to have a lower excess mortality than municipalities classified as rural.

Socio-economic position

The Swiss neighbourhood index of socio-economic position provides an estimate of socio-economic position (SEP) based on census data for 1.5 million buildings (Panczak et al. 2023). We consider the median index of each municipality, then group municipalities in quintiles before adding to the model (reference is 5th quintile with highest SEP).

model1.4b = INLA::inla(munici_observed ~ - 1 + offset(E) +
                         sex:age_group +
                         f(id_space, model = "bym2", graph = "data/nb/gg_wm_q.adj", scale.model = TRUE, 
                           hyper = hyper.bym2, constr=TRUE) +
                         sep1 + sep2 + sep3 + sep4,
                       data = data1,
                       family = "Poisson",
                       control.compute = list(config = TRUE, waic = TRUE),
                       quantiles = c(0.025, 0.5, 0.975),
                       num.threads = threads,
                       safe = TRUE)
summary(model1.4b)


Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, 
   quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, 
   Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = 
   link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, 
   control.compute = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, control.fixed = 
   control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " 
   control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = 
   control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg 
   = inla.arg, num.threads = num.threads, ", " blas.num.threads = blas.num.threads, keep = 
   keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, 
   safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
Time used:
    Pre = 8.73, Running = 4.36, Post = 0.456, Total = 13.5 
Fixed effects:
                           mean    sd 0.025quant 0.5quant 0.975quant   mode kld
sep1                      0.072 0.016      0.040    0.072      0.103  0.072   0
sep2                      0.034 0.013      0.008    0.034      0.060  0.034   0
sep3                      0.025 0.012      0.000    0.025      0.049  0.025   0
sep4                      0.002 0.011     -0.020    0.002      0.024  0.002   0
sexFemale:age_group40-59 -0.006 0.025     -0.056   -0.006      0.044 -0.006   0
sexMale:age_group40-59    0.340 0.020      0.301    0.340      0.379  0.340   0
sexFemale:age_group60-69  0.023 0.021     -0.018    0.023      0.064  0.023   0
sexMale:age_group60-69    0.405 0.017      0.372    0.405      0.439  0.405   0
sexFemale:age_group70-79  0.222 0.015      0.192    0.222      0.251  0.222   0
sexMale:age_group70-79    0.461 0.013      0.435    0.461      0.487  0.461   0
sexFemale:age_group80+    0.480 0.010      0.461    0.480      0.500  0.480   0
sexMale:age_group80+      0.136 0.010      0.116    0.136      0.157  0.136   0

Random effects:
  Name    Model
    id_space BYM2 model

Model hyperparameters:
                          mean      sd 0.025quant 0.5quant 0.975quant    mode
Precision for id_space 955.536 242.463    577.544  921.998   1525.727 855.831
Phi for id_space         0.961   0.039      0.854    0.973      0.998   0.993

Watanabe-Akaike information criterion (WAIC) ...: 46974.24
Effective number of parameters .................: 12.71

Marginal log-Likelihood:  -22754.14 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.4b$summary.fixed)[c(1,3,5)]

                              mean 0.025quant 0.975quant
sep1                     1.0746661  1.0412310   1.108975
sep2                     1.0343216  1.0078085   1.061363
sep3                     1.0251608  1.0003805   1.050414
sep4                     1.0020228  0.9801422   1.024357
sexFemale:age_group40-59 0.9940482  0.9456350   1.044964
sexMale:age_group40-59   1.4051921  1.3511187   1.461484
sexFemale:age_group60-69 1.0232876  0.9818185   1.066542
sexMale:age_group60-69   1.4999336  1.4509266   1.550677
sexFemale:age_group70-79 1.2480166  1.2121342   1.285045
sexMale:age_group70-79   1.5859160  1.5456989   1.627313
sexFemale:age_group80+   1.6167121  1.5860118   1.648265
sexMale:age_group80+     1.1462144  1.1230406   1.170020

model1.4b$waic$waic - model1.3$waic$waic

[1] -6.858679

drivers_plot(model1.4b,data1)

A gradient appears clearly, so that we can conclude that municipalities of lowest median SEP had a higher relative excess mortality in 2020 compared to municipalities of highest median SEP (5th quintile).

International borders

We now consider whether the municipality belongs to a cross-border labor region (https://www.bfs.admin.ch/bfs/de/home/grundlagen/raumgliederungen.assetdetail.8706500.html). This identifies municipalities with a high level of connectivity with neighboring countries France, Germany and Italy.

model1.4c = INLA::inla(munici_observed ~ - 1 + offset(E) +
                         sex:age_group +
                         f(id_space, model = "bym2", graph = "data/nb/gg_wm_q.adj", scale.model = TRUE, 
                           hyper = hyper.bym2, constr=TRUE) +
                         border,
                       data = data1,
                       family = "Poisson",
                       control.compute = list(config = TRUE, waic = TRUE),
                       quantiles = c(0.025, 0.5, 0.975),
                       num.threads = threads,
                       safe = TRUE)
summary(model1.4c)


Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, 
   quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, 
   Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = 
   link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, 
   control.compute = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, control.fixed = 
   control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " 
   control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = 
   control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg 
   = inla.arg, num.threads = num.threads, ", " blas.num.threads = blas.num.threads, keep = 
   keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, 
   safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
Time used:
    Pre = 8.93, Running = 3.34, Post = 0.471, Total = 12.7 
Fixed effects:
                          mean    sd 0.025quant 0.5quant 0.975quant  mode kld
border                   0.041 0.015      0.011    0.041      0.070 0.041   0
sexFemale:age_group40-59 0.008 0.024     -0.040    0.008      0.056 0.008   0
sexMale:age_group40-59   0.355 0.019      0.319    0.355      0.392 0.355   0
sexFemale:age_group60-69 0.037 0.020     -0.002    0.037      0.076 0.037   0
sexMale:age_group60-69   0.420 0.015      0.390    0.420      0.450 0.420   0
sexFemale:age_group70-79 0.236 0.013      0.210    0.236      0.262 0.236   0
sexMale:age_group70-79   0.475 0.011      0.453    0.475      0.497 0.475   0
sexFemale:age_group80+   0.494 0.007      0.480    0.494      0.508 0.494   0
sexMale:age_group80+     0.150 0.008      0.135    0.150      0.166 0.150   0

Random effects:
  Name    Model
    id_space BYM2 model

Model hyperparameters:
                          mean      sd 0.025quant 0.5quant 0.975quant    mode
Precision for id_space 733.110 161.170    472.590  713.787   1104.173 674.729
Phi for id_space         0.965   0.035      0.871    0.976      0.998   0.994

Watanabe-Akaike information criterion (WAIC) ...: 46977.16
Effective number of parameters .................: 14.46

Marginal log-Likelihood:  -22738.69 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.4c$summary.fixed)[c(1,3,5)]

                             mean 0.025quant 0.975quant
border                   1.041472  1.0115449   1.072380
sexFemale:age_group40-59 1.007930  0.9608114   1.057363
sexMale:age_group40-59   1.426601  1.3753495   1.479769
sexFemale:age_group60-69 1.037441  0.9978187   1.078642
sexMale:age_group60-69   1.522199  1.4771990   1.568580
sexFemale:age_group70-79 1.266075  1.2338749   1.299129
sexMale:age_group70-79   1.607786  1.5731915   1.643163
sexFemale:age_group80+   1.639002  1.6168112   1.661561
sexMale:age_group80+     1.162313  1.1446724   1.180257

model1.4c$waic$waic - model1.3$waic$waic

[1] -3.930987

drivers_plot(model1.4c,data1)

We observe higher relative excess mortality in municipalities in the vicinity of international borders.

Language

Most Swiss municipalities have one official language: German, French or Italian. A few municipalities have several official languages, but given the relatively low numbers, we consider only the majority language. The difficulty is the colinearity between language regions and the first COVID-19 wave of 2020, that primarily affected Ticino (Italian) and Southwestern Switzerland (French), mostly because of how the initial global spread of COVID-19 occurred (with large early epidemics in Italy and then France). These effects are much larger than any effect that could be attributed to cultural differences between language regions, so it is quite difficult to estimate the latter. We still attempt to do so by adding the language of each municipality (reference is German) to our model.

model1.4d = INLA::inla(munici_observed ~ - 1 + offset(E) +
                         sex:age_group +
                         f(id_space, model = "bym2", graph = "data/nb/gg_wm_q.adj", scale.model = TRUE, 
                           hyper = hyper.bym2, constr=TRUE) +
                         lang_fr + lang_it,
                       data = data1,
                       family = "Poisson",
                       control.compute = list(config = TRUE, waic = TRUE),
                       quantiles = c(0.025, 0.5, 0.975),
                       num.threads = threads,
                       safe = TRUE)
summary(model1.4d)


Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, 
   quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, 
   Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = 
   link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, 
   control.compute = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, control.fixed = 
   control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " 
   control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = 
   control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg 
   = inla.arg, num.threads = num.threads, ", " blas.num.threads = blas.num.threads, keep = 
   keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, 
   safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
Time used:
    Pre = 9.25, Running = 3.93, Post = 0.47, Total = 13.6 
Fixed effects:
                           mean    sd 0.025quant 0.5quant 0.975quant   mode kld
lang_fr                   0.090 0.013      0.063    0.090      0.116  0.091   0
lang_it                   0.153 0.022      0.110    0.154      0.196  0.154   0
sexFemale:age_group40-59 -0.021 0.025     -0.069   -0.021      0.028 -0.021   0
sexMale:age_group40-59    0.326 0.019      0.289    0.326      0.363  0.326   0
sexFemale:age_group60-69  0.007 0.020     -0.032    0.007      0.047  0.007   0
sexMale:age_group60-69    0.392 0.016      0.361    0.392      0.422  0.392   0
sexFemale:age_group70-79  0.206 0.013      0.180    0.206      0.233  0.206   0
sexMale:age_group70-79    0.446 0.011      0.423    0.446      0.468  0.446   0
sexFemale:age_group80+    0.464 0.007      0.449    0.464      0.479  0.464   0
sexMale:age_group80+      0.121 0.008      0.105    0.121      0.137  0.121   0

Random effects:
  Name    Model
    id_space BYM2 model

Model hyperparameters:
                           mean      sd 0.025quant 0.5quant 0.975quant     mode
Precision for id_space 3107.867 1870.36   1064.845 2618.835   8051.973 1937.674
Phi for id_space          0.856    0.14      0.463    0.903      0.992    0.981

Watanabe-Akaike information criterion (WAIC) ...: 46967.71
Effective number of parameters .................: 8.91

Marginal log-Likelihood:  -22719.47 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.4d$summary.fixed)[c(1,3,5)]

                             mean 0.025quant 0.975quant
lang_fr                  1.094028  1.0650015   1.122499
lang_it                  1.165839  1.1158181   1.216709
sexFemale:age_group40-59 0.979651  0.9336435   1.027950
sexMale:age_group40-59   1.385655  1.3353479   1.437917
sexFemale:age_group60-69 1.007406  0.9685883   1.047820
sexMale:age_group60-69   1.479263  1.4349687   1.525020
sexFemale:age_group70-79 1.228978  1.1970985   1.261826
sexMale:age_group70-79   1.561569  1.5271499   1.596961
sexFemale:age_group80+   1.590278  1.5673892   1.614042
sexMale:age_group80+     1.128295  1.1103044   1.146865

model1.4d$waic$waic - model1.3$waic$waic

[1] -13.38201

drivers_plot(model1.4d,data1) + theme(legend.position=c(.3,.87))

At first sight, we may thing that there is a large effect of language region on excess mortality, with around 7-12% more deaths than expected in French-speaking municipalities and 12-22% more in Italian-speaking municipalities compared to German. However, as expected this association is likely confounded by the regional variability associated with COVID-19 waves in 2020. Indeed, if we now look at the geographically-structured municipality effect for this model, which can be interpreted as residual effects, we see that the higher excess in South and Southwestern Switzerland is now more evenly distributed (captured by the language effect), while French-speaking regions that were comparatively less impacted during the first wave (such as Neuchâtel and Jura) now have a negative municipality effect to compensate. These nonsensical results highlight the difficulty to estimate the effect of language regions. For this reason, in the following we rely upon observing the residual municipality effects to draw conclusion about the association with language rather than using the language as a fixed effect, as shown in the next map based on model 1.3.

In this last map, we can make two observations. First, French-speaking and Italian-speaking municipalities are not systematically more affected by excess mortality than German-speaking municipalities, with exceptions like the area around Neuchâtel and the Italian-speaking municipalities of Graubunden. Second, there is a clear separation between the French- and German-speaking municipalities, so that differences could rather be attributed to a lower level of connectivity between populations of different language rather than intrinsic differences favoring SARS-CoV-2 transmission or severity.

Referendums on COVID-19 measures

We now focus on results from two referendums about COVID-19 control measures held in June and November 2020. The point here is not to look at causality one way or the other, as we look at overall excess for 2020, and the voting took place at two separated points. A preliminary analysis has reported a negative association between the proportion of “yes” vote at the November referendum at the cantonal level and 7-day incidence on December 7, 2021 (https://smw.ch/index.php/smw/announcement/view/50). We classify municipalities according to the proportion of “yes” vote (expressing support of government-issued measures aimed at controlling COVID-19) at each vote, in quintiles (taking the 5th quintile - with highest support - as a reference).

model1.4e = INLA::inla(munici_observed ~ - 1 + offset(E) +
                         sex:age_group +
                         f(id_space, model = "bym2", graph = "data/nb/gg_wm_q.adj", scale.model = TRUE, 
                           hyper = hyper.bym2, constr=TRUE) +
                         vote_jun_q1 + vote_jun_q2 + vote_jun_q3 + vote_jun_q4,
                       data = data1,
                       family = "Poisson",
                       control.compute = list(config = TRUE, waic = TRUE),
                       quantiles = c(0.025, 0.5, 0.975),
                       num.threads = threads,
                       safe = TRUE)
summary(model1.4e)


Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, 
   quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, 
   Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = 
   link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, 
   control.compute = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, control.fixed = 
   control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " 
   control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = 
   control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg 
   = inla.arg, num.threads = num.threads, ", " blas.num.threads = blas.num.threads, keep = 
   keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, 
   safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
Time used:
    Pre = 9.32, Running = 6.54, Post = 0.459, Total = 16.3 
Fixed effects:
                          mean    sd 0.025quant 0.5quant 0.975quant  mode kld
vote_jun_q1              0.045 0.015      0.015    0.045      0.075 0.045   0
vote_jun_q2              0.033 0.014      0.005    0.033      0.061 0.033   0
vote_jun_q3              0.021 0.013     -0.003    0.021      0.046 0.021   0
vote_jun_q4              0.011 0.011     -0.011    0.011      0.032 0.011   0
sexFemale:age_group40-59 0.000 0.025     -0.049    0.000      0.049 0.000   0
sexMale:age_group40-59   0.347 0.019      0.309    0.347      0.385 0.347   0
sexFemale:age_group60-69 0.029 0.021     -0.012    0.029      0.069 0.029   0
sexMale:age_group60-69   0.412 0.016      0.380    0.412      0.444 0.412   0
sexFemale:age_group70-79 0.229 0.014      0.201    0.229      0.256 0.229   0
sexMale:age_group70-79   0.467 0.012      0.443    0.467      0.492 0.467   0
sexFemale:age_group80+   0.487 0.009      0.471    0.487      0.504 0.487   0
sexMale:age_group80+     0.143 0.009      0.125    0.143      0.162 0.143   0

Random effects:
  Name    Model
    id_space BYM2 model

Model hyperparameters:
                          mean      sd 0.025quant 0.5quant 0.975quant    mode
Precision for id_space 716.412 154.948    467.669  697.183   1074.761 657.558
Phi for id_space         0.979   0.024      0.911    0.986      0.999   0.997

Watanabe-Akaike information criterion (WAIC) ...: 46974.92
Effective number of parameters .................: 14.05

Marginal log-Likelihood:  -22760.47 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.4e$summary.fixed)[c(1,3,5)]

                             mean 0.025quant 0.975quant
vote_jun_q1              1.045806  1.0148236   1.077765
vote_jun_q2              1.033527  1.0054030   1.062455
vote_jun_q3              1.021409  0.9966143   1.046817
vote_jun_q4              1.010717  0.9889794   1.032934
sexFemale:age_group40-59 1.000141  0.9522514   1.050445
sexMale:age_group40-59   1.414452  1.3615214   1.469455
sexFemale:age_group60-69 1.029322  0.9885622   1.071772
sexMale:age_group60-69   1.509576  1.4621740   1.558542
sexFemale:age_group70-79 1.256899  1.2225230   1.292268
sexMale:age_group70-79   1.595678  1.5575369   1.634799
sexFemale:age_group80+   1.628040  1.6009003   1.655742
sexMale:age_group80+     1.154025  1.1330425   1.175452

model1.4e$waic$waic - model1.3$waic$waic

[1] -6.174953

drivers_plot(model1.4e,data1)

model1.4f = INLA::inla(munici_observed ~ - 1 + offset(E) +
                         sex:age_group +
                         f(id_space, model = "bym2", graph = "data/nb/gg_wm_q.adj", scale.model = TRUE, 
                           hyper = hyper.bym2, constr=TRUE) +
                         vote_nov_q1 + vote_nov_q2 + vote_nov_q3 + vote_nov_q4,
                       data = data1,
                       family = "Poisson",
                       control.compute = list(config = TRUE, waic = TRUE),
                       quantiles = c(0.025, 0.5, 0.975),
                       num.threads = threads,
                       safe = TRUE)
summary(model1.4f)


Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, 
   quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, 
   Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = 
   link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, 
   control.compute = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, control.fixed = 
   control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " 
   control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = 
   control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg 
   = inla.arg, num.threads = num.threads, ", " blas.num.threads = blas.num.threads, keep = 
   keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, 
   safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
Time used:
    Pre = 9.22, Running = 5.14, Post = 0.488, Total = 14.8 
Fixed effects:
                           mean    sd 0.025quant 0.5quant 0.975quant   mode kld
vote_nov_q1               0.044 0.015      0.015    0.044      0.074  0.044   0
vote_nov_q2               0.043 0.013      0.017    0.043      0.070  0.043   0
vote_nov_q3               0.022 0.012      0.000    0.022      0.045  0.023   0
vote_nov_q4               0.017 0.011     -0.005    0.017      0.038  0.017   0
sexFemale:age_group40-59 -0.004 0.025     -0.053   -0.004      0.045 -0.004   0
sexMale:age_group40-59    0.343 0.019      0.305    0.343      0.381  0.343   0
sexFemale:age_group60-69  0.025 0.021     -0.015    0.025      0.065  0.025   0
sexMale:age_group60-69    0.408 0.016      0.376    0.408      0.440  0.408   0
sexFemale:age_group70-79  0.225 0.014      0.197    0.225      0.253  0.225   0
sexMale:age_group70-79    0.464 0.012      0.440    0.464      0.488  0.464   0
sexFemale:age_group80+    0.484 0.008      0.467    0.484      0.500  0.484   0
sexMale:age_group80+      0.140 0.009      0.122    0.140      0.158  0.140   0

Random effects:
  Name    Model
    id_space BYM2 model

Model hyperparameters:
                          mean      sd 0.025quant 0.5quant 0.975quant    mode
Precision for id_space 779.888 175.200    497.667  758.586   1184.072 715.767
Phi for id_space         0.968   0.031      0.883    0.978      0.998   0.995

Watanabe-Akaike information criterion (WAIC) ...: 46974.75
Effective number of parameters .................: 13.60

Marginal log-Likelihood:  -22758.78 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.4f$summary.fixed)[c(1,3,5)]

                              mean 0.025quant 0.975quant
vote_nov_q1              1.0454678  1.0149771   1.076855
vote_nov_q2              1.0444218  1.0174558   1.072064
vote_nov_q3              1.0227422  0.9996515   1.046301
vote_nov_q4              1.0167237  0.9950556   1.038837
sexFemale:age_group40-59 0.9961336  0.9485066   1.046162
sexMale:age_group40-59   1.4091147  1.3564863   1.463808
sexFemale:age_group60-69 1.0254934  0.9850296   1.067633
sexMale:age_group60-69   1.5040277  1.4570849   1.552518
sexFemale:age_group70-79 1.2522786  1.2182269   1.287318
sexMale:age_group70-79   1.5897221  1.5521168   1.628301
sexFemale:age_group80+   1.6222125  1.5956871   1.649315
sexMale:age_group80+     1.1498773  1.1293344   1.170870

model1.4f$waic$waic - model1.3$waic$waic

[1] -6.347276

drivers_plot(model1.4f,data1)

In both cases it appears that excess mortality is consistently about 5% higher in municipalities expressing lowest support to control measures (first quantile) in both referendums, with a clear gradient.

Multivariable model

We now jointly estimate the effects of interest identified in the univariable analysis: rural or urban status, border, median SEP quintile and results from the COVID-19 referendums (using only the June referendum to limit complexity). We leave out language regions for the reasons explained above.

model1.5 = INLA::inla(munici_observed ~ - 1 + offset(E) +
                        sex:age_group +
                        f(id_space, model = "bym2", graph = "data/nb/gg_wm_q.adj", scale.model = TRUE, 
                          hyper = hyper.bym2, constr=TRUE) +
                        border +
                        type_periurban + type_urban +
                        sep1 + sep2 + sep3 + sep4 +
                        vote_jun_q1 + vote_jun_q2 + vote_jun_q3 + vote_jun_q4,
                      data = data1,
                      family = "Poisson",
                      control.compute = list(config = TRUE, waic = TRUE),
                      quantiles = c(0.025, 0.5, 0.975),
                      num.threads = threads,
                      safe = TRUE)
summary(model1.5)


Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, 
   quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, 
   Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = 
   link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, 
   control.compute = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, control.fixed = 
   control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " 
   control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = 
   control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg 
   = inla.arg, num.threads = num.threads, ", " blas.num.threads = blas.num.threads, keep = 
   keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, 
   safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
Time used:
    Pre = 9.66, Running = 4.71, Post = 0.544, Total = 14.9 
Fixed effects:
                           mean    sd 0.025quant 0.5quant 0.975quant   mode kld
border                    0.043 0.014      0.016    0.043      0.071  0.043   0
type_periurban           -0.023 0.013     -0.049   -0.023      0.003 -0.023   0
type_urban               -0.048 0.013     -0.073   -0.048     -0.022 -0.048   0
sep1                      0.056 0.018      0.021    0.056      0.090  0.056   0
sep2                      0.026 0.014     -0.003    0.026      0.054  0.026   0
sep3                      0.019 0.013     -0.007    0.019      0.045  0.019   0
sep4                      0.000 0.011     -0.022    0.000      0.022  0.000   0
vote_jun_q1              -0.008 0.018     -0.044   -0.008      0.028 -0.008   0
vote_jun_q2              -0.007 0.016     -0.039   -0.008      0.024 -0.008   0
vote_jun_q3              -0.005 0.014     -0.032   -0.005      0.022 -0.005   0
vote_jun_q4              -0.005 0.011     -0.027   -0.005      0.017 -0.005   0
sexFemale:age_group40-59  0.027 0.029     -0.029    0.027      0.084  0.027   0
sexMale:age_group40-59    0.374 0.024      0.327    0.374      0.422  0.374   0
sexFemale:age_group60-69  0.057 0.025      0.008    0.057      0.106  0.057   0
sexMale:age_group60-69    0.439 0.022      0.397    0.439      0.482  0.439   0
sexFemale:age_group70-79  0.256 0.020      0.216    0.256      0.295  0.256   0
sexMale:age_group70-79    0.495 0.019      0.458    0.495      0.532  0.495   0
sexFemale:age_group80+    0.515 0.017      0.482    0.515      0.547  0.515   0
sexMale:age_group80+      0.170 0.017      0.137    0.170      0.203  0.170   0

Random effects:
  Name    Model
    id_space BYM2 model

Model hyperparameters:
                         mean     sd 0.025quant 0.5quant 0.975quant    mode
Precision for id_space 998.47 256.76     599.66  962.499   1603.503 891.712
Phi for id_space         0.96   0.04       0.85    0.972      0.997   0.993

Watanabe-Akaike information criterion (WAIC) ...: 46958.35
Effective number of parameters .................: 13.80

Marginal log-Likelihood:  -22796.10 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

exp(model1.5$summary.fixed)[c(1,3,5)]

                              mean 0.025quant 0.975quant
border                   1.0443881  1.0156991  1.0739306
type_periurban           0.9771043  0.9519695  1.0028949
type_urban               0.9534305  0.9294926  0.9779838
sep1                     1.0572317  1.0210039  1.0945124
sep2                     1.0258589  0.9970622  1.0552796
sep3                     1.0193663  0.9931238  1.0461396
sep4                     0.9998385  0.9778917  1.0222444
vote_jun_q1              0.9919884  0.9569852  1.0284108
vote_jun_q2              0.9925411  0.9616492  1.0245207
vote_jun_q3              0.9947892  0.9686157  1.0217391
vote_jun_q4              0.9950804  0.9732749  1.0174220
sexFemale:age_group40-59 1.0276277  0.9713554  1.0871858
sexMale:age_group40-59   1.4542422  1.3873319  1.5244301
sexFemale:age_group60-69 1.0585346  1.0079405  1.1117016
sexMale:age_group60-69   1.5515998  1.4874597  1.6185796
sexFemale:age_group70-79 1.2911573  1.2416425  1.3427143
sexMale:age_group70-79   1.6404637  1.5812474  1.7019932
sexFemale:age_group80+   1.6732965  1.6199343  1.7285486
sexMale:age_group80+     1.1853712  1.1468813  1.2252393

model1.5$waic$waic - model1.3$waic$waic

[1] -22.74136

drivers_plot(model1.5,data1)

This multivariate analysis confirms the association between high relative excess mortality and low median SEP, with a clear gradient. We also see that the association between high excess mortality and border municipalities remains, and than urban and peri-urban areas are still associated with comparatively higher excess mortality in 2020, although uncertainty increases. Estimates of an association with voting results are faint after adjusting for the other covariates. This decrease in the association can be attributed to the collinearity with other covariates, in particular with SEP and urbanicity (see correlogram in section 1). Therefore, we cannot disentangle between the associations with voting results, SEP and urbanicity.

There are also some interesting patterns in the residual effects at the level of the municipality (adjusting for all aforementioned covariates), with in particular, expected higher excesses in Ticino and Southwestern Switzerland, a visible barrier between French-speaking and German-speaking regions, lower excess in the large cities of the German-speaking part (Zurich, Basel, Bern) and in relatively isolated valleys of Graubunden.

save(list=ls(pattern="model1."),data1,
     file = "results_inla/local_corr_models.Rdata")

Step 2: multivariate model with full uncertainty propagation

In the previous section we modeled the variation of the median excess mortality over 2020 by municipality. This approach underestimates the uncertainty from two sources, first from the prediction error in the expected mortality at the cantonal level, second from the downscaling to the municipal level. At this stage, we bring back these two sources of uncertainty in the final estimates by repeatedly fitting model 1.5 to 50 different sets of posterior draws of excess mortality by municipality, then combining with equal weights.

model1.5_merg = readRDS("results_inla/model1.5_merg.rds")
summary(model1.5_merg)


Call:
   c("inla.core(formula = formula, family = family, contrasts = contrasts, ", " data = data, 
   quantiles = quantiles, E = E, offset = offset, ", " scale = scale, weights = weights, 
   Ntrials = Ntrials, strata = strata, ", " lp.scale = lp.scale, link.covariates = 
   link.covariates, verbose = verbose, ", " lincomb = lincomb, selection = selection, 
   control.compute = control.compute, ", " control.predictor = control.predictor, 
   control.family = control.family, ", " control.inla = control.inla, control.fixed = 
   control.fixed, ", " control.mode = control.mode, control.expert = control.expert, ", " 
   control.hazard = control.hazard, control.lincomb = control.lincomb, ", " control.update = 
   control.update, control.lp.scale = control.lp.scale, ", " control.pardiso = 
   control.pardiso, only.hyperparam = only.hyperparam, ", " inla.call = inla.call, inla.arg 
   = inla.arg, num.threads = num.threads, ", " blas.num.threads = blas.num.threads, keep = 
   keep, working.directory = working.directory, ", " silent = silent, inla.mode = inla.mode, 
   safe = FALSE, debug = debug, ", " .parent.frame = .parent.frame)") 
Time used:
    Pre = 568, Running = 147, Post = 34.6, Total = 750 
Fixed effects:
                           mean    sd
border                    0.038 0.029
type_periurban            0.003 0.019
type_urban                0.009 0.020
sep1                      0.041 0.031
sep2                      0.023 0.025
sep3                      0.018 0.025
sep4                      0.004 0.022
vote_jun_q1              -0.010 0.032
vote_jun_q2              -0.013 0.030
vote_jun_q3              -0.008 0.025
vote_jun_q4              -0.008 0.024
sexFemale:age_group40-59 -0.044 0.050
sexMale:age_group40-59    0.274 0.053
sexFemale:age_group60-69 -0.013 0.045
sexMale:age_group60-69    0.326 0.047
sexFemale:age_group70-79  0.201 0.052
sexMale:age_group70-79    0.396 0.044
sexFemale:age_group80+    0.451 0.041
sexMale:age_group80+      0.125 0.044

Random effects:
  Name    Model
    id_space BYM2 model

Model hyperparameters:
                         mean     sd
Precision for id_space 93.520 24.652
Phi for id_space        0.433  0.176

Marginal log-Likelihood:  -45658.28 
 is computed 
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

drivers_plot(model1.5_merg,data1,compute_cri = TRUE)

As expected, this approach leads to a dilution of the observed associations between relative excess mortality and local covariates. We still observe a linear gradient in the association between excess mortality and median SEP at the municipal level in age groups 40 to 79, and a likely association between excess mortality and border municipalities, although in both cases with higher uncertainty. We also observe similar patterns in the residual effects at the level of the municipality, again with higher uncertainty.

References

Konstantinoudis, Garyfallos, Michela Cameletti, Virgilio Gómez-Rubio, Inmaculada León Gómez, Monica Pirani, Gianluca Baio, Amparo Larrauri, et al. 2022. “Regional Excess Mortality During the 2020 COVID-19 Pandemic in Five European Countries.” Nature Communications 13 (1): 482.

Panczak, Radoslaw, Claudia Berlin, Marieke Voorpostel, Marcel Zwahlen, and Matthias Egger. 2023. “The Swiss Neighbourhood Index of Socioeconomic Position: Update and Re-Validation.” Swiss Medical Weekly 153 (40028): 40028.