Sum to Zero vs Dummy Coding

class: center, middle, inverse, title-slide

# <b>Sum to Zero vs Dummy Coding </b>
## Data Analysis for Psychology in R 2<br><br>
### dapR2 Team
### Department of Psychology<br>The University of Edinburgh
### AY 2020-2021

---

# Weeks Learning Objectives
1. Interpret the output from a model using dummy coding and sum-to-zero coding

2. Create specific contrast matrices to test specific effects.

3. Recognise other forms of contrasts.

4. Construct models to test factorial designs.

---
# Topics for today
+ Last time we looked at the `$F$`-test in one-way designs and linear models

+ This time we are going to consider contrasts and `$\beta$` coefficients

---
# Looking beneath the F-test
+ The `$F$`-test gives us an overall test of the model, or the effect of an experimental condition.

+ But we may want to know something more specific.
  + Differences between specific groups or sets of groups.

+ In such cases we talk about...
  + contrasts & planned comparisons
  + post-hoc test

+ So how do we approach these from the linear model perspective?

---
# Contrasts and Planned comparisons
+ Sometimes we want to make comparisons between pairs of things.
  + Treatment A vs Treatment B
  + Treatment A vs (Treatment B & Treatment C) etc.

+ Such comparisons can be... 
  + Specified a priori (confirmatory)
  + For all possible comparisons (exploratory)
  
+ We achieve these comparisons via assigning weights to groups.

+ May sound complicated, but we have already seen this practice in action this year
  
---
# Dummy coding
+ (Also called reference group coding.)

+ Create `$k$`-1 dummy variables/contrasts
  + where `$k$` is the number of levels of the categorical predictor.

+ Assign reference group 0 on all dummies.

+ Assign 1 to the focal group.

+ Enter the dummies into the linear model and they code the difference in means between the focal group/level and the reference.

+ We are going to use `$g$` from here on to be explicit these are experimental groups.
  + But if we think in terms of the categorical variables produced by a design `$g = k$`

---
# Why do we need a reference group?
+ Consider our example.

+ We have three groups each given a specific Treatment A, B or C

+ We want a model that represents our data (observations), but all we "know" are what groups an observation belongs to. So;

`$$y_{ij} = \mu_i + \epsilon_{ij}$$`

+ Where 
  + `$y_{ij}$` are the individual observations
  + `$\mu_i$` is the mean of group `$i$` and
  + `$\epsilon_{ij}$` is the individual deviation from that mean.

???
+ And this hopefully makes sense.
  + Given we know someone's group, our best guess is the mean
  + But people wont all score the mean, so there is some deviation for every person.

---
# Why do we need a reference group?
+ An alternative way to present this idea looks much more like our linear model:

`$$y_{ij} = \beta_0 + \underbrace{(\mu_{i} - \beta_0)}_{\beta_i} + \epsilon_{ij}$$`
+ Where 
  + `$y_{ij}$` are the individual observations
  + `$\beta_0$` is an estimate of reference/overall average
  + `$\mu_i$` is the mean of group `$i$` 
  + `$\beta_1$` is the difference between the reference and the mean of group `$i$`, and
  + `$\epsilon_{ij}$` is the individual deviation from that mean.

---
# Why do we need a reference group?
+ We can write this idea equation more generally:

$$\mu_i = \beta_0 + \beta_i $$

+ or for the specific groups (in our case 3):

`$$\mu_{treatmentA} = \beta_0 + \beta_{1A}$$`

`$$\mu_{treatmentB} = \beta_0 + \beta_{2B}$$`

`$$\mu_{treatmentC} = \beta_0 + \beta_{3C}$$`

+ **The problem**: we have four parameters ( `$\beta_0$` , `$\beta_{1A}$` , `$\beta_{2B}$` , `$\beta_{3C}$` ) to model three group means ( `$\mu_{TreatmentA}$` , `$\mu_{TreatmentB}$` , `$\mu_{TreatmentC}$` )
  + This means our model is under-identified.
  + We are trying to estimate too much with too little.

---
# Constraints fix identification
+ Consider dummy coding.

+ Suppose we make Treatment A the reference. Then,

`$$\mu_{treatmentA} = \beta_0$$`

`$$\mu_{treatmentB} = \beta_0 + \beta_{2B}$$`

`$$\mu_{treatmentC} = \beta_0 + \beta_{3C}$$`
+ Fixed!

+ We now only have three parameters ( `$\beta_0$` , `$\beta_{2B}$` , `$\beta_{3C}$` ) for the three group means ( `$\mu_{TreatmentA}$` , `$\mu_{TreatmentB}$` , `$\mu_{TreatmentC}$` ).

---
# Why not always use dummy coding?

+ We might not always want to compare against a reference group.

+ We might want to compare to:
  + The overall or grand mean
  + Group 1 vs groups 2, 3, 4 combined
  + and on we go!

