lme() 函数中截距的两个随机效应

问题描述 投票：0回答：1

我正在尝试使用 

lme()
包中的 
nlme
函数将两个随机效应引入到截距中。这个想法是在 
barleyprogeny1.lmer
中重写 
nlme()
模型。

我在这篇文章中找到了一些信息：nlme中的多重随机效应和https://stat.ethz.ch/pipermail/r-sig-mixed-models/2018q2/026750.html， 然后我决定根据信息调整模型，它们看起来像这样：

library(nlme)
library(lme4)
barleyprogeny1.lme <- lme(yield_g_m2 ~ 1,
                           random = list(fblock = ~ 1, fline = ~1),
                           method="REML",
                           dataexperiment)
barleyprogeny1.lmer <- lmer(yield_g_m2 ~ 1 + 
                            (1|fblock) + (1|fline), 
                            REML=TRUE,
                            data = dataexperiment) 
barleyprogeny2.lme <- lme(yield_g_m2 ~ 1,
                          random = ~1|fblock,
                          method="REML",
                          data = dataexperiment)                         
barleyprogeny2.lmer <- lmer(yield_g_m2 ~ 1 + 
                          (1|fblock), 
                          REML=TRUE,
                          data = dataexperiment)

但是，当比较这些模型是否相同时，看到有显着差异：

all.equal(REMLcrit(barleyprogeny1.lmer), c(-2*logLik(barleyprogeny1.lme))) 
"Mean relative difference: 0.021192723770312"

all.equal(fixef(barleyprogeny1.lme), fixef(barleyprogeny1.lmer))
"Mean relative difference: 0.016793324438639"

此外，例如，当我比较 

barleyprogeny1.
中的两个模型 
nlme
与 
lme4
中的对数似然之间的差异时，它们是非常不同的， 在我看来，
fline = ~1
中的效果
barleyprogeny1.lme
被忽略了：

-2*logLik(barleyprogeny2.lmer)
'log Lik.' 1912.0625636471 (df=3)   #In this case they are equal
-2*logLik(barleyprogeny2.lme)
'log Lik.' 1912.0625636472 (df=3)
#################################
-2*logLik(barleyprogeny1.lmer)
'log Lik.' 1872.3816955801 (df=4)   # In this case they are not equal.
-2*logLik(barleyprogeny1.lme)
'log Lik.' 1912.0625636471 (df=4)

我对 

lme()
函数感兴趣的事实是因为与 
lmer()
相比，该函数允许一些可能性。我在本帖开头提到的帖子是 2019 年的。此修复有任何实现或想法吗？

数据在这里：

