From john.maindonald at anu.edu.au  Wed May 11 01:06:33 2011
From: john.maindonald at anu.edu.au (John Maindonald)
Date: Wed, 11 May 2011 09:06:33 +1000
Subject: [Rd] anova.lm fails with test="Cp"
In-Reply-To: <AA81EAC1-22CD-4550-BFCA-C04F42F8FB3A@gmail.com>
References: <0E28A214-9F72-45A2-AAA2-E1D7EC46C6B1@anu.edu.au>
	<AA81EAC1-22CD-4550-BFCA-C04F42F8FB3A@gmail.com>
Message-ID: <13D70C43-99FF-45E6-B4EE-6D9951DE00DB@anu.edu.au>

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
> 



[Rd] anova.lm fails with test="Cp"
John Maindonald