1
$\begingroup$

Using the ChickWeight data provided by R,

> head(ChickWeight)
  weight Time Chick Diet
1     42    0     1    1
2     51    2     1    1
3     59    4     1    1
4     64    6     1    1
5     76    8     1    1
6     93   10     1    1

using linear regression methods I am trying to answer the question "Is there a significant difference in weight gain per day between the chicks on diets 2 compared to 4"? In the regression model, the gradient will be the chick weight gain per day, so I want to compare the gradient of chicks on diet 2 to those on diet 4.

Using sub-setting I get the gradient of Diet 4 is not significantly more with a p-value of 0.0822.

> ChickWeight_D2v4 = subset(ChickWeight,Diet==2|Diet==4)
> ChickWeight.lm.D2v4 <- lm(weight~Time*Diet,data=ChickWeight_D2v4)
> summary(ChickWeight.lm.D2v4)

Call:
lm(formula = weight ~ Time * Diet, data = ChickWeight_D2v4)

Residuals:
     Min       1Q   Median       3Q      Max 
-135.425  -10.900   -2.255   10.208  123.402 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  28.6336     5.7023   5.021 1.02e-06 ***
Time          8.6091     0.4438  19.398  < 2e-16 ***
Diet4         2.1585     8.0805   0.267   0.7896    
Time:Diet4    1.1052     0.6331   1.746   0.0822 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 32.94 on 234 degrees of freedom
Multiple R-squared:  0.7837,    Adjusted R-squared:  0.7809 
F-statistic: 282.6 on 3 and 234 DF,  p-value: < 2.2e-16

Using planed contrasts I get the same result but a p-value of 0.150. Which method is correct?

> ChickWeight$Diet2 <- ChickWeight$Diet
> t2.vs.4 <- c(0,-1,0,1)
> dummy <- c(0,0,0,0)
> contrasts(ChickWeight$Diet2) <- cbind(t2.vs.4, dummy, dummy) 
> m3 <- lm(weight~Diet2*Time, data = ChickWeight)
> summary(m3)

Call:
lm(formula = weight ~ Diet2 * Time, data = ChickWeight)

Residuals:
     Min       1Q   Median       3Q      Max 
-126.034  -13.972   -0.055   13.451  160.534 

Coefficients: (4 not defined because of singularities)
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        27.4126     3.0160   9.089   <2e-16 ***
Diet2t2.vs.4        1.1829     4.7403   0.250    0.803    
Diet2dummy              NA         NA      NA       NA    
Diet2dummy              NA         NA      NA       NA    
Time                8.8121     0.2381  37.009   <2e-16 ***
Diet2t2.vs.4:Time   0.5357     0.3714   1.442    0.150    
Diet2dummy:Time         NA         NA      NA       NA    
Diet2dummy:Time         NA         NA      NA       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 38.65 on 574 degrees of freedom
Multiple R-squared:  0.7058,    Adjusted R-squared:  0.7043 
F-statistic:   459 on 3 and 574 DF,  p-value: < 2.2e-16

I worked out a better answer: By changing the order of the factors I no longer loose the information lost by sub-setting. And now have a p-value of 0.09197:

> ChickWeight$Diet <- factor(ChickWeight$Diet,levels=c(2,4,1,3))
> ChickWeight$Diet2 <- ChickWeight$Diet
> CW.lm <- lm(weight~Diet*Time,data=ChickWeight)
> summary(CW.lm)

Call:
lm(formula = weight ~ Diet * Time, data = ChickWeight)

Residuals:
     Min       1Q   Median       3Q      Max 
-135.425  -13.757   -1.311   11.069  130.391 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  28.6336     5.8972   4.855 1.55e-06 ***
Diet4         2.1585     8.3567   0.258  0.79627    
Diet1         2.2974     7.2672   0.316  0.75202    
Diet3       -10.3833     8.3399  -1.245  0.21364    
Time          8.6091     0.4590  18.757  < 2e-16 ***
Diet4:Time    1.1052     0.6548   1.688  0.09197 .  
Diet1:Time   -1.7673     0.5717  -3.092  0.00209 ** 
Diet3:Time    2.8137     0.6491   4.335 1.72e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 34.07 on 570 degrees of freedom
Multiple R-squared:  0.773,    Adjusted R-squared:  0.7702 
F-statistic: 277.3 on 7 and 570 DF,  p-value: < 2.2e-16

And following an example in Aho, Ken A.. Foundational and Applied Statistics for Biologists Using R (Page 433) my planned contrasts give the same answer:

> contrasts(ChickWeight$Diet2) <- contr.helmert(4)
> contrasts(ChickWeight$Diet2)
  [,1] [,2] [,3]
2   -1   -1   -1
4    1   -1   -1
1    0    2   -1
3    0    0    3
> CW.pc.lm <- lm(weight~Diet2*Time, data = ChickWeight)
> summary(CW.pc.lm)

Call:
lm(formula = weight ~ Diet2 * Time, data = ChickWeight)

Residuals:
     Min       1Q   Median       3Q      Max 
-135.425  -13.757   -1.311   11.069  130.391 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  27.1518     2.7686   9.807  < 2e-16 ***
Diet21        1.0793     4.1783   0.258   0.7963    
Diet22        0.4060     1.9859   0.204   0.8381    
Diet23       -2.9671     1.6685  -1.778   0.0759 .  
Time          9.1470     0.2173  42.093  < 2e-16 ***
Diet21:Time   0.5526     0.3274   1.688   0.0920 .  
Diet22:Time  -0.7733     0.1575  -4.909 1.20e-06 ***
Diet23:Time   0.7586     0.1302   5.827 9.46e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 34.07 on 570 degrees of freedom
Multiple R-squared:  0.773,    Adjusted R-squared:  0.7702 
F-statistic: 277.3 on 7 and 570 DF,  p-value: < 2.2e-16

