[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:
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
More information about the R-sig-mixed-models
mailing list