--- category: literaturenote citekey: abdellapenalizedconwaymaxwellpoisson2019 title: "Penalized Conway-Maxwell-Poisson regression for modelling dispersed discrete data: The case study of motor vehicle crash frequency" authors: "Abdella, Galal M.; Kim, Jinho; Al-Khalifa, Khalifa N.; Hamouda, Abdel Magid" year: 2019 date: 2019-00-00 2019 doi: 10.1016/j.ssci.2019.06.036 publication: Safety Science url: "https://doi.org/10.1016/j.ssci.2019.06.036" zotero_key: DEHVSRAU zotero_storage: 8IUMGZV6 collections: transpordi_planeerimine folder: 001_artiklid/Road Safety Indicators firstAuthor: "Abdella, Galal M." status: converted --- ![](_page_0_Picture_1.jpeg) ### Contents lists available at [ScienceDirect](http://www.sciencedirect.com/science/journal/09257535) ## Safety Science journal homepage: [www.elsevier.com/locate/safety](https://www.elsevier.com/locate/safety) ![](_page_0_Picture_5.jpeg) # Penalized Conway-Maxwell-Poisson regression for modelling dispersed discrete data: The case study of motor vehicle crash frequency ![](_page_0_Picture_7.jpeg) Galal M. Abdella[a,](#page-0-0)[⁎](#page-0-1) , Jinho Kim[b](#page-0-2) , Khalifa N. Al-Khalifa[a](#page-0-0) , Abdel Magid Hamoud[aa](#page-0-0),[c](#page-0-3) - a *Department of Mechanical and Industrial Engineering, Qatar University, Doha, Qatar* - b *Rutgers Business School, Rutgers University, NJ, USA* - c *Qatar Transportation and Traffic Safety Center, Qatar University, Doha, Qatar* #### ARTICLE INFO #### *Keywords:* COM-Poisson regression Motor vehicle crash rate Poisson regression Generalized linear models ABSTRACT Statistical modelling of road crashes has been of extreme interest to researchers over the last decades. Such models are necessary for the investigation of the opportunities for road safety improvement. The motor vehicle crash frequency (MVC-F) is probably the most important count of road crashes. In practice, like many of other discrete variables, this count is often diagnosed with over- or underdispersion, i.e. the variance is greater or less than the mean. The traditional regression models, especially those based on the Poisson distribution, are inefficient in modelling dispersed count data. On the contrary, the Conway-Maxwell-Poisson (COM-Poisson) distribution has been proven powerful in modelling count data with a wide range of dispersion. In crash data modelling, many situations may give rise to collinearity between contributory crash factors. Under this situation, the maximum likelihood estimates of the coefficients of the COM-Poisson GLM become increasingly unreliable as the collinearity among the model predictors increases. This paper addresses this issue and proposes a penalized likelihood scheme to be used with the COM-Poisson GLM regression for improving its prediction performance. For better GLM regression output, we suggest implementing the penalized COM-Poisson GLM regression under a K- fold cross-validation framework. A real-world crash example is provided, showing the performance of the penalized COM-Poisson GLM regression compared to the Poisson and the classical COM-Poisson GLM regressions. ### **1. Introduction** Statistical models for road crash data have become of extreme importance for all aspects of transport and traffic systems (TTS). These models, in particular, are helpful in establishing an evidence-based framework through which transportation authority can develop innovative road safety strategies and interventions. Over the past years, there has been a growing interest in literature on the development of statistical models for road safety studies ([Lord and Persaud, 2000; Lord](#page-6-0) [et al., 2005; Lord and Mannering, 2010; Kokonendji, 2014; Frank and](#page-6-0) [Samia, 2015; Paraskevi et al., 2015; Abdella et al., 2017; Abdur Rouf](#page-6-0) [et al., 2018](#page-6-0)). The Poisson distributions are the most popular for road safety management and analysis. The motor vehicle crash frequency (MVC-F) is probably the most common crash count variable ([Nazha et al., 2018](#page-6-1)). Since the MVC-F is usually characterized by over- and underdispersion, the use of the traditional Poisson-based models becomes inappropriate. The over-dispersion phenomenon usually occurs when the variance of the crash count exceeds the expected mean, i.e. *V Y E Y* ( ) ( ) > . On the contrary, the under-dispersion occurs when the variance is less than the expected mean, i.e. *V Y E Y* ( ) ( ) < . The over-dispersion in the MVC-F count usually occurs when the dataset on some road segments or timeintervals exhibits a large number of zero values ([Vangala et al., 2015](#page-6-2))*.* Several significant contributions have been conducted over the last years on the applications of the Conway-Maxwell-Poisson (COM-Poisson) in modeling over- and underdispersion count variables. [Guikema and Coffelt \(2008\)](#page-6-3) developed a COM-Poisson GLM formulation using two different links (dual-link) functions for modeling discrete counts. These two functions are = + = *ln µ*( ) *x i i* 0 *i* 1 and = + = *ln*( ) *z j q i* 0 *i* 1 , where *p* and *q* are the covariates used in the location and the shape link functions. [Lord et al. \(2008\)](#page-6-4) modified the dual-link formulation by replacing the second link function by the value of the shape parameter . The authors used log( ) 1 to link the response variable with the model predictors and then suggested using the *E-mail addresses:* [gmg5005@qu.edu.qa](mailto:gmg5005@qu.edu.qa) (G.M. Abdella), [jinho.kim@rutgers.edu](mailto:jinho.kim@rutgers.edu) (J. Kim), [Alkhalifa@qu.edu.qa](mailto:Alkhalifa@qu.edu.qa) (K.N. Al-Khalifa), [Hamouda@qu.edu.qa](mailto:Hamouda@qu.edu.qa) (A.M. Hamouda). ⁎ Corresponding author. Bayesian method for estimating the coefficients of the resulted model. The selection of the Bayesian estimation here is mainly attributed to its favourable mathematical and practical properties. The study has shown that the COM-Poisson GLMs perform the same as the NB models. One advantage of the COM-Poisson GLM over the NB model is that it can handle both the over- and underdispersed count data [\(Lord et al.,](#page-6-4) [2008\)](#page-6-4). This feature would, for sure, expands the domain of usage of the COM-Poisson-based models to include a wide range of industrial and service sectors. [Sellers and Shmueli, \(2010\)](#page-6-5) considered the log linear form of the Poisson regression by [McCullagh and Nelder, \(1997\)](#page-6-6) and proposed a COM-Poisson GLM formulation using a ( ( )) ( ) *E Y log* = as a link function. This function is often preferable over others for it would lead to a direct estimation of the coefficients of the GLM regressors. Both of the above studies are mainly based on the assumption of having a fixed value of the shape parameter, i.e. the shape parameter is independent of the covariates of the model. [Sellers and Andrew \(2016\)](#page-6-7) developed another format of the COM-Poisson, named as a zero-inflated COM–Poisson (ZICMP), having the capability of modelling both overand under-dispersion dataset. The ZICMP has revealed a performance relatively close to the zero-inflated negative binomial (ZINB) with respect to the Akaike Information Criterion (AIC). In the practice of road crash modeling, there are many situations may give rise to multicollinearity –or simply collinearity between contributory crash factors. In statistics, the collinearity refers to the existence of a significant linear relationship between the model factors. The collinearity has a severe impact on the reliability of the GLM estimates as well as the interpretation of the results. Despite the attractive prediction performance of the COM-Poisson-based regression models, we believe that there is still a need for developing COM-Poisson regression models that are capable of not only accommodating the dispersion issue but also being insensitive to the potential collinearity between road crash factors. This work develops a penalized likelihood formulation to be used in the COM-Poisson GLM regression to stabilizing the estimates of the regression parameters. We develop the penalized likelihood function by adding an *L*2 norm to the squared norm of the GLM coefficients in the likelihood function by [Sellers and Shmueli \(2010\)](#page-6-5). Considering the effectiveness of the COM-Poisson distribution in modelling dispersed count data, our penalized COM-Poisson GLM regression would provide a practitioner with a powerful tool for simultaneously overcoming both the "dispersion" and "collinearity" problems. In order to achieve a better predictive performance of our model, we implement the penalized COM-Poisson GLM regression under a K *fold*-cross-validation framework. We organize the rest of this research paper as follows: [Section 2](#page-1-0) introduces the COM-Poisson distribution. [Section 3](#page-1-1) describes the penalized COM-Poisson GLM regression model. [Section 4](#page-2-0) describes the cross-validation approach. [Section 5](#page-3-0) illustrates the implementation of the penalized COM-Poisson GLM for modeling a real-world example. [Section 6](#page-5-0) presents a summary of the findings and outlines extensions for future research work. #### **2. The COM-Poisson distribution** The COM-Poisson is usually considered as a general form of the Poisson distribution. [Conway and Maxwell, \(1962\)](#page-6-8) firstly develop the COM-Poisson for analyzing queuing systems. Later, its statistical properties are studied and reported by [Shmueli et al., \(2005\)](#page-6-9). The probability mass function (pmf) of the COM-Poisson distribution is as follows: $$Pr(Y_i = y_i \mid \gamma_i, \theta) = f(y_i; \gamma_i, \theta) = \frac{\gamma_i^{y_i}}{Q(\gamma_i, \theta)(y_i!)^{\theta}};$$ $$y_i = 0, 1, \dots : i = 1, 2, \dots n$$ (1) where *yi n* is the *i* th discrete response value, *i* is the mean value of the observations ( >0) *i* , and is the shape parameter. When < 1 , the problem is called over-dispersion and under-dispersion when > 1 . The function *Q*( , ) *i* is defined as the normalization constant and found as follow [\(Shmueli et al., 2005; Sellers and Shmueli, 2010](#page-6-9)) $$Q(\gamma, \theta) = \sum_{j=0}^{\infty} \frac{\gamma^{j}}{(j!)^{\theta}}$$ (2) and, $$Pr(Y_i = y_i \mid Y_i, \ \theta) = \frac{Y_i^{y_i}}{\sum_{j=0}^{m} \frac{y_j^{j}}{(j)^{\beta}}(y_i!)^{\theta}}; \quad y_i = 0, \ 1, \ \cdots; i = 1, \ 2, \ \cdots n$$ (3) The sum of COM-Poisson independent variables may reduce to three different probability distributions. These are (1) Poisson ( = 1) ( , *Q* ) = *e i i*, (2) Geometric ( *i*< = 1, 0) ( , *Q i* ) 1 (1 ) = *i* , and (3) Bernoulli ( [\(Shmueli et al., 2005; Sellers and Shmueli,](#page-6-9) [2010\)](#page-6-9). The moments of the COM-Poisson distribution are estimated as derivations of *Q*( , ) *i* ; these are: $$E(Y_i) = \frac{\partial logQ(\gamma_i, \vartheta)}{\partial log\gamma_i} \approx \gamma_i^{\frac{1}{9}} - \frac{\vartheta - 1}{2\vartheta}; \tag{4}$$ $$V(Y_i) = \frac{\partial E(Y_i)}{\partial \log \gamma_i} = \frac{\partial \log^2 Q(\gamma_i, \, \theta)}{\partial \log^2 \gamma_i};$$ (5) The approximations in Eqs. [\(4\)](#page-1-2) and [\(5\)](#page-1-3) are accurate when 1 or > 10 *i* ([Minka et al., 2003\)](#page-6-10). ## **3. Methodology: the penalized COM-Poisson GLM** This section introduces the COM-Poisson GLM reported by [Sellers &](#page-6-5) [Shmueli, \(2010\)](#page-6-5) and shows the steps followed in this work to develop its penalized version. The GLM using ( ( )) ( ) *E Y log* = is denoted as follows: $$logE(Y_i) = log(\gamma) = X_i \beta = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \dots + \beta_p X_{ip},$$ (6) where *Xij n* is the *i th* observation of the *j th*model predictor. The values are the coefficients of the model predictor in the GLM. There are two well-known techniques for estimating the coefficients of the GLM shown in Eq. [\(6\).](#page-1-4) These are the weighted least squares and the maximum likelihood. As the COM-Poisson is usually classified as one of the exponential family distributions, the maximum likelihood method is more preferable as a way to achieve an accurate estimation. The loglikelihood function for a set of *n* independent readings with constant is as follows ([Sellers and Shmueli, 2010\)](#page-6-5): $$\log \mathcal{L}(\gamma_{h} \vartheta | y_{1}, y_{2}, \dots, y_{n}) = \Lambda_{\gamma_{l}} = \sum_{i=1}^{n} y_{i} \mathbf{X}_{i} \boldsymbol{\beta} - \sum_{i=1}^{n} \vartheta log y_{i}! - \sum_{i=1}^{n} log Q(\mathbf{X}_{i} \boldsymbol{\beta}, \vartheta).$$ $$\tag{7}$$ The *L*2 penalized regression shrinks the coefficient estimates of the COM-Poisson GLM towards zeros, but will never become exactly zeros. The purpose of this shrinkage is mainly to prevent overfitting arising due to the collinearity between the model predictors. The penalized scheme of the maximum likelihood is shown in Eq. [\(8\)](#page-1-5) and reformatted in Eq. [\(9\):](#page-1-6) $$\lambda_{\lambda}^{Ridge} = \sum_{i=1}^{n} y_{i} \mathbf{X}_{i} \dot{\boldsymbol{\beta}} - \sum_{i=1}^{n} \theta \log y_{i}! - \sum_{i=1}^{n} \log Q(\mathbf{X}_{i} \dot{\boldsymbol{\beta}}, \theta) - \lambda \sum_{j=1}^{p} \beta_{j}^{2}.$$ (8) Or; $$\Lambda_{\lambda}^{Ridge} = \sum_{i=1}^{n} y_{i} X_{i} \beta - \sum_{i=1}^{n} \theta \log y_{i}! - \sum_{i=1}^{n} \log Q(X_{i} \beta, \theta) - \lambda ||\beta||_{2}^{2}.$$ (9) Here, 0 is a pre-chosen penalization (shrinkage) parameter that can be chosen to achieve certain penalty strength and . 2 2 is the squared Euclidean norm. Note that, larger means more shrinkage – and vice versa. The *L*2 penalization performs particularly well when the true values of are small or even zero. Let ' and ' be the solution of the penalized log-maximum likelihood in Eq. [\(9\),](#page-1-6) then the following two properties are true: - Property 1: As 0, ' , where is the solution of the COM-Poisson GLM model by [Sellers and Shmueli \(2010\).](#page-6-5) - Property 2: As , 0 ' Our regression model and the classical COM-Poisson GLM regression by [Sellers & Shmueli \(2010\)](#page-6-5) are similar in that both use the same link function, but they are dissimilar in that they use different maximum likelihood formulations. The classical COM-Poisson GLM regression using a constant dispersion parameter is a special case of our regression model when no shrinkage on the coefficients of the GLM estimates is applied, i.e. =0. In the next section, we refer to our model as the penalized COM-Poisson GLM. ## **4. Penalized likelihood estimation based on the cross-validation** The *L*2 penalization may cause some bias. The proper selection of the penalization parameter ( ) reduces the negative effect of large bias and results in a better mean squared error. [Hoerl and Kennard, \(1970\)](#page-6-11) firstly proposed this parameter when the ridge regression is introduced as a solution for the collinearity problem. Since then, several methods have been developed to find the most efficient estimate of this parameter ([Kibria, 2003; Khalaf and Shukur, 2005; Alkhamisi et al., 2006;](#page-6-12) [Alkhamisi and Shukur, 2008; Muniz and Kibria 2009\)](#page-6-12). The ridge trace plotted by [Hoerl and Kennard, \(1970\)](#page-6-11) is a graphical method to determine an appropriate value of . This method suggests visualizing the changes in the GLM coefficients over a wide range of values and to select the value of when the changes in these coefficients are insignificant. This method is a kind of trade-off between the variance and the bias of the coefficient estimates. To maintain a small bias, the minimum value of the corresponding range of would be an appropriate choice. This method is applied in this paper to specify the practical range of . Later, K values of are selected from this range, including extremes, and arbitrarily assigned to the K cross-validation cycles. The proposed approach for solving the penalized likelihood estimation under the cross-validation is described below: - **Step 1**: Split the original dataset into K equally sized sub-datasets. - **Step 2**: Define *w n* = K as the sub-dataset size, where *n* is the size of the original dataset. - **Step 3:** Let // the subset containing indices of K sub-datasets. - **Step 4**: Let *j* // the order index of a sub-dataset, where *j* = 1, 2, , K - **Step 5:** Let be a subset of size K containing all possible values of . That is, = { 1, 2, , *K*}. - **Step 6**: For *j* = 1 to K - (i) Name the *j* th sub-dataset in as a "training" sub-dataset and the (K 1) remaining sub-datasets as the "validating" sub-datasets. - (ii) Solve Equation [\(9\)](#page-1-6) and find ' *j* and ' . Note that, is the tuning parameter applied when the *j th* sub-dataset is selected as the "training" sub-dataset. - (iii) Calculate the Mean Absolute Error (MAE), which is one of the methods available in the literature to measure the prediction accuracy. The Mean Absolute Error for any sub-dataset in is calculated as follows: $$MAE_j = \frac{1}{w} \sum_{i=1}^{w} |y_i - \hat{y}_i| = \frac{1}{w} \sum_{i=1}^{w} |e_i|,$$ (10) where the fitted value *yi* can be calculated by applying Eq. [\(4\)](#page-1-2) after setting $$\hat{y}_i = exp\left(x_i^T \hat{\beta}_{\lambda_j}\right)$$ . The error, $e_i = y_i - \hat{y}_i$ , is calculated as follows: $e_i = y_i - \hat{\beta}_0 - \hat{\beta}_1 X_{i1} - \hat{\beta}_2 X_{i2} - \cdots - \hat{\beta}_p X_{ip}$ (11) (iv) Compare the Median Absolute Deviation (MAD) for each cycle as follows: $$MAD_j = \text{median } |MAE_t - \text{median}(\mathbf{M}_j)|, t = 1, 2, \dots, (K - 1),$$ (12) where **M***j* is a list of MAE values defined as $$\mathbf{M}_{j} = [\mathrm{MAE}_{1}, \, \mathrm{MAE}_{2}, \cdots, \mathrm{MAE}_{j}] \tag{13}$$ - **Step 7**: End of the "For" loop - **Step 8**: Let **D** the matrix contains the MAD*j* values of all validation cycles. - **Step 9**: Find MADmin to specify the optimal penalization parameter as follows: $$MAD_{min} = argmin\{\mathbf{D}\} = arg \min_{j=1,...,K} \{MAD_j\}$$ (14) - **Step 10**: Let *opt* be the value of used in the cross-validation cycle having the minimum MAD. This cycle is already defined in **Step 9**. - **Step 11**: By using *opt* , fit the penalized COM-Poisson GLM model to the whole dataset and find *opt* and *opt* . - **Step 12**: Let and *opt* be the optimal solution of the penalized maximum likelihood (PML) problem Eq. [\(9\)](#page-1-6) could be directly solved by using the constrained nonlinear optimization tool (nlminb) available in R software. In this paper, the R Software available at [https://cran.r-project.org/web/packages/](https://cran.r-project.org/web/packages/COMPoissonReg) [COMPoissonReg](https://cran.r-project.org/web/packages/COMPoissonReg) was modified to capture the structure of our penalized likelihood estimation. Further methods for solving such problem are available in [Sellers and Shmueli, \(2010\).](#page-6-5) Another way to find y*i* is by obtaining the inverse-CDF of y*i* and [\(Sellers and Shmueli, 2010\)](#page-6-5). In order to specify the optimal penalization parameters, we calculate the MAD for each cycle using the MAE values, which are used to measure the prediction accuracy, and select the parameters for the minimum value of MAD, that is, MADmin, as the optimal penalization parameters. This mechanism for the optimal parameters is described in the steps 6, 9 and 10. The penalized likelihood estimation is widely used to reduce overfitting in estimating the regression parameters, including linear, logistic, semi-parametric and proportional hazards regression [\(Verweij](#page-6-13) [and Van Houwelingen, 1994; Ciuperca et al., 2003; Rondeau et al.,](#page-6-13) [2003; Goeman, 2010; Göbl et al., 2015\)](#page-6-13). The overfitting problem may arise due to the collinearity among predictors or in the condition of the high dimensionality. [Goeman \(2010\)](#page-6-14) provided several GLMs in R software, named as "*penalized*" package, with *L*1 – LASSO and/or *L*2 – ridge penalty terms, using likelihood-cross validations. The regression models supported by the "*penalized*" package are the linear, logistic, Poisson regression, and the Cox Proportional Hazards model. In this paper, we are specifically interested in tackling the over- and under-dispersion issues arising during the crash count regression for **Table 1** Summary statistics for the *K* sub-datasets. | Statistic | Sub-dataset | | | | | | | | | |-----------|-------------|-------|-------|-------|--|--|--|--|--| | | 1 | 2 | 3 | 4 | | | | | | | Min | 0 | 0 | 0 | 0 | | | | | | | Max | 54 | 51 | 53 | 48 | | | | | | | Mean | 13.76 | 11.25 | 10.95 | 10.25 | | | | | | | Median | 11 | 7 | 7 | 8 | | | | | | | AIC | 1401 | 1335 | 1503 | 1332 | | | | | | | R2 | 0.30 | 0.21 | 0.20 | 0.17 | | | | | | **Table 2** The VIF and Tolerance values of X1 and X2 for the 4-tested sub-datasets. | | Regressors | | | | | | | | | | | | |----------------------|--------------|--------------|--------------|--------------|--------------|--------------|--------------|--------------|--|--|--|--| | | X1 | | | | X2 | | | | | | | | | | Sample 1 | Sample 2 | Sample 3 | Sample 4 | Sample 1 | Sample 2 | Sample 3 | Sample 4 | | | | | | VIF
Tolerance (T) | 1.15
0.86 | 1.14
0.88 | 1.11
0.90 | 1.08
0.92 | 1.15
0.86 | 1.14
0.88 | 1.11
0.90 | 1.08
0.92 | | | | | **Table 3** Comparison of estimated coefficients in the first cycle. | Parameter | Regression Model | | | | | | | | | | |-----------|------------------|-----------------|---------------------------|--|--|--|--|--|--|--| | | Poisson | COM-Poisson GLM | Penalized COM-Poisson GLM | | | | | | | | | | – | 0.39360 | 0.37940 | | | | | | | | | 0 | −9.9467 | −10.988 | −10.761 | | | | | | | | | 1 | 0.5515 | 0.5819 | 0.5561 | | | | | | | | | 2 | 0.7361 | 0.8051 | 0.8086 | | | | | | | | which we use the ridge COM-Poisson GLM. Note that the ridge Poisson GLM is a special case of our ridge COM-Poisson GLM when *Q*( , ) = *e i i*. Even under this special case of similarity, the estimated regression parameters resulted by using the ridge COM-Poisson GLM are expected to be different from [Goeman, \(2010\),](#page-6-14) because the outcomes of the cross-validation vary according to the run-up settings and the criteria of a decision. ### **5. An illustrative example: the MVC-F in Toronto, Ontario** In this section, we compare the prediction accuracy of the penalized COM-Poisson GLM with both the Poisson and the classical COM-Poisson GLM. For a fair comparison, this paper uses the overdispersed crash dataset collected in 1995 at 868 signalized intersections located in Toronto, Ontario. Several researchers have used this dataset [\(Lord,](#page-6-15) [2000; Miaou and Lord, 2003; Miranda-Moreno and Fu 2007; Lord et al.,](#page-6-15) [2008; Sellers and Shmueli, 2010](#page-6-15)). For a further detailed description of this dataset, see [Lord et al. \(2008\)](#page-6-4). There are 868 observations in the dataset which consist of the MVC-F, as a response variable ( Y ), and two log-transformed traffic flow, as predictors ( X1 and X2). In this example, the original dataset is divided into four sub-datasets for applying the systematic K -fold CV approach. Each of them contains 217 observations. According to this division, the parameters K and *w* are set at 4 and 217, respectively. [Table 1](#page-2-1) summarizes the statistics of the sub-datasets of one cycle solution. Many diagnostics have been proposed in the previous literature to assess the degree of collinearity[\(Abdella et al. 2019\)](#page-5-1). The Pearson correlation coefficient ( ) **R**2 is the typical among them. The **R**2 takes a value between −1 and +1, where a negative value indicates a negative correlation, 0 means no correlation between the two variables and a positive value means positive correlation. The **R**2 in the Toronto crash dataset is 0.31. For each sub-dataset, the **R**2 of the two traffic flow regressors is calculated and reported in [Table 1](#page-2-1). The collinearity, can also be indicated by using the Variance Inflation Factor (VIF) of regression coefficients, *i* . The VIF shows the degree of influence of the multicollinearity on the variance of regression estimates, which is defined as $$VIF_{i} = \frac{1}{(1 - r_{i}^{2})}; \qquad i = 1, 2, \dots, p,$$ (15) where *p* is the number of model variables and *ri* 2 is the **R**2 value when the regressor X*i* is regressed on the other factors. An equivalent indicator to the VIF is the Tolerance (T = VIF−1), which measures the amount of information whether the regressor X*i* can deliver under a particular level of collinearity. The Tolerance takes a maximum value of 1. If the regressor X*i* has no correlation (*ri* 2 = 0) with the other ( *p* 1 regressors, then the T is equal to 1. The VIF and T values of all subdatasets for the classical Poisson regression are calculated in [Table 2](#page-3-1). When VIF >10 , the model predictors are considered to be highly correlated with each other. The multicollinearity level in the Toronto crash dataset is relatively low; see [Table 2.](#page-3-1) The effective range of λ of the Toronto crash dataset is found to be from 0.1 to 0.4. Four values of , including extremes, are selected. These are = { 1, 2 3 4 , , } = { 0.1, 0.2, 0.3, 0.4 }. [Table 3](#page-3-2) displays the results of fitting the penalized COM-Poisson GLM to the first sample of the crash data. Two more regression models are applied in this cycle and compared with the proposed model. These are the Poisson model and the classical COM-Poisson GLM by [Sellers and Shmueli, \(2010\)](#page-6-5). Note that, the dispersion parameters ( of both the penalized and classical COM-Poisson GLM regression models are assumed constant and unknown across observations. [Table 3](#page-3-2) indicates that the first "training" sub-dataset is characterized with over-dispersion ( = < 0. 3794 1) . This finding is used here to justify the usage of the COM-Poisson instead of the ordinary Poisson model. In order to create a scale consistency of the estimated coefficients of the three compared models, we follow [Sellers](#page-6-5) [and Shmueli, \(2010\)](#page-6-5) and report the estimated coefficients of the penalized and classical COM-Poisson models in terms of . The comparison in [Table 3](#page-3-2) is only to show the extent to which the compared models are close to each other. When developing prediction models, it is essential to test how well the data-set meets the assumptions of the model. To check the assumption that the residuals are uncorrelated with each other, we ![](_page_3_Figure_15.jpeg) ![](_page_3_Figure_16.jpeg) ![](_page_3_Figure_17.jpeg) **Fig. 1.** Plotting the deviance residuals produced by the fitted penalized COM-Poisson GLM at *j*=1 versus the observations order. ![](_page_4_Figure_0.jpeg) ![](_page_4_Figure_1.jpeg) **Fig. 2.** The Histogram plot for the deviance residuals of the three "validating" sub-datasets of the first cycle. ![](_page_4_Figure_3.jpeg) **Fig. 3.** Plotting the Mean Absolute Error (MAE) under different levels of shrinkage in the coefficients of the penalized COM-Poisson GLM regression. **Table 4** Calculating the MAD of the first cycle of validation. | Sample | MAE | M )1
Median( | MAE
median(M ) 1 | MAD1 | |--------|------|-----------------|----------------------|--------| | 2 | 3.96 | 3.96 | 0.0000 | 0.1887 | | 3 | 4.74 | | 0.7843 | | | 4 | 3.77 | | 0.1887 | | **Table 5** Summary of the last three cycles of cross-validation. | Cycle,
j | j | | Penalized COM-Poisson GLM Coefficients | MADj | | | |-------------|-----|--------|----------------------------------------|--------|--------|--------| | | | | 0 | 1 | 2 | | | 2 | 0.2 | 0.3590 | −12.502 | 0.7269 | 0.7970 | 0.1866 | | 3 | 0.3 | 0.2614 | −11.902 | 0.6501 | 0.8099 | 0.2136 | | 4 | 0.4 | 0.3825 | −10.313 | 0.5886 | 0.7146 | 0.0289 | **Table 6** Estimated GLM coefficients of the Toronto, Ontario crash dataset. | Parameter | Regression Model | | | | | | | |-----------|------------------|-----------------|---------------------------|--|--|--|--| | | Poisson | COM-Poisson GLM | Penalized COM-Poisson GLM | | | | | | | – | 0.34917 | 0.33714 | | | | | | 0 | −10.234 | −11.702 | −11.507 | | | | | | 1 | 0.6029 | 0.6558 | 0.6333 | | | | | | 2 | 0.7037 | 0.7911 | 0.7945 | | | | | | MSE | 23.60 | 32.57 | 32.54 | | | | | calculated the residual values of the three "validating" sub-datasets and plotted these values versus the order of the observations ([Fig. 1](#page-3-3)). The residual plots are supposed to show random patterns. Non-random patterns may violate the assumption that the residuals are independent of each other. Fortunately, here, the independence assumption of residuals ![](_page_4_Figure_12.jpeg) **Fig. 4.** The scatter plot of MSE of fitting the penalized COM-Poisson GLM, classical COM-Poisson GLM and Poisson models for the crash dataset of Toronto. is confirmed as the three "validating" sub-dataset reveal a random behaviour of residuals. Numerous techniques are available to test the normality assumption of residuals, i.e. *e N µ i* ( , )2 . The simplest one is the Histogram, which is also an efficient to present whether the residuals are skewed or whether the dataset includes outliers. [Fig. 2](#page-4-0) shows an approximately normal distribution of the residuals produced by the three "validating" sub-datasets. To build the **M**1 matrix**,** the first cycle solution and 0.1 is used to compute the MAE of each "validating" sub-dataset. The **M**1 matrix is found to equal [3.96 4.74 3.77] *T***.** In our cross-validation approach, a K -number of values are selected and arbitrarily assigned to the cycles of cross-validation. The best results of our approach should be reached when the value MAE from each "validating" sub-datasets at a cycle *j* is as low as possible which can be achieved over the range of . To investigate the convenience of choosing 1 = 0.1 for the first cycle, we plot the MAE values of the three "validating" sub-datasets versus several values of ; see [Fig. 3](#page-4-1). Graphically, for all "validating" sub-datasets, setting 1=0.1 is the best choice among others. Note that, the importance of performing the "suitability check" step decreases as the slope of the practical range of decreases – and vice versa. The MAD1 value is calculated as shown in **Table 7** Summary statistics for the *K* = {2,6,8} sub-datasets. | Statistic | K = 2 | | K = 6 | | | | | | K = 8 | | | | | | | | |-----------|-------|------|-------|------|------|------|------|------|-------|------|------|------|-----|------|------|------| | | 1 | 2 | 1 | 2 | 3 | 4 | 5 | 6 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | | Min | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | | Max | 54 | 53 | 51 | 53 | 47 | 46 | 54 | 50 | 51 | 53 | 47 | 46 | 51 | 50 | 54 | 50 | | Mean | 10.6 | 12.4 | 10.4 | 13.9 | 10.4 | 11.3 | 10.9 | 12.0 | 10.5 | 14.7 | 10.6 | 12.1 | 9.4 | 10.9 | 11.8 | 12.7 | | Median | 8 | 18 | 7 | 10.5 | 8 | 7 | 9 | 10 | 7 | 11 | 8 | 8 | 7.5 | 9 | 10 | 9 | **Table 8** The coefficients of estimation using K = {2,6,8} | Parameter | | Penalized COM-Poisson GLM | | | | | | | |-----------|---------|---------------------------|---------|--|--|--|--|--| | | K = 2 | K = 6 | K = 8 | | | | | | | 0 | −3.6828 | −4.737 | −4.9727 | | | | | | | 1 | 0.2034 | 0.3168 | 0.3377 | | | | | | | 2 | 0.2533 | 0.248 | 0.2429 | | | | | | | MSE | 32.95 | 31.57 | 29.49 | | | | | | [Table 4.](#page-4-2) Having the MAD1 value is calculated, we repeat the above procedures with *j* = 2, 3 and 4. [Table 5](#page-4-3) summarizes the results of the last three cycles of the cross-validation approach. Comparing the MAD*j* values in [Table 5,](#page-4-3) the fourth cycle in which 4=0.4 has the lower MAD. Hence, we set *opt* = 0.4 as the optimal shrinkage parameter and fit the penalized COM-Poisson GLM model to the whole crash dataset; see [Table 6](#page-4-4). The estimated GLM coefficients of both COM-Poisson GLM model are very close ([Table 6\)](#page-4-4). The results of fitting the classical COM-Poisson GLM are identical to those reported by [Sellers and Shmueli \(2010\)](#page-6-5). How well the fitted model predicts the response variable is necessary to assess the model adequacy. In this example, we use the mean squared error (MSE) as a measure of adequacy. Since the well-fitted model leads to the minimum MSE value, the penalized GLM COM-Poisson is the best to fit the Toronto crash dataset (MSE = 32.54) shown in [Table 6](#page-4-4). The last column in [Table 6](#page-4-4) represents the optimal solution of Eq. [\(9\).](#page-1-6) For further understanding of the performance of the penalized COM-Poisson GLM regression model, one can compare its performance with the Negative Binomial (NB) and the restricted generalized Poisson regression (RGPR) models available in [Sellers and Shmueli \(2010\).](#page-6-5) To draw inference about the sensitivity of the penalized COM-Poisson GLM to changes in , we proposed different levels of *λ* and computed the MSE at each level; see [Fig. 4.](#page-4-5) The MSE values of both the Poisson and the classical COM-Poisson GLM models were plotted for the purpose of comparison. Since the last two models are not penalizationbased models, their MSE values over the proposed range of *λ* are set fixed at their magnitudes reported in [Table 6.](#page-4-4) These are 32.60 and 32.57, respectively. [Fig. 4](#page-4-5) shows that the penalized COM-Poisson GLM provides the lowest MSE value when it is compared with other models. The one exception here is the case when *λ* = 1.2. In this case, the classical COM Poisson GLM outperforms our model. Moreover, as the range of *λ*, in this example, is set from 0.1 to 0.4, the penalized COM Poisson GLM would deliver the lowest MSE regardless of how the values are assigned to the cycles of the cross-validation. ## **6. Optimizing the K-parameter** The selection of K is critical to the quality of the cross-validation. In [Section 5](#page-3-0), we used K = 4. The results have shown the advantage of the penalized-COM Poisson over its counterparts the classical COM-Poisson and the Poisson regression. In this section, we conduct an analytical study to investigate the performance of the Penalized-COM-Poisson under several values of K. Three different values of K = {2,6,8} were suggested. [Table 7](#page-5-2) shows the descriptive analysis of the sub-dataset under each value of K. [Table 8](#page-5-3) reports the fitting results using the three values of K. The comparison between the MSE values of the three models shows a significant advantage of the third model using K = 8. The case of K = 2 has the worst MSE values comparing with all K values considered in this study. The three proposed models, similar to the case of K = 4, have a negative intercept. ## **7. Conclusion** This paper integrates the penalized likelihood estimation, more specifically the ridge penalty function, with the COM-Poisson GLM regression to enhance its prediction accuracy under the condition of the collinearity. Using the penalized likelihood estimation makes the COM-Poisson GLM regression more rigorous to inflation in the standard errors of GLM estimates. The real-world example illustrated the excellent performance of the penalized COM-Poisson GLM regression in simultaneously accommodating both the dispersion and collinearity problems. To specify the value of the penalization parameter , we applied the the penalized likelihood estimation under the K-*fold* cross-validation. A limitation associated with this approach is the case of a small dataset in which the estimate of the penalization parameter *opt* is unreliable. To handle this situation, we recommend to select K = 1 or K n = . The second option ( K n = ) is known in literatures as leave-one-out (LOO) cross-validation. The preference of the penalized COM-Poisson GLM should significantly arise in high-dimensional count data regression in which the collinearity is more likely to occur. Comparing the COM-Poisson GLM with those in the "*penalized*" package, in modeling high-dimensional crash data, could be one of our future works. Under the high-dimensional road safety, collecting large-scale observations is an expensive task in terms of time and effort. According to [Sellers and Shmueli, \(2010\),](#page-6-5) the COM-Poisson has the flexibility to reduce to a logistics regression and to produce similar estimates. Motivated by the importance of the logistic regression in several real-world disciplines, this finding would also be interesting for future investigations. #### **Acknowledgment** Authors are grateful to Dominique Lord (Texas A&M University) and Srinivas Geedipally (Texas A&M Transportation Institute) for providing us the dataset of the real-example of Toronto, Ontario. #### **References** [Abdella, G.M., Alhajyaseen, W., Al-Khalifa, N.K., Hamouda, S.A., 2016. Usage of non](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0005)[linear regression for modeling the behavior of motor vehicle crash fatality \(MVF\)](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0005) [rate. In: The Proceedings of the International Conference on Industrial Engineering](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0005) [and Operations Management, Rabat, Morocco, April 11–13, 2017, pp. 1827–1834](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0005). [Abdella, G.M., Al-Khalifa, N.K., Maha, A., Tayseer, Hamouda, S.A., 2019. Modelling](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0010) [trends in road crash frequency in Qatar State. Int. J. Operat. Res. 34 \(4\), 507–523](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0010). [Abdur Rouf, K.B., Abdella, G.M., Al-Khalifa, N.K., Alhajyaseen, W., 2018. Ridge](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0015) - [penalization-based generalized linear model \(GzLM\) for predicting risky-driving](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0015) [index. In: The Proceedings of the International Conference on Industrial Engineering](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0015) [and Operations Management, Washington DC, USA, September 27–29, 2018, pp.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0015) [1462–1473](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0015). - [Alkhamisi, M.A., Shukur, G., 2008. Developing ridge parameters for sur model. Commun.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0020) [Stat. Theory Meth. 37, 544–564.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0020) - [Alkhamisi, M.A., Khalaf, G., Shukur, G., 2006. Some modifications for choosing ridge](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0025) [parameter. Commun. Stat. Theory Meth. 35, 1–16](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0025). - [Ciuperca, G., Ridolfi, A., Idier, J., 2003. Penalized maximum likelihood estimator for](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0030) [normal mixtures. Board Found. Scand. J. Stat. 30, 45–59.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0030) - [Conway, R.W., Maxwell, W.L., 1962. A queuing model with state-dependent service rates.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0035) [J. Ind. Eng. 12, 132–136](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0035). - [Frank, G., Samia, I., 2015. Determinants of seat belt use: a regression analysis with FARS](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0040) [data corrected for self-selection. J. Saf. Res. 55, 7–15](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0040). - [Ghadban, Nazha R., Abdella, Galal M., Alhajyaseen, Wael, Alkhalifa, Khlaifa N., 2018.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0125) [Analyzing the impact of human characteristics on the comprehensibility of road](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0125) [traffic signs. The Proceedings of the International Conference on Industrial](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0125) [Engineering and Operations Management, Bandung, Indonesia, March 6–8, 2018](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0125). - [Göbl, C.S., Bozkurt, L., Tura, A., Pacini, G., Kautzky-Willer, A., Mittlböck, M., 2015.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0045) [Application of penalized regression techniques in modelling insulin sensitivity by](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0045) [correlated metabolic parameters. PLoS ONE 10 \(11\).](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0045) - Goeman, J.J., 2010. L1 [penalized estimation in the cox proportional hazards model.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0050) [Biometr. J. 52, 70–84.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0050) - [Guikema, S.D., Coffelt, J.P., 2007. A flexible count data regression model for risk analysis.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0055) [Risk Anal. 28, 213–223](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0055). - [Hoerl, A.E., Kennard, R.W., 1970. Ridge regression: biased estimation for non- orthogonal](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0060) [problems. Technometrics 12, 55–67.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0060) - [Khalaf, G., Shukur, G., 2005. Choosing ridge parameters for regression problems.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0065) [Commun. Stat. Theory Meth. 34, 1177–1182](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0065). - [Kibria, B.M., 2003. Performance of some new ridge regression estimators. Commun. Stat.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0070) [Simul. Comput. 32, 419–435.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0070) - [Kokonendji, C.C., 2014. Over- and Underdispersion Models, Methods and Applications of](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0075) [Statistics in Clinical Trials: Planning, Analysis, and Inferential Methods. John Wiley](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0075) [and Sons Inc.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0075) - [Lord, D., 2000. The Prediction of Accidents on Digital Networks: Characteristics and](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0080) [Issues Related to the Application of Accident Prediction Models. Ph.D. Dissertation.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0080) [Department of Civil Engineering, University of Toronto, Toronto, Ontario](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0080). - [Lord, D., Mannering, F., 2010. The statistical analysis of crash-frequency data: a review](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0085) - [and assessment of methodological alternatives. Transp. Res. Part A 44, 291–305](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0085). - [Lord, D., Persaud, B.N., 2000. Accident prediction models with and without trend: ap](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0090)[plication of the generalized estimating equation procedure. Transp. Res. Rec. 1717,](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0090) [102–108](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0090). - [Lord, D., Guikema, S.D., Geedipally, S.R., 2008. Application of the](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0095) [Conway–Maxwell–Poisson generalized linear model for analyzing motor vehicle](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0095) [crashes. Accid. Anal. Prev. 40, 1123–1134](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0095). - [McCullagh, P., Nelder, J.A., 1997. Generalized Linear Models, second ed. Chapman and](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0100) [Hall/CRC, London, pp. R0727836.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0100) - [Miaou, S.P., Lord, D., 2003. Modelling traffic crash-flow relationships for intersections:](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0105) [dispersion parameter, functional form, functional form, and Bayes versus empirical](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0105) [Bayes methods. Transp. Res. Rec. J. Transp. Res. Board 1840, 31–40](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0105). - [Minka, T.P., Shmueli, G., Kadane, J.B., Borle, S., Boatwright, P., 2003. Computing with](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0110) [the COM-Poisson Distribution. Technical Report 776. Dept. Statistics, Carnegie](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0110) [Mellon Univ., Pittsburgh, PA](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0110). - [Miranda-Moreno, L.F., Fu, L., 2007. Traffic safety study: empirical Bayes or Full Bayes?](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0115) [The 84th Annual Meeting of the Transportation Research Board, Washington, D.C.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0115) [Paper 07-1680](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0115). - [Muniz, G., Kibria, B.M.G., 2009. On some ridge regression estimators: an empirical](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0120) [comparisons. Commun. Stat. Simul. Comput. 38, 621–630](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0120). - [Paraskevi, M., Quddus, M.A., David, P., Andrew, H., 2015. Exploring the factors affecting](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0130) [motorway accident severity in England using the generalized ordered logistic re](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0130)[gression model. J. Saf. Res. 55, 89–97](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0130). - [Rondeau, V., Commenges, D., Joly, P., 2003. Maximum penalized likelihood estimation in](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0135) [a gamma-frailty model. Lifetime Data Anal. 9, 139–153](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0135). - [Sellers, K.F., Andrew, R., 2016. A flexible zero-inflated model to address data dispersion.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0140) [Comput. Stat. Data Anal. 99, 68–80](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0140). - [Sellers, K.F., Shmueli, G., 2010. A flexible regression model for count data. Ann. Appl.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0145) [Stat. 4, 943–961.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0145) - [Shmueli, G., Minka, T.P., Kadane, J.B., Borle, S., Boatwright, P., 2005. A useful dis](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0150)[tribution for fitting discrete data: revival of the Conway-Maxwell-Poisson distribu](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0150)[tion. J. Roy. Stat. Soc.: Ser. C \(Appl. Stat.\) 54, 127–142](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0150). - [Vangala, P., Lord, D., Geedipally, S.R., 2015. Exploring the application of the Negative](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0155) [Binomial – Generalized Exponential model for analyzing traffic crash data with ex](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0155)[cess zeros. Anal. Meth. Accid. Res. 7, 29–36](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0155). - [Verweij, P.J.M., Van Houwelingen, J.C., 1994. Penalized likelihood in cox regression.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0160) [Stat. Med. 13, 2427–2436.](http://refhub.elsevier.com/S0925-7535(18)31395-X/h0160)