heterogeneous variance models

Also, the results don't seem to agree with the results from the old nlme package.

Example (using a slightly more sensible model than from the help file): Here's what the data looks like:

  xyplot(distance~age|Sex, groups=Subject, data=Orthodont, type="l", col="darkblue")

Let's fit lme4 and nlme models:

  l_new1 = lmer(distance ~ age + Sex + (age|Subject), data = Orthodont)
  l_new2 = update(l_new1, . ~ . + (0+dummy(Sex, "Female")|Subject))
  summary(l_new1)
  summary(l_new2)

  l_old1 = lme(distance ~ age + Sex, random=~age|Subject, data = Orthodont)
  l_old2 = update(l_old1, weights=varIdent(form=~1|Sex))
  summary(l_old1)
  summary(l_old2)

While the fixed effects estimates for the two models assuming homoscedasticity are very ~identical, the estimates from the heteroscedasticity models differ:

fixef(l_new1)

(Intercept)         age   SexFemale 
 17.6351989   0.6601852  -2.1454883

fixef(l_new2)

(Intercept)         age   SexFemale 
 17.6487393   0.6601852  -2.1787238

fixef(l_old1)

(Intercept)         age   SexFemale 
 17.6352002   0.6601852  -2.1454915

fixef(l_old2)

(Intercept)         age   SexFemale 
 18.4385143   0.5795744  -1.9407706

It's a clever trick, but doesn't seem to work properly (or is it nlme that doesn't work properly?). (Also, it assumes that one *knows* which group has the lowest residual variance.)


Regards,
Karl Ove Hufthammer


-----Opphavleg melding-----
From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces at r-project.org] P? vegner av Aaron Mackey
Sent: 11. mai 2015 15:50
To: Ben Bolker
Kopi: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] heterogeneous variance models

Thanks for this. What's confusing to me in the example is that the residual std. dev is still 1.3100 whether the dummy term is included or not; I would have thought that the residual sd without the dummy term would be intermediate between the two sex-specific sd's seen in the full model (reported as 1.0521 and 1.311).  The reported residuals also don't seem to change between the two models.

thanks for any clarification you can provide, -Aaron

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1


   Someone sent me an e-mail recently asking how/whether heterogeneous 
variances can be implemented in lme4.  I can't find it any more, so 
I'm answering here in hopes that they'll see it. The short answer is 
that they can, at least for heterogeneity among categorical groups 
(similar to what varIdent does in the nlme package), although it's a 
bit of a nuisance for categorical variables with lots of levels, as 
the variance terms all need to be specified separately.

   The trick is to set up a dummy variable for each level: there is a 
helper function dummy() for this purpose, and example("dummy") gives 
an example showing how to estimate heterogeneous residual variances 
for females vs. males in the standard sleep study example.

  Ben Bolker


-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.11 (GNU/Linux)

iQEcBAEBAgAGBQJVT73oAAoJEOCV5YRblxUHlzQIAMuoYv9K06NGG+3mEcvHcmPW
Xp7rLyxRfiTtRkmr/UieR52WpFDapxNgMurLT1aESD5TWeSBTXmeev2ltajPM5It
A5Mz1cNob1Tl67zsBCUrKMVVclRicO/Ymbpab1QJYcXxGAVZ03gkaguPWJpw3aHQ
kt2GZW6SGidbI17PwZXTyVb69cJ8p/mq459LwMAFNusD4MCiIt1/K2HTL8NnekVn
G7XPMDBg/JoyHmnqWpJbhzccjxZetNMtK6GyXaamlnyPJaJI0i9j6t/X0jKiAqO0
k1BdwNkLo03sneHodyhdRqAS5KMzeNSo8FOinDIUCLtVikz5k0IbQhHO6JWChz4=
=KPfA
-----END PGP SIGNATURE-----

_______________________________________________
R-sig-mixed-models at r-project.org mailing list 
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

_______________________________________________
R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models