Abstract

Introduction

The world is facing a pandemic due to a novel coronavirus (Sars-Cov-2) since January 30, 2020, when the World Health Organization (WHO) declared the outbreak a public health emergency of international concern (¹1. WHO Coronavirus Disease (COVID-19) [Internet]. Covid19.who.int. 2020 [cited 11 August 2020]. Available from: https://covid19.who.int/.
https://covid19.who.int/... ). The first reported case of the new coronavirus disease, named COVID-19, was in early December in Wuhan, China (²2. Huang C, Wang Y, Li X, Ren L, Zhao J, Hu Y, et al. Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China. Lancet 2020; 395: 497–506, doi: 10.1016/S0140-6736(20)30183-5.
https://doi.org/10.1016/S0140-6736(20)30... ), and as of August 11, 2020, the disease has killed more than 732,000 people and infected almost 20 million worldwide (¹1. WHO Coronavirus Disease (COVID-19) [Internet]. Covid19.who.int. 2020 [cited 11 August 2020]. Available from: https://covid19.who.int/.
https://covid19.who.int/... ). As of December 3, 2020, Brazil was the second country with most deaths and infections due to COVID-19 with more than 3 million reported cases and 101,049 deaths, only behind the United States of America that has almost 5 million reported cases and 161,547 deaths (¹1. WHO Coronavirus Disease (COVID-19) [Internet]. Covid19.who.int. 2020 [cited 11 August 2020]. Available from: https://covid19.who.int/.
https://covid19.who.int/... ). The first death by COVID-19 in Brazil occurred in the city of São Paulo on March 16, 2020.

Most initial cases of COVID-19 in western countries were of individuals returning from China (³3. Holshue M, DeBolt C, Lindquist S, Lofy K, Wiesman J, Bruce H, et al. First case of 2019 novel coronavirus in the United States. N Eng J Med 2020; 382: 929–936, doi: 10.1056/NEJMoa2001191.
https://doi.org/10.1056/NEJMoa2001191... ) and Italy (⁴4. Escalera-Antezana J, Lizon-Ferrufino N, Maldonado-Alanoca A, Alarcón-De-la-Vega G, Alvarado-Arnez L, Balderrama-Saavedra M, et al. Clinical features of the first cases and a cluster of Coronavirus Disease 2019 (COVID-19) in Bolivia imported from Italy and Spain. Travel Med Infect Dis 2020; 35: 101653, doi: 10.1016/j.tmaid.2020.101653.
https://doi.org/10.1016/j.tmaid.2020.101... ), but many countries still have an issue for establishing when the first case of local transmission occurred. This is a consequence of many problems; one of them is the existence of a fraction of asymptomatic individuals (due to mild or lack of symptoms) that facilitated the dissemination of the coronavirus before its detection by the health system (⁵5. Li R, Pei S, Chen B, Song Y, Zhang T, Yang W, et al. Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2). Science 2020; 368: 489–493, doi: 10.1126/science.abb3221.
https://doi.org/10.1126/science.abb3221... ), combined with insufficient testing, which limited the capacity to count the number of infections and, consequently, establish the real timeline of the coronavirus infection. On the other hand, the excess mortality rate can provide valuable information about the dynamics of the pandemic (⁶6. Subbaraman N. Why daily death tolls have become unusually important in understanding the coronavirus pandemic. Nature 2020.). Recently, this approach was used in several studies. Li et al. established that a great number of infections in China happened long before social restrictions and control policies took place, revealing a large number of undocumented cases (⁵5. Li R, Pei S, Chen B, Song Y, Zhang T, Yang W, et al. Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2). Science 2020; 368: 489–493, doi: 10.1126/science.abb3221.
https://doi.org/10.1126/science.abb3221... ). Also, Delatorre et al. (⁷7. Delatorre E, Mir D, Gräf T, Bello G. Tracking the onset date of the community spread of SARS- CoV-2 in western countries. Mem Inst Oswaldo Cruz 2020; 115: e200183, doi: 10.1590/0074-02760200183.
https://doi.org/10.1590/0074-02760200183... ) confirmed such findings for multiple countries, including the USA, UK, and Italy. For Brazil, they found that COVID-19 was already circulating in the Brazilian population in early February (between January 31st and February 7th), 2020.

Following the trend of using the excess number of deaths, many time-series models can be used to forecast new deaths and understand the coronavirus dynamics. One of them is the auto regressive integrated moving average (ARIMA) model. The ARIMA model, or Box-Jenkins model, was developed in 1976 by George Box and Gwilym Jenkins, who introduced the concept of transforming non-stationary data into stationary (⁸8. Box G, Jenkins G. Time series analysis forecasting and control. San Francisco: Holden-Day; 1976.). Data is considered stationary when its parameters, mean and variance, do not depend on time (i.e., are constant throughout time) (⁹9. Sato RC. Disease management with ARIMA model in time series. Einstein (São Paulo) 2013; 11: 128–131, doi: 10.1590/S1679-45082013000100024.
https://doi.org/10.1590/S1679-4508201300... ). The transformation introduced by Box and Jenkins is called seasonal differencing and consists of computing the differences between consecutive observations.

Using the excess number of deaths is preferred over the number of notified infected cases, mainly because of the small number of tested people leading to sub-notification in the number of cases (⁷7. Delatorre E, Mir D, Gräf T, Bello G. Tracking the onset date of the community spread of SARS- CoV-2 in western countries. Mem Inst Oswaldo Cruz 2020; 115: e200183, doi: 10.1590/0074-02760200183.
https://doi.org/10.1590/0074-02760200183... ). Therefore, the objective of this study was to explore the number of autopsies performed at the São Paulo Autopsy Service (SPAS) to find a possible excess in the number of deaths before the first official death by Sars-CoV-2 in the city of São Paulo. This analysis can help us to establish if deaths related to the coronavirus occurred previous to March 16, 2020, implying a possible early community transmission of the virus and the spread of the disease in Brazil for some time before the first official case.

Material and Methods

Participants

São Paulo is the largest city in Brazil, with an estimated population near 12 million in 2019 (¹⁰10. Fundação Sistema Estadual De Análise De Dados (SEADE). 2020. Available from: http://produtos.seade.gov.br/produtos/projpop/index.php.
http://produtos.seade.gov.br/produtos/pr... ). Data on mortality were collected from the SPAS, which since 1931 is responsible for clarifying the cause of non-traumatic deaths that were due to ill-defined disease or of people without medical assistance in the municipality of São Paulo. SPAS performs more than 14,000 full-body autopsies yearly, accounting for almost 20% of the autopsies in the city of São Paulo (¹¹11. SVO - Serviço De Verificação De =bitos Da Capital. Introdução: Principais atividades. 2020. Available from: http://www.svoc.usp.br/Introducao.htm.
http://www.svoc.usp.br/Introducao.htm... ).

Socioeconomic and clinical data are collected by the SPAS staff from the next of kin. The causes of death determined with the information from the autopsy examination were coded according to the 10th revision of the International Statistical Classification of Diseases and Related Health Problems (ICD-10) (¹²12. ICD-10. Geneva: World Health Organization; 2011.).

The SPAS stopped the conventional autopsy procedure on March 21st, 2020 following São Paulo State regulations based on the high risk of staff contamination during the COVID-19 pandemic, which caused an interruption in data collection. Therefore, this study included death data from January 1st, 2019 to March 20th, 2020. We excluded individuals younger than 18 years (n=1,641), cases for which the date of autopsy was after March 20, 2020, and a few cases with data processing errors (n=63). Family members signed an informed consent after a SPAS employee explained the autopsy procedure and the possible participation in a research project.

Autopsy evaluation

Individuals who die from an undefined cause are referred to the SPAS, where a full-body autopsy is performed to determine the cause of death. Upon the arrival of the deceased body at the SPAS, the next of kin is requested to answer a few questions regarding the moment of death, previous diseases, and lifestyle factors. The average postmortem interval between death and autopsy exam in our sample was 24.77±5.64 h.

Autopsies are performed by certified pathologists that examine major organs macroscopically. When necessary, tissue samples are collected for additional microscopic analysis. Afterward, the pathologist determines the cause of death based on the findings from the autopsy exam. The basic and immediate causes of death are determined, as well as additional pathologies that could have contributed to the death. For the all-cause mortality analysis, all ICD-10s (¹²12. ICD-10. Geneva: World Health Organization; 2011.) are included. Also, specific-cause mortality is determined by selecting ICD-10s that were possibly related to the coronavirus: pneumonia (J18), respiratory-related (ICD-10s in Table 1), and coronary artery disease (I21 and I25) cases (¹³13. Vetter P, Vu D, L'Huillier A, Schibler M, Kaiser L, Jacquerioz F. Clinical features of covid-19. BMJ 2020; 369: 1470, doi: 10.1136/bmj.m1470.
https://doi.org/10.1136/bmj.m1470... ,¹⁴14. Zhang Y, Xiao M, Zhang S, Xia P, Cao W, Jiang W. et al. Coagulopathy and antiphospholipid antibodies in patients with Covid-19. N Engl J Med 2020; 382: e38, doi: 10.1056/NEJMc2007575.
https://doi.org/10.1056/NEJMc2007575... ).

Model training and testing

Because the main focus of this study was to understand the death pattern at the beginning of the pandemic in São Paulo, we divided the data into two sets: training and test. The model was trained using data from the training set and the test set was later used to estimate the model performance because we were interested in estimating the model accuracy on genuine forecasts, which refers to data that was not used when fitting the model. To validate and choose our models, we performed the classical holdout technique splitting our data following the 80-20 proportion: the training set consisted of death records from January 1, 2019 to October 19, 2019, and the test set from October 20, 2019 to January 1, 2020. We chose the best model considering its performance on the test set. Thereafter, we forecasted daily aggregated values from January 1, 2020 to March 20, 2020, and looked for abnormalities comparing the forecasts with the real values from the series. We chose the best model considering its performance on the test set.

Model fitting

The ARIMA model can be decomposed into two models: 1) autoregression (AR): in a linear regression model, a target variable is forecasted using a linear combination of covariates. In an AR, we forecast the target variable using a linear combination of past values of the target variable (¹⁵15. National Institute of Standards and Technology (U.S.). Croarkin C, Tobias P, Zey C. Engineering Statistics Handbook. Gaithersburg: The Institute; 2001.). The term p indicates the order of the AR model or the number of past values that are used; 2) moving average (MA): instead of using past values of the target variable, the MA model uses past forecast errors. In other words, the regression accounts for the current value of the series against a Gaussian white noise distribution (¹⁵15. National Institute of Standards and Technology (U.S.). Croarkin C, Tobias P, Zey C. Engineering Statistics Handbook. Gaithersburg: The Institute; 2001.). The term m indicates the order of the MA model or the number of error components used.

To incorporate both of these models, Box and Jenkins proposed a methodology that used stationary transformations to integrate both models. Thus, in the integration component in the ARIMA model, the term i indicates the number of differentiations performed to transform the non-stationary series into a stationary one. The task of estimating the parameters of the ARIMA or Box-Jenkins model can sometimes be subjective and difficult. To solve this, we first used the Hyndman-Khandakar algorithm (¹⁶16. Hyndman R, Khandakar Y. Automatic time series forecasting: the forecast package for R. J Stat Software 2008; 27.) to estimate the order of the model and we later performed a residual analysis to validate the model chosen by the algorithm.

The default behavior for the algorithm can be summarized in the following steps: a) because ARIMA modeling requires stationary time-series, the number of differences (i.e., obtained subtracting consecutive observations) is determined using repeated Kwiatkowski-Phillips-Schmidt-Shin (KPSS) tests; b) the values of p (the autoregressive component of the model) and q (the moving-average component of the model) are chosen by minimizing the small-sample corrected Akaike information criterion (AICc) after differencing the data, which includes a penalty term for complex models. The algorithm uses a stepwise search to transverse the solution space to find the optimal p and q parameters. This step consists of searching for the best model defined by the lowest AICc among the simple ARIMA models and making variations to the best model until no lower AICc can be found. To perform the residual analysis, checking if the model had adequately captured all the information in the data and can be used for prediction, we visually checked for correlations in the residuals using autocorrelation plots (ACF) and partial autocorrelation plots (PCF). We have also performed portmanteau tests to assert that there was no correlation using the Ljung-Box test. Furthermore, we calculated the mean value for the residuals, looking for a zero mean value, to assert that the model did not have any bias.

There are some challenges when dealing with long daily time-series data, such as dealing with multiple seasonal patterns and long-term seasonal patterns. To account for these effects, we tested two approaches: a) pure ARIMA modeling: components were estimated using the previously mentioned Hyndman-Khandakar algorithm; b) ARIMA modeling with time series decomposition: we performed decomposition of the series with the Seasonal-trend decomposition procedure with loess (¹⁷17. Cleveland R, Cleveland W, McRae J, Terpenning I. STL: A seasonal-trend decomposition procedure based on loess. J Offic Statist 1990; 6: 3–33.) and used an ARIMA modeling in the seasonal-adjusted component. The ARIMA estimated components were obtained using the Hyndman-Khandakar algorithm.