structure(list(block = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L), line = c(2L, 39L, 41L, 4L, 79L, 76L, 78L, 
35L, 25L, 67L, 17L, 30L, 73L, 64L, 72L, 44L, 18L, 82L, 29L, 23L, 
38L, 31L, 40L, 57L, 63L, 83L, 36L, 66L, 10L, 13L, 62L, 81L, 43L, 
26L, 45L, 12L, 24L, 51L, 6L, 14L, 80L, 46L, 42L, 47L, 37L, 27L, 
54L, 32L, 70L, 9L, 8L, 77L, 22L, 11L, 58L, 19L, 49L, 28L, 50L, 
74L, 7L, 68L, 59L, 16L, 33L, 53L, 75L, 20L, 61L, 1L, 71L, 55L, 
69L, 60L, 21L, 65L, 34L, 5L, 48L, 56L, 3L, 15L, 52L, 64L, 27L, 
29L, 68L, 16L, 33L, 67L, 51L, 66L, 50L, 36L, 80L, 17L, 26L, 63L, 
74L, 56L, 81L, 43L, 79L, 14L, 32L, 70L, 41L, 31L, 71L, 65L, 7L, 
58L, 69L, 34L, 11L, 44L, 49L, 54L, 72L, 9L, 23L, 48L, 59L, 78L, 
61L, 45L, 18L, 19L, 35L, 6L, 40L, 8L, 1L, 10L, 42L, 46L, 37L, 
60L, 38L, 13L, 24L, 39L), yield_g_m2 = c(483.33, 145.84, 321.84, 
719.14, 317.63, 344.48, 260.02, 374.28, 428.61, 407.25, 551.84, 
353.29, 355.3, 647.92, 165.76, 517.52, 366.24, 251.3, 606.37, 
605.75, 641.42, 166.75, 410.87, 181.97, 562.9, 280.44, 800.35, 
687.92, 764.88, 541.15, 730.48, 315.63, 678.46, 580.22, 519.88, 
436.59, 671.22, 692.55, 849.66, 910.76, 487.86, 724.01, 793.43, 
192.43, 895.3, 731.87, 809.41, 669.16, 996.19, 774.84, 636.45, 
357.94, 340.65, 644.83, 521.67, 622.72, 830.57, 679.92, 721.13, 
489.31, 907.38, 325.96, 553.46, 210.71, 770.23, 559.14, 617.33, 
632.46, 611.52, 717.78, 595.86, 555.29, 467.24, 572.9, 514.62, 
818.74, 673.43, 798.99, 786.06, 522.61, 873.04, 600.06, 603.04, 
681.64, 762.42, 932.33, 385.47, 240, 846.85, 702.58, 746.11, 
846.05, 885.67, 1054.7, 478.12, 959.25, 639.39, 755.9, 551.41, 
435.62, 303.72, 836.82, 439.17, 934.72, 836.95, 904.9, 538, 226.12, 
569.61, 713.43, 820.08, 435.34, 378.89, 639.11, 516.84, 873.18, 
823.25, 859.36, 258.59, 587.07, 817.51, 645.1, 634.58, 260.26, 
472.44, 575.76, 265.37, 423.76, 554.69, 755.05, 568.31, 299.92, 
591.19, 756.63, 552.53, 627.25, 552.7, 284.72, 540.68, 475.1, 
463.22, 212.66), fblock = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("1", "2"), class = "factor"), 
    fline = structure(c(2L, 39L, 41L, 4L, 79L, 76L, 78L, 35L, 
    25L, 67L, 17L, 30L, 73L, 64L, 72L, 44L, 18L, 82L, 29L, 23L, 
    38L, 31L, 40L, 57L, 63L, 83L, 36L, 66L, 10L, 13L, 62L, 81L, 
    43L, 26L, 45L, 12L, 24L, 51L, 6L, 14L, 80L, 46L, 42L, 47L, 
    37L, 27L, 54L, 32L, 70L, 9L, 8L, 77L, 22L, 11L, 58L, 19L, 
    49L, 28L, 50L, 74L, 7L, 68L, 59L, 16L, 33L, 53L, 75L, 20L, 
    61L, 1L, 71L, 55L, 69L, 60L, 21L, 65L, 34L, 5L, 48L, 56L, 
    3L, 15L, 52L, 64L, 27L, 29L, 68L, 16L, 33L, 67L, 51L, 66L, 
    50L, 36L, 80L, 17L, 26L, 63L, 74L, 56L, 81L, 43L, 79L, 14L, 
    32L, 70L, 41L, 31L, 71L, 65L, 7L, 58L, 69L, 34L, 11L, 44L, 
    49L, 54L, 72L, 9L, 23L, 48L, 59L, 78L, 61L, 45L, 18L, 19L, 
    35L, 6L, 40L, 8L, 1L, 10L, 42L, 46L, 37L, 60L, 38L, 13L, 
    24L, 39L), .Label = c("1", "2", "3", "4", "5", "6", "7", 
    "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", 
    "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", 
    "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", 
    "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", 
    "48", "49", "50", "51", "52", "53", "54", "55", "56", "57", 
    "58", "59", "60", "61", "62", "63", "64", "65", "66", "67", 
    "68", "69", "70", "71", "72", "73", "74", "75", "76", "77", 
    "78", "79", "80", "81", "82", "83"), class = "factor"), plot = structure(c(1L, 
    2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 
    15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 
    27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L, 38L, 
    39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 50L, 
    51L, 52L, 53L, 54L, 55L, 56L, 57L, 58L, 59L, 60L, 61L, 62L, 
    63L, 64L, 65L, 66L, 67L, 68L, 69L, 70L, 71L, 72L, 73L, 74L, 
    75L, 76L, 77L, 78L, 79L, 80L, 81L, 82L, 83L, 1L, 2L, 3L, 
    4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 
    17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 
    29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 
    41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 
    53L, 54L, 55L, 56L, 57L, 58L, 59L), .Label = c("1", "2", 
    "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", 
    "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", 
    "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", 
    "34", "35", "36", "37", "38", "39", "40", "41", "42", "43", 
    "44", "45", "46", "47", "48", "49", "50", "51", "52", "53", 
    "54", "55", "56", "57", "58", "59", "60", "61", "62", "63", 
    "64", "65", "66", "67", "68", "69", "70", "71", "72", "73", 
    "74", "75", "76", "77", "78", "79", "80", "81", "82", "83"
    ), class = "factor")), row.names = c(NA, -142L), class = "data.frame")
r regression lme4 mixed-models nlme
1个回答
1
投票

您想要的标准术语是交叉随机效应，众所周知，在

lme
中实现起来很痛苦。

改编自这个简历问题，最终改编自 2003 年 Doug Bates 和 Peter Dalgaard 之间关于 r-help 的对话：

dataexperiment$int <- 1  ## create dummy variable
barleyprogeny1.lme <- lme(yield_g_m2 ~ 1,
                          random = list(int = pdIdent(form = ~ 0 + factor(fblock)), 
                                        fline = ~1),
                          method="REML",
                          dataexperiment)

barleyprogeny1.lmer <- lmer(yield_g_m2 ~ 1 + 
                            (1|fblock) + (1|fline), 
                            REML=TRUE,
                            data = dataexperiment)

结果比较：

VarCorr(barleyprogeny1.lme)

                Variance                    StdDev   
int =           pdIdent(0 + factor(fblock))          
factor(fblock)1    14.61884                   3.82346
factor(fblock)2    14.61884                   3.82346
fline =         pdLogChol(1)                         
(Intercept)     30645.45039                 175.05842
Residual        13221.74805                 114.98586

VarCorr(barleyprogeny1.lmer)

 Groups   Name        Std.Dev.
 fline    (Intercept) 175.0585
 fblock   (Intercept)   3.8228
 Residual             114.9858

不相同，但相对差异仅约2e-6；对数似然相当于标准容差 

all.equal()
。

从之前的答案来看：因为

lme
只能做嵌套随机效果，所以你必须

...[技巧] lme 通过创建一个具有单个级别的组。该模型仍然是嵌套的，但现在它是作为嵌套一部分的单级别组，这没有问题。

通过将一组的随机效应视为斜率而不是截距的参数来合并交叉随机效应。

pdIdent
规范规定，
fblock
随机效应的协方差矩阵应该是齐次且对角的（根据文档，“身份正定矩阵的倍数”）；这模仿了如果您可以指定截距项会发生什么
1|fline
。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.