8 Extending the Fundamental Principles of Moderation Analysis

As Hayes opened, “in this chapter, [we’ll see] how [the] principles of moderation analysis are applied when the moderator is dichotomous (rather than a continuum, as in the previous chapter) as well as when both focal antecedent and moderator are continuous (p. 267).”

8.1 Moderation with a dichotomous moderator

Here we load a couple necessary packages, load the data, and take a glimpse().

## Observations: 211
## Variables: 5
## $ id      <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23…
## $ frame   <dbl> 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, …
## $ donate  <dbl> 5.6, 4.2, 4.2, 4.6, 3.0, 5.0, 4.8, 6.0, 4.2, 4.4, 5.8, 6.2, 6.0, 4.2, 4.4, 5.8, 5…
## $ justify <dbl> 2.95, 2.85, 3.00, 3.30, 5.00, 3.20, 2.90, 1.40, 3.25, 3.55, 1.55, 1.60, 1.65, 2.6…
## $ skeptic <dbl> 1.8, 5.2, 3.2, 1.0, 7.6, 4.2, 4.2, 1.2, 1.8, 8.8, 1.0, 5.4, 2.2, 3.6, 7.8, 1.6, 1…

Regardless of whether the antecedent variables are continuous or binary, the equation for the simple moderation is still

\[Y = i_Y + b_1 X + b_2 W + b_3 XW + e_Y.\]

We can use that equation to fit our first moderation model with a binary \(W\) (i.e., frame) like so.

Check the summary.

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: justify ~ 1 + skeptic + frame + frame:skeptic 
##    Data: disaster (Number of observations: 211) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##               Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
## Intercept        2.455     0.148    2.157    2.740 1.000     2727     3010
## skeptic          0.104     0.037    0.031    0.180 1.001     2631     2711
## frame           -0.567     0.216   -0.986   -0.151 1.001     2332     2631
## skeptic:frame    0.202     0.054    0.096    0.310 1.001     2205     2528
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
## sigma    0.819     0.041    0.744    0.906 1.001     3270     2730
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

We’ll compute our Bayeisan \(R^2\) in the typical way.

##    Estimate Est.Error  Q2.5 Q97.5
## R2    0.249     0.044 0.164 0.333

This model should look familiar to you because it is exactly the same model estimated in the analysis presented in Chapter 7 (see Table 7.4, model 3). The only differences between these two analyses are how the corresponding question is framed, meaning which variable is deemed the focal antecedent and which is the moderator, and how these variables are symbolically labeled as \(X\) and \(W\). In the analysis in Chapter 7, the focal antecedent variable was a dichotomous variable coding the framing of the cause of the disaster (labeled \(X\) then, but \(W\) now), whereas in this analysis, the focal antecedent is a continuous variable placing each person on a continuum of climate change skepticism (labeled \(W\) then, but \(X\) now), with the moderator being a dichotomous variable coding experimental condition.

So this example illustrates the symmetry property of interactions introduced in section 7.1. (p. 272)

8.1.1 Visualizing and probing the interaction.

For the plots in this chapter, we’ll take our color palette from the ochRe package, which provides Australia-inspired colors. We’ll also use a few theme settings from good-old ggthemes. As in the last chapter, we’ll save our adjusted theme settings as an object, theme_08.

Happily, the ochRe package has a handy convenience function, viz_palette(), that makes it easy to preview the colors available in a given palette. We’ll be using “olsen_qual” and “olsen_seq”.

Behold our Figure 8.3.

# these will come in handy with `geom_text()`, below
green_slope <- (fixef(model8.1)["skeptic", 1] + fixef(model8.1)[4, 1]) %>% round(digits = 3)
blue_slope  <-  fixef(model8.1)["skeptic", 1] %>% round(digits = 3)

# define the new data
nd <-
  crossing(frame   = 0:1,
           skeptic = c(0, 7))

