A predictive model for daily cumulative COVID-19 cases in Ghana

2.1. Polynomial regression model

Given the plot of the cumulative cases of the COVID-19 in Figure 1A, it is obvious that linear regression models

do not provide accurate statistical inferences since the relationship between the observed cumulative cases and time (in days) is non-linear. There is the need to modify the linear model to account for the non-linear relationship, by using polynomials of higher degree (i.e. degrees greater than one). In general, non-linear effects can be modelled by using polynomials of degree p defined as follows:²⁷

In the regression model (1), the parameters β_j are independent of the predictor variable x_i, however, the basis functions depend non-linearly on the predictor variable. Consequently, the parameters can be estimated using ordinary least squares approaches. In general, we can express model (1) in terms of a smooth function as follows:

where f represents a function or a transformation of the predictor variable x_i.²⁷

Polynomial models are easy to implement, however, their non-local property (i.e. the fitted function at any given value x₀ depends on data values that are far from x₀) is their major disadvantage. This issue can be avoided by dividing the domain of x into smaller intervals, fitting accurate polynomials in each interval and then finally combining the piecewise polynomial into a global one.²⁷ The domain of x is divided into smaller intervals using an arbitrary number/position of points τ known as knots.²⁷ A piecewise continuous model is fitted by specifying the following functions:

with + as a function defined by:

The combination of these sets of functions give rise to a composite function defined as f(x).

2.2. Spline regression

Polynomial regression does not capture the complete non-linear relationship. An alternative, and often superior approach for modeling non-linear relationships is the use of splines.⁷ A spline can be perceived as a flexible thin strip of wood or metal that can be used to draw smooth curves.²⁷ They require several weights to be placed at certain positions so that the strip of wood would bend according to the number/position of the weights.²⁷ Statistically, splines are used to reproduce flexible smooth curves.²⁷ That is, splines enable smooth interpolation between fixed points, called knots. They are series of polynomial segments strung together.⁷

Assume that the curve f(X) evaluates to a single value y for each set of predictors x, where x can be univariate or multivariate. If the set of knots is defined by τ₁ < τ₂ < … < τ_u in the domain of X,X ∈ R, then f(X) is a special polynomial of degree p, called a spline.

In modeling studies a smoothness criterion, which states that all derivatives of order less than p are continuous, is usually imposed.²⁷ A physical spline is linear beyond the last knot, thus, more constraints are imposed on derivatives of order 2 or greater at the leftmost and rightmost knots.²⁷ Splines which have these extra constraints are known as restricted or natural splines. Flexibility of the curves can be achieved by increasing the number of knots or the degree of the polynomial. However, it is worth noting that increasing the number of knots may lead to over-fitting due to associated high variances. Furthermore, decreasing the number of knots may lead to a rigid and restrictive function that has more bias.²⁷

Let f denote any spline function with a fixed knot sequence and a fixed degree p. Since the spline functions are objects in a vector space V, then f can be expressed as follows:

where the B_k are a set of basis functions spanning V and β_k are the associated spline coefficients.²⁷ For any k knots, there are k + 1 polynomials of degree p and p × k constraints. This leads to (d + 1)(k + 1) − p × k = d + k − 1 free parameters.^13,34 For natural or restricted splines, there are k free parameters. Since βB = (βA)(A⁻¹B) = δB* and for any non-singular matrix, there are an infinite number of possible basis sets for the spline.

The advantage of the equation (3) is that the estimation of f reduces to the estimation of the regression coefficients β_k. Specifically, the specification of Model (3) indicates that f is non-linear in the predictor but linear in the vector of regression coefficient β = (β₁,β₂,…,β_K+p+1). One can view the estimation of f as an optimization problem that is linear in the transformed variables B₁(X),…,B_K+p+1(X). Consequently, a framework is established for the estimation approaches to be adapted for splines in a wide range of generalized or multivariate regression.²⁷ A more appealing property of spline models is their ability to reduce the estimates to a few regression coefficients.²⁷

Although the flexibility property of splines makes them a better choice for fitting datasets, there are challenges associated with the number of tuning parameters.^13,27,34 That is, the choice of the basis functions B and the degree of the polynomial eventually have little impact. Sauerbrei and colleagues²⁷ noted that spline models are robust to the degree p of the polynomial. Polynomials with degree p = 3 (cubic polynomial) are standard because they are smooth curves. If the derivatives of the fitted curves are required, then a higher order polynomial is appropriate. However, the authors in²⁷ have observed that polynomial models with degree p > 3 are “effectively indistinguishable”.

