[R-sig-ME] Fitted values in lmer
Douglas Bates
bates at stat.wisc.edu
Fri Mar 14 23:11:34 CET 2008
Patch is now committed. I took the liberty of shortening it to a
one-liner and making similar fixes for the residuals and resid
methods.
Does anyone know enough about the definitions of the different types
of residuals for generalized linear models to implement a "type"
argument for the residuals method?
On Fri, Mar 14, 2008 at 4:38 PM, Douglas Bates <bates at stat.wisc.edu> wrote:
> Thank you, Henric. That's exactly what I wanted. I will change one
> glich, which is my glich, not yours. The original method should have
> returned the mu slot, not the eta slot. These are the same numerical
> values for linear mixed models but not for generalized linear mixed
> models.
>
>
>
> On Fri, Mar 14, 2008 at 4:28 PM, Henric Nilsson (Public)
> <nilsson.henric at gmail.com> wrote:
> >
> > Douglas Bates wrote:
> >
> > > What I was referring to is the na.action attribute of the model frame
> > >
> > >> library(lme4)
> > > Loading required package: Matrix
> > > Loading required package: lattice
> > >
> > > Attaching package: 'Matrix'
> > >
> > >
> > > The following object(s) are masked from package:stats :
> > >
> > > xtabs
> > >
> > >> slps <- sleepstudy
> > >> slps$Reaction[c(14, 34, 92)] <- NA
> > >> fm1 <- lmer(Reaction ~ Days + (Days|Subject), slps)
> > >> str(fm1 at frame)
> > > 'data.frame': 177 obs. of 3 variables:
> > > $ Reaction: num 250 259 251 321 357 ...
> > > $ Days : num 0 1 2 3 4 5 6 7 8 9 ...
> > > $ Subject : Factor w/ 18 levels "308","309","310",..: 1 1 1 1 1 1 1 1 1 1 ...
> > > - attr(*, "terms")=Classes 'terms', 'formula' length 3 Reaction ~ Days
> > > .. ..- attr(*, "variables")= language list(Reaction, Days)
> > > .. ..- attr(*, "factors")= int [1:2, 1] 0 1
> > > .. .. ..- attr(*, "dimnames")=List of 2
> > > .. .. .. ..$ : chr [1:2] "Reaction" "Days"
> > > .. .. .. ..$ : chr "Days"
> > > .. ..- attr(*, "term.labels")= chr "Days"
> > > .. ..- attr(*, "order")= int 1
> > > .. ..- attr(*, "intercept")= int 1
> > > .. ..- attr(*, "response")= int 1
> > > .. ..- attr(*, ".Environment")=<R_GlobalEnv>
> > > .. ..- attr(*, "predvars")= language list(Reaction, Days)
> > > .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
> > > .. .. ..- attr(*, "names")= chr [1:2] "Reaction" "Days"
> > > - attr(*, "na.action")=Class 'omit' Named int [1:3] 14 34 92
> > > .. ..- attr(*, "names")= chr [1:3] "14" "34" "92"
> > >
> > > At this point one can get the fitted values with the NA's in the
> > > correct positions via the thoroughly unobvious
> > >
> > >> stats:::naresid.exclude(attr(fm1 at frame, "na.action"), fm1 at mu)
> > > [1] 253.4849 273.1808 292.8766 312.5725 332.2684 351.9642 371.6601 391.3559
> > > [9] 411.0518 430.7476 213.1921 214.8561 216.5201 NA 219.8481 221.5121
> > > [17] 223.1761 224.8401 226.5041 228.1681 212.8984 217.8564 222.8143 227.7723
> > > [25] 232.7302 237.6882 242.6462 247.6041 252.5621 257.5200 275.9163 281.5164
> > > [33] 287.1166 NA 298.3168 303.9170 309.5171 315.1173 320.7174 326.3175
> > > [41] 273.5320 280.9433 288.3545 295.7658 303.1771 310.5884 317.9997 325.4110
> > > [49] 332.8223 340.2336 260.3778 270.5809 280.7841 290.9872 301.1903 311.3934
> > > [57] 321.5965 331.7996 342.0027 352.2058 268.1086 278.3689 288.6293 298.8896
> > > [65] 309.1500 319.4103 329.6707 339.9310 350.1914 360.4517 244.2257 255.7619
> > > [73] 267.2981 278.8342 290.3704 301.9066 313.4427 324.9789 336.5151 348.0512
> > > [81] 251.2783 250.9591 250.6398 250.3205 250.0013 249.6820 249.3627 249.0435
> > > [89] 248.7242 248.4050 283.8511 NA 322.6857 342.1031 361.5204 380.9377
> > > [97] 400.3550 419.7723 439.1897 458.6070 226.4061 238.0217 249.6373 261.2529
> > > [105] 272.8685 284.4841 296.0997 307.7153 319.3310 330.9466 238.3395 255.4259
> > > [113] 272.5123 289.5986 306.6850 323.7714 340.8578 357.9442 375.0305 392.1169
> > > [121] 256.0056 263.4524 270.8991 278.3458 285.7925 293.2393 300.6860 308.1327
> > > [129] 315.5795 323.0262 272.0277 286.0636 300.0995 314.1354 328.1714 342.2073
> > > [137] 356.2432 370.2791 384.3150 398.3510 254.6441 265.9888 277.3335 288.6782
> > > [145] 300.0229 311.3676 322.7123 334.0570 345.4017 356.7464 225.9405 241.2159
> > > [153] 256.4914 271.7668 287.0422 302.3176 317.5930 332.8684 348.1438 363.4192
> > > [161] 252.2315 261.7074 271.1833 280.6592 290.1351 299.6110 309.0870 318.5629
> > > [169] 328.0388 337.5147 263.5955 275.3631 287.1308 298.8985 310.6662 322.4338
> > > [177] 334.2015 345.9692 357.7368 369.5045
> > >
> > > There is some subtlety about the difference between na.omit and
> > > na.exclude that eludes me. I'm sure it is documented and I could,
> > > with time, unravel the mysteries but I really would prefer not to
> > > spend that time. I was hoping that someone else knew what the
> > > distinction was and how one would construct an extractor that would
> > > return the desired value and not be such an abominable kludge as the
> > > above.
> >
> > Yes. I've prepared and tested the following patch:
> >
> > Index: C:/Users/hennil/Desktop/R_packages/source/lme4_0.999375-8/R/lmer.R
> > ===================================================================
> > --- C:/Users/hennil/Desktop/R_packages/source/lme4_0.999375-8/R/lmer.R
> > (revision 99)
> > +++ C:/Users/hennil/Desktop/R_packages/source/lme4_0.999375-8/R/lmer.R
> > (working copy)
> > @@ -882,7 +882,14 @@
> > })
> >
> > setMethod("fitted", signature(object = "mer"),
> > - function(object, ...) object at eta)
> > + function(object, ...)
> > + {
> > + eta <- object at eta
> > + na.action <- attr(object at frame, "na.action")
> > + if(!is.null(na.action))
> > + napredict(na.action, eta)
> > + else eta
> > + })
> >
> > setMethod("formula", signature(x = "mer"),
> > function(x, ...)
> >
> >
> > After applying the patch, the distinction between `na.omit' and
> > `na.exclude' is easily shown:
> >
> >
> > > fm1 <- lmer(Reaction ~ Days + (Days|Subject), slps)
> > > fitted(fm1)
> >
> > [1] 253.4849 273.1808 292.8766 312.5725 332.2684 351.9642
> > [7] 371.6601 391.3559 411.0518 430.7476 213.1921 214.8561
> > [13] 216.5201 219.8481 221.5121 223.1761 224.8401 226.5041
> > [19] 228.1681 212.8984 217.8564 222.8143 227.7723 232.7302
> > [25] 237.6882 242.6462 247.6041 252.5621 257.5200 275.9163
> > [31] 281.5164 287.1166 298.3168 303.9170 309.5171 315.1173
> > [37] 320.7174 326.3175 273.5320 280.9433 288.3545 295.7658
> > [43] 303.1771 310.5884 317.9997 325.4110 332.8223 340.2336
> > [49] 260.3778 270.5809 280.7841 290.9872 301.1903 311.3934
> > [55] 321.5965 331.7996 342.0027 352.2058 268.1086 278.3689
> > [61] 288.6293 298.8896 309.1500 319.4103 329.6707 339.9310
> > [67] 350.1914 360.4517 244.2257 255.7619 267.2981 278.8342
> > [73] 290.3704 301.9066 313.4427 324.9789 336.5151 348.0512
> > [79] 251.2783 250.9591 250.6398 250.3205 250.0013 249.6820
> > [85] 249.3627 249.0435 248.7242 248.4050 283.8511 322.6857
> > [91] 342.1031 361.5204 380.9377 400.3550 419.7723 439.1897
> > [97] 458.6070 226.4061 238.0217 249.6373 261.2529 272.8685
> > [103] 284.4841 296.0997 307.7153 319.3310 330.9466 238.3395
> > [109] 255.4259 272.5123 289.5986 306.6850 323.7714 340.8578
> > [115] 357.9442 375.0305 392.1169 256.0056 263.4524 270.8991
> > [121] 278.3458 285.7925 293.2393 300.6860 308.1327 315.5795
> > [127] 323.0262 272.0277 286.0636 300.0995 314.1354 328.1714
> > [133] 342.2073 356.2432 370.2791 384.3150 398.3510 254.6441
> > [139] 265.9888 277.3335 288.6782 300.0229 311.3676 322.7123
> > [145] 334.0570 345.4017 356.7464 225.9405 241.2159 256.4914
> > [151] 271.7668 287.0422 302.3176 317.5930 332.8684 348.1438
> > [157] 363.4192 252.2315 261.7074 271.1833 280.6592 290.1351
> > [163] 299.6110 309.0870 318.5629 328.0388 337.5147 263.5955
> > [169] 275.3631 287.1308 298.8985 310.6662 322.4338 334.2015
> > [175] 345.9692 357.7368 369.5045
> > >
> > > fm2 <- lmer(Reaction ~ Days + (Days|Subject), slps,
> > + na.action = "na.exclude")
> > > fitted(fm2)
> >
> > [1] 253.4849 273.1808 292.8766 312.5725 332.2684 351.9642
> > [7] 371.6601 391.3559 411.0518 430.7476 213.1921 214.8561
> > [13] 216.5201 NA 219.8481 221.5121 223.1761 224.8401
> > [19] 226.5041 228.1681 212.8984 217.8564 222.8143 227.7723
> >
> > [25] 232.7302 237.6882 242.6462 247.6041 252.5621 257.5200
> > [31] 275.9163 281.5164 287.1166 NA 298.3168 303.9170
> > [37] 309.5171 315.1173 320.7174 326.3175 273.5320 280.9433
> > [43] 288.3545 295.7658 303.1771 310.5884 317.9997 325.4110
> >
> > [49] 332.8223 340.2336 260.3778 270.5809 280.7841 290.9872
> > [55] 301.1903 311.3934 321.5965 331.7996 342.0027 352.2058
> > [61] 268.1086 278.3689 288.6293 298.8896 309.1500 319.4103
> > [67] 329.6707 339.9310 350.1914 360.4517 244.2257 255.7619
> >
> > [73] 267.2981 278.8342 290.3704 301.9066 313.4427 324.9789
> > [79] 336.5151 348.0512 251.2783 250.9591 250.6398 250.3205
> > [85] 250.0013 249.6820 249.3627 249.0435 248.7242 248.4050
> > [91] 283.8511 NA 322.6857 342.1031 361.5204 380.9377
> >
> > [97] 400.3550 419.7723 439.1897 458.6070 226.4061 238.0217
> > [103] 249.6373 261.2529 272.8685 284.4841 296.0997 307.7153
> > [109] 319.3310 330.9466 238.3395 255.4259 272.5123 289.5986
> > [115] 306.6850 323.7714 340.8578 357.9442 375.0305 392.1169
> >
> > [121] 256.0056 263.4524 270.8991 278.3458 285.7925 293.2393
> > [127] 300.6860 308.1327 315.5795 323.0262 272.0277 286.0636
> > [133] 300.0995 314.1354 328.1714 342.2073 356.2432 370.2791
> > [139] 384.3150 398.3510 254.6441 265.9888 277.3335 288.6782
> >
> > [145] 300.0229 311.3676 322.7123 334.0570 345.4017 356.7464
> > [151] 225.9405 241.2159 256.4914 271.7668 287.0422 302.3176
> > [157] 317.5930 332.8684 348.1438 363.4192 252.2315 261.7074
> > [163] 271.1833 280.6592 290.1351 299.6110 309.0870 318.5629
> >
> > [169] 328.0388 337.5147 263.5955 275.3631 287.1308 298.8985
> > [175] 310.6662 322.4338 334.2015 345.9692 357.7368 369.5045
> > >
> >
> > I hope you find this useful.
> >
> >
> > Henric
> >
> >
> >
> >
> > >
> > >
> > > On Fri, Mar 14, 2008 at 9:43 AM, Douglas Bates <bates at stat.wisc.edu> wrote:
> > >> On Fri, Mar 14, 2008 at 7:31 AM, Andy Fugard <a.fugard at ed.ac.uk> wrote:
> > >>
> > >> > Is there any way to get fitted values out of lmer so that any NAs in the
> > >> > predictors/outcome pop out in the correct position?
> > >>
> > >> Yes.
> > >>
> > >>
> > >> > From what I can tell, the "fitted" function doesn't have access to
> > >> > missingness in the inputs.
> > >>
> > >> It does not, in its current implementation. I believe there is
> > >> information regarding the pattern of missingness in the model frame
> > >> (stored in the "frame" slot of the fitted model) and that can be used
> > >> to generate an appropriate fitted response.
> > >>
> > >> Does anyone know exactly how this is done for linear models? If so,
> > >> are you willing to provide a patch to the fitted method for mer
> > >> objects? I can put it on my ToDo list but when the first item on the
> > >> list is "finish the book" you almost never get to the second and
> > >> subsequent items.
> > >>
> > >
> >
> >
> > > _______________________________________________
> > > R-sig-mixed-models at r-project.org mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> > >
> >
>
More information about the R-sig-mixed-models
mailing list