[R] Simulate binary data for a logistic regression Monte Carlo
Rolf Turner
r.turner at auckland.ac.nz
Tue Apr 7 22:58:36 CEST 2009
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}}
More information about the R-help
mailing list