A linear mixed model to estimate COVID-19-induced excess mortality.

The Corona Virus Disease (COVID-19) pandemic has increased mortality in countries worldwide. To evaluate the impact of the pandemic on mortality, the use of excess mortality rather than reported COVID-19 deaths has been suggested. Excess mortality, however, requires estimation of mortality under nonpandemic conditions. Although many methods exist to forecast mortality, they are either complex to apply, require many sources of information, ignore serial correlation, and/or are influenced by historical excess mortality. We propose a linear mixed model that is easy to apply, requires only historical mortality data, allows for serial correlation, and down-weighs the influence of historical excess mortality. Appropriateness of the linear mixed model is evaluated with fit statistics and forecasting accuracy measures for Belgium and the Netherlands. Unlike the commonly used 5-year weekly average, the linear mixed model is forecasting the year-specific mortality, and as a result improves the estimation of excess mortality for Belgium and the Netherlands.

INTRODUCTION

During the Corona Virus Disease (COVID‐19) pandemic, most countries have reported the number of COVID‐19 deaths as an essential part of their monitoring strategy (Giattino et al., 2021). However, reported COVID‐19 deaths depend on the completeness and strategy of counting deaths. Variation in counting exists through testing strategy, availability of test material, in‐ or excluding nursing home deaths, or by variation in coding and registration. Hence, reported COVID‐19 mortality is prone to mis‐reporting. Therefore, excess mortality has been suggested to assess the overall impact on mortality of the SARS‐CoV‐2 virus (Aron et al., 2020; Morgan et al., 2020; Beaney et al., 2020). Excess mortality is obtained by subtracting the expected deaths based on the pre‐pandemic period from the registered all‐cause deaths in the pandemic period. Reported all‐cause deaths are not only more reliable across countries, but excess death captures also both direct and indirect effects of this pandemic on mortality. Lower death counts due to mitigation measures or higher counts due to COVID‐19 or an overloaded health system will both be reflected in excess mortality.

The critical part in determining excess mortality is a reliable estimate of baseline mortality, that is, the mortality that is expected under nonpandemic conditions. A simple method to determine baseline mortality is the average mortality based on historical data, most commonly the weekly average of the past 5 years (Beaney et al., 2020; Modig et al., 2020; Stang et al., 2020; Giattino et al., 2021; Michelozzi et al., 2020). However, this methodology treats mortality in subsequent weeks as independent observations. Moreover, by predicting the average, this methodology ignores year‐specific trends in mortality. For example, in Belgium, mortality was below average during the first weeks of 2020, that is, right before the start of the pandemic, due to a mild influenza season (Molenberghs et al., 2020; Bustos Sierra et al., 2020). Additionally, the weekly average may be influenced by peaks of increased mortality due to heat waves or seasonal influenza epidemics in recent history. Two popular methods to minimize past influence of excess mortality in the forecasting of baseline mortality are the Farrington (Farrington et al., 1996) and Euro‐MOMO models (Vestergaard et al., 2020; Fouillet et al., 2020). The Farrington model uses residuals to down‐weigh the influence of outbreaks in the past, while the Euro‐MOMO model takes only historical periods without excess mortality into account to forecast future baseline mortality. However, excluding Winter and Summer seasons due to influenza or heat waves may not be sufficient to eliminate the influence of these events on mortality (Aron et al., 2020). The number of deaths in Spring may, for example, be below average after a severe seasonal influenza season.

Variations of the Euro‐MOMO model exist that do not exclude historical data for mortality forecasting and add, similarly as for Serfling's models (Serfling, 1963), a cyclic term to model seasonality (Cox et al., 2010; Nielsen et al., 2018, 2021). However, these extensions often require additional information, such as historical influenza data, temperature, and ozone concentration, which may not always be easy to access, especially when several countries are considered in a common analysis. Time series models, such as Dynamic Harmonic Regression (Chen et al., 2020) and ARIMA models (Faust et al., 2021), exploit the serial correlation in the historical mortality data. The latter models, however, require stationarity of the time series and are well suited for short‐term forecasting, but may be of limited use in the long run (Harvey, 1989; Harvey and Todd, 1983). Others have suggested simple linear models to determine COVID‐19 excess mortality, including a yearly time trend in combination with fixed or smoothing‐spline‐based weekly effects (The Economist, 2020; The New York Times, 2020; Karlinsky and Kobak, 2021). That said, the longitudinal mortality data likely violate the independent error assumption of such linear models.

