[R] lme - problems with model
Douglas Bates
bates at stat.wisc.edu
Mon Feb 23 18:30:26 CET 2004
CG Pettersson <cg.pettersson at evp.slu.se> writes:
> Thanks a lot for the answer!
> Now, I only have the last one left - How do I get round it?
> I knew about the missing cells in the design, but didn´t know how lme
> would react on them.
> In this case, I can remove the water:temp term, but how can I be sure
> that this is the right thing to do?
Others may be able to come up with inventive ways of creating the
model matrix to do this but I would simply compare the additive model
to the cell means model.
> milk = read.table("/tmp/milk.txt", header = TRUE)
> milk$water = factor(milk$water)
> milk$temp = factor(milk$temp)
> milk$WaterTemp = factor(paste(milk$water, milk$temp, sep = '/'))
> 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
> xtabs(~ WaterTemp, milk)
WaterTemp
1/100 1/110 1/120 2/100 2/120 2/140 3/100 3/110 3/140
3 3 3 3 3 3 3 3 3
> library(lme4)
Loading required package: stats4
Loading required package: lattice
> summary(fm1 <- lme(maill6 ~ water + temp, milk, ~ 1 | rep))
Linear mixed-effects model fit by REML
Data: milk
AIC BIC logLik
6.155644 14.51182 4.922178
Random effects:
Groups Name Variance Std.Dev.
rep (Intercept) 0.021838 0.14778
Residual 0.017505 0.13231
Fixed effects: maill6 ~ water + temp
Estimate Std. Error DF t value Pr(>|t|)
(Intercept) 2.194667 0.103828 19 21.1376 1.162e-14 ***
water2 0.328000 0.068322 19 4.8008 0.0001243 ***
water3 -0.045333 0.068322 19 -0.6635 0.5149678
temp110 0.078000 0.072467 19 1.0764 0.2952465
temp120 0.321333 0.072467 19 4.4342 0.0002847 ***
temp140 0.350667 0.072467 19 4.8390 0.0001140 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Correlation of Fixed Effects:
(Intr) water2 water3 tmp110 tmp120
water2 -0.329
water3 -0.329 0.500
temp110 -0.310 0.236 0.000
temp120 -0.310 0.000 0.236 0.333
temp140 -0.155 -0.236 -0.236 0.333 0.333
Number of Observations: 27
Number of Groups: 3
> summary(fm2 <- lme(maill6 ~ WaterTemp, milk, ~ 1 | rep))
Linear mixed-effects model fit by REML
Data: milk
AIC BIC logLik
14.96052 24.75461 3.519738
Random effects:
Groups Name Variance Std.Dev.
rep (Intercept) 0.021862 0.14786
Residual 0.017286 0.13148
Fixed effects: maill6 ~ WaterTemp
Estimate Std. Error DF t value Pr(>|t|)
(Intercept) 2.1766667 0.1142339 16 19.0545 2.016e-12 ***
WaterTemp1/110 0.1400000 0.1073502 16 1.3041 0.21064
WaterTemp1/120 0.3133333 0.1073502 16 2.9188 0.01004 *
WaterTemp2/100 0.2833333 0.1073502 16 2.6393 0.01785 *
WaterTemp2/120 0.6933333 0.1073502 16 6.4586 7.897e-06 ***
WaterTemp2/140 0.7333333 0.1073502 16 6.8312 4.035e-06 ***
WaterTemp3/100 0.0533333 0.1073502 16 0.4968 0.62608
WaterTemp3/110 0.0066667 0.1073502 16 0.0621 0.95125
WaterTemp3/140 0.2866667 0.1073502 16 2.6704 0.01676 *
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Correlation of Fixed Effects:
(Intr) WT1/11 WT1/12 WT2/10 WT2/12 WT2/14 WT3/10 WT3/11
WtrTmp1/110 -0.470
WtrTmp1/120 -0.470 0.500
WtrTmp2/100 -0.470 0.500 0.500
WtrTmp2/120 -0.470 0.500 0.500 0.500
WtrTmp2/140 -0.470 0.500 0.500 0.500 0.500
WtrTmp3/100 -0.470 0.500 0.500 0.500 0.500 0.500
WtrTmp3/110 -0.470 0.500 0.500 0.500 0.500 0.500 0.500
WtrTmp3/140 -0.470 0.500 0.500 0.500 0.500 0.500 0.500 0.500
Number of Observations: 27
Number of Groups: 3
Both AIC and BIC, which are on the scale of "smaller is better",
indicate strong preference for the additive model. A likelihood ratio
test would not show significant improvement in the fit of the cell
means model relative to the additive model.
In general one should be cautious when using likelihood ratio tests on
the fixed-effects terms but in this case it is probably ok because the
estimates for the random effects are nearly identical for the two
models.
> Is the lm run without the random term enough for removing water:temp
> from the lme model, or do I have to do a PROC MIXED run with the
> random term to make that decision in a case like this?
I would use this analysis instead.
> Is it possible (for me) to understand why MIXED accepts the design
> but not lme? They ought to get the same sort of problems, or have I
> missed something?
Because of the way that model matrices are created in SAS, the
computational methods *must* detect rank deficiencies and compensate
for them. Whenever there is a categorical factor in the model the SAS
model matrix will be rank deficient. You may be aware that SAS always
uses a parameterization of a factor where the last level of the factor
is the "reference" level. That is not a choice - it is a necessary
consequence of the way in which the computation is carried out. It is
the detection of the rank deficiency and elimination of the column
where that is detected that causes SAS to eliminate the coefficient
for the last level in a factor.
The approach used in the S language is to use a set of k-1 "contrasts"
to generate the columns for terms involving a factor with k levels.
This automatically creates a full rank model matrix unless you have
missing cells. The lm code check for rank deficiency and compensates
for it. However we did not build that capability into lme (from the
nlme package or from the lme4 package) because there are two possible
explanations for the rank deficiency and in the one case we want to
circumvent it and in the other case we don't. It is difficult to
distinguish between those two cases so we don't even try.
> -------------------
> > 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/
> >
> CG Pettersson, MSci, PhD Stud.
> Swedish University of Agricultural Sciences
> Dep. of Ecology and Crop Production. Box 7043
> SE-750 07 Uppsala
--
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