Matthieu Stigler
2009-Nov-28 15:14 UTC
[R] Include manually an intercept in lm without breaking it?
Hi Say I want to add manually an intercept in the function lm. Even if almost all results will be identical, few stats are different as DF counting will be different as intercept will not be included in "automatic" case, while it will be in "manual" case. See: ###usual lm on freeny fr<-lm(freeny.y~freeny.x) ###manual lm on freeny man<-cbind(1,freeny.x) colnames(man)<-c("const",colnames(freeny.x)) fr_man<-lm(freeny.y~man-1) ###coef are the same cbind(coef(fr), coef(fr_man)) ###but summary output is different (but should be the same!). #Difference comes from fact that in the "automatic case", DF are 4 (intercept not included) summary(fr) summary(fr_man) ###Workaround: attr(fr_man$terms, "intercept") <- 1 ##so now: summary(fr) summary(fr_man) ###but have negative effect that intercept is used twice in model.matrix, see: model.matrix(fr_man) So I could not find a good way to add manually an intercept and preserving the right output... any idea? Thanks a lot!! Matthieu Stigler
Duncan Murdoch
2009-Nov-28 15:31 UTC
[R] Include manually an intercept in lm without breaking it?
On 28/11/2009 10:14 AM, Matthieu Stigler wrote:> Hi > > Say I want to add manually an intercept in the function lm. Even if > almost all results will be identical, few stats are different as DF > counting will be different as intercept will not be included in > "automatic" case, while it will be in "manual" case. See: > > ###usual lm on freeny > fr<-lm(freeny.y~freeny.x) > ###manual lm on freeny > man<-cbind(1,freeny.x) > colnames(man)<-c("const",colnames(freeny.x)) > fr_man<-lm(freeny.y~man-1) > > ###coef are the same > cbind(coef(fr), coef(fr_man)) > > ###but summary output is different (but should be the same!).Why do you say it should be the same? The residual sum of squares needs to be calculated against some reference model. When you have the intercept included automatically, that's an indication that you want the reference model to include the intercept. When you put "-1" in your model spec, that's an indication that you don't want to include it, you want to compare against 0. If you want to change the conventions, you'll need to write your own summary.lm function. But it would be easier to just use the standard convention. Duncan Murdoch> #Difference comes from fact that in the "automatic case", DF are 4 > (intercept not included) > summary(fr) > summary(fr_man) > > ###Workaround: > attr(fr_man$terms, "intercept") <- 1 > > ##so now: > summary(fr) > summary(fr_man) > > > ###but have negative effect that intercept is used twice in > model.matrix, see: > model.matrix(fr_man) > > So I could not find a good way to add manually an intercept and > preserving the right output... any idea? > > Thanks a lot!! > > Matthieu Stigler > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code.