We propose a linear mixed model based on uninterrupted historical mortality data to forecast the year‐specific baseline mortality for the year 2020 in order to estimate the COVID‐19 excess mortality. Although linear mixed models offer a versatile modeling family, which can incorporate many mean and (co)variability structures for longitudinal data, including serial correlation (Verbeke and Molenberghs, 2000; Verbeke et al., 1998; Chi and Reinsel, 1989), we show that they are particularly well designed for modeling baseline mortality patterns. The influence of historical excess mortality is downsized by two distinct strategies, more specifically by down‐weighing the residuals, similar to Farrington et al. (1996), and by down‐weighing the historical excess mortality data. Although marginal population‐averaged predictions and conditional year‐specific predictions cannot be compared directly, we will demonstrate the advantage of using the year‐specific predictions of our proposed linear mixed model over the commonly used 5‐year weekly population average in forecasting the mortality in 2020 for Belgium and the Netherlands.

DATA

Open source daily all‐cause Belgian mortality data from the year 2009 until 2020 are available from the National Statistical Institute, Statistics Belgium (STATBEL, 2020). These data were downloaded on January 25, 2021 and temporally aggregated in weekly periods. The weeks are defined according to the International Standard ISO 8601 definition, that is, Monday is the first day of the week and the first week of the year is the week that contains the first Thursday of January. The first, incomplete, week of the year 2009 was excluded from the data. Also, the weeks numbered 53, present in some years, are excluded.

In Belgium, daily COVID‐19 mortality data are registered by the Belgian institute for public health, Sciensano (EPISTAT: COVID‐19, 2020). These open source data were extracted on January 26, 2021 and aggregated in weeks using the same week definition as for the all‐cause mortality. Registered COVID‐19 related deaths in Belgium include confirmed and possible COVID‐19 deaths (Bustos Sierra et al., 2020).

For the Netherlands, weekly historical mortality data are available from The Human Mortality Database (2021), reported by the National Statistical Office, Central Bureau for Statistics (2021b), for the year 1995 until 2020, while daily reported COVID‐19 mortality is available from the National Institute for Public Health and the Environment (National Institute for Public Health and the Environment, 2021). The same aggregation methods as used for the Belgian data are applied to the data, which was downloaded on April 7, 2021.

METHODOLOGY

To forecast the year‐specific mortality in the year 2020 based on historical mortality data, in Section 3.1 a general linear mixed model is presented, allowing for a possible temporal correlation structure and down‐weighting of past excess mortality, such as heat waves and influenza outbreaks. With these models, the excess mortality and its 95% prediction interval (PI) is estimated, as explained in Section 3.2. Finally, in Section 3.3, methods to compare several choices of serial temporal correlation structures and two down‐weighting methods for past excess mortality are presented.

The data analyses were performed and figures produced using SAS 9.4 Software and R Studio 4.0.3.

Model proposal for excess mortality

General linear mixed model

We model the weekly mortality with week by year . The number of deaths is usually modeled with a Poisson distribution. However, since the mean of the weekly deaths is sufficiently high for the central limit theorem to be invoked, we will use a Gaussian model. The Gaussian model comes with the additional advantage that it has less convergence issues and a more straightforward connection between marginal and conditional year‐specific interpretation.

Mortality often shows a cyclic pattern within a year, captured here by Fourier series. The number of Fourier terms m is determined via a correlogram of the historical mortality. As mortality may fluctuate year by year due to increasing population sizes or changing age distributions, a random intercept is added, resulting in the following model: with , , and and mutually independent.

The variation in the cyclic pattern from year to year can be modeled by including random effects of the Fourier terms, resulting in a random‐effects vector , with and mutually independent. The need to include additional random effects in the model is evaluated via likelihood ratio tests. Importantly, it is known that random effects are often able to represent the serial correlation among the measurements (Chi and Reinsel, 1989; Verbeke and Molenberghs, 2000), thus including random effects may be sufficient to also capture the serial correlation between mortality of consecutive weeks. Alternatively, in addition to the random effects, serial correlation may be introduced in the model. However, as this strategy sometimes overparameterizes the covariance structure, we will carefully evaluate whether splitting the error in a serial correlation and a measurement error is required. In any case, it has been shown that, if serial correlation is present in addition to the random‐effects correlation, the inclusion of a serial correlation structure is preferable over correctly specifying the model (Verbeke et al., 1998). Model parameters in our proposed model are estimated via restricted maximum likelihood (REML) (Molenberghs et al., 2020).

Reducing the influence of historical excess mortality