Furthermore, modeling with splines involves deciding the number/spacing of knots and whether to use or not use a penalty function (the integrated second derivative of the spline). The absence of a penalty term in the spline model implies the generation of transformed variables which are added to the standard model. Such a procedure where the flexibility of the resulting non-linear function is entirely based on the number of knots is referred to as regression splines.²⁷ If the penalty term is added to the spline modeling, modification of the procedure is required to take into account the penalty term. In that case, each regression function has to be modified separately to obtain smooth splines that exhibit several desirable properties.

Moreover, a discussion on choices of basis B_k functions for splines can be found in.^{27,34, Chp. 5} The discussion here will involve B-splines and bases that are based on a special parametrisation of a cubic spline. These set of bases depend on the sequence of knots.^9,27 An advantage of the B-basis is that the bases have a local support. That is, the B-bases are larger than zero in intervals spanned by p + 1 knots and zero elsewhere.⁹ This property of the B-bases makes them numerically stable as well as present an efficient algorithm for building the basis functions.³⁴ Detailed information on different types of basis for splines and guidelines for the use of splines can be found in.²⁷

Futher, the selection and placement of knots is challenging due to the arbitrary nature of the task. That is, whenever a non-linear relationship is detected in data, the polynomial terms are not flexible enough to capture the relationship, however, splines require specification of the knots. GAMs provide a tool to automatically fit a spline regression.^16,17,34,35 GAMs will be discussed in the section that follows immediately.

2.3. GAM

The purpose of this section is to discuss GLMs and their extension to GAMs. The linear models (LMs) are used to model response variables that follow normal distributions whereas GLMs are used to model either normal or non-normal responses.²⁵ The general form of GLMs is:

where μ_i = E[Y_i], g is a smooth monotonic “link function”, X is n × p design matrix of covariates, X_i is the covariate associated with the i^th subject or item, β is a 1 × p vector of unknown parameters describing the effects of the covariates on the 1 × n matrix of responses Y_i, and n is the number of observations. The GLMs assumes that the responses Y_i are independent and follow some exponential family of distributions. The exponential family of distributions include Poisson, binomial, gamma, and normal distributions.²⁵ For a detailed discussion of GLMs, see.^{10, 23} Under the generalized linear mixed (GLMM) effects model, random effect components Z_ib_i are added to the fixed effect components X_iβ, where b_i is 1 × q is a vector of random effects and Z_i is a p × q design matrix of the random effects and $b_i\sim N\left(0,{\sigma }^2\right),{\sigma }^2$ is the variance of the random effect. So the general form of the GLMM is defined as:

GLMs are specified in terms of the linear predictor η = X_iβ which is the same as in the linear models. Hence, most of the concepts of linear modeling are maintained under the GLM framework with little modification. The formulation of the model is the same except that one has to choose the link function and the distributional assumption of the data. When data distribution is assumed to follow the normal distribution, the identity-link function is used and the GLM becomes the linear model for normal data. When data are counts such as number of new cases or number of cumulative cases, the appropriate distribution is the Poisson distribution with the log-link function option. When the outcome or response variable is binary, such as whether one is infected with the disease or not, then the appropriate distribution to assume is the binomial distribution with logit-link function.^10,23 A detailed discussion of exponential family of distributions and link functions can be found in.^{34, Section 3.1.1} Estimations of parameters and statistical inferences under the GLMs are based on the theory of maximum likelihood estimation. However maximization of the likelihood requires an iterative least squares approach discussed in.^{34, Section 1.8.8 (p. 54)} Also see^{34, Section 3.1.2} for detailed theory on fitting of the Generalized Linear Models.

A GAM is a GLM ith a linear predictor involving a sum of smooth functions of covariates.^16,17 The general form of the GAM is:

