Class Outline

GPT AI Tutor for ESRM64103
Planned Contrast
Example: Group Comparison - STEM vs. Non-STEM Groups
- ANOVA-style t statistics and regression-style coding in R
- Effect sizes

2 Planned Contrasts

2.1 Definition: Planned Contrasts

Pre-defined:
- Unlike post hoc tests, planned contrasts are determined before looking at the data, meaning the researcher has a specific hypothesis about which groups to compare.
- Definition: Planned contrasts are hypothesis-driven comparisons made before data collection.
Weights assigned:
- To perform a planned contrast, each group is assigned a numerical “weight” which reflects its role in the comparison, with the weights usually summing to zero.
$D = weights * means$
- D: unscaled group differences given coding scheme

2.2 Example of Planned Contrasts

Imagine a study comparing the effects of three different study methods (A, B, C) on test scores.
- One planned contrast might be to compare the average score of method A (considered the “experimental” method) against the combined average of methods B and C (considered the “control” conditions),
- Testing the hypothesis that method A leads to significantly higher scores than the traditional methods.
- $H_0: \mu_{A} = \frac{\mu_B+\mu_C}{2}$ ; we also call this a complex contrast
When to use planned contrasts:
- When you have a clear theoretical basis for predicting specific differences between groups in your study.
- When you are only interested in a few specific comparisons, not all possible pairwise comparisons.

2.3 What Does Each Contrast Tell Us?

Each contrast is a mean comparison (via t-test).
Simple contrast (pairwise) compares two individual group means.
Complex contrast compares a combination of group means.
Must be theoretically justified for meaningful interpretation.

2.4 Simple vs. Complex Comparisons

Simple Comparison: Two groups directly compared.
- Example: $H_0: \mu_2 = \mu_3$
Complex Comparison: Combines means of multiple groups.
- Example: $H_0: \frac{(\mu_1 + \mu_2)}{2} = \mu_3$
- Example: $H_0: \frac{(\mu_1 + \mu_2 + \mu_3)}{3} = \frac{(\mu_4 + \mu_5)}{2}$

Note

We should not test all possible combinations of groups. Instead, justify your comparison plan before performing statistic analysis.

2.5 Today’s focus: Complex Comparisons

We performed omnibus tests in the last lecture, which provide all pairwise group comparisons (simple contrasts)
Today we focus more on complex contrasts.
- Helmert contrast: Compares each mean to the mean of subsequent groups.
- Sum (deviation) contrast: each group compared to the grand mean
- Polynomial contrast: Tests for trends in ordered data.
By default, R uses treatment contrasts: each group compared to the reference group
- Question: is “treatment contrast” a simple or complex contrast?
- G1 vs. Treatment (reference)
- G2 vs. Treatment (reference)
- ….

2.6 Orthogonal vs. Non-Orthogonal Contrasts

Orthogonal Contrasts: Independent from each other, sum of product of weights equals zero.
Non-Orthogonal Contrasts: Not independent, lead to inflated Type I error rates.

Note

Orthogonal contrasts allow clear interpretation without redundancy.
Orthogonal contrasts follow a series of group comparisons that do not overlap variances.

2.7 Orthogonal contrasts from variances: no redundancy

Helmert contrast example

With a logical control group, a good first contrast compares all treatment groups to the one control group.
To get each level of the IV alone, you should have one fewer contrast than your number of IV levels (3 levels = 2 contrasts)
Once an IV level appears by itself, it shouldn’t reappear in subsequent contrasts

2.8 Example of Orthogonal Contrasts

Contrast 1: g3 vs. (g1, g2)

Contrast 2: g1 vs. g2

2.9 Orthogonal Planned Contrasts

If the same exact combination of means is not found in more than one contrast, the contrasts are independent (orthogonal)
- Check this by ensuring that the product of the weights across all contrasts sums to zero
For an orthogonal comparison, contrasts are independent of each other:
- We weight the means included on each side of the contrast
- Each contrast has weights that sum to 0
- Groups not in the contrast get a weight of 0
Why does independence matter?
- Type I error rate is unaffected by independent (orthogonal) contrasts
- Interpretation of contrasts is cleaner because contrasts aren’t related (you’ve isolated effects)

Group	Contrast 1	Contrast 2	Product
G1	+1	-1	-1
G2	+1	+1	+1
G3	-2	0	0
Sum	0	0	0

2.10 Contrasts’ Independence checking in R

⌘+C

cat("## Constrasts and Coding \n")
contras <- matrix(
  c(1, 1, -2,
    -1, 1, 0), ncol = 2
)
contras

## Constrasts and Coding 
     [,1] [,2]
[1,]    1   -1
[2,]    1    1
[3,]   -2    0

To understand if these contrasts are orthogonal, one could compute the dot product between the contrast vectors. The dot product of first and second contrast is:

$(1 * -1) + (1 * 1) + (-2 * 0) = -1 + 1 + 0 = 0$

Since the dot product of the contrasts is equal to zero, these contrasts are indeed orthogonal.

Orthogonal contrasts have the advantage that the comparison of the means of one contrast does not affect the comparison of the means of any other contrast.

⌘+C

cat("## the dot-product of two contrasts should be zero\n")
t(contras[,1]) %*% contras[,2] 
cat("## if non-diagonal elements are zero and diagnoal elements = 0 \n")
crossprod(contras)

## the dot-product of two contrasts should be zero
     [,1]
[1,]    0
## if non-diagonal elements are zero and diagnoal elements = 0 
     [,1] [,2]
[1,]    6    0
[2,]    0    2

new_code <- matrix(c(
0, 0, 0, 0,
1, 0, 0, 0,
0, 1, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1
), byrow = TRUE, nrow = 5)

cat("## the cross-product of two contrasts should be zero\n")
t(new_code[,1]) %*% new_code[,2]
t(new_code[,1]) %*% new_code[,3]
t(new_code[,1]) %*% new_code[,4]
t(new_code[,2]) %*% new_code[,3]
cat("## if non-diagonal elements are zero and diagnoal elements = 0 \n")
crossprod(new_code)

2.11 Computing Planned Contrasts

Formula for contrast value: $C = c_1\mu_1 + c_2\mu_2 + \dots + c_k\mu_k$
Test statistic: $t = \frac{C}{\sqrt{MSE \sum \frac{c_i^2}{n_i}}}$
- $MSE$ : Mean Square Error from ANOVA
- $c_i$ : Contrast coefficients
- $n_i$ : Sample size per group

2.12 Planned Contrasts and Coding Scheme

The relationship between planned contrasts in ANOVA and coding in regression lies in how categorical variables are represented and interpreted in statistical models.
Both approaches aim to test specific hypotheses about group differences, but their implementation varies based on the framework:
- ANOVA focuses on partitioning variance,
- Regression interprets categorical predictors through coding schemes.

2.12.1 Planned Contrast and Linear Regression

Planned contrast can be done using linear regression + contrasts
Let’s look at the default contrasts plan: treatment contrasts == dummy coding

library(tidyverse)
library(kableExtra)
library(here)
set.seed(42)
dt <- read.csv(here("teaching/2025-01-13-Experiment-Design/Lecture05","week5_example.csv"))
## Treatment contrast matrix 
attributes(C(dt$group, treatment, 4))$contrasts

   g2 g3 g4 g5
g1  0  0  0  0
g2  1  0  0  0
g3  0  1  0  0
g4  0  0  1  0
g5  0  0  0  1

## Sum contrast matrix 
attributes(C(dt$group, sum, 4))$contrasts

   [,1] [,2] [,3] [,4]
g1    1    0    0    0
g2    0    1    0    0
g3    0    0    1    0
g4    0    0    0    1
g5   -1   -1   -1   -1

attributes(C(dt$group, helmert, 4))$contrasts

   [,1] [,2] [,3] [,4]
g1   -1   -1   -1   -1
g2    1   -1   -1   -1
g3    0    2   -1   -1
g4    0    0    3   -1
g5    0    0    0    4

crossprod(attributes(C(dt$group, treatment, 4))$contrasts)

   g2 g3 g4 g5
g2  1  0  0  0
g3  0  1  0  0
g4  0  0  1  0
g5  0  0  0  1

