Probably the most obvious difference is that ML estimates of variance components are biased. The bias is equivalent to estimating the variance of a population by dividing the sum of squared deviations from the sample mean of N observations by N rather than N-1.
To check for other differences, I just did library(mle4) in R 1.9.1 and ran the example in help("lme") then ran the same lme call with method="ML". The results (see below) show that the parameter estimates for the fixed effects are the same for ML as for REML, but the estimates of standard errors of those parameters seems to perpetuate the bias corrected by REML.
If stepAIC works with ML but not with REML, I would use stepAIC with ML then reestimate the final model with REML. However, I would also do some simulations to estimate the false positive (Type I) error rate. AIC includes a penalty for the complexity of the model but NOT for the number of alternatives considered, and if you present enough random alternatives to the procedure, it will find some that are statistically significant. The nlme package contains a function "simulate.lme" to facilitate this.
hope this helps. spencer graves
> library(lme4)
> data(bdf)
> fm <- lme(langPOST ~ IQ.ver.cen + avg.IQ.ver.cen, data = bdf,
+ random = ~ IQ.ver.cen | schoolNR)
> summary(fm)
Linear mixed-effects model fit by REML
Fixed: langPOST ~ IQ.ver.cen + avg.IQ.ver.cen
Data: bdf
AIC BIC logLik
15231.87 15272.02 -7608.937Random effects:
Groups Name Variance Std.Dev. Corr schoolNR (Intercept) 8.07563 2.84177 IQ.ver.cen 0.20801 0.45608 -0.642
Residual 41.34968 6.43037 # of obs: 2287, groups: schoolNR, 131
Fixed effects:
Estimate Std. Error DF t value Pr(>|t|) (Intercept) 4.0750e+01 2.8808e-01 2284 141.452 < 2.2e-16 ***
IQ.ver.cen 2.4598e+00 8.3638e-02 2284 29.410 < 2.2e-16 ***
avg.IQ.ver.cen 1.4089e+00 3.2374e-01 2284 4.352 1.409e-05 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Correlation of Fixed Effects:
(Intr) IQ.vr.
IQ.ver.cen -0.274 avg.IQ.vr.c 0.029 -0.213
> fm.ml <- lme(langPOST ~ IQ.ver.cen + avg.IQ.ver.cen, data = bdf,
+ random = ~ IQ.ver.cen | schoolNR, method="ML")
> summary(fm.ml)
Linear mixed-effects model fit by maximum likelihood
Fixed: langPOST ~ IQ.ver.cen + avg.IQ.ver.cen
Data: bdf
AIC BIC logLik
15227.53 15267.68 -7606.767
Random effects:
Groups Name Variance Std.Dev. Corr schoolNR (Intercept) 7.91904 2.81408 IQ.ver.cen 0.20006 0.44728 -0.652
Residual 41.35055 6.43044 # of obs: 2287, groups: schoolNR, 131
Fixed effects:
Estimate Std. Error DF t value Pr(>|t|) (Intercept) 4.0750e+01 2.8610e-01 2284 142.4342 < 2.2e-16 ***
IQ.ver.cen 2.4589e+00 8.3237e-02 2284 29.5406 < 2.2e-16 ***
avg.IQ.ver.cen 1.4052e+00 3.2168e-01 2284 4.3683 1.308e-05 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Correlation of Fixed Effects:
(Intr) IQ.vr.
IQ.ver.cen -0.274 avg.IQ.vr.c 0.028 -0.214
#######################
Klaus Thul wrote:
Dear all,
I am planning to use nlme library for analysis of experiments in semiconductor industry. Currently I am using "lm" but plan to move to "lme" to handle within wafer / wafer-to-wafer and lot-to-lot variation correctly.
So far everything is working well, but I have a fundamentel question:
NLME offers "maximum likelihood" and "restricted maximum likelihood" for parameter estimation. ML has the advantage, that likelihood ratios can be computed even with changes in model structure. In addition, ML works with the stepAIC function from MASS-library which I am currently using for model-building.
I am wondering, why REML is the default setting in NLME and therefore somehow preferred by the authors. What is the main reason to use REML?
Maybe I am lacking here statistical knowledge. Any hint, including reference to literature would be very helpful.
Best regards, Klaus Thul
______________________________________________
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
-- Spencer Graves, PhD, Senior Development Engineer O: (408)938-4420; mobile: (408)655-4567
______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
