Greetings. Here is sample code, with some comments. It shows how I can simulate data and estimate glm with binomial family when there is no individual level random error, but when I add random error into the linear predictor, I have a difficult time getting reasonable estimates of the model parameters or the variance component.
There are no clusters here, just individual level responses, so perhaps I am misunderstanding the translation from my simple mixed model to the syntax of glmmML and lmer. I get roughly the same (inaccurate) fixed effect parameter estimates from glmmML and lmer, but the information they give on the variance components is quite different. Thanks in advance. Now I paste in the example code ### Paul Johnson <[EMAIL PROTECTED]> ### 2006-03-08 N <- 1000 A <- -1 B <- 0.3 x <- 1 + 10 * rnorm(N) eta <- A + B * x pi <- exp(eta)/(1+exp(eta)) myunif <- runif(N) y <- ifelse(myunif < pi, 1, 0) plot(x,y, main=bquote( eta[i] == .(A) + .(B) * x[i] )) text ( 0.5*max(x), 0.5, expression( Prob( y[i] == 1) == frac( 1 , 1 + exp(-eta[i] )))) myglm1 <- glm ( y ~ x, family=binomial(link="logit") ) summary(myglm1) ## Just for fun.... myglm2 <- glm(y~x, family=quasibinomial) summary(myglm2) ### Mixed model: random intercept with large variance eta <- A + B * x + 5 * rnorm(N) pi <- exp(eta)/(1+exp(eta)) myunif <- runif(N) y <- ifelse(myunif < pi, 1, 0) plot(x,y, main=bquote( eta[i] == .(A) + .(B) * x[i] + u[i])) text ( 0.5*max(x), 0.5, expression( Prob( y[i] == 1) == frac( 1 , 1 + exp(-eta[i] )))) ### Parameter estimates go to hell, as expected myglm3 <- glm ( y ~ x, family=binomial(link="logit") ) summary(myglm3) ### Why doesn't quasibinomial show more evidence of the random intercept? myglm4 <- glm(y~x, family=quasibinomial) summary(myglm4) # Can I estimate with lmer? library(lme4) ### With no repeated observations, what does lmer want? ### Clusters of size 1 ? ### In lme, I'd need random= ~ 1 idnumber <- 1: length(y) mylmer1 <- lmer( y ~ x + ( 1 | idnumber ), family=binomial, method="Laplace" ) summary(mylmer1) ### Glmm wants clusters, and I don't have any clusters with more than 1 observation ### library(glmmML) myglmm1 <- glmmML(y~x, family=binomial, cluster = idnumber ) summary(myglmm1) -- Paul E. Johnson Professor, Political Science 1541 Lilac Lane, Room 504 University of Kansas ______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
