I am having difficulty finding a correct glm model to fit my data.

The outcome is the length of time in months a person will spend in prison (sentence length). It's technically a count, all positive integer values, although there are no zeros (I modeled those separately). It is censored at 600 - 50 years, a "life" sentence. The data are badly overdispersed - mean is 53 and the variance is 7338.

I have approached the data using a heckman two-step model where I first fit a probit model to predict who gets prison and who does not and then calculated the inverse mills ratio which I included in the sentence length model to adjust for selection effects. The predictors are almost but not entirely the same across the two models. The differences in the models are meaningful, I didn't just "leave out" some variables from the selection model.

The data are hierarchical and I have hypotheses that relate to cross level interactions, but we can ignore this complexity for simplification purposes.

I am inclined to fit a negative binomial regression with a log link but the data are stripped of all the zeros, the lowest value possible is 12 (1 year in prison). I'm having a hell of a time getting the mixed effects negative binomial model to converge even using start values from the fixed effects model but that's another problem which is something I'll have to deal with if the negative binomial approach is the right one.

When I plot the deviance against the fitted value (mean response, not the linear prediction) from the negative binomial model with the log link I get this: glm negative binomial, log link deviance residuals vs. mean response

I think the large positive residuals at the low end of the predicted values are worrisome but am not sure what the remedy would be.

Now, the situation looks worse when modeled with a gaussian family and the log link. Here's the deviance vs. mean response for that model. These downward diagonal lines which I believe are indicative of the data being discrete are more pronounced and the residuals tend to be positive and low values and negative at high values. glm with gauss family and log link

These plots are from simplified fixed effects models without the interactions. I used Stata 14 to produce these results. I originally considered a gamma model with log link but that model is the least cooperative in converging and the deviance vs. fitted plot looks pretty much identical to the plot from the negative binomial model. A gaussian model with the identity link converges and produces coefficients that are literally 0 which is bizarre (I've never seen it) and must indicate incredibly poor model fit.

A negative binomial model on the whole set of data is not so great either, there are just so many zeros - about 80% of full the sample are zeros. And yes, it's a big dataset if you can't tell from the plots.

EDITED TO ADD: In hopes of getting a response, I'll add another plot that I think might be insightful. There is a variable I include in the analysis that reflects the minimum sentence a person is supposed to receive based on their overall "score" for the case (crime seriousness, victim injury, etc.). Here is a scatter plot of that score against the observed values and a regression line. You can see there's a clear linear relationship in the data and the fitted tracks it exactly parallel but the intercept places the line a bit high. You can also see the horizontal bars which I believe correspond to the diagonal lines in the residual vs. fitted plots. The bars correspond to widely used or common sentence terms. For example, the bar between 300 and 400 is 360 months - a 30 year sentence. There is also a pronounced line at 240 which would be a 20 year sentence. I'm guessing that these horizontal bars are the result of mandatory minimum sentences but they don't seem to be identified in the data, I'm not sure what to do about it.

Given that this is obviously one of the most important predictors in the model I'm thinking that I need to model the relationship as linear (because it looks linear) and not with a log link. This would rule out the negative binomial model at least as implemented in Stata's meglm command. But then I'm back to this issue with the data being censored at 12, strictly positive, and overdispersed. expected minimum sentence vs observed

Edited: I desperately need some advice about what direction to go. It seems like every glm specification is giving me biased estimates as indicated by the residuals being positive at low fitted values and negative at high fitted values. QQ plots likewise indicate non normality of the residuals which, I think, means that p values and confidence intervals are not accurate. The models just don't fit!

I've thought about boxcox transforming the predictor but am not sure how to go about that. I've thought about quantile regression which should be a bit more robust to non normal residuals. There's also hurdle models although I do not believe these will accommodate hierarchical data. There's also tobit, I could set 600 as an upper censor point but this doesn't resolve issues with the overdispersion of the data and the discrete nature of the outcome.

  • $\begingroup$ The latest issue of Stata Journal has a nice paper by Hardin and Hilbe on regression models for count data from truncated distributions that may prove helpful. $\endgroup$
    – dimitriy
    Commented Apr 25, 2015 at 0:32

1 Answer 1


Convergence problems aside, in any glm, you have two issues:

  1. How does your mean $E(Y)$ depend on covariates $X$? This is what should determine your link function.
  2. How does your variance $V(Y)$ depend on $E(Y)$. This is what should determine your family.

1.is more important than 2. because you can fix 2. by using robust variance (sandwich estimates), or bootstrapping. If you get 1. wrong, then your coefficient estimates become hard to interpret.

How about a quasi-Poisson model, with an identity link, and perhaps make your model a bit more flexible by including a regression spline for your prison length "score"?

  • $\begingroup$ Is there any merit in looking at mean and variance of the outcome across some of the predictors? Would this be helpful in determining what error family best fits the data? The CV looks to be about constant although it drops a bit for some groups who tend to get long sentences because of the upper limit censoring. This was the first thing I tried and it led me to the gamma family but the models don't like to converge and the residuals look just like those from the negative binomial model. I'll look at quasi poisson but I don't believe stata's glm will fit that family. $\endgroup$
    – whauser
    Commented Apr 25, 2015 at 13:26
  • $\begingroup$ I also wanted to add that I have experimented with linear splines. I haven't tried a restricted cubic spline on account of the difficult interpretation. In any case, one would think that the number of months a person gets would be a 1 to 1 identity link with the expected minimum sentence, but then there are those with these ridiculous point scores and resulting expected minimum sentence on account of various sentence "enhancements" and their sentences can't possibly be higher than 600 because of the censoring. $\endgroup$
    – whauser
    Commented Apr 25, 2015 at 13:39
  • $\begingroup$ Yes, I think that looking for trends in 1. residuals vs fitted, 2. residuals vs predictors and 3. squared residuals vs fitted (you will probably need to do some binning by fitted values here to get a graph that makes sense) is a sensible way to ID a link (items 1 and 2) and family (item 3). There are also some non-parametric ways to estimate links, at least, I believe under the heading of "projection pursuit" $\endgroup$
    – Andrew M
    Commented Apr 28, 2015 at 20:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.