46

我有一个具有固定和随机效果的 mer 对象。如何提取随机效应的方差估计？这是我的问题的简化版本。

study <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy)
study

这会产生很长的输出 - 在这种情况下不会太长。无论如何，我如何明确选择

Random effects:
Groups   Name        Variance Std.Dev.
Subject  (Intercept) 1378.18  37.124  
Residual              960.46  30.991

输出的一部分？我想要价值观本身。

我看了很久

str(study)

那里什么都没有！还检查了 lme4 包中的任何提取器功能，但无济于事。请帮忙！

4

7 回答 7

81

其他一些答案是可行的，但我声称最好的答案是使用为此设计的访问器方法—— VarCorr（这与lme4的前身nlme包中的相同）。

最新版本lme4（版本 1.1-7，但以下所有内容可能适用于 >= 1.0 版本）中的更新，VarCorr比以前更灵活，并且应该做你想做的一切，而无需在拟合的模型对象内四处寻找。

library(lme4)
study <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy)
VarCorr(study)
##  Groups   Name        Std.Dev.
##  Subject  (Intercept) 37.124  
##  Residual             30.991

默认情况下VarCorr()打印标准偏差，但如果您愿意，您可以获取方差：

print(VarCorr(study),comp="Variance")
##  Groups   Name        Variance
##  Subject  (Intercept) 1378.18 
##  Residual              960.46

（comp=c("Variance","Std.Dev.")将打印两者）。

为了获得更大的灵活性，您可以使用该as.data.frame方法转换VarCorr对象，该方法给出了分组变量、效果变量以及方差/协方差或标准差/相关性：

as.data.frame(VarCorr(study))
##        grp        var1 var2      vcov    sdcor
## 1  Subject (Intercept) <NA> 1378.1785 37.12383
## 2 Residual        <NA> <NA>  960.4566 30.99123

最后，VarCorr对象的原始形式（如果您不必这样做，您可能不应该与您混淆）是方差-协方差矩阵的列表，其中包含编码标准偏差和相关性的附加（冗余）信息，以及属性 ( "sc") 给出残差标准偏差并指定模型是否具有估计的尺度参数 ( "useSc")。

unclass(VarCorr(fm1))
## $Subject
##             (Intercept)      Days
## (Intercept)  612.089748  9.604335
## Days           9.604335 35.071662
## attr(,"stddev")
## (Intercept)        Days 
##   24.740448    5.922133 
## attr(,"correlation")
##             (Intercept)       Days
## (Intercept)  1.00000000 0.06555134
## Days         0.06555134 1.00000000
## 
## attr(,"sc")
## [1] 25.59182
## attr(,"useSc")
## [1] TRUE
##
于 2011-12-15T23:46:52.610 回答
14

lmer返回一个 S4 对象，所以这应该工作：

remat <- summary(study)@REmat
print(remat, quote=FALSE)

哪个打印：

 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 1378.18  37.124  
 Residual              960.46  30.991

...通常，您可以查看“mer”对象的print和summary方法的来源：

class(study) # mer
selectMethod("print", "mer")
selectMethod("summary", "mer")
于 2011-12-15T21:42:12.670 回答
5 
> attributes(summary(study))$REmat
 Groups     Name          Variance  Std.Dev.
 "Subject"  "(Intercept)" "1378.18" "37.124"
 "Residual" ""            " 960.46" "30.991"
于 2011-12-15T21:41:38.583 回答
1

这个答案很大程度上基于@Ben Bolker 的答案，但是如果人们对此不熟悉并且想要自己的值，而不仅仅是值的打印输出（正如 OP 似乎想要的那样），那么您可以按如下方式提取值：

将VarCorr对象转换为数据框。

re_dat = as.data.frame(VarCorr(study))

然后访问每个单独的值：

int_vcov = re_dat[1,'vcov']
resid_vcov = re_dat[2,'vcov']

使用此方法（在您创建的日期框架中指定行和列），您可以访问您想要的任何值。

于 2019-03-20T16:38:59.573 回答
1

另一种可能是

sum <- summary (study)
var <- data.frame (sum$varcor)
于 2021-06-26T21:51:46.693 回答
0

这个包对这样的事情很有用（https://easystats.github.io/insight/reference/index.html）

library("insight")

get_variance_random(study) #Where study is your fit mixed model
于 2020-04-26T15:46:16.467 回答
-7

尝试

attributes(study)

举个例子：

> women
   height weight
1      58    115
2      59    117
3      60    120
4      61    123
5      62    126
6      63    129
7      64    132
8      65    135
9      66    139
10     67    142
11     68    146
12     69    150
13     70    154
14     71    159
15     72    164

> lm1 <- lm(height ~ weight, data=women)
> attributes(lm1)
$names
 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "xlevels"       "call"          "terms"         "model"        

$class
[1] "lm"

> lm1$coefficients
(Intercept)      weight 
 25.7234557   0.2872492 

> lm1$coefficients[[1]]

[1] 25.72346


> lm1$coefficients[[2]]

[1] 0.2872492

结束。

于 2011-12-15T21:38:52.660 回答