[R-sig-ME] Fitted values in lmer
Douglas Bates
bates at stat.wisc.edu
Fri Mar 14 22:38:53 CET 2008
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