anova.lm fails with test="Cp"
For unknown sigma^2, a version that is a modification of AIC may be preferred, i.e. n log(RSS/n) + 2p - n I notice that this is what is given in Maindonald and Braun (2010) Data Analysis & Graphics Using R, 3rd edition. Cf: Venables and Ripley, MASS, 4th edn, p.174. V&R do however stop short of actually saying that Cp should be modified in the same way as AIC when sigma^2 has to be estimated. Better still, perhaps, give the AIC statistic. This would make the output consistent with dropterm(), drop1() and add1(). Or if Cp is to stay, allow AIC as a a further "test". John Maindonald email: john.maindonald at anu.edu.au phone : +61 2 (6125)3473 fax : +61 2(6125)5549 Centre for Mathematics & Its Applications, Room 1194, John Dedman Mathematical Sciences Building (Building 27) Australian National University, Canberra ACT 0200. http://www.maths.anu.edu.au/~johnm
On 08/05/2011, at 6:15 PM, peter dalgaard wrote:
On May 8, 2011, at 09:25 , John Maindonald wrote:
Here is an example, modified from the help page to use test="Cp": --------------------------------------------------------------------------------
fit0 <- lm(sr ~ 1, data = LifeCycleSavings) fit1 <- update(fit0, . ~ . + pop15) fit2 <- update(fit1, . ~ . + pop75) anova(fit0, fit1, fit2, test="Cp")
Error in `[.data.frame`(table, , "Resid. Dev") : undefined columns selected
Yes, the "Resid. Dev" column is only there in analysis of deviance tables. For the lm() case, it looks like you should have "RSS".
This has probably been there "forever". Just goes to show how often people use these things...
Also, now that I'm looking at it, are we calculating it correctly in any case? We have
cbind(table, Cp = table[, "Resid. Dev"] + 2 * scale *
(n - table[, "Resid. Df"]))
whereas all the references I can find have Cp=RSS/MS-N+2P, so the above would actually be scale*Cp+N.
--
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com