# extract the `fitted()` draws
fitted(model8.1, 
       newdata = nd, 
       summary = F) %>% 
  data.frame() %>% 
  gather() %>% 
  # bind the new data
  bind_cols(
    nd %>% 
      expand(nesting(frame, skeptic),
             iter = 1:4000)
  ) %>% 
  
  # plot!
  ggplot(aes(x = skeptic, y = value, 
             group = interaction(frame, iter),
             color = frame %>% as.character())) +
  geom_line(aes(color = frame %>% as.character()),
            size = 1/6, alpha = 1/25) +
  geom_text(data = tibble(skeptic = c(4, 4.6),
                          value   = c(3.5, 2.6),
                          frame   = 1:0,
                          iter    = 0,
                          label   = c(paste("the slope when frame = 1 is about", green_slope),
                                      paste("the slope when frame = 0 is about", blue_slope)),
                          angle   = c(28, 5)),
            aes(label = label, angle = angle)) +
  scale_color_manual(NULL, values = ochre_palettes[["olsen_qual"]][(5:6)]) +
  scale_x_continuous(breaks = 1:6) +
  coord_cartesian(xlim = 1:6,
                  ylim = 2:4) +
  labs(title    = "The moderation effect in a spaghetti plot",
       subtitle = "Using brms defaults, we have 4000 posterior draws. Much like we have\n4000 draws for each model parameter, when we combine those draws\nacross parameters, we end up with 4000 model equations. The\nconsequence is we can also express the line plot as 4000\nsemitransparent lines, one for each level of our moderator, frame.",
       x = expression(paste("Climate Change Skepticism (", italic("X"), ")")),
       y = "Strength of Justification for Withholding Aid") +
  theme_08 +
  theme(legend.position = "none")

In addition to our fancy Australia-inspired colors, we’ll also play around a bit with spaghetti plots in this chapter. To my knowledge, this use of spaghetti plots is uniquely Bayesian. If you’re trying to wrap your head around what on earth we just did, take a look at the first few rows from posterior_samples() object, post.

##   b_Intercept  b_skeptic    b_frame b_skeptic:frame     sigma      lp__
## 1    2.171268 0.18122884 -0.6213774       0.1814370 0.8144831 -264.4922
## 2    2.112918 0.19286032 -0.5346854       0.1749145 0.8169612 -264.5730
## 3    2.580092 0.04726581 -0.3385162       0.1657682 0.8691147 -264.6696
## 4    2.609206 0.04635323 -0.3941262       0.1594145 0.8658938 -264.7612
## 5    2.408386 0.12717550 -0.5910781       0.1956912 0.8121655 -260.1720
## 6    2.383672 0.11459057 -0.5955617       0.2302462 0.8024598 -260.5167

The head() function returned six rows, each one corresponding to the credible parameter values from a given posterior draw. The lp__ is uniquely Bayesian and beyond the scope of this project. You might think of sigma as the Bayesian analogue to what the OLS folks often refer to as error or the residual variance. Hayes doesn’t tend to emphasize it in this text, but it’s something you’ll want to pay increasing attention to as you move along in your Bayesian career. All the columns starting with b_ are the regression parameters, the model coefficients or the fixed effects. But anyways, notice that those b_ columns correspond to the four parameter values in Formula 8.2 on page 270. Here they are, but reformatted to more closely mimic the text:

  1. \(\hat{Y}\) = 2.171 + 0.181\(X\) + -0.621\(W\) + 0.181XW
  2. \(\hat{Y}\) = 2.113 + 0.193\(X\) + -0.535\(W\) + 0.175XW
  3. \(\hat{Y}\) = 2.58 + 0.047\(X\) + -0.339\(W\) + 0.166XW
  4. \(\hat{Y}\) = 2.609 + 0.046\(X\) + -0.394\(W\) + 0.159XW
  5. \(\hat{Y}\) = 2.408 + 0.127\(X\) + -0.591\(W\) + 0.196XW
  6. \(\hat{Y}\) = 2.384 + 0.115\(X\) + -0.596\(W\) + 0.23XW

Each row of post, each iteration or posterior draw, yields a full model equation that is a credible description of the data–or at least as credible as we can get within the limits of the model we have specified, our priors (which we typically cop out on and just use defaults in this project), and how well those fit when applied to the data at hand. So when we use brms convenience functions like fitted(), we pass specific predictor values through those 4000 unique model equations, which returns 4000 similar but distinct expected \(Y\)-values. So although a nice way to summarize those 4000 values is with summaries such as the posterior mean/median and 95% intervals, another way is to just plot an individual regression line for each of the iterations. That is what’s going on when we depict out models with a spaghetti plot.

