Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
82,335 characters · 9 sections · 111 citation commands
Bayesian panel quantile regression for binary outcomes with correlated random effects: An application on crime recidivism in Canada
The concept of quantile regression introduced in Koenker-Basset-1978 has captured the attention of both statisticians and econometricians, theorists as well as applied researchers, and across school of thoughts i.e., Classicals (or Frequentists) and Bayesians. Quantile regression offers several advantages over mean regression (such as robustness against outliers, desirable equivariance properties, etc.) and estimation methods, particularly for cross-section data, are also well developed\footnote{Some Classical techniques include simplex method Dantzig-1963, Dantzig-Thapa-1997,Dantzig-Thapa-2003, Barrodale-Roberts-1973, Koenker-dOrey-1987, interior point algorithm Karmarkar-1984,Mehrotra-1992 and smoothing algorithm Madsen-Nielsen-1993, Chen-2007. Bayesian methods using Markov chain Monte Carlo (MCMC) algorithms for estimating quantile regression was introduced in Yu-Moyeed-2001 and refined, amongst others, in Kozumi-Kobayashi-2011. A non-Markovian simulation based algorithm was proposed in Rahman-2013. See also Soares-Fagundes-2018 for interval quantile regression using swarm intelligence.}. The method has been employed in various disciplines including economics, finance, and the social sciences KoenkerBook-2005, Davino-etal-2013. However, the development of quantile regression for panel data witnessed noticeable delay (more than two decades) because of complexities in estimation. The primary challenge was that quantiles, unlike means, are not linear operators and hence standard differencing (or demeaning) methods are not applicable to estimation of quantile regression. The challenges in estimation increases, if, for example, the outcome variable is discrete (such as binary or ordinal) because quantiles for such variables are not readily defined. Besides, modeling of panel data brings in consideration of unobserved individual-specific heterogeneity and the related debate on the choice of “random-effects” versus “fixed-effects”. Motivated by these challenges in modeling and estimation, this paper considers a quantile regression model for panel data in the presence of correlated-random effects (CRE) and introduces two Markov chain Monte Carlo (MCMC) algorithms for its estimation. The proposed framework is applied to study crime recidivism in the Province of Quebec, Canada, using a novel data constructed from the administrative correctional files.
The current paper touches on at least two growing econometric/statistic literatures -- quantile regression for panel data and panel quantile regression for discrete outcomes. In reference to the former, Koenker-2004 was first to suggest a penalization based approach to estimate quantile regression model with unobserved individual-specific effects\footnote{For other development in quantile regression on panel data see, amongst others, Lamarche-2010, Canay-2011, Chernozhukov-etal-2013, Galvao-etal-2013, Galvao-Kato-2017, Graham-etal-2018, and Galvao-Poirier-2019 to mention a few.}. Geraci-Bottai-2007 adopted the likelihood based approach of Yu-Moyeed-2001 and constructed a working likelihood using the asymmetric Laplace (AL) distribution. They proposed a Monte Carlo expectation-maximization (EM) algorithm to estimate the panel quantile regression model and apply it to study labor pain data reported in Davis-1991. Later, Geraci-Bottai-2014 extended the panel quantile regression model of Geraci-Bottai-2007 to accommodate multiple individual-specific effects and suggested strategies to reduce the computational burden of the Monte Carlo EM algorithm. A Bayesian approach to estimate the panel quantile regression was presented in Luo-Lian-Tian-2012, where they propose a Gibbs sampling algorithm by exploiting the normal-exponential mixture representation of the AL distribution Kozumi-Kobayashi-2011. Wang-2012 also utilized the AL density to develop a Bayesian estimation method for quantile regression in a parametric nonlinear mixed-effects model.
The papers on quantile regression mentioned in the previous paragraph have assumed that the unobserved individual-specific effects are uncorrelated with the regressors -- also known as “random-effects” in the Classical econometrics literature. In contrast, when the individual-specific effects are assumed to be correlated with the regressors, the models have been termed as “fixed-effects” model. Fixed-effects models suffer from the limitation that it cannot estimate the coefficient for time-invariant regressors. So, when most of the variation in a regressor is located in the individual dimension (rather than in the time dimension), estimation of coefficients of time varying regressors may be imprecise. Most disciplines in applied statistics, other than econometrics, use the random-effects model Cameron-Trivedi-2005. However, as shown in Baltagi-2013, most applied work in economics have settled the choice between the two specifications using the specification test proposed in Hausman-1978.
Between the questionable orthogonality assumption of the random-effects model and the limitations of the fixed-effects specification, lies the idea of correlated random-effects (CRE). This concept is utilized in the current paper to soften the assertion of unobserved individual heterogeneity being uncorrelated with regressors. The CRE was introduced in Mundlak-1978, where he models the individual-specific effects as a linear function of the time averages of all the regressors. Hausman-Taylor-1981 proposed an alternative specification in which some of the time-varying and time-invariant regressors are related to the unobserved individual-specific effects.\footnote{Baltagi-etal-2003 suggested an alternative pretest estimator based on the Hausman-Taylor (HT) model. This pretest alternative considers an HT model in which some of the variables, but not all, may be correlated with the individual effects. The pretest estimator becomes the random-effects estimator if the standard Hausman test is not rejected. The pretest estimator becomes the HT estimator if a second Hausman test (based on the difference between the FE and HT estimators) does not reject the choice of strictly exogenous regressors. Otherwise, the pretest estimator is the FE estimator.} Later, Chamberlain1982, Chamberlain-1984 considered a richer model and defined the individual-specific effects as a weighted sum of the regressors. These CRE models lead to an estimator of the coefficients of the regressors that equals the fixed-effects estimator. The literature has numerous publications on the Hausman tests or the CRE models in a linear or non-linear framework. We refer the reader to Baltagi-2013, Wooldridge-2010, Arellano-1993, Burda-Harding-2013, Greene-2015 and references therein. Most recently, Joshi-Wooldridge-2019 extended the CRE approach to linear panel data models when instrumental variables are needed and the panel is unbalanced.
Within the quantile regression for panel data literature, Abrevaya-Dahl-2008 incorporated the CRE to the quantile panel regression model and utilized it to study birth weight using a balanced panel data from Arizona and Washington. They make certain simplifying assumptions which allows them to estimate the model using pooled linear quantile regression. Following the quantile regression framework of Abrevaya-Dahl-2008, Bache-etal-2013 considers a more restricted specification to model birth weight using an unbalanced panel data from Denmark. Arellano-Bonhomme-2016 introduced a class of QR estimators for short panels, where the conditional quantile response function of the unobserved heterogeneity is also specified as a function of observables. The literature on Bayesian panel quantile regression with CRE is limited to Kobayashi-Kozumi-2012, where they develop Bayesian quantile regression for censored dynamic panel data and proposed a Gibbs sampling algorithm to estimate the model. The initial condition problem arising due to the dynamic nature of the model was successfully managed using correlated random effects. In addition, they implement the framework to study
The literature on panel quantile regression for discrete outcomes is quite sparse and most of the work has only come recently\footnote{A body of work related to quantile regression for discrete outcomes include, but is not limited to, Kordas-2006, Benoit-Poel-2010, Alhamzawi-2016, Omata-etal-2017, Alhamzawi-Ali-2018 and Rahman-Karnawat-2019}. Alhamzawi-Ali-2018 extended the Bayesian ordinal quantile regression introduced in Rahman-2016 to panel data and use it to analyze treatment related changes in illness severity using data from the National Institute of Mental Health Schizophrenia Collaborative (NIMHSC), and previously analyzed in Gibbons-Hedeker-1994. Ghasemzadeh-etal-2018-METRON proposed a Gibbs sampling algorithm to estimate Bayesian quantile regression for ordinal longitudinal response in the presence of non-ignorable missingness and use it to analyze the Schizophrenia data of Gibbons-Hedeker-1994. Ghasemzadeh-etal-2018-Comm developed a Bayesian quantile regression model for bivariate longitudinal mixed ordinal and continuous responses to study the relationship between reading ability and antisocial behavior amongst children using the Peabody Individual Achievement Test (PIAT) data. Most recently, Rahman-Vossmeyer-2019 considered a panel quantile regression model with binary outcomes and develop an efficient blocked sampling algorithm. They apply the framework to study female labor force participation and home ownership using data from the Panel Study of Income Dynamics (PSID).
This article contributes to the two literatures by incorporating the CRE concept into the panel quantile regression model for binary outcomes. Our proposed framework is more general and can accommodate the binary panel quantile regression model of Rahman-Vossmeyer-2019 as a special case. We present two MCMC algorithms -- a simple (non-blocked) Gibbs sampling algorithm and another blocked Gibbs sampling algorithm that exploits the block sampling of parameters to reduce the autocorrelation in MCMC draws. We also explain how to calculate the marginal effects, relative risk and the odds ratio using the MCMC draws. The performance of the blocked algorithm is thoroughly tested in multiple simulation studies and shown to perform extremely well. Lastly, we implement the model to study crime recidivism in the Province of Quebec, Canada, using data from the administrative correction files for the period 2007$-$2017. The results provide strong support for including the CRE into the binary panel quantile regression framework. On the applied side, we find that the recently implemented “tough-on-crime” policy has been successful in reducing the probability of repeat offenses and this is most pronounced at the lower quantiles. Besides, our results confirm existing findings from recent studies on crime recidivism, such as, schooling (unemployment rate) is negatively (positively) associated with crime recidivism. Moreover, the marginal effects and relative risk show considerable variability across the considered quantiles.
The remainder of the paper is organized as follows. Section (ref) introduces the binary panel regression model with correlated random-effects and the two MCMC algorithms. Section (ref) presents the simulation studies and discusses the performance of the algorithm. Section (ref) discusses how to compute the marginal effects, relative risk and odds ratio using the MCMC draws. Section (ref) implements the proposed framework to study crime recidivism in Quebec, a Canadian Province. Section (ref) presents concluding remarks.
We propose a binary quantile regression framework for panel data where the individual-specific effects are correlated with the covariates giving rise to correlated random effects. The resulting binary panel quantile regression with correlated random effects (BPQRCRE) model can be conveniently expressed in the latent variable formulation of Albert-Chib-2001 as follows,
where $z_{it}$ is a continuous latent variable associated with the binary outcome $y_{it}$, $x_{it}^{\prime} = (x_{it,1}, x_{it,2}, \cdots , x_{it,k})$ is a $\left( 1\times k\right) $ vector of explanatory variables including the intercept, $\beta $ is the $\left( k\times 1\right) $ vector of common parameters, and $\alpha _{i}$ is the individual-specific effect assumed to be independently distributed as a normal distribution, i.e., $\alpha_{i} \sim N\left( \overline{m}^{\prime}_{i} \zeta, \sigma _{\alpha }^{2}\right) $. Here $\overline{m}_{i,j}=\mathop{\textstyle \sum }_{t=1}^{T_{i}}x_{it,j}/T_{i}$ (for $j=2,...,k$) and $\overline{m}^{\prime}_{i} = (\overline{m}_{i,2},\cdots, \overline{m}_{i,k})$ is a $\left( 1\times (k-1)\right) $ vector of individual means of explanatory variables excluding the intercept. The dependence of $\alpha$ on the covariates $(x)$ yields a correlated random effects model Mundlak-1978. The error term $\varepsilon_{it}$, conditional on $\alpha_{i}$, is assumed to be independently and identically distributed (iid) as an Asymmetric Laplace (AL) distribution i.e., $\varepsilon_{it}|\alpha_{i} \overset{iid}{\sim} AL\left(0,1,p\right)$, where $p$ denotes the quantile. The AL error distribution is used to create a working likelihood and has been utilized in previous studies on longitudinal data models such as Luo-Lian-Tian-2012 and Rahman-Vossmeyer-2019.
In the proposed BPQRCRE framework, the modeling of correlated random effects as a function of the means of the covariates is inspired from Mundlak-1978. Utilizing $\overline{m}^{\prime}_{i}$ as a set of controls for unobserved heterogeneity is both intuitive and advantageous. It is intuitive because it estimates the effect of the covariates holding the time average fixed, and advantageous because it serves a compromise between the questionable orthogonality assumptions of the random effects model and the limitation of the fixed effects specification which leads to the incidental parameters problem. The considered model reduces to the standard uncorrelated random effects case, if we set $\zeta=0$, i.e., assume $\alpha_{i}$ is independent of the covariates Rahman-Vossmeyer-2019. Here, we note that Chamberlain1982,Chamberlain-1984 allowed for correlation between $\alpha _{i}$ and the covariates $x_{it}^{\prime }$ (excluding the intercept) through a more general formulation: $\alpha_{i}\sim N\left( \mathop{\textstyle \sum }_{t=1}^{T_{i}}x_{it}^{\prime }\zeta _{t}, \sigma _{\alpha }^{2}\right)$. However, this approach is more involved for an unbalanced panel, particularly if endogeneity attrition is the reason for the panel to be unbalanced Wooldridge-2010. Besides, the correlated random effects specification has a number of virtues for nonlinear panel data models as underlined in Burda-Harding-2013 and Greene-2015. Hence, we prefer the approach presented in Mundlak-1978 compared to the method in Chamberlain1980,Chamberlain1982, Chamberlain-1984.
The BPQRCRE model as presented in equation (ref) can be directly estimated using MCMC algorithms, but the resulting posterior will not yield the full set of tractable conditional posteriors necessary for a Gibbs sampler. Therefore, as done in Luo-Lian-Tian-2012 and Rahman-Vossmeyer-2019, we utilize the normal-exponential mixture representation of the AL distribution to facilitate Gibbs sampling Kozumi-Kobayashi-2011. The mixture representation for $\varepsilon_{it}$ can be written as follows,
where $u_{it}\sim N\left( 0,1\right) $ is mutually independent of $ w_{it}\sim \mathcal{E}\left( 1\right) $ with $\mathcal{E}$ representing the exponential distribution and the constants are $\theta =\frac{1-2p}{p(1-p)}$ and $\tau ^{2}=\frac{2}{p(1-p)}$. The mixture representation gives access to the appealing properties of the normal distribution.
To implement the Bayesian approach, we stack the model across $i$. Define $z_{i}=(z_{i1},...,z_{iT_{i}})^{\prime }$, $ y_{i}=(y_{i1},\cdots,y_{iT_{i}})^{\prime }$, $X_{i}=(x'_{i1},\cdots,x'_{iT_{i}})^{\prime}$, $w_{i}=(w_{i1},\cdots,w_{iT_{i}})^{ \prime }$, $D_{\tau \sqrt{w_{i}}}=\tau \; \mathrm{diag}(\sqrt{w_{i1}},\cdots,\sqrt{ w_{iT_{i}}})^{\prime }$ and $u_{i}=(u_{i1},\cdots,u_{iT_{i}})^{\prime }$. The resulting hierarchical model can be written as,
where $\iota _{T_{i}}$ is a $\left( T_{i}\times 1\right) $ vector of ones and the last line in equation (ref) presents the prior distribution on the parameters. The notation $N_{k}(\cdot)$ denotes a multivariate normal distribution of dimension $k$ and $IG(\cdot)$ denotes an inverse-gamma distribution. We note that the form of the prior distribution on $\beta $ holds a penalty interpretation on the quantile loss function Koenker-2004. A normal prior on $\beta$ implies an $\ell_{2}$ penalty and has been used in Geraci-Bottai-2007, yuan2010, Luo-Lian-Tian-2012 and Rahman-Vossmeyer-2019.
By Bayes' theorem, we express the \textquotedblleft complete joint posterior\textquotedblright\ density as proportional to the product of complete likelihood function and the prior distributions as follows,
where the first line assumes independence between prior distributions and second line follows from the fact that given $z_{it}$, the observed $y_{it}$ is independent of all parameters because the second line of equation (ref) determines $y_{it}$ given $z_{it}$ with probability $1$. Substituting the distribution of the variables associated with the likelihood and the prior distributions in equation (ref) yields the following expression,
The complete joint posterior density in equation (ref) does not have a tractable form, and thus simulation techniques are necessary for estimation. Similar to Rahman-Vossmeyer-2019, we adopt a Bayesian approach due to the following two reasons.. First, the likelihood function of a discrete panel data model is analytically intractable which makes optimization difficult using standard hill-climbing techniques. Second, numerical simulation methods for discrete panel data models are often slow and difficult to implement as noted in Burda-Harding-2013 and others. The complete joint posterior distribution (equation (ref)) readily yields a full set of conditional distributions (outlined below) which can be readily employed to estimate the model using Gibbs sampling.
We can derive the conditional posteriors of the parameters and latent variables from the joint posterior density (ref) by a straightforward extension of the non-blocked sampling method presented in Rahman-Vossmeyer-2019. This is presented in Algorithm (ref), and the derivations of the conditional posterior densities can be found in the supplementary material. The parameters $\beta$ are sampled from an updated multivariate normal distribution. Similarly, the parameters $\alpha_i$ are sampled from an updated multivariate normal distribution. The latent weights $w_{it}$ are sampled element wise from a generalized inverse Gaussian ($GIG$) distribution Devroye-2014. The variance $\sigma^{2}_{\alpha}$ is sampled from an updated inverse-gamma ($IG$) distribution. The parameters $\zeta$ are sampled from an updated multivariate normal distribution. Last, the latent variable $z_{it}$ is sampled element wise from an univariate truncated normal ($TN$) distribution. Note that while drawing each of the parameters or latent variables, we hold the remaining quantities fixed as presented in Algorithm (ref).
The MCMC procedure presented in Algorithm (ref) exhibits the conditional posterior distributions for the parameters and latent variables necessary for a Gibbs sampler. While this Gibbs sampler is straightforward, there is potential for poor mixing of the MCMC draws due to correlation between ($\beta $, $\alpha_{i}$) and ($z_{i}$, $\alpha _{i}$). This correlation arises because the variables corresponding to the parameters in $\alpha _{i}$ are often a subset of those in $x_{it}^{\prime }$. Thus conditioning these items on one another leads to high autocorrelation in MCMC draws as demonstrated in Chib-Carlin-1999 and noted in Rahman-Vossmeyer-2019.
To avoid the high autocorrelation in MCMC draws, we present an alternative algorithm that jointly samples ($\beta $, $z$) in one block within the Gibbs sampler Rahman-Vossmeyer-2019. The details of our blocked sampler are described in Algorithm (ref), and the derivations of the conditional posterior densities are presented in the supplementary file. Specifically, $\beta$ is sampled marginally of $\alpha_i$ from a multivariate normal distribution. Then the latent variable $z_i$ is sampled marginally of $\alpha_i$ from a truncated multivariate normal distribution denoted by $TMVN_{B_i} $, where $B_i$ is the truncation region given by $B_i=\left( B_{i1} \times B_{i2} \times ... \times B_{iT_i} \right)$ such that $B_{it}$ is the interval $\left(0, \infty \right)$ if $y_{it}=1$ and the interval $\left(-\infty, 0 \right]$ if $y_{it}=0$. To draw from a truncated multivariate normal distribution, we utilize the method proposed in geweke1991,geweke2005; as done in Rahman-Vossmeyer-2019. This involves drawing from a series of conditional posteriors which are univariate truncated normal distributions. The parameter $\alpha_i$ is sampled conditional on $(\beta, z, w, \sigma^{2}_{\alpha},\zeta)$ from an updated multivariate normal distribution. The latent weights $w_{it}$ are sampled element wise from a generalized inverse Gaussian ($GIG$) distribution Devroye-2014. The variance $\sigma^{2}_{\alpha}$ is sampled from an updated inverse-gamma ($IG$) distribution. Lastly, the parameters $\zeta$ are sampled from an updated multivariate normal distribution. Once again, while sampling each quantity of interest, we hold the remaining parameters or latent variables fixed as exhibited in Algorithm (ref).
In this section, we present two simulation studies to demonstrate the performance of the blocked algorithm for the BPQRCRE model. The simulation data are generated from the following model,
where $x^{\prime}_{it}=\left[ 1, \, x_{it,2}, \, x_{it,3}, \, x_{it,4} \right]$, $\overline{m}^{\prime}_{i} =\left[ \overline{m}_{i,3}, \, \overline{m}_{i,4} \right]$, $\overline{m}_{i,j}=\sum_{t=1}^{T_{i}}x_{it,j}/T_{i}$, $j=3,4$, $\beta=\left(\beta_{1}, \, \beta_{2}, \, \beta_{3}, \, \beta_{4} \right)^{\prime}=\left(0.5, \, 1, \, 0.6, \, -0.8\right)^{\prime}$, $\zeta=\left(\zeta_3, \, \zeta_4 \right)^{\prime} =\left(-1, 1\right)^{\prime}$. The covariates are generated as $x_{it,2} \sim U(-2,2)$, $x_{it,3} \sim U(-2,2)$, $x_{it,4} \sim U(-2,2)$, where $U$ denotes a uniform distribution, and $\sigma_{\alpha }^{2}=1$. Our first sample is unbalanced with $n=1,000$ and $T_i \sim U(5,15)$, leading to $T=\sum^{n}_{i=1} T_{i} = 9,989$ observations. In a second exercise, we increase the number of individuals $n = 2,000$ leading to $T=19,985$ observations. The error term is generated from a standard AL distribution, i.e., $\varepsilon_{it} \sim AL(0, 1, p)$ for $i=1,\cdots,n$, and $t=1,\cdots,T_i$ at three different quantiles $p$ = $0.25$, $0.5$, $0.75$.
The binary outcome variable $y$ is constructed from the continuous variable $z$, by assigning $y_{it}=1$ whenever $z_{it}>0$ and $y_{it}=0$ whenever $z_{it} \le 0$ for all of $i=1,\cdots,n$ and $t=1,\cdots,T_{i}$. We note that the binary response values of $0$s and $1$s are different at each quantile, because the error values generated from an AL distribution are different for each quantile. In the first simulation exercise with $n=1,000$, the number of observations corresponding to $0$s and $1$s for the $25$th, $50$th and $75$th quantiles are $(2283, 7706)$, $(4217, 5772)$ and $(6442, 3547)$, respectively. In the second simulation exercise with $n=2,000$, the number of observations corresponding to $0$s and $1$s for the $25$th, $50$th and $75$th quantiles are $(4640, 15345)$, $(8691, 11294)$ and $(13234, 6751)$, respectively. To complete the Bayesian setup for estimation, we use the following independent prior distributions: $\beta \sim N_{k}\left( 0_{k} , 10^3 I_{k} \right)$, $\zeta \sim N_{k-1}\left( 0_{k-1} , 10^3 I_{k-2} \right)$, $\sigma_{\alpha }^{2} \sim IG \left( 10/2 , 9/2 \right)$. For each exercise, we generate $16,000$ MCMC samples where the first $1,000$ values are discarded as burn-ins. The posterior estimates are reported based on the remaining $15,000$ MCMC iterations with a thinning factor of $10$. The mixing of the MCMC chain is extremely good as illustrated in Figure (ref), which reports the trace and autocorrelation plots of the parameters from the second simulation exercise at the 75th quantile. The figure shows that, as desired, the chains mix well and the autocorrelation of the MCMC draws are close to zero. The plots from the first simulation exercise and the remaining quantiles in the second simulation exercise are extremely similar and not presented to avoid repetition and keep the paper within reasonable length. To supplement the plots in Figure 1, Table (ref) presents the autocorrelation in MCMC draws at lag $1$, lag $5$, and lag $10$ confirming the good mixing across simulation exercises and at all quantiles.
The results from the two simulation exercises are presented in Table (ref). Specifically, the table reports the true values of the parameters used to generate the data, along with the posterior mean, standard deviation and inefficiency factor Greenberg-2012 of the MCMC draws. In general, the results show that the posterior means for $(\beta, \zeta)$ are near to their respective true values, $\beta=\left(0.5,1, 0.6, -0.8\right)^{\prime}$ and $\zeta=\left(-1,1\right)^{\prime}$ across all considered quantiles. The posterior standard deviations for all the parameters are small and all the coefficients are statistically different from zero. So, the proposed MCMC algorithm is successful in correctly estimating all the model parameters across all quantiles. This is especially important because the number of $0$s and $1$s were different for each quantile. Moreover, the inefficiency factor for all the parameters is close to $1$, suggesting a good sampling performance and a nice mixing of the Markov chain. Comparing the results from the first and second simulation exercise, we see that when the sample size is increased from ($n=1,000$, $T=9,989$) to ($n=2,000$, $T=19,985$), the results improve and the posterior means of the coefficients are closer to their true values. In particular, some small observed biases for $\beta_1$, $\zeta_3$, and $\zeta_4$ at the 25th quantile are reduced to a large extent. To summarize, the proposed algorithm for estimating BQQRCRE model does well in both the simulations, but the advantages of having a larger data is clearly evident in the posterior results.
Our proposed binary panel quantile model is nonlinear, as such the coefficients by themselves do not give the marginal effects Rahman-2016,Rahman-Vossmeyer-2019. However, marginal effects are important to understand the effect of a covariate on the probability of success. For example, in our current application one may be interested in seeing how the probability of recidivism is affected due to an additional year of schooling, decreasing regional unemployment rate by 1 percentage, or involvement in violent crime. These may be useful to policy makers and researchers alike.
To formally derive the marginal effects, we rewrite the BPQRCRE model presented in Equation (1) as follows,
where $\varepsilon_{it} = w_{it}\theta + \tau \sqrt{w_{it}} u_{it}$. We know $\varepsilon_{it} \overset{iid}{\sim} AL(0,1,p)$ for $i=1,\cdots,n$ and $t=1,\cdots,T_{i}$, which implies $z_{it}|\alpha_{i} \overset{ind} \sim AL(x'_{it} \beta + \alpha_{i},1,p)$, where $ind$ denotes independently distributed.
Given the model framework, the probability of success can be calculated as,
for $i=1,\cdots,n$ and $t=1,\cdots,T_{i}$, where $F_{AL}(x,0,1,p)$ denotes the cumulative distribution function (cdf) of an AL distribution evaluated at $x$, with location 0, scale 1 and quantile $p$.
Marginal effect (i.e., the derivative of the probability of success with respect to a covariate) is often computed at the average covariate values or by averaging the marginal effects over the sample, alias average partial effects Wooldridge-2010,Greene-2017. However, Jeliazkov-Vossmeyer-2018 show that both these quantities can be clearly inadequate in nonlinear settings (e.g., binary, ordinal and Poisson models) because they employ point estimates rather than their full distribution. To account for the uncertainty in parameters, we need another layer of integration over the model parameters. This idea of calculating the marginal effect that accounts for uncertainty in parameters and the covariates has been previously considered, amongst others, by Chib-Jeliazkov-2006 in the context of semiparametric dynamic binary longitudinal models, and Jeliazkov-etal-2008 and Jeliazkov-Rahman-2012 in relation to ordinal and binary models. Within the quantile literature, this has been mentioned by Rahman-2016 in the context of ordinal models and discussed by Rahman-Vossmeyer-2019 in connection to binary longitudinal outcome models.
Suppose, we are interested in the average marginal effect i.e., average difference between probabilities of success when the $j$-th covariate $\{x_{it,j}\}_{t=1}^{T_{i}}$ is set to the values $a$ and $b$, denoted as $\{x_{it,j}^{a}\}_{t=1}^{T_{i}}$ and $ \{ x_{it,j}^{b}\}_{t=1}^{T_{i}}$, respectively. To proceed, we split the covariate and parameter vectors as follows: $x_{it}^{a} = ( x_{it,j}^{a}, x_{it,-j})$, $x_{it}^{b} = ( x_{it,j}^{b}, x_{it,-j})$, and $\beta = (\beta_{j}, \beta_{-j})$, where $-j$ in the subscript denotes all covariates/parameters except the $j$-th covariate/parameter. We are interested in the distribution of the difference $\{\Pr(y_{it}=1|x_{it,j}^{b}) - \Pr(y_{it}=1|x_{it,j}^{a} ) \}$, marginalized over $\{x_{it,-j}\}$ and the parameters $(\beta,\alpha)$, given the data $y=(y_{1}, \cdots, y_{n})'$. As done in Chib-Jeliazkov-2006 and Rahman-Vossmeyer-2019, we marginalize the covariates using their empirical distribution and integrate the parameters using their posterior distribution.
To obtain a sample of draws from the distribution of the difference in probabilities of success, marginalized over $\{x_{it,-j}\}$ and $(\beta,\alpha)$, we express it as follows,
Drawing a sample from the above predictive distribution (i.e., equation (ref)) utilizes the method of composition. This involves randomly drawing an individual, extracting the corresponding sequence of covariate values, drawing a value $(\beta, \alpha)$ from the posterior distribution and finally evaluating $\{\Pr(y_{it}=1|x_{it,j}^{b}) - \Pr(y_{it}=1|x_{it,j}^{a} ) \}$. This is repeated for all other individuals and other draws from the posterior distribution. Finally, the average marginal effect $(AME_{Bayes})$ is calculated as the average of the difference in pointwise probabilities of success as follows,
where the expression for probability of success follows from equation (ref), $T=\sum_{i=1}^{n} T_{i}$ is the total number of observations, and M is the number of MCMC draws. Here, $(\beta^{(m)}, \alpha^{(m)})$ is an MCMC draw of $(\beta, \alpha)$ for $m=1,...,M$. The quantity in equation (ref) provides estimate that integrates out the variability in the sample and the uncertainty in parameter estimation.
Relative risk ($RR$) can be calculated to demonstrate the association between the risk factor or exposure ($x_j$) and the event ($y$) being studied. It is the ratio of the probability of the outcome with the risk factor ($x_j=b$) to the probability of the outcome with the risk factor ($x_j=a$) (e.g., exposed ($b=1$) /non-exposed ($a=0$)). Following equation (ref), the relative risk is given by,
where $H_{AL}^{r} = 1 - F_{AL}(-x_{it,j}^{r} \beta_{j}^{(m)} - x'_{it,-j} \beta_{-j}^{(m)} - \alpha_{i}^{m},0,1,p)$ for $r = a, b$, is the complement of the cdf of the AL distribution. If there is a causal effect between the exposure and the outcome, values of $RR$ can be interpreted as follows: if $RR>1$ (resp. $RR<1$), the risk of outcome is increased (resp. decreased) by the exposure and if $RR=1$, the exposure does not affect the outcome.
The odds ratio is the ratio of the odds of the event occurring with the risk factor ($x_j=b$) to the odds of it occurring with the risk factor ($x_j=a$). It is given by:
The odds ratio, for a given exposure $x_j$, does not have an intuitive interpretation as the relative risk. OR are often interpreted as if they were equivalent to relative risks while ignoring their meaning as a ratio of odds. Two main factors influence the discrepancies between $RR$ and $OR$: the initial risk of an event $y_{it}$, and the strength of the association between exposure $x_{it,j}$ and the event $y_{it}$. When the event $y_{it}=1$ is rare, then $OR(b/a)$ $\approx $ $RR(b/a)$, but the odds ratio generally overestimates the relative risk, and this overestimation becomes larger with increasing incidence of the outcome.
Crime has been extensively studied by economists both theoretically and empirically (see, e.g.chalfin2017 for a recent survey). Many empirical analyses have used panel data either at the state Cornwell1994,baltagi2006,baltagi2017 or at the individual level bhuller2019. The vast majority of the published papers focus on the situation in the U.S. Here, we study crime recidivism in Canada between 2007-2017 for two reasons. First, the Canadian government implemented a “tough-on-crime" policy in 2012 which marked a shift from rehabilitating to warehousing people. Our proposed estimator is well suited to measure the sensitivity of recidivism to this new policy.\footnote{Starting in 2012, the government enacted a series of legislations that made prison conditions more austere; imposed lengthier incarceration periods; significantly expanded the scope of mandatory minimum penalties; and reduced opportunities for conditional release, parole, and alternatives to incarceration.} Second, offenders who are sentenced to less than two years serve their sentence in a provincial correctional institution while offenders sentenced to two years or more serve their's in a federal penitentiary. The former have committed less serious crimes and are more likely to reoffend over the time span of our panel. Because our analysis focuses on this population, the impact of the “tough-on-crime” policy may be more easily unearthed from the data than if it focused on detainees serving long sentences.
We utilize a sample data drawn from the administrative correctional files for the Province of Quebec. The files are used by corrections personnel to manage activities and interventions related to housing offenders and contain detailed information on inmates' characteristics, correctional facilities, and sentence administration. While they offer a wealth of information, the files have never been used for research purposes. For illustrative purpose, we have drawn a random sample of 8,974 detainees out of a population of 148,441. Each detainee is observed upon release and up until 2017. The earliest releases occur in 2007 and the latest in 2016. Overall, our unbalanced panel includes 61,880 observations. Of the 8,974 detainees, as many as 3,466 had at least one repeat offense over our sample period.
Table (ref) presents the main characteristics of our sample. Detainees are 41 years of age on average, have a level of schooling corresponding to a high-school degree, and few are married. Aboriginal detainees represent 4.5% of our sample and most are incarcerated in a correctional institution suited to their needs and specificities. Approximately 7% of inmates do not have French or English, Canada's two official languages, as their mother tongue. These include some Aboriginal residents as well as recent immigrants. Crimes have been aggregated into 4 distinct categories. By far the most common concerns property crime. Traffic related and infractions to the criminal code usually entail shorter sentences. Violent crimes receive the longest sentences in our data but necessarily less than two years. As mentioned above, major crimes fall under the federal jurisdiction. The yearly unemployment rate is measured at the regional level where a detainee is released. Over our sample period, it varies between 4.4% and 17.5%. The “Post 2012” variable is equal to one if a detainee entered the panel at any time during or after 2012 while the “Pre-Post 2012” variable is equal to one if a detainee entered at any time before 2012. In the latter case, repeat offenses are observed over the entire duration of the panel, i.e. 2007-2017. In the former, they are only observed over 2012-2017. Roughly a quarter of our sample belongs to the period post the implementation of “tough-on-crime” policy. The remaining observations (74.8%) were sanctioned prior to 2012 and may or may not have reoffended in the Post 2012 period. The next 3 lines of the table provide information on the rates of recidivism for distinct periods.\footnote{Recidivism is a yearly dummy variable equal to one the year at which the new incarceration begins and zero otherwise. Recidivism may be equal to one in consecutive years so long as the repeat offenses occurred after the end of the previous sentence. Reincarcerations while on parole or on conditional release are not considered repeat offenses. } Thus, the overall rate of recidivism is equal to 11.4%. The next line focuses on individuals who are present both before and after the implementation of the “tough-on-crime” policy. Their recidivism rate is approximately 9%. The last line focuses on individuals who entered the panel on or after 2012. Naturally, as they are observed for a shorter period of time, their recidivism rate is relatively smaller at 2.3%.
Figure (ref) depicts the proportions of repeat offenses for the entire sample period and for those who entered the panel in 2012 or later. The figure provides prima facie evidence on the impact of the policy. Indeed, the proportion of detainees who do not reoffend upon release in the post-policy period is 15 percentage points larger (74.1%) than the proportion for the whole sample period (51.5%). Likewise, the proportion of repeat offenders is between 3 to 6 percentage points lower in the post-policy period for any given number of repeat offenses.\footnote{Obviously, detainees who entered the sample on or after 2012 have had less time to reoffend. Yet, in our sample as many as 34% of detainees are reincarcerated within 12 months upon release, and as many as 43% within two years. Hence, the sharp decline in repeat offenses in the post-2012 period is unlikely due to the sampling frame. See lalande2015.} Naturally, such differences may results from factors other than the “tough-on-crime" policy, such as, but not limited to, better economic opportunities, and demographic compositional changes. In order to net these out, we now turn to formal econometric modelling.\footnote{To the extent the new legislation has indeed lowered the recidivism rates, it not clear whether it did so through deterrent or incapacitative effects. Yet, see bhuller2019 for U.S. evidence according to which deterrence dominates incapacitation.}
\vskip-220pt
The dependent variable $y$ is an indicator variable that equals 1 if an individual commits a repeat offense and 0 otherwise. We regress the probability of recidivism on time-varying covariates (age, schooling, unemployment rate), on time-invariant policy variables (Pre-Post 2012 and Post 2012) and on other time-invariant control variables.\\ Our Bayesian setup uses the same independent prior distributions as in the simulation exercise: $\beta \sim N_{k}\left( 0_{k} , 10^3 I_{k} \right)$, $\zeta \sim N_{k-1}\left( 0_{k-1} , 10^3 I_{k-2} \right)$, $\sigma_{\alpha }^{2} \sim IG \left( 10/2 , 9/2 \right)$. We generate $60,000$ MCMC samples of which the first $10,000$ are discarded as burn-ins. The posterior estimates are reported using a thinning factor of $50$, optimized following the approach in Owen-2017.\footnote{Thinning has been criticized by some MacEachern-Berliner-1994,Link-Eaton-2012 while others acknowledge that it can increase statistical efficiency Geyer-1991. See Owen-2017 who claims that the arguments against thinning may be misleading.}
The mixing of the MCMC chain is extremely good as illustrated in Figure (ref) which exhibits the trace plots of the parameters at the 75th quantile.\footnote{Note that the time-varying covariates (Age, Schooling and Unemployment rate) have been “demeaned” and that Age has been divided by 10. The parameter estimates must thus be interpreted accordingly.} Trace plots at other quantiles are similar and not reported for the sake of brevity but they are available upon request. Figure (ref) provides additional information on the performance of the MCMC chain. The left-hand-side figure depicts the boxplots of the inefficiency factors of the parameters ($\beta$s, $\zeta$s and $\sigma^2_{\alpha}$) for each of the five different quantiles used in estimating the model. Except perhaps for the 10th quantile, all are reasonably close to one. Consistent with the simulation results, the parameter with the largest inefficiency factor at the 10th quantile is $\sigma^2_{\alpha}$ (not shown, see Table (ref)). The right-hand-side figure reports the boxplots of the convergence diagnostics of the parameter estimates for the same five specifications based on the first 10% and the last 40% values of the Markov chain geweke1992. As depicted, all parameters have $Z$-scores within 2 standard deviation of the mean at the $5\%$ level or within $2.58$ standard deviation at $1\%$ level. All in all, the Markov chains behave satisfactorily and thus lend themselves to statistical inference.
Table (ref) reports the posterior means and standard deviations at five different quantiles separately. To ease interpretation, the quantile-specific estimates are reported column-wise in increasing order. Row-wise, we distinguish the time-varying covariates from the time-invariant and the correlated random effects variables. Note that the correlated random effects specification does not include an intercept. This is to allow the identification of the two time-invariant policy variables, Pre-Post 2012 and Post 2012. The former, is equal to one if the detainee was incarcerated prior to 2012 and thus observed both before and after the implementation of the “tough-on-crime” policy. The latter is equal to one if a detainee's first incarceration occurred during or after 2012, and thus always exposed to the policy. All other time-invariant variables are measured at first entry in the panel.\footnote{Recall from Table (ref) that very few men are married. In addition, next to none report a change in their marital status in between incarcerations. Further, since the marital status of non-repeaters is not observed in the data we are constrained to use the information at entry in the panel.} The estimates of the correlated random components associated with the individual mean Age, Schooling and Unemployment, $\widehat{\zeta}$, are all statistically different from zero regardless of the quantile. The individual-specific effects, $\alpha_i$, are thus highly correlated with the individual means of the time-varying variables. Omitting this correlation may therefore bias the model estimates and hence their intrinsic marginal effects and relative risks. This provides empirical support to the worthiness of incorporating correlated random effects within a quantile regression.
The first noteworthy feature of the table is that all parameter estimates are statistically different from zero, except for the parameter associated with Other Mother Tongue. Thus detainees who report speaking a language other than English or French at home are no more and no less likely to eventually reoffend. A second interesting feature concerns the sign of the parameter estimates. Indeed, all are consistent with recent research on crime recidivism. For instance, Age and Schooling are associated with lower rates of recidivism bhuller2019 whereas being released during a period of high unemployment has been found to favour recidivism siwach2018,rege2019. Likewise, married men are less likely to reoffend whereas Aboriginal detainees are more likely to do so justice2017. The type of crime is also associated with recidivism. The estimates must be interpreted relative to traffic related crimes, which is the base or omitted category in our analysis. Clearly, sentences for Violent Crimes will be harsher and so the large parameter estimate presumably reflects an incapacitative effect. Finally, the parameter estimates of Post 2012 is larger than that of Pre-Post 2012 which suggests that the implementation of the “tough-on-crime” policy may have had a detrimental effect on recidivism.
As stated in Section (ref), the parameter estimates such as those reported in Table (ref) do not give the marginal effects. Yet, the latter are important from a policy perspective. Thus, while the parameter estimates vary considerably across quantiles, it is not clear that the marginal effects are equally sensitive since they depend both on the time-varying variables and the correlated random components. Figure (ref) reports the average marginal effects computed according to equation ((ref)), along with their highest posterior density intervals (HPDI).\footnote{The marginal effects for Age correspond to 1/10 of an additional year relative to the mean. Those for Unemployment and Schooling correspond to one additional year and one additional percentage point relative to their individual means, respectively. The remaining marginal effects correspond to a change in the indicator variables.} Note that most marginal effects have a relatively flat profile between $p10$ and $p75$ and then exhibit a small kink between $p75$ and $p90$. For instance, increasing Age by 1/10th reduces the probability of reoffending by 1% at the 10th quantile and by 1.6% at the 90th quantile. Similar results hold for Schooling (1% vs 2.0%), \texttt{Married} (0.3% \textit{vs} 0.45%), and \texttt{Violent Crime} (5% \textit{vs} 6.5%). Thus, for all three time-varying covariates the marginal effects increase by one half as we move from $p10$ to $p90$. As for the time-invariant variables, their marginal effects all increase by at least 50% as we move from $p10$ to $p90$. In particular, the marginal effects associated to \texttt{First Nation}, \texttt{Property Crime} and \texttt{Other Crime} exhibit a twofold increase. More importantly, the marginal effects of the two “tough-on-crime” variables increase manifold and in a steady fashion between $p10$ and $p90$. Furthermore, the HPDI is relatively narrow in both cases. Hence, according to the parameter estimates associated with \texttt{Pre-Post 2012}, the probability of reoffending decreases from 78% at the 10th quantile to as little as 10% at the 90th. Likewise, the parameters of \texttt{Post 2012} imply that the probability decreases from 79% to 14% at both extremes. These results are important from a policy perspective for two reasons. First, they imply that detainees from both groups are sensitive to the “tough-on-crime” policy, and even more so for those in the \texttt{Post 2012} group. Consequently long-run recidivism (\textit{i.e.} recidivism by the \texttt{Pre-Post 2012} group between 2012-2017) can be addressed just as well as short-run recidivism (\textit{i.e.} recidivism by the \texttt{Post 2012} group between 2012-2017) by such policies. Second, the policy does not impact all detainees alike. Those in the lower quantiles are much more responsive than those in the upper quantiles.
In order to gain further insight into the sensitivity of recidivism to various covariates, we report the corresponding relative risks in Figure (ref) (see equation ((ref))) along with their HDPI. Not surprisingly given the marginal effects, the relative risks are fairly constant for the first two or three quantiles ($p=10\%,25\%, 50\%$), with a few exceptions. Beyond the second or third quantiles, most increase or decrease sharply. The figure also shows which covariates influence recidivism most. Thus, while Age, Schooling and Unemployment Rate are associated with slightly different rates of repeat offenses, only those in the highest quantiles exhibit significantly different recidivism rates. On the other hand, marital status (Married), First Nation and types of crime (Violent, Property, Other) all have significantly higher or lower relative risks of reoffending as the case may be, and all exhibit a sharp change between the last two quantiles. Here, as with the previous figure, the results concerning the “tough-on-crime” variables are particularly interesting. Indeed, according to the figure all detainees were much less likely to reoffend in the post 2012 period, irrespective of whether they where first convicted prior to 2012 or after. As with the marginal effects, the policy appears to have had a larger impact on those in the lower quantiles. Thus for every quantile the risk of recidivism is much lower (and significantly different) for those who were exposed to the “tough-on-crime” policy. For instance, the 95% HPDI at quantile $p10$ is $\left[0.087;0.094\right]$ for the \texttt{Pre-Post 2012} group and $\left[0.066; 0.074\right]$ for the \texttt{Post 2012} group. On the other hand, the 95% HPDI at quantile $p90$ for the two groups are $\left[0.323;0.385\right]$ and $\left[0.189; 0.240\right]$, respectively. In other words, for the lowest quantile ($p10$), exposure to the policy decreases recidivism by as much as $\left[90;91\right]\%$ and $\left[92;93\right]\%$ for the \texttt{Pre-Post 2012} and \texttt{Post 2012} groups, respectively. In contrast, for those in the highest quantile, $p90$, the \texttt{Post 2012} group decreases its recidivism rate more than that of the \texttt{Pre-Post 2012} ($\left[76;81\right]\%$ vs $\left[61;67\right]\%$).
This paper presents a panel quantile regression model for binary outcomes with correlated random-effects (CRE) and proposes two MCMC algorithms for its estimation. By incorporating the CRE into the panel quantile regression for discrete outcomes, we move beyond the random-effects framework typically considered in the Bayesian quantile regression literature. The paper makes an important contribution to the literature on quantile regression for panel data and panel quantile regression for discrete outcomes. The two proposed MCMC algorithms are simpler to implement, but we prefer the algorithm that exploits block sampling of parameters to reduce the autocorrelation in MCMC draws. This blocked algorithm is tested in multiple simulation studies and shown to perform extremely well. We also emphasize the calculation of marginal effects in models with discrete outcome and explain its computation, along with those of relative risk and odds ratio, using the MCMC draws. Finally, we implement the proposed quantile framework to analyze crime recidivism in Quebec (a Canadian Province) for the period 2007$-$2017 using a novel data from the administrative correctional files. Amongst other things, we investigate the effect of the recently implemented “tough-on-crime” policy on the probability of repeat offense. Our results show that the policy negatively affects the probability of repeat offenses across quantiles and hence has been largely successful in achieving its objective. Besides, the results suggest that the CRE structure is relevant in modeling the probability of repeat offenses across quantiles.
This paper opens avenues for future research in several directions. The proposed framework can be readily extended to panel quantile regression models with continuous and other discrete response variables (e.g., count and ordinal outcomes). One may also consider the Hausman-Taylor version of CRE, where the individual-specific effects are related to only some of the time-varying and time-invariant regressors, and merge it with the panel quantile regression model for continuous or discrete outcomes. Besides, a dynamic relationship can be introduced to panel quantile regression models (with continuous or discrete outcomes) and the initial condition problems can be tackled using the CRE structure.
The authors declare that they have no conflict of interest.