+ Let's consider the example of the grand mean...

---
# Sum to zero constraint

+ With dummy coding we had a reference group, and the mean of that group was equal to the value of `$\beta_0$`, or

`$$\mu_{reference} = \beta_0$$`

+ Alternately, we can apply what is referred to as the sum to zero constraint (again using example of three levels).

`$$\beta_1 + \beta_2 + \beta_3 = 0$$`

+ There are two consequences of this constraint (see practical exercises for full explanation):

`$$\beta_0 = \frac{\mu_1 + \mu_2 + \mu_3}{3}$$`

+ And

`$$\mu_1 = \beta_0 + \beta_1$$`

`$$\mu_2 = \beta_0 + \beta_2$$`

`$$\mu_3 = \beta_0 - (\beta_1 + \beta_2)$$`

---
class: center, middle
# Time for a break

---
class: center, middle
# Welcome Back!

---
# OK, but how do we apply the constraint?

+ Answer, in the same way as we did with dummy coding.

+ We can create a set of sum to zero (sometimes called effect, or deviation) variables
  + Or the equivalent contrast matrix.
  
+ For effect code variables we:
  + Create `$g-1$` variables
  + For observations in the focal group, assign 1
  + For observations in the last group, assign -1
  + For all other groups assign 0

---
# Comparing coding matrices

.pull-left[
<table class="table" style="width: auto !important; margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> Level </th>
   <th style="text-align:right;"> D1 </th>
   <th style="text-align:right;"> D2 </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> Treatment A </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Treatment B </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 0 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Treatment C </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 1 </td>
  </tr>
</tbody>
</table>

`$$y_{ij} = \beta_0 + \beta_1D_1 + \beta_2D_2 + \epsilon_{ij}$$`

]

.pull-right[

<table class="table" style="width: auto !important; margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> Level </th>
   <th style="text-align:right;"> E1 </th>
   <th style="text-align:right;"> E2 </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> Treatment A </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 0 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Treatment B </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 1 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Treatment C </td>
   <td style="text-align:right;"> -1 </td>
   <td style="text-align:right;"> -1 </td>
  </tr>
</tbody>
</table>

`$$y_{ij} = \beta_0 + \beta_1E_1 + \beta_2E_2 + \epsilon_{ij}$$`

]

---
# Sum to zero/effects for group means

`$$\mu_1 = \beta_0 + 1*\beta_1 + 0*\beta_2 = \beta_0 + \beta_1$$`

`$$\mu_2 = \beta_0 + 0*\beta_1 + 1*\beta_2 = \beta_0 + \beta_2$$`

`$$\mu_3 = \beta_0 -1*\beta_1 -1*\beta_2 = \beta_0 - \beta_1 -\beta_2$$`

+ Now we will look practically at the implementation and differences

---
# The data

```r
hosp_tbl <- read_csv("hospital.csv", col_types = "dff")
hosp_tbl %>%
  select(SWB, Treatment) %>%
  slice(1:10)
```

```
## # A tibble: 10 x 2
##      SWB Treatment
##    <dbl> <fct>    
##  1   6.2 TreatA   
##  2  15.9 TreatA   
##  3   7.2 TreatA   
##  4  11.3 TreatA   
##  5  11.2 TreatA   
##  6   9   TreatA   
##  7  14.5 TreatA   
##  8   7.3 TreatA   
##  9  13.7 TreatA   
## 10  12.6 TreatA
```

---
# Group Means

```r
hosp_tbl %>%
  select(1:2) %>%
  group_by(Treatment) %>%
  summarise(
    mean = round(mean(SWB),3),
    sd = round(sd(SWB),1),
    N = n()
  )
```

```
## `summarise()` ungrouping output (override with `.groups` argument)
```

```
## # A tibble: 3 x 4
##   Treatment  mean    sd     N
##   <fct>     <dbl> <dbl> <int>
## 1 TreatA     9.33   2.9    60
## 2 TreatB    11.3    2.5    60
## 3 TreatC     9.04   2      60
```

---
# Dummy (reference) model

```r
summary(lm(SWB ~ Treatment, data = hosp_tbl))
```

```
## 
## Call:
## lm(formula = SWB ~ Treatment, data = hosp_tbl)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.373 -1.987 -0.300  1.838  7.173 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       9.3267     0.3242  28.770  < 2e-16 ***
## TreatmentTreatB   1.9467     0.4585   4.246 3.51e-05 ***
## TreatmentTreatC  -0.2850     0.4585  -0.622    0.535    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.511 on 177 degrees of freedom
## Multiple R-squared:  0.1369,	Adjusted R-squared:  0.1271 
## F-statistic: 14.04 on 2 and 177 DF,  p-value: 2.196e-06
```

---
# Dummy (reference) model