When evaluating these two approaches, we selected the best model in each one and compared the performance on the test set. The winning model was the one that had the best performance (i.e., the smaller error). When the errors were too close, we chose the model with a smaller number of degrees of freedom (i.e., the simplest model).

Statistical analysis

To estimate the accuracy of the model in genuine forecasts, we used the mean absolute error [MAE = mean(|et|)], a scale-dependent error that refers to the differences between the real and the forecasted death count. This was calculated using the forecasts obtained in the test set.

We used mean and standard deviation to describe continuous variables and percentages for categorical ones. Characteristics of individuals were compared using unpaired Student's t-test and chi-square tests, when appropriate. To calculate the annual mortality rates, we performed a population correction using the estimated total number of inhabitants in the city of São Paulo according to the State System of Data Analysis Foundation (¹⁰10. Fundação Sistema Estadual De Análise De Dados (SEADE). 2020. Available from: http://produtos.seade.gov.br/produtos/projpop/index.php.
http://produtos.seade.gov.br/produtos/pr... ). The alpha level was set to 0.05. Data analysis was performed using the R 3.6.3 software (¹⁸18. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. 2020. Available from: https://www.R-project.org/.
https://www.R-project.org/... ).

Results

Characteristics of the sample

In total, there were 17,011 autopsies conducted at the SPAS from 2019 to 2020 (54.5% men) with an average age of 68.7±15.6 years. Table 2 describes the baseline characteristics of the training, test, and 2020 datasets. Individuals were older adults, predominantly male, and white with elementary education. Diabetes was one of the most common risk factors followed by smoking habit. In our sample, coronary artery disease (ICD I21 and I25) was one of the most common causes of death comprising almost 25% of the cases for all three datasets (training, test, and 2020 dataset).