where μ_i = E[Y_i] and the response variable $Y_i\sim expf({\mu }_i,\varphi )$ where $expf({\mu }_i,\varphi )$ denotes an exponential family distribution with mean μ_i and scale parameter ϕ. The variable H_i represents a design matrix of covariates for any strict parametric model components, θ is a vector of parameter estimates describing the effects of the covariates on the response, f are the smooth functions of the covariates x_k. This model introduces flexibility in the specification of the response variable on the covariates.³⁴ However, complications are avoided when the model is specified in terms of “smooth functions” rather than detailed parametric relationships.³⁴ Simon N.³⁴ showed how GAMs can be represented using basis expansions for smooth functions, where each smooth function has an associated penalty controlling function smoothness. Estimation of parameters can be achieved by using penalized regression approaches. The appropriate degree of smoothness for f_j can be estimated from data using cross validation or marginal likelihood maximization.³⁴ For univariate smoothing, the representation and estimation of component functions of a model are best introduced taking into account a model consisting of a function of one covariate defined as:

where y_i is the response variable, x_i is the covariate, f is the smooth function and ϵ_i is random variable defined as ϵ_i ∼ N(0,σ²). Given equation (7), it is possible to represent a function with basis expansions. To estimate f, using the approaches applied to linear models,^{34, Chp. 1 and 3} it is required that f be represented such that the function (7) becomes a linear model. This can be achieved by selecting a basis that spans the space of functions of f or a close approximation to it. The chosen basis functions will be considered as completely known. That is, if B_j(x) is the j^th basis function, then f is assumed to have a representation defined by:

where β_j is a vector of unknown parameters. Substituting (8) into (7) yields:

which is a linear model.

Suppose that f is in the space of fourth order polynomials, then it follows that a basis for this space is

and the equation (8) becomes:

and the equation (7) becomes the following model:

In the case of additive models, suppose that there are two covariates, x and v, describing the changes in a response variable, y, then an additive model is defined as:

where β₀ is the intercept, f_j are the smooth functions, and ϵ_i are independent and identically normally distributed random variable with mean zero and variance σ². A notable issue is that the model now contains more than one function which leads to identifiability issue. It requires identifiability constraints to be imposed on the model before fitting. If the identifiability problem is addressed, then the additive model can be represented using penalized regression splines.^{34, P. 175, Section 4.3.1} The degree of smoothing is selected by cross validation or (RE)ML as done under the univariate model. Here the basis functions for f₁ are defined by using a sequence of k₁ knots with $x_j^{*}$ equally spaced over the domain of x and unknown γ_j coefficients. Also, the basis functions for f₂ are defined by using a sequence of k₂ knots with $v_j^{*}$ equally spaced over the range of v and unknown δ_j coefficients. It follows that

and hence

where the basis b_j(x_i) is the i,j elements of X₁. On the other hand,

and hence

where the basis B_j(x_i) is the i,j elements of X₂. For the identifiability problem, the best constraints according to Simon N.³⁴ are the sum-to-zero constraints:

where 1^′ is a 1 × n vector of 1s. This constraint does not change the shape of the smooth function f₁ but shifts f₁ vertically so that the mean value of f₁ is zero. For details on how this constraint can be applied and how additive models can be fitted using penalized least squares, see Simon N.^{34, P. 176-177}

The GAMs are extensions of additive models. Under the GAM framework, the linear predictors predict the known smooth monotonic function of the expected value of the response variable, where the response may follow any exponential family distribution. The linear predictor may simply have a known mean variance relationship which allows for the use of a quasi-likelihood methods. The GAM has the form described in equation (6), i.e.:

Detailed information on GAM theory can be found in.^{34, Chp. 6} The GAM is fitted using penalized likelihood maximization, which practically can be achieved by using penalized iterative least squares (PIRLS) described in.^{34, P. 180}

3.1. Data: COVID-19 Ghana cases

The data used in this study was obtained from the Ghana Health Service and the global cases from the Center for Systems Science and Engineering at Johns Hopkins University.³¹ The data shows that, as of February 24, 2021, the number of COVID-19 cases registered is about 80,759 with 582 deaths and 73,365 recoveries.³¹ The left panel of Figure 1 shows the monthly new cases of COVID-19 and the right panel of Figure 1 shows the trend of the number of cumulative cases from March 14, 2020 to February 28 2021. In general, the cumulative number of cases increased over the study period. The new cases registered peaked in July 2020 and then decreased until October 2020. The new cases continued to increase from November 2020 to February 2021 with a sharp increase in January 2021 and a slight decrease in December 2020. This continuous increase in the number of new cases is captured by the curve of cumulative cases.

