# Heteroskedasticity underlying data generation process in logistic regression

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.

• @DimitriyV.Masterov thanks for the link. I came across the blog post at some point. I went into a little detail in my question because there are people who are only familiar with the GLM formulation for the logistic regression model. Oct 10, 2018 at 2:57
• If you really only care about $\beta$, note that it is possible to estimate this up to scale without assuming that $e$ and $x$ are independent, and without assuming a distribution on $e$ using the maximum score estimator - in fact even the variance $Var(e|x)$ is unrestricted. What you need to assume is that the median of $e$ and $x$ is 0 and that at least of on the estimated parameters is not zero and continuous Oct 10, 2018 at 6:40
• @Repmat Although I have accepted an answer, can you develop this further in an answer below? I or anyone who comes across this question could benefit from it. Oct 10, 2018 at 15:41

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).

• Thank you very much for your reply. In the past day, I've come across these models and found the glmx package. I've started testing its behaviour under simulation. I was also able to replicate it using standard MLE in R - I originally kept the intercept in the equation for the scale, that didn't go well. Thanks again for the explanation and the links you provided. Oct 9, 2018 at 22:28
• Models with identity link in the scale equation are always a bit tricky. In the cases we considered the log link either worked better or was on par with identity or square root link. But the main problem is really to check what kind of mean function you get on the manifest scale from the heteroscedastic function on the latent scale. It's easy to specify something that collapses or is very close to collapsing to something simpler... Oct 9, 2018 at 23:08
• I found the manuscript by Keele & Park very informative! Thanks again. Oct 10, 2018 at 4:01

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)

• I'm a DHARMa user and had not considered this. I'll see if it can pick this up in the data and report back, thanks. And my interest is not related to any application, more theoretical. One thing I'm almost sure of is that glmmTMB will not let you model dispersion if your variable is Bernoulli because of the standard relationship between mean and variance for a Bernoulli variable. Oct 9, 2018 at 11:57
• Jim, you can adjust the dispersion in glmmTMB with the dispersion parameter, as long as you take a distribution with variable dispersion. For a binomial, the betabinomial should work ... for anything more complicated, one would probably have to move to a general framework, e.g. Jags or Stan. Oct 9, 2018 at 17:01
• I tried out DHARMa on data with a heteroskedastic data generating process for the logistic regression model; the resulting plots appeared adequate. I do not think it's a type of model misspecification that it can pick up. Oct 10, 2018 at 3:59
• The issue with 0/1 data is that one always has to group them to see dispersion problems. I have added an example in my answer. I guess Achim's comments are more to the point regarding this specific example, but in general I think the workflow described above works well for logistic regression. Oct 10, 2018 at 9:54
• Thanks! Tested it for false positive and false negative and it worked about it. Oct 10, 2018 at 13:50

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.

• In empirical data analysis, one has no idea of the data generation process. So it is no use deriving an MLE for this specific process. I'm more interested in general approaches. It is acceptable if there are none. I see how my last sentence could have been misleading in that regard. Oct 9, 2018 at 4:49
• @HeteroskedasticJim thats why god invented robust standard errors, under the right assumptions this will give you the (nominally) right inference Oct 9, 2018 at 5:28
• @Repmat lol, but my question is about bias not inference. Oct 9, 2018 at 5:42
• @Repmat So-called robust standard errors only help you if you have the estimating function correctly specified and the rest of the likelihood misspecified. Here, this does not help: The mean function of the manifest variable is misspecified! Hence coefficient estimates are inconsistent and no standard errors can save you save you from that... Oct 9, 2018 at 23:48
• @AchimZeileis There are certain situations why it is still usefull. Hence my "under the rightassumptions" - of course I should have qualified this (its not a general case in any sense). Regardless I think it largely misses the point, there is not particular reason to care about the estimate of $\beta$. I have yet to see a single application where $E(y^* \mid x)$, i.e. the effect on the latent response, was of any interest. The interest is almost always on $P(Y=1 \mid x)$, and probit/logit might still do a pretty good job here even with heteroskedasticity Oct 10, 2018 at 6:24