On Jun 27, 2013, at 3:47 PM, Shaun Jackman wrote:
> Hi Jean,
>
> contr.treatment(4) shows what the default contrast matrix looks like
> for a factor with 4 levels. What function do I use to create a
> contrast matrix to compare each level with the global mean (four
> comparisons in total), and produce a table similar to `summary.lm`?
>
I believe you asking for "contr.sum" although I think there might be some
differences between how it operates and what you are expressing as your
expectations.
> contrasts(ChickWeight$Diet) <- contr.sum(4)
> model <- lm(weight ~ Diet, ChickWeight)
> summary(model)
Call:
lm(formula = weight ~ Diet, data = ChickWeight)
Residuals:
Min 1Q Median 3Q Max
-103.95 -53.65 -13.64 40.38 230.05
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 125.869 2.986 42.150 < 2e-16 ***
Diet1 -23.223 4.454 -5.214 2.59e-07 ***
Diet2 -3.252 5.380 -0.604 0.54576
Diet3 17.081 5.380 3.175 0.00158 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 69.33 on 574 degrees of freedom
Multiple R-squared: 0.05348, Adjusted R-squared: 0.04853
F-statistic: 10.81 on 3 and 574 DF, p-value: 6.433e-07
> mean(ChickWeight$weight)
[1] 121.8183
> table(ChickWeight$Diet)
1 2 3 4
220 120 120 118
So in an unbalanced data situation, the Intercept is only approximately the
grand mean.
To see what you are requesting in the summary you can an offset from the mean
and use the Intercept suppression syntax:
> model <- lm(weight ~ Diet+0+offset(rep(mean(ChickWeight$weight),
> nrow(ChickWeight) )), ChickWeight)
> summary(model)
Call:
lm(formula = weight ~ Diet + 0 + offset(rep(mean(ChickWeight$weight),
nrow(ChickWeight))), data = ChickWeight)
Residuals:
Min 1Q Median 3Q Max
-103.95 -53.65 -13.64 40.38 230.05
Coefficients:
Estimate Std. Error t value Pr(>|t|)
Diet1 -19.1729 4.6740 -4.102 4.69e-05 ***
Diet2 0.7983 6.3286 0.126 0.899660
Diet3 21.1317 6.3286 3.339 0.000895 ***
Diet4 13.4444 6.3820 2.107 0.035584 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 69.33 on 574 degrees of freedom
Multiple R-squared: 0.7599, Adjusted R-squared: 0.7583
F-statistic: 454.3 on 4 and 574 DF, p-value: < 2.2e-16
Notice this does estimate waht you requested, but I think it is more due to the
use of an offset than to the choice of contrasts.
> with(ChickWeight, tapply(weight, Diet, function(categ) mean(categ)-
> mean(weight) ) )
1 2 3 4
-19.1728846 0.7983276 21.1316609 13.4443728
I'm very worried this might be inferentially suspect, since the degrees of
freedom and the anava F statistic are different than the usual methods.
--
David.
> Thanks,
> Shaun
>
>
> On 26 June 2013 05:50, Adams, Jean <jvad...@usgs.gov> wrote:
>> Shaun,
>>
>> See the help on contrasts ...
>> ?contr.treatment
>>
>> Jean
>>
>>
>> On Tue, Jun 25, 2013 at 7:07 PM, Shaun Jackman <sjack...@gmail.com> wrote:
>>>
>>> Hi,
>>>
>>> I've used `lm` to create a linear model of a continuous variable
>>> against a factor variable with four levels using an example R data set
>>> (see below). By default, it uses a treatment contrast matrix that
>>> compares each level of the factor variable with the first reference
>>> level (three comparisons in total). I'd like to compare each level
>>> with the global mean (four comparisons in total), and produce a table
>>> similar to `summary.lm`. How do I go about this?
>>>
>>> ```r
>>> model <- lm(weight ~ Diet, ChickWeight)
>>> summary(model)
>>> ```
>>>
>>> Thanks,
>>> Shaun
>>>
>>> ______________________________________________
>>> R-help@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.
>>
>>
>
> ______________________________________________
> R-help@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.
David Winsemius
Alameda, CA, USA
______________________________________________
R-help@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.
