Estimating the error term in a logistic regression This simulation depends on a parameter called epsilon = 2. How do I estimate epsilon given the binary outcome ($y$) and the predictors ($x1$, $x2$)? The latent variable z is not observable.
If I'm able to estimate epsilon, how do I modify my confidence intervals to account for the uncertainty caused by epsilon?
# Load required library
library(tidyverse)

nobs <- 100000

epsilon <- 2

# Generate data
set.seed(123)
x1 <- rnorm(nobs, mean = 0, sd = 1)
x2 <- rnorm(nobs, mean = 0, sd = 1)

z <- -5 + 2*x1 + -1*x2 + rnorm(nobs,mean=0,sd=epsilon)

y <- rbinom(nobs, size = 1, prob = boot::inv.logit(z))

mynewdata <- data.frame(y=y, x1=x1, x2=x2)

summary(mynewdata)

# Fit logistic regression model
fit <- glm(y ~ x1 + x2, family = "binomial", data=mynewdata)


## confidence intervals for log odds

predictions <- predict(fit, type = "link", se.fit = TRUE)

predicted_log_odds <- predictions$fit

ci_lower <- predictions$fit - 1.96 * predictions$se.fit
ci_upper <- predictions$fit + 1.96 * predictions$se.fit

results <- data.frame(pred_log_odds=predicted_log_odds,
                      lower=ci_lower,
                      upper=ci_upper)

head(results)


Writing your model mathematically helps us see what is going on.  Let's track the code step by step.  Start with the generation of the intermediate variable z conditional on the values x1 and x2:
z <- -5 + 2*x1 + -1*x2 + rnorm(nobs,mean=0,sd=epsilon)

Writing $\beta=(-5,2,-1)^\prime,$ $X = \pmatrix{\mathbf 1 & x_1 & x_2}$ for the design matrix, and $\varepsilon$ for the iid Normal disturbances (whose common standard deviation is given by the variable epsilon), that code implements the model
$$Z = X\beta + \varepsilon\tag{1}$$
for the random response $Z$ conditional on $X.$
Next,
y <- rbinom(nobs, size = 1, prob = boot::inv.logit(z))

The inv.logit function inverts the logit link, $y(z) = 1/(1 + \exp(-z)).$  The rbinom function generates a $1$ with probability prob and otherwise generates a $0.$  Thus, y represents a random variable $Y$ with
$$\Pr(Y = 1\mid Z) = \frac{1}{1 + e^{-Z}}.\tag{2}$$
Putting $(1)$ and $(2)$ together and taking expectations shows
$$E(Y = 1\mid X, \varepsilon) = \frac{1}{1 + e^{-X\beta - \varepsilon}}.$$
The regression you are trying to fit models the expectation of $Y$ conditional on $X$ alone.  To find out what this model is, we have to "marginalize" $X$ by integrating out the disturbance terms, which we can do one observation at a time because the disturbances are independent:
$$E(Y=1\mid X) = \int_{\mathbb R} \frac{1}{\sigma}\frac{\phi\left(t/\sigma\right)}{1 + e^{-X\beta - t}}\,\mathrm dt$$
where $\phi$ is the standard Normal density function.
This integral is a nonlinear function of the parameters $\beta.$  How can we tell?  If it were linear, its derivative with respect to any $\beta_j$ would be constant.  The integrand is bounded and smooth, permitting us to differentiate under the integral sign.  Here, I write $x_i$ for the values of the explanatory variables associated with observation $i:$
$$\frac{\partial}{\partial \beta_j} E(Y_i=1\mid x_i) = \frac{x_{ij}}{\sigma}\int_{\mathbb R} \frac{e^{-x_i\beta-t}}{\left(1 + e^{-x_i\beta - t}\right)^2} \phi\left(t/\sigma\right)\,\mathrm dt,$$
which is not a constant function of $\beta$ unless $x_i$ is the zero vector.  (Simple demonstration: take a sequence of values of $\beta$ for which $x_i\beta$ grows arbitrarily large.  The limit clearly is zero, but the integral just as clearly cannot be constantly zero for all $\beta.$)
Thus, you are generating data according to an unusually complicated model that is not linear in the parameters.  However, with the code
fit <- glm(y ~ x1 + x2, family = "binomial", data=mynewdata)

you are fitting a model that is linear in its parameters.  That is, you are fitting a model that differs from the one used to generate the data.
Something has to give.  Because you are using a GLM to fit the model, the applicability of its confidence intervals to the parameter $\beta$ is questionable due to the model mis-specification.  The solution is not to try to alter those intervals, but to fit a model that is appropriate for the data-generation process.
You could do that using Maximum Likelihood, for instance.  But would it be worth it?  Do you want to fit this complicated model to your data, or would it be acceptable to use a simpler model for your regression?
