Re: [R] lme - problems with model
As Spencer Graves suggested, I tried this with continuous variables. Seems to work ok: lme(maill6 ~ water * temp , random= ~1|rep, data = milk) Linear mixed-effects model fit by REML Data: milk Log-restricted-likelihood: -10.57237 Fixed: maill6 ~ water * temp (Intercept)water temp water:temp -1.107227891 0.928965420 0.032507653 -0.008792517 Random effects: Formula: ~1 | rep (Intercept) Residual StdDev: 0.1358565 0.2189339 Number of Observations: 27 Number of Groups: 3 For the smaller model, I get: lme(maill6 ~ water + temp , random= ~1|rep, data = milk) Linear mixed-effects model fit by REML Data: milk Log-restricted-likelihood: -8.068963 Fixed: maill6 ~ water + temp (Intercept) watertemp 1.1708 -0.05819444 0.01212500 Random effects: Formula: ~1 | rep (Intercept) Residual StdDev: 0.1328748 0.2348303 Number of Observations: 27 Number of Groups: 3 Cheers, Dick *** Richard P. Beyer, Ph.D. University of Washington Tel.:(206) 616 7378 Env. Occ. Health Sci. , Box 354695 Fax: (206) 685 4696 4225 Roosevelt Way NE, # 100 Seattle, WA 98105-6099 http://depts.washington.edu/ceeh/ServiceCores/FC5/FC5.html *** Message: 14 Date: Mon, 23 Feb 2004 07:41:07 -0800 From: Spencer Graves [EMAIL PROTECTED] Subject: Re: [R] lme - problems with model To: CG Pettersson [EMAIL PROTECTED] Cc: Douglas Bates [EMAIL PROTECTED], [EMAIL PROTECTED] Message-ID: [EMAIL PROTECTED] Content-Type: text/plain; charset=ISO-8859-1; format=flowed If you want to try to get the same answers as PROC MIXED, I suggest you try to figure out how SAS codes interactions and which ones it retains. Then you can try code those manually and include them as separate explanatory variables, e.g., I((water==2)(temp==110)). You could work this out in lm then try the result on lme. An alternative would be to convert temp from a factor to a continuous variable. I would make plots of the response variables vs. temp with different lines and symbols for water and rep to make sure I had something that was mostly linear in some transformation of temp. hope this helps. spencer graves CG Pettersson wrote: 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? 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? 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? /CG --- CG Pettersson [EMAIL PROTECTED] 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. Errort value Pr(|t|) (Intercept) 2.1767 0.1142339 19.0544730 2.218661e-13 water2 0.2833 0.1615511 1.7538308 9.647013e-02 water3 0.0533 0.1615511 0.3301329 7.451108e-01 temp110 0.1400 0.1615511 0.8665987 3.975669e-01 temp120 0.3133 0.1615511 1.9395305 6.827304e-02 temp140 0.2333 0.1615511 1.4443312 1.658280e-01 water3:temp110 -0.1867 0.2284678 -0.8170371 4.245898e-01 water2:temp120 0.0967 0.2284678 0.4231085 6.772282e-01 water2:temp140 0.2167 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
Re: [R] lme - problems with model
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? 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? 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? /CG --- CG Pettersson [EMAIL PROTECTED] 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. Errort value Pr(|t|) (Intercept) 2.1767 0.1142339 19.0544730 2.218661e-13 water2 0.2833 0.1615511 1.7538308 9.647013e-02 water3 0.0533 0.1615511 0.3301329 7.451108e-01 temp110 0.1400 0.1615511 0.8665987 3.975669e-01 temp120 0.3133 0.1615511 1.9395305 6.827304e-02 temp140 0.2333 0.1615511 1.4443312 1.658280e-01 water3:temp110 -0.1867 0.2284678 -0.8170371 4.245898e-01 water2:temp120 0.0967 0.2284678 0.4231085 6.772282e-01 water2:temp140 0.2167 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.09.6 2 1 100 2 2.19 2.20 2.27 11.09.3 11.0 3 1 100 3 2.13 2.20 2.41 10.17.09.6 4 1 110 1 2.13 2.34 2.41 11.0 10.59.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.429.4 10.79.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.79.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.29.2 12 2 100 3 2.80 2.63 3.338.39.77.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.99.09.4 16 2 140 1 2.56 2.84 3.10 10.7 11.29.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.759.29.69.6 19 3 100 1 2.60 2.24 2.32 10.88.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.49.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.819.39.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 ~
Re: [R] lme - problems with model
If you want to try to get the same answers as PROC MIXED, I suggest you try to figure out how SAS codes interactions and which ones it retains. Then you can try code those manually and include them as separate explanatory variables, e.g., I((water==2)(temp==110)). You could work this out in lm then try the result on lme. An alternative would be to convert temp from a factor to a continuous variable. I would make plots of the response variables vs. temp with different lines and symbols for water and rep to make sure I had something that was mostly linear in some transformation of temp. hope this helps. spencer graves CG Pettersson wrote: 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? 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? 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? /CG --- CG Pettersson [EMAIL PROTECTED] 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. Errort value Pr(|t|) (Intercept) 2.1767 0.1142339 19.0544730 2.218661e-13 water2 0.2833 0.1615511 1.7538308 9.647013e-02 water3 0.0533 0.1615511 0.3301329 7.451108e-01 temp110 0.1400 0.1615511 0.8665987 3.975669e-01 temp120 0.3133 0.1615511 1.9395305 6.827304e-02 temp140 0.2333 0.1615511 1.4443312 1.658280e-01 water3:temp110 -0.1867 0.2284678 -0.8170371 4.245898e-01 water2:temp120 0.0967 0.2284678 0.4231085 6.772282e-01 water2:temp140 0.2167 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.09.6 2 1 100 2 2.19 2.20 2.27 11.09.3 11.0 3 1 100 3 2.13 2.20 2.41 10.17.09.6 4 1 110 1 2.13 2.34 2.41 11.0 10.59.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.429.4 10.79.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.79.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.29.2 12 2 100 3 2.80 2.63 3.338.39.77.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.99.09.4 16 2 140 1 2.56 2.84 3.10 10.7 11.29.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.759.29.69.6 19 3 100 1 2.60 2.24 2.32 10.88.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
Re: [R] lme - problems with model
CG Pettersson [EMAIL PROTECTED] 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 NameVariance 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 NameVariance 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.177 0.1142339 16 19.0545 2.016e-12 *** WaterTemp1/110 0.140 0.1073502 16 1.3041 0.21064 WaterTemp1/120 0.313 0.1073502 16 2.9188 0.01004 * WaterTemp2/100 0.283 0.1073502 16 2.6393 0.01785 * WaterTemp2/120 0.693 0.1073502 16 6.4586 7.897e-06 *** WaterTemp2/140 0.733 0.1073502 16 6.8312 4.035e-06 *** WaterTemp3/100 0.053 0.1073502 16 0.4968 0.62608 WaterTemp3/110 0.007 0.1073502 16 0.0621 0.95125 WaterTemp3/140 0.287 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.
Re: [R] lme - problems with model
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 [EMAIL PROTECTED] 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 NameVariance 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 NameVariance 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.177 0.1142339 16 19.0545 2.016e-12 *** WaterTemp1/110 0.140 0.1073502 16 1.3041 0.21064 WaterTemp1/120 0.313 0.1073502 16 2.9188 0.01004 * WaterTemp2/100 0.283 0.1073502 16 2.6393 0.01785 * WaterTemp2/120 0.693 0.1073502 16 6.4586 7.897e-06 *** WaterTemp2/140 0.733 0.1073502 16 6.8312 4.035e-06 *** WaterTemp3/100 0.053 0.1073502 16 0.4968 0.62608 WaterTemp3/110 0.007 0.1073502 16 0.0621 0.95125 WaterTemp3/140 0.287 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)