The thing I like about spaghetti plots is that they give a three-dimensional sense of the posterior. Note that each individual line is very skinny and semitransparent. When you pile a whole bunch of them atop each other, the peaked or most credible regions of the posterior are the most saturated in color. Less credible posterior regions almost seamlessly merge into the background. Also, note how the combination of many similar but distinct straight lines results in a bowtie shape. Hopefully this clarifies where that shape’s been coming from when we use geom_ribbon() to plot the 95% intervals.

Back to the text, on the bottom of page 274, Hayes pointed out the conditional effect of skeptic when frame == 1 is \(b_1 + b_3 = 0.306\). We can show that with a little arithmetic followed up with tidybayes::mean_qi().

## # A tibble: 1 x 6
##   `b1 + b3` .lower .upper .width .point .interval
##       <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1     0.306   0.23  0.384   0.95 mean   qi

But anyways, you could recode frame in a number of ways, including if_else() or, in this case, by simple arithmetic.

With frame_ep in hand, we’re ready to refit the model.

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: justify ~ skeptic + frame_ep + skeptic:frame_ep 
##    Data: disaster (Number of observations: 211) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept            1.88      0.16     1.56     2.19 1.00     2193     2986
## skeptic              0.31      0.04     0.23     0.39 1.00     2194     2480
## frame_ep             0.57      0.21     0.14     0.98 1.00     1891     2273
## skeptic:frame_ep    -0.20      0.05    -0.31    -0.10 1.00     1907     2081
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.82      0.04     0.74     0.90 1.00     3017     2732
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Our results match nicely with the formula on page 275.

If you want to follow along with Hayes on page 276 and isolate the 95% credible intervals for the skeptic parameter, you can use the posterior_interval() function.

##  2.5% 97.5% 
## 0.231 0.389

8.2 Interaction between two quantitative variables

Here’s the glbwarm data.

## Observations: 815
## Variables: 7
## $ govact   <dbl> 3.6, 5.0, 6.6, 1.0, 4.0, 7.0, 6.8, 5.6, 6.0, 2.6, 1.4, 5.6, 7.0, 3.8, 3.4, 4.2, …
## $ posemot  <dbl> 3.67, 2.00, 2.33, 5.00, 2.33, 1.00, 2.33, 4.00, 5.00, 5.00, 1.00, 4.00, 1.00, 5.…
## $ negemot  <dbl> 4.67, 2.33, 3.67, 5.00, 1.67, 6.00, 4.00, 5.33, 6.00, 2.00, 1.00, 4.00, 5.00, 4.…
## $ ideology <dbl> 6, 2, 1, 1, 4, 3, 4, 5, 4, 7, 6, 4, 2, 4, 5, 2, 6, 4, 2, 4, 4, 2, 6, 4, 4, 3, 4,…
## $ age      <dbl> 61, 55, 85, 59, 22, 34, 47, 65, 50, 60, 71, 60, 71, 59, 32, 36, 69, 70, 41, 48, …
## $ sex      <dbl> 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1,…
## $ partyid  <dbl> 2, 1, 1, 1, 1, 2, 1, 1, 2, 3, 2, 1, 1, 1, 1, 1, 2, 3, 1, 3, 2, 1, 3, 2, 1, 1, 1,…

In this section we add three covariates (i.e., \(C\) variables) to the basic moderation model. Although Hayes made a distinction between the \(X\), \(M\), and \(C\) variables in the text, that distinction is conceptual and doesn’t impact the way we enter them into brm(). Rather, the brm() formula clarifies they’re all just predictors.

