logo标签列表

从`lmermod`对象中提取结果

rlme4
7
我想从一个 lmerMod 对象中提取一些结果。
require(lme4)
(fm1 <- lmer(Yield ~ 1|Batch, Dyestuff))

这会产生
Linear mixed model fit by REML ['lmerMod']
Formula: Yield ~ 1 | Batch
   Data: Dyestuff
REML criterion at convergence: 319.6543
Random effects:
 Groups   Name        Std.Dev.
 Batch    (Intercept) 42.00   
 Residual             49.51   
Number of obs: 30, groups:  Batch, 6
Fixed Effects:
(Intercept)  
       1527

特别是，我想提取出两个随机效应的标准差，即42.00和49.51。我猜可能有一种内置方法可以做到这一点，但我没有快速找到它，所以我认为通过使用str(fm1)应该能够找到它，它会产生以下结果：
Formal class 'lmerMod' [package "lme4"] with 13 slots
  ..@ resp   :Reference class 'lmerResp' [package "lme4"] with 9 fields
  .. ..$ Ptr    :<externalptr> 
  .. ..$ mu     : num [1:30] 1510 1510 1510 1510 1510 ...
  .. ..$ offset : num [1:30] 0 0 0 0 0 0 0 0 0 0 ...
  .. ..$ sqrtXwt: num [1:30] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ sqrtrwt: num [1:30] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ weights: num [1:30] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ wtres  : num [1:30] 35.1 -69.9 -69.9 10.1 70.1 ...
  .. ..$ y      : num [1:30] 1545 1440 1440 1520 1580 ...
  .. ..$ REML   : int 1
  .. ..and 26 methods, of which 14 are possibly relevant:
  .. ..  allInfo, copy#envRefClass, initialize, initialize#lmResp,
  .. ..  initializePtr, initializePtr#lmResp, objective, ptr, ptr#lmResp,
  .. ..  setOffset, setResp, setWeights, updateMu, wrss
  ..@ Gp     : int [1:2] 0 6
  ..@ call   : language lmer(formula = Yield ~ 1 | Batch, data = Dyestuff)
  ..@ frame  :'data.frame':     30 obs. of  2 variables:
  .. ..$ Yield: num [1:30] 1545 1440 1440 1520 1580 ...
  .. ..$ Batch: Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 2 2 2 2 2 ...
  .. ..- attr(*, "terms")=Classes 'terms', 'formula' length 3 Yield ~ 1 + Batch
  .. .. .. ..- attr(*, "variables")= language list(Yield, Batch)
  .. .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. ..$ : chr [1:2] "Yield" "Batch"
  .. .. .. .. .. ..$ : chr "Batch"
  .. .. .. ..- attr(*, "term.labels")= chr "Batch"
  .. .. .. ..- attr(*, "order")= int 1
  .. .. .. ..- attr(*, "intercept")= int 1
  .. .. .. ..- attr(*, "response")= int 1
  .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. .. ..- attr(*, "predvars")= language list(Yield, Batch)
  .. .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
  .. .. .. .. ..- attr(*, "names")= chr [1:2] "Yield" "Batch"
  .. .. .. ..- attr(*, "predvars.fixed")= language list(Yield)
  .. ..- attr(*, "formula")=Class 'formula' length 3 Yield ~ 1 | Batch
  .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  ..@ flist  :List of 1
  .. ..$ Batch: Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 2 2 2 2 2 ...
  .. ..- attr(*, "assign")= int 1
  ..@ cnms   :List of 1
  .. ..$ Batch: chr "(Intercept)"
  ..@ lower  : num 0
  ..@ theta  : num 0.848
  ..@ beta   : num 1527
  ..@ u      : num [1:6] -20.755 0.461 33.669 -27.212 66.877 ...
  ..@ devcomp:List of 2
  .. ..$ cmp : Named num [1:10] 9.15 1.88 61495.41 9590.84 71086.25 ...
  .. .. ..- attr(*, "names")= chr [1:10] "ldL2" "ldRX2" "wrss" "ussq" ...
  .. ..$ dims: Named int [1:12] 30 30 1 29 1 6 1 1 0 1 ...
  .. .. ..- attr(*, "names")= chr [1:12] "N" "n" "p" "nmp" ...
  ..@ pp     :Reference class 'merPredD' [package "lme4"] with 18 fields
  .. ..$ Lambdat:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. ..@ i       : int [1:6] 0 1 2 3 4 5
  .. .. .. ..@ p       : int [1:7] 0 1 2 3 4 5 6
  .. .. .. ..@ Dim     : int [1:2] 6 6
  .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. ..$ : NULL
  .. .. .. .. ..$ : NULL
  .. .. .. ..@ x       : num [1:6] 0.848 0.848 0.848 0.848 0.848 ...
  .. .. .. ..@ factors : list()
  .. ..$ LamtUt :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. ..@ i       : int [1:30] 0 0 0 0 0 1 1 1 1 1 ...
  .. .. .. ..@ p       : int [1:31] 0 1 2 3 4 5 6 7 8 9 ...
  .. .. .. ..@ Dim     : int [1:2] 6 30
  .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. ..$ : NULL
  .. .. .. .. ..$ : NULL
  .. .. .. ..@ x       : num [1:30] 0.848 0.848 0.848 0.848 0.848 ...
  .. .. .. ..@ factors : list()
  .. ..$ Lind   : int [1:6] 1 1 1 1 1 1
  .. ..$ Ptr    :<externalptr> 
  .. ..$ RZX    : num [1:6, 1] 1.98 1.98 1.98 1.98 1.98 ...
  .. ..$ Ut     :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. ..@ i       : int [1:30] 0 0 0 0 0 1 1 1 1 1 ...
  .. .. .. ..@ p       : int [1:31] 0 1 2 3 4 5 6 7 8 9 ...
  .. .. .. ..@ Dim     : int [1:2] 6 30
  .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. ..$ : chr [1:6] "A" "B" "C" "D" ...
  .. .. .. .. ..$ : NULL
  .. .. .. ..@ x       : num [1:30] 1 1 1 1 1 1 1 1 1 1 ...
  .. .. .. ..@ factors : list()
  .. ..$ Utr    : num [1:6] 6384 6481 6634 6354 6787 ...
  .. ..$ V      : num [1:30, 1] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ VtV    : num [1, 1] 30
  .. ..$ Vtr    : num 45825
  .. ..$ X      : num [1:30, 1] 1 1 1 1 1 1 1 1 1 1 ...
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:30] "1" "2" "3" "4" ...
  .. .. .. ..$ : chr "(Intercept)"
  .. .. ..- attr(*, "assign")= int 0
  .. ..$ Xwts   : num [1:30] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ Zt     :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. ..@ i       : int [1:30] 0 0 0 0 0 1 1 1 1 1 ...
  .. .. .. ..@ p       : int [1:31] 0 1 2 3 4 5 6 7 8 9 ...
  .. .. .. ..@ Dim     : int [1:2] 6 30
  .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. ..$ : chr [1:6] "A" "B" "C" "D" ...
  .. .. .. .. ..$ : NULL
  .. .. .. ..@ x       : num [1:30] 1 1 1 1 1 1 1 1 1 1 ...
  .. .. .. ..@ factors : list()
  .. ..$ beta0  : num 0
  .. ..$ delb   : num 1527
  .. ..$ delu   : num [1:6] -20.755 0.461 33.669 -27.212 66.877 ...
  .. ..$ theta  : num 0.848
  .. ..$ u0     : num [1:6] 0 0 0 0 0 0
  .. ..and 42 methods, of which 30 are possibly relevant:
  .. ..  b, beta, CcNumer, copy#envRefClass, initialize, initializePtr,
  .. ..  installPars, L, ldL2, ldRX2, linPred, P, ptr, RX, RXdiag, RXi,
  .. ..  setBeta0, setDelb, setDelu, setTheta, solve, solveU, sqrL, u, unsc,
  .. ..  updateDecomp, updateL, updateLamtUt, updateRes, updateXwts
  ..@ optinfo:List of 7
  .. ..$ optimizer: chr "bobyqa"
  .. ..$ control  :List of 1
  .. .. ..$ iprint: int 0
  .. ..$ derivs   :List of 2
  .. .. ..$ gradient: num 1.61e-07
  .. .. ..$ Hessian : num [1, 1] 14.1
  .. ..$ conv     :List of 2
  .. .. ..$ opt : int 0
  .. .. ..$ lme4: list()
  .. ..$ feval    : int 16
  .. ..$ warnings : list()
  .. ..$ val      : num 0.848

