Heteroskedasticity underlying data generation process in logistic regression

Question

We can develop the logistic regression model using the latent variable approach:

\begin{equation} y = \begin{cases} 1, & \mathrm{if}\ X\beta + \epsilon > 0\\ 0, & \mathrm{otherwise} \end{cases} \end{equation}

where $\epsilon \sim \mathcal{L}(0,1)$. This development is equivalent to the generalized linear model formulation which makes no mention of logistic errors:

\begin{align} \begin{split} P(y = 1 \mid X) &{}= P(X\beta + \epsilon > 0 \mid X) = P(\epsilon > -X\beta \mid X) \\ &{}= P(\epsilon < X\beta \mid X) \quad \epsilon \text{ is symmetric about about zero}\\ &{}= \frac{1}{1+e^{-X\beta}} \quad \text{equivalent to GLM formulation} \end{split} \end{align}

My question is how can one obtain unbiased estimates of $\beta$ when the error is heteroskedastic, $\epsilon_i \sim \mathcal{L}(0,\sigma_i)$? I report a simple simulation to demonstrate the problem. The data generation process for the logistic regression model under heteroskedasticity was:

\begin{equation} y = \begin{cases} 1, & \mathrm{if}\ x + \epsilon > 0\\ 0, & \mathrm{otherwise} \end{cases} \quad \text{where }x\sim\mathcal{N}(0, 1)\text{ and }\epsilon\sim\mathcal{L}(0,\lvert x\rvert) \end{equation}

R syntax:

set.seed(12345)
res <- t(replicate(10000, { # 10,000 replications
  dat <- data.frame(x = rnorm(500))
  # Use Bernoulli RNG to determine homoskedasticity or heteroskedasticity
  flip <- rbinom(1, 1, .5)
  # Assign logistic error to data frame
  ifelse(flip == 0, dat$error <- rlogis(nrow(dat)),
         dat$error <- rlogis(nrow(dat), scale = abs(dat$x)))
  # Create binary outcome
  dat$y <- (rowSums(dat) > 0) + 0
  # Return coin flip and logistic regression coefficient
  c(flip, unname(coef(glm(y ~ x, binomial, dat))["x"]))
}))
# Summarize the results
par(mfrow = c(1, 2))
hist(res[res[, 1] == 0, 2], main = paste(
  "homosked, mean:", round(mean(res[res[, 1] == 0, 2]), 3)))
abline(v = mean(res[res[, 1] == 0, 2]))
hist(res[res[, 1] == 1, 2], main = paste(
  "heterosked, mean:", round(mean(res[res[, 1] == 1, 2]), 3)))
abline(v = mean(res[res[, 1] == 1, 2]))
par(mfrow = c(1, 1))

As is evident from the plot above, when the logistic error is heteroskedastic, the coefficient estimated using ML is biased downwards. Does anyone know an approach for unbiased estimation of the model coefficients under this form of heteroskedasticity? Not the specific process here, but under heteroskedastic logistic errors.

Answer 1

The model you are describing is known in econometrics as the heteroscedastic logit model and if you specify the regression relationship for the latent scale parameter, you can estimate it by maximum likelihood. Of course, this does not have to be unbiased but it will be consistent and asymptotically normal provided that you have specified everything correctly.

If you look at the literature, I suggest that you look for heteroscedastic probit rather than logit because there seems to be more material on the Gaussian case. However, clearly the same principles apply to logit. A simple package in R that allows to estimate these models is glmx with function hetglm().

Having said that: A common misunderstanding about these models is, though, that they are similar to linear regression with heteroscedastic errors. However, the heteroscedasticity is only on the latent scale. On the manifest scale, the heteroscedasticity will look like nonlinearity in the mean! Specifically, in the case you designed you simply get a step function.