The focus of this study is to determine an appropriate model that can be used to explain the dynamics or trend of cumulative cases and then predict/forecast cumulative cases of the virus for better management decisions. This requires the researcher to find a model that can fit the blue line the data points in the black curve. Statistical models in this work will be implemented for the number of cumulative COVID-19 cases. About 80% of the data was used as training data and the remaining 20% as test data to validate the models. The left panel of the Figure 2 represents the number of cumulative COVID-19 cases for the training dataset and test dataset are presented in the right panel of the Figure 2.

Figure 2. Plots of the training dataset (left panel) and test dataset (right panel) for cumulative coronavirus disease 2019 (COVID-19) cases.

3.2. Polynomial modeling of COVID-19 data

Firstly, a naive linear regression model to the cumulative COVID-19 cases in the left panel of Figure 1. The left panel of Figure 1 shows the curve of the linear regression model compared with the real data. This model provides the worst fit with highest RMSE = 6023.14, MAE = 5292.04, MAPE = 28.80654 and R² = 0.93. The R² indicates that 93% of the dynamics in the COVID-19 cases have been explained by time because of the general increase in the number of cases. However, the linear model does not capture non-linearity in the data leading to a very high RMSE. This is evident from the left panel of Figure 1, where the predictions of the fitted model (in the blue line) do not follow the observed trend of the COVID-19 cases. Best fit should approximately follow the observed trend shown by the black curve.

Next, a polynomial model with appropriate degree p is fitted on the cumulative COVID-19 training dataset in Figure 2 (left panel) and then applied on the test datasets, shown on the right panel of the Figure 2, to validate the model. Various polynomials defined by different degrees p were fitted and the polynomial model with degree p = 11 proves to produce the highest R² and lowest RMSE. The polynomial degrees beyond or below 11 are not significant. That is, polynomials with degree p < 11 produce the highest RMSE and lowest R² relative to polynomial with degree p = 11. On the other hand, polynomials with degree p > 12 lead to prediction with a rank-deficient fits. The curves of polynomial with degrees 3,7, and 11 are respectively shown in the top-right, bottom-left, and bottom-right panels of the Figure 3. The polynomial with degree 3 has the highest RMSE = 5297.00, MAE = 4914.50, and MAPE = 73.6702 giving the worst fit similar to that of the linear regression model. On the other hand, polynomials with degrees of 7 and 11 appear to provide accurate fits for the cumulative cases but polynomials with degree 7 have a very high RMSE = 1547.25, MAE = 1236.03 and MAPE = 19.67334 relative to RMSE = 591.2077, MAE = 484.9354, MAPE = 7.189698 of the polynomial with degree 11. The polynomial with degree 11 has the highest R² = 0.999 followed by the degree 7 polynomial (R² = 0.996) and degree 3 polynomial (R² = 0.947). The best fitting models from these models are the polynomial with degree 11 since it has the lowest RMSE and the highest R² value (see Table 1). In addition, this model also has the lowest AIC of 4383.862, whereas the polynomials with degree 3 and 7 have AICs of 5687.345 and 4955.772 respectively. Although the polynomial with degree 11 appears to capture the non-linearity in the data, it gives a very poor prediction. This is exhibited in the forecasts in Table 2 and the top-right panel of Figure 6. Although forecasts from the linear model suggest increasing cases (see the top-left panel of the Figure 6), the forecasts from day 1 to day 14 of March 2021 compared with the real data indicate that the linear models are inaccurate for the COVID-19 Ghana data (see Table 1).

Figure 3. Fits of the linear regression model (top-left panel), polynomial of degree 3 (top-right panel), polynomial of degree 7 (bottom-left panel), and polynomial of degree 11 (bottom-right panel) to the cumulative cases.

3.3. Spline modeling of COVID-19 Data

Again we fit a spline model with appropriate knots and degree of polynomial to the cumulative COVID-19 training dataset in in the left panel of Figure 2 and then use the test datasets in the right panel of Figure 2 to evaluate the fitted spline model. This checks the ability of the fitted spline model to capture and explain the non-linearity in the COVID-19 cases. This means that we have specify two parameters include the degree of polynomial and the location of the knots.⁷ Following⁷ example, we have to chose values between 0.20 and 0.95 quantiles as the knots. Choosing and placing three knots at the lower, median, and upper quartiles produced a very bad fit of the data. In fact, we need to identify at least 14 knots between 0.20 and 0.95 quantiles for placement rather knots at the lower, median, and upper quartiles in Bruce and Bruce’s example.⁷

