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

[Nick Isaac]

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:

w<-rep(1,nrow(sleepstudy)
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))}
a<-sapply(c(fm1,fm2,fm3,fm4),myfunc)
colnames(a)<-c("unweighted","w=1","sum(w)=1","w=2")
rownames(a)<-c("sigmaML^2","pwrss/n","pwrss/sum(w)","pwrss/sum(w^2)")

a

               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
#1).

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