[R-sig-ME] Desperately seeking help with simple lmer fit (lmer vs Proc Mixed)
Douglas Bates
bates at stat.wisc.edu
Sat Feb 21 00:01:28 CET 2009
On Wed, Feb 18, 2009 at 2:20 PM, Douglas Bates <bates at stat.wisc.edu> wrote:
> On Wed, Feb 18, 2009 at 1:56 PM, Douglas Bates <bates at stat.wisc.edu> wrote:
>> Thanks for the example.
>>
>> The estimates of the variance components are the result of an
>> optimization and it appears that there are two local optima for this
>> problem. I'm sorry to say that the optimum determined by SAS PROC
>> MIXED represents a better fit (a lower deviance) than the one
>> determined by lmer.
>>
>> If you plot the response versus run, say
>>
>> dotplot(reorder(Run, Y) ~ Y)
>>
>> you will see that there is one run (number 13) with unusually large Y,
>> which may have abnormal influence. I was able to reproduce the SAS
>> results by modifying the starting estimates. This is not optimal
>> because it means you need to know the correct answer before you can
>> calculate it.
>>
>>> (fm1 <- lmer(Y ~ (1|Run) + (1|Prep), pr, verbose = TRUE))
>> 0: 1536.6629: 0.666667 0.384900
>> 1: 1441.4487: 1.44044 1.01837
>> 2: 1421.4662: 2.35838 0.621655
>> 3: 1408.1762: 3.27848 0.229969
>> 4: 1401.9874: 4.04888 0.00000
>> 5: 1399.8678: 4.59188 0.00000
>> 6: 1398.9477: 5.12145 0.00000
>> 7: 1398.7722: 5.41633 0.00000
>> 8: 1398.7543: 5.53332 0.00000
>> 9: 1398.7539: 5.55421 0.00000
>> 10: 1398.7539: 5.55527 1.68436e-07
>> 11: 1398.7539: 5.55527 1.68436e-07
>> Linear mixed model fit by REML
>> Formula: Y ~ (1 | Run) + (1 | Prep)
>> Data: pr
>> AIC BIC logLik deviance REMLdev
>> 1407 1417 -699.4 1410 1399
>> Random effects:
>> Groups Name Variance Std.Dev.
>> Run (Intercept) 3.5859e+05 5.9882e+02
>> Prep (Intercept) 3.2965e-10 1.8156e-05
>> Residual 1.1619e+04 1.0779e+02
>> Number of obs: 108, groups: Run, 18; Prep, 6
>>
>> Fixed effects:
>> Estimate Std. Error t value
>> (Intercept) 28216.5 141.5 199.4
>>> (fm1a <- lmer(Y ~ (1|Run) + (1|Prep), pr, verbose = TRUE, start = c(4,4)))
>> 0: 1393.6144: 4.00000 4.00000
>> 1: 1393.5078: 3.74366 4.09884
>> 2: 1393.4819: 3.72800 4.32664
>> 3: 1393.4810: 3.73535 4.37362
>> 4: 1393.4810: 3.73373 4.37885
>> 5: 1393.4810: 3.73404 4.37877
>> 6: 1393.4810: 3.73404 4.37875
>> 7: 1393.4810: 3.73404 4.37875
>> Linear mixed model fit by REML
>> Formula: Y ~ (1 | Run) + (1 | Prep)
>> Data: pr
>> AIC BIC logLik deviance REMLdev
>> 1401 1412 -696.7 1406 1393
>> Random effects:
>> Groups Name Variance Std.Dev.
>> Run (Intercept) 162010 402.50
>> Prep (Intercept) 222784 472.00
>> Residual 11619 107.79
>> Number of obs: 108, groups: Run, 18; Prep, 6
>>
>> Fixed effects:
>> Estimate Std. Error t value
>> (Intercept) 28216 215 131.2
>>
>> If you remove Run 13 you do get results like those from SAS using lmer
>>
>>> (fm2 <- lmer(Y ~ (1|Run) + (1|Prep), pr, verbose = TRUE, subset = Run != 13))
>> 0: 1426.8110: 0.666667 0.396059
>> 1: 1342.0591: 1.47666 0.982502
>> 2: 1325.5670: 2.35364 0.501989
>> 3: 1311.1286: 3.60322 0.00000
>> 4: 1309.3047: 3.96698 0.00000
>> 5: 1307.8040: 4.56388 1.65053e-05
>> 6: 1307.5200: 4.88366 6.39529e-05
>> 7: 1307.4763: 5.04660 0.000170199
>> 8: 1307.4746: 5.08457 0.000341660
>> 9: 1307.4745: 5.08786 0.000627531
>> 10: 1307.4745: 5.08806 0.00113168
>> 11: 1307.4743: 5.09250 0.0198356
>> 12: 1306.7915: 5.55085 2.00076
>> 13: 1306.7202: 5.57778 2.29410
>> 14: 1304.5962: 3.26016 2.72067
>> 15: 1304.1219: 3.90630 4.98691
>> 16: 1303.9334: 3.70822 4.85924
>> 17: 1303.7839: 3.39994 4.44476
>> 18: 1303.7039: 3.78966 4.10571
>> 19: 1303.6539: 3.50693 3.67340
>> 20: 1303.6333: 3.63992 3.74762
>> 21: 1303.6285: 3.55042 3.87085
>> 22: 1303.6259: 3.59760 3.87143
>> 23: 1303.6258: 3.59381 3.86414
>> 24: 1303.6258: 3.59257 3.86456
>> 25: 1303.6258: 3.59308 3.86577
>> 26: 1303.6258: 3.59318 3.86504
>> 27: 1303.6258: 3.59312 3.86514
>> Linear mixed model fit by REML
>> Formula: Y ~ (1 | Run) + (1 | Prep)
>> Data: pr
>> Subset: Run != 13
>> AIC BIC logLik deviance REMLdev
>> 1312 1322 -651.8 1316 1304
>> Random effects:
>> Groups Name Variance Std.Dev.
>> Run (Intercept) 135857 368.59
>> Prep (Intercept) 157206 396.49
>> Residual 10523 102.58
>> Number of obs: 102, groups: Run, 17; Prep, 6
>>
>> Fixed effects:
>> Estimate Std. Error t value
>> (Intercept) 28161.8 185.5 151.8
>>
>>
>> This example may be motivation for me to try out a modification in the
>> optimization method that I have been contemplating.
>
> Well, the good news is that the modification was successful on this example
>
>> (fm1 <- lmer(Y ~ 1 + (1|Run) + (1|Prep), pr, verbose = TRUE))
> 0: 1536.6629: 0.816497 0.620403
> 1: 1405.8868: 1.66562 1.14860
> 2: 1399.0326: 1.71590 1.44355
> 3: 1394.7782: 1.84561 1.71319
> 4: 1393.9976: 1.79780 2.00855
> 5: 1393.9690: 2.09290 2.05797
> 6: 1393.5053: 1.96491 2.07162
> 7: 1393.4829: 1.92326 2.08972
> 8: 1393.4810: 1.93265 2.09030
> 9: 1393.4810: 1.93240 2.09189
> 10: 1393.4810: 1.93236 2.09254
> 11: 1393.4810: 1.93237 2.09255
> Linear mixed model fit by REML
> Formula: Y ~ 1 + (1 | Run) + (1 | Prep)
> Data: pr
> AIC BIC logLik deviance REMLdev
> 1401 1412 -696.7 1406 1393
> Random effects:
> Groups Name Variance Std.Dev.
> Run (Intercept) 162010 402.50
> Prep (Intercept) 222784 472.00
> Residual 11619 107.79
> Number of obs: 108, groups: Run, 18; Prep, 6
>
> Fixed effects:
> Estimate Std. Error t value
> (Intercept) 28216 215 131.2
>
> The bad news is that I made the modification in a development version
> of the lme4 package and that version won't be ready for prime time for
> a while. I still have a lot of work to do on the generalized linear
> mixed models.
I usually want to see a picture to help me understand what is going on
with a model so I created a contour plot of the profiled deviance as a
function of the two parameters in the new version of the optimization.
These parameters represent the square roots of the standard
deviations of the two sets of random effects relative to the residual
standard deviation. The conditional estimates of all the other
parameters, given these two, can be evaluated directly. The contours
are the joint confidence regions for this pair of parameters. You can
see that even on this square root scale the contours are somewhat
elongated to the right, and the value of 0 for the second variance
component (the prep component) is within confidence regions at about
85%. This is an indication that the variance component for Prep is
ill-determined.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: contours.pdf
Type: application/pdf
Size: 18787 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20090220/68619abe/attachment.pdf>
More information about the R-sig-mixed-models
mailing list