The cause of death distribution was statistically different for all three datasets, with an increase in the ‘Other' category and a decrease in most of the selected causes of death. When comparing our training and test datasets, we only found significantly different frequencies for education. For the training set and the 2020 dataset, we also found different distributions for education and for AIDS, epilepsy, and smoking habits (Table 2).

Temporal analysis

The average monthly mortality rate during the period was 110.17 per 100,000 individuals. In the period of one year (January 2019 to January 2020), the monthly all-cause mortality rate increased from 8.84/100,000 in 2019 to 9.41/100,000 in 2020 (a 5.99% increase).

For coronary artery disease and pneumonia, we found evidence of a 7-day seasonal pattern, but for respiratory-related issues the seasonal pattern was less identifiable resembling a 6-day seasonal pattern; for the all-cause mortality, we found evidence of a 365-day seasonal shape pattern.

When looking at daily data, the average number of daily deaths during the period (2019-2020) was 38.22±7.29 with 4.65±2.88 pneumonia deaths, 5.22±3.02 deaths due to respiratory diseases, and 8.83±4.25 by coronary artery disease.

Figure 1 shows the time-series plots and the autocorrelation functions (showing the seasonal patterns described above) for all-cause, pneumonia, respiratory-related, and coronary artery disease death cases.

Figure 1
Time-series and autocorrelation function (ACF) plot for all-cause mortality and specific causes (coronary artery disease, respiratory diseases, and pneumonia).

