I'm trying to introduce two random effects into the intercept using the lme() function from the nlme package. The idea is to rewrite the barleyprogeny1.lmer model in nlme().
I found some information in this post: multiple random effects in nlme and https://stat.ethz.ch/pipermail/r-sig-mixed-models/2018q2/026750.html, and then I decided to adjust the models based on the information and they looked like this:
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)
However, when comparing whether these models are the same, see that there is a significant difference:
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"
Furthermore, when I compare, for example, the difference between the Log-Likelihood of both models barleyprogeny1. in nlme vs. lme4, they are quite different,
what seems to me is that the effect fline = ~1 in barleyprogeny1.lme is ignored:
-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)
The fact that I am interested in the lme() function is due to some possibilities that this function allows in contrast to lmer(). The post I mention at the beginning of this thread is from 2019. Is there any implementation or idea for this fix?
The data is here:
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")
The standard term for what you want is crossed random effects, which are known to be a pain to implement in lme.
Adapting from this CV question and ultimately from a 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)
Comparing results:
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
Not identical, but the relative difference is only about 2e-6; the log-likelihoods are equivalent up to the standard tolerance of all.equal().
From the previous answer: because lme can only do nested random effects you have to
... [trick] lme by creating a group with a single level. The model is still nested but now it is the single level group that is part of the nesting which is no problem.
The crossed random effect is incorporated by treating the random effect for the one group as a parameter for a slope instead of an intercept.
The pdIdent specification says that the covariance matrix for the fblock random effect should be homogeneous and diagonal ("a multiple of the identity positive-definite matrix", according to the documentation); this mimics what would happen if you could specify an intercept term 1|fline.