Data-driven tools for assessing and combating COVID-19 outbreaks in Brazil based on analytics and statistical methods

2.1 Materials

For smart analytics, it is needed to employ and optimize various freely accessible data sources to update and select databases [15, 16]. Initial data about the care of COVID-19 suspected patients, or patients who presented moderate symptoms of the disease, are recorded in the Family Health Units. These units consist of a public health network (Brazilian Unique Health System) composed of private clinics, ambulatory centers, and public hospitals located in different neighborhoods of the Brazilian cities, such as Recife. Also, healthcare systems for patients who use supplemental services (Brazilian Supplementary Health) provide data about severe disease symptoms. Another way of tracking and recording data are the call centers created to answer questions about the symptoms and guide the population in their care. Therefore, the Brazilian Ministry of Health implemented the Research Electronic Data Capture (REDCap), e-SUS Epidemiological Surveillance (e-SUS-VE), and Influenza Epidemiological Surveillance Information System (SIVEP-Gripe) platforms to report prospective suspected, probable, and confirmed COVID-19 cases as informed by public and private health services (primary and emergency care). The above-mentioned systems are interrelated on the COVID-19 website (infoms.saude.gov.br/extensions/covid-19_html/covid-19_html.html, accessed on 10 August 2021), which summarizes daily the aggregated counts from both platforms. In the case of Recife, Brazil, governmental and official data sources come from:

• The Center for Strategic Information on Health Surveillance, whose website is available at cievsrecife.wordpress.com (accessed on 10 August 2021).

• The Secretariat of Planning and Management from the Pernambuco State, whose website is available at www.seplag.ce.gov.br (accessed on 10 August 2021).

The Brazilian government and public data sources are available at:

• brasil.io/dataset/covid19/caso (accessed on 10 August 2021).

• dados.gov.br/dataset (accessed on 10 August 2021).

• github.com/wcota/covid19br (accessed on 10 August 2021).

The mobility indexes of the Google Community Mobility Reports (GCMR) are available at www.google.com/covid19/mobility (accessed on 10 August 2021). The GCMR collates data from those accessing Google applications with smartphones or handheld devices, allowing us to record the location history or traceability. In addition, data based on scientific articles reviewed by peers that report primary data as the gold standard for data inclusion provide further information cited in the references of this article.

To find additional details for each case or patient, the data can be augmented with medical reports, registry data, and captured reports (web scraping), primarily through news websites. The primary databases contain epidemiological records used for surveillance, evaluation, and research to address public questions. The administrative records (such as outpatient and hospital data anonymized and privatized) are employed for accounting and controlling the production of the services provided. The medical records contain privatized and anonymized clinical data on patients (in this work, these data were not used because they would have to be approved by the respective ethics committees before they could be used for any kind of implementation). The databases of the Brazilian Ministry of Health have demographic and socio-economic characteristics such as conditional cash transfer and housing programs. Combining these databases is employed to study social determinants of health and evaluate policies for epidemic control. Among the recorded variables to be used, we can cite: (i) daily cases reported; (ii) daily deaths reported; (iii) age; (iv) sex; (v) complete blood analyses (which is a type of blood test) on patients tested and other coronavirus or influenza pathogen diagnoses; (vi) isolation index of confirmed and suspected cases; (vii) number of hospitalizations; (viii) number of hospital beds; (ix) population density for the year 2020 by age groups and geographical mesoregions; (x) age distribution of confirmed COVID-19 infection rates; (xi) specific rates by age group; and (xii) intensive care unit (ICU) admissions for COVID-19 cases; among others.

The CASTLab, DADO, and MODINTERV platforms integrate several interconnected modules, including tools for data collection, data preprocessing, feature extraction, learning, classifying, making decisions, and data augmentation. These tools and their inter-dependencies are summarized in a simplified workflow sketched in Fig. 1, representing a framework for producing insights using smart analytics to combat outbreaks (SACO). From learning modules, it is possible, for example, to evaluate the epidemic curves via statistical methods, trends models, susceptible, exposed, infectious, recovered (SEIR) epidemic models, and Shewhart cohort chart, among others. Also, rates, averages, and another quantity are provided such that, if correlated with vaccination rates, hospital beds, population density, and under-reported cases, permit us to capture trends, produce epidemiological/clinical profiles, and identify possible risk factors for mortality of Brazilian COVID-19 patients under different aggregated levels, such as population segment, or region.

