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.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
