Generalized additive models to capture the death rates in Canada COVID-19

07/05/2020 ∙ by Farzali Izadi, et al. ∙ 0

To capture the death rates and strong weekly, biweekly and probably monthly patterns in the Canada COVID-19, we utilize the generalized additive models in the absence of direct statistically based measurement of infection rates. By examining the death rates of Canada in general and Quebec, Ontario and Alberta in particular, one can easily figured out that there are substantial overdispersion relative to the Poisson so that the negative binomial distribution is an appropriate choice for the analysis. Generalized additive models (GAMs) are one of the main modeling tools for data analysis. GAMs can efficiently combine different types of fixed, random and smooth terms in the linear predictor of a regression model to account for different types of effects. GAMs are a semi-parametric extension of the generalized linear models (GLMs), used often for the case when there is no a priori reason for choosing a particular response function such as linear, quadratic, etc. and need the data to 'speak for themselves'. GAMs do this via the smoothing functions and take each predictor variable in the model and separate it into sections delimited by 'knots', and then fit polynomial functions to each section separately, with the constraint that there are no links at the knots - second derivatives of the separate functions are equal at the knots.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

Code Repositories

Generalized-additive-models-to-capture-the-death-rates-in-Canada-COVID-19

None


view repo
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

In late December 2019, some severe pneumonia cases of unknown cause were reported in Wuhan, Hubei province, China. These cases epidemiologically linked to a seafood wholesale market in Wuhan, although many of the first 41 cases were later reported to have no known exposure to the market. In early January 2020, the World Health Organisation (WHO) named this novel coronavirus as Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) and the associated disease as COVID-19 [14]. Following this report, there has been a rapid increase in the number of cases as on 24th June 2020, there were over 9.2 million confirmed cases and almost 475,000 deaths worldwide, while the confirmed cases in Canada was 103,000, with 8,500 deaths mostly distributed in three provinces of Quebec, Ontario, and Alberta. The death rates for these cases which are available from https://www.canada.ca/en/public-health/services/diseases/2019-novel-coronavirus-infection.html

are shown in figure 1.

Figure 1. Reported deaths with COVID-19 in the Canada, Quebec, Ontario and Alberta since January 31, 2020 The Canada lock down started on March 17, 2020.

From the box plot one can see that there are a few outliers in the Alberta data only, see figure 2.

Figure 2. Note that there are some outliers in Alberta deaths.

Also the trend, seasonal and random components of each case can be obtained by decomposition function of time series where in figure 3 we plotted only the Canada deaths.

Figure 3. Decomposition of Canada death series into trend, seasonal and random components.

To see the normality of the data, we can plot the histogram as well as the density of the deaths, figure 4. The plots clearly show that there are substantial overdispersion relative to the Poisson for all cases except possibly for the Alberta data.

Figure 4. Histogram and the density of the deaths series. Note that the substantial overdispersion in all cases.

2. Generalized additive models

In general, generalized additive models (GAMs), see e.g. S. Wood, 2017, [10] are a semi-parametric extension of the generalized linear models (GLMs), used often for the case when there is no a priori reason for choosing a particular response function (such as linear, quadratic, etc.) and need the data to speak for themselves. GAMs first introduced by Hastie and Tibshirani, 1986 [7] and Hastie and Tibshirani, 1990 [8] and are widely used in practice (Ruppert et al. 2003) [6]; (Fahrmeir, L., Lang, S. 2001) [1]; (Fahrmeir, et al. 2004) [2] (Fahrmeir et al. 2013) [3]

. On the other hand, GLMs themselves are extension of the linear models (LMs). To understand the thing better, let us start with LMs. For a univariate response variable of multiple predictors one may simply write

(2.1)

It is clear that the response variable

is normally distributed with mean

, and variance

and the linearity of the model is apparent from the equation. One of the issues with this model is that, the assumptions about the data generating process are limited. One remedy for this is to consider other types of distribution, and include a link function relating the mean , i.e.,

(2.2)

In fact the typical linear regression model is a generalized linear model with a Gaussian distribution and identity link function. To further illustrate the generalization, we may consider a distribution other than the Gaussian for example Poisson or a negative binomial distribution for a count data where the link function is natural log function.

