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

[Rolf Turner]

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}}