Skip to content

Instantly share code, notes, and snippets.

@schnarrd
Last active February 8, 2019 19:39
Show Gist options
  • Save schnarrd/f9b217a4eac09f3729d4e2011e04391d to your computer and use it in GitHub Desktop.
Save schnarrd/f9b217a4eac09f3729d4e2011e04391d to your computer and use it in GitHub Desktop.
#### 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment