代码之家 › 专栏 › 技术社区 › user55546

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

nlme mixed-models lme4 regression r

user55546 · 技术社区 · 2 年前

我正试图使用 lme() 的函数 nlme 包裹这个想法是重写 barleyprogeny1.lmer 中的模型 nlme() .

我在这篇文章中找到了一些信息: multiple random effects in 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 vs。 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")

0 回复 | 直到 2 年前

Ben Bolker 2 年前

你想要什么的标准术语是 交叉随机效应 ,这是一个很难实现的问题 lme .

适应于 this CV question 最终从 2003 conversation on r-help between Doug Bates and Peter Dalgaard :

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 .