lmer 输出的方差分量的标准误差
Standard Error of variance component from the output of lmer
我需要从 lmer 的输出中提取 standard error 方差分量。
library(lme4)
model <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)
以下产生方差分量的估计：
s2 <- VarCorr(model)$Subject[1]
它是不是方差的标准误差。我想要标准错误。我怎样才能拥有它？
编辑：
也许我无法让你理解我所说的“方差分量的标准误差”是什么意思。所以我正在编辑我的 post .
在道格拉斯·C·蒙哥马利 (Douglas C. Montgomery) 着的 Design and Analysis of Experiments 一书中的第 12 章随机因素实验中，在本章末尾，示例 12-2 由 SAS 完成。在示例 12-2 中，模型是双因子因子随机效应模型。输出在 Table 12-17
中给出
我正在尝试通过 lmer 在 R 中拟合模型。
library(lme4)
fit <- lmer(y~(1|operator)+(1|part),data=dat)
用于提取 Estimate 的 R 代码，在 table 12-17 中由 4 注释：
est_ope=VarCorr(fit)$operator[1]
est_part = VarCorr(fit)$part[1]
sig = summary(fit)$sigma
est_res = sig^2
现在我想提取 Std Errors 的结果，由 lmer 输出中的 table 12-17 中的 5 注释。
非常感谢！
我不太确定你所说的 "standard error of variance component" 是什么意思。我最好的猜测（基于您的代码）是您想要随机效应的标准误差。您可以使用软件包 arm 获得它：
library(arm)
se.ranef(model)
#$Subject
# (Intercept)
#308 9.475668
#309 9.475668
#310 9.475668
#330 9.475668
#331 9.475668
#332 9.475668
#333 9.475668
#334 9.475668
#335 9.475668
#337 9.475668
#349 9.475668
#350 9.475668
#351 9.475668
#352 9.475668
#369 9.475668
#370 9.475668
#371 9.475668
#372 9.475668
这其实就是随机效应的条件方差-协方差矩阵的平方根：
sqrt(attr(ranef(model, condVar = TRUE)$Subject, "postVar"))
我认为您正在寻找方差估计值的 Wald 标准误差。请注意，这些（正如 Doug Bates 经常指出的那样）Wald 标准误差通常 非常差 对方差不确定性的估计，因为似然曲线通常远非方差的二次方规模......我假设你知道你在做什么并且对这些数字有一些很好的用途......
这可以（现在）使用 merDeriv 包完成。
library(lme4)
library(merDeriv)
m1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
sqrt(diag(vcov(m1, full = TRUE)))
vv <- vcov(m1, full = TRUE)
colnames(vv)
## [1] "(Intercept)" "Days"
## [3] "cov_Subject.(Intercept)" "cov_Subject.Days.(Intercept)"
## [5] "cov_Subject.Days" "residual"
如果我们想要方差分量的标准误差，我们取对角线的平方根并只保留最后三个元素：
sqrt(diag(vv)[3:5])
## [1] 288.78602 46.67876 14.78208
旧答案
library("lme4")
model <- lmer(Reaction ~ Days + (1|Subject), sleepstudy, REML=FALSE)
（目前对于 REML 估计要做到这一点相当困难......）
提取根据标准偏差和相关性而不是 Cholesky 因子参数化的偏差函数（注意这是一个内部函数，因此不能保证它在未来会以相同的方式继续工作...... )
dd.ML <- lme4:::devfun2(model,useSc=TRUE,signames=FALSE)
提取参数作为原始尺度的标准差：
vv <- as.data.frame(VarCorr(model)) ## need ML estimates!
pars <- vv[,"sdcor"]
## will need to be careful about order if using this for
## a random-slopes model ...
现在计算二阶导数（Hessian）矩阵：
library("numDeriv")
hh1 <- hessian(dd.ML,pars)
vv2 <- 2*solve(hh1) ## 2* converts from log-likelihood to deviance scale
sqrt(diag(vv2)) ## get standard errors
这些是标准差的标准误差：将它们加倍得到方差的标准误差（当你转换一个值时，它的标准误差根据转换的导数缩放）。
我认为这应该可以，但您可能需要仔细检查一下...
mn2=lmer(pun~ pre + (pre|pro), REML = TRUE, data = pro)
summary(mn2)
coe2=coef(mn2)
coe2
# Matriz de varianza-covarianza (covarianza)
as.data.frame(VarCorr(mn2))
# Extraer coeficientes fijos
fixef(mn2)
# Extraer desvios de a - alfa y b - beta
re=as.data.frame(ranef(mn2))
我需要从 lmer 的输出中提取 standard error 方差分量。
library(lme4)
model <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)
以下产生方差分量的估计：
s2 <- VarCorr(model)$Subject[1]
它是不是方差的标准误差。我想要标准错误。我怎样才能拥有它？
编辑：
也许我无法让你理解我所说的“方差分量的标准误差”是什么意思。所以我正在编辑我的 post .
在道格拉斯·C·蒙哥马利 (Douglas C. Montgomery) 着的 Design and Analysis of Experiments 一书中的第 12 章随机因素实验中，在本章末尾，示例 12-2 由 SAS 完成。在示例 12-2 中，模型是双因子因子随机效应模型。输出在 Table 12-17
中给出我正在尝试通过 lmer 在 R 中拟合模型。
library(lme4)
fit <- lmer(y~(1|operator)+(1|part),data=dat)
用于提取 Estimate 的 R 代码，在 table 12-17 中由 4 注释：
est_ope=VarCorr(fit)$operator[1]
est_part = VarCorr(fit)$part[1]
sig = summary(fit)$sigma
est_res = sig^2
现在我想提取 Std Errors 的结果，由 lmer 输出中的 table 12-17 中的 5 注释。
非常感谢！
我不太确定你所说的 "standard error of variance component" 是什么意思。我最好的猜测（基于您的代码）是您想要随机效应的标准误差。您可以使用软件包 arm 获得它：
library(arm)
se.ranef(model)
#$Subject
# (Intercept)
#308 9.475668
#309 9.475668
#310 9.475668
#330 9.475668
#331 9.475668
#332 9.475668
#333 9.475668
#334 9.475668
#335 9.475668
#337 9.475668
#349 9.475668
#350 9.475668
#351 9.475668
#352 9.475668
#369 9.475668
#370 9.475668
#371 9.475668
#372 9.475668
这其实就是随机效应的条件方差-协方差矩阵的平方根：
sqrt(attr(ranef(model, condVar = TRUE)$Subject, "postVar"))
我认为您正在寻找方差估计值的 Wald 标准误差。请注意，这些（正如 Doug Bates 经常指出的那样）Wald 标准误差通常 非常差 对方差不确定性的估计，因为似然曲线通常远非方差的二次方规模......我假设你知道你在做什么并且对这些数字有一些很好的用途......
这可以（现在）使用 merDeriv 包完成。
library(lme4)
library(merDeriv)
m1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
sqrt(diag(vcov(m1, full = TRUE)))
vv <- vcov(m1, full = TRUE)
colnames(vv)
## [1] "(Intercept)" "Days"
## [3] "cov_Subject.(Intercept)" "cov_Subject.Days.(Intercept)"
## [5] "cov_Subject.Days" "residual"
如果我们想要方差分量的标准误差，我们取对角线的平方根并只保留最后三个元素：
sqrt(diag(vv)[3:5])
## [1] 288.78602 46.67876 14.78208
旧答案
library("lme4")
model <- lmer(Reaction ~ Days + (1|Subject), sleepstudy, REML=FALSE)
（目前对于 REML 估计要做到这一点相当困难......）
提取根据标准偏差和相关性而不是 Cholesky 因子参数化的偏差函数（注意这是一个内部函数，因此不能保证它在未来会以相同的方式继续工作...... )
dd.ML <- lme4:::devfun2(model,useSc=TRUE,signames=FALSE)
提取参数作为原始尺度的标准差：
vv <- as.data.frame(VarCorr(model)) ## need ML estimates!
pars <- vv[,"sdcor"]
## will need to be careful about order if using this for
## a random-slopes model ...
现在计算二阶导数（Hessian）矩阵：
library("numDeriv")
hh1 <- hessian(dd.ML,pars)
vv2 <- 2*solve(hh1) ## 2* converts from log-likelihood to deviance scale
sqrt(diag(vv2)) ## get standard errors
这些是标准差的标准误差：将它们加倍得到方差的标准误差（当你转换一个值时，它的标准误差根据转换的导数缩放）。
我认为这应该可以，但您可能需要仔细检查一下...
mn2=lmer(pun~ pre + (pre|pro), REML = TRUE, data = pro)
summary(mn2)
coe2=coef(mn2)
coe2
# Matriz de varianza-covarianza (covarianza)
as.data.frame(VarCorr(mn2))
# Extraer coeficientes fijos
fixef(mn2)
# Extraer desvios de a - alfa y b - beta
re=as.data.frame(ranef(mn2))