Let's look at this with some more notation. I'll discuss the Gaussian case here as I'm re-using materials from presentation slides (ERCIM-2013.pdf) on this topic. We assume a latent variable $y^*$ and corresponding manifest response $y$: $$ \begin{eqnarray*} y_i^* & = & x_i^\top \beta + \varepsilon_i \\ y_i & = & I(y_i^* > 0) \end{eqnarray*} $$ For homoscedastic Gaussian errors $\varepsilon_i \sim \mathcal{N}(0, \sigma^2)$: $$ \begin{eqnarray*} \text{Prob}(y_i = 1 ~|~ x_i) & = & \text{Prob}(y_i^* > 0 ~|~ x_i) \\ %% & = & \text{Prob}(-\varepsilon_i > x_i^\top \beta ~|~ x_i) \\ & = & \Phi \left( \frac{x_i^\top \beta - 0}{\sigma} \right) \end{eqnarray*} $$ In this case the scale parameter $\sigma$ is not identified and hence taken to be $\sigma = 1$. However, while we cannot identify both $\beta$ and $\sigma$ in absolute terms, we can identify scale differences via: $$ \begin{equation*} \log(\sigma_i) ~=~ z_i^\top \gamma \end{equation*} $$ where $z_i$ must not include a constant term for identifiability. Instead of using a log-link as above to assure positivity of $\sigma_i$ you can also assume an identity link. In that case, one should really add a constant though to avoid $\sigma_i$ becoming zero: $$ \begin{equation*} \sigma_i ~=~ 1 + z_i^\top \gamma \end{equation*} $$ so that when $z_i = 0$ or $\gamma = 0$ you drop back to the standard probit model. The resulting mean function of the model is then: $$ \begin{equation*} E(y_i | x_i) ~=~ \text{Prob}(y_i = 1 ~|~ x_i) ~=~ \Phi \left( \frac{x_i^\top \beta}{\sigma_i} \right) \end{equation*} $$ Thus, on top of the nonlinearity due to the probit link function you get a rational specification instead of the usual linear predictor. So even on the link scale, this is nonlinear.

With your specification of $\eta_i = 0 + 1 \cdot x_i$ and $\sigma_i = |x_i|$ you get: $$ \Phi^{-1}(E(y_i | x_i)) ~=~ \frac{0 + x_i}{|x_i|} ~=~ \text{sign}(x_i) $$ Exactly the same applies in the logit case. Thus, your conditional mean aka the conditional "success" probability simply switches from low to high when $x_i$ changes sign. (Note that, strictly speaking, for $x_i = 0$ the above is not defined, though.) Thus, if you want to estimate that model you could also simply use glm(y ~ factor(x > 0), family = binomial).

The set of slides linked above point out a few other potential problems and pitfalls. All in all, I would recommend to only use these models if there are convincing assumptions that heteroscedasticity on a latent scale is really the issue. Otherwise I would simply use a logit GAM instead of a logit GLM with smooth nonlinear effects. Such convincing cases do exist, see e.g., our precipitation forecasting example in doi:10.1175/MWR-D-13-00271.1. But in economics and the social sciences this is more difficult, e.g., as pointed out by Alvarez & Brehm (1995, American Journal of Political Science) or Keele & Park (2006, Unpublished manuscript).

Answer 2

I'm not sure if this is what you're asking, but from the practical side, you can check for heteroscedasticity of the residuals in logistic regression models with DHARMa, and you can correct it, e.g., with glmmTMB.

The fact that you will never be 100% sure if you have the right model is inherent to any applied statistical problem. It is a problem for sure, but I don't see that the issue of heteroscedasticity is more problematic than any other misspecification problem (though it is probably less frequently checked).

EDIT - due to the comments, here an example how to diagnose the model misspecification in DHARMa

set.seed(12345)

dat <- data.frame(x = rnorm(1000))
# switch here to use one or the other option 
#dat$error <- rlogis(nrow(dat))
dat$error <- rlogis(nrow(dat), scale = abs(dat$x))
dat$y <- (rowSums(dat) > 0) + 0

fit <- glm(y ~ x, binomial, dat)

library(DHARMa)

# standard plots for 0/1  data don't show the problem 

res <- simulateResiduals(fit)
plot(res)
testDispersion(res)

# as explained at the end of the DHARMa vignette, one has to 
# group the data to see dispersion problems

qnt <- quantile(dat$x,seq(0,1,0.02))
dat$sizegroup = cut(dat$x,unique(qnt),include.lowest=TRUE)


res2 = recalculateResiduals(res , group = dat$sizegroup)
testDispersion(res2)
plot(res2)  
plotResiduals(1:nlevels(dat$sizegroup), res2$scaledResiduals)

Answer 3

The MLE was derived based on the assumption $ϵ∼L(0,1)$. Then it was used on the totally different situation, $ϵ∼L(0,|x|)$. Of course, the original MLE does not work anymore. The basic principle for MLE is: MLE is based on the distribution of the random variables; if the assumptions are changed, new MLE based on new distribution should be derived separately, and the old MLE cannot be used.

So for your situation, get $P(\epsilon < X\beta \mid X)$ under the assumption that $ϵ∼L(0,|x|)$, then derive MLE. You will get unbiased estimate.