severe_contrasts.R

#### OVERVIEW ####

set.seed(432432) # For reproducibility

# For all examples, I'm using the below sampling error sd and n per cell
err <- 1
n_per_cell <- 200
# I'm also assuming the following:
# (1) our smallest effect of substantive interest is .4
# (2) all contrasts are Helment (successively compare one group to the average of the others)
# (3) all contrasts are unit-weighted

# This method should generalize to other orthogonal contrasts
# The parameter estimates for unit-weighted Helmert contrasts, however, have a 
# nice interpretation as a series of mean differences and are therefore easy to understand

library(tidyverse) # For mutate()
library(car) # For lht()


#### THREE GROUP ####

# Ceofficients for this case. Modify as desired
# First is the intercept, second is the focal contrast, remainder are the residual contrasts
coefs <- c(0, .7, .2)

# Create the data
# The last line creates the outcome using the coefficients above and the desired sampling error, err
# %*% is matrix multiplication
dat_three <- data.frame(group = rep(paste("group", 1:3), n_per_cell))
dat_three <- mutate(dat_three,
                    c1 = case_when(group=="group 1" ~ 2/3, TRUE ~ -1/3),
                    c2 = case_when(group=="group 1" ~ 0, group=="group 2" ~ 1/2, TRUE ~ -1/2),
                    outcome = as.vector(cbind(1, c1, c2) %*% coefs + rnorm(nrow(dat_three), sd=err)))

m_three <- lm(outcome ~ c1 + c2, data=dat_three)

### Test of the focal contrast against 0 ###
summary(m_three)
### We also get the test of the residual contrast against 0 from this model ###

### Equivalence tests ###
# This method allows for some "crud" in the residual contrast
lht(m_three, "c2=.4")
lht(m_three, "c2=-.4")
# We are permitted to divide the p-values by 2 because these are one-sided tests 
lht(m_three, "c2=.4")[2, 6]/2
lht(m_three, "c2=-.4")[2, 6]/2


#### FOUR GROUP ####

# Ceofficients for this case. Modify as desired
# First is the intercept, second is the focal contrast, remainder are the residual contrasts
coefs <- c(0, .7, .2, .1)

# Create the data
# The last line creates the outcome using the coefficients above and the desired sampling error, err
# %*% is matrix multiplication
dat_four <- data.frame(group = rep(paste("group", 1:4), n_per_cell))
dat_four <- mutate(dat_four,
                    c1 = case_when(group=="group 1" ~ -3/4, TRUE ~ 1/4),
                    c2 = case_when(group=="group 1" ~ 0, group=="group 2" ~ 2/3, TRUE ~ -1/3),
                    c3 = case_when(group %in% c("group 1", "group 2") ~ 0, group=="group 3" ~ 1/2, TRUE ~ -1/2),
                    outcome = as.vector(cbind(1, c1, c2, c3) %*% coefs + rnorm(nrow(dat_four), sd=err)))

m_four <- lm(outcome ~ c1 + c2 + c3, data=dat_four)

### Test of the focal contrast against 0 ###
summary(m_four)

### Test of the residual contrasts against 0 ###
# It probably makes sense to do a joint test of all of them in this case
# because there are multiple residual contrasts rather than one
lht(m_four, c("c2=0", "c3=0"))

### Equivalence tests ###
# This method allows for some "crud" in the residual contrast
# It probably makes sense to do these tests jointly, as above
lht(m_four, c("c2=.4", "c3=.4"))
lht(m_four, c("c2=-.4", "c3=-.4"))
# We are permitted to divide the p-values by 2 because these are one-sided tests 
lht(m_four, c("c2=.4", "c3=.4"))[2, 6]/2
lht(m_four, c("c2=-.4", "c3=-.4"))[2, 6]/2


#### FIVE GROUP ####

# Ceofficients for this case. Modify as desired
# First is the intercept, second is the focal contrast, remainder are the residual contrasts
coefs <- c(0, .7, .2, .1, -.6)

# Create the data
# The last line creates the outcome using the coefficients above and the desired sampling error, err
# %*% is matrix multiplication
dat_five <- data.frame(group = rep(paste("group", 1:5), n_per_cell))
dat_five <- mutate(dat_five,
                   c1 = case_when(group=="group 1" ~ -4/5, TRUE ~ 1/5),
                   c2 = case_when(group=="group 1" ~ 0, group=="group 2" ~ 3/4, TRUE ~ -1/4),
                   c3 = case_when(group %in% c("group 1", "group 2") ~ 0, group=="group 3" ~ 2/3, TRUE ~ -1/3),
                   c4 = case_when(group %in% c("group 1", "group 2", "group 3") ~ 0, group=="group 4" ~ 1/2, TRUE ~ -1/2),
                   outcome = as.vector(cbind(1, c1, c2, c3, c4) %*% coefs + rnorm(nrow(dat_five), sd=err)))

m_five <- lm(outcome ~ c1 + c2 + c3 + c4, data=dat_five)

### Test of the focal contrast against 0 ###
summary(m_five)

### Test of the residual contrasts against 0 ###
# It probably makes sense to do a joint test of all of them in this case
# because there are multiple residual contrasts rather than one
lht(m_five, c("c2=0", "c3=0", "c4=0"))

### Equivalence tests ###
# This method allows for some "crud" in the residual contrast
# It probably makes sense to do these tests jointly, as above
lht(m_five, c("c2=.4", "c3=.4", "c4=.4"))
lht(m_five, c("c2=-.4", "c3=-.4", "c4=-.4"))
# We are permitted to divide the p-values by 2 because these are one-sided tests 
lht(m_five, c("c2=.4", "c3=.4", "c4=.4"))[2, 6]/2
lht(m_five, c("c2=-.4", "c3=-.4", "c4=-.4"))[2, 6]/2