The study focused on the recorded deaths and the causes of death of all citizens currently living in the territory of Belgrade for the period from 2015 to 2021. The age-specific and sex-stratified midyear population estimates for the territory of Belgrade were collected from the electronic databases of the Statistical Office of the Republic of Serbia for the given period [18]. These midyear population estimates were calculated considering the natural population change and net migration rates. Overall mortality data, as well as data on the mortality from COVID-19 as the main cause of death, were collected. The mortality data for COVID-19 as the primary cause of death were coded per the WHO guidelines [19]. These guidelines define COVID-19 as a cause of death if clinical death occurs, which is compatible with COVID-19 disease, in a case with confirmed or probable COVID-19 disease, and no alternative, more probable causes exist. For this analysis, we used the Codes for special purposes (U00-U89), Chapter XXII, ICD-10, where only deaths with the main cause of death–COVID-19 (U07.1, U07.2) were exclusively attributed for 2020 and 2021. This coding scheme was used to derive crude COVID-19 and non-COVID-19 mortality estimates. The data were also sourced from the electronic databases of the Statistical Office of the Republic of Serbia [18].

We first calculated the standardized direct death rates for the populations from 2015 to 2019 to estimate the excess deaths [20]. The 2020 and 2021 population estimates used this standardized death rate to generate age-specific and sex-stratified mortality. This total was then compared to the actual crude death rate. Direct standardization is not without its limits. It can severely underestimate excess deaths due to COVID-19. The 5-year average estimate does not capture seasonal variation and year-to-year mortality. To account for these year-to-year trends, we use a mixed-effects hierarchical Poisson generalized linear regression model (GLM) to predict mortality for 2020 and 2021 using a yearly time series of sex, age group membership, and the logarithm of absolute population count in the region. The model was fit to mortality data observed from 2015 to 2019. To fit the Bayesian model, we estimated the parameter posterior distributions using Hamiltonian Markov Chain Monte Carlo, as implemented in Stan (version 2.26.1). The posterior distributions quantify which parameter values are most consistent with the data and are subsequently used for generating predictions. We assessed convergence by the inspection of trace plots and R-hat values, which should be below 1.01 [21], and the effective sample size, which should be greater than 1000 [22]. These quality assurance measures indicate that the sampler has achieved robust convergence and that the tails of the parameter distributions are reliable when estimating confidence intervals of excess mortality. Model out-of-sample predictive performance was assessed using leave-one-out cross-validation approximated via Pareto smoothed importance sampling [23]. Leave-one-out cross-validation metrics estimate model performance on unseen data and are useful measures of out-of-sample predictive performance. The leave-one-out cross-validation indicated that the predictors contain valuable predictive information. Priors were chosen through prior predictive simulation so that predata predictions span the range of scientifically plausible outcomes. Using weakly informative priors ensures that the model is partially regularized and protects against implausible parameter values. We constructed posterior predictive distributions to assess the model fit adequacy to evaluate predictions at the observed data’s mean (Multimedia Appendix 1) [24]. Posterior predictive distributions evaluated the models’ predictions for consistency against the observed data before any inference could be made. Furthermore, to assess model calibration across the whole observed distribution, we compared the probability integral transformed values of the predictive distribution to our data’s empirical (observed) cumulative distribution function (Multimedia Appendix 2) [25]. The model convergence diagnostic is shown in Multimedia Appendix 3. After successful model fitting, we generated predictions for each age group for the next two pandemic years and compared the expected and the observed mortalities.

The raw number of excess deaths was not comparable across countries due to significant differences in the population. To enable easier comparison across countries, we calculated the P-score. The P-Score measures the percentage difference between the reported and projected deaths.

All calculations were performed using R software (version 4.3.1; R Foundation for Statistical Computing).

This study was the secondary analysis of publicly available data collected by the Statistical Office of the Republic of Serbia. Data collection adhered to strict ethical guidelines ensuring the confidentiality and privacy of all participants. All collected data were fully anonymized to prevent the identification of individuals.

From 2015 to 2019, the territory of Belgrade had an average population of 1,688,185, of which 797,080 (47.22%) citizens were male and 891,105 (52.78%) were female. The largest number of people were included in the 40‐59 year age group (459,345), closely followed by the 20‐39 year age group (459,032). The number of people aged >80 years was the lowest (72,757). These demographics remained relatively stable throughout 2021 (Multimedia Appendix 4). The raw population characteristics are provided in Multimedia Appendix 5.