However, I'd love to know what was wrong with my first attempt at planned contrasts. Thanks

Edit: After reading 'The R book' by Crawley, I believe the planned contrast (aka posteriori contrasts) need to be "orthogonal" (by Crawley's definition) and non-zero. Putting the Diet factors back in order, I got the same result as above a third time.

> t2.vs.4 <- c(0,-1,0,1)
> t2.4.vs.3 <- c(0,-0.5,1,-0.5)
> ChickWeight$Diet2 <- ChickWeight$Diet
> dummy <- c(-3,1,1,1)
> contrasts(ChickWeight$Diet2) <- cbind(t2.vs.4,t2.4.vs.3, dummy) 
> contrasts(ChickWeight$Diet2)
  t2.vs.4 t2.4.vs.3 dummy
1       0       0.0    -3
2      -1      -0.5     1
3       0       1.0     1
4       1      -0.5     1
> CW.pc.lm <- lm(weight~Diet2*Time, data = ChickWeight)
> summary(CW.pc.lm)

Call:
lm(formula = weight ~ Diet2 * Time, data = ChickWeight)

Residuals:
     Min       1Q   Median       3Q      Max 
-135.425  -13.757   -1.311   11.069  130.391 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          27.1518     2.7686   9.807  < 2e-16 ***
Diet2t2.vs.4          1.0793     4.1783   0.258    0.796    
Diet2t2.4.vs.3       -7.6417     4.8183  -1.586    0.113    
Diet2dummy           -1.2597     1.3615  -0.925    0.355    
Time                  9.1470     0.2173  42.093  < 2e-16 ***
Diet2t2.vs.4:Time     0.5526     0.3274   1.688    0.092 .  
Diet2t2.4.vs.3:Time   1.5074     0.3759   4.011 6.86e-05 ***
Diet2dummy:Time       0.7684     0.1082   7.104 3.62e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 34.07 on 570 degrees of freedom
Multiple R-squared:  0.773, Adjusted R-squared:  0.7702 
F-statistic: 277.3 on 7 and 570 DF,  p-value: < 2.2e-16
$\endgroup$
4
  • $\begingroup$ Your models include the interaction between time and diet. It means you assume that weight gain due to diet changes according time. Today diet 2 is better than diet 4; tomorrow diet 2 is worse than diet 4. So under this situation, How to compare diet 2 vs diet 4? $\endgroup$
    – user158565
    Oct 8, 2018 at 2:18
  • $\begingroup$ @a_statistician, I have edited the question for clarity, a need the interaction term so gradient can differ between diets. $\endgroup$
    – Greg F.
    Oct 8, 2018 at 2:38
  • $\begingroup$ The first part is OK, but some information was lost because you used subset of data. I cannot figure out how to get that point estimate 0.5337 in the second part, which is totally different from 1.1052 in the first part. I verified that 1.1052 is correct. $\endgroup$
    – user158565
    Oct 8, 2018 at 3:32
  • $\begingroup$ Thanks @a_statistician, I think I worked it out... $\endgroup$
    – Greg F.
    Oct 8, 2018 at 8:23

1 Answer 1

1
$\begingroup$

I have discussed this at length and came to the follow conclusions:

The dummy contrast (0,0,0,0) is incorrect and gives biased estimates of the parameters in the model.

To illustrate this, I’ve used four different methods to estimate the intercepts (at Time = 0) for each of the four diets. The methods are: · 'subset': fit simple regressions separately for each Diet · 'default': fit to whole dataset without messing with contrasts · 'crawley': fit to whole dataset using the contrasts plus c(-3,1,1,1) · 'dummy': fit to whole dataset using the contrasts plus c(0,0,0,0) Here are the estimates of the intercepts:

Diet  subset default crawley dummy
1      30.9    30.9   30.9    27.5
2      28.6    28.6   28.6    30.1
3      18.3    18.3   18.3    19.7
4      30.8    30.8   30.8    32.7

As you can see, the dummy method gives incorrect estimates. This occurs because the two planned contrasts alone only allow Diet 2 to differ from Diet 4, and Diet 3 to differ from the mean of Diets 2 and 4. Crawley’s final contrast allows Diet 1 to differ from the mean of Diets2-4. Put another way, omitting this contrast (using instead c(0,0,0,0)) essentially forces Diet1 to be equal to the mean of Diets2-4. This can be seen in the estimates in the table above.

> mean(c(30.1, 19.7, 32.7))
[1] 27.5

Interestingly, if you supply just the two planned contrasts, R will automatically provide the last one, which is equivalent to Crawley’s c(-3,1,1,1).

> ChickWeight$Diet2 <- ChickWeight$Diet
> t2.vs.4 <- c(0,-1,0,1) 
> t2.4.vs.3 <- c(0,-0.5,1,-0.5) 
> contrasts(ChickWeight$Diet2) <- cbind(t2.vs.4, t2.4.vs.3) 
> contrasts(ChickWeight$Diet2)
  t2.vs.4 t2.4.vs.3          
1       0       0.0 -0.8660254
2      -1      -0.5  0.2886751
3       0       1.0  0.2886751
4       1      -0.5  0.2886751

The issues with sub-setting are (1) it can open the door to multiple testing errors, although this doesn't appear to be a problem in this example; and (2) it excludes a whole lot of data that can be used to get a better estimate of the error variance and thus more accurate tests of hypotheses (loss of information).

Crawley's method (or supplying just two planned contrasts) appears to be the best method to answer this question.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.