2016-05-02 59 views
0

我試圖從由混合效果模型的概要調用創建的對象中包含的隨機效果表中提取單個元素。具體而言,我想提取每個2級隨機效果。使用nlme從混合效果模型的輸出中提取元素

玩具數據:

set.seed(1234) 
score <- c(rnorm(8, 20, 3), rnorm(8, 35, 5)) 
rep <- rep(c(0,1,2,3), each = 8) 
group <- rep(0:1, times = 16) 
id <- factor(rep(1:8, times = 4)) 

df <- data.frame(id, group, rep, score) 

現在創建一個模型

require(nlme) 

modelLME <- summary(lme(score ~ group*rep, data = df, random = ~ rep|id)) 

modelLME 

當我們調用它,我們得到的輸出

Linear mixed-effects model fit by REML 
Data: df 
     AIC  BIC logLik 
    219.6569 230.3146 -101.8285 

Random effects: 
Formula: ~rep | id 
Structure: General positive-definite, Log-Cholesky parametrization 
      StdDev  Corr 
(Intercept) 2.664083e-04 (Intr) 
rep   2.484345e-05 0  
Residual 7.476621e+00  

Fixed effects: score ~ group * rep 
       Value Std.Error DF t-value p-value 
(Intercept) 22.624455 3.127695 22 7.233587 0.0000 
group  -1.373324 4.423229 6 -0.310480 0.7667 
rep   2.825635 1.671823 22 1.690152 0.1051 
group:rep 0.007129 2.364315 22 0.003015 0.9976 
Correlation: 
      (Intr) group rep 
group  -0.707    
rep  -0.802 0.567  
group:rep 0.567 -0.802 -0.707 

Standardized Within-Group Residuals: 
     Min   Q1   Med   Q3   Max 
-1.86631781 -0.74498367 0.03515508 0.76672652 1.91896578 

Number of Observations: 32 
Number of Groups: 8 

現在我可以提取隨機效應殘差以上表格通過

modelLME$sigma 

但我不能在這個輸出發現在(Intercept)rep行隨機效應表StdDev列的值(2.664083e-042.484345e-05分別)它必須是在某處,我通過搜索str(modelLME),但我看了找不到它。

+0

你是對的@beetroot。謝謝你幫我找到答案,它完美地回答了我的問題。 – llewmills

回答

0

你想要這樣的東西嗎?

library(nlme) 
library(broom) 

modelLME = lme(score ~ group*rep, data = df, random = ~ rep|id) 
tidy(modelLME) 
+0

嗨@bramtayl,謝謝你的回答,但我正在尋找的答案是在@beetroot發現的重複問題中。 'nlme'和'lme4'中的'VarCorr'函數爲您提供了整個殘差表,然後您可以使用'[,]'subsetter來提取您需要的元素。 – llewmills

相關問題