[R] lme - problems with model
Spencer Graves
spencer.graves at pdf.com
Mon Feb 23 18:55:41 CET 2004
Doug's "xtabs" suggests to me that the following might be
estimable, with data.$Temp <- as.numeric(as.character(data.$temp))
water*(Temp+I(Temp^2))
It looks like it should be estimable in "lm", and depending on the
noise model, it should also be estimable in "lme". ???
hope this helps. spencer graves
Douglas Bates wrote:
>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
>>
>>
>
>
>
More information about the R-help
mailing list