Our results cohere nicely with the Hayes’s formula in the middle of page 278 or with the results he displayed in Table 8.2.

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: govact ~ 1 + negemot + age + negemot:age + posemot + ideology + sex 
##    Data: glbwarm (Number of observations: 815) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##             Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
## Intercept      5.176     0.340    4.528    5.838 1.000     2494     2796
## negemot        0.119     0.083   -0.042    0.282 1.002     2320     2502
## age           -0.024     0.006   -0.036   -0.012 1.001     2279     2394
## posemot       -0.021     0.028   -0.078    0.034 1.002     4074     2772
## ideology      -0.211     0.027   -0.264   -0.159 1.000     4268     2995
## sex           -0.012     0.075   -0.163    0.134 1.002     4037     2925
## negemot:age    0.006     0.002    0.003    0.009 1.001     2364     2564
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
## sigma    1.059     0.026    1.009    1.112 1.000     3783     2748
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Here’s the \(R^2\) summary.

##    Estimate Est.Error Q2.5 Q97.5
## R2    0.401      0.02 0.36  0.44

As the \(R^2\) is a good bit away from the boundaries, it’s nicely Gaussian.

8.2.1 Visualizing and probing the interaction.

For our version of Figure 8.5, we’ll need to adjust our nd data for fitted().

## # A tibble: 6 x 5
##   negemot   age posemot ideology   sex
##     <dbl> <dbl>   <dbl>    <dbl> <dbl>
## 1       1    30    3.13     4.08 0.488
## 2       1    50    3.13     4.08 0.488
## 3       1    70    3.13     4.08 0.488
## 4       6    30    3.13     4.08 0.488
## 5       6    50    3.13     4.08 0.488
## 6       6    70    3.13     4.08 0.488

Our fitted() and ggplot2 code will be quite similar to the last spaghetti plot. Only this time we’ll use filter() to reduce the number of posterior draws we show in the plot.

When we reduce the number of lines depicted in the plot, we lose some of the three-dimensional illusion. It’s nice, however, to get a closer look to each individual line. To each their own.

We’ll continue with our spaghetti plot approach for Figure 8.7. Again, when we made the JN technique plot for Chapter 7, we computed values for the posterior mean and the 95% intervals. Because the intervals follow a bowtie shape, we had to compute the \(Y\)-values for many values across the x-axis in order to make the curve look smooth. But as long as we stick with the spaghetti plot approach, all we need are the values at the endpoints of each iteration. Although each line is straight, the combination of many lines is what produces the bowtie effect.

In other words, each of those orange lines is a credible expression of \(\theta_{X \rightarrow Y}\) (i.e., \(b_1 + b_3 W\)) across a continuous range of \(W\) values.

8.3 Hierarchical versus simultaneous entry

Many investigators test a moderation hypothesis in regression analysis using a method that on the surface seems different than the procedure described thus far. This alternative approach is to build a regression model by adding the product of \(X\) and \(W\) to a model already containing \(X\) and \(W\). This procedure is sometimes called hierarchical regression or hierarchical variable entry (and easily confused by name with hierarchical linear modeling, which is an entirely different thing). The goal using this method is to determine whether allowing \(X\)’s effect to be contingent on \(W\) produces a better fitting model than one in which the effect of \(X\) is constrained to be unconditional on \(W\). According to the logic of hierarchical entry, if the contingent model accounts for more of the variation in Y than the model that forces \(X\)’s effect to be independent of \(W\), then the better model is one in which \(W\) is allowed to moderate \(X\)’s effect. Although this approach works, it is a widely believed myth that it is necessary to use this approach in order to test a moderation hypothesis. (p. 289, emphasis in the original)

Although this method is not necessary, it can be handy to slowly build your model. This method can also serve nice rhetorical purposes in a paper. Anyway, here’s our multivariable but non-moderation model, model8.4.

Here we’ll compute the corresponding \(R^2\) and compare it with the one for the original interaction model with a difference score.

## # A tibble: 3 x 7
##   R2                  value .lower .upper .width .point .interval
##   <chr>               <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 difference           0.05 -0.073  0.169   0.95 median qi       
## 2 moderation_model     0.25  0.164  0.333   0.95 median qi       
## 3 multivariable_model  0.2   0.113  0.284   0.95 median qi

Note that the Bayesian \(R^2\) performed differently than the \(F\)-test in the text.

We can also compare these with the LOO, which, as is typical of information criteria, corrects for model complexity. First, we compute them and attach the results to the model fit objects.

Now use the loo_compare() function to compare them directly.

##          elpd_diff se_diff elpd_loo se_elpd_loo p_loo  se_p_loo looic  se_looic
## model8.1    0.0       0.0  -259.1     10.8         5.5    0.9    518.2   21.6  
## model8.4   -5.8       4.1  -264.8     11.2         4.9    0.9    529.7   22.4

As a reminder, we generally prefer models with lower information criteria, which in this case is clearly the moderation model (i.e., model8.1). However, the standard error value (i.e., se_diff) for the difference (i.e., elpd_diff) is quite large, which suggests that the model with the lowest value isn’t the clear winner. Happily, these results match nicely with the Bayesian \(R^2\) difference score. The moderation model appears somewhat better than the multivariable model, but its superiority is hardly decisive.

8.4 The equivalence between moderated regression analysis and a 2 X 2 factorial analysis of variance

I’m just not going to encourage ANOVA \(F\)-testing methodology. However, I will show the Bayesian regression model. First, here are the data.

## Observations: 541
## Variables: 7
## $ policy   <dbl> 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0,…
## $ interest <dbl> 4.0, 2.0, 3.0, 1.0, 1.0, 2.0, 1.0, 2.5, 3.0, 1.0, 2.0, 3.5, 1.0, 1.0, 1.5, 3.0, …
## $ age      <dbl> 39, 57, 63, 56, 50, 87, 33, 64, 82, 28, 18, 52, 42, 39, 64, 72, 54, 84, 55, 27, …
## $ educ     <dbl> 3, 3, 2, 5, 3, 2, 7, 2, 3, 3, 1, 1, 5, 4, 3, 2, 3, 4, 7, 2, 3, 5, 4, 5, 5, 3, 3,…
## $ male     <dbl> 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1,…
## $ conserv  <dbl> 4, 3, 6, 3, 3, 5, 6, 3, 6, 7, 4, 2, 7, 6, 5, 6, 6, 3, 7, 6, 5, 5, 3, 4, 6, 2, 7,…
## $ kerry    <dbl> 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0,…

Fit the moderation model.

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: interest ~ 1 + policy + kerry + policy:kerry 
##    Data: caskets (Number of observations: 541) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept        1.78      0.09     1.61     1.96 1.00     2374     2968
## policy          -0.38      0.12    -0.63    -0.14 1.00     2096     2639
## kerry            0.60      0.13     0.34     0.84 1.00     2176     2434
## policy:kerry     0.36      0.17     0.03     0.71 1.00     1926     2375
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.04      0.03     0.98     1.10 1.00     3741     2853
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Those results don’t look anything like what Hayes reported in Tables 8.3 or 8.4. However, a little deft manipulation of the posterior samples can yield equivalent results to Hayes’s Table 8.3.

Here are the cell-specific means in Table 8.3.

## # A tibble: 4 x 7
##   name    value .lower .upper .width .point .interval
##   <chr>   <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 Y_bar_1  1.78   1.61   1.96   0.95 median qi       
## 2 Y_bar_2  1.40   1.23   1.58   0.95 median qi       
## 3 Y_bar_3  2.38   2.21   2.56   0.95 median qi       
## 4 Y_bar_4  2.36   2.19   2.53   0.95 median qi

And here are the marginal means from Table 8.3.

## # A tibble: 4 x 7
##   name     value .lower .upper .width .point .interval
##   <chr>    <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 Y_bar_12  1.59   1.47   1.72   0.95 median qi       
## 2 Y_bar_13  2.08   1.96   2.21   0.95 median qi       
## 3 Y_bar_24  1.88   1.76   2.00   0.95 median qi       
## 4 Y_bar_34  2.37   2.24   2.50   0.95 median qi

For kicks and giggles, here are what the cell-specific means look like in box plots.

And here are the same for the marginal means. This time we’ll show the shapes of the posteriors with violin plots with horizontal lines depicting the median and interquartile ranges.

On page 294, Hayes used point estimates to compute the simple effect of policy information among Kerry supporters and then the same thing among Bush supporters. Here’s how we’d do that when working with the full vector of posterior draws.

