[R] Include manually an intercept in lm without breaking it?

[Matthieu Stigler]

Thanks for your answer!

Actually, the reason to include it manually was that the function I 
wanted to create should have args like:
type=c("const", "trend", "both", "none") and so was easier to create a 
matrix and simply put
lm(y~-1+datamat)

But following your advices I think rather change that approach and will use:

reg<-if(any(c("const", "both") %in% type)){
  lm(y ~ 1 + . -const, data = datamat) #there is const, use internal one 
and remove from datamat
}else {
  lm(y ~ -1 + ., data = datamat)}#there is no const, remove default one

So this should avoid the problem mentionned.

Thanks again!

Matthieu


Duncan Murdoch a écrit :
> 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.
>