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.