Categorical Data Analysis Using a Skewed Weibull Regression Model

In this paper, we present a Weibull link (skewed) model for categorical response data arising from binomial as well as multinomial model. We show that, for such types of categorical data, the most commonly used models (logit, probit and complementary log–log) can be obtained as limiting cases. We further compare the proposed model with some other asymmetrical models. The Bayesian as well as frequentist estimation procedures for binomial and multinomial data responses are presented in detail. The analysis of two datasets to show the efficiency of the proposed model is performed.


Introduction
The statistical problem of estimating binary response variables is very important in many areas including social science, biology and economics [1]. The vast bibliography of categorical data presents the big evolution of the methods that handle appropriately binary and polychotomous data. More details can be found in Agresti [2]. Generalized linear model (GLM) has a wide range of tools in regression for count data [3]. Two important and commonly used symmetric link functions in GLM are the logit and probit links [4]. Many studies have investigated the limitations of these symmetric link functions. It is well accepted that when the probability of the binary response approaches 0 at a different rate from the rate (as a function of covariate) it approaches 1, symmetric link functions cannot be appropriate [5]. Many parameteric classes of link functions are in the literature, including the power transform of logit link by Aranda-Ordaz [6] and the a general link class of Chen et al. [5]. Other works with one-parameter class include Guerrero and Johnson [7], Morgan [8], Whittmore [9] and a host of others. Existing models for two-parameter families include Stukel [10], Prentice [11], Pregibon [12], Czado [13] and Czado [14].
Stukel's model with transformation of both tails of logit link is very general and can approximate many important links including probit, logit and complementary log-log. However, the Bayesian analysis of Stukel's model is not straightforward to implement, particularly in presence of multiple covariates and noninformative improper priors. The model proposed by Chen et al. [5], which includes the skew-probit model, uses a latent variable approach [15] that is convenient for sampling from the posterior distribution. Using the Albert and Chib [15] technique, Kim et al. [16] proposed the generalized t-link models, Naranjo et al. [17] proposed the asymmetric exponential power (AEP) model, and Rubio and Liseo [18] discuss the Jeffreys prior for skew-symmetric models. However the frequentist analysis for these models are not trivial. For the skew-probit model, The existence of the maximum likelihood estimator (MLE) of the linear regression parameters (β) can be proved only under the restrictive condition that the skewness parameter of the link function is known [19].

