Skip to content
{{ message }}

Instantly share code, notes, and snippets.

# schnarrd/severe_contrasts.R

Last active Feb 8, 2019
 #### 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
to join this conversation on GitHub. Already have an account? Sign in to comment