In the territory of Belgrade, there were 27,722 and 32,593 deaths in 2020 and 2021, respectively. The average number of deaths in the previous 5-year period was 23,792. In 2020 and 2021, there were 3028 (12%) and 7669 (24%) deaths when COVID-19 was coded as the primary cause.

The crude death rate increased significantly in 2020 and 2021 from its average of 14.1 per 1000 people in the previous 5-year period to an average of 17.85 per 1000 people. Those in their seventh and eighth decade were the most severely impacted age groups. Age-specific crude death rates for citizens in their seventh decade of life were 33.90 per 1000 people and 40.08 per 1000 people in 2020 and 2021, respectively. Those in their eight-decade had crude death rates of 153.24 per 1000 people and 183.91 per 1000 people for 2020 and 2021, respectively. Raw death rate trends are shown in Figure 1.

The total number of excess deaths, calculated from model estimates, was 3175 deaths (99% CI 1715-4094) for 2020 and 8321 deaths (99% CI 6975-9197) for 2021. The ratio of estimated excess deaths to reported COVID-19 deaths was 1.07. Excess mortality, stratified by age group and sex, is shown in Figures 2-4. The estimated increase in mortality during 2020 was 12.93% (99% CI 15.74-17.33%). During 2021, mortality increased by 39.32% (99% CI 35.91-39.32%) from the expected value. The average increase in mortality during 2020 and 2021 was 23.55%, compared to the expected mortality extrapolated from the previous 5-year period. There seemed to be fewer deaths than expected based on the predicted trends in children and young people, which aligned with the raw observed data. To be precise, between 40 and 106 and between 44 and 110 fewer deaths occurred among children and young people in 2020 and 2021, respectively; this amounted to an average decrease in mortality of 22.43% and 23.71% for this age group. The most impacted groups were people in their seventh and eighth decades. People aged 60‐79 years experienced an estimated 1279 to 2143 and 3594 to 4423 excess deaths in 2020 and 2021, respectively. People aged >80 years experienced an estimated 504 to 1625 and 2895 to 3858 excess deaths in 2020 and 2021, respectively. In people aged 60 to 79 years, there was an estimated increase in mortality of 16.95% and 38.50% in 2020 and 2021, respectively. For those aged >80 years, the increase was estimated at 11.50% and 34.14% in 2020 and 2021, respectively. For those aged up to 39 years old, there was little increase in mortality (4.72%) during 2020. However, in 2021, even those aged 20‐39 had an estimated increase in mortality of 32.95%. The crude, standardized, and model-predicted death rates are directly compared in Figure 3 and Tables 1 and 2. The model-predicted deaths matched the number of non-COVID deaths recorded in the territory of Belgrade (Figure 3). This concordance between the predicted and recorded non-COVID deaths provided evidence that the cause-of-death misattribution did not occur in the territory of Belgrade. The error whiskers represent 99% credible intervals.

When revised death statistics were released in Serbia, some authors noted significant discrepancies between the preliminary and finalized COVID-19 dataset [1,2]. Specifically, considerable mortality underreporting occurred in the initial dataset compared to the finalized COVID-19 mortality dataset. We collected data on the number of deaths in Belgrade, coding the SARS-CoV-2 as the primary cause, to validate the accuracy of this finalized dataset. We opted for a modeling-based approach because classical standardization does not address the hypothesized problem of death misattribution, which involves coding another primary cause of death as a COVID-19 death [20]. We fitted a Poisson GLM model to the acquired data for 2015‐2019, and then used the model to generate mortality predictions for 2020 and 2021. Given the expected age and comorbidity structure within the population over the studied time frame, the model’s mortality estimates reflect the expected deaths in the absence of the pandemic. When we removed COVID-19 as a cause of death from the dataset, the model’s predictions matched the reported deaths, indicating that there were no significant cause-of-death misattributions in the final published dataset for Belgrade. Furthermore, this establishes a firm link between excess deaths and COVID-19 in Belgrade.

Possible mechanisms for this undercount discrepancy have yet to be fully investigated. Differences in the case definition were a common explanation for reporting heterogeneity between countries [5,26,27]. However, according to official reports, Serbia adopted the case definition proposed by the WHO early on during the pandemic, which remained unchanged throughout [1]. The ratio of excess mortality to the officially reported COVID-19 death count (undercount ratio) for the amended dataset was 1.07, indicating slight undercounting. This is a substantial improvement compared to the ratio of greater than 3.00 reported in the preliminary dataset [1]. The undercount ratio was larger in 2020 than in 2021, with a slight decrease in divergence. This improvement is consistent with a worldwide improvement in undercount ratios throughout the pandemic, attributed to improvements in accuracy, case definition revisions, and better reporting [3]. However, substantial undercounting remained in the preliminary dataset. Our analysis indicates that the revised datasets accurately reflect actual COVID-19 mortality, at least in Belgrade. This suggests that accurate death certificate data were collected, and that the problem was localized to the reporting system. No prior reporting system of this scale had been used in Serbia; thus, a learning curve could have existed initially, causing delays. However, the divergence between the actual and reported deaths continues into 2021, making any attribution to a learning curve unlikely.

