我試圖從由混合效果模型的概要調用創建的對象中包含的隨機效果表中提取單個元素。具體而言,我想提取每個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-04
和2.484345e-05
分別)它必須是在某處,我通過搜索str(modelLME)
,但我看了找不到它。
你是對的@beetroot。謝謝你幫我找到答案,它完美地回答了我的問題。 – llewmills