When estimating the baseline mortality, the influence on the parameter estimates of historical excess mortality, mainly due to heat waves and seasonal influenza epidemics, needs to be reduced. Two strategies are proposed, both requiring a three‐step analysis, where model (1) is fitted twice.

The first method follows the weighted regression of Farrington et al. (1996) and downweighs historical excess mortality for standardized conditional residuals (Nobre and da Motta Singer, 2007), 1. After fitting model (1), for the first time, a weight based on the standardized residuals is obtained for . Next, a weighted regression model (1) with weights is fitted a second time.

The second method also uses the standardized residuals obtained after fitting model (1) a first time, but down‐weighs the observations by multiplication with the weight: for . Finally, model (1) is fitted a second time, but now on the weighted observations. Observations in historical excess mortality weeks will have higher standardized residuals, which results in a larger reduction of the observation by the weight .

Estimating excess mortality

The weekly prediction and 95% prediction interval (PI) of the baseline mortality during the pandemic (week 11 to 52) are based on the year 2020 year‐specific conditional empirical best linear unbiased predictions of the linear mixed models. The advantage of the year‐specific predictions is that they predict mortality in a specific year, while a marginal prediction, such as the 5‐year weekly average method, predicts mortality in an average year. For the prediction of the weekly mortality by the 5‐year weekly average, we use the prediction interval with unknown mean and variance: with n the number of years ( in our case), and the standard deviation.

The weekly excess mortality results from subtracting the predicted mortality from the observed pandemic weekly mortality. Finally, the weekly estimated excess mortality, and similarly its lower and upper bounds, is summed over all pandemic weeks to result in an estimate and 95% PI of the excess mortality of the year 2020.

Evaluation of models

As discussed in Sections 3.1.1 and 3.1.2, model (1) can be extended with additional correlation structures and in each of these models past excess mortality can be down‐weighted by two methods. To evaluate the appropriateness of the models, we use several statistics. The likelihood ratio test (LRT) compares the −2 log likelihood difference between two nested model with a mixture χ2 distribution. Additionally, the root mean square error percentage (RMSE%) evaluates the forecasting accuracy of the models. If the forecasting error is the difference between the forecasted death and observed death , then Because the proposed models down‐weigh the historical excess mortality, the forecasting accuracy measure can only sensibly be evaluated in years where there has been no or little excess mortality. Historical years with substantial excess mortality will by definition have a large deviation between the observed and forecast deaths for the weeks with excess mortality.

APPLICATION

Linear mixed model (1) is fitted to historical mortality data from week 2 of year 2009 to week 10 of year 2020 for Belgium and the Netherlands and forecast the remaining weeks of the year 2020. In Belgium, the first COVID‐19 related death was reported in week 11 (Bustos Sierra et al., 2020), while four COVID‐19 related deaths occurred in week 10 in the Netherlands. After a relatively mild influenza season during the Winter of 2019–2020, Belgian and Dutch mortality counts during the first 10 weeks of 2020 were lower than average (Figures 1 and 3). This below‐average mortality will influence expected mortality in the following weeks and thus should be taken into account when forecasting.

FIGURE 1

Weekly all‐cause mortality in Belgium from year 2015 to 2020, with the 5‐year average (years 2015–2019) and the sum of the 5‐year average with the reported COVID‐19 deaths, including the 95% prediction interval. This figure appears in color in the electronic version of this article, and any mention of color refers to that version

FIGURE 3

Weekly all‐cause mortality in the Netherlands from year 2015 to 2020, with the 5‐year average (year 2015–2019) and the sum of the 5‐year average with the reported COVID‐19 deaths, including the 95% prediction interval. This figure appears in color in the electronic version of this article, and any mention of color refers to that version

Mortality in Belgium and the Netherlands clearly shows a cyclic pattern (Figures 1 and 3). A correlogram indicates that a yearly cycle is strongly present with a less pronounced half‐yearly cycle. Therefore, both yearly and half‐yearly Fourier series are included into model (1).

The need for modeling additional correlation is evaluated by adding a random effect to each of the Fourier terms in model (1) in turn. In the Supporting Information, it is shown that the model with a random effect on the yearly sine wave fits the data best in most models for both Belgium and the Netherlands. Further modeling of the correlation structure is evaluated by comparing the −2 log likelihood difference of the model including and the expanded models with a distribution for additional serial correlation and mixture distribution for an additional Fourier term random effect.

Belgium