...但这并没有让我更明白。
你尝试过使用 attributes(fm1) 吗？ - David Arenburg
@DavidArenburg 我之前没有找到，所以我试了一下，发现第二个项目是 fm1@devcomp$cmp[10]，但我仍然找不到第一个。 - Joe King
该类的拼写为 lmerMod。 - IRTFM
1
紧密相关：https://dev59.com/HYLba4cB1Zd3GeqPZhnu#25085736 - Ben Bolker
attr(summary(fm1)$varcor, "sc")：返回 fm1 模型的方差相关性的摘要属性中的 "sc" 值。attr(summary(fm1)$varcor$Batch, "stddev")：返回 fm1 模型中 "Batch" 方差相关性的摘要属性中的 "stddev" 值。 - Vlo
1个回答
5

lme4提供VarCorr函数来提取方差和相关系数:

varcor <- VarCorr(fm1)
as.data.frame(varcor)
#       grp        var1 var2    vcov   sdcor
#1    Batch (Intercept) <NA> 1764.05 42.0006
#2 Residual        <NA> <NA> 2451.25 49.5101
2
如果提问者能够建议在文档的哪个位置查找/期望找到此信息，我可以尝试改进文档，使其更易于找到... - Ben Bolker
好的。我按顺序查看了 ?lmer，在索引中查看了 ?'lmerMod-class'，并尝试了 ?summary.lmerMod（尽管最后一个不存在）。我还考虑了 vcov，但没有去那里，因为上下文表明它仅适用于固定效应。看来我也许应该去那里，因为将来可能会有一些有帮助的注释者在那里放置链接。 - IRTFM
1
我认为不方便的地方是 class(fm1) 是 lmerMod，但你想要的方法是 VarCorr.merMod 类。当你执行 methods(class="lm") 时，会返回一堆有用的函数，但是 methods(class="lmerMod") 却没有返回任何东西。除非用户知道 VarCorr 函数（通过阅读文档），或者在查看 str 后提取了 summary(fm1)$varcor，否则用户很难找到他们想要的信息。当我在去年遇到这个问题时，我使用了超级繁琐的 capture.output 到文本处理的方法。编辑以反映 lmerMod。 - Vlo
@Vlo 尝试使用 methods(class='merMod')，虽然我不确定没有这个讨论的天真用户（比如我自己）怎么能知道 lmerMod 是 merMod 的子集。 - IRTFM
做到了，谢谢。但是提取一个值真是太麻烦了！ - Phil
显示剩余2条评论
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接