.pull-left[

```
##     (Intercept) TreatmentTreatB TreatmentTreatC 
##           9.327           1.947          -0.285
```

+ Recall the equations for the group means:

`$$\mu_{treatmentA} = \beta_0$$`
`$$\mu_{treatmentB} = \beta_0 + \beta_1$$`

`$$\mu_{treatmentC} = \beta_0 + \beta_2$$`

]

.pull-right[

```
## `summarise()` ungrouping output (override with `.groups` argument)
```

<table class="table" style="width: auto !important; margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> Treatment </th>
   <th style="text-align:right;"> mean </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> TreatA </td>
   <td style="text-align:right;"> 9.327 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> TreatB </td>
   <td style="text-align:right;"> 11.273 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> TreatC </td>
   <td style="text-align:right;"> 9.042 </td>
  </tr>
</tbody>
</table>

]

---
# Effects (sum to zero) model

+ We need to change the contrast scheme from default.

```r
contrasts(hosp_tbl$Treatment) <- contr.sum 
contrasts(hosp_tbl$Treatment)
```

```
##        [,1] [,2]
## TreatA    1    0
## TreatB    0    1
## TreatC   -1   -1
```

---
# Effects (sum to zero) model

```r
summary(lm(SWB ~ Treatment, data = hosp_tbl))
```

```
## 
## Call:
## lm(formula = SWB ~ Treatment, data = hosp_tbl)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.373 -1.987 -0.300  1.838  7.173 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.8806     0.1872  52.791  < 2e-16 ***
## Treatment1   -0.5539     0.2647  -2.093   0.0378 *  
## Treatment2    1.3928     0.2647   5.262 4.09e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.511 on 177 degrees of freedom
## Multiple R-squared:  0.1369,	Adjusted R-squared:  0.1271 
## F-statistic: 14.04 on 2 and 177 DF,  p-value: 2.196e-06
```

---
# Effects (sum to zero) model

.pull-left[

```
## (Intercept)  Treatment1  Treatment2 
##       9.881      -0.554       1.393
```

+ Recall the equations for the group means:

`$$\beta_0 = \frac{\mu_1 + \mu_2 + \mu_3}{3}$$`

`$$\mu_1 = \beta_0 + \beta_1$$`

`$$\mu_2 = \beta_0 + \beta_2$$`

`$$\mu_3 = \beta_0 - (\beta_1 + \beta_2)$$`
]

.pull-right[

```
## `summarise()` ungrouping output (override with `.groups` argument)
```

<table class="table" style="width: auto !important; margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;"> Treatment </th>
   <th style="text-align:right;"> mean </th>
   <th style="text-align:right;"> Gmean </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> TreatA </td>
   <td style="text-align:right;"> 9.327 </td>
   <td style="text-align:right;"> 9.881 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> TreatB </td>
   <td style="text-align:right;"> 11.273 </td>
   <td style="text-align:right;"> 9.881 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> TreatC </td>
   <td style="text-align:right;"> 9.042 </td>
   <td style="text-align:right;"> 9.881 </td>
  </tr>
</tbody>
</table>

]

---
# Manual contrast testing
+ We can structure a wide variety of contrasts so long as they can be written

1. A as a linear combination of population means
2. The associated coefficients sum to zero

+ So

$$H_0: c_1\mu_1 + c_1\mu_2 + c_3\mu_3 $$

+ With

`$$c_1 + c_2 + c_3 = 0$$`

---
# Manual contrast testing
+ For example, say we wanted to compare TreatmentA to TreatmentB and TreatmentC combined.

+ Our null becomes:

`$$H_0: \mu_{treatmentA} = \frac{1}{2}(\mu_{treatmentB} + \mu_{treatmentC})$$`

+ Or, is the average SWB of TreatmentA group equal to the average SWB of the combined TreatmentB and TreatmentC

+ We can set the contrasts here:
  + `$c_1 = 1$`
  + `$c_2 = -1/2$`
  + `$c_3 = -1/2$`
  
  
---
# The wide world of contrasts 
+ We have now seen two examples of coding schemes (dummy and effect).

+ We have also seen that so long as we apply some set of constraints, we are able to do different things.

+ This means there are **lots** of different coding's we can use for categorical variables to make different comparisons.
  + If you are interested, see the excellent resource on [UCLA website](https://stats.idre.ucla.edu/r/library/r-library-contrast-coding-systems-for-categorical-variables/)

+ These can include some custom contrasts (planned comparisons).
  + Next lecture.

---
# Summary of today

+ We have considered different ways in which we can code categorical predictors.

+ Take home:
  + Use of coding matrices allows us to compare groups (or levels) in lots of ways.
  + Our `$\beta$`'s will represent differences in group means.
  + The scheme we use determines which groups.
  + This makes it very flexible for testing hypotheses.
  + We will come to see this is very useful for testing specific hypotheses.
  + Effects (sum to zero, or deviation) coding = traditional ANOVA