For Belgium, adding either Gaussian serial correlation or an additional random effect for both the weighted regression as the weighted observations strategy, significantly improves the model (Table 1), although the estimated excess mortality from week 11 to 52 for the year 2020 is not very different between all models. In addition to the random effect of the yearly sine wave, only a random effect of the half‐yearly sine wave converged. Other serial correlation structures either do not converge or fit the data significantly worse. For the weighted regression models, the model including Gaussian serial correlation fits the data best, while for the weighted observations, the model with the two random sine wave effects fits the data slightly better. Comparing the log‐likelihood between the weighted regression and weighted observation models is inappropriate given that the observations between the two models are distinct.

TABLE 1

Model fit, forecasting accuracy, and excess mortality estimation (95% PI) comparing linear mixed models, the 5‐year weekly average method and published alternative models for Belgium

	Weighted regression					Weighted observations
Model	−2LL	LRT	RMSE 2014	RMSE 2016	Excess 2020	−2LL	LRT	RMSE 2014	RMSE 2016	Excess 2020
1 + b_1j	6952	/	5.17	5.34	20,586	6682	/	4.87	5.95	20,893
					(18,437;22,738)					(18,843;22,934)
1+ b_1j+ ε_(1)tj	6802	<0.001	4.81	5.14	20,693	6671	<0.001	4.71	5.42	21,008
					(13,137;28,212)					(14,409;27,610)
1+ b_1j+b_3j	6919	<0.001	4.38	4.88	20,467	6666	<0.001	4.25	5.38	20,982
					(18,041;22,900)					(18,639;23,319)
5‐y average	/	/	5.21	5.73	18,989
					(6852;31,122)
BE‐MOMO	/	/	/	/	19,110
					NA
Karlinksy	/	/	/	/	17,421
					(14,799;20,043)
The Economist	/	/	/	/	19,863
					NA
Reported	/	/	/	/	19,288
COVID‐19 deaths

Abbreviations: LL, log likelihood; LRT, Likelihood ratio test; NA, not available; RMSE, root mean square error.

Neither of the years 2014 and 2016 had marked episodes of higher mortality than expected, since influenza related mortality was very low in the year 2014, while 2016 had a heat wave with little excess mortality. As such, the forecasting accuracy of the models for the years 2014 and 2016 is evaluated by excluding deaths from week 11 onward from these years during estimation and forecast the baseline mortality using the first 10 weeks. The forecasting accuracy of the year‐specific linear mixed model in Belgium for both years 2014 and 2016 is better than the 5‐year weekly population average (Table 1). Note that for the 5‐year weekly average, the years 2009–2013 were used to forecast mortality for 2014 and the years 2011–2015 to forecast 2016. Although the differences are small, the model with two random sine wave effects has a slightly better forecasting accuracy. Comparing the weighted regression models with the weighted observation models, there is little difference in forecasting accuracy. For the year 2014, the weighted observation models have slightly better forecasting accuracy, with the reverse holding for 2016. Using the Mean Absolute Error (MAE) or Mean Absolute Percentage Error (MAPE) as forecasting accuracy measure, the conclusions remain the same.

A clear advantage of the linear mixed models over the 5‐year weekly average is the accuracy in estimation of the baseline and excess mortality. The 95% prediction interval of the mortality forecasting is much wider for the 5‐year weekly average (Table 1, Figures 1 and 2). Using the years 2009–2019, rather than only the last 5 years, decreases the variability for the 5‐year average (excess mortality estimate: 19,957 with 95% PI 10,531;29,381), but it is still wider than the linear mixed model variance. When using only the last 5 years to fit the linear mixed models, convergence issues arise because of insufficient data to estimate the correlation structure.

FIGURE 2

Weekly all‐cause mortality in Belgium from year 2015 to 2020, with the mortality forecast by the linear mixed model with two random sine wave effects and the weighted regression methodology and its sum with the reported COVID‐19 deaths, including the 95% prediction interval. This figure appears in color in the electronic version of this article, and any mention of color refers to that version

From week 11 to 52 in 2020, 19,288 COVID‐19 deaths were reported in Belgium (EPISTAT: COVID‐19, 2020). Adding the excess mortality of 1460 deaths from the heat wave during the Summer of 2020 (Bustos Sierra et al., 2020), it seems that the excess mortality forecast of the linear mixed model coincides well with this sum of 20,748 (Figure 2 and Table 1). The excess mortality forecast by BE‐MOMO (Leroy et al., 2021), accessed August 4th 2021 and the simple linear regression models (The Economist, 2020; Karlinsky and Kobak, 2021), accessed August 5th 2021, are clearly lower. Note that the model by Karlinsky excludes the weeks during the heat wave in August 2020 when estimating the excess mortality.