## # A tibble: 2 x 7
##   name                 value .lower .upper .width .point .interval
##   <chr>                <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 simple_effect_Bush  -0.382 -0.631 -0.142   0.95 mean   qi       
## 2 simple_effect_Kerry -0.023 -0.265  0.221   0.95 mean   qi

So then computing the main effect for policy information using the simple effects is little more than an extension of those steps.

## # A tibble: 1 x 6
##   main_effect .lower .upper .width .point .interval
##         <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1      -0.203 -0.381 -0.029   0.95 mean   qi

And we get the same results by strategically subtracting the marginal means.

## # A tibble: 1 x 6
##   main_effect .lower .upper .width .point .interval
##         <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1      -0.203 -0.381 -0.029   0.95 mean   qi

The main effect of for candidate is similarly computed using either approach.

## # A tibble: 1 x 6
##   main_effect .lower .upper .width .point .interval
##         <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1        0.78  0.606  0.958   0.95 mean   qi
## # A tibble: 1 x 6
##   main_effect .lower .upper .width .point .interval
##         <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1        0.78  0.606  0.958   0.95 mean   qi

We don’t have an \(F\)-test for our Bayesian moderation model. But we do have an interaction term. Here’s its distribution.

Following Hayes’s work on the bottom of page 295, here’s how you’d reproduce that by manipulating our \(\overline Y\) vectors.

## # A tibble: 1 x 6
##   reproduced_interaction_term .lower .upper .width .point .interval
##                         <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1                        0.36   0.03   0.71   0.95 mean   qi

Extending that logic, we also get the answer this way.

## # A tibble: 1 x 6
##   reproduced_interaction_term .lower .upper .width .point .interval
##                         <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1                        0.36   0.03   0.71   0.95 mean   qi

8.4.1 Simple effects parameterization.

We might reacquaint ourselves with the formula from model8.5.

## interest ~ 1 + policy + kerry + policy:kerry

The results cohere nicely with the “Model 1” results at the top of Table 8.5.

##              Estimate Est.Error   Q2.5  Q97.5
## Intercept       1.781     0.089  1.611  1.961
## policy         -0.382     0.124 -0.631 -0.142
## kerry           0.601     0.127  0.345  0.845
## policy:kerry    0.359     0.174  0.026  0.712

The Bayesian \(R^2\) portion looks on point, too.

##    Estimate Est.Error  Q2.5 Q97.5
## R2     0.14     0.025 0.092  0.19

Our various Y_bar transformations from before continue to cohere with the coefficients, above, just like in the text. E.g., the policy coefficient may be returned like so.

## # A tibble: 2 x 3
##   name                mean    sd
##   <chr>              <dbl> <dbl>
## 1 b1                -0.382 0.124
## 2 Y_bar_2 - Y_bar_1 -0.382 0.124

We can continue to use Hayes’s Y_bar transformations to return the kerry coefficient, too.

## # A tibble: 2 x 3
##   name               mean    sd
##   <chr>             <dbl> <dbl>
## 1 b2                0.601 0.127
## 2 Y_bar_3 - Y_bar_1 0.601 0.127

Here we compute \(b_3\) with the difference between the simple effects of \(X\) at levels of \(W\).

## # A tibble: 2 x 3
##   name                                       mean    sd
##   <chr>                                     <dbl> <dbl>
## 1 (Y_bar_4 - Y_bar_3) - (Y_bar_2 - Y_bar_1) 0.359 0.174
## 2 b3                                        0.359 0.174

And now \(b_{3}\) with the difference between the simple effects of \(W\) at levels of \(X\).

## # A tibble: 2 x 3
##   name                                       mean    sd
##   <chr>                                     <dbl> <dbl>
## 1 (Y_bar_4 - Y_bar_2) - (Y_bar_3 - Y_bar_1) 0.359 0.174
## 2 b3                                        0.359 0.174

8.4.2 Main effects parameterization.

A nice feature of brms is you can transform your data right within the brm() or update() functions. Here we’ll make our two new main-effects-coded variables, policy_me and kerry_me, with the mutate() function right within update().

Transforming your data within the brms functions won’t change the original data structure. However, brms will save the data used to fit the model within the brm() object. You can access that data like so.

