explaining lme variance component results
Mike, I find asinh(x/2) to be an interesting (more palatable?) alternative to log(x+1). A couple of other comments inserted below: Kevin Wright
On 9/3/07, Mike Dunbar <mdu at ceh.ac.uk> wrote:
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 th! e 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.
In my experience with biological data, factors that have levels that
are widely (temporal/spatial) separated are often more variable than
factors with levels that are closer together. In your data, I would
(a priori) expect MONTH is more likely to be significant than POLE.
Have you plotted this data? I often find dotplots very useful for
exploring/understanding data that have a few factors. For example:
levels(temp$MONTH) <- paste("Mo",levels(temp$MONTH),sep="")
levels(temp$POLE) <- paste("Po",levels(temp$POLE),sep="")
levels(temp$TRANSECT) <- paste("Tr",levels(temp$TRANSECT),sep="")
levels(temp$TIME) <- paste("Ti",levels(temp$TIME),sep="")
# Time2 and Time3 have big differences between months
dotplot(TRANSECT~y|TIME*MONTH, data=temp)
# Pole-to-pole variation is often small
dotplot(POLE~y|MONTH*TIME*TRANSECT, data=temp,
scales=list(y=list(cex=.5,relation="free")))
I often try many different arrangements of the panels/lines to find
the most informative view.
The zero values of insectdens look odd, given the
not-so-unlike-Gaussian distribution of the other values.
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}}
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models