On 8/04/2009, at 1:26 AM, jjh21 wrote:
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.
Well, it's the epsilon term that's doing you in. The model
that you are fitting takes no cognizance of this random effect.
You can *probably* fit a model which does take cognizance of it
using one of the variants of glm() (available in R) which include
random effects, but don't ask me how!
I don't know if it is possible to work out analytically what the
bias should be if the epsilon term is ignored when fitting the
model, but it might be. Or an approximation might be achieved.
Good luck.
cheers,
Rolf Turner
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)
######################################################################
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
______________________________________________
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.
