[R] Simulate binary data for a logistic regression Monte Carlo

[jjh21]

Hello,

I am trying to simulate binary outcome data for a logistic regression Monte
Carlo study. I need to eventually be able to manipulate the structure of the
error term to give groups of observations a random effect. Right now I am
just doing a very basic set up to make sure I can recover the parameters
properly. I am running into trouble with the code below. It works if you
take out the object 'epsilon,' but I need that object because it is where I
plan to add in the random effects eventually. Any suggestions? Thanks.

set.seed(1)
alpha <- numeric(100) # Vector to store coefficient estimates
X1 <- numeric(100) 
X2 <- numeric(100)

a <- 0 # True parameters
B1 <- .5
B2 <- .85

for (i in 1:100){
 x1 <- rnorm(1000)
 x2 <- rnorm(1000)
 epsilon <- rnorm(1000) 
  
 LP <- a + B1*x1 + B2*x2 + epsilon # Linear predictor
 y <- rbinom(length(LP), 1, plogis(LP)) # Binomial link
 
 model <- glm(y ~ x1 + x2, family = binomial)
 alpha[i] <- coef(model)[1]
 X1[i] <- coef(model)[2]
 X2[i] <- coef(model)[3]
}

mean(X1) # Shows a bias downward (should be about .5)
mean(X2) # Shows a bias downward (should be about .85)
mean(alpha) # Pretty close (should be about 0)
-- 
View this message in context: http://www.nabble.com/Simulate-binary-data-for-a-logistic-regression-Monte-Carlo-tp22928872p22928872.html
Sent from the R help mailing list archive at Nabble.com.