(2.3)

where is any exponential family distribution. We can still generalize more to add nonlinear terms in the above equation namely

where and an exponential or non-exponential family distribution and

are any univariate or multivariate functions of independent variables called smooth and nonparametric part of the model, which mean that the shape of predictor functions are fully characterized by the data as opposed to parametric terms that are defined by a set of parameters like the parameter vector

in the linear part in which both to be estimated.

Having said that, the ordinary least square problem generalizes into

(2.4)

where , controlling the extent of penalization and establishes a trade-off between the goodness of fit and the smoothness of the model.

On the other hand, these functions are represented using appropriate intermediate rank spline basis expansions of modest rank , as

(2.5)

Substituting the function under the integral sign with the above equation we get

(2.6)

where the right hand side is a quadratic form with respect to the known matrix .

Collecting both the parametric and non-parametric coefficients into a double parameter , we obtain

(2.7)

Writing we get the more compact form

(2.8)

Let be the maximizer of and the negative Hessian of at . From a Bayesian viewpoint is a posterior mode for . The Bayesian approach views the smooth functions as intrinsic Gaussian random fields with prior given by , where is a Moore–Penrose orpseudoinverse of . Furthermore in the large sample limit

(2.9)

Writing the density in (2.8) as , and the joint density of and as , the Laplace approximation to the marginal likelihood for the smoothing parameters and is . Finally, nested Newton iterations are used to find the values of and ; maximizing ) and the corresponding . Wood, S. N. 2016 [11]

3. Application of GAMs for fatal cases of Canada

Let denotes the deaths reported on day . We will study the Canada deaths in general and the provinces of Quebec, Ontario, and Alberta in particular cases. Due to overdispersion seen in the deaths, we assume that follows a negative binomial distribution with mean and variance

. Note that , where .

Then let

(3.1)

where is a link function log in our discussion, is a smooth function of time, , measured in days, is a zero mean cyclic smooth function of day of the week , set up so that , is a zero mean cyclic smooth function of day of the biweek , set up so that , is a zero mean cyclic smooth function of day of the month , set up so that , and denotes the order of the derivative in , and , see S. Wood 2020 [13]. Based on the discussion and notations in section.1, , , and are basis functions representing the underlying death rate and the strong weekly, biweekly and monthly cycles seen in the data respectively.

From (2.8) we can easily compute the confidence intervals for each f and make inferences about when the peak in f occurs. This can be done by executing gam function of mgcv library in R code. The fitted models to the reported deaths in Canada, Quebec, Ontario, and Alberta are shown in figures 5, 6, 7 and 8 respectively. In the Canada case, we model the deaths with respect to day, day of week and day of month. As the results of figure 5 shows, all the variables are statistically significant with

as R-squared.

Figure 5. Results of the fitted model to the Canada deaths. Note that all variable are statistically significant.

For the deaths of Quebec, we see that the day, day of week is an appropriate choice for the predicted variables. As the results of figure 6 shows, all the variables are statistically significant with as R-squared.

Figure 6. Results of the fitted model to the Quebec deaths. Note that all variable are statistically significant.

In the case of the Ontario we see that also the day, day of week and day of month is an appropriate choice for the predicted variables. As the results of Figure 7 shows, all the variables are highly significant with as R-squared.

Figure 7. Results of the fitted model to the Ontario deaths. Note that all variable are highly significant.

We figured out that the day, day of biweek as predictors for the Alberta case is an appropriate choice for the model. As the results of figure 8 shows, all the variables are statistically significant with as R-squared.

Figure 8. Results of the fitted model to the Alberta deaths. Note that all variable are statistically significant.

The next statistics is the results of the gam.check consists of Q-Q plot where for a good fit the data should lie on the red line. The second is the histogram of the residuals. In this plot, the histogram must be symmetric with respect to the line . Third plot is the Residuals vs. Linear predictors. This plot must be symmetric with respect to the line . Finally, the last plot is the Response vs. Fitted values. The more closer the data to the line , the better the fitted model. These plots for the corresponding data of Canada, Quebec, Ontario and Alberta are given in figures 9, 10, 11 and 12 respectively.