##   interest policy_me kerry_me
## 1        4       0.5      0.5
## 2        2      -0.5      0.5
## 3        3       0.5      0.5
## 4        1       0.5      0.5
## 5        1       0.5      0.5
## 6        2      -0.5     -0.5

But we digress. Here’s our analogue to the “Model 2” portion of Table 8.5.

##                    Estimate Est.Error   Q2.5  Q97.5
## Intercept             1.981     0.045  1.893  2.070
## policy_me            -0.208     0.090 -0.387 -0.032
## kerry_me              0.780     0.090  0.607  0.955
## policy_me:kerry_me    0.354     0.180  0.007  0.707
## R2                    0.140     0.025  0.093  0.191

Like with model8.6, above, we’ll need a bit of algebra to compute our \(\overline Y_i\) vectors.

With our post for fit5 in hand, we’ll follow the formulas at the top of page 298 to compute our \(b_1\) and \(b_2\) distributions.

## # A tibble: 2 x 3
##   name    mean    sd
##   <chr>  <dbl> <dbl>
## 1 b1    -0.208  0.09
## 2 b2     0.78   0.09

Hayes pointed out that the interaction effect, \(b_3\), is the same across models his OLS Models 1 and 2. This is largely true for our Bayesian HMC model8.5 and model8.6 models.

##  Estimate Est.Error      Q2.5     Q97.5 
##     0.359     0.174     0.026     0.712
##  Estimate Est.Error      Q2.5     Q97.5 
##     0.354     0.180     0.007     0.707

However, the results aren’t exactly the same because of simulation error. If you were working on a project requiring high precision, increase the number of posterior iterations. To demonstrate, here we’ll increase each chain’s post-warmup iteration count by an order of magnitude, resulting in 80,000 post-warmup iterations rather than the default 4,000.

Now they’re quite a bit closer.

##  Estimate Est.Error      Q2.5     Q97.5 
##     0.360     0.178     0.011     0.711
##  Estimate Est.Error      Q2.5     Q97.5 
##     0.360     0.178     0.011     0.709

And before you get fixate on how there are still differences after 80,000 iterations, each, consider comparing the two density plots.

8.4.3 Conducting a \(2 \times 2\) between-participants factorial ANOVA using PROCESS another regression model with brms.

Since we’re square in single-level regression land with our brms approach, there’s no direct analogue for us, here. However, notice the post-ANOVA \(t\)-tests Hayes presented on page 300. If we just want to consider the \(2 \times 2\) structure of our two dummy variables as indicative of four groups, we have one more coding system available for the job. With the handy str_c() function, we’ll concatenate the policy and kerry values into a nominal variable, policy_kerry. Here’s what that looks like:

## # A tibble: 6 x 8
##   policy interest   age  educ  male conserv kerry policy_kerry
##    <dbl>    <dbl> <dbl> <dbl> <dbl>   <dbl> <dbl> <chr>       
## 1      1        4    39     3     1       4     1 11          
## 2      0        2    57     3     1       3     1 01          
## 3      1        3    63     2     0       6     1 11          
## 4      1        1    56     5     1       3     1 11          
## 5      1        1    50     3     0       3     1 11          
## 6      0        2    87     2     1       5     0 00

Now check out what happens if we reformat our formula to interest ~ 0 + policy_kerry.

The brm() function recnognized policy_kerry was a character vector and treated it as a nominal variable. The 0 + part of the function removed the model intercept. Here’s how that effects the output.

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: interest ~ 0 + policy_kerry 
##    Data: caskets (Number of observations: 541) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## policy_kerry00     1.78      0.09     1.61     1.96 1.00     4728     3033
## policy_kerry01     2.38      0.09     2.21     2.56 1.00     4673     3337
## policy_kerry10     1.40      0.09     1.22     1.57 1.00     5379     3114
## policy_kerry11     2.36      0.09     2.19     2.52 1.00     4375     2931
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.04      0.03     0.98     1.10 1.00     4792     2939
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Without the typical intercept, brm() estimated the means for each of the four policy_kerry groups. It’s kinda like an intercept-only model, but with four intercepts. Here’s what their densities look like:

