### Quantification of excess mortality

The study focuses on a retrospective analysis and comparison of the (standardized) excess mortality of Belarus with neighboring countries. Basically, the *excess mortality* over a period of time (e.g. a month) is defined as the actual number of deaths minus the number of deaths that would have been “normally expected”^{26}, where “normally expected” can be interpreted quite broadly. In this work, the formal definition of “excess deaths” over a period of time (i.e. one month) is the difference between the actual number of deaths minus the “expected number” that are predicted from modeling statistical. Such modeling can incorporate auxiliary data such as age characteristics of the population over time. The approach is to formulate and fit models to historical mortality data on *not* periods (i.e. months) that directly precede the start date of the outbreak. Then, the same models are used to produce the “expected” predicted mortality during the epidemic period for comparison with the reported mortality data. The differences between reported and predicted produce the formal quantification of “excess mortality”.

### Standardized Mortality Scores

Although modeling efforts provide a formal quantification of “excess” mortality numbers, they lack a meaningful way to compare these numbers between countries or regions. Subsequent comparison may be highly desirable to compare the effectiveness of applied NPIs.

Formal comparisons between countries (or regions) can be made via the *standardized scores*^{40}. In particular, the *P*-scores allow comparison of mortality rates for different time periods and populations of different regions with different sizes^{41}. Scores are defined non-parametrically or parametrically^{26} and both are evaluated in this work. The nonparametric version of *I* periods studied has the form:

$$begin{aligned} P(t_i) = frac{x(t_i) – bar{x}(t_i)}{bar{x}(t_i)}, text{ for } , i= 1,2, dots , I, end{aligned}$$

(1)

or (bar{x}(t_i)) is the average mortality of all (t_i) periods of the past *not* years preceding the epidemic period studied.

The parametric version of *I* study periods has the form:

$$begin{aligned} mathscr {P}(t_i) = frac{x(t_i) – hat{mu }(t_i)}{hat{mu }(t_i)}, text{ for } , i=1,2, dots , I end{aligned}$$

(2)

or (hat{mu }(t_i)) is the population mortality predicted by the model for (t_i) period studied based on fitting the model using *not* years of data preceding the epidemic period studied. The differences between nonparametric (1) and parametric (2) and supposedly minor^{26}. Sheet music (P(t_i)) and (mathscr {P}(t_i)) are usually converted to a percentage for easier interpretation. For example, the converted score value of 120.47 indicates that mortality is exceeded by 120.47%. As the mortality scores are on the same scale regardless of the size of the population studied, they can be easily interpreted and compared. *between regions*.

Non-parametric (1) and parametric (2) scores are not limited to mortality data and corresponding time series. The same scores can be calculated and investigated for any time series such as Google Trends search queries, which is also done in this work.

### Predictive models for parametric estimates

The model chosen to produce the expected estimates (hat{mu }(t_i)) for (i=1,2, dots , I) because formula (2) must take into account the (competing) processes that occur simultaneously within the population under study and must incorporate the ability to disentangle the effects of each process on mortality. In particular, the model should at least incorporate the annual seasonality of mortality^{42}potential time lags as well as changes in the structure and size of the population being studied^{43}. For COVID-19, demographic change is crucial, as the Belarusian population is aging and the number of high-risk people over 65 is increasing^{44}. These demographic characteristics have not been incorporated into previous excess mortality studies^{25}, while such changes can potentially affect mortality dynamics. Two time series models were used: the Prophet forecast model^{45} and an autoregressive integrated moving average (ARIMA) model^{46}.

The model of the Prophet^{45} is implemented in R wrap Prophet and consists of three main elements: trend, seasonality and holidays. The model is formulated as follows:

$$begin{aligned} y

(3)

