在 R 中提取 glmer() 函数的结果
Extracting result of glmer() function in R
我想模拟来自多级物流分布的数据 1000 次,每次估计参数并计算估计值的平均值。但似乎在 glmer()
函数中无法像 lm()
函数那样提取结果,例如 lm(y~x)$coef
。如何从 glmer()
函数中提取结果?
这是 R 代码:
#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
如何从 glmer()
函数中提取结果?
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
在混合效应模型中,您有两种类型的系数(因此 "mixed"):固定系数和随机系数。两者都可以使用专用函数从 lmer/glmer 对象中提取。例如:
lmer_obj = glmer(Y ~ X1 + X2 + (1|Subj), data = D, family = binomial)
fixef(lmer_obj) ## returns fixed effects
ranef(lmer_obj) ## returns random effects
我想模拟来自多级物流分布的数据 1000 次,每次估计参数并计算估计值的平均值。但似乎在 glmer()
函数中无法像 lm()
函数那样提取结果,例如 lm(y~x)$coef
。如何从 glmer()
函数中提取结果?
这是 R 代码:
#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
如何从 glmer()
函数中提取结果?
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
在混合效应模型中,您有两种类型的系数(因此 "mixed"):固定系数和随机系数。两者都可以使用专用函数从 lmer/glmer 对象中提取。例如:
lmer_obj = glmer(Y ~ X1 + X2 + (1|Subj), data = D, family = binomial)
fixef(lmer_obj) ## returns fixed effects
ranef(lmer_obj) ## returns random effects