Since each of the four primary vectors in our post object is of a group mean, it’s trivial to compute difference scores. To compute the difference score analogous to Hayes’s two \(t\)-tests, we’d do the following.

## # A tibble: 2 x 7
##   name          value .lower .upper .width .point .interval
##   <chr>         <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 difference_1 -0.386 -0.639 -0.137   0.95 median qi       
## 2 difference_2 -0.026 -0.273  0.213   0.95 median qi

Session info

## R version 3.6.0 (2019-04-26)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] tidybayes_1.1.0 ochRe_1.0.0     ggthemes_4.2.0  brms_2.10.3     Rcpp_1.0.2      forcats_0.4.0  
##  [7] stringr_1.4.0   dplyr_0.8.3     purrr_0.3.3     readr_1.3.1     tidyr_1.0.0     tibble_2.1.3   
## [13] ggplot2_3.2.1   tidyverse_1.2.1
## 
## loaded via a namespace (and not attached):
##   [1] colorspace_1.4-1          ellipsis_0.3.0            ggridges_0.5.1           
##   [4] rsconnect_0.8.15          ggstance_0.3.2            markdown_1.1             
##   [7] base64enc_0.1-3           rstudioapi_0.10           rstan_2.19.2             
##  [10] svUnit_0.7-12             DT_0.9                    fansi_0.4.0              
##  [13] lubridate_1.7.4           xml2_1.2.0                codetools_0.2-16         
##  [16] bridgesampling_0.7-2      knitr_1.23                shinythemes_1.1.2        
##  [19] zeallot_0.1.0             bayesplot_1.7.0           jsonlite_1.6             
##  [22] broom_0.5.2               shiny_1.3.2               compiler_3.6.0           
##  [25] httr_1.4.0                backports_1.1.5           assertthat_0.2.1         
##  [28] Matrix_1.2-17             lazyeval_0.2.2            cli_1.1.0                
##  [31] later_1.0.0               htmltools_0.4.0           prettyunits_1.0.2        
##  [34] tools_3.6.0               igraph_1.2.4.1            coda_0.19-3              
##  [37] gtable_0.3.0              glue_1.3.1.9000           reshape2_1.4.3           
##  [40] raster_2.9-5              cellranger_1.1.0          vctrs_0.2.0              
##  [43] nlme_3.1-139              crosstalk_1.0.0           xfun_0.10                
##  [46] ps_1.3.0                  rvest_0.3.4               mime_0.7                 
##  [49] miniUI_0.1.1.1            lifecycle_0.1.0           gtools_3.8.1             
##  [52] zoo_1.8-6                 scales_1.0.0              colourpicker_1.0         
##  [55] hms_0.4.2                 promises_1.1.0            Brobdingnag_1.2-6        
##  [58] parallel_3.6.0            inline_0.3.15             shinystan_2.5.0          
##  [61] gridExtra_2.3             loo_2.1.0                 StanHeaders_2.19.0       
##  [64] stringi_1.4.3             dygraphs_1.1.1.6          pkgbuild_1.0.5           
##  [67] rlang_0.4.1               pkgconfig_2.0.3           matrixStats_0.55.0       
##  [70] evaluate_0.14             lattice_0.20-38           rstantools_2.0.0         
##  [73] htmlwidgets_1.5           labeling_0.3              tidyselect_0.2.5         
##  [76] processx_3.4.1            plyr_1.8.4                magrittr_1.5             
##  [79] R6_2.4.0                  generics_0.0.2            pillar_1.4.2             
##  [82] haven_2.1.0               withr_2.1.2               xts_0.11-2               
##  [85] sp_1.3-1                  abind_1.4-5               modelr_0.1.4             
##  [88] crayon_1.3.4              arrayhelpers_1.0-20160527 utf8_1.1.4               
##  [91] rmarkdown_1.13            grid_3.6.0                readxl_1.3.1             
##  [94] callr_3.3.2               threejs_0.3.1             digest_0.6.21            
##  [97] xtable_1.8-4              httpuv_1.5.2              stats4_3.6.0             
## [100] munsell_0.5.0             shinyjs_1.0