[R-sig-ME] Bug in weights in lmer

Nick Isaac njbisaac at googlemail.com
Thu Apr 24 19:46:57 CEST 2008

Dear Prof Bates,

Thanks for your reply. My sympathies re your computer ailments.

Following your suggestion, I explored the relationship between sigma^2
and pwrss:

fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
fm2 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy, weights = w)
fm3 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy, weights = w/sum(w))
fm4 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy, weights = 2*w)

myfunc<-function(x){c(x at deviance[5]^2,x at deviance[7]/nrow(x at frame),x at deviance[7]/sum(x at pWt),x at deviance[7]/sum(x at pWt^2))}


               unweighted     w=1     sum(w)=1      w=2
sigmaML^2        647.6672 7026169 9.992688e+05 11571700
pwrss/n          647.6672 7026169 5.551493e+03 23143401
pwrss/sum(w)          Inf 7026169 9.992688e+05 11571700
pwrss/sum(w^2)        Inf 7026169 1.798684e+08  5785850

>  It is likely that this estimate should be  the pwrss divided by either
>  the sum of the elements in pWt or the sum of the squares of the
>  elements in pWt when we don't have unit weights.  Do either of those
>  numbers seem reasonable?

The unweighted model (fm1) behaves as it should: sigma^2 equals
pwrss/n and the empty pWt slot means that the other two entries return
Inf. When all observations receive weight=1 (fm2), all four estimates
are equal because n=sum(w)=sum(w^2)=180. However, these numbers are
far higher than the comparable values in fm1. The third and fourth
models, in which the observations are weighted equally (1/180 & 2
respectively) show that sigma^2 equals pwrss/sum(w) (your suggestion

But this doesn't explain why sigma should differ between fm1 and fm2
(they should be identical, right?). So I compared the contents of the
u slot:

plot(fm1 at u, fm2 at u)

Whilst fm1 at u is approximately normally distributed around zero, fm2 is
bimodal and all values are positive. It turns out that the even
elements of the two are perfectly correlated: fm2 at u = 158 + 2.4 *
fm1 at u:

summary(lm(fm2 at u[(1:18)*2] ~ fm1 at u[(1:18)*2]))

However, the relationship between the odd elements is much weaker:

summary(lm(fm2 at u[(1:18)*2-1] ~ fm1 at u[(1:18)*2-1]))

My understanding is not sufficient to know what this means, but I hope
it provides you with some useful clues.

Best wishes, Nick

More information about the R-sig-mixed-models mailing list