Link function
Let X = (1, X 1 , . . . , X r ) be the design matrix, where 1 is a vector with all values equal to 1, j = 1, . . . , r. We denote the vector of binary response variable as Y. Similar to GLM, our interest lies in modeling the probability Pr . . , n, where η = βX, β = (β 0 , β 1 , . . . , β r ) are the linear coefficients, and g(·) is the link function. The link function relates the covariates X with the mean response µ = E(Y | X). In this case, the g −1 is a cumulative distribution function (cdf) on the real line. Our interest is a link function that can accommodate symmetric and asymmetric tails which has a simple parameteric functional form, and can be easily tractable. To obtain these goals, we use the cdf of Weibull distribution for g −1 , where α ∈ R is the location/threshold parameter, γ > 0 is the shape parameter, and I (η>α) is the indicator of η > α. I A is the indicator function of event A, that is I A = 1 and I A c = 0. Alternatively, the Weibull link function is defined as where µ(η) = E(Y | η) ≥ 0, γ > 0, and η > 0.
Note that, in the above parameterization, the restriction of η > 0 is not a problem because the parameter β 0 plays the role of both the location/threshold parameter α and the intercept of linear predictor η = βX. By doing this, we avoid the identifiability problem in estimation of β 0 , also we have a more parsimonious model. The skewness of the Weibull link depends only on the parameter γ, and can be evaluated by ( is the Gamma function. The skewness lies in the interval (−1.1395, ∞). We also evaluated the Arnold-Groeneveld (AG) skewness measure [26], which is a skewness measure related to the mode of a distribution. Again, the AG skewness depends only on the parameter γ, and can be evaluated as 2 exp{(1 − γ)/γ} − 1, and lies in the interval (−0.26424, ∞). However, sometimes, a model with skewness lower than −1.1395 is desired; in this case, we can use the reflected Weibull distribution to define the link as µ = g −1 (ω) = exp{−η γ }, and the skewness lies in the interval (−∞, 1.1395). The different forms of Weibull link are shown in Figure 1 with solid line for the Weibull link and dashed line for the reflected Weibull link.

Special Cases
The choice of the Weibull distribution as link function is due to its flexible properties. Rinne [27] (Chapter 3) discusses the various properties of Weibull along with Weibull distribution as approximation to some symmetrical distributions. We highlight the relations of Weibull with the normal and logistic distributions, because they explain the relations of Weibull link with probit and logit link functions. Based on results of Rinne [27], we have where Φ is the distribution function of the standard normal distribution. These results show that Weibull link can approximate the probit link and the logit link. The degrees of these approximations are illustrated in Figure 2. We have the following proposition for another important case of link, the complementary log-log link [4]. Proposition 1. The complementary log-log link defined by g −1 (η) = 1 − exp{− exp(η)} is a limiting case of the Weibull link because Proof. Taking α = −1 in Equation (1) and dividing η by γ, without loss of generality, we can rewrite the Weibull link given in Equation (2) as: Now, taking the limit γ → ∞ of g −1 (η) completes the proof.
Given this result, we can say that for a dataset when the estimated value of γ is large then the complementary log-log link should be appropriate. Using the reflected Weibull link, we have a similar result with the log-log link, defined as g −1 (η) = exp{− exp(−η)} [4]. The complementary log-log and log-log link as limiting cases are illustrated in Figure 3.

Binomial Data
Consider a sample of size n from the binary variable/response Y, with Pr[Y i = 1] = p i for i = 1, . . . , n. We denote the observed data as D = {n, Y = y, X = x}, where y = (y 1 , . . . , y n ) is the observed vector of Y = (Y 1 , . . . , Y n ), and x = (1, x 1 , . . . , x r ) , is the observed covariate matrix of X = (1, X 1 , . . . , X r ) . The likelihood function for the Weibull link can be written as and the log-likelihood as where η i is the i-th element of the vector η = βX, and β, γ are the parameters to be estimated.
A numerical method such as Nelder and Mead [28] can be used to obtain the MLE for (β, γ). The expression of the gradient vector and Hessian matrix are given in Appendix A. Using the gradient vector and Hessian matrix, it is simple to implement a Newton-Raphson algorithm to obtain the MLE. As initial guesses for the numerical algorithm, we suggest to use the estimator β i,probit under probit model for β i (i = 0), β 0,guess = − min( β probit x) + 0.001 for β 0 , and 3.60235 for γ. The initial guesses (β, γ) can be interpreted as the Weibull link being an approximate probit link.
For the Bayesian analysis, the posterior density is where p(β, γ) is the joint prior. We suggest using the hierarchical Bayes model. Assuming the parameters are a priori independent, the first level of hierarchy has γ following a gamma distribution with mean m γ and variance v γ , and β with multivariate normal distribution with mean vector m β and covariance matrix v β I, where I is the identity matrix. The values of v γ and v β are fixed, and for m γ and m β we consider a prior, that is p(m γ ) = p(m β ) ∝ 1. Arguably, we can use the mode of the integrated likelihood of (m γ , m β ) to determine a prior distribution [23]. The hyper-parameters v γ and v β are viewed as prior precision parameters. The EM (Expectation-Maximization) algorithm [29] can be used to obtain the estimates of m γ and m β . The MCMC procedure is used to generate a sample from the posterior distribution. For the MCMC procedure, we used a Gibbs sampler with Metropolis-Hasting. The convergence of the chain was monitored using ergodic means. We omit the details about these computational tools because they are already well known tools and are not the main subject of the present paper. In addition, it was not necessary to develop any special scheme to sample from the posterior chain. Another advantage of the Weibull link is that the posterior distributions are proper even when we use a wide range of non-informative priors. The Jeffreys' prior for the parameter β has the form p(β | γ) ∝ |I(β | γ)| 1/2 , where the Fisher information matrix I(β | γ) can be obtained by taking the expectation of the Hessian matrix given in Appendix A.
Considering the improper prior p(β) ∝ 1, and the non-informative prior p(γ) ∝ 1/γ c , for γ > 1 and c > 1 a known constant [30], we have the non-informative prior distribution With this constraint (in the parameter γ of the Weibull link), the skewness lies in the interval (−1.1395, 2], which is still a flexible link. For the improper prior of Equation (7), the propriety of the resulting posterior distribution in Equation (6) is stated in Theorem 1. Theorem 1. Let z i = −1 when y i = 0 and z i = 1 when y i = 1, and X * be the matrix with rows z i x i . Suppose that the design matrix X is of full rank, and there exists a positive vector a = (a 1 , . . . , a n ) ∈ R n , with a i > 0, for i = 1, . . . , n, such that X * a = 0, under the non-informative prior of Equation (7), then the posterior density Equation (6) is proper.
Proof. Let u, u 1 , . . . , u n be independent random variables with common Weibull distribution with shape parameter γ.
From Lemma 4.1 of Chen and Shao [31] there exists a constant K depending only on X * such that This prior give a constraint in the parameter γ. However, any proper prior can be used with the proposed model, avoiding any constraint problem in the parameter γ.

Multinomial Data
For multinomial responses, we have that Y i ∈ {1, . . . , K}, and p k = Pr(Y i = k), for k = 1, . . . , K and ∑ K j=1 p j = 1. The logistic multinomial regression model consider a reference category, generally the category K, and have a link function where , and I A is the indicator function of event A, that is I A = 1 and I A c = 0.
Note that Using a reparameterization [32] of p as (1 − θ ) the likelihood function in Equation (8) can be rewritten as This shows that the estimation for multinomial data is equivalent to estimating K − 1 binomial response models. We can consider any link function for binary data, taking θ k = g −1 (η k ). For the MLE, . For Bayesian estimation, we generate a sample from the posterior distribution of each θ k , then we can do the transformation to obtain the estimators of p. Considering the Weibull link function, we need to generate a sample from the posterior of γ k and β k , for each k = 1, . . . , K − 1, and then perform the proper transformation to obtain the sample from the posterior of θ k . In this case, the prior of θ k can be viewed as a transformation of the priors of γ k and β k . Thus, for both MLE and Bayesian estimator, we can use the procedures described in Section 3.1. The partition scheme presented to solve the multinomial model estimation is intuitive. For more details about the reparameterization used here, see Pereira and Stern [32].

Model Selection and Diagnostics
In the case of binomial data, to compare models within frequentist set up, we use the Akaike Information Criterion (AIC) [33] and the Bayesian Information Criterion (BIC) [34]. For Bayesian analysis, we use long established tool of Deviance Information Criterion (DIC) [35]. We omit the details of these popular tools for the sake of brevity. In addition, for Bayesian analysis we use the Pr(D|M) [36], where D is the observed data and M is the used model. Pr(D|M) is approximated by where θ i is the i-th sample from the posterior distribution of θ under model M, given the data D. This measure is directly related to the Bayes Factor (BF). If the interest is to evaluate the BF 10 of the model M 1 against M 0 , considering Pr(M 0 ) = Pr(M 1 ) = 0.5, we have that BF 10 = Pr(D|M 1 )/Pr(D|M 0 ). For both Bayesian and frequentist paradigms, we use: a version of Kolmogorov-Smirnov statistics (KS) as measure of goodness of fit (KS is defined as KS = max i |y i − y i |, the maximum absolute error of the predicted and the observed frequencies, where y i is the predicted value of y i ); the Mean Absolute Error (MAE) defined as MAE = 1 n ∑ n i=1 |y i − y i |; and the Brier Score (B-S) defined by Brier [37].

Data Example
In the data examples, we compare the proposed link function with some others links. Table 1 presents the link functions considered.

Link Function Parameteric Space
Weibull AEP is the asymmetric exponential power link from Naranjo et al. [17]; Γ(·) is the mathematical gamma function; F G (·, a, b) is the distribution function of a random variable with distribution Gamma with shape a and scale b; F B (·, λ 1 , λ 1 ) is the distribution function of a random variable with distribution Beta(λ 1 , λ 2 ); F N (·) is the distribution function of a random variable with distribution normal, with mean zero and variance 1; and F SN (·, δ) is the distribution function of a random variable with distribution skew-normal, with mean zero, variance 1 and asymmetric parameter δ [38].

Binary Data Example
We analyze the study of relative potency of three different poisons: Rotenone, Deguelin and Mixture [24]. The experiment was to test the different poisons with different doses, with objective to understand the potency of the poisons. Five doses for rotenone, six doses for deguelin and six doses for the mixture were considered. For each dose and poison, around 50 insects were considered by observing how many insects were killed. The data are presented in Table 2. We consider that the response variable is binary with Y = 1 representing the insect killed, and as covariates: X 1 as the log (Dose), X 2 as an indicator of Rotonone, and X 3 as an indicator of Deguelin. The mixture of poisons is considered as the reference poison (that is, X 2 = 0 and X 3 = 0). Our main objective is to find the model that better represents (fits) these Data. We are not looking for the "best" poison or dose. We obtain the MLE for Weibull parameters and for comparison we also estimated the parameters of complementary log-log, Stukel, probit, logit, Aranda-Ordaz, and Prentice models. Table 3, presents some statistics of each model to compare them. The best models, based on AIC, are complementary log-log and Weibull. The advantage of the complementary log-log is that this model has one fewer parameter. However,γ = 114.5084 (Table 4 Table 4. Another important models are the skew-probit proposed by Chen et al. [5] and AEP proposed by Naranjo et al. [17]. To compare with the Weibull model we perform a Bayesian analysis for these models. The priors for the parameters of the Weibull model are the same as that described in Section 3.1. We use the values v γ = 100 and v β = 25. The estimated values for the hyper-parameters of first hierarchical level are m γ = 9.1089 and m β = (0.1588, 0.8879, 0.1261, −0.1717). For the priors of skew-probit model, we used a uniform distribution over the interval (−1, 1) for the asymmetry parameter and independent normal distribution with mean 0 and variance 25 for each β j , j = 0, . . . , 3. For the priors of AEP model, we used independent priors: gamma distribution with mean 1 and variance 100 for the parameters θ 1 and θ 2 , and normal distribution with mean 0 and variance 25 for each β j , j = 0, . . . , 3. Table 5 presents the model selection criteria (DIC, KS, MAE, B-S and Pr(D|M)) for the three competing models. For criteria DIC, KS, MAE and B-S, a smaller value indicates a better agreement between the model and observed data. For criterion Pr(D|M), a higher value indicates a better agreement between the model and observed data. We note that the Weibull model has better KS; AEP has the better DIC and MAE; and skew-probit has better Pr(D|M). The B-S was very similar for all models. For the three models, the posterior mean of relevant parameters are given in Table 6.  Table 6. Bayesian estimates for binomial example.

Multinomial Data Example
Grazeffe et al. [25] reported a study of DNA mutation of the cells of adult snails, each irradiated with a single dose of gamma radiation. They recorded four categories of DNA mutation with Y = 1, 2, 3 and 4 representing no mutation, low, intermediate, and high DNA mutation respectively. The snails are randomized into five different dose levels with X ∈ {0, 2.5, 5, 10, 20}. The data are presented in Table 1 of Grazeffe et al. [25]. The objective is to compare effects of different dose levels on DNA mutation Y (Y = 1 for C0, Y = 2 for C1, Y = 3 for C2 and Y = 4 for C3). We illustrate the use of Weibull link model, under frequentist approach, for analysis of this study with multinomial responses. Further, in Table 10, we compare our estimates of Pr[Y = k | x] with those obtained by Grazeffe et al. [25] based on the logit link model, and the other models discussed here.
For a proper comparison with previous method of Grazeffe et al. [25], we obtain the MLE with only X and X 2 as covariates. We use the reflected Weibull link, because this model has lowest values of KS and MAE than those for Weibull link. To obtain the estimation of the reflected Weibull model we first estimate the values of θ 1 , θ 2 , θ 3 . To simplify, consider the three binary variables Z 1 , Z 2 and Z 3 , where θ k = Pr(Z k = k), k = 1, 2, 3. Then, using the results in Section 3.2, we construct Table 7 with the observed values of Zs, and estimate the models for Zs. The parameter estimates for the three binary models are presented in Table 8, and we have θ 1 (x) = e −(0.0234−1.6395x+0.6748x 2 ) 0.1742 , θ 2 (x) = e −(1.0930−0.0368x+0.0030x 2 ) 2.3604 and θ 3 (x) = e −(1.2429−0.0866x+0.0047x 2 ) 1.7562 .   Table 9 presents the inferential statistics for model comparisons. All statistics indicate a preference for Weibull link model. The main difference for the Stukel model was because Weibull model has three parameters fewer than the Stukel model in this multinomial example. The estimated frequencies, under Weibull, Stukel and logit models, of DNA mutation for each class is presented in Table 10, and illustrated in Figure 5. This figure shows that the Weibull link model has a better fit for categories Y = 1 and 4, when compared with the logit link model. For categories Y = 2 and 3, both models have comparable performances. Weibull and Stukel models have similar values of estimated frequencies.

Final Comments
In this paper, we have presented a Weibull model to estimate the problem of binary and multinomial regression analysis. The model is very flexible and capable to handle with many different types of data. The comparison with other skew-link model, in binomial data example (Section 4.1), shows that the performance of the Weibull link was good when compared to the others models. The model with worst measures was the Prentice model. All others had an equivalent result. We are convinced that our proposed model is a good option. A good feature of the model is that the logit, probit, complementary log-log, and log-log link functions are approximations of Weibull link. Then, the proposed model can accommodate even symmetric link function. For the flexibility of the Weibull link model, we are comfortable to suggest its use in practice.
Other aspect of the proposed Weibull model is that the associated numerical procedure of MLE is very simple to implement, particularly in comparison to other competing. For Bayesian estimates, we also suggest an Empirical Bayes approach to determine the prior. Under full Bayesian estimation, we compare the model with the skew-probit model [5] and AEP model [17], in Section 4.1. Again, all models had similar results, however the KS of Weibull model were the measures with the greatest differences among all models. The performance of our model was good, even under full Bayesian framework, in binomial data example (Section 4.1).
We also develop a partition scheme for the multinomial regression model simplifying the problem to K − 1 binomial regression analysis. This is a general scheme that can be used for other link functions, which opens a vast options to estimate multinomial data. In Section 4.2, we analyze a multinomial data problem, where the Weibull model had the best measure values when compared with all other models. Our perspective is that the Weibull model is a good option for binary/multinomial regression, mainly due to its simplicity. We have analytic form for the link function, as well as for the gradient and Hessian matrix.