2
$\begingroup$

I am trying to understand how contrasts work, so I ran a small simulation using the following R code:

hsb2 <- read_csv("http://www.ats.ucla.edu/stat/data/hsb2.csv")
hsb2$race.f <- factor(hsb2$race, labels=c("Hispanic", "Asian", "African-Am", "Caucasian"))

# method 1
(contrasts(hsb2$race.f) <- matrix(c(0,1,0,0,0,0,1,0,0,0,0,1),ncol = 3, byrow=FALSE))
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    1    0    0
[3,]    0    1    0
[4,]    0    0    1
m1 <- lm(write ~ race.f, hsb2)


# method 2
(contrasts(hsb2$race.f) <- matrix(c(-1,1,0,0,-1,0,1,0,-1,0,0,1),ncol = 3, byrow=FALSE))
     [,1] [,2] [,3]
[1,]   -1   -1   -1
[2,]    1    0    0
[3,]    0    1    0
[4,]    0    0    1
m2 <- lm(write ~ 0 + race.f, hsb2)

I was expecting that method 1 and method 2 return the same set of coefficients because:

  • In method 1, the three groups are compared to the reference group, and the contrast reflects that
  • In method 2, I planned to take out the intercept in the linear model, so that all factors are in the model. To emulate the scenario in the method1, I encoded that information in the contrast.

But I didn't get the same set of coefficients:

summary(m1) returns:

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   46.458      1.842  25.218  < 2e-16 ***
race.f1       11.542      3.286   3.512 0.000552 ***
race.f2        1.742      2.732   0.637 0.524613    
race.f3        7.597      1.989   3.820 0.000179 ***

summary(m2) returns:

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
race.fHispanic    46.4583     1.8422   25.22   <2e-16 ***
race.fAsian       58.0000     2.7212   21.31   <2e-16 ***
race.fAfrican-Am  48.2000     2.0181   23.88   <2e-16 ***
race.fCaucasian   54.0552     0.7495   72.12   <2e-16 ***

Clearly, I am not understanding contrasts correctly. Where do my reasoning fail and what methods would have made the contrast 1 and contrast 2 produce the same results(or is it just a complete nonsense?)

$\endgroup$
1
$\begingroup$

Using the 0+ syntax in the model specification will bypass the defined contrasts, and give you separately estimated coefficients for each level. Consider the examples below - you get the same estimates when you remove the intercept (0+) no matter what contrasts you use:

#Dummy coding
> contr.treatment(4)
  2 3 4
1 0 0 0
2 1 0 0
3 0 1 0
4 0 0 1

> contrasts(hsb2$race.f) = contr.treatment(4)
#model with intercept
> summary(lm(write ~ race.f, data=hsb2))

Call:
lm(formula = write ~ race.f, data = hsb2)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.0552  -5.4583   0.9724   7.0000  18.8000 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   46.458      1.842  25.218  < 2e-16 ***
race.f2       11.542      3.286   3.512 0.000552 ***
race.f3        1.742      2.732   0.637 0.524613    
race.f4        7.597      1.989   3.820 0.000179 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.025 on 196 degrees of freedom
Multiple R-squared:  0.1071,    Adjusted R-squared:  0.0934 
F-statistic: 7.833 on 3 and 196 DF,  p-value: 5.785e-05

# removing intercept
> summary(lm(write ~ 0 + race.f, data=hsb2))

Call:
lm(formula = write ~ 0 + race.f, data = hsb2)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.0552  -5.4583   0.9724   7.0000  18.8000 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
race.fHispanic    46.4583     1.8422   25.22   <2e-16 ***
race.fAsian       58.0000     2.7212   21.31   <2e-16 ***
race.fAfrican-Am  48.2000     2.0181   23.88   <2e-16 ***
race.fCaucasian   54.0552     0.7495   72.12   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.025 on 196 degrees of freedom
Multiple R-squared:  0.9722,    Adjusted R-squared:  0.9717 
F-statistic:  1716 on 4 and 196 DF,  p-value: < 2.2e-16

# Deviation coding
> contr.sum(4)
  [,1] [,2] [,3]
1    1    0    0
2    0    1    0
3    0    0    1
4   -1   -1   -1

> contrasts(hsb2$race.f) = contr.sum(4)
# model with intercept
> summary(lm(write ~ race.f, data=hsb2))

Call:
lm(formula = write ~ race.f, data = hsb2)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.0552  -5.4583   0.9724   7.0000  18.8000 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  51.6784     0.9821  52.619  < 2e-16 ***
race.f1      -5.2200     1.6314  -3.200  0.00160 ** 
race.f2       6.3216     2.1603   2.926  0.00384 ** 
race.f3      -3.4784     1.7323  -2.008  0.04602 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.025 on 196 degrees of freedom
Multiple R-squared:  0.1071,    Adjusted R-squared:  0.0934 
F-statistic: 7.833 on 3 and 196 DF,  p-value: 5.785e-05

#model without intercept
> summary(lm(write ~ 0 + race.f, data=hsb2))

Call:
lm(formula = write ~ 0 + race.f, data = hsb2)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.0552  -5.4583   0.9724   7.0000  18.8000 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
race.fHispanic    46.4583     1.8422   25.22   <2e-16 ***
race.fAsian       58.0000     2.7212   21.31   <2e-16 ***
race.fAfrican-Am  48.2000     2.0181   23.88   <2e-16 ***
race.fCaucasian   54.0552     0.7495   72.12   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.025 on 196 degrees of freedom
Multiple R-squared:  0.9722,    Adjusted R-squared:  0.9717 
F-statistic:  1716 on 4 and 196 DF,  p-value: < 2.2e-16
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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