for a given time interval *you*. Formula (3) consists of a dependent (response) model variable *there*(*you*) that represents the number of deaths or search trend values, independent (exploratory) variables *g*(*you*), *s*(*you*) and *h*(*you*) which represent the trends of the model, and the error of the model (epsilon _t). Variable *g*(*you*) represents the non-periodic modeling trend component, *s*(*you*)—periodic modeling component, and *h*(*you*)—extreme events (or holidays) that may occur with an irregular schedule. The model error (epsilon _t) represents processes that are not explained by the model. All model variables and error depend on time *you*. Model (3) can also incorporate a covariate, i.e. an independent predictor for *there*(*you*)^{47}. Therefore, the number of individuals in the population of Belarus aged 65 or older was considered as such a covariate in the model fitting process.

Autoregressive Integrated Moving Average (ARIMA) model is implemented in R wrap provide^{46}. The ARIMA model is determined by the set of three parameters (*p*, *D*, *q*). The parameter *p* determines the model offset i.e. the number of previous observations (Y(t-1), Y(t-2), dots , Y(tp)) which are used to adjust the current observation *Yes*(*you*). The parameter *D* represents the “degree” of the model, i.e. how many times the previous response must be subtracted to guarantee the required stationarity assumption of the model^{48}. For example, for (j=0) the newly defined response variable is (y

(4)

Formula (4) contains the response pattern variable *there*(*you*) which represents the number of deaths (or search trend values) and is parameterized by the global mean parameter (mu)autoregressive parameters (varphi _i) or (i=1,2,points , p)moving average settings (theta _j) and model errors (epsilon _{tj}) or (j=0,1,2,points , q).

### Score calculations

The non-parametric *P*– January to June 2020 mortality data scores were calculated based on monthly averages (bar{x}(t_i)) for the first six months of the year, i.e. for (t_1, t_2, dots , t_6). The averages were based on *five* (2015-2019) and *nine* (2011-2019) years of data. The analogous calculations of *P*-scores were preformed for three Google trends based on *five* (2015-2019) years of data.

The parametric (mathscr {P})-scores for mortality data were calculated based on predicted means (hat{mu }(t_i)) which came from the Prophet forecasting model and the ARIMA model. Both models were fitted to two time series with *five* (January 2015-February 2020) and *nine* (January 2011-February 2020) years of monthly mortality. For all adapts to the *four* the following epidemic months (March 2020-June 2020) were predicted by these models. The Prophet and ARIMA models were equipped *both* without covariate and with a covariate which was the number of individuals aged 65 or over. For the Prophet model, the default package settings were used while for the ARIMA model auto.arima function was used for *nine* years to determine (*p*, *D*, *q*) settings. The corresponding ARIMA parameters were used to *five* years adapts since those auto.arima options were found to be the most reliable.

For the model Prophet the (95%) the upper bound of the model fits (and the forecasts for the last four months) were used to produce the expected mortality (hat{mu }(t_i)) and the correspondent (mathscr {P}(t_i))-notes (2). These estimates are conservative since they use the upper bound of expected expected mortality, thus producing the lower bound estimate of (mathscr {P}(t_i))-scores. the provide The package only provides ARIMA upper uncertainty bounds for model predictions and not for model fits. Therefore, instead of an upper limit of 95%, these values produced by the model (hat{mu }(t_i)) *both* for the adjustment and for the planned four months were used to produce the (mathscr {P}(t_i))-scores.

Analog scores (2) based on *five* (January 2015-February 2020) years of monthly trend data and the Prophet model were used to produce the expected averages (hat{mu }(t_i)) and the correspondent (mathscr {P}(t_i))-scores from three Google Trends. Only the Prophet model was fitted for Google Trends since the results of the Prophet and ARIMA models for mortality were very similar.

### Coupling Mortality Data and Google Trends

To investigate associations and (potential) relationships between mortality and Google Trends, the corresponding time series were compared in pairs. The comparisons were made for the period when they were all available, i.e. for the monthly readings from January 2015 to June 2020. To facilitate visual comparison, the series have been normalized to be on the same scale. Specifically, the values of each time series were divided by the corresponding median of the values in that series and the resulting values were multiplied by 100. In addition to this, the four standardized series were smoothed using loess function in R with span = 0.25 parameter to remove potential noise. The resulting standardized and standardized values were compared using Granger’s dependency test for time series *with dependent* recordings^{49}. Graphic summaries were also produced.