Figure 9. Canada: Q-Q plot, Histogram of residual, Residuals vs. Linear predictors and Response vs. Fitted values.
Figure 10. Quebec: Q-Q plot, Histogram of residual, Residuals vs. Linear predictors and Response vs. Fitted values.
Figure 11. Ontario: Q-Q plot, Histogram of residual, Residuals vs. Linear predictors and Response vs. Fitted values.
Figure 12. Alberta: Q-Q plot, Histogram of residual, Residuals vs. Linear predictors and Response vs. Fitted values.

The other statistics of the model fits to the reported deaths in all cases are shown in figures 13, 14, 15 and 16 respectively. The posterior modes (solid) and 95% confidence intervals for the model functions as well as auto-correlation functions and the deviance residuals against day for the Quebec and Alberta.

Figure 13. Posterior modes (solid) and 95% confidence intervals for the model functions for the Canada as well as auto-correlation functions of the deviance residuals. Lag 1 is just significantly negative.
Figure 14. Posterior modes (solid) and 95% confidence intervals for the model functions for the Quebec as well as auto-correlation functions of the deviance residuals and deviance residuals against day. Lag 7 is just significantly negative.
Figure 15. Posterior modes (solid) and 95% confidence intervals for the model functions for the Ontario as well as auto-correlation functions of the deviance residuals. Lag 7 is just significantly negative.
Figure 16. Posterior modes (solid) and 95% confidence intervals for the model functions for the Alberta as well as auto-correlation functions of the deviance residuals against day. No significant lag.

4. Inference about the peak of the day of each four cases

With gam model it is also straightforward to make inferences about when the peak in f occurs. To this end, it is enough to use the model matrix by removing cyclic part, the coefficients and variance-covariance matrix of the model to estimate the model functions and confidence intervals. To find the day of occurrence of the peak for each corresponding underlying death rate function, , simply we generates multivariate normal random deviates using ’rmvn’ function from ’mgcv’ package in which it takes 3 arguments as number of simulations, the coefficients and variance-covariance matrix of the model. The results for all 4 cases are shown on the figures 17 and 18 respectively.

Figure 17. Underlying death rates for the Canada(left) and Quebec (right) with confidence intervals. Day 0 is March 17th 2020 the Canada lock down. The scaled red bar chart illustrates the distribution of the location of the peak of the death rate curve, obtained by simulation from the approximate posterior for each model.
Figure 18. Underlying death rates for the Ontario(left) and Alberta (right) with confidence intervals. Day 0 is March 17th 2020 the Canada lock down. The scaled red bar chart illustrates the distribution of the location of the peak of the death rate curve, obtained by simulation from the approximate posterior for each model.

5. Inference about the past fatal infection cases

To obtain the model priors with auto-generated code and associated data is to simulate jagam from rjags library in GAMs. We also load glm to improve samplers for GLMs. This is useful for inference about models with complex random effects structure. The new mgcv function jagam is designed to be called in the same way that the modeling function gam would be called. That is, a model formula and family object specify the required model structure, while the required data are supplied in a data frame or list or on the search path. However, unlike gam, jagam does no model fitting. Rather it writes JAGS code to specify the model as a Bayesian graphical model for simulation with JAGS, and produces a list containing the data objects referred to in the JAGS code, suitable for passing to JAGS via the rjags function jags.model (Plummer 2016)[5]. To infer the sequence of past fatal infections one needs to produce the observed sequence of deaths. Verity et al. 2020 [9] show that the distribution of time from onset of symptoms to death for fatal cases can be modelled by a gamma density with mean 17.8 and variance 71.2 (s.d. 8.44). Let be the function describing the variation in the number of eventually fatal cases over time. Let be the lower triangular square matrix of order , describing the onset-to-death gamma density, given above where is the number of day of pandemic under consideration. Then we have if and 0 otherwise. Now , where is the expected number of deaths on day . Then can be represented using an intermediate rank spline, again with a cubic spline penalty. We can then employ exactly the model of the previous section but setting . The only difference is that we need to infer over a considerable period before the first death occurs where 15 days is clearly sufficient given the form of deaths data. After executing the jagam code for the data frame of deaths, day, day of week, day of biweek and day of month variables - note that we have used different combination of the above variables to see which one is appropriate of Canada deaths in general and Quebec, Ontario and Alberta in particular cases. Jagam code with Poisson family distribution produces a model containing all the priors. Then by adding the matrix