Figures with the mortality forecast of all models in Table 1 are available in the Supporting Information.

The Netherlands

In the Netherlands, the first COVID‐19‐related death was reported on March 6, 2020 (National Institute for Public Health and the Environment, 2021). Since only four deaths were reported in week 10, we chose also for the Netherlands to initiate forecasting mortality from week 11 onward. In the Netherlands, the year 2020 started with a lower‐than‐average mortality during the first 10 weeks (Figure 3). The likelihood ratio tests show that adding a correlation structure to model (1) significantly improves the model fit for the historical mortality data in the Netherlands (Table 2). Although the variability between the estimation of the excess mortality of the different models is larger than in the Belgian analysis, the precision of each estimation is much better than the 5‐year weekly average (Table 2, Figures 3, 4). Similar to Belgium, only data from 2009 to 2019 have been used to fit the linear mixed models. Using data from 1995 onward for the Netherlands does not reduce the variance of the excess estimation by the linear mixed models in any relevant way.

TABLE 2

Model fit and excess mortality estimation (95% PI) comparing linear mixed models, the 5‐year weekly average method and published alternative models for the Netherlands

	Weighted regression			Weighted observations
Model	−2LL	LRT	Excess 2020	−2LL	LRT	Excess 2020
1+b_1j	7109	/	20,025	7020	/	21,125
			(15,634;24,414)			(16,663;25,589)
1+b_1j+ε_(1)tj	6935	<0.001	20,698	7001	<0.001	20,727
			(11,534;29,858)			(11,337;30,113)
1+b_1j+b_3j	7044	<0.001	20,585	7004	<0.001	22,796
			(15,023;26,145)			(17,308;28,277)
5‐y average	/	/	19,024
			(5324;32,726)
EURO‐MOMO	/	/	15,807
			NA
Karlinksy	/	/	15,739
			(13,003;18,475)
The Economist	/	/	16,700
			NA
Reported	/	/	11,527
COVID‐19 deaths

Abbreviations: LL, log likelihood; LRT, likelihood ratio test.

FIGURE 4

Weekly all‐cause mortality in the Netherlands from year 2015 to 2020, with the mortality forecast by the linear mixed model with two random sine wave effects and the weighted regression methodology and its sum with the reported COVID‐19 deaths, including the 95% prediction interval. This figure appears in color in the electronic version of this article, and any mention of color refers to that version

Contrary to Belgium, in the Netherlands the reported COVID‐19 mortality (11,527 across weeks 11–52 in 2020) is well below the estimated excess mortality in all models (Figures 3, 4 and Table 2). Taking account of the heat wave in the Netherlands in week 33, which resulted in an estimated 400 excess deaths (National Institute for Public Health and the Environment, 2021), the linear mixed models estimate that between 51% and 56% of the COVID‐19 mortality has been reported in 2020. The excess mortality forecast by EURO‐MOMO (Central Bureau for Statistics, 2021c), accessed August 4, 2021 and the simple linear regression models (The Economist, 2020; Karlinsky and Kobak, 2021), accessed August 5, 2021, are lower than both the linear mixed models and the 5‐year average. Again, the model by Karlinsky excludes the weeks during the heat wave in August 2020 when estimating the excess mortality.

Figures with the mortality forecast of all models in Table 2 are available in the Supporting Information.

DISCUSSION

In 2020, the world was confronted with the most lethal pandemic in 100 years, the severity of which is underscored by the all‐cause mortality from Belgium and the Netherlands. They show how hard we were hit, but also how important the mitigation measures against the spread of the virus are. To understand the complete picture of the pandemic, for individual countries and for the comparison between countries, it is useful to estimate the excess mortality, since both direct and indirect effects of the pandemic on mortality are captured by excess mortality.

Determining excess mortality requires the estimation of mortality under nonpandemic conditions. The often used method of averaging the 5‐year historical weekly mortality, however, ignores the trend in mortality of the first weeks of 2020 by estimating the population average and may be influenced by recent excess mortality in the past. The advantage of the weekly average method is that it is easy to apply and only information about mortality is required. We propose linear mixed models to forecast year‐specific baseline mortality, which address the shortcomings of the 5‐year weekly averaging method, while maintaining simplicity in application and limited data requirements.

