[R] [OFF] stepwise using REML???

[Ronaldo Reis Jr.]

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