3 Example - STEM vs. Non-STEM Groups

3.1 Background

Hypothesis: STEM students have different growth mindset scores than non-STEM students.
Weights assigned:
- STEM (Engineering, Chemistry): $+\frac{1}{2}$
- Non-STEM (Education, Political Science, Psychology): $-\frac{1}{3}$
Compute contrast value and test using t-statistic.

3.2 Set Contrasts in R

⌘+C

# Set seed for reproducibility
options(digits = 5)
summary_tbl <- dt |> 
  group_by(group) |> 
  summarise(
    N = n(),
    Mean = mean(score),
    SD = sd(score),
    shapiro.test.p.values = shapiro.test(score)$p.value
  ) |> 
  mutate(department = c("Engineering", "Education", "Chemistry", "Political", "Psychology")) |> 
  relocate(group, department)
summary_tbl

group	department	N	Mean	SD	shapiro.test.p.values
g1	Engineering	28	4.2500	3.15054	0.07759
g2	Education	28	2.7589	2.19478	0.07605
g3	Chemistry	28	3.5446	2.86506	0.00623
g4	Political	28	3.8568	0.58325	0.03023
g5	Psychology	28	2.0243	1.30911	0.06147

Homogeneity of variance assumption: Levene’s test

3.2.1 Stem vs. Non-Stem

⌘+C

aov_fit <- aov(score ~ group, data = dt)
car::leveneTest(aov_fit)

Levene's Test for Homogeneity of Variance (center = median)
       Df F value  Pr(>F)    
group   4      13 5.8e-09 ***
      135                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Even though the assumption checks did not pass using the original categorical levels, we may still be interested in different group contrasts.

⌘+C

coef(lm(score ~ group, data = dt))

(Intercept)     groupg2     groupg3     groupg4     groupg5 
    4.25000    -1.49107    -0.70536    -0.39321    -2.22571

$(\beta_0 + \beta_0 + \beta_2) / 2 - (\beta_0 + \beta_1 + \beta_0 + \beta_3 + \beta_0 + \beta_4) / 3 \\ = \beta_2 / 2 - (\beta_1 + \beta_3 + \beta_3) /3$

⌘+C

contrast <- matrix(
  c(0, -1/3, 1/2, -1/3, -1/3), 
  nrow = 1
)
rownames(contrast) <- "Stem vs. Non-Stem"
summary(multcomp::glht(aov_fit, linfct = contrast))


     Simultaneous Tests for General Linear Hypotheses

Fit: aov(formula = score ~ group, data = dt)

Linear Hypotheses:
                       Estimate Std. Error t value Pr(>|t|)   
Stem vs. Non-Stem == 0    1.017      0.386    2.64   0.0093 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

means <- summary_tbl$Mean
mean(means[c(1,3)]) - mean(means[-c(1,3)])

[1] 1.0173

3.2.2 Complex Contrast Matrix

There are multiple “canned” contrasts: Helmert, sum (effect coding), treatment

For example, Helmert: four contrasts

g1 vs. g2: $\mu_{Engineering} = \mu_{Education}$
$\frac{g1+g2}{2}$ vs. g3: $\mu_{non-Chemistry} = \mu_{Chemistry}$
$\frac{g1+g2+g3}{3}$ vs. g4: $\mu_{non-Political} = \mu_{Political}$
$\frac{g1+g2+g3+g4}{4}$ vs. g5: $\mu_{non-Psychology} = \mu_{Psychology}$

Summary Statistics:

⌘+C

dt$group <- factor(dt$group, levels = c("g1", "g2", "g3", "g4", "g5"))
groups <- levels(dt$group)
cH <- contr.helmert(groups) # pre-defined four contrasts
colnames(cH) <- paste0("Ctras", 1:4)
summary_ctras_tbl <- cbind(summary_tbl, cH)
summary_ctras_tbl

   group  department  N   Mean      SD shapiro.test.p.values Ctras1 Ctras2
g1    g1 Engineering 28 4.2500 3.15054             0.0775874     -1     -1
g2    g2   Education 28 2.7589 2.19478             0.0760542      1     -1
g3    g3   Chemistry 28 3.5446 2.86506             0.0062253      0      2
g4    g4   Political 28 3.8568 0.58325             0.0302312      0      0
g5    g5  Psychology 28 2.0243 1.30911             0.0614743      0      0
   Ctras3 Ctras4
g1     -1     -1
g2     -1     -1
g3     -1     -1
g4      3     -1
g5      0      4

3.2.3 Helmert contrasts are Orthogonal

apply(cH, 2, sum)

Ctras1 Ctras2 Ctras3 Ctras4 
     0      0      0      0

crossprod(cH) # diagonal -- columns are orthogonal

       Ctras1 Ctras2 Ctras3 Ctras4
Ctras1      2      0      0      0
Ctras2      0      6      0      0
Ctras3      0      0     12      0
Ctras4      0      0      0     20

summary(aov(score ~ group, dt))

             Df Sum Sq Mean Sq F value Pr(>F)   
group         4     89    22.3    4.47  0.002 **
Residuals   135    675     5.0                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.3 ANOVA: t-value formula for Defined Contrast Matrix

$t = \frac{C}{\sqrt{MSE \sum \frac{c_i^2}{n_i}}}$

Sum_C2_n <- colSums(cH^2 / summary_tbl$N)
C <- crossprod(summary_tbl$Mean, cH)
MSE <- 5
t <- as.numeric(C / sqrt(MSE * Sum_C2_n))
tibble(
  t_value = t,
  p_value = pt(t, df = 135) ## p-values
)

# A tibble: 4 × 2
  t_value  p_value
    <dbl>    <dbl>
1 -2.50   0.00690 
2  0.0776 0.531   
3  0.695  0.756   
4 -3.34   0.000541

g1 vs. g2: We reject the null and determine that the mean of Education is different from the mean of Engineering in their growth mindset scores (p = 0.531).
$\frac{g1+g2}{2}$ vs. g3: We retain the null and determine that the mean of Chemistry is not significantly different from the mean of Education and Engineering in their growth mindset scores (p = 0.531).

3.3.1 Helmert Contrast

Remember the planned contrast: g1 vs. g2 from the Helmert contrast:

t-value: -2.495
p-value: 0.0069
df: 134

⌘+C

contrasts(dt$group) <- "contr.helmert"
fit_helmert <- lm(score ~ group, dt)
contr.helmert(levels(dt$group))
summary(fit_helmert)$coefficients

   [,1] [,2] [,3] [,4]
g1   -1   -1   -1   -1
g2    1   -1   -1   -1
g3    0    2   -1   -1
g4    0    0    3   -1
g5    0    0    0    4
             Estimate Std. Error   t value   Pr(>|t|)
(Intercept)  3.286929   0.189003 17.390874 2.8188e-36
group1      -0.745536   0.298840 -2.494765 1.3810e-02
group2       0.013393   0.172535  0.077624 9.3824e-01
group3       0.084732   0.122001  0.694520 4.8855e-01
group4      -0.315661   0.094502 -3.340271 1.0825e-03

mean(summary_tbl$Mean)

[1] 3.2869

unique(cbind(model.matrix(fit_helmert), group = dt$group))

    (Intercept) group1 group2 group3 group4 group
1             1     -1     -1     -1     -1     1
29            1      1     -1     -1     -1     2
57            1      0      2     -1     -1     3
85            1      0      0      3     -1     4
113           1      0      0      0      4     5

summary_tbl$Mean %*% contrasts(dt$group) / c(2, 6, 12, 20)

         [,1]     [,2]     [,3]     [,4]
[1,] -0.74554 0.013393 0.084732 -0.31566

3.3.2 Treatment contrasts

For treatment contrasts, four dummy variables are created to compare:
- G1 (ref) vs. G2
- G1 (ref) vs. G3
- G1 (ref) vs. G4
- G1 (ref) vs. G5

Intercept: G1’s mean
group2: G2 vs. G1
group3: G3 vs. G1
group4: G4 vs. G1
group5: G5 vs. G1

library(multcomp)
contrasts(dt$group) <- "contr.treatment"
fit <- lm(score ~ group, dt)
unique(cbind(model.matrix(fit), group = dt$group))

    (Intercept) groupg2 groupg3 groupg4 groupg5 group
