Is there a binomial regression model that captures data with fat tails?

Specifically, are there any binomial regression models that use a kernel with heavier tails and higher kurtosis than the standard kernels (logistic/probit/cloglog)?

As a function of the linear predictor $\textbf{x}'\mathbf{\hat{\beta}}$, the logistic distribution

• Underestimates the probability of my data being in the tails of the distribution
• Underestimates the kurtosis, or clustering of data, in the middle of the distribution:

This can be seen from a diagnostic plot of my fit: • The red line is the logistic CDF, representing a perfect fit
• The black line represents the fitted probabilities from my dataset (calculated by binning observations into 0.1 intervals of $\textbf{x}'\mathbf{\hat{\beta}}$, where $\mathbf{\hat{\beta}}$ is obtained from my fit)
• The grey bars in the background represent number of observations on which the true probabilities are based upon
• The grey areas are where the tail 10% of the data lie (5% each side).

Ideally, any solution would use R.

Edit

Why am I talking about CDFs? Our GLM equation is:

$$\mathbb{P}(Y = 1) = \mathbb{E}[Y] = g^{-1}(\textbf{x}'\mathbf{\beta})$$

Where $g$ is the link function.

Further, if $g^{-1}$ is a valid probability distribution (i.e. monotonically increasing from 0 to 1, indeed the case with probit, logit, cloglog), then consider a latent (not directly observed) continuous random variable $Y^{*}$ whose distribution (CDF) is given by $g^{-1}$. Then by definition

$$\mathbb{P}(Y^{*} \leq \textbf{x}'\mathbf{\beta}) = g^{-1}(\textbf{x}'\mathbf{\beta})$$

Equating the two equations above, we see the probability of $Y=1$ is exactly equal to the CDF of $Y^{*}$

$$\mathbb{P}(Y = 1) = \mathbb{P}(Y^{*} \leq \textbf{x}'\mathbf{\beta})$$

Hence I talk interchangeably about the expected response $\mathbb{E}[Y]$ and CDF of $Y^{*}$ over linear-predictor ($\textbf{x}'\mathbf{\hat{\beta}}$) space.

• I believe you're confusing the CDF with the relationship between your predictor variable and the probability of a 1 in your binomial data. I don't think any of your lines represent a CDF. The black line definitely can't tell you anything directly about the CDF.
– John
Sep 3 '14 at 15:00
• You are confused; a binomial random variable cannot have fat tails! (it doesn't have tails at all) Sep 3 '14 at 15:17
• As John suggests, you're confounding possible lack of fit with the appearance of the marginal distribution. Investigate it instead as what it is - possible lack of fit. Sep 4 '14 at 3:52
• I see you edit, if your model is equivalent to a CDF then how is the line ever descending? CDF's cannot do that. I'm not saying you're not trying to model a CDF. I'm suggesting your casual equivalence to one isn't warranted so you should word your question more specifically.
– John
Sep 4 '14 at 12:23

You can generalise the logistic regression model so that the latent distribution is something other than logistic. Using the t distribution lets you capture relationships where the data are contaminated, meaning observations from the "wrong" class appear unexpectedly far away from the decision boundary -- this would be the binary equivalent of fat tails in regression with a continuous response. IOW, you model

\begin{align} y_i & = \begin{cases} 1 & \text{if } z_i > 0 \\ 0 & \text{if } z_i < 0 \end{cases}\\ z_i & = X_i\beta + \epsilon_i \end{align}

where the latent errors $\epsilon$ are distributed as $$\epsilon_i \sim t_\nu \left(0, \frac{\nu-2}{\nu} \right)$$

with $\nu > 2$ estimated from the data.

This is also called robit regression. In particular see pp.124-125 of Gelman & Hill's Data Analysis Using Regression And Multilevel/Hierarchical Models (what a mouthful); the above is eqn 6.15 in the book. You can fit it in R using package stan (which is also by Andrew Gelman and his collaborators and goes with the book), list and maybe others. You could also write code to maximise the log-likelihood directly using optim or nlminb; I haven't tried it, but it shouldn't be hard.

• Another thing Gelman suggests is trying to estimate a model like P(y=1) = 0.01 + .98 logit^-1(XB), which bounds the minimum and maximum guessed probabilities. This will mechanically put more weight in the tails but of course you can also make the .01 and .98 parameters to estimate.
– BKay
Sep 4 '14 at 18:29