# Logistic regression simulation in order to show that intercept is biased when Y=1 is rare [duplicate]

I'm trying to simulate a logistic regression. My goal is showing that if Y=1 is rare, than the intercept is biased. In my R script I define the logistic regression model through the latent variable's approach (see for example pp. 140 http://gking.harvard.edu/files/abs/0s-abs.shtml):

x   <- rnorm(10000)

b0h <- numeric(1000)
b1h <- numeric(1000)

for(i in 1:1000){
eps <- rlogis(10000)
eta <- 1+2*x+eps
y   <-numeric(10000)
y   <- ifelse (eta>0,1,0)

m      <- glm(y~x,family=binomial)
b0h[i] <- coef(m)
b1h[i] <- coef(m)
}

mean(b0h)
mean(b1h)
hist(b0h)
hist(b1h)


The problem here is that I don't know how to force the observations y to be balanced before (50:50), then unbalanced (90:10). As we can see with the function table(), in my script the proportion of ones is random.

table(y)


How to solve this problem?

• If you are truly simulating from a logistic regression model then the proportion has to be random. – Xi'an Dec 9 '14 at 10:57
• Thanks for your answer @Xi'an. I know that the proportion has to be random, but I just want to show that if the class is skewed (one class is rare), than if I use logistic regression the intercept is biased (this is what I think to have understood from gking.harvard.edu/files/abs/0s-abs.shtml. How to force the class to be skewed in my Rscript? – Luca Dibo Dec 9 '14 at 11:15
• Then you have to sample from another model. – Xi'an Dec 9 '14 at 11:46
• – gung - Reinstate Monica Dec 9 '14 at 19:23
• You're never going to show this with a simulation! All you could possibly discover is that a certain amount of bias persists in the largest simulations your computer can handle. I am having a hard time understanding why the intercept should be biased asymptotically--you actually contradict yourself with this statement (since "asymptotically unbiased" implies the bias becomes unmeasurable with sufficiently large samples)--so I would guess it might have something to do with the nature of the asymptotics. – whuber Dec 9 '14 at 23:02

Logistic regression doesn't really have an error term. Alternatively, you can think of the response distribution (the binomial) as having its random component intrinsically 'built-in' (for more, it may help to read my answer here: Difference between logit and probit models). As a result, I think it is conceptually clearer to generate data for simulations directly from a binomial parameterized as the logistic transformation of the structural component of the model, rather than use the logistic as a sort of error term.

From there, if you want the long run probability that $Y = 1$ to be $.5$ or $.1$, you just need your structural component to be balanced around $0$ (for $.5$), or $-2.197225$ (for $.1$). I got those values by converting the response probability to the log odds:
$$\log(\text{odds}(Y=1)) = \frac{\exp(Pr(Y = 1))}{(1+\exp(Pr(Y = 1))}$$ The most convenient way to do this will be to use those values for your intercept ($\beta_0$) and have your slope be $0$. (Alternatively, you can use any two parameter values, $\beta_0$ and $\beta_1$, that you like such that, given your $X$ values, the log mean odds equals, e.g., $-2.197225$.) Here is an example with R code:

lo2p = function(lo){      # this function will perform the logistic transformation
odds = exp(lo)          #   of the structural component of the data generating
p    = odds / (1+odds)  #   process
return(p)
}

N     = 1000              # these are the true values of the DGP
beta0 = -2.197225
beta1 = 0

set.seed(8361)            # this makes the simulation exactly reproducible
x     = rnorm(N)
lo    = beta0 + beta1*x
p     = lo2p(lo)          # these will be the parameters of the binomial response

b0h   = vector(length=N)  # these will store the results
b1h   = vector(length=N)
y1prp = vector(length=N)  # (for the proportion of y=1)

for(i in 1:1000){         # here is the simulation
y        = rbinom(n=N, size=1, prob=p)
m        = glm(y~x, family=binomial)
b0h[i]   = coef(m)
b1h[i]   = coef(m)
y1prp[i] = mean(y)
}

mean(b0h)                 # these are the results
#  -2.205844
mean(b1h)
#  -0.0003422177
mean(y1prp)
#  0.100036
hist(b0h)
hist(b1h)  • You don't want to control the mean log odds to achieve a given proportion of successes; you want to control the (weighted) mean probability. The two means are not quite equivalent. – whuber Dec 9 '14 at 19:32
• @whuber, hmmm I think I see your point: because of the intervening non-linear transformation, the mean of p(Y=1) for a large sample of simulated datasets could diverge from lo(Y=1) because of the asymmetry in the transformations of the datasets in the tails. Is that what you're getting at? (It does seem to work awfully well w/ lo(Y=1)=.1, but maybe that's not extreme enough.) – gung - Reinstate Monica Dec 9 '14 at 19:44
• Yes, that's what I mean. Because the logistic transformation is nearly linear for values close to zero (around $-1$ to $1$), the difference really shows up only for more extreme values--precisely the case you are discussing in your answer! For instance, let $p_1=0.01$ and $p_2=0.19$. Although their mean is $0.1$, their mean log odds corresponds to a probability of only $0.046$. Rare events have more "leverage" in the averaging when it is done in terms of log odds. – whuber Dec 9 '14 at 19:55
• @whuber, do you think I should delete this answer? – gung - Reinstate Monica Dec 9 '14 at 19:59
• @whuber, I changed that phrasing, let me know if you think it needs more. BTW, I am aware of replicate, but I don't tend to use such here; my goal is to make my R code as simple & transparent as possible, even if clunky, in the hopes that even those who don't know R could get the gist. – gung - Reinstate Monica Dec 9 '14 at 21:24