1             1       0       0       0       0     1
29            1       1       0       0       0     2
57            1       0       1       0       0     3
85            1       0       0       1       0     4
113           1       0       0       0       1     5

summary(fit)$coefficients

            Estimate Std. Error t value   Pr(>|t|)
(Intercept)  4.25000    0.42262 10.0562 4.2275e-18
groupg2     -1.49107    0.59768 -2.4948 1.3810e-02
groupg3     -0.70536    0.59768 -1.1802 2.4001e-01
groupg4     -0.39321    0.59768 -0.6579 5.1172e-01
groupg5     -2.22571    0.59768 -3.7239 2.8718e-04

3.3.3 Sum Contrasts

Another type of coding is effect coding. In R, the corresponding contrast type is the so-called sum contrasts.
A detailed post about sum contrasts can be found here
With sum contrasts, the reference level is the grand mean.
- $\frac{g1+g2+g3+g4+g5}{5}$ vs. g1/g2/g3/g4: the difference between mean score of g1 with grand mean across all five groups

contrasts(dt$group) <- "contr.sum"
fit2 <- lm(score ~ group, dt)
contr.sum(levels(dt$group))

   [,1] [,2] [,3] [,4]
g1    1    0    0    0
g2    0    1    0    0
g3    0    0    1    0
g4    0    0    0    1
g5   -1   -1   -1   -1

summary(fit2)$coefficients

            Estimate Std. Error  t value   Pr(>|t|)
(Intercept)  3.28693    0.18900 17.39087 2.8188e-36
group1       0.96307    0.37801  2.54777 1.1962e-02
group2      -0.52800    0.37801 -1.39680 1.6476e-01
group3       0.25771    0.37801  0.68177 4.9655e-01
group4       0.56986    0.37801  1.50753 1.3401e-01

mean(dt$score) # (Intercept) grand mean

[1] 3.2869

tibble(
  Label = paste0("group", 1:4),
  Estimate = summary_tbl$Mean[1:4] - mean(dt$score) 
)

# A tibble: 4 × 2
  Label  Estimate
  <chr>     <dbl>
1 group1    0.963
2 group2   -0.528
3 group3    0.258
4 group4    0.570

3.4 Effect Coding (Deviation Coding)

In modern statistics, regression-style coding is statistically equivalent to an ANOVA-style contrast matrix.
- Equivalent to ANOVA-style contrasts. (we will use this in R to reproduce ANOVA-style contrast matrix)
Compares each level to the grand mean.

Note

Effect coding is a method of encoding categorical variables in regression models, similar to dummy coding, but with a different interpretation of the resulting coefficients. It is particularly useful when researchers want to compare each level of a categorical variable to the overall mean rather than to a specific reference category.

3.4.1 Definition and Representation

In effect coding, categorical variables are transformed into numerical variables, typically using values of -1, 0, and 1. The key difference from dummy coding is that the reference category is represented by -1 instead of 0, and the coefficients indicate deviations from the grand mean.

For a categorical variable with k levels, effect coding requires k-1 coded variables. If we have a categorical variable X with three levels: $A, B, C$ , the effect coding scheme could be:

Category	$X_1$	$X_2$
A	1	0
B	0	1
C (reference)	-1	-1

The last category ( $C$ ) is the reference group, coded as -1 for all indicator

In effect coding, one level of the categorical variable (usually the last one or the one that is considered a ‘reference’ category) is coded as -1, and the others are coded as +1 or 0, suggesting their relation to the overall mean, rather than to a specific category.

Here’s how you could effect code a categorical variable with three levels (e.g., groups):

    [,1] [,2]
g1   1    0
g2   0    1
g3  -1   -1

This coding scheme shows that:

Group1 (g1) is compared to the overall effect by coding it as 1 in the first column and 0 in the second, suggesting it is above or below the overall mean.
Conversely, Group2 (g2) is also contrasted with the overall effect.
Group3 (g3), coded as -1 in both columns, serves as the reference category against which the other two are compared.

In your regression output, the coefficients for the first two groups will show how the mean of these groups differs from the overall mean. The intercept will represent the overall mean across all groups.

3.4.2 Interpretation of Coefficients in effective coding

When effect coding is used in a regression model:

$Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \epsilon$

$X_1$ and $X_2$ are coded variables. They have no much meaning, but their coefficients are important
$\beta_0$ represents the grand mean of $Y$ across all categories.
$\beta_1$ and $\beta_2$ represent the deviation of categories $A$ and $B$ from the grand mean.
The reference group ( $C$ ) does not have a separate coefficient; instead, its deviation can be inferred as $-(\beta_1 + \beta_2)$ .

⌘+C

library(ggplot2)
# Create a data frame for text labels
text_data <- data.frame(
  x = rep(0.25, 3),  # Repeating the same x-coordinate
  y = c(0.3, 0.7, 0.9),  # Different y-coordinates
  label = c("C: beta[0] - beta[1] - beta[2]", 
            "A: beta[0] + 1*'×'*beta[1] + 0*'×'*beta[2]", 
            "B: beta[0] + 0*'×'*beta[1] + 1*'×'*beta[2]")  # Labels
)

# Create an empty ggplot with defined limits
ggplot() +
  geom_text(data = text_data, aes(x = x, y = y, label = label), parse = TRUE, size = 11) +
  # Add a vertical line at x = 0.5
  # geom_vline(xintercept = 0.5, color = "blue", linetype = "dashed", linewidth = 1) +
  # Add two horizontal lines at y = 0.3 and y = 0.7
  geom_hline(yintercept = c(0.35, 0.75, 0.95), color = "red", linetype = "solid", linewidth = 1) +
  geom_hline(yintercept = 0.5, color = "grey", linetype = "solid", linewidth = 1) +
  geom_text(aes(x = .25, y = .45, label = "grand mean of Y"), color = "grey", size = 11) +
  # Set axis limits
  xlim(0, 1) + ylim(0, 1) +
  labs(y = "Y", x = "") +
  # Theme adjustments
  theme_minimal() +
  theme(text = element_text(size = 20))

3.4.3 Comparison to Dummy Coding

Dummy Coding: Compares each category to a specific reference category (e.g., comparing A and B to C).

Category	$X_1$	$X_2$
A	1	0
B	0	1
C (reference)	0	0

Effect Coding: Compares each category to the grand mean rather than a single reference category.

3.4.4 Use Cases

Effect coding is beneficial when:

There is no natural baseline category, and comparisons to the overall mean are more meaningful.
Researchers want to maintain sum-to-zero constraints for categorical variables in linear models.
In ANOVA-style analyses, where main effects and interaction effects are tested under an equal-weight assumption.

3.4.5 Implementation in R

Effect coding can be set in R using the contr.sum function:

X <- factor(c("A", "B", "C"))
contrasts(X) <- contr.sum(3) # set up effect coding in R
model <- lm(Y ~ X, data = mydata) # use linear regression to mimic ANOVA-style results
summary(model)

3.5 Exercise

library(ESRM64103)
head(ESRM64103::exp_political_attitude)

       party scores
1   Democrat      4
2   Democrat      3
3   Democrat      5
4   Democrat      4
5   Democrat      4
6 Republican      6

Under treatment coding, which level is the reference and how do you interpret the coefficients?
Under effect coding, what does the intercept represent and how do you interpret each group coefficient relative to the grand mean?

library(ESRM64103)
dat <- exp_political_attitude
dat$party <- factor(dat$party, levels = c("Democrat", 
"Republican", "Independent"))
# 1) Fit with treatment (dummy) coding
contrasts(dat$party) <- __________
fit_tr <- lm(__________)
# 2) Fit with effect (sum) coding
contrasts(dat$party) <- __________
fit_eff <- lm(__________)
# Compare ANOVA tables and coefficients
anova_tr <- anova(fit_tr)
coef_tr  <- summary(fit_tr)$coef
anova_eff <- anova(fit_eff)
coef_eff  <- summary(fit_eff)$coef
    
list(
  Treatment_ANOVA = anova_tr,
  Treatment_Coef  = round(coef_tr, 3),
  Effect_ANOVA    = anova_eff,
  Effect_Coef     = round(coef_eff, 3)
)