Fig. 1.Simplified flow diagram on smart analytics for Brazilian COVID-19 data. The process involves feedback of augmented data by using data cleaning, data transformation, and feature extraction by applying learning algorithms that can be employed to prioritize the risks factors according to the quantified judgment of a decision-maker.

The classification of epidemic curves by including covariate data, such as vaccination rates, hospital beds, population density, and under-reported cases, implies a certain risk level (high or low). Thus, a decision-maker can choose to implement interventions to reduce the number of fatalities while softening the economic impact of the epidemic. Implementing epidemic curves by using purely statistical and mathematical techniques that evaluate quantities, such as rates and averages, to capture the trends is quite challenging. In this sense, SACO aligns with the emerging data science domain that employs scientific methods, processes, algorithms, and systems to extract helpful knowledge to respond to the COVID-19 epidemic. The SACO intersects health, economy, public policy knowledge, planning, epidemiology, methodological development, and information technologies to collect, select, analyze, visualize, model, simulate, optimize, and report outbreak data of an epidemic in a smart and interconnected way.

The SACO used in Brazil is based on a wide range of approaches including, among others: (i) gathering data, database design and mobile technology [17]; (ii) statistical estimates with the maximum likelihood method [18, 19]; (iii) interactive data visualization [20, 21, 22]; (iv) geostatistics [23, 24]; (v) graph theory [25, 26, 27]; (vi) Bayesian statistics [28, 29]; (vii) mathematical and computational modeling [30, 31, 32, 33]; (viii) genetic analyses [34, 35]; (ix) evidence synthesis methods [36, 37, 38, 39]; and (x) statistical quality control charts adapted to monitor COVID-19 deaths and other indexes [40, 41].

After the initial entry of records in the platform, the database should undergo a consistency check and analysis utilizing complementary methodologies such as data anonymization, imputation, and data linkage to identify possible duplicate records and missing values effects. Using the number of confirmed COVID-19 cases and deaths reported up to 10 August 2021, we apply data-driven innovative epidemiology tools. Among them, we can mention the manual verification performed by data curators or the use of statistical and artificial intelligence tools [42], such as probabilistic record linkage, feature hashing, and principal component analysis [43, 44], among others [45, 46]. The epidemiological outbreak of COVID-19 can be investigated with the verified data, including descriptive mapping of occurrences over time and estimation of the main epidemiological parameters employing mathematical models. Specifically, we may gain intelligent insights from epidemic curves that represent the number of new cases or deaths per unit of time, based on the date of symptom onset [47].

We put here only some of the smart tools and statistical methods used by the CASTLab, DADO, and MODINTERV platforms as illustration. We utilize the moving average model to estimate the non-stationarity in the number of confirmed COVID-19 cases and deaths to identify the moments when there was a significant change in the time-series trend. Note the importance of determining the rhythm of contagion after a policy intervention (such as social distancing, quarantine, closing schools). Also, we use georeferenced standardized data on specific settlements or areas, typically obtained from administrative units such as neighborhoods or census sectors. We utilize choropleth charts as a thematic map in which areas are colored proportionally to a statistical variable that indicates an aggregate summary of a geographic characteristic within each area. These geographical records help to assess the spread of contagion in a macro view and make massive contact tracing between regions easier to monitor. Thus, when implementing successful protocols for containing the epidemic, such as social isolation and containment barriers (for example, vehicles and individuals’ flux where are not allowed circulation between different geographical locations), does the handling of the pandemic to be more efficient, that is, the number of deaths and contagions decreases.

The quantification of transmissibility during epidemics is essential to design and adjust public health responses. Epidemics can be measured by the effective or instantaneous reproduction number, that is, the average number of secondary cases caused by an infected individual. Knowledge of the reproduction number is key to understanding the dynamics of any infectious disease, and these should be reevaluated as the epidemic progresses in space and time. The instantaneous reproduction number, which considers the speed with which a disease spreads in a population, is an epidemiological parameter that quantifies the average of contagions that an infected person causes, assuming that disease transmissibility is constant in a window of time [6]. Thus, an effective quarantine restricts the free spread of the disease and maintains this parameter at values less than one to achieve a controlled spread [6]. With the implementation of COVID-19 control measures, the disease transmissibility is likely to change as well. Based on the surveillance data of COVID-19, we calculate the time-varying reproduction number, instantaneous reproduction number, and peak date across basic reproduction number, R${}_0$, namely, after interventions. This is conducted to evaluate its changing dynamics and determine the effectiveness of possible intervention strategies as vaccination policies. Note that R${}_0$ is the expected number of cases generated by one case in a population where all individuals are susceptible to infection.