The spline model with 3 knots or degrees of freedom (df) which poorly fit the data are shown in the top-left panel of Figure 4. We observed that knots of less than 14 do not provide a best fit for the data with relatively high RMSE. For instance, a spline fit with 3 knots in the top-right panel of Figure 4 and 8 knots in the top-left panel of Figure 4B poorly fit the data. However, knots greater than or equal to 14 provide the best fit of the data with relatively low RMSE and AIC, MAE, and MAPE as shown in Table 1. For example, the bottom-left panel of Figure 4 and the bottom-right panel of Figure 4, with knots 14 and 50 respectively, appear to provide the best fit for the data. The spline model provides predictions almost exactly the same as the original data. This is exhibited in the forecast values in Table 2, where forecasts and observed cases for 1 to 8 March are almost the same and showed a general increase in the covid-19 cases from 1 March to 31 March 2021 in the bottom-left panel of Figure 6. The green dots in Figure 7 show an increasing trend of the observed COVID-19 cumulative cases in March which support the forecasts produced by this model.

Figure 4. Fit of the spline regression models with 3 knots (top-left panel), 8 knots (top-right panel), 14 knots (bottom-left panel), and 50 knots (bottom-right panel) to the cumulative cases.

3.4. GAM for COVID-19

The gam function from mgcv package in R software was used to implement the GAM. The gam model formulation allows for the inclusion of smooth terms such as splines s() and tensor products te(). In the gam function, there are a number of options available for controlling automatic smoothing parameter estimation.³⁴

The left panel of Figure 5 presents the plot of the fitted GAM to the COVID-19 cumulative cases. It can be observed that the GAM is able to capture the non-linearity exhibited by the COVID-19 cases. The effect of time (Days) is estimated as a smooth curve with 8.98 degrees of freedom and the p −value associated spline term s(Days) is less than 0.05 which gives an indication that time in days has significant effect on COVID-19 cases. The effective degrees of freedom (edf) is approximately 9 indicating that polynomial of degree 9 can be used for predicting. The total degrees of freedom is 9.98. The right panel of Figure 5 shows the plot of partial residuals:

versus time (days). The right panel of Figure 5 shows that the estimated effect of days with a corresponding 95% confidence intervals is strictly Bayesian credible intervals^{34, P. 293} shown as dashed lines. The points where the confidence limits and the fitted curve pass through zero on the vertical axis are due to the identifiability constraints imposed to smoothen the time (Days) term. From the right panel of Figure 5, it can be observed that the partial residuals are uniformly scatted round the fitted curve. This gives an indication that the model describes the data well.

Figure 5. Fit of the generalized additive model (left panel) with the effect of time (days) estimated as a smooth curve with 8.98 degrees of freedom and a plot of days versus partial residuals (right panel).

The GAM model provides predictions similar to the original data. This observation is shown in the forecast values in Table 2, where forecasts and real data values from day 1 to day 8 of March 2021 are almost the same and show an increase in the COVID-19 cases (see the bottom-right panel of the Figure 6). The green dots in Figure 7 shows an increasing trend of cumulative COVID-19 cases in March which supports the forecasts produced by this model.

Figure 6. Plots of the observed data values (blue marker dots) and forecasts (green marker dots) of cumulative coronavirus disease 2019 (COVID-19) cases using linear regression (top-left panel), polynomial regression (top-right panel), spline regression (bottom-left panel), and generalized additive model (bottom-right panel).

Figure 7. Observed COVID-19 data used in analysis from March 14, 2020 to February 28, 2021 (blue marker dots) with forecast data from day 1 to day 31 of March, 2021 (green marker dots).

3.5. Forecasting of cumulative COVID-19 Cases

In this section, the most accurate polynomial, spline and GAM regression models are applied to forecast the number of cumulative COVID-19 cases for one month (from 1 March 2021 to 31 March 2021). Figure 6 presents plots of the forecasted cumulative COVID-19 cases from 1 March (353 days) to 31 March (383 days) 2021.