Is it possible to simulate logistic regression without randomness?

Question

We can simulate linear regression without randomness, which means we make $y=X\beta$ instead of $y=X\beta+\epsilon$. Then if we fit a linear model the coefficients will be identical to the "ground truth". Here is an example.

set.seed(0)
n    <- 1e5
p    <- 3
X    <- matrix(rnorm(n*p), ncol=p)
beta <- runif(p)
# y  <- X %*% beta + rnorm(n)*0.5
# remove the randomness
y    <- X %*% beta
dat  <- data.frame(y=y, x=X)
lm.res = lm(y ~ .-1, data=dat)
norm(as.matrix(lm.res$coefficients - beta))
[1] 2.176037e-14

My question is can we do the similar simulation with logistic regression? From this question I get the point of remove randomness can be done by using deterministic statement instead of sample from binomial distribution.

y <- ifelse(plogis(X %*% beta)>0.5,1,0) 

instead of

y <- rbinom(n,1,prob=plogis(X %*% beta))

But if we do that, complete separation will happen, and we cannot get the coefficients. On the other hand, if we add regularization, then the coefficients will not be the ones generated data.

So, what can I do to "remove the randomness in logistic regression" and solve for the exact "ground truth" coefficients like the linear regression case?

I feel I have some fundamental misunderstanding of the concept, what am I missing?

Answer 1

Logistic regression does not have an "error" term as with classical linear regression. The exception to this might be thresholded linear regression with a logistic error term, but this isn't a commonly accepted probability model which results in a logistic regression model. This is because logistic models have a mean-variance relationship. The analogue to "adding an error term" to a linear regression model is actually a quasibinomial model in which the variance is merely proportional to p*(1-p).

A related question may be how to obtain regression model results which are identical over various designs or replications. This can be done with a "trick" in regression modeling software. You can generate non-integral $Y$ outcomes from the predicted risk which result in the same logistic regression results independent of the design of $X$. For instance: x1 <- seq(-3, 3, 0.1) and x2 <- rnorm(61) as two different designs. As in your case, y1 <- plogis(0.3*x1) and y2 <- plogis(0.3*x2) both result in the same logistic regression model results with 0.3 as the log odds ratio and 0.0 as the log odds for $x=0$.

> glm(y1 ~ x1, family=binomial)

Call:  glm(formula = y1 ~ x1, family = binomial)

Coefficients:
(Intercept)           x1  
 -2.528e-16    3.000e-01  

This relates to your question because the parameter estimates are exactly as defined in your probability model, independent of the design of $x$, and without separation (e.g. log odds ratios, $\beta = \pm \infty$).

Modeling fractional results in a logistic model is an accepted way of analyzing ecological data, where the outcome may indeed be fractional. Not coincidentally, this is also a type of modeling when quasibinomial models are of the most use. Also not coincidentally, I think dispersion is proportional to a scale parameter for the logistic error term when doing "latent logistic regression".

Answer 2

I always like to think of logistic regression as what happens if you apply a binary decision to a linear model. That is, let's assume there is some underlying relationship that follows the linear model: $$ y = X\beta+\varepsilon $$ where $X$ is your independent variable and $\beta$ the coefficient (or slope) on that variable, and $\varepsilon$ is random noise. And then let's say we apply a function to the continuous variable $y$ that maps it onto a binary outcome: $$ f(y) = \left\{\begin{matrix}0, ~\operatorname{if}~ y \operatorname{\leqslant \theta} \\1, ~\operatorname{if}~ y \operatorname{>\theta} \end{matrix}\right. $$ where $\theta$ is a threshold. What is the probability that this function returns $1$, given a certain value of $X$? If we assume that $\varepsilon$ is Normally distributed with mean $0$ and variance $\sigma^2$, then we can calculate this probability as:

$$ p(f(y)=1|X)=p(y>\theta|X)=\int_\theta^\infty N\left(y; X\beta, \sigma^2 \right)dy $$

In other words, this is computing the area under the Normal distribution that is to the right of the threshold. Note that this probability is essentially what a logistic regression model tries to describe. And indeed, if you plot this probability as a function of $X$, you get something pretty close in shape to the logistic function (in fact the logistic function is often used as a convenient approximation to the cumulative Normal distribution).

For values of $X\beta$ near the threshold, the probability that $y$ will be above threshold is near $0.5$, because the noise $\varepsilon$ can sway the outcome either way. As you increase $X$, $X\beta$ will get further away from $\theta$ and $f(y)=1$ becomes more likely. Crucially, how quickly $p(f(y)=1|X)$ increases with $X$ depends on two things: the slope $\beta$ and the noise variance $\sigma^2$. More precisely, it depends on the ratio $\frac{\beta}{\sigma}$. It is this (signal-to-noise) ratio which determines the (expected) coefficient you get from a logistic regression. In other words, you can think of the coefficients in a logistic regression as controlling how much each independent variable needs to change relative to the noise in the data in order to increase the probability of a certain outcome by some amount.

Now to come to your question: you're asking if it is possible to eliminate all randomness, i.e. to have no noise. This would mean that $\sigma$ equals $0$, and therefore $\frac{\beta}{\sigma}$ would be undefined (or "infinite"). This explains what you found, that you cannot estimate the coefficients when there is no noise. Indeed, you can think of the perfect separation you achieve without noise as corresponding to an infinite coefficient on your independent variable, since (for $X\beta$ near the threshold $\theta$) you only need to change $X$ an infinitesimal amount in order to go all the way from $p(y>\theta|X)=0$ to $p(y>\theta|X)=1$.

Edit: actually one thing you could do is instead of drawing samples from a binomial distribution to simulate your data, replace these samples by their expectation, i.e. the probability predicted by the simulated logistic function. That way, you're removing the randomness that derives from simulating a limited sample (i.e. the sampling variability), and thus your coefficient estimates should then equal the ground truth (since there is one logistic function that exactly fits these values).