[R-sig-ME] lme: predictions variance collapses when one more level is added

[Dieter Menne]

I have a simple mixed-model, with predictive factor treat (levels M1,M2,M3,
M4), continuous par, and a grouping variable subj from a cross-over experiment.

Everything works as expected when I only use M1, M2, M3; see subset.lme
below. The residuals are well distributed; resid(.,type="p")~fitted(.)|treat

When I add level M4 (all.lme below), the variance of the predictions shrinks
to almost zero. I know that level M4 adds heteroscedasticity, so I tried
with varPower(); this corrects for the residual, but the fitted() appear
nonsensical. 


d = structure(list(subj = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 
2L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 
7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 
11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 14L, 
14L, 14L, 14L, 15L, 15L, 15L, 15L, 16L, 17L, 17L, 17L, 17L), .Label = c("1", 
"3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", 
"15", "16", "17", "18"), class = "factor"), treat = structure(c(4L, 
1L, 3L, 2L, 2L, 3L, 1L, 4L, 3L, 2L, 4L, 1L, 4L, 3L, 2L, 3L, 4L, 
2L, 1L, 1L, 4L, 2L, 3L, 4L, 3L, 1L, 2L, 1L, 4L, 3L, 2L, 2L, 3L, 
4L, 1L, 2L, 4L, 3L, 4L, 3L, 1L, 2L, 2L, 3L, 4L, 1L, 2L, 1L, 3L, 
4L, 4L, 1L, 3L, 2L, 1L, 2L, 4L, 3L, 1L, 3L, 2L, 4L, 1L), .Label = c("M1", 
"M2", "M3", "M4"), class = "factor"), par = c(162.1, 105.2, 146.1, 
119.7, 129.9, 152.2, 156.8, 235.4, 122.6, 107.6, 199, 82.5, 165.8, 
115.3, 128.4, 96.8, 334.3, 67.9, 84.3, 152.3, 63.7, 103.6, 97.7, 
335.2, 130.5, 113.6, 123.9, 79.7, 172.1, 94.7, 97, 101.9, 125.4, 
177.1, 96, 90.9, 174.2, 96.3, 343.6, 106.3, 94.9, 112.5, 111.6, 
110.4, 197.9, 116, 93, 85.1, 123.1, 222, 28.4, 123.3, 111.6, 
106.3, 81.8, 101.5, 301.1, 102, 59.8, 131.1, 99.2, 104.6, 107.9
)), .Names = c("subj", "treat", "par"), row.names = c(1L, 2L, 
3L, 5L, 7L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 18L, 19L, 
20L, 21L, 22L, 23L, 25L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 35L, 
36L, 38L, 39L, 40L, 42L, 43L, 45L, 46L, 47L, 48L, 49L, 50L, 51L, 
52L, 54L, 55L, 57L, 58L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 68L, 
70L, 71L, 72L, 73L, 74L, 75L, 76L, 77L, 79L), class = "data.frame")

library(nlme)
str(d)
# Drop one level first to see reasonable result
subset.lme = lme(par~treat,random=~1|subj,data=d[d$treat!="M4",])
# This plot looks good (scaling chosen for comparison below)
plot(subset.lme,  resid(.,type="p")~fitted(.)|treat,pch=16,xlim=c(70,250))
summary(subset.lme)

# Use all levels: variance of predictions per group shrinks to nothing
all.lme = lme(par~treat,random=~1|subj,data=d)
plot(all.lme,  resid(.,type="p")~fitted(.)|treat,pch=16,xlim=c(70,250))
summary(all.lme)

# Taking into account heteroscedastic M4 does not change the picture
allpower.lme = lme(par~treat,random=~1|subj,data=d,weights=varPower())
plot(allpower.lme,  resid(.,type="p")~fitted(.)|treat,pch=16,xlim=c(70,250))
summary(allpower.lme)


sessionInfo()


#R version 3.0.2 (2013-09-25)
#Platform: i386-w64-mingw32/i386 (32-bit)

#locale:
#  [1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252   
LC_MONETARY=German_Germany.1252
#[4] LC_NUMERIC=C                    LC_TIME=German_Germany.1252    
#
#attached base packages:
#  [1] stats     graphics  grDevices utils     datasets  methods   base     

#other attached packages:
#  [1] nlme_3.1-111

...