Despite intense government measures to curb the spread of the SARS-CoV-2, the territory of Belgrade had the most considerable death rates in the country [2]. There is a complex interplay of contributing causes to this increase in mortality. Belgrade is home to one of the world’s largest medical health care centres. The University Clinical Centre of Serbia possesses a hospital floor capacity of 3150 beds [28]. The health care system in Serbia is heavily centralized, with most medical services provided by the state. Overreliance on the centralized health care system may have contributed to delays in reporting death, particularly given worldwide reports of barriers to access in rural areas [29-31]. Health care disparities between communities of differing socioeconomic status and mobility are a well-recognized moderator of pandemic impact [31-34]. The government’s response included the incorporation of many hospitals into a ‘COVID-19 Response System,’ which had limited capacity to treat non-COVID-19 diseases. Thus, treatment priorities could have contributed to limited health care access and might have increased morbidity and mortality due to other causes. However, as our model’s death predictions account for the non-COVID deaths almost perfectly, especially in the lower age subgroups, deaths from other causes could not have been acute—or at least, could not have occurred during 2020 and 2021. Increased mortality and morbidity due to the delayed diagnosis of oncological and other chronic diseases in the following years can be an expected consequence, which could manifest as a delayed increases in disease burden. Further studies in this area would shed light on the delayed impact of COVID-19.

On assessing the stringency index, a composite measure of 9 government response metrics, we found that Serbia quickly reached the maximal (strictest response) by March 22, 2020. Both stringency as well as containment and health indexes (an extension of the stringency index) settled into a 40‐70 range during 2020‐2021, briefly dropping below 40 for the period leading up to and including the June 21, 2020 elections [35]. Increased deaths led to more stringent government responses, while every reduction led to waning compliance [36]. Over time, people become less risk-averse, and government policies lose effectiveness, especially after initiating the vaccination campaign. The relative percentage of people who received all doses prescribed by the initial vaccination protocol was 46.87% by January 1, 2022, and showed significant plateauing [37]. The relative effect of nonpharmaceutical interventions in general declined by 12% from May 2021 due to less stringent methods and the introduction of vaccination strategies [38].

Mortality from COVID-19 is more significant among those with advanced age, mimicking overall mortality [39-41]. The territory of Belgrade observed a reduction in the mortality rate of children and young people (aged 0‐19 years) during 2020 and 2021. COVID-19 deaths in this age group were very low (1 case in 2020 for males and 0 for females). Most of the mortality risk in children and young people comes from injuries and infectious diseases [42]. Government intervention and the mobility restrictions during the pandemic probably contributed to lower mortality in children and young people by indirectly lowering their overall risk. For example, traffic accident fatalities decreased in the Western Balkans during lockdowns [43]. The evidence suggests a decrease in pediatric intensive care unit admissions between 2020 and 2021. Additionally, research has found no increase in adolescent self-harm or suicide [44]. These effects contributed to a decrease in the death rate among the child and adolescent population. However, the changes in health care functioning during the same period could have had long-term consequences in this age group, which have yet to manifest.

We acknowledge that the estimates obtained are Bayesian and are posterior distributions of expected deaths. As some of the raw non-COVID-19 deaths are on the edges of our estimates, where the predicted response parameter is least likely (given our data), there is room to argue that deaths from other causes might have increased for specific age groups. In contrast, other age groups showed fewer deaths than expected. Additionally, the resolution of our data limited our modeling choice, as the most accurate dataset only includes measurements from yearly reports.

The finalized mortality dataset for Belgrade can be safely used in COVID-19 impact analysis. In 2020 and 2021, Belgrade experienced a significant increase in mortality, with the majority of the excess mortality attributable to the SARS-CoV-2. Concerns about increased mortality from causes other than COVID-19 seem misplaced, as their impact appears negligible. However, we do caution that the delayed effects of the SARS-CoV-2, barriers to health care, and the psychological impact of the COVID-19 pandemic on Belgrade territory remain unexplored.