I want to simulate data from Multilevel Logistic Distribution 1000 time and each time estimating parameter and compute the average of the estimates . But it seems , in glmer() function results cannot be extracted like lm() function , say , lm(y~x)$coef . How can I extract the results from glmer() function ?
Here is the R code :
#Simulating data from multilevel logistic distribution
library(mvtnorm)
library(lme4)
set.seed(1234)
## J = number of groups
## n = group size
## g00,g10,g01,g11 = fixed effect parameters .
## s2_0,s2_1,s01 = variance values for the group level random effect .
simu <- function(J,n,g00,g10,g01,g11,s2_0,s2_1,s01){
n_j <- rep(n,J) ## number of individuals in jth group
N <- sum(n_j) ## sample size
#Simulate the covariate value for this sample size .
z <- rnorm(J)
x <- rnorm(N)
#Generate (u_0j,u_1j) from a bivariate normal .
mu <- c(0,0)
sig <- matrix(c(s2_0,s01,s01,s2_1),ncol=2)
u <- rmvnorm(J,mean=mu,sigma=sig,method="chol")
#Now form the linear predictor .
pi_0 <- g00 +g01*z + u[,1]
pi_1 <- g10 + g11*z + u[,2]
eta <- rep(pi_0,n_j)+rep(pi_1,n_j)*x
#Transform back to the probability scale .
p <- exp(eta)/(1+exp(eta))
#Simulate a bernoulli from each individual distribution .
y <- rbinom(N,1,p)
# estimating parameters
sim_data_mat <- matrix(c(y,x,rep(z,n_j),rep(1:30,n_j)),ncol=4)
sim_data <- data.frame(sim_data_mat)
colnames(sim_data) <- c("Y","X","Z","cluster")
res <-summary(glmer(Y~X*Z+(1|cluster),family=binomial,data=sim_data,nAGQ=10))
res$coef
}
out <- replicate(10,simu(30,5,-1,.3,.3,.3,.13,1,0))
##Error in res$coef : $ operator not defined for this S4 class
How can I extract the results from glmer() function ?
sessionInfo()
R version 2.14.0 (2011-10-31)
Platform: i386-pc-mingw32/i386 (32-bit)
locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4_0.999999-0 Matrix_1.0-1 lattice_0.20-10 mvtnorm_0.9-9994
loaded via a namespace (and not attached):
[1] grid_2.14.0 nlme_3.1-102 stats4_2.14.0 tools_2.14.0
In mixed-effects models you have two types of coefficients (hence "mixed"): fixed and random. Both can be extracted from lmer/glmer objects using the dedicated functions. For example:
lmer_obj = glmer(Y ~ X1 + X2 + (1|Subj), data = D, family = binomial)
fixef(lmer_obj) ## returns fixed effects
ranef(lmer_obj) ## returns random effects