[R] lme - problems with model

Douglas Bates bates at stat.wisc.edu
Mon Feb 23 14:26:39 CET 2004


CG Pettersson <cg.pettersson at evp.slu.se> writes:

> Hello all!
> 
> I´m working with some training datasets in a SAS-based course, trying
> to do the same things in lme that I do in PROC MIXED. 
> 
> Why don´t lme do an analysis on this dataset when I use the model
> water*temp?
> The trouble comes from the water:temp term, as it works with
> water+temp.
> The data are, indeed, assymetric but lm accepts the water:temp term
> giving results in the F-test near what PROC MIXED produces. MIXED
> accepts the model.
> 
> The water:temp term can be removed from the model according to the
> F-test in SAS (and to the lm model without any random term). Doing so
> in both MIXED and lme gives reasonably similar results for both
> systems.
> 
> What do the error message mean, and how can I get around this?

Because of missing cells in the design

> xtabs(~water + temp, milk)
     temp
water 100 110 120 140
    1 3   3   3   0  
    2 3   0   3   3  
    3 3   3   0   3  

the model matrix for the fixed effects is rank deficient.  In lm the
rank deficiency is detected and appropriate adjustments made

> coef(summary(lm(maill6 ~ water * temp, milk)))
                  Estimate Std. Error    t value     Pr(>|t|)
(Intercept)     2.17666667  0.1142339 19.0544730 2.218661e-13
water2          0.28333333  0.1615511  1.7538308 9.647013e-02
water3          0.05333333  0.1615511  0.3301329 7.451108e-01
temp110         0.14000000  0.1615511  0.8665987 3.975669e-01
temp120         0.31333333  0.1615511  1.9395305 6.827304e-02
temp140         0.23333333  0.1615511  1.4443312 1.658280e-01
water3:temp110 -0.18666667  0.2284678 -0.8170371 4.245898e-01
water2:temp120  0.09666667  0.2284678  0.4231085 6.772282e-01
water2:temp140  0.21666667  0.2284678  0.9483467 3.555125e-01

Notice that you would expect 6 degrees of freedom for the interaction
term but only three coefficients are estimated.

In lme it is much more difficult to compensate for such rank
deficiencies because they could be systematic, like this, or they
could be due to relative precision parameters approaching zero during
the iterations.  Because of this we just report the error (although
admittedly we could be a bit more explicit about the nature of the
problem - we are reporting the symptom that we detect, not the
probable cause).


> The dataset:
> > milk
>    water temp rep maill4 maill6 maill8 taste4 taste6 taste8
> 1      1  100   1   2.90   2.13   2.39   10.1   10.0    9.6
> 2      1  100   2   2.19   2.20   2.27   11.0    9.3   11.0
> 3      1  100   3   2.13   2.20   2.41   10.1    7.0    9.6
> 4      1  110   1   2.13   2.34   2.41   11.0   10.5    9.8
> 5      1  110   2   2.32   2.27   2.25   11.0   11.3   11.2
> 6      1  110   3   2.13   2.34   2.42    9.4   10.7    9.0
> 7      1  120   1   2.00   2.49   2.71   11.1   11.2   11.4
> 8      1  120   2   2.41   2.49   2.46   11.6   11.7    9.6
> 9      1  120   3   2.22   2.49   2.73   10.7   10.3   10.2
> 10     2  100   1   2.13   2.41   2.49   11.1   10.8   11.2
> 11     2  100   2   2.49   2.34   2.53   11.1   11.2    9.2
> 12     2  100   3   2.80   2.63   3.33    8.3    9.7    7.8
> 13     2  120   1   2.38   2.85   2.06   11.9   11.2   11.2
> 14     2  120   2   2.61   2.70   2.70   11.7   10.8   11.0
> 15     2  120   3   2.77   3.06   3.25   10.9    9.0    9.4
> 16     2  140   1   2.56   2.84   3.10   10.7   11.2    9.8
> 17     2  140   2   2.63   2.61   2.81   10.8   11.0   11.6
> 18     2  140   3   2.99   3.28   3.75    9.2    9.6    9.6
> 19     3  100   1   2.60   2.24   2.32   10.8    8.4   10.8
> 20     3  100   2   2.06   2.11   2.20   11.0   11.2   11.8
> 21     3  100   3   1.98   2.34   2.80   10.3   10.2   10.6
> 22     3  110   1   1.91   2.06   2.29   11.0   11.4    9.4
> 23     3  110   2   1.98   1.98   2.15   10.0   11.8   10.6
> 24     3  110   3   1.98   2.51   2.81    9.3    9.2   10.2
> 25     3  140   1   2.27   2.42   2.72   10.8   11.6   12.0
> 26     3  140   2   2.27   2.20   2.41   11.2   11.0   11.4
> 27     3  140   3   2.20   2.77   3.06   10.5   10.2   10.0
> 
> The failing model:
> > lme(maill6 ~ water * temp  , random= ~1|rep, data = milk)
> Error in MEEM(object, conLin, control$niterEM) : 
>         Singularity in backsolve at level 0, block 1
> 
> The smaller (working) model:
> > lme(maill6 ~ water + temp  , random= ~1|rep, data = milk)
> Linear mixed-effects model fit by REML
>   Data: milk 
>   Log-restricted-likelihood: 4.922178
>   Fixed: maill6 ~ water + temp 
> (Intercept)      water2      water3     temp110     temp120    
> temp140 
>  2.19466667  0.32800000 -0.04533333  0.07800000  0.32133333 
> 0.35066667 
> 
> Random effects:
>  Formula: ~1 | rep
>         (Intercept)  Residual
> StdDev:   0.1477760 0.1323057
> 
> Number of Observations: 27
> Number of Groups: 3 
> > 
> 
> 
> 
> 
> CG Pettersson, MSci, PhD Stud.
> Swedish University of Agricultural Sciences
> Dep. of Ecology and Crop Production. Box 7043
> SE-750 07 Uppsala
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

-- 
Douglas Bates                            bates at stat.wisc.edu
Statistics Department                    608/262-2598
University of Wisconsin - Madison        http://www.stat.wisc.edu/~bates/




More information about the R-help mailing list