6.3 Using a different group coding system
Here we’ll make our alternative dummies, what we’ll call D_1 and D_2, with orthogonal contrast coding.
protest <-
protest %>%
mutate(D_1 = ifelse(protest == 0, -2/3, 1/3),
D_2 = ifelse(protest == 0, 0,
ifelse(protest == 1, -1/2, 1/2)))Here are the sub-models.
m_model <- bf(respappr ~ 1 + D_1 + D_2)
y_model <- bf(liking ~ 1 + D_1 + D_2 + respappr)Now we fit.
model4 <-
brm(data = protest, family = gaussian,
m_model + y_model + set_rescor(FALSE),
chains = 4, cores = 4)Here are our intercepts and regression coefficient summaries.
fixef(model4)## Estimate Est.Error Q2.5 Q97.5
## respappr_Intercept 4.8399962 0.10578752 4.6301115 5.0460324
## liking_Intercept 3.6303107 0.35919791 2.9137984 4.3166260
## respappr_D_1 1.4358633 0.21557632 1.0162848 1.8717145
## respappr_D_2 0.3497242 0.25111687 -0.1406033 0.8324690
## liking_D_1 -0.1146844 0.20498317 -0.5065575 0.2993028
## liking_D_2 -0.2168985 0.20270585 -0.6162726 0.1886207
## liking_respappr 0.4129279 0.07201573 0.2767746 0.5557400It’s important to note that these will not correspond to the “TOTAL EFFECT MODEL” section of the PROCESS output of Figure 6.3. Hayes’s PROCESS has the mcx=3 command which tells the program to reparametrize the orthogonal contrasts. brms doesn’t have such a command.
For now, we’ll have to jump to equation 6.8 towards the bottom of page 207. Those parameters are evident in our output.
fixef(model4)[c(1, 3:4) , ] %>% round(digits = 3)## Estimate Est.Error Q2.5 Q97.5
## respappr_Intercept 4.840 0.106 4.630 5.046
## respappr_D_1 1.436 0.216 1.016 1.872
## respappr_D_2 0.350 0.251 -0.141 0.832Thus it’s easy to get the \(\overline{M}_{ij}\) means with a little posterior manipulation.
post <- posterior_samples(model4)
post <-
post %>%
mutate(M_np = b_respappr_Intercept + b_respappr_D_1*-2/3 + b_respappr_D_2*0,
M_ip = b_respappr_Intercept + b_respappr_D_1*1/3 + b_respappr_D_2*-1/2,
M_cp = b_respappr_Intercept + b_respappr_D_1*1/3 + b_respappr_D_2*1/2)
post %>%
select(starts_with("M_")) %>%
gather() %>%
mutate(key = factor(key, levels = c("M_np", "M_ip", "M_cp"))) %>%
group_by(key) %>%
summarize(mean = mean(value),
sd = sd(value))## # A tibble: 3 x 3
## key mean sd
## <fct> <dbl> <dbl>
## 1 M_np 3.88 0.182
## 2 M_ip 5.14 0.181
## 3 M_cp 5.49 0.173With these in hand, we can compute \(a_{1}\) and \(a_{2}\).
post <-
post %>%
mutate(a1 = (M_ip + M_cp)/2 - M_np,
a2 = M_cp - M_ip)
post %>%
select(a1:a2) %>%
gather() %>%
group_by(key) %>%
summarize(mean = mean(value),
sd = sd(value))## # A tibble: 2 x 3
## key mean sd
## <chr> <dbl> <dbl>
## 1 a1 1.44 0.216
## 2 a2 0.350 0.251Happily, our model output will allow us to work with Hayes’s \(\overline{Y}^*_{ij}\) equations in the middle of page 210.
post <-
post %>%
mutate(Y_np = b_liking_Intercept + b_liking_D_1*-2/3 + b_liking_D_2*0 + b_liking_respappr*mean(protest$respappr),
Y_ip = b_liking_Intercept + b_liking_D_1*1/3 + b_liking_D_2*-1/2 + b_liking_respappr*mean(protest$respappr),
Y_cp = b_liking_Intercept + b_liking_D_1*1/3 + b_liking_D_2*1/2 + b_liking_respappr*mean(protest$respappr))
post %>%
select(starts_with("Y_")) %>%
gather() %>%
mutate(key = factor(key, levels = c("Y_np", "Y_ip", "Y_cp"))) %>%
group_by(key) %>%
summarize(mean = mean(value),
sd = sd(value))## # A tibble: 3 x 3
## key mean sd
## <fct> <dbl> <dbl>
## 1 Y_np 5.72 0.160
## 2 Y_ip 5.71 0.144
## 3 Y_cp 5.49 0.149And with these in hand, we can compute \(c'_{1}\) and \(c'_{2}\).
post <-
post %>%
mutate(c1_prime = (Y_ip + Y_cp)/2 - Y_np,
c2_prime = Y_cp - Y_ip)
post %>%
select(c1_prime:c2_prime) %>%
gather() %>%
group_by(key) %>%
summarize(mean = mean(value),
sd = sd(value))## # A tibble: 2 x 3
## key mean sd
## <chr> <dbl> <dbl>
## 1 c1_prime -0.115 0.205
## 2 c2_prime -0.217 0.203It appears Hayes has a typo in the formula for \(c'_{2}\) on page 211. The value he has down for \(\overline{Y}^*_{IP}\), 5.145, is incorrect. It’s not the one he displayed at the bottom of the previous page and it also contradicts the analyses herein. So it goes… These things happen.
We haven’t spelled out, but the \(b\) parameter is currently labeled b_liking_respappr in our post object. Here we’ll make a b column to make things easier. While we’re at it, we’ll compute the indirect effects, too.
post <-
post %>%
mutate(b = b_liking_respappr) %>%
mutate(a1b = a1*b,
a2b = a2*b)
post %>%
select(a1b:a2b) %>%
gather() %>%
group_by(key) %>%
summarize(mean = mean(value),
sd = sd(value),
ll = quantile(value, probs = .025),
ul = quantile(value, probs = .975)) %>%
mutate_if(is.double, round, digits = 3)## # A tibble: 2 x 5
## key mean sd ll ul
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 a1b 0.593 0.137 0.345 0.882
## 2 a2b 0.144 0.108 -0.061 0.366Now we can compute and summarize() our \(c_{1}\) and \(c_{1}\).
post <-
post %>%
mutate(c1 = c1_prime + a1b,
c2 = c2_prime + a2b)
post %>%
select(c1:c2) %>%
gather() %>%
group_by(key) %>%
summarize(mean = mean(value),
sd = sd(value))## # A tibble: 2 x 3
## key mean sd
## <chr> <dbl> <dbl>
## 1 c1 0.478 0.199
## 2 c2 -0.0726 0.225