3.6 Self-defined contrast

Extended example 2: Assume now that the average of the STEM groups is different from the average of the non-STEM groups

3.6.1 Method 1: Calculation by Hand

  group  department  N   Mean      SD shapiro.test.p.values Contrasts
1    g1 Engineering 28 4.2500 3.15054             0.0775874   0.50000
2    g2   Education 28 2.7589 2.19478             0.0760542  -0.33333
3    g3   Chemistry 28 3.5446 2.86506             0.0062253   0.50000
4    g4   Political 28 3.8568 0.58325             0.0302312  -0.33333
5    g5  Psychology 28 2.0243 1.30911             0.0614743  -0.33333

$H_0: \frac{\mu_{Engineering}+\mu_{Chemistry}}{2} = \frac{\mu_{Education}+\mu_{PoliSci}+\mu_{Psychology}}{3}$

Weighted mean difference:

$C = c_1\mu_{Eng}+c_2\mu_{Edu}+c_3\mu_{Chem}+c_4\mu_{PoliSci}+c_5\mu_{Psych}\\ = \frac{1}{2}*4.25+(-\frac13)*2.75+(\frac12)*3.54+(-\frac13)*3.85+(-\frac13)*2.02\\ = 1.0173$

(C <- sum(summary_tbl_ext$Contrasts*summary_tbl_ext$Mean))

[1] 1.0173

$\sum\frac{c^2}{n} = \frac{(\frac12)^2}{28}+\frac{(-\frac13)^2}{28}+\frac{(\frac12)^2}{28}+\frac{(-\frac13)^2}{28}+\frac{(-\frac13)^2}{28}$

Sum_C2_n <- sum(summary_tbl_ext$Contrasts^2 / summary_tbl$N)
MSE = sum((residuals(aov(score ~ group, dt)))^2) / (nrow(dt) - 5)
t = as.numeric(C / sqrt(MSE * Sum_C2_n))
p.value = pt(t, df = 135, lower.tail = FALSE) * 2
data.frame(t, p.value)

       t   p.value
1 2.6369 0.0093476

$t = \frac{C}{\sqrt{MSE*\sum\frac{c^2}{n} }} = \frac{1.0173}{\sqrt{5.0011*0.029762}}=2.6368$

3.6.2 Method 2: Linear regression contrasts in R

# set first contrast
contrasts(dt$group) <- matrix(
  c(1/2, -1/3, 1/2, -1/3, -1/3)
)
fit_extended <- lm(score ~ group, dt)
unique(model.matrix(fit_extended))

    (Intercept)   group1    group2   group3   group4
1             1  0.50000 -0.692619 -0.10070 -0.10070
29            1 -0.33333  0.164436 -0.56552 -0.56552
57            1  0.50000  0.692619  0.10070  0.10070
85            1 -0.33333 -0.082218  0.78276 -0.21724
113           1 -0.33333 -0.082218 -0.21724  0.78276

summary(fit_extended)$coefficient[2,] |> round(3)

  Estimate Std. Error    t value   Pr(>|t|) 
     1.221      0.463      2.637      0.009

3.7 Exercise

(g1, g2, g3) vs (g4, g5)

contrasts(dt$group) <- matrix(
  __________________
)
fit <- lm(score ~ group, dt)
unique(_____)
summary(________)$coefficient[_,] |> round(3)

contrasts(dt$group) <- matrix(
  c(-1/3, -1/3, -1/3, 1/2, 1/2)
)
fit <- lm(score ~ group, dt)
unique(model.matrix(fit))

    (Intercept)   group1      group2    group3    group4
1             1 -0.33333 -4.3068e-01 -0.490501 -0.490501
29            1 -0.33333 -3.8540e-01  0.508987  0.508987
57            1 -0.33333  8.1608e-01 -0.018486 -0.018486
85            1  0.50000 -1.6593e-18  0.500000 -0.500000
113           1  0.50000 -1.6593e-18 -0.500000  0.500000

summary(fit)$coefficient[2,] |> round(3)

  Estimate Std. Error    t value   Pr(>|t|) 
    -0.693      0.463     -1.496      0.137

4 Effect Sizes

4.1 What Are Effect Sizes?

Effect size measures the magnitude of an effect beyond statistical significance.
- Put simply: a p-value is partially dependent on sample size and does not give us any insight into the strength of the relationship
- Lower p-value → can result from increasing sample size
Provides context for interpreting practical significance.
- In scientific experiments, it is often useful to know not only whether an experiment has a statistically significant effect, but also the size (magnitude) of any observed effects.
Common measures: Eta squared ( $\eta^2$ ), Omega squared ( $\omega^2$ ), Cohen’s d.

Note

Many psychology journals require the reporting of effect sizes

4.2 Eta squared

$\eta^2$ : Proportion of total variance explained by the independent variable.
Formula: $\eta^2 = \frac{SS_{Model}}{SS_{Total}}$
Interpretation:
- Small: 0.01, Medium: 0.06, Large: 0.14

(F_table <- as.data.frame(anova(fit)))

           Df  Sum Sq Mean Sq F value    Pr(>F)
group       4  89.368 22.3420  4.4674 0.0020173
Residuals 135 675.149  5.0011      NA        NA

(eta_2 <- F_table$`Sum Sq`[1] / sum(F_table$`Sum Sq`))

[1] 0.11689

Interpretation: 11.69% of variance in the DV is due to group differences.

4.3 Drawbacks of eta squared

As you add more variables to the model, the proportion explained by any one variable will automatically decrease.
- This makes it hard to compare the effect of a single variable in different studies.
- Partial Eta Squared solves this problem. There, the denominator is not the total variation in Y, but the unexplained variation in Y plus the variation explained just by that IV.
  - Any variation explained by other IVs is removed from the denominator.
- In a one-way ANOVA, Eta Squared and Partial Eta Squared will be equal, but this isn’t true in models with more than one independent variable (factorial ANOVA).
Eta Squared is a biased measure of population variance explained (although it is accurate for the sample).
- It always overestimates it. This bias gets very small as sample size increases, but for small samples an unbiased effect size measure is Omega Squared.

4.4 Omega squared

Omega Squared ( $\omega^2$ ) has the same basic interpretation but uses unbiased measures of the variance components.
- Because it is an unbiased estimate of population variances, omega squared is always smaller than eta squared.
Unbiased estimate of effect size, preferred for small samples.
Formula: $\omega^2 = \frac{SS_{Model} - df_{Model} \cdot MSE}{SS_{Total} + MSE}$
Interpretation follows $\eta^2$ scale but slightly smaller values.

F_table

           Df  Sum Sq Mean Sq F value    Pr(>F)
group       4  89.368 22.3420  4.4674 0.0020173
Residuals 135 675.149  5.0011      NA        NA

1attach(F_table)
2(Omega_2 <- (`Sum Sq`[1] - Df[1] * MSE) / (sum(`Sum Sq`) + MSE))
detach(F_table)

1: Attach data set so that you can directly call the columns without “$”
2: The formula of Omega square

[1] 0.090139

4.5 Effect Size for Planned Contrasts

Correlation-based effect size: $r = \sqrt{\frac{t^2}{t^2 + df}} = \sqrt{\frac{F}{F + df}}$
Example: For $t = 2.49, df = 135$ : $r = \sqrt{\frac{2.49^2}{2.49^2 + 135}} = 0.21$
- Small to moderate effect.

(coef_tbl <- as.data.frame(summary(fit)$coefficients))

               Estimate Std. Error    t value   Pr(>|t|)
(Intercept)  3.28692857    0.18900 17.3908736 2.8188e-36
group1      -0.69278571    0.46296 -1.4964232 1.3688e-01
group2      -0.00096962    0.42262 -0.0022943 9.9817e-01
group3       0.17035380    0.42262  0.4030862 6.8752e-01
group4      -1.66214620    0.42262 -3.9329222 1.3359e-04

attach(coef_tbl)
round(sqrt(`t value`^2 / (`t value`^2 + 135)), 3)
detach(coef_tbl)

[1] 0.831 0.128 0.000 0.035 0.321

