[OFF] stepwise using REML???

Em Dom 22 Jun 2003 11:18, Douglas Bates escreveu:

For nested design it may be very dangerous due the difference in
variance structure, mainly in a splitplot design. ML make
significative variables that REML dont make.

It would be good to quote an example that shows this.  I'm not sure
that this occurs in general.

Look this example:

Using stepwise with a ML estimation:
--------------------------------------
m0.ml <- lme(response~1,random=~1|plot1/plot2,method="ML")
mfull.ml <- update(m0.ml,.~.+v1*v2+v1*v3)

stepAIC(mfull.ml)

Start:  AIC= 250.23 
 response ~ v1 + v2 + v3 + v1:v2 + v1:v3 

                   Df    AIC
- v1:v3		   11 249.82
<none>                250.23
- v1:v2		    2 253.65

...

Linear mixed-effects model fit by maximum likelihood
  Data: NULL 
  Log-likelihood: -112.9370
  Fixed: response ~ v1 + v2 + v1:v2 
            (Intercept)                   v1       		   v2l2 
             12.5936495              -0.3327049             -15.9920774 
                   v2l3 	        v1:v2l2   		v1:v2l3 
            -12.5285727               0.3750894               0.3014936 

Random effects:
 Formula: ~1 | plot1
        (Intercept)
StdDev: 0.004214203

 Formula: ~1 | plot2 %in% plot1
        (Intercept)  Residual
StdDev:   0.3051747 0.5556307

Number of Observations: 124
Number of Groups: 
            plot1  plot2 %in% plot1 
                6                17 
--------------------------------------
in this case the selected model is v1*v2, but in this case it use the same 
denominator DF for all variables, and it is not true.

In anova made by REML:
--------------------------------------
m0 <- lme(response~1,random=~1|plot1/plot2)
mfull <- update(m0,.~.+v1*v2+v1*v3)
anova(mfull)
                 numDF denDF   F-value p-value
(Intercept)          1    85 126.08414  <.0001
v1                   1     8   1.86229  0.2095
v2                   2     3   1.90189  0.2928
v3                  11    85   1.58872  0.1167
v1:v2                2     8   3.53897  0.0792
v1:v3               11    85   1.71976  0.0825
---------------------------------------

In this case all variables are not significative.

I read an article that is made a stepwise procedure using GENSTAT.

from article:
"Terms were dropped from a model in a stepwise procedure by assessing the
change in deviance between the full model and the submodel."

All are made using REML.

It is possible?! I dont know GENSTAT.

You would need to be more specific about how the comparisons are made.
I assume that you plan to keep the random effects structure constant
and compare two nested models that differ only in the fixed effects
terms.  I can think of four ways of doing this:

1) Use the F-test obtained by fitting the full model and conditioning
on the estimates of the random effects parameters.  This is what the
anova function applied to an model fit by lme gives.

2) Fit both models and compare the values of the REML criterion in a
likelihood ratio test.

3) Fit both models by REML and compare the values of the
log-likelihood (i.e. the ML criterion) in a likelihood ratio test.
You can obtain that value with logLik(fm, REML=FALSE) if fm is your
fitted model.

4)Fit both models and evaluate the REML criterion for the full model
at the two sets of estimates.  Compare these values in a likelihood
ratio test.

I feel that 1) is appropriate, 2) is inappropriate, 3) may be
appropriate and 4) looks interesting.  4) is based on recent work by
Greg Reinsel.

In some simulations reported in chapter 3 of Pinheiro and Bates (2000)
3) fared badly compared to 1).

Thanks for all
 ______________________________________________

R-help at stat.math.ethz.ch mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help

Nothing is but what is not.
--
|   // | \\   [***********************************]
|> ( ?   ? )  [Ronaldo Reis J?nior                ]
|      V      [UFV/DBA-Entomologia                ]
|>  /     \   [36571-000 Vi?osa - MG              ]
|  /(.''`.)\  [Fone: 31-3899-2532                 ]
|>/(: :'  :)\ [chrysopa at insecta.ufv.br            ]
|/ (`. `'` ) \[ICQ#: 5692561 | LinuxUser#: 205366 ]
|>  ( `-  )   [***********************************]
|>> _/   \_Powered by GNU/Debian Woody/Sarge