[R] something confusing about stepAIC
Steve Cumming
stevec at berl.ab.ca
Tue Apr 2 00:44:47 CEST 2002
Folks, I'm using stepAIC(MASS) to do some automated, exploratory, model
selection for binomial and Poisson glm models in R 1.3. Because I wanted to
experiment with the small-sample correction AICc, I dug around in the code
for the functions
glm.fit
stepAIC
dropterm.glm
addterm.glm
extractAIC.glm
and came across something I just don't understand.
stepAIC() passes dropterm.glm() a model object. dropterm.glm() then fits a
number of submodels, computing for each some measure DeltaFit of the
relative change in goodness of fit. It then returns to stepAIC() with some
indication of which submodel is best. addterm.glm behaves similarly. The
problem is, both functions use the submodel deviances to compute the
DeltaFits, not the submodel AICs, even though these are available for at
least the poisson, binomial and negbin families. The step() function works
much the same way.
Is this behaviour correct? It is not obvious to me that the deviance-based
criteria will produce the same sequence of steps as the AIC criteria, in
general. Surely that depends on the exact deviance statistic used, and/or
the error model? It seems to me that the appropriate change to the inner
loop in dropterm.glm() and addterm.glm() would not break anything, and would
make stepAIC() safe to use on models returned by glm.nb().
Am I missing something?
I note that v.1.4 has some new AIC methods. Should I migrate?
Any comments or suggestions would be much appreciated.
Steve Cumming
Boreal Ecosystems Research Ltd.
P.S. I don't know how to express my gratitude and amazement to the people
responsible for R, so I'll just say "thanks".
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list