Model fitting and forecasting

For all subseries, all-cause and specific causes of death (coronary artery disease, respiratory disease, and pneumonia), the best model was the pure ARIMA model, giving the best out-of-sample mean absolute error (MAE) metrics. The results of portmanteau tests on the residuals showed that the pure ARIMA models had no correlations left on the residuals and, because of that, these were suitable for making predictions (Table 3).

The model for coronary artery disease indicated an excess of deaths for four dates in 2020: January 23rd, February 29th, March 8th, and March 9th. For pneumonia, the chosen model only indicated an excess of deaths on 2 dates: March 16th and March 20th. For respiratory diseases, the chosen model also only indicated an excess of deaths on 2 dates: March 16th and March 20th. Finally, the all-cause mortality series model revealed an excess of deaths on March 16th and March 20th. Figure 2 shows the forecast for each one of the series and their confidence intervals for the period of January 1 to March 20, 2020.

Figure 2
Left panel: Expected (forecast) and observed (real) number of deaths from January to March 2020. Right panel: Expected (forecast) and observed (real) number of deaths in March 2020.

Discussion

In this study, we investigated whether an excess number of deaths occurred in the SPAS during the first trimester of 2020 using forecasts from a time-series model trained on daily data from January 2019 to mid-October 2019. We focused on the all-cause mortality and selected coronary artery disease and respiratory ICDs related to COVID-19 according to multiple reports (¹³13. Vetter P, Vu D, L'Huillier A, Schibler M, Kaiser L, Jacquerioz F. Clinical features of covid-19. BMJ 2020; 369: 1470, doi: 10.1136/bmj.m1470.
https://doi.org/10.1136/bmj.m1470... ,¹⁴14. Zhang Y, Xiao M, Zhang S, Xia P, Cao W, Jiang W. et al. Coagulopathy and antiphospholipid antibodies in patients with Covid-19. N Engl J Med 2020; 382: e38, doi: 10.1056/NEJMc2007575.
https://doi.org/10.1056/NEJMc2007575... ).

For the all-cause mortality, we only found two dates with an excess of deaths (with 3 and 6 more deaths than expected for March 16th and 20th, respectively) from our model forecasted from previously observed values. In São Paulo, the health crisis due to the pandemic started around the end of April, when the ICU hospitalization rate for COVID-19 reached 67% and the number of deaths due to the virus reached 1,124 confirmed cases and 1,603 suspected cases (¹⁹19. Secretary of Health of the City of São Paulo. Boletim diário Covid-19 no município de São Paulo de 23 de abril de 2020 [Internet]. 2020. Available from: https://www.prefeitura.sp.gov.br/cidade/secretarias/saude/noticias/index.php?p=296672
https://www.prefeitura.sp.gov.br/cidade/... ). Therefore, we did not expect to find an excess in the number of deaths for the all-cause mortality related to the pandemic.

InfoGripe, a system from Fundação Oswaldo Cruz that monitors cases of hospitalization and presents alert levels for cases of severe acute respiratory syndrome (SARS) in Brazil, recorded a ten-fold increase in the historical average of hospitalizations on April 31, 2020 (²⁰20. Carvalho G. Casos de Síndrome Respiratória Aguda têm recorde histórico. Agência Fiocruz de Notícias [Internet; in Portuguese]. 2020. Available from: https://agencia.fiocruz.br/casos-de-sindrome-respiratoria-aguda-tem-recorde-historico.
https://agencia.fiocruz.br/casos-de-sind... ), that is, after the notification of the first case in Brazil on February 25, 2020. Because of that, we expected to find possible abnormalities or excess deaths due to pulmonary causes, especially for pneumonia and SARS. Nevertheless, for pneumonia, we only found an excess in the number of deaths on March 15th (1 more than expected) and March 20th (5 more than expected), the latter being the last day of our daily time-series and when the city of São Paulo started to apply public policies to contain the pandemic because of the first death on the 17th of March. The increase above the forecasted confidence intervals matched the overall health crisis that started in São Paulo, with social distancing policies being enforced to prevent more infections. Furthermore, we found an excess of deaths on the same dates due to respiratory disease, mainly driven by pneumonia. Finally, for coronary artery disease death cases, we found a few days with an excess in the predicted number of deaths by our model, but they lacked a constant increasing pattern that could be related to the COVID-19 pandemic.

Although we expected to see an abnormal increase in deaths on these series, especially for pneumonia and the respiratory ICDs, our results suggested that only in mid-March did deaths begin to increase above the confidence intervals of the models. This did not contradict evidence that COVID-19 reached Brazil in early-February, as reported in several studies (⁷7. Delatorre E, Mir D, Gräf T, Bello G. Tracking the onset date of the community spread of SARS- CoV-2 in western countries. Mem Inst Oswaldo Cruz 2020; 115: e200183, doi: 10.1590/0074-02760200183.
https://doi.org/10.1590/0074-02760200183... ,²¹21. Candido DS, Claro IM, de Jesus JG, Souza WM, Moreira FRR, Dellicour S, et al. Evolution and epidemic spread of SARS-CoV-2 in Brazil. Science 2020; 369: 1255–1260, doi: 10.1126/science.abd2161.
https://doi.org/10.1126/science.abd2161... ), but suggested that the mortality dynamics for the virus in Brazil were slower than other Western countries probably because of the warm months in the first trimester, which was also reported previously (²²22. Backer A. Slower COVID-19 case and death count growth at higher irradiances and temperatures. SSRN Electronic J 2020.).

There are some limitations to this study. First, our data (truncated at March 20, 2020) correspond to the beginning of the pandemic in Brazil, making it hard to make conclusions from March onwards. Additional observations after this date would allow more significant findings regarding the COVID-19 evolution in São Paulo. Second, due to the context of the pandemic, many autopsy protocols suffered changes to reduce contamination risks. In São Paulo, the first official resolution by the Mayor’s office of the City of São Paulo took place on March 20, 2020 (²³23. PRO-AIM/CEInfo. Resolução SS-32, 20-03-2020. São Paulo: City Hall of São Paulo; 2020. Available from: https://www.prefeitura.sp.gov.br/cidade/secretarias/upload/saude/arquivos/mortalidade/Resolucao-SS_32.pdf.
https://www.prefeitura.sp.gov.br/cidade/... ), establishing rules for the investigation of death cases including indirect and verbal autopsy procedures for suspected cases and expansion in the use of less invasive tools. It is worth mentioning that due to this resolution the SPAS stopped the conventional autopsy procedure on March 21, 2020. Third, the ARIMA modeling approach only accounts for one type of seasonality, which means that other seasonality patterns may not have been incorporated by the model and that could affect its performance since we are not using any covariates.

The analysis conducted in this article provides interesting insights into the dynamics of the pandemic during the first trimester of 2020, especially because data comes from a heterogeneous source, representative of the city of São Paulo (the largest city in Brazil). Also, because data come from autopsy reports, the cause of death is much more reliable than clinical data. Our study suggests that mortality and morbidity rates in warm regions had lower growth than in cold regions.

Future work could use hierarchical time-series modeling to accumulate effects from different death causes or using covariates accounting for different events that could have an impact on mortality (like the Brazilian Carnival). In conclusion, our findings are in line with official information that reported that the first death due to the coronavirus in Brazil occurred in mid-March since we did not observe an excess in the number of deaths before this period.

Acknowledgments

We would like to thank the São Paulo Autopsy Service from University of São Paulo for providing data access for this study.

References

Publication Dates

History