Shows a small to moderate positive relationship between g1 and g5.

4.6 Cohen’s d and Hedges’ g

Used for simple mean comparisons.
Cohen’s d formula: $d = \frac{M_1 - M_2}{SD_{pooled}}$
Hedges’ g corrects for small sample bias.
Guidelines:
- Small: 0.2, Medium: 0.5, Large: 0.8

4.7 Guidelines for Effect Size

For our example: there is a significant effect of academic program on growth mindset scores (F(4, 135) = 4.47).
Academic program explains 11.69% of the variance in growth mindset scores. This is a medium to large effect ( $\eta^2$ = 0.1169).

4.8 Summary

Planned contrasts allow hypothesis-driven mean comparisons.
Orthogonal contrasts maintain Type I error control.
Effect sizes help interpret the importance of results.
Combining planned contrasts with effect size measures enhances statistical analysis.

--- title: "Lecture 05: ANOVA Comparisons and Contrasts" subtitle: "Experimental Design in Education" date: "2025-02-17" execute: eval: true echo: true message: false warning: false format: html: code-tools: true code-line-numbers: false code-fold: false number-offset: 1 fig.width: 10 fig-align: center message: false grid: sidebar-width: 350px uark-revealjs: chalkboard: true embed-resources: false code-fold: false number-sections: true number-depth: 1 footer: "ESRM 64503" slide-number: c/t tbl-colwidths: auto scrollable: true output-file: slides-index.html mermaid: theme: forest filters: - webr --- ::: objectives [Class Outline]{.redcolor .bigger} - Planned Contrast - Example: Group Comparison - STEM vs. Non-STEM Groups - ANOVA-style t statistics and regression-style coding in R - Effect sizes ::: # Planned Contrasts ## Definition: Planned Contrasts - Pre-defined: - Unlike post hoc tests, planned contrasts are determined **before** looking at the data, meaning the researcher has a specific [[hypothesis about which groups to compare]{.underline}]{.bigger}. - **Definition**: Planned contrasts are hypothesis-driven comparisons made before data collection. - Weights assigned: - To perform a planned contrast, each group is assigned a numerical "weight" which reflects its role in the comparison, with the [weights usually summing to zero]{.underline}. $$ D = weights * means $$ - **D**: unscaled group differences given coding scheme ## Example of Planned Contrasts - **Imagine** a study comparing the effects of three different study methods (**A, B, C**) on test scores. - One planned contrast might be to compare the average score of method A (considered the "experimental" method) against the combined average of methods B and C (considered the "control" conditions), - Testing the hypothesis that method A leads to significantly higher scores than the traditional methods. - $H_0: \mu_{A} = \frac{\mu_B+\mu_C}{2}$; we also call this a **complex contrast** - **When** to use planned contrasts: - When you have a clear theoretical basis for predicting specific differences between groups in your study. - When you are only interested in a few specific comparisons, not all possible pairwise comparisons. ## What Does Each Contrast Tell Us? - Each contrast is a mean comparison (via t-test). - **Simple** contrast (pairwise) compares two individual group means. - **Complex** contrast compares a combination of group means. - Must be theoretically justified for meaningful interpretation. ## Simple vs. Complex Comparisons - **Simple Comparison:** Two groups directly compared. - Example: $H_0: \mu_2 = \mu_3$ - **Complex Comparison:** Combines means of multiple groups. - Example: $H_0: \frac{(\mu_1 + \mu_2)}{2} = \mu_3$ - Example: $H_0: \frac{(\mu_1 + \mu_2 + \mu_3)}{3} = \frac{(\mu_4 + \mu_5)}{2}$ ::: callout-note We should not test all possible combinations of groups. Instead, justify your comparison plan before performing statistic analysis. ::: ## Today's focus: Complex Comparisons - We performed omnibus tests in the last lecture, which provide all pairwise group comparisons (simple contrasts) - Today we focus more on **complex contrasts**. - **Helmert** contrast: Compares each mean to the mean of subsequent groups. - **Sum (deviation)** contrast: each group compared to the grand mean - **Polynomial** contrast: Tests for trends in ordered data. - By default, R uses **treatment** contrasts: each group compared to the reference group - Question: is "treatment contrast" a simple or complex contrast? - G1 vs. Treatment (reference) - G2 vs. Treatment (reference) - .... ## Orthogonal vs. Non-Orthogonal Contrasts - **Orthogonal Contrasts:** Independent from each other, sum of product of weights equals zero. - **Non-Orthogonal Contrasts:** Not independent, lead to inflated Type I error rates. ::: callout-note Orthogonal contrasts allow clear interpretation without redundancy. ::: - Orthogonal contrasts follow a series of group comparisons that do not overlap variances. ## Orthogonal contrasts from variances: no redundancy [Helmert contrast example]{.redcolor} - With a logical control group, a good first contrast compares all treatment groups to the one control group. - To get each level of the IV alone, you should have one fewer contrast than your number of IV levels (3 levels = 2 contrasts) - Once an IV level appears by itself, it shouldn’t reappear in subsequent contrasts ```{dot} //| echo: FALSE //| fig-width: 15 digraph VarianceDiagram { rankdir=LR; A [label="Total Variance Explained", shape=box, style=filled, fillcolor=tomato, fontcolor=blue, fontsize=20, width=2.5, height=1.5]; B [label="Variance for G1 and G2", shape=box, style=filled, fillcolor=lightgrey, fontcolor=red, fontsize=16, width=2.0, height=1]; C [label="Variance for G1", shape=box, style=filled, fillcolor=lightblue, fontcolor=red]; D [label="Variance for G2", shape=box, style=filled, fillcolor=lightblue, fontcolor=red]; E [label="Variance for G3", shape=box, style=filled, fillcolor=lightgrey, fontcolor=red]; A -> B; B -> C; B -> D; A -> E; } ``` ## Example of Orthogonal Contrasts ::::: columns ::: {.column width="50%"} - Contrast 1: g3 vs. (g1, g2) ![](image/g3vsg12.jpg){width="90%"} ::: ::: {.column width="50%"} - Contrast 2: g1 vs. g2 ![](image/g1vsg2.jpg){width="90%"} ::: ::::: ## Orthogonal Planned Contrasts - If the same exact combination of means is not found in more than one contrast, the contrasts are independent (orthogonal) - Check this by ensuring that the product of the weights across all contrasts sums to zero - For an orthogonal comparison, **contrasts are independent of each other:** - We weight the means included on each side of the contrast - Each contrast has weights that sum to 0 - Groups not in the contrast get a weight of 0 - Why does independence matter? - [Type I error rate is unaffected by independent (orthogonal) contrasts]{.underline} - Interpretation of contrasts is cleaner because contrasts aren’t related (you’ve isolated effects) | Group | Contrast 1 | Contrast 2 | Product | |-------|------------|------------|---------| | G1 | +1 | -1 | -1 | | G2 | +1 | +1 | +1 | | G3 | -2 | 0 | 0 | | Sum | 0 | 0 | 0 | ## Contrasts' Independence checking in R ```{r} #| results: hold #| code-fold: true cat("## Constrasts and Coding \n") contras <- matrix( c(1, 1, -2, -1, 1, 0), ncol = 2 ) contras ``` To understand if these contrasts are orthogonal, one could compute the dot product between the contrast vectors. The dot product of first and second contrast is: $(1 * -1) + (1 * 1) + (-2 * 0) = -1 + 1 + 0 = 0$ Since the dot product of the contrasts is equal to zero, these contrasts are indeed orthogonal. Orthogonal contrasts have the advantage that the comparison of the means of one contrast does not affect the comparison of the means of any other contrast. ```{r} #| results: hold #| code-fold: true cat("## the dot-product of two contrasts should be zero\n") t(contras[,1]) %*% contras[,2] cat("## if non-diagonal elements are zero and diagnoal elements = 0 \n") crossprod(contras) ``` ::: panel-tabset ### Exercise ```{webr-r} new_code <- matrix(c( 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1 ), byrow = TRUE, nrow = 5) new_code ``` ### Answer ```{webr-r answer} new_code <- matrix(c( 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1 ), byrow = TRUE, nrow = 5) cat("## the cross-product of two contrasts should be zero\n") t(new_code[,1]) %*% new_code[,2] t(new_code[,1]) %*% new_code[,3] t(new_code[,1]) %*% new_code[,4] t(new_code[,2]) %*% new_code[,3] cat("## if non-diagonal elements are zero and diagnoal elements = 0 \n") crossprod(new_code) ``` ::: ## Computing Planned Contrasts - Formula for contrast value: $C = c_1\mu_1 + c_2\mu_2 + \dots + c_k\mu_k$ - Test statistic: $t = \frac{C}{\sqrt{MSE \sum \frac{c_i^2}{n_i}}}$ - $MSE$: Mean Square Error from ANOVA - $c_i$: Contrast coefficients - $n_i$: Sample size per group ## Planned Contrasts and Coding Scheme - The relationship between planned contrasts in ANOVA and coding in regression lies in how categorical variables are represented and interpreted in statistical models. - Both approaches aim to test specific hypotheses about group differences, but their implementation varies based on the framework: - ANOVA focuses on partitioning variance, - Regression interprets categorical predictors through coding schemes. ------------------------------------------------------------------------ ### Planned Contrast and Linear Regression - Planned contrast can be done using `linear regression` + `contrasts` - Let's look at the default contrasts plan: treatment contrasts == dummy coding ```{r} #| output-location: column library(tidyverse) library(kableExtra) library(here) set.seed(42) dt <- read.csv(here("teaching/2025-01-13-Experiment-Design/Lecture05","week5_example.csv")) ## Treatment contrast matrix attributes(C(dt$group, treatment, 4))$contrasts ``` ```{r} #| output-location: column ## Sum contrast matrix attributes(C(dt$group, sum, 4))$contrasts ``` ```{r} #| output-location: column attributes(C(dt$group, helmert, 4))$contrasts ``` ```{r} #| output-location: column crossprod(attributes(C(dt$group, treatment, 4))$contrasts) ``` # Example - STEM vs. Non-STEM Groups ## Background - Hypothesis: STEM students have different growth mindset scores than non-STEM students. - Weights assigned: - STEM (Engineering, Chemistry): $+\frac{1}{2}$ - Non-STEM (Education, Political Science, Psychology): $-\frac{1}{3}$ - Compute contrast value and test using t-statistic. ## Set Contrasts in R ```{r} #| code-fold: true #| results: hide # Set seed for reproducibility options(digits = 5) summary_tbl <- dt |> group_by(group) |> summarise( N = n(), Mean = mean(score), SD = sd(score), shapiro.test.p.values = shapiro.test(score)$p.value ) |> mutate(department = c("Engineering", "Education", "Chemistry", "Political", "Psychology")) |> relocate(group, department) summary_tbl ``` ```{r} #| echo: false kable(summary_tbl) ``` - Homogeneity of variance assumption: Levene's test ----------------------- ### Stem vs. Non-Stem ```{r} #| code-fold: show aov_fit <- aov(score ~ group, data = dt) car::leveneTest(aov_fit) ``` Even though the assumption checks did not pass using the original categorical levels, we may still be interested in different group contrasts. ----------------------- ```{r} #| code-fold: true coef(lm(score ~ group, data = dt)) ``` $$ (\beta_0 + \beta_0 + \beta_2) / 2 - (\beta_0 + \beta_1 + \beta_0 + \beta_3 + \beta_0 + \beta_4) / 3 \\ = \beta_2 / 2 - (\beta_1 + \beta_3 + \beta_3) /3 $$ ```{r} #| code-fold: true contrast <- matrix( c(0, -1/3, 1/2, -1/3, -1/3), nrow = 1 ) rownames(contrast) <- "Stem vs. Non-Stem" summary(multcomp::glht(aov_fit, linfct = contrast)) ``` ```{r} means <- summary_tbl$Mean mean(means[c(1,3)]) - mean(means[-c(1,3)]) ``` ------------------------------------------------------------------------ ### Complex Contrast Matrix - There are multiple "canned" contrasts: Helmert, sum (effect coding), treatment [For example, Helmert: four contrasts]{.redcolor} 1. g1 vs. g2: $\mu_{Engineering} = \mu_{Education}$ 2. $\frac{g1+g2}{2}$ vs. g3: $\mu_{non-Chemistry} = \mu_{Chemistry}$ 3. $\frac{g1+g2+g3}{3}$ vs. g4: $\mu_{non-Political} = \mu_{Political}$ 4. $\frac{g1+g2+g3+g4}{4}$ vs. g5: $\mu_{non-Psychology} = \mu_{Psychology}$ Summary Statistics: ```{r} #| code-fold: true #| results: hold dt$group <- factor(dt$group, levels = c("g1", "g2", "g3", "g4", "g5")) groups <- levels(dt$group) cH <- contr.helmert(groups) # pre-defined four contrasts colnames(cH) <- paste0("Ctras", 1:4) summary_ctras_tbl <- cbind(summary_tbl, cH) summary_ctras_tbl ``` ------------------------------------------------------------------------ ### Helmert contrasts are Orthogonal ```{r} apply(cH, 2, sum) crossprod(cH) # diagonal -- columns are orthogonal ``` ```{r} summary(aov(score ~ group, dt)) ``` ## ANOVA: t-value formula for Defined Contrast Matrix $t = \frac{C}{\sqrt{MSE \sum \frac{c_i^2}{n_i}}}$ ```{r} #| results: hold Sum_C2_n <- colSums(cH^2 / summary_tbl$N) C <- crossprod(summary_tbl$Mean, cH) MSE <- 5 t <- as.numeric(C / sqrt(MSE * Sum_C2_n)) tibble( t_value = t, p_value = pt(t, df = 135) ## p-values ) ``` - g1 vs. g2: We reject the null and determine that the mean of Education is different from the mean of Engineering in their growth mindset scores (p = 0.531). - $\frac{g1+g2}{2}$ vs. g3: We retain the null and determine that the mean of Chemistry is not significantly different from the mean of Education and Engineering in their growth mindset scores (p = 0.531). ------------------------------------------------------------------------ ### Helmert Contrast Remember the planned contrast: g1 vs. g2 from the Helmert contrast: - t-value: -2.495 - p-value: 0.0069 - df: 134 ```{r} #| code-fold: true #| results: hold contrasts(dt$group) <- "contr.helmert" fit_helmert <- lm(score ~ group, dt) contr.helmert(levels(dt$group)) summary(fit_helmert)$coefficients ``` ```{r} mean(summary_tbl$Mean) unique(cbind(model.matrix(fit_helmert), group = dt$group)) summary_tbl$Mean %*% contrasts(dt$group) / c(2, 6, 12, 20) ``` ------------------------------------------------------------------------ ### Treatment contrasts ::::: columns ::: {.column width="50%"} - For treatment contrasts, four dummy variables are created to compare: - G1 (ref) vs. G2 - G1 (ref) vs. G3 - G1 (ref) vs. G4 - G1 (ref) vs. G5 ::: ::: {.column width="\"50%"} - `Intercept`: G1's mean - `group2`: G2 vs. G1 - `group3`: G3 vs. G1 - `group4`: G4 vs. G1 - `group5`: G5 vs. G1 ::: ::::: ```{r} #| output-location: column library(multcomp) contrasts(dt$group) <- "contr.treatment" fit <- lm(score ~ group, dt) unique(cbind(model.matrix(fit), group = dt$group)) ``` ```{r} #| output-location: column summary(fit)$coefficients ``` ------------------------------------------------------------------------ ### Sum Contrasts - Another type of coding is **effect coding**. In R, the corresponding contrast type is the so-called **sum contrasts**. - A detailed post about sum contrasts can be found [here](https://learnb4ss.github.io/learnB4SS/articles/contrasts.html) - With sum contrasts, the reference level is the grand mean. - $\frac{g1+g2+g3+g4+g5}{5}$ vs. g1/g2/g3/g4: the difference between mean score of g1 with grand mean across all five groups ```{r} #| output-location: column contrasts(dt$group) <- "contr.sum" fit2 <- lm(score ~ group, dt) contr.sum(levels(dt$group)) ``` ```{r} #| output-location: column summary(fit2)$coefficients ``` ```{r} #| output-location: column mean(dt$score) # (Intercept) grand mean ``` ```{r} #| output-location: column tibble( Label = paste0("group", 1:4), Estimate = summary_tbl$Mean[1:4] - mean(dt$score) ) ``` ## Effect Coding (Deviation Coding) - In modern statistics, regression-style coding is statistically equivalent to an ANOVA-style contrast matrix. - Equivalent to ANOVA-style contrasts. ([we will use this in R to reproduce ANOVA-style contrast matrix]{.redcolor}) - Compares each level to the grand mean. ::: callout-note Effect coding is a method of encoding categorical variables in regression models, similar to dummy coding, but with a different interpretation of the resulting coefficients. It is particularly useful when researchers want to compare each level of a categorical variable to the overall mean rather than to a specific reference category. ::: ------------------------------------------------------------------------ ### Definition and Representation In effect coding, categorical variables are transformed into numerical variables, typically using values of -1, 0, and 1. The key difference from dummy coding is that the reference category is represented by -1 instead of 0, and the coefficients indicate deviations from the grand mean. For a categorical variable with `k` levels, effect coding requires `k-1` coded variables. If we have a categorical variable X with three levels: $A, B, C$, the effect coding scheme could be: | Category | $X_1$ | $X_2$ | |---------------|-------|-------| | A | 1 | 0 | | B | 0 | 1 | | C (reference) | -1 | -1 | The last category ($C$) is the reference group, coded as -1 for all indicator ------------------------------------------------ - In effect coding, one level of the categorical variable (usually the last one or the one that is considered a 'reference' category) is coded as -1, and the others are coded as +1 or 0, suggesting their relation to the overall mean, rather than to a specific category. Here's how you could effect code a categorical variable with three levels (e.g., groups): ``` [,1] [,2] g1 1 0 g2 0 1 g3 -1 -1 ``` This coding scheme shows that: - Group1 (`g1`) is compared to the overall effect by coding it as `1` in the first column and `0` in the second, suggesting it is above or below the overall mean. - Conversely, Group2 (`g2`) is also contrasted with the overall effect. - Group3 (`g3`), coded as `-1` in both columns, serves as the reference category against which the other two are compared. In your regression output, the coefficients for the first two groups will show how the mean of these groups differs from the overall mean. The intercept will represent the overall mean across all groups. ------------------------------------------------------------------------ ### Interpretation of Coefficients in effective coding ::::: columns ::: {.column width="40%"} When effect coding is used in a regression model: $$ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \epsilon $$ - $X_1$ and $X_2$ are coded variables. They have no much meaning, but their coefficients are **important** - $\beta_0$ represents the **grand mean** of $Y$ across all categories. - $\beta_1$ and $\beta_2$ represent the **deviation** of categories $A$ and $B$ from the grand mean. - The reference group ($C$) does not have a separate coefficient; instead, its deviation can be inferred as $-(\beta_1 + \beta_2)$. ::: ::: {.column width="60%"} ```{r} #| code-fold: true #| eval: true #| fig-height: 8 library(ggplot2) # Create a data frame for text labels text_data <- data.frame( x = rep(0.25, 3), # Repeating the same x-coordinate y = c(0.3, 0.7, 0.9), # Different y-coordinates label = c("C: beta[0] - beta[1] - beta[2]", "A: beta[0] + 1*'×'*beta[1] + 0*'×'*beta[2]", "B: beta[0] + 0*'×'*beta[1] + 1*'×'*beta[2]") # Labels ) # Create an empty ggplot with defined limits ggplot() + geom_text(data = text_data, aes(x = x, y = y, label = label), parse = TRUE, size = 11) + # Add a vertical line at x = 0.5 # geom_vline(xintercept = 0.5, color = "blue", linetype = "dashed", linewidth = 1) + # Add two horizontal lines at y = 0.3 and y = 0.7 geom_hline(yintercept = c(0.35, 0.75, 0.95), color = "red", linetype = "solid", linewidth = 1) + geom_hline(yintercept = 0.5, color = "grey", linetype = "solid", linewidth = 1) + geom_text(aes(x = .25, y = .45, label = "grand mean of Y"), color = "grey", size = 11) + # Set axis limits xlim(0, 1) + ylim(0, 1) + labs(y = "Y", x = "") + # Theme adjustments theme_minimal() + theme(text = element_text(size = 20)) ``` ::: ::::: ------------------------------------------------------------------------ ### Comparison to Dummy Coding - **Dummy Coding**: Compares each category to a **specific reference category** (e.g., comparing A and B to C). | Category | $X_1$ | $X_2$ | |---------------|-------|-------| | A | 1 | 0 | | B | 0 | 1 | | C (reference) | 0 | 0 | - **Effect Coding**: Compares each category to the **grand mean** rather than a single reference category. ------------------------------------------------------------------------ ### Use Cases Effect coding is beneficial when: - There is no natural baseline category, and comparisons to the **overall mean** are more meaningful. - Researchers want to maintain **sum-to-zero constraints** for categorical variables in linear models. - In ANOVA-style analyses, where main effects and interaction effects are tested under an equal-weight assumption. ------------------------------------------------------------------------ ### Implementation in R Effect coding can be set in R using the `contr.sum` function: ```{r} #| eval: false X <- factor(c("A", "B", "C")) contrasts(X) <- contr.sum(3) # set up effect coding in R model <- lm(Y ~ X, data = mydata) # use linear regression to mimic ANOVA-style results summary(model) ``` ## Exercise ```{r} library(ESRM64103) head(ESRM64103::exp_political_attitude) ``` ::: panel-tabset ### Quiz 1. Under treatment coding, which level is the reference and how do you interpret the coefficients? 2. Under effect coding, what does the intercept represent and how do you interpret each group coefficient relative to the grand mean? ```{webr-r} library(ESRM64103) dat <- exp_political_attitude dat$party <- factor(dat$party, levels = c("Democrat", "Republican", "Independent")) # 1) Fit with treatment (dummy) coding contrasts(dat$party) <- __________ fit_tr <- lm(__________) # 2) Fit with effect (sum) coding contrasts(dat$party) <- __________ fit_eff <- lm(__________) # Compare ANOVA tables and coefficients anova_tr <- anova(fit_tr) coef_tr <- summary(fit_tr)$coef anova_eff <- anova(fit_eff) coef_eff <- summary(fit_eff)$coef list( Treatment_ANOVA = anova_tr, Treatment_Coef = round(coef_tr, 3), Effect_ANOVA = anova_eff, Effect_Coef = round(coef_eff, 3) ) ``` ### Answer ```{webr-r} library(ESRM64103) dat <- exp_political_attitude dat$party <- factor(dat$party, levels = c("Democrat", "Republican", "Independent")) # 1) Fit with treatment (dummy) coding contrasts(dat$party) <- contr.treatment(3) fit_tr <- lm(scores ~ party, data = dat) # 2) Fit with effect (sum) coding contrasts(dat$party) <- contr.sum(3) fit_eff <- lm(scores ~ party, data = dat) # Compare ANOVA tables and coefficients anova_tr <- anova(fit_tr) coef_tr <- summary(fit_tr)$coef anova_eff <- anova(fit_eff) coef_eff <- summary(fit_eff)$coef list( Treatment_ANOVA = anova_tr, Treatment_Coef = round(coef_tr, 3), Effect_ANOVA = anova_eff, Effect_Coef = round(coef_eff, 3) ) ``` ::: ## Self-defined contrast - Extended example 2: Assume now that the average of the STEM groups is different from the average of the non-STEM groups ### Method 1: Calculation by Hand ```{r} #| echo: false (summary_tbl_ext <- cbind(summary_tbl, Contrasts = c(1/2, -1/3, 1/2, -1/3, -1/3))) ``` $$ H_0: \frac{\mu_{Engineering}+\mu_{Chemistry}}{2} = \frac{\mu_{Education}+\mu_{PoliSci}+\mu_{Psychology}}{3} $$ Weighted mean difference: $$ C = c_1\mu_{Eng}+c_2\mu_{Edu}+c_3\mu_{Chem}+c_4\mu_{PoliSci}+c_5\mu_{Psych}\\ = \frac{1}{2}*4.25+(-\frac13)*2.75+(\frac12)*3.54+(-\frac13)*3.85+(-\frac13)*2.02\\ = 1.0173 $$ ```{r} (C <- sum(summary_tbl_ext$Contrasts*summary_tbl_ext$Mean)) ``` $$ \sum\frac{c^2}{n} = \frac{(\frac12)^2}{28}+\frac{(-\frac13)^2}{28}+\frac{(\frac12)^2}{28}+\frac{(-\frac13)^2}{28}+\frac{(-\frac13)^2}{28} $$ ```{r} Sum_C2_n <- sum(summary_tbl_ext$Contrasts^2 / summary_tbl$N) MSE = sum((residuals(aov(score ~ group, dt)))^2) / (nrow(dt) - 5) t = as.numeric(C / sqrt(MSE * Sum_C2_n)) p.value = pt(t, df = 135, lower.tail = FALSE) * 2 data.frame(t, p.value) ``` $$ t = \frac{C}{\sqrt{MSE*\sum\frac{c^2}{n} }} = \frac{1.0173}{\sqrt{5.0011*0.029762}}=2.6368 $$ ------------------------------------------------------------------------ ### Method 2: Linear regression contrasts in R ```{r} # set first contrast contrasts(dt$group) <- matrix( c(1/2, -1/3, 1/2, -1/3, -1/3) ) fit_extended <- lm(score ~ group, dt) unique(model.matrix(fit_extended)) ``` ```{r} #| output-line-numbers: "3" #| class-output: highlight summary(fit_extended)$coefficient[2,] |> round(3) ``` ## Exercise ::: panel-tabset ### Quiz (g1, g2, g3) vs (g4, g5) ```{r} #| eval: false contrasts(dt$group) <- matrix( __________________ ) fit <- lm(score ~ group, dt) unique(_____) summary(________)$coefficient[_,] |> round(3) ``` ### Anwser ```{r} contrasts(dt$group) <- matrix( c(-1/3, -1/3, -1/3, 1/2, 1/2) ) fit <- lm(score ~ group, dt) unique(model.matrix(fit)) summary(fit)$coefficient[2,] |> round(3) ``` ::: # Effect Sizes ## What Are Effect Sizes? - Effect size measures the magnitude of an effect beyond statistical significance. - Put simply: a p-value is partially dependent on sample size and does not give us any insight into the strength of the relationship - Lower p-value → can result from increasing sample size - Provides context for interpreting practical significance. - In scientific experiments, it is often useful to know not only whether an experiment has a statistically significant effect, but also the size (magnitude) of any observed effects. - Common measures: Eta squared ($\eta^2$), Omega squared ($\omega^2$), Cohen’s d. ::: callout-note Many psychology journals require the reporting of effect sizes ::: ## Eta squared - $\eta^2$: Proportion of total variance explained by the independent variable. - Formula: $\eta^2 = \frac{SS_{Model}}{SS_{Total}}$ - Interpretation: - Small: 0.01, Medium: 0.06, Large: 0.14 ```{r} (F_table <- as.data.frame(anova(fit))) ``` ```{r} (eta_2 <- F_table$`Sum Sq`[1] / sum(F_table$`Sum Sq`)) ``` [Interpretation: 11.69% of variance in the DV is due to group differences.]{.redcolor} ## Drawbacks of eta squared 1. **As you add more variables to the model, the proportion explained by any one variable will automatically decrease.** - This makes it hard to compare the effect of a single variable in different studies. - Partial Eta Squared solves this problem. There, the denominator is not the total variation in Y, but the unexplained variation in Y plus the variation explained just by that IV. - Any variation explained by other IVs is removed from the denominator. - In a one-way ANOVA, Eta Squared and Partial Eta Squared will be equal, but this isn’t true in models with more than one independent variable (factorial ANOVA). 2. **Eta Squared is a biased measure of population variance explained (although it is accurate for the sample).** - It always overestimates it. This bias gets very small as sample size increases, but for small samples an unbiased effect size measure is Omega Squared. ## Omega squared - Omega Squared ($\omega^2$) has the same basic interpretation but uses unbiased measures of the variance components. - Because it is an unbiased estimate of population variances, omega squared is always smaller than eta squared. - Unbiased estimate of effect size, preferred for small samples. - Formula: $\omega^2 = \frac{SS_{Model} - df_{Model} \cdot MSE}{SS_{Total} + MSE}$ - Interpretation follows $\eta^2$ scale but slightly smaller values. ```{r} F_table ``` ```{r} #| results: hold attach(F_table) #<1> (Omega_2 <- (`Sum Sq`[1] - Df[1] * MSE) / (sum(`Sum Sq`) + MSE)) #<2> detach(F_table) ``` 1. Attach data set so that you can directly call the columns without "\$" 2. The formula of Omega square ## Effect Size for Planned Contrasts - Correlation-based effect size: $r = \sqrt{\frac{t^2}{t^2 + df}} = \sqrt{\frac{F}{F + df}}$ - Example: For $t = 2.49, df = 135$: $r = \sqrt{\frac{2.49^2}{2.49^2 + 135}} = 0.21$ - Small to moderate effect. ```{r} (coef_tbl <- as.data.frame(summary(fit)$coefficients)) ``` ```{r} #| results: hold attach(coef_tbl) round(sqrt(`t value`^2 / (`t value`^2 + 135)), 3) detach(coef_tbl) ``` - Shows a small to moderate positive relationship between g1 and g5. ## Cohen’s d and Hedges’ g - Used for simple mean comparisons. - Cohen’s d formula: $d = \frac{M_1 - M_2}{SD_{pooled}}$ - Hedges’ g corrects for small sample bias. - Guidelines: - Small: 0.2, Medium: 0.5, Large: 0.8 ## Guidelines for Effect Size ![](images/clipboard-1482285190.png) - For our example: there is a significant effect of academic program on growth mindset scores (F(4, 135) = 4.47). - Academic program explains 11.69% of the variance in growth mindset scores. This is a medium to large effect ($\eta^2$ = 0.1169). ## Summary - Planned contrasts allow hypothesis-driven mean comparisons. - Orthogonal contrasts maintain Type I error control. - Effect sizes help interpret the importance of results. - Combining planned contrasts with effect size measures enhances statistical analysis.