We modeled the correlation between mortality in consecutive weeks and between weeks in a year via several correlation structures in the linear mixed models, but could not identify a single best model. This underscores that the inclusion of a serial correlation structure is more important than its precise parametric form (Verbeke et al., 1998). Also between the two methods for down‐weighting historical excess mortality there is no clear winner. Therefore, we choose not to recommend any of the variations to our model, but rather recommend to down‐weigh past excess mortality and include a serial correlation structure.

For Belgium and the Netherlands, we have shown that the linear mixed models not only fit the mortality data better, but also that the prediction for the years 2014 and 2016 are superior to the 5‐year weekly average. For Belgium, the excess mortality in 2020 is estimated by the linear mixed models to lie between 20,467 and 21,008. Taking account of 1460 excess deaths during the Summer heat wave in 2020 (Bustos Sierra et al., 2020), between 19,007 and 19,548 deaths can then be attributed to direct and indirect effects of the COVID‐19 pandemic. Since 19,288 COVID‐19 related deaths have been reported in Belgium by the Belgian Institute for Public Health (EPISTAT: COVID‐19, 2020), the reported COVID‐19 mortality coincides fairly well with the excess mortality. A priori, this coincidence is not a solid proof for COVID‐19 deaths being well reported in Belgium. That said, it is independently known that the country reports not only deaths from COVID‐19 laboratory test or chest CT scan confirmed cases, but also deaths in clinically confirmed COVID‐19 cases (Renard et al., 2020). This was done regardless of the setting (with hospitals and nursing homes the most important ones). Also, COVID‐19 has been a reportable disease since the pandemic's outbreak. It is therefore not surprising that around 90% of deaths suspected to be COVID‐19 related are caused by COVID‐19 (Grewelle and De Leo, 2020). For the Netherlands, the excess mortality in 2020 is estimated by the linear mixed models to lie between 20,585 and 22,796, although only 11,527 COVID‐19 deaths have been reported by the National Institute for Public Health and the Environment (2021). Likely, the Netherlands have under‐reported COVID‐19 deaths to an estimated 51–56%. Recently, the Dutch Central Bureau for Statistics (2021b) has attributed a little more than 20,000 deaths in 2020 to COVID‐19 (Central Bureau for Statistics, 2021a), which is in line with the estimation by the linear mixed models.

The excess mortality in 2020 for Belgium and the Netherlands, estimated by the EURO‐MOMO models or variations thereof, as well as by simple linear models, is estimated lower in both countries than the estimates by the linear mixed models. The EURO‐MOMO models forecast mortality by ignoring past historical excess mortality, which has been criticized (Aron et al., 2020). BE‐MOMO, a variation of these models, uses all historical mortality data, but requires additional data on climate and influenza, which may be difficult to obtain in general. Simple linear models do not down‐weigh past excess mortality and the independent error assumption is likely violated in longitudinal mortality data. The linear mixed models achieve a compromise between flexibility and simplicity. It is a simple approach that captures seasonal and yearly variation, while reducing the influence of past excess mortality.

Both 2020 waves of SARS‐CoV‐2 infections in both Belgium and the Netherlands have led to thousands of COVID‐related deaths. Although the first wave was shorter and more intense, the second wave, which started roughly at week 30 and continued after the study period investigated here, was longer, resulting in more or less equal numbers of COVID‐19 deaths in both periods. As the timing of the deaths are parallel to the excess mortality, it shows that likely not the nonpharmaceutical interventions, but the virus itself is responsible for the majority of the excess mortality.

It is of course artificial to evaluate the effect of the pandemic on mortality by calendar year. The full impact of COVID‐19 on mortality will be evident when sufficient individuals in the population will be vaccinated and endemicity reached. This would also allow for international comparison between countries with a different timing of the pandemic or with different mitigation measures. In the meantime, using the linear mixed models allows for intermediate evaluation.

Finally, excess mortality over an entire country does not capture the entire effect of the COVID‐19 pandemic on mortality. Age groups, gender, and regional differences within a country are smoothed out in an overall excess‐mortality evaluation. To account for these differences, the linear mixed models can be extended to include these variables or can be applied to the different subgroups. They can also be applied to historical data that are interrupted, for example, if only seasonal historical mortality data are available. When applying the linear mixed models to sub‐populations, it is important to evaluate whether the death count is still sufficiently high for the central limit theorem to be invoked.

Web Appendices, Tables, and Figures referenced in Section 4 are available with this paper at the Biometrics website on Wiley Online Library. These include all technical details, and additional simulation results and data analyses. SAS code to reproduce all results are available online at https://ibiostat.be/online‐resources/longitudinal#SASLLM_Mort.

Click here for additional data file.