[R-sig-ME] explaining lme variance component results
Mike Dunbar
mdu at ceh.ac.uk
Mon Sep 3 14:45:54 CEST 2007
Dear All
Ben's recent question "boneheaded (?) question about SS in anova (lmer vs lme)", prompted me to post this very marginally related question. Time and time again I come across problems presenting results of multilevel models fitted with lme to either scientists not familiar with random effects, or only familiar with random effects in the traditional balanced, linear anova context. It's the latter I'm having particular troubles with at the moment.
Recently I've had a manuscript come back from review from an ecology journal where the reviewers and also some of my co-workers fall into the latter category. The analysis in brief is a variance components analysis, with 4 nested random effects, hence no fixed effects. I illustrated the results showing the variance component at each level, including the residual. I then undertook likelihood ratio tests for each random effect by deletion. The results are a little strange in that in one instance a higher level component (MONTH: 25% of the total variance) is significant, but a lower level component (POLE: 29% of the total variance) is not. The results are also a little strange as for the highest level factor, MONTH, there are only 4 levels, hence I'd have thought that its variance component would be fairly imprecise and hence less likely to be significant. I'm prepared to believe these results as I think the deletion tests are entirely sound, but I'm not sure if I understand the reasons. I think its inherent in the nested design and the degrees of freedom, but I'm having a hell of a job explaining it without alot of arm waving. One way of course would be to provide an explanation using F-tests and expected mean squares, but I don't know how to do that for such a complex design.
Data and analysis below. If the data get mangled somehow then I'm happy to email as an attachment.
I posted a similar question on the multilevel list a few weeks ago, but had no replies, hence this is another try to a different audience with a slightly different question. I'm aware this isn't a comparison of lmer and lme (hope that doesn't matter), and I'm also aware that a GLMM analysis may be more appropriate, but I want to understand fully this more basic analysis first!
regards, and thanks in advance for any replies
Mike Dunbar
# analysis (run data below first)
varcor.insects <- lme(log(insectdens+1) ~ 1, random=~1|MONTH/TIME/TRANSECT/POLE, data=temp)
VarCorr(varcor.insects)
# variance components: MONTH: 0.757, TIME: 1.074, TRANSECT: 0.000, POLE: 0.898, resid: 0.307
varcor.insects.nopole <- lme(log(insectdens+1) ~ 1, random=~1|MONTH/TIME/TRANSECT, data=temp)
varcor.insects.nomonth <- lme(log(insectdens+1) ~ 1, random=~1|TIME/TRANSECT/POLE, data=temp)
varcor.insects.notime <- lme(log(insectdens+1) ~ 1, random=~1|MONTH/TRANSECT/POLE, data=temp)
anova(varcor.insects.nomonth, varcor.insects)
anova(varcor.insects.notime, varcor.insects)
anova(varcor.insects.nopole,varcor.insects)
# variance component for pole is greater than that for month, yet month is significant, pole is not
# data: run this first
temp <-
structure(list(MONTH = 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,
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, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("4", "5", "6", "7"), class = "factor"),
TRANSECT = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 3L, 3L, 3L,
3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 2L, 2L, 3L,
3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 1L,
2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 1L,
1L, 1L, 1L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L,
5L, 5L, 1L, 1L, 1L, 1L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L,
4L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 1L, 2L, 2L, 3L, 3L, 3L, 3L,
4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 1L, 2L, 2L, 3L,
3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 1L,
2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 1L, 1L,
1L, 1L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L,
5L, 1L, 1L, 1L, 1L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L,
5L, 5L, 5L, 5L, 1L, 1L, 1L, 1L, 2L, 2L, 3L, 3L, 3L, 3L, 4L,
4L, 4L, 4L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 1L, 2L, 2L, 3L, 3L,
3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 1L, 2L,
2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 1L, 1L,
1L, 1L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L,
5L, 1L, 1L, 1L, 1L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L,
5L, 5L, 5L, 5L, 1L, 1L, 1L, 1L, 2L, 2L, 3L, 3L, 3L, 3L, 4L,
4L, 4L, 4L, 5L, 5L, 5L, 5L), .Label = c("1", "2", "3", "4",
"5"), class = "factor"), POLE = structure(c(1L, 2L, 3L, 4L,
5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L,
18L, 2L, 3L, 4L, 5L, 5L, 7L, 8L, 9L, 10L, 11L, 12L, 13L,
14L, 15L, 16L, 17L, 18L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L,
9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 1L, 2L,
3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L,
16L, 17L, 18L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L,
12L, 13L, 14L, 15L, 16L, 17L, 18L, 1L, 2L, 3L, 4L, 5L, 6L,
7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L,
1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L,
15L, 16L, 17L, 18L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L,
11L, 13L, 14L, 15L, 16L, 17L, 18L, 1L, 2L, 3L, 4L, 5L, 6L,
7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L,
1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L,
15L, 16L, 17L, 18L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L,
11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 1L, 2L, 3L, 4L, 5L,
6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L,
1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L,
15L, 16L, 17L, 18L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L,
11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 1L, 2L, 3L, 4L, 5L,
6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L,
1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L,
15L, 16L, 17L, 18L), .Label = c("11", "12", "13", "14", "23",
"24", "31", "32", "33", "34", "41", "42", "43", "44", "51",
"52", "53", "54"), class = "factor"), TIME = structure(c(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, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 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, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 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, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 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, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L
), .Label = c("1", "2", "3", "4"), class = "factor"), insectdens = c(0,
63.64, 14.57, 22.5, 0, 107.6, 87.16, 22.24, 51.92, 42.1,
0, 15.78, 43.02, 9.98, 7.47, 38.92, 11.78, 0, 44.61, 45.9,
26.78, 80.2, 37.4, 17.58, 8.94, 54.7, 141.89, 10.39, 14.6,
158.34, 9.52, 14.36, 47.94, 12.26, 31.44, 104.37, 74.4, 73.38,
94.36, 120.26, 48.12, 264.87, 129.36, 107.03, 145.53, 95.2,
116.55, 492, 140.4, 111.28, 104.8, 59.76, 91.92, 12.92, 22.6,
53.6, 102.6, 151.92, 0, 0, 34.96, 21.38, 60.4, 16.9, 24.96,
263.76, 37.12, 30.76, 26.24, 25.88, 7.2, 48.87, 28.1, 44.28,
67.26, 0, 50.49, 0, 0, 10.82, 27.06, 9.09, 13.56, 13.25,
31.74, 11.98, 12.08, 30.55, 31.16, 16.35, 78.15, 13.54, 80.44,
104.55, 37.32, 107.7, 21.52, 22.28, 55.6, 40.72, 15.12, 25.44,
87.36, 19.94, 78.32, 13.04, 40.55, 14.37, 34.68, 69.4, 62.28,
117.96, 18.27, 34.74, 69.86, 100.89, 95.7, 0, 129.6, 147.4,
92.16, 76.1, 52.78, 52.32, 46.7, 49.41, 54.8, 59.2, 0, 0,
15.89, 35.42, 8.54, 52.98, 48.81, 49.85, 64.6, 20.2, 80.29,
38.52, 17.28, 41.55, 237.25, 38.88, 25.69, 0, 0, 0, 467.64,
36, 112.05, 42.08, 26.86, 0, 17.48, 34.95, 43.88, 33.76,
23.24, 16.29, 189.99, 365.6, 329.29, 228, 140.91, 433.84,
228.9, 130.3, 609.9, 371.88, 360.36, 338.91, 139.15, 285.9,
253.68, 353.35, 839.16, 717.42, 2081.2, 1052.03, 1276.65,
512.25, 614.46, 479.52, 3020.4, 885.96, 932.49, 1476.09,
576.46, 568.8, 712.53, 1864.56, 792.05, 1807.52, 899.25,
1487.7, 166.5, 78.7, 169.2, 99.35, 176.85, 104.6, 0, 193.05,
204.49, 201.9, 49.71, 93.06, 75.72, 207.79, 73.26, 96.4,
227.63, 141.7, 98.25, 58.4, 30.84, 141.72, 277.16, 534.19,
508.04, 68.44, 215.37, 0, 48.68, 81.36, 65.07, 133.5, 70.28,
79.62, 107.73, 30.14, 407.76, 422.24, 317.7, 241.5, 644.49,
1104.24, 268.24, 449.25, 475.57, 311.63, 461.3, 247.69, 395.25,
399.84, 529.48, 440.82, 394.56, 322, 353.5, 452.4, 699.72,
89.04, 347.6, 563.67, 456.88, 513, 440.37, 398.86, 428, 398.06,
438.75, 367.9, 501.48, 522.6, 616.11, 309.96, 232.08, 198.06,
109.59, 49.59, 152.08, 617.83, 372.75, 66.81, 68.28, 157.48,
46.24, 99.18, 158.05, 91.68, 112.62, 98.21, 117.54, 77.3)), .Names = c("MONTH",
"TRANSECT", "POLE", "TIME", "insectdens"), class = "data.frame", row.names = c(1L,
4L, 7L, 9L, 13L, 16L, 19L, 22L, 25L, 28L, 31L, 34L, 37L, 40L,
43L, 44L, 46L, 49L, 55L, 58L, 60L, 64L, 67L, 70L, 73L, 76L, 79L,
82L, 85L, 88L, 91L, 94L, 95L, 97L, 100L, 103L, 106L, 109L, 111L,
115L, 118L, 121L, 124L, 127L, 130L, 133L, 136L, 139L, 142L, 145L,
146L, 148L, 151L, 154L, 157L, 160L, 162L, 166L, 169L, 172L, 175L,
178L, 181L, 184L, 187L, 190L, 193L, 196L, 197L, 199L, 202L, 205L,
208L, 211L, 213L, 217L, 220L, 223L, 226L, 229L, 232L, 235L, 238L,
241L, 244L, 247L, 248L, 250L, 253L, 256L, 259L, 262L, 264L, 268L,
271L, 274L, 277L, 280L, 283L, 286L, 289L, 292L, 295L, 298L, 299L,
301L, 304L, 307L, 310L, 313L, 315L, 319L, 322L, 325L, 328L, 331L,
334L, 337L, 340L, 343L, 346L, 349L, 350L, 352L, 355L, 358L, 361L,
364L, 366L, 370L, 373L, 376L, 379L, 382L, 385L, 388L, 394L, 397L,
400L, 401L, 403L, 406L, 409L, 412L, 415L, 417L, 421L, 424L, 427L,
430L, 433L, 436L, 439L, 442L, 445L, 448L, 451L, 452L, 454L, 457L,
460L, 463L, 466L, 468L, 472L, 475L, 478L, 481L, 484L, 487L, 490L,
493L, 496L, 499L, 502L, 503L, 505L, 508L, 511L, 514L, 517L, 519L,
523L, 526L, 529L, 532L, 535L, 538L, 541L, 544L, 547L, 550L, 553L,
554L, 556L, 559L, 562L, 565L, 568L, 570L, 574L, 577L, 580L, 583L,
586L, 589L, 592L, 595L, 598L, 601L, 604L, 605L, 607L, 610L, 613L,
616L, 619L, 621L, 625L, 628L, 631L, 634L, 637L, 640L, 643L, 646L,
649L, 652L, 655L, 656L, 658L, 661L, 664L, 667L, 670L, 672L, 676L,
679L, 682L, 685L, 688L, 691L, 694L, 697L, 700L, 703L, 706L, 707L,
709L, 712L, 715L, 718L, 721L, 723L, 727L, 730L, 733L, 736L, 739L,
742L, 745L, 748L, 751L, 754L, 757L, 758L, 760L, 763L, 766L, 769L,
772L, 774L, 778L, 781L, 784L, 787L, 790L, 793L, 796L, 799L, 802L,
805L, 808L, 809L, 811L, 814L))
--
This message (and any attachments) is for the recipient only...{{dropped}}
More information about the R-sig-mixed-models
mailing list