and a bit of extra regularization to the output we pass the model to jags.model function JAGS is Just Another Gibbs Sampler. It is a program for the analysis of Bayesian models using Markov Chain Monte Carlo (MCMC) simulation. Then by passing the result to the jags.samples function with the parameters of thin=300 and iteration=1000000 we get the values of

, , and the Monte carlo array of 3 dimension with the first dimension equals to the number of parameters in the gam model. As in the gam model, we take the data in the first component of the jags model as the model.matrix and the data in first 2 dimension of the to simulate the fatal infection profiles and get all the necessary statistics such as median with confidence intervals, the peak points of the median profile, the squared second difference of the median infection profile on the log scale which is proportional to the smoothing penalty and the absolute gradient of the infection profile. See Wood, S. N. 2016 [12] These are depicted in figure 19.

Figure 19. The inferred Canada fatal infection profile over time (day 0 is March 17th 2020). The red curve is the posterior median, the dashed curves delimit an confidence interval and the dotted curves a confidence interval. The scaled black bar chart shows the posterior distribution of the day of peak fatal infection. The dashed grey curve is proportional to the squared second difference of the median infection profile on the log scale which is proportional to the smoothing penalty. The green curve is proportional to the absolute gradient of the infection profile.

Acknowledgement. I am very grateful to Simon S. Wood as his paper [13] was very inspirational to my work.

References

  • [1]

    Fahrmeir, L., and Lang, S. 2001, Bayesian Inference for Generalized Additive Mixed Models Based on Markov Random Field Priors, Applied Statistics, 50, 201–220.

  • [2] Fahrmeir, L., Kneib, T., and Lang, S. 2004, Penalized Structured Additive Regression for Space-TimeData: ABayesian Perspective, Statistica Sinica, 14, 731–761.
  • [3] Fahrmeir, L., Kneib, T., Lang, S., and Marx, B. 2013, Regression Models, New York: Springer.
  • [4] Green, P. J. and B. W. Silverman 1994. Nonparametric Regression and Generalized Linear Models. Chapman & Hall.
  • [5] Plummer M 2016. rjags: Bayesian Graphical Models Using MCMC. R package version 4-6, URL https://CRAN.R-project.org/package=rjags.
  • [6] Ruppert, D.,Wand, M. P., and Carroll, R. J. 2003, Semiparametric Regression, New York: Cambridge University Press.
  • [7] Hastie, T., and Tibshirani, R. 1986, Generalized Additive Models (with discussion), Statistical Science, 1, 297–318.
  • [8] Hastie, T., and Tibshirani, R. 1990, Generalized AdditiveModels, London: Chapman & Hall.
  • [9] Verity, R., L. C. Okell, I. Dorigatti, P. Winskill, C. Whittaker, N. Imai, G. Cuomo-Dannenburg, H. Thompson, P. G. Walker, H. Fu, et al. 2020. Estimates of the severity of coronavirus disease 2019: a model-based analysis. The Lancet Infectious Diseases.
  • [10] Wood, S. N. 2017. Generalized Additive Models: An Introduction with R (2 ed.). Boca Raton, FL: CRC press.
  • [11] Wood, S. N., N. Pya, and B. 2016. Smoothing parameter and model selection for general smooth models (with discussion). Journal of the American Statistical Association 111, 1548–1575.
  • [12] Wood, S. N. 2016. Just another Gibbs additive modeller: Interfacing JAGS and mgcv. Journal of Statistical Software 75(7).
  • [13] Wood, S. N. Simple models for COVID-19 death and fatal infection profiles arXiv:2005.02090v1 [stat.AP] 5 May 2020.
  • [14] World Health Organization. www.who.int/emergencies/diseases/novel-coronavirus-2019. Accessed 15 Mar. 2020.