webR Code Links

2 Planned Contrasts

2.1 Definition: Planned Contrasts

2.2 Example of Planned Contrasts

2.3 What Does Each Contrast Tell Us?

2.4 Simple vs. Complex Comparisons

2.5 Today’s focus: Complex Comparisons

2.6 Orthogonal vs. Non-Orthogonal Contrasts

2.7 Orthogonal contrasts from variances: no redundancy

2.8 Example of Orthogonal Contrasts

2.9 Orthogonal Planned Contrasts

2.10 Contrasts’ Independence checking in R

2.11 Computing Planned Contrasts

2.12 Planned Contrasts and Coding Scheme

2.12.1 Planned Contrast and Linear Regression

3 Example - STEM vs. Non-STEM Groups

3.1 Background

3.2 Set Contrasts in R

3.2.1 Stem vs. Non-Stem

3.2.2 Complex Contrast Matrix

3.2.3 Helmert contrasts are Orthogonal

3.3 ANOVA: t-value formula for Defined Contrast Matrix

3.3.1 Helmert Contrast

3.3.2 Treatment contrasts

3.3.3 Sum Contrasts

3.4 Effect Coding (Deviation Coding)

3.4.1 Definition and Representation

3.4.2 Interpretation of Coefficients in effective coding

3.4.3 Comparison to Dummy Coding

3.4.4 Use Cases

3.4.5 Implementation in R

3.5 Exercise

3.6 Self-defined contrast

3.6.1 Method 1: Calculation by Hand

3.6.2 Method 2: Linear regression contrasts in R

3.7 Exercise

4 Effect Sizes

4.1 What Are Effect Sizes?

4.2 Eta squared

4.3 Drawbacks of eta squared

4.4 Omega squared

4.5 Effect Size for Planned Contrasts

4.6 Cohen’s d and Hedges’ g

4.7 Guidelines for Effect Size

4.8 Summary

R History Command Contents