9Some Myths and Additional Extensions of Moderation Analysis
In this chapter, our Bayesian approach complicates some of Hayes’s flow. We’ll see how Bayesian HMC estimation can make us reconsider the value in mean centering and we’ll also slip in some missing data talk.
9.1 Truths and myths about mean-centering
Mean-centering has been recommended in a few highly regarded books on regression analysis (e.g., Aiken & West, 1991; Cohen et al., 2003), and several explanations have been offered for why mean-centering should be undertaken prior to computation of the product and model estimation. The explanation that seems to have resulted in the most misunderstanding is that \(X\) and \(W\) are likely to be highly correlated with \(XW\) and this will produce estimation problems caused by collinearity and result in poor or “strange” estimates of regression coefficients, large standard errors, and reduced power of the statistical test of the interaction. But this is, in large part, simply a myth. As I describe later, there are some reasons that mean-centering the focal antecedent or moderator variables can be a beneficial thing to do, which is why it has been recommended by some. However, it is incorrect to claim that it is necessary, that a failure to do so will lead one to incorrect inferences about moderation, or that the resulting regression coefficients are somehow strange or inherently uninterpretable. (Hayes, 2018, pp. 304–305, emphasis in the original)
Here we load a couple necessary packages, load the data, and take a glimpse().
Before we fit our models, we’ll go ahead and make our mean-centered predictors, negemot_c and age_c.
glbwarm <- glbwarm %>%mutate(negemot_c = negemot -mean(negemot),age_c = age -mean(age))
Now we’re ready to fit Models 1 and 2. But before we do, it’s worth repeating part of the text:
Let’s load brms.
library(brms)
As we’ll see in just a bit, there are some important reasons for Bayesians using HMC to mean center that wouldn’t pop up within the OLS paradigm. First let’s fit model9.1 and model9.2. model9.1 follows the conventional moderation equation
\[Y = i_Y + b_1 X + b_2 W + b_3 XW + e_Y.\]
model9.2 is our mean-centered model, which we can express formally as
\[\begin{align*}
Y & = i_Y + b_1 \left (X - \overline X \right ) + b_2 \left (W - \overline W \right ) + b_3 \left (X - \overline X \right ) \left (W - \overline W \right ) + e_Y, \; \text{or more simply} \\
& = i_Y + b_1 X' + b_2 W' + b_3 X'W' + e_Y,
\end{align*}\]
where \(X' = \left (X - \overline X \right)\) and so on.
model9.1<-brm(data = glbwarm, family = gaussian, govact ~1+ negemot + age + negemot:age,cores =4,file ="fits/model09.01")model9.2<-brm(data = glbwarm, family = gaussian, govact ~1+ negemot_c + age_c + negemot_c:age_c,cores =4,file ="fits/model09.02")
As with Hayes’s OLS models, our HMC models yield the same Bayesian \(R^2\) distributions, within simulation error.
Our model summaries also correspond nicely with those at the top of Table 9.1.
print(model9.1, digits =3)
Family: gaussian
Links: mu = identity
Formula: govact ~ 1 + negemot + age + negemot:age
Data: glbwarm (Number of observations: 815)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.343 0.336 3.690 4.989 1.003 1504 2102
negemot 0.145 0.086 -0.021 0.315 1.003 1429 1945
age -0.031 0.006 -0.043 -0.019 1.002 1513 1982
negemot:age 0.007 0.002 0.004 0.010 1.003 1466 2004
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 1.097 0.026 1.046 1.149 1.003 2464 2419
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
print(model9.2, digits =3)
Family: gaussian
Links: mu = identity
Formula: govact ~ 1 + negemot_c + age_c + negemot_c:age_c
Data: glbwarm (Number of observations: 815)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.598 0.039 4.523 4.673 1.001 3383 2605
negemot_c 0.500 0.026 0.448 0.550 1.003 3153 2405
age_c -0.005 0.002 -0.010 -0.001 1.000 5299 2965
negemot_c:age_c 0.007 0.002 0.004 0.010 1.000 5879 3421
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 1.098 0.027 1.047 1.152 1.002 3581 2970
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
However, notice the ‘Bulk_ESS’ and ‘Tail_ESS’ columns. The values for model9.2 were substantially larger than those for model9.1. These columns denote the number of effective samples (a.k.a. the effective sample size). Versions of brms before 2.10.0 returned a single effective sample size (ESS) value per parameter. We will cover where the change came from in a bit. But first, recall that we’ve been using brms defaults, which results in 4 HMC chains, each of which contains 2,000 draws (iterations), the first 1,000 of which are warmup values. After we discard the warmup values, that leaves 1,000 draws from each chain–4,000 total. As it turns out, Markov chains, and thus HMC chains, are typically autocorrelated, which means that each draw is partially dependent on the previous draw. Ideally, the autocorrelations are near zero. That’s often not the case.
The bayesplot package offers a variety of diagnostic plots. Here we’ll use the mcmc_acf() function to make autocorrelation plots for all model parameters.
But before we get there, we’ll be using an ggthemes-based theme for our plots in this chapter.
# loadlibrary(bayesplot)library(ggthemes)# adjust the global settingstheme_set(theme_tufte(base_size =13) +theme(panel.background =element_rect(fill ="grey98", linewidth =0)))# to reset the default ggplot2 theme to its default parameters, execute:# theme_set(theme_gray())
The theme_tufte() function will let us follow some of Edward Tufte’s preferred aesthecic principles (Tufte, 2001), and the theme_set() function will make these the default plot settings. To undo this, just execute theme_set(theme_gray()). But anyways, here are our mcmc_acf() plots.
draws1 <-as_draws_df(model9.1)draws2 <-as_draws_df(model9.2)mcmc_acf(draws1, pars =c("b_Intercept", "b_negemot", "b_age", "b_negemot:age", "sigma"),lags =4)
mcmc_acf(draws2, pars =c("b_Intercept", "b_negemot_c", "b_age_c", "b_negemot_c:age_c", "sigma"),lags =4)
The effective sample size is an estimate of the number of independent draws from the posterior distribution of the estimand of interest. Because the draws within a Markov chain are not independent if there is autocorrelation, the effective sample size, \(n_\textit{eff}\), will be smaller than the total sample size, \(N\). The larger the ratio of \(n_\textit{eff}\) to \(N\) the better.
In that quote, they spoke as if there was only one measure of ESS. Though this has been the case for some time, times have changed. In a (2019)paper, Stan-team allstars Vehtari, Gelman, Simpson, Carpenter, and Bürkner proposed two measures of ESS: “bulk-ESS” and “tail-ESS.” From their paper, we read:
if you plan to report quantile estimates or posterior intervals, we strongly suggest assessing the convergence of the chains for these quantiles. In Section 4.3 we show that convergence of Markov chains is not uniform across the parameter space and propose diagnostics and effective sample sizes specifically for extreme quantiles. This is different from the standard ESS estimate (which we refer to as the “bulk-ESS”), which mainly assesses how well the centre of the distribution is resolved. Instead, these “tail-ESS” measures allow the user to estimate the MCSE for interval estimates. (p. 5, emphasis in the original)
You can read the paper for technical details. In short, the Bulk_ESS in the output of brms 2.10.0 and up is what was previously referred to as Eff.Sample, the effective sample size. This indexed the number of effective samples in “the centre of the” posterior distribution (i.e., the posterior mean or median). But since we also care about uncertainty in our models, we also care about how well we have characterized the tails of the posterior distribution in 95% intervals and such. The new Tail_ESS in brms output allows us to gauge the effective sample size for those intervals.
The Bulk_ESS and Tail_ESS values were all well above 2,000 with model9.2 and the autocorrelations were very low, too. model9.1 had higher autocorrelations and lower ESS values. The upshot is that even though we have 4,000 samples for each parameter, those samples don’t necessarily give us the same quality of information fully independent samples would. Bulk_ESS and Tail_ESS values help you determine how concerned you should be. And, as it turns out, things like centering can help increase a models Bulk_ESS and Tail_ESS values.
Wading in further, we can use the neff_ratio() function to collect the \(N_\textit{eff}\) to \(N\) ratio for each model parameter and then use mcmc_neff() to make a visual diagnostic. Here we do so for model9.1 and model9.2.
Although none of the \(N_\textit{eff}\) to \(N\) ratios were in the shockingly-low range for either model, there were substantially higher for model9.2.
In addition to autocorrelations and \(N_\textit{eff}\) to \(N\) ratios, there is also the issue that the parameters in the model can themselves be correlated. If you like a visual approach, you can use brms::pairs() to make quick and dirty histograms for each parameter along with scatter plots showing the shape of their correlations.
Here we’ll take a similar approach, but with functions from the GGally package(Schloerke et al., 2021), which will also compute the point estimates for the bivariate correlations.
When fitting models with HMC, centering can make a difference for the parameter correlations. If you prefer a more numeric approach, vcov() will yield the variance/covariance matrix–or correlation matrix when using correlation = T–for the parameters in a model.
And so wait, what does that even mean for a parameter to correlate with another parameter? you might ask. Fair enough. Let’s compute a parameter correlation step by step. The first step requires as_draws_df().
Now we’ve put our posterior draws into a data object, draws, we can make a scatter plot of two parameters. Here we’ll choose b_negemot and the interaction coefficient, b_negemot:age.
draws %>%ggplot(aes(x = b_negemot, y =`b_negemot:age`)) +geom_point(size =1/10, alpha =1/5) +labs(subtitle ="Each dot is the parameter pair from a\nsingle draw. Looking across the 4,000\ntotal posterior draws, it becomes clear\nthe two parameters are strongly\nnegatively correlated.")
And indeed, the Pearson’s correlation coefficient is strong.
draws %>%summarise(rho =cor(b_negemot, `b_negemot:age`))
# A tibble: 1 × 1
rho
<dbl>
1 -0.956
And what was that part from the vcov() output, again?
This entire topic of HMC diagnostics can seem baffling, especially when compared to the simplicity of OLS. If this is your first introduction, you might want to watch lectures 10 and 11 from McElreath’s Statistical Rethinking Fall 2017 lecture series. Accordingly, you might check out Chapter 8 of his (2015) text, Chapter 9 of his (2020) text, and either of my ebooks translating his work into brms and tidyverse code (Kurz, 2026b, 2026a).
9.1.1 The effect of mean-centering on multicollinearity and the standard error of \(b_3\)
This can be difficult to keep track of, but what we just looked at were the correlations among model parameters. These are not the same as correlations among variables. As such, those correlations are not the same as those in Table 9.2. But we can get those, too. First we’ll have to do a little more data processing to get all the necessary mean-centered variables and standardized variables.
Recall that to get our sweet Bayesian correlations, we use the multivariate mvbind() syntax to fit an intercepts-only model. Here we do that for all three of the Table 9.2 sections.
For the sake of space, I’ll let you check that out for model9.4 and model9.5. If you’re tricky with your VarCorr() indexing, you can also get the model-implied variances.
And if you’re like totally lost with all this indexing, you might execute VarCorr(correlations1) %>% str() and spend a little time looking at what VarCorr() returns.
On page 309, Hayes explained why the OLS variance for \(b_3\) is unaffected by mean centering. The story was similar for our HMC model, too.
9.1.2 The effect of mean-centering on \(b_1\), \(b_2\), and their standard errors posterior \(\textit{SD}\)’s
A second explanation given for why mean-centering is preferred is that it makes \(b_1\) and \(b_2\), the regression coefficients for \(X\) and \(W\), more meaningful. This is generally true and thus not a myth, although it is not necessarily true in all circumstances… .
Mean-centering \(X\) and \(W\) prior to computation of the product and estimation of the model will produce \(b_1\) and \(b_2\) that are always meaningful, rather than meaningful only when \(X\) and/or \(W\) are meaningful when equal to zero… [After mean centering,] \(b_1\) estimates the difference in \(Y\) between two cases that differ by one unit on \(X\) among cases that are average on \(W\). Similarly, \(b_2\) estimates the difference in \(Y\) between two cases that differ by one unit on \(W\) among cases that are average on \(X\). These will always estimate conditional effects of \(X\) on \(Y\) within the range of the data, and they can always be interpreted. (p. 310. emphasis in the original)
If you only care about posterior means, you can use model9.1 to reproduce the results at the bottom of page 310 like this.
But we’re proper Bayesians and like a summary of the spread in the posterior. So we’ll evoke as_draws_df() and the other usual steps, this time just focusing on model9.1.
Note how the standard error Hayes computed at the top of page 311 corresponds nicely with the posterior \(\textit{SD}\) we just computed. Hayes employed a fancy formula; we just used sd(). At any rate, the main message is centering did not effect our estimate of the conditional effect of \(X\). It turns out \((\theta_{X \rightarrow Y} \mid W) = \left (\theta_{\overline X \rightarrow Y} \mid \overline W \right)\).
9.1.3 The centering option in PROCESS
I’m not aware of a similar function in brms. You’ll have to use your data wrangling skills.
9.2 The estimation and interpretation of standardized regression coefficients in a moderation analysis
Mean-centering does nothing to change the scaling of regression coefficients. Whether or not mean-centering is used when estimating a model of the form \(\hat Y = i_Y + b_1 X + b_2 W + b_3 XW\), \(b_1\), \(b_2\), and \(b_3\) are interpreted with respect to the measured metrics of \(X\), \(W\), and \(Y\) (i.e., in unstandardized form). Although [Hayes] generally prefer[s] to report and interpret regression analyses based on unstandardized coefficients, it is possible to generate regression coefficients that are analogous to standardized regression coefficients in regression models without a product term as a predictor. However, one must be careful when doing so. (p. 313, emphasis in the original)
9.2.1 Variant 1
We’ve already computed standardized predictors. Now we just need to standardize the criterion, govact.
Family: gaussian
Links: mu = identity
Formula: govact_z ~ negemot_z + age_z + negemot_x_age_z
Data: glbwarm (Number of observations: 815)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.002 0.028 -0.054 0.056 1.003 2209 2230
negemot_z 0.165 0.092 -0.017 0.344 1.003 1474 1898
age_z -0.369 0.071 -0.506 -0.226 1.003 1492 1853
negemot_x_age_z 0.511 0.111 0.296 0.731 1.004 1429 1770
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.807 0.020 0.769 0.847 1.001 2543 2297
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
The results correspond nicely to those in Table 9.1, too.
9.3 A caution on manual centering and standardization [because of missing data]
It’s worthwhile considering the issue of listwise deletion when data are partially missing. The brms default is to delete rows with missingness, NA in R, for the predictors. However, brms allows users to perform one-step Bayesian imputation for missing values using the mi() syntax. You can learn all the details in Bürkner’s (2022b) vignette, Handle missing values with brms. First we’ll fit see what happens when you fit a model in brms when some of the negemot_z values are missing, but without using the mi() syntax. And of course before we do that, we’ll make a negemot_z_missing variable, which is identical to negemot_z, but about 10% of the values are missing.
If you’ve never used rbinom() before, execute ?rbinom or look it up in your favorite web search engine. Here’s our listwise deletion model, which corresponds to what you’d get from a typical OLS-based program.
Warning: Rows containing NAs were excluded from the model.
Warning: Rows containing NAs were excluded from the model.
Note the warning. Let’s compare the listwise deletion results with the model based on all the data.
print(model9.6)
Family: gaussian
Links: mu = identity
Formula: govact_z ~ negemot_z + age_z + negemot_z:age_z
Data: glbwarm (Number of observations: 815)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.01 0.03 -0.05 0.06 1.00 5380 2903
negemot_z 0.56 0.03 0.51 0.62 1.00 4320 3068
age_z -0.06 0.03 -0.12 -0.01 1.00 4786 3109
negemot_z:age_z 0.13 0.03 0.07 0.19 1.00 4625 3222
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.81 0.02 0.77 0.85 1.00 4859 2512
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
print(model9.8)
Family: gaussian
Links: mu = identity
Formula: govact_z ~ negemot_z_missing + age_z + negemot_z_missing:age_z
Data: glbwarm (Number of observations: 731)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.00 0.03 -0.05 0.06 1.00 4834 3093
negemot_z_missing 0.57 0.03 0.51 0.63 1.00 4589 3116
age_z -0.05 0.03 -0.10 0.01 1.00 5043 2980
negemot_z_missing:age_z 0.12 0.03 0.06 0.18 1.00 5132 3114
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.81 0.02 0.77 0.85 1.00 5265 3238
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
In this case, the model results were similar to those based on all the data because we used rbinom() to delete the predictor values completely at random. With real data and real-live missing data mechanisms, the situation isn’t often so rosy. But anyway, the real story, here, is the Data: glbwarm (Number of observations: n) line at the top of the print() outputs. The number, \(n\), was 815 in the model using all the data and 731 for the one based on listwise deletion. That’s a lot of missing information.
The mi() syntax will allow us to use all the rows in a model, even if one or more of the predictors contain missing values. The syntax makes the model a multivariate model in that now we’ll be modeling both govact_zandnegemot_z_missing. As we’ve covered in prior chapters, there are multiple ways to write a multivariate model in brms(see Bürkner, 2022a). One nice way is to write the model for each criterion separately in a bf() statement. You combine the bf() statements together with the + operator. And for models like the ones in Hayes’s text, you’ll also want to tack on set_rescor(FALSE). You can do this within the brm() function, as usual. But I find that this clutters the code up more than I like. So another approach is to save the combination of bf() statements as an object.
Family: MV(gaussian, gaussian)
Links: mu = identity
mu = identity
Formula: govact_z ~ 1 + mi(negemot_z_missing) + age_z + mi(negemot_z_missing):age_z
negemot_z_missing | mi() ~ 1
Data: glbwarm (Number of observations: 815)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
govactz_Intercept 0.01 0.03 -0.05 0.06 1.00 5688 3130
negemotzmissing_Intercept -0.00 0.04 -0.07 0.07 1.00 5091 3310
govactz_age_z -0.06 0.03 -0.12 -0.00 1.00 6276 2884
govactz_minegemot_z_missing 0.57 0.03 0.51 0.63 1.00 5251 2770
govactz_minegemot_z_missing:age_z 0.13 0.03 0.07 0.18 1.00 5527 3058
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_govactz 0.81 0.02 0.77 0.85 1.00 5216 2937
sigma_negemotzmissing 1.00 0.03 0.95 1.05 1.00 5228 3052
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
When using the multivariate mi() syntax, your print() output becomes more complicated. We got regression sub-models for both govact_z and negemot_z_missing. But look at the line at the top of the output that reads “Data: glbwarm (Number of observations: 815)”. The multivariate mi() syntax used all the available data. No listwise deletion necessary.
The print() output for our model obscured some of the results. To clarify what the mi() syntax did, let’s peek at the first columns returned by as_draws_df().
draws <-as_draws_df(model9.9)draws[, 1:20] %>%glimpse()
Columns b_govactz_Intercept through sigma_negemotzmissing were business as usual. But notice all the Ymi_negemotzmissing[i] columns. In each of these we see 4,000 posterior draws for the missing negemot_z_missing values. The [i] part of the column names indexes which row number the iterations correspond to. Summarizing these columns might help us get a sense of the results.
That’s a lot of output. Here’s what those summaries look like in a coefficient plot.
# summarize just like beforedraws %>%pivot_longer(starts_with("Ymi"),names_to ="row") %>%group_by(row) %>%mean_qi(value) %>%# plot!ggplot(aes(x = row, y = value, ymin = .lower, ymax = .upper)) +geom_hline(yintercept =0, color ="grey75") +geom_pointinterval(linewidth =1/3) +scale_x_discrete("rank-ordered row number", breaks =NULL) +ylab("imputed value")
Each missing negemot_z value got an entire posterior distribution. Just as the model is uncertain about what those values might have been, that uncertainty was baked right into the primary submodel predicting govact_z. That’s a good thing. We want our Bayesian models to use as much information as they can and yield results with as much certainty as possible. But we don’t want our models to be more certain than the data–and priors–allow. When you use listwise deletion methods, you leave information on the table, which we don’t want. But when you use old fashioned ad hock methods like mean imputation, you underestimate the uncertainty in the model, which we also don’t want. We want the middle path.
Here’s a focused look at two other important new parameters.
Our model has estimated the mean and standard deviations for our negemot_z_missing variable. Hopefully that isn’t a surprise. This is exactly what we asked brms to do with the negemot_z_missing | mi() ~ 1 part of the model formula. Since that submodel had no predictors, the intercept was just the mean. Correspondingly, the residual variance was the entire variance–but expressed in the usual \(\sigma\) metric since we’re using brms. And since our negemot_z_missing variable is a subset of the standardized negemot_z variable, naturally the estimates for the mean and standard deviation are about 0 and 1, respectively.
Another method we could have used is multiple imputation. With multiple imputation, you create a small number of alternative data sets, typically 5, into which you impute plausible values into the missing value slots. You then fit the model separately for each of the imputed data sets. Because the imputed values will vary across the data sets, that uncertainty will get appropriately transmitted to the model. I know this might sound crazy, but it typically leads to much lower model bias when compared to mean imputation or listwise deletion.
But we didn’t quite use multiple imputation. With one-step Bayesian imputation using the mi() syntax, you get an entire posterior distribution for each missing value. And if you have variables in the data set that might help predict what those missing values are, you’d just plug that into the missing data submodel. Improving the imputation model can improve the subsequent substantive model.
For more on the mi() approach, see Bürkner’s (2022b)vignette. McElreath lectured on this topic (here and here) within the context of his Statistical rethinking texts (McElreath, 2015, 2020), and you can find links to my translations of both of his texts here. For a more general introduction to missing data theory, check out van Burren’s (2018) book, Flexible imputation of missing data, or Enders’s great (2022) text, Applied missing data analysis. You can also find Enders lecturing on missing data here.
The take home message is there is no need to ignore missing data or use outdated procedures like listwise deletion. Be a champion and model your missing data with brms.
9.4 More than one moderator
None of this is a problem for brms. But instead of using the model=i syntax in Hayes’s PROCESS, you just have to specify your model formula in brm().
9.4.1 Additive multiple moderation
It’s trivial to add sex, its interaction with negemot, and the two covariates (i.e., posemot and ideology) to the model. We can even do it within update().
Our output matches nicely with the formula at the bottom of page 232 and the PROCESS output in Figure 9.2.
print(model9.10, digits =3)
Family: gaussian
Links: mu = identity
Formula: govact ~ negemot + sex + age + posemot + ideology + negemot:sex + negemot:age
Data: glbwarm (Number of observations: 815)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 5.273 0.333 4.635 5.941 1.000 2492 2414
negemot 0.092 0.082 -0.073 0.252 1.002 2328 2799
sex -0.750 0.195 -1.130 -0.366 1.003 2127 2693
age -0.018 0.006 -0.030 -0.006 1.000 2310 2432
posemot -0.023 0.027 -0.078 0.030 1.001 4039 2594
ideology -0.207 0.026 -0.257 -0.155 1.000 3825 2496
negemot:sex 0.206 0.050 0.109 0.301 1.003 2133 2624
negemot:age 0.005 0.002 0.002 0.008 1.001 2290 2513
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 1.048 0.026 0.998 1.100 1.000 3710 2876
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
On page 325, Hayes discussed the unique variance each of the two moderation terms accounted for after controlling for the other covariates. In order to get our Bayesian version of these, we’ll have to fit two additional models, one after removing each of the interaction terms.
model9.11<-update( model9.10, newdata = glbwarm, govact ~1+ negemot + sex + age + posemot + ideology + negemot:sex,cores =4,file ="fits/model09.11")model9.12<-update( model9.10, newdata = glbwarm, govact ~1+ negemot + sex + age + posemot + ideology + negemot:age,cores =4,file ="fits/model09.12")
Here we’ll extract the bayes_R2() iterations for each of the three models, place them all in a single tibble, and then do a little arithmetic to get the difference scores. After all that data wrangling, we’ll summarize() as usual.
r2 <-tibble(r2_without_age_interaction =bayes_R2(model9.11, summary = F)[, 1],r2_without_sex_interaction =bayes_R2(model9.12, summary = F)[, 1],r2_with_both_interactions =bayes_R2(model9.10, summary = F)[, 1]) %>%mutate(`delta R2 due to age interaction`= r2_with_both_interactions - r2_without_age_interaction,`delta R2 due to sex interaction`= r2_with_both_interactions - r2_without_sex_interaction)r2 %>%pivot_longer(`delta R2 due to age interaction`:`delta R2 due to sex interaction`) %>%group_by(name) %>%mean_qi(value) %>%mutate_if(is.double, round, digits =3)
# A tibble: 2 × 7
name value .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 delta R2 due to age interaction 0.007 -0.049 0.063 0.95 mean qi
2 delta R2 due to sex interaction 0.012 -0.044 0.068 0.95 mean qi
Recall that \(R^2\) is in a 0-to-1 metric. It’s a proportion. If you want to convert that to a percentage, as in percent of variance explained, you’d just multiply by 100. To make it explicit, let’s do that.
r2 %>%pivot_longer(`delta R2 due to age interaction`:`delta R2 due to sex interaction`) %>%group_by(name) %>%mutate(percent = value *100) %>%mean_qi(percent) %>%mutate_if(is.double, round, digits =3)
# A tibble: 2 × 7
name percent .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 delta R2 due to age interaction 0.681 -4.94 6.26 0.95 mean qi
2 delta R2 due to sex interaction 1.16 -4.40 6.85 0.95 mean qi
Hopefully it’s clear how our proportions turned percentages correspond to the numbers on page 325. However, note how our 95% credible intervals do not cohere with the \(p\)-values from Hayes’s \(F\)-tests.
If we want to prep for our version of Figure 9.3, we’ll need to carefully specify the predictor values we’ll pass through the fitted() function. Here we do so and save them in nd.
tibble [180 × 5] (S3: tbl_df/tbl/data.frame)
$ negemot : num [1:180] 0.5 0.5 0.5 0.5 0.5 ...
$ sex : int [1:180] 0 0 0 1 1 1 0 0 0 1 ...
$ age : num [1:180] 30 50 70 30 50 70 30 50 70 30 ...
$ posemot : num [1:180] 3.13 3.13 3.13 3.13 3.13 ...
$ ideology: num [1:180] 4.08 4.08 4.08 4.08 4.08 ...
With our nd values in hand, we’re ready to make our version of Figure 9.3.
fitted(model9.10, newdata = nd) %>%data.frame() %>%bind_cols(nd) %>%# these lines will make the strip text match with those with Hayes's Figuremutate(sex =if_else(sex ==0, str_c("Females, W = ", sex),str_c("Males, W = ", sex)),age =str_c("Age, Z, = ", age)) %>%# finally, we plot!ggplot(aes(x = negemot, group = sex)) +geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5, fill = sex),alpha =1/3, color ="transparent") +geom_line(aes(y = Estimate, color = sex),linewidth =1) +scale_fill_manual(values =colorblind_pal()(8)[2:3]) +scale_color_manual(values =colorblind_pal()(8)[2:3]) +scale_x_continuous(breaks =1:6) +coord_cartesian(xlim =c(1, 6),ylim =c(3, 6)) +labs(x =expression("Negative Emotions about Climate Change, "*italic(X)),y =expression("Support for Government Action to Mitigate Climate Change, "*italic(Y))) +theme(legend.position ="top",legend.title =element_blank()) +facet_grid(age ~ .)
Note how we’ve colored the lines and ribbons with help from the handy ggthemes::colorblind_pal() function, which provides colorblind-friendly colors for categorical variables.
Anyway, recall that the conditional effect of \(X\) for various values of \(W\) and \(Z\) is
This is easiest to show with as_draws_df() and a little algebra. As in the middle of page 329, here we solve for both sexes and age == 50.
draws <-as_draws_df(model9.10)draws %>%# algebramutate(men = b_negemot +`b_negemot:sex`*1+`b_negemot:age`*50,women = b_negemot +`b_negemot:sex`*0+`b_negemot:age`*50) %>%# more algebramutate(`men - women`= men - women) %>%pivot_longer(men:`men - women`) %>%# this just orders the outputmutate(name =factor(name, levels =c("men", "women", "men - women"))) %>%group_by(name) %>%mean_qi(value) %>%mutate_if(is.double, round, digits =3)
# A tibble: 3 × 7
name value .lower .upper .width .point .interval
<fct> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 men 0.536 0.467 0.606 0.95 mean qi
2 women 0.33 0.259 0.399 0.95 mean qi
3 men - women 0.206 0.109 0.301 0.95 mean qi
Switching our discussion to \(b_5\) (i.e., \(b_{\text{negemot} \times \text{age}}\)), Hayes showed its value for two groups 10 years apart. Here it is for model9.10.
draws %>%# algebratransmute(`10 * b_5`=`b_negemot:age`*10) %>%mean_qi() %>%mutate_if(is.double, round, digits =3)
Because we haven’t changed the predictor variables in the model–just added interactions among them–there’s no need to redo our nd values. Rather, all we need to do is pass them through fitted() based on our new model9.13 and plot. Without further ado, here’s our Figure 9.6.
fitted(model9.13, newdata = nd) %>%data.frame() %>%bind_cols(nd) %>%# these lines will make the strip text match with those with Hayes's Figuremutate(sex =if_else(sex ==0, str_c("Females, W = ", sex),str_c("Males, W = ", sex)),age =str_c("Age, Z, = ", age)) %>%# behold, Figure 9.6!ggplot(aes(x = negemot, group = sex)) +geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5, fill = sex),alpha =1/3, color ="transparent") +geom_line(aes(y = Estimate, color = sex),linewidth =1) +scale_fill_manual(values =colorblind_pal()(8)[2:3]) +scale_color_manual(values =colorblind_pal()(8)[2:3]) +scale_x_continuous(breaks =1:6) +coord_cartesian(xlim =c(1, 6),ylim =c(3, 6)) +labs(x =expression("Negative Emotions about Climate Change, "*italic(X)),y =expression("Support for Government Action to Mitigate Climate Change, "*italic(Y))) +theme(legend.position ="top",legend.title =element_blank()) +facet_grid(age ~ .)
For the pick-a-point values Hayes covered on page 338, recall that when using as_draws_df(), our \(b_4\) is b_negemot:sex and our \(b_7\) is b_negemot:sex:age.
# A tibble: 3 × 7
theta_XW_on_Y_given value .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 age = 30 0.07 -0.092 0.227 0.95 mean qi
2 age = 50 0.204 0.105 0.305 0.95 mean qi
3 age = 70 0.338 0.177 0.504 0.95 mean qi
Our method for making a JN technique plot with fitted() way back in Chapter 7 isn’t going to work, here. At least not as far as I can see. Rather, we’re going to have to skillfully manipulate our draws object. For those new to R, this might be a little confusing at first. So I’m going to make a crude attempt first and then get more sophisticated.
Crude attempt:
draws %>%mutate(`age = 30`=`b_negemot:sex`+`b_negemot:sex:age`*30, `age = 50`=`b_negemot:sex`+`b_negemot:sex:age`*50, `age = 70`=`b_negemot:sex`+`b_negemot:sex:age`*70) %>%pivot_longer(contains("="), names_to ="theta_XW_on_Y_given") %>%mutate(`theta XW on Y given`=str_extract(theta_XW_on_Y_given, "\\d+") %>%as.double()) %>%group_by(`theta XW on Y given`) %>%mean_qi(value) %>%# the plotggplot(aes(x =`theta XW on Y given`)) +geom_hline(yintercept =0) +geom_vline(xintercept =38.114) +geom_ribbon(aes(ymin = .lower, ymax = .upper),alpha =1/2) +geom_line(aes(y = value), linewidth =1) +coord_cartesian(xlim =c(20, 85),ylim =c(-0.25, 0.75))
Notice how we just took the code from our pick-a-point analysis and dumped it into a plot. So one obvious approach would be to pick like 30 or 50 age values to plug into transmute(), and just do the same thing. If you’re super afraid of coding, that’d be one intuitive but extremely verbose attempt. I have done stuff like that earlier in my R career. There’s no shame in being extremely verbose and redundant if that’s what makes sense. Another way is to think in terms of functions. When we made age = 30 within transmute(), we took a specific age value (i.e., 30) and plugged it into the formula b_negemot:sex + b_negemot:sex:age * i where \(i\) = 30. And when we made age = 50 we did exactly the same thing but switched out the 30 for a 50. So what we need is a function that will take a range of values for \(i\), plug them into our b_negemot:sex + b_negemot:sex:age * i formula, and then neatly return the output. A nice base R function for that is sapply().
'data.frame': 4000 obs. of 76 variables:
$ X1 : num 0.17222 -0.10668 0.01881 0.00871 0.22339 ...
$ X2 : num 0.1748 -0.0971 0.0261 0.0141 0.2225 ...
$ X3 : num 0.1774 -0.0875 0.0333 0.0196 0.2217 ...
$ X4 : num 0.18 -0.0779 0.0406 0.025 0.2208 ...
$ X5 : num 0.1826 -0.0682 0.0478 0.0304 0.22 ...
$ X6 : num 0.1852 -0.0586 0.0551 0.0358 0.2191 ...
$ X7 : num 0.1878 -0.049 0.0623 0.0413 0.2182 ...
$ X8 : num 0.1904 -0.0394 0.0696 0.0467 0.2174 ...
$ X9 : num 0.1931 -0.0298 0.0769 0.0521 0.2165 ...
$ X10: num 0.1957 -0.0202 0.0841 0.0576 0.2157 ...
$ X11: num 0.1983 -0.0106 0.0914 0.063 0.2148 ...
$ X12: num 0.200861 -0.000999 0.098618 0.06842 0.213938 ...
$ X13: num 0.20346 0.00861 0.10587 0.07385 0.21308 ...
$ X14: num 0.2061 0.0182 0.1131 0.0793 0.2122 ...
$ X15: num 0.2087 0.0278 0.1204 0.0847 0.2114 ...
$ X16: num 0.2113 0.0374 0.1276 0.0901 0.2105 ...
$ X17: num 0.2139 0.047 0.1349 0.0956 0.2096 ...
$ X18: num 0.2165 0.0566 0.1422 0.101 0.2088 ...
$ X19: num 0.2191 0.0663 0.1494 0.1064 0.2079 ...
$ X20: num 0.2217 0.0759 0.1567 0.1118 0.2071 ...
$ X21: num 0.2243 0.0855 0.1639 0.1173 0.2062 ...
$ X22: num 0.2269 0.0951 0.1712 0.1227 0.2053 ...
$ X23: num 0.229 0.105 0.178 0.128 0.204 ...
$ X24: num 0.232 0.114 0.186 0.134 0.204 ...
$ X25: num 0.235 0.124 0.193 0.139 0.203 ...
$ X26: num 0.237 0.134 0.2 0.144 0.202 ...
$ X27: num 0.24 0.143 0.207 0.15 0.201 ...
$ X28: num 0.243 0.153 0.215 0.155 0.2 ...
$ X29: num 0.245 0.162 0.222 0.161 0.199 ...
$ X30: num 0.248 0.172 0.229 0.166 0.198 ...
$ X31: num 0.25 0.182 0.236 0.172 0.198 ...
$ X32: num 0.253 0.191 0.244 0.177 0.197 ...
$ X33: num 0.256 0.201 0.251 0.182 0.196 ...
$ X34: num 0.258 0.21 0.258 0.188 0.195 ...
$ X35: num 0.261 0.22 0.265 0.193 0.194 ...
$ X36: num 0.263 0.23 0.273 0.199 0.193 ...
$ X37: num 0.266 0.239 0.28 0.204 0.192 ...
$ X38: num 0.269 0.249 0.287 0.21 0.192 ...
$ X39: num 0.271 0.258 0.295 0.215 0.191 ...
$ X40: num 0.274 0.268 0.302 0.22 0.19 ...
$ X41: num 0.276 0.278 0.309 0.226 0.189 ...
$ X42: num 0.279 0.287 0.316 0.231 0.188 ...
$ X43: num 0.282 0.297 0.324 0.237 0.187 ...
$ X44: num 0.284 0.306 0.331 0.242 0.186 ...
$ X45: num 0.287 0.316 0.338 0.248 0.186 ...
$ X46: num 0.289 0.326 0.345 0.253 0.185 ...
$ X47: num 0.292 0.335 0.353 0.258 0.184 ...
$ X48: num 0.295 0.345 0.36 0.264 0.183 ...
$ X49: num 0.297 0.354 0.367 0.269 0.182 ...
$ X50: num 0.3 0.364 0.374 0.275 0.181 ...
$ X51: num 0.302 0.374 0.382 0.28 0.18 ...
$ X52: num 0.305 0.383 0.389 0.286 0.18 ...
$ X53: num 0.308 0.393 0.396 0.291 0.179 ...
$ X54: num 0.31 0.403 0.403 0.296 0.178 ...
$ X55: num 0.313 0.412 0.411 0.302 0.177 ...
$ X56: num 0.315 0.422 0.418 0.307 0.176 ...
$ X57: num 0.318 0.431 0.425 0.313 0.175 ...
$ X58: num 0.321 0.441 0.432 0.318 0.174 ...
$ X59: num 0.323 0.451 0.44 0.324 0.174 ...
$ X60: num 0.326 0.46 0.447 0.329 0.173 ...
$ X61: num 0.328 0.47 0.454 0.334 0.172 ...
$ X62: num 0.331 0.479 0.461 0.34 0.171 ...
$ X63: num 0.334 0.489 0.469 0.345 0.17 ...
$ X64: num 0.336 0.499 0.476 0.351 0.169 ...
$ X65: num 0.339 0.508 0.483 0.356 0.168 ...
$ X66: num 0.341 0.518 0.49 0.362 0.168 ...
$ X67: num 0.344 0.527 0.498 0.367 0.167 ...
$ X68: num 0.347 0.537 0.505 0.372 0.166 ...
$ X69: num 0.349 0.547 0.512 0.378 0.165 ...
$ X70: num 0.352 0.556 0.519 0.383 0.164 ...
$ X71: num 0.354 0.566 0.527 0.389 0.163 ...
$ X72: num 0.357 0.575 0.534 0.394 0.162 ...
$ X73: num 0.36 0.585 0.541 0.4 0.162 ...
$ X74: num 0.362 0.595 0.548 0.405 0.161 ...
$ X75: num 0.365 0.604 0.556 0.41 0.16 ...
$ X76: num 0.367 0.614 0.563 0.416 0.159 ...
Okay, so that looks a little monstrous. What we did in the first argument was tell sapply() which values we’d like to use in some function. We chose each integer ranging from 15 to 90–which, if you do the math, is 76 values. We then told sapply() to plug those values into a custom function, which we defined as function(i) {draws$b_negemot:sex + draws$b_negemot:sex:age * i}. In our custom function, i was a placeholder for each of those 76 integers. But remember that draws has 4,000 rows, each one corresponding to one of the 4,000 posterior iterations. Thus, for each of our 76 i-values, we got 4,000 results. After all that sapply() returned a matrix. Since we like to work within the tidyverse and use ggplot2, we just went ahead and put those results in a tibble.
With our sapply() output in hand, all we need to do is a little more indexing and summarizing and we’re ready to plot. The result is our very own version of Figure 9.7.
Or for kicks and giggles, another way to get a clearer sense of how our data informed the shape of the plot, here we replace our geom_ribbon() + geom_line() code with geom_pointrange().
Although I probably wouldn’t try to use a plot like this in a manuscript, I hope it makes clear how the way we’ve been implementing the JN technique is just the pick-a-point approach in bulk. No magic, here.
For all you tidyverse fanatics out there, don’t worry. There are more tidyverse-centric ways to get the plot values than with sapply(). We’ll get to them soon enough. It’s advantageous to have good old base Rsapply() up your sleeve, too. And new R users, it’s helpful to know that sapply() is one part of the apply() family of base R functions, about which you might learn more here or here.
Now the conditional effect of \(X\) on \(Y\) given \(W\) and \(Z\) is
\[\theta_{X \rightarrow Y} = b_1 + b_4 W + b_5 Z + b_7 WZ.\]
In the terms of model9.13 where sex = \(W\) and age = \(Z\), we can restate that as
Following Hayes at the bottom of page 341, here is the difference in the effect of negative emotions between men and women among 30-year-olds, \(\theta_{XW \rightarrow Y} \mid (Z = 30)\).
draws %>%# algebramutate(men = b_negemot +`b_negemot:sex`*1+`b_negemot:age`*30+`b_negemot:sex:age`*1*30,women = b_negemot +`b_negemot:sex`*0+`b_negemot:age`*30+`b_negemot:sex:age`*0*30) %>%# more algebramutate(`men - women`= men - women) %>%pivot_longer(men:`men - women`) %>%# this just orders the outputmutate(name =factor(name, levels =c("men", "women", "men - women"))) %>%group_by(name) %>%mean_qi(value) %>%mutate_if(is.double, round, digits =3)
# A tibble: 3 × 7
name value .lower .upper .width .point .interval
<fct> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 men 0.371 0.243 0.493 0.95 mean qi
2 women 0.301 0.195 0.407 0.95 mean qi
3 men - women 0.07 -0.092 0.227 0.95 mean qi
In contrast, here is \(\theta_{XW \rightarrow Y} \mid (Z = 50)\).
# A tibble: 3 × 7
name value .lower .upper .width .point .interval
<fct> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 men 0.522 0.452 0.593 0.95 mean qi
2 women 0.318 0.244 0.393 0.95 mean qi
3 men - women 0.204 0.105 0.305 0.95 mean qi
Finally, here we compute \(\theta_{XW \rightarrow Y} \mid (Z = 70)\).
# A tibble: 3 × 7
name value .lower .upper .width .point .interval
<fct> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 men 0.674 0.582 0.77 0.95 mean qi
2 women 0.336 0.202 0.464 0.95 mean qi
3 men - women 0.338 0.177 0.504 0.95 mean qi
9.5 Comparing conditional effects
In this section, [Hayes discussed] a test of the difference between two conditional effects in a moderation model. The method [Hayes described] is called a “slope difference test” by Dawson (Dawson, 2014; Dawson & Richter, 2006), who offers some guidelines on how to conduct this test for the moderated moderation model. (p. 342)
As usual, we will follow along, but with alternative methods to the statistical testing paradigm. We will steadfastly continue summarizing and plotting the posterior distribution from various vantage points.
9.5.1 Comparing conditional effects in the additive multiple moderation model
When we’re interested in estimating the difference in the conditional effect of \(X\) on \(Y\) (i.e., \(\theta_{X \rightarrow Y}\)) for \(W = w_1\) and \(Z = z_1\) versus when \(W = w_2\) and \(Z = z_2\), that follows the formula
As in other cases, we don’t have to worry about special considerations for computing the standard errors for out Bayesian models. All we need to do is follow the simple algebraic manipulations of the posterior distribution. Because of the correlation structure within the parameters, the uncertainty in the conditional distribution will work itself out.
9.5.2 Comparing conditional effects in the moderated moderation model
Since we don’t have the contrast feature automated like in PROCESS, we’ll have to carefully follow the equations above to specify the values properly in R. Here we’ll use the equation in the first line,
# A tibble: 3 × 5
name mean sd ll ul
<chr> <dbl> <dbl> <dbl> <dbl>
1 30-year-old men 0.371 0.0623 0.244 0.493
2 50-year-old women 0.318 0.0372 0.244 0.393
3 contrast 0.0523 0.0713 -0.0882 0.191
Notice how our posterior \(\textit{SD}\) corresponded nicely to the standard error in Hayes’s contrast test. And we didn’t even have to worry about using the frightening Formula 9.21 on page 345. That information was contained in the posterior distribution all along. All we had to do was combine the parameter iterations with a little algebra and then summarize().
For good measure, we’ll compute using the equation in the second line,
Cohen, J., Cohen, P., West, S. G., & Aiken, L. S. (2003). Applied multiple regression/correlation analysis for the behavioral sciences (3rd edition). Lawrence Erlbaum Associates. https://doi.org/10.4324/9780203774441
Dawson, J. F. (2014). Moderation in management research: What, why, when, and how. Journal of Business and Psychology, 29(1), 1–19. https://doi.org/10.1007/s10869-013-9308-7
Dawson, J. F., & Richter, A. W. (2006). Probing three-way interactions in moderated multiple regression: Development and application of a slope difference test. Journal of Applied Psychology, 91(4), 917–926. https://doi.org/10.1037/0021-9010.91.4.917
Schloerke, B., Crowley, J., Di Cook, Briatte, F., Marbach, M., Thoen, E., Elberg, A., & Larmarange, J. (2021). GGally: Extension to ’ggplot2’. https://CRAN.R-project.org/package=GGally
Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. (2019). Rank-normalization, folding, and localization: An improved for assessing convergence of MCMC. https://arxiv.org/abs/1903.08008?
Comments