The epidemic curves are helpful in many ways. They provide a simple visual sketch of demographic dynamics that may be used to assess the growth or decline of an outbreak [48] and evaluate the effect of some measures, such as non-pharmacological interventions. In addition, epidemic curves often form the raw material employed by various modeling techniques for monitoring and forecasting [49, 50, 51].

2.2 Statistical methods

Time-series model

We estimate the trend of the number of COVID-19 cases per day in Recife, Brazil, with a time-series model as follows. Let $\mathit{y}\left(\mathit{t}\right)$ be the number of confirmed cases at time $\mathit{t}$. Then, the trend of this series, $\mathit{m}\left(\mathit{t}\right)$ say, can be estimated through the moving average on day $\mathit{t}$ given by

with $\mathit{d}$ defined in Eqn. 1 being a constant that regulates the width of the data interval through which the average is calculated, controlling the degree of smoothing. Observe that, as d decreases, the moving average is more sensitive, so it may capture recent trends very well, although it might produce false alerts. Nevertheless, when $\mathit{d}$ is large, these false alerts are avoided, but the identification of current trends is more transient [6]. The choice of $\mathit{d}$ can be made according to the seasonality of the minor scale detected and considering that the confinement measures for the observation period were the weekends.

Epidemic model

When intervention measures (such as lockdown and social displacement) are put in place and a certain proportion of the population gains immunity, interest switches to knowing the time-varying effective reproduction number. Consider an individual who turns infectious on day $\mathit{t}$. We denote by ${\mathit{R}}_e\left(\mathit{t}\right)$ the expected number of secondary cases that this infectious individual causes. We assume that the times of infecting others and detecting this infectivity (for example, by symptoms or by a test) coincides, that is, on day $\mathit{t}$, this person also appears in the incidence time-series. The time interval between symptom onset in the primary case and symptom onset in the secondary case is called the serial interval (SI), whose associated statistical model is called the SI distribution. Note that the SI is different from the generation time (GT), which is the period between exposure of the primary and secondary cases. However, since exposure is rarely observable, we must employ the time series of incident symptom onsets as a basis for the statistical inference [52]. Using a renewable equation, we can estimate the time-varying effective reproduction number [53]. Here, we consider a simple growth model and denote by ${\mathit{y}}_{\mathit{e}}\left(\mathit{t}\right)$ the expected number of new symptom onsets that we observe on day $\mathit{t}$.

Let $({\mathit{g}}_1,\dots ,{\mathit{g}}_M)$ be the probability mass functions of the SI distribution, that is, $\mathit{P}\left(\mathit{i}\right)={\mathit{g}}_{\mathit{i}}$, for $\mathit{i}\in \{1,\dots ,\mathit{M}\}$. The expected number of cases can be described by the homogeneous linear difference equation stated as

where $\mathit{t}\in \{2,3,\dots \}$}, and we ignore the terms when $\mathit{t}-\mathit{M}\le 0$.

We can estimate ${\mathit{R}}_{\mathit{e}}\left(\mathit{t}\right)$ defined in Eqn. 2 by the method proposed by Wallinga and Teunis (W&T) [54, 55, 56]. The probability that case $\mathit{i}$ infected by case $\mathit{j}$, ${\mathit{p}}_{\mathit{i}\mathit{j}}$ namely, given its difference in time of symptom onset at $({\mathit{t}}_{\mathit{i}}-{\mathit{t}}_{\mathit{j}})$, may be expressed in terms of the SI distribution, denoted by $\mathit{w}$, which is formulated as

Thus, ${\mathit{p}}_{\mathit{i}\mathit{j}}$ stated in Eqn. 3 is the probability that case $\mathit{i}$ has been infected by case $\mathit{j}$, normalized by the probability that case $\mathit{i}$ has been infected by any other case k. The effective reproduction number for case $\mathit{j}$ is the sum of all cases $\mathit{i}$, weighted by the relative probability that case $\mathit{i}$ has been infected by case $\mathit{j}$, that is, ${\mathit{R}}_{\mathit{j}}=\sum _i^{}{\mathit{p}}_{\mathit{i}\mathit{j}}$. Therefore, ${\mathit{R}}_{\mathit{e}}\left(\mathit{t}\right)$ is obtained with the W&T method by averaging the individual reproduction number over all cases detected on the same day j, where t${}_{\mathit{j}}$ = t. The analysis can be carried out by utilizing the R0 package of the R software. Another way to evaluate ${\mathit{R}}_{\mathit{e}}\left(\mathit{t}\right)$ is to use the Robert Koch Institute (RKI) method, where a constant generation time of four days is assumed [57]. The RKI method estimates the time from infection to: (i) first symptoms to be around five days; (ii) being infectious around three days (which results in two days during which an affected person is infectious but might not suspect he/she is sick him/herself); and (iii) infecting other people around four days.

The instantaneous reproduction number $\mathit{R}\left(\mathit{t}\right)$ is interpreted as the average number of secondary cases that each symptomatic individual at time $\mathit{t}$ would infect if the conditions remained as they were at time $\mathit{t}$. Note that $\mathit{R}\left(\mathit{t}\right)$ is evaluated by comparing the number of new infections on day $\mathit{t}$ with the infection pressure from the days prior to $\mathit{t}$, that is, by the expression established as

Note that $\mathit{R}\left(\mathit{t}\right)$ given in Eqn. 4 is different from the basic reproduction number, which is presented next. Unlike the W&T forward-looking method, the estimate defined in Eqn. 4 is backward in time, whereas the timing of $\mathit{R}\left(\mathit{t}\right)$ can make a big difference when comparing it with intervention measures. The estimator of the instantaneous reproduction number is implemented in an R package named EpiEstim. If $\mathit{R}\left(\mathit{t}\right)$$>$ 1, it means that a case infects more than one person, resulting in a spread of the virus, whereas $\mathit{R}\left(\mathit{t}\right)$$<$ 1 translates into the containment of the virus, as one case infects less than one person.

We use a deterministic SEIR model (one-wave) to describe the spread of the disease because of the uncertainty due to the incomplete identification of the infected population. Note that, at the beginning of the pandemic, there was no clear policy on testing the population for SARS-CoV-2, and there were no test kits that would allow the correct identification of the actual infected cases. As mentioned, the SEIR model classifies the population into four classes: susceptible (S), exposed (E, infected but not yet infectious), infected (and infectious, I), and recovered (R). The dynamics of the transitions between the four different compartments (S, E, I, and R) is described by the non-linear system of ordinary differential equations (ODEs) given by

where the elements defined in Eqn. 5 correspond to $\mathit{N}=\mathit{S}+\mathit{E}+\mathit{I}+\mathit{R}$, $\sigma$ is the loss of latency rate, $\gamma$ is the recovery rate, $\mu$ is the disease-induced mortality rate, and $\beta$ is the transmission rate. We assume that those non-pharmacological interventions (for example, the wearing of face masks, social distancing, or self-isolation when sick) can decrease the number of contacts per individual per unit of time during the epidemic. To accommodate such effect, we consider that $\beta$ change in time $\mathit{t}$ by the dynamic formulated as

where ${\mathit{t}}_0$ is the starting time of the intervention (initial day) and $\tau$ is the average duration of the interventions. Furthermore, ${\beta }_0$ is the initial transmission rate and the product ${\beta }_0{\beta }_1$ represents the transmission rate at the end of the epidemic. The explicit solution to the equation given in Eqn. 6, subject to the constraint to $\beta \left({\mathit{t}}_{\mathit{0}}\right)={\beta }_0$, is stated as

where the elements defined in Eqn. 7 are established in the previous expressions. Once the SEIR model without the time-dependent transmission rate is fitted, the basic reproduction number is estimated [58, 59] as ${\mathit{R}}_0=\mathit{r}/(\gamma +\mu )$. Here, r is the exponential rate of increase of cases at the beginning of the epidemic, which, in turn, is estimated by fitting an exponential function employing only the first data. The estimation of parameters is based on a nonlinear curve fitted by minimization of the sum of square errors (by utilizing the least squares method). To solve the ODEs given in Eqn. 5, we use the ode function of an R package named dSolve. This package has implemented stiff and nonstiff integration routines utilizing the ODE-PACK Fortran codes (LSODE, LSODES, LSODA, LSODAR), DVODE, and DASPK2.0. A suite of Runge-Kutta integrators and special-purpose solvers to efficiently integrate 1-, 2- and 3-dimensional partial differential equations are available. The dSolve package includes fixed and adaptive time-step explicit Runge-Kutta solvers and the Euler method. We employ the optim function of the R software that has implemented the Nelder-Mead, quasi-Newton and conjugate-gradient algorithms.

We use curves for the cumulative number of COVID-19 deaths as a function of time. Note that they can show three distinct regions of behavior in the first wave of epidemics. There is an initial phase in which the epidemic curve is typically exponential, based on the number of deaths observed on the first day that a death occurred. This rapid rise is characterized by the growth rate $\mathit{r}$. Note that contagion may result in a complex network of human contacts or when mitigation measures are adopted at the beginning of the epidemic. However, depending on the complexity of the contagion dynamic, the growth may happen more slowly than an exponential curve characterized by a parameter $\mathit{q}\in (0,1)$, which interpolates between the linear regime $\mathit{q}=0$ and the exponential regime $\mathit{q}=1$. This type of epidemic curve presents an inflection point, denoted by ${\mathit{t}}_{\mathit{c}}$ and shown in Fig. 2, which corresponds to the instant in which the accumulated curve changes its concavity, that is, the growth rate reaches its maximum value and begins to decrease after that. The final part of the curve, after the inflection point, can be characterized by an $\alpha$ parameter that controls how quickly the epidemic curve moves away from the linear trend and bends towards the plateau. The plateau value, denoted by the parameter $\mathit{K}$, represents the total number of deaths at the end of the epidemic.

Fig. 2.Qualitative outline of an epidemic curve (black) for the accumulated number of deaths, with an indication of their different phases in the first wave of the pandemic, with the red circle indicating the change point of the curve.

We may model the epidemic curves in a deterministic way through a generalized Richards model (GRM) [52, 53]. Here, the GRM is defined by the ODE formulated as

where $\mathit{C}\left(\mathit{t}\right)$ is the cumulative number of deaths at $\mathit{t}$, $\mathit{r}$ is the growth rate in the initial phase, as mentioned, $0\le \mathit{q}\le 1$ is the parameter that enables us to interpolate among linear ($\mathit{q}=0$), sub-exponential $\mathit{q}<1$, and exponential ($\mathit{q}=1$) growth. The ODE stated in Eqn. 8 must be supplemented with an initial condition given as

for a given value of ${\mathit{C}}_0$. The exact solution of the expression established in Eqn. 8 for $0\le \mathit{q}<1$, subject to the constraint shown in Eqn. 9, can be written implicitly as

where ${}_2^{}{\mathit{F}}_1\left(\mathit{a}\mathit{,}\mathit{b}\mathit{;}\mathit{c}\mathit{;}\mathit{x}\right)$ is the hypergeometric Gauss function. The constant ${\mathit{t}}_{\mathit{i}}$ given in Eqn. 10 is determined by the initial condition ${\mathit{C}}_{\mathit{0}}$ through the relation defined as

The inflection point, ${\mathit{t}}_{\mathit{c}}$ namely, of the curve $\mathit{C}\left(\mathit{t}\right)$ is defined as the time when the second derivative of $\mathit{C}\left(\mathit{t}\right)$ with respect to $\mathit{t}$ is zero, that is, $\ddot{C}\left({\mathit{t}}_{\mathit{c}}\right)$, where the two points on the top of $\mathit{C}\left(\mathit{t}\right)$ indicate the second derivative with respect to the time. Using the expression formulated in Eqn. 11, we find that

where the elements defined in Eqn. 12 are established in the previous expressions.

The solution of the GRM given in Eqn. 8 is attained in an implicit form. It does not represent any numerical difficulty since the solution can easily be obtained, for empirical curve fitting purposes, by simply treating the data in the same implicit form, that is, as a curve of type $\mathit{C}\left(\mathit{t}\right)$. In the case where $\mathit{q}=1$, the GRM defined in Eqn. 8 reduces to the standard Richards model (RM), which has an explicit formula in terms of elementary functions stated as

For fitting the formula defined in Eqn. 13 to empirical data, it is convenient to set $\mathit{C}\left(0\right)={\mathit{C}}_0$, where ${\mathit{C}}_0$ is the number of deaths recorded on the first day that a death was reported. Using that $\mathit{C}\left(0\right)=\mathit{K}/{(1+\alpha {\mathit{e}}^{-\alpha \mathit{r}{\mathit{t}}_{\mathit{c}}})}^{1/\alpha }$, we can eliminate ${\mathit{t}}_{\mathit{c}}$ in favor of the other parameters, so we work with only three free parameters, namely $\mathit{r},\alpha$, and $\mathit{K}$, to be numerically determined. All statistical fittings may be performed using the Levenberg-Marquardt algorithm to solve the corresponding non-linear least-square optimization problem.

When estimating the model parameters from short time-series, we have overfitting problems due to the redundancy of parameters. This may lead, for example, to the estimation of specific parameters that are outside of biologically reasonable ranges or of other types. For instance, when applied to the number of infected cases in an epidemic, the parameter $\alpha$ should be constrained to the interval $(0,1)$. We use the GRM instead of the number of deaths, but we assume that the same constraint should be considered. For our purposes, we need to consider the constraints $0<\mathit{r}<1$;and $0<\alpha <1$ as empirical criteria for validating the GRM, which is unsuitable when the available data do not encompass the inflection point ${\mathit{t}}_{\mathit{c}}$. Nonetheless, as more data are considered, the model is expected to become more accurate. As an empirical criterion, we consider here that the GRM is only acceptable if ${\mathit{t}}_{\mathit{c}}$ is less than the time of the last data point. Hence, an intervention strategy (for example, adding covariate data, such as vaccination rates, hospital beds, population density, under-reported cases) may be used to model the data by assuming that its next result is captured by a change in the values of the parameters of the GRM after a given time.

Control charts

Looking at monitoring the number of deaths, adapted Shewhart control charts can be strategically utilized to differentiate variation due to the disease effect and variation due to common causes. Adaptation of these control charts is needed given the exponential behavior of indicators such as COVID-19 death rates. The use of a control chart for monitoring COVID-19 was introduced in [40, 41].

In the context of an epidemic, such as COVID-19, it is possible to analyze the evolution of indicators in a manner analogous to that of a production process. The idea behind the method assumes that observations are made periodically. This is equivalent to having sample sizes ${\mathit{n}}_{\mathit{t}}=1$ analyzed utilizing $\mathit{Q}\left(\mathit{t}\right)$ or some statistic indicator as the mean, standard deviation, or percentage. Its choice depends on the characteristic measured in each item and is calculated to summarize the values measured in the sample, $\mathit{S}\left(\mathit{t}\right)$ namely, representing the behavior of the process at time $\mathit{t}$. These statistics are compared to the monitoring limits such as central line (CL) as well as lower (LCL) and upper (UCL) control limits plotted over time. In the case of the indicator of COVID-19 deaths per day, the signal is not constant, and, in general, its growth is exponential. Then, it is possible to model the signal and what must be checked, in terms of stability, is the behavior of the residuals of the fitted model. If the growth is modeled by an exponential function described as

where $\mathit{y}\left(\mathit{t}\right)$ stated in Eqn. 14 is the number of deaths registered on day $\mathit{t}\in \{1,\dots ,\mathit{M}\}$, $\mathit{a}$ is the initial value of the growth model, and $\mathit{b}$ is the growth factor. Here, $\mathit{a}$ and $\mathit{b}$ are constants to be estimated. Hence, ${\log }_{}\left(\mathit{y}\left(\mathit{t}\right)\right)$ is a linear function on $\mathit{t}$ formulated as

The procedure for adapting the Shewhart chart to our context of COVID-19 is described as follows in Algorithm 1.

Algorithm 1: Adapted Shewhart control chart

Step 1. Fit ${\log }_{}\left(\mathit{y}\left(\mathit{t}\right)\right)$ defined in Eqn. 15 as a function of time $\mathit{t}$ via the linear model stated as

where ${\beta }_0$, ${\beta }_1$ are regression coefficients, and ${\epsilon }_{\mathit{t}}$ is a random error of zero expected value and constant variance.

Step 2. Calculate the residuals of the model stated in Eqn. 16 as

Step 3. Obtain the limit of the Shewhart chart drawn up for the residuals defined in Eqn. 17 as

where ${\hat{\sigma }}_{\mathit{Q}}$ is a measure of the variance of the statistic $\mathit{Q}\left(\mathit{t}\right)$.

Step 4. Calculate the LCL and UCL of the Shewhart adding and subtract ${\mathit{L}}_{\mathit{r}\mathit{e}\mathit{s}\mathit{i}\mathit{d}\mathit{u}\mathit{a}\mathit{l}}$ stated in Eqn. 18 to the values fitted by the regression considering

Step 5. The CL and the respective control limits stated in Eqn. 19 of the fitted Shewhart chart are calculated by exponentiating the fitted values and the limit calculated in step 4.

The implementation of Algorithm 1 must be used the logarithm in base 10. In addition, it is possible that, in some situations, null values are observed, that is, $\mathit{y}\left(\mathit{t}\right)=0$, for some $\mathit{t}\in \{1,\dots ,\mathit{M}\}$. Therefore, for implementing the adapted Shewhart control chart, $(1+\mathit{y}\left(\mathit{t}\right))$ must be employed as the response variable.