Centering in Moderation Analysis: A Guide

R
statistics
code
datawaizard
moderation
interactions
Author

Mattan S. Ben-Shachar

Published

January 21, 2024

This blog post will cover what centering is, what a sum-to-zero contrast is, why and when you should use them, and how you can do them in R. Specifically, we will focus on their use and how they affect coefficient interpretation in regression tables and affect interpretation in ANOVA tables.

I will be using the following packages and versions:

library(datawizard) # 0.9.1
library(parameters) # 0.21.3

We will be using the following toy data set:

Generate Toy Dataset 
data <- data.frame(
  Y = c(21, 21, 22.8, 21.4, 18.7, 18.1, 14.3, 24.4, 22.8, 19.2), 
  G = factor(c("g2", "g2", "g1", "g2", "g3", "g2", "g3", "g1", "g1", "g2")), 
  X = c(160, 160, 108, 258, 360, 225, 360, 146.7, 140.8, 167.6),
  Q = c(3.9, 3.9, 3.85, 3.08, 3.15, 2.76, 3.21, 3.69, 3.92, 3.92), 
  A = factor(c("a1", "a1", "a2", "a2", "a1", "a2", "a1", "a2", "a2", "a2"))
)

(This is really a re-labeled version of mtcars[1:10, c(1:3, 5, 8)].)

str(data)
#> 'data.frame':    10 obs. of  5 variables:
#>  $ Y: num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2
#>  $ G: Factor w/ 3 levels "g1","g2","g3": 2 2 1 2 3 2 3 1 1 2
#>  $ X: num  160 160 108 258 360 ...
#>  $ Q: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92
#>  $ A: Factor w/ 2 levels "a1","a2": 1 1 2 2 1 2 1 2 2 2

Simple Regression

We will start with not centering anything.

Continuous Predictor

model1 <- lm(Y ~ X, 
             data = data)

Something that is always true - the (Intercept) represents the predicted values when all the predictors are set to 0.

model_parameters(model1)
#> Parameter   | Coefficient |       SE |         95% CI |  t(8) |      p
#> ----------------------------------------------------------------------
#> (Intercept) |       25.56 |     1.62 | [21.82, 29.31] | 15.74 | < .001
#> X           |       -0.02 | 7.20e-03 | [-0.04, -0.01] | -3.46 | 0.009

# Predict for X=0
pred_grid <- data.frame(X = 0)
predict(model1, newdata = pred_grid)
#>       1 
#> 25.5638

But say we wanted the intercept to represent the predicted value when \(x_i=\bar{X}\) - how would we do that? By centering!

So what is centering?

Don’t just be a zero…

Centering is the processes of assigning meaning to the value of 0 of some variable. We do this by:

  1. defining a meaningful value
  2. Subtracting that value from all values of the variable.

For example, I can center X so that 0 is 258.

rbind(
  X = data$X,
  "X - 258" = data$X - 258
)
#>         [,1] [,2] [,3] [,4] [,5] [,6] [,7]   [,8]   [,9] [,10]
#> X        160  160  108  258  360  225  360  146.7  140.8 167.6
#> X - 258  -98  -98 -150    0  102  -33  102 -111.3 -117.2 -90.4

We can see that by subtracting 258 from X, values that were above 258 are positive, values that were below 258 are negative, and the values that were 258 are now 0. Thus, for this new variable, when we speak of 0, we actually mean “258”.

model2 <- lm(Y ~ I(X - 258),
             data = data)

model_parameters(model2)
#> Parameter   | Coefficient |       SE |         95% CI |  t(8) |      p
#> ----------------------------------------------------------------------
#> (Intercept) |       19.14 |     0.71 | [17.50, 20.78] | 26.86 | < .001
#> X - 258     |       -0.02 | 7.20e-03 | [-0.04, -0.01] | -3.46 | 0.009

# Predict for X=258
pred_grid <- data.frame(X = 258)
predict(model1, newdata = pred_grid)
#>        1 
#> 19.14033

Typically, centering is done around the mean (mean-centering) - thus giving 0 the meaning of “mean”.

model3 <- lm(Y ~ I(X - mean(X)),
             data = data)

model_parameters(model3)
#> Parameter   | Coefficient |       SE |         95% CI |  t(8) |      p
#> ----------------------------------------------------------------------
#> (Intercept) |       20.37 |     0.62 | [18.95, 21.79] | 32.99 | < .001
#> X - mean(X) |       -0.02 | 7.20e-03 | [-0.04, -0.01] | -3.46 | 0.009

# Predict for X=mean(X)
pred_grid <- data.frame(X = mean(data$X))
predict(model1, newdata = pred_grid)
#>     1 
#> 20.37

We can also do this with the scale() function, but since this function also standardizes (re-scales to and sd of 1), we need to set scale(scale = FALSE):

model4 <- lm(Y ~ scale(X, scale = FALSE),
             data = data)

model_parameters(model4)
#> Parameter   | Coefficient |       SE |         95% CI |  t(8) |      p
#> ----------------------------------------------------------------------
#> (Intercept) |       20.37 |     0.62 | [18.95, 21.79] | 32.99 | < .001
#> X           |       -0.02 | 7.20e-03 | [-0.04, -0.01] | -3.46 | 0.009

# Predict for X=mean(X)
pred_grid <- data.frame(X = mean(data$X))
predict(model1, newdata = pred_grid)
#>     1 
#> 20.37

Alternatively, we can just use datawizard::center():1

model5 <- lm(Y ~ center(X),
             data = data)

model_parameters(model5)
#> Parameter   | Coefficient |       SE |         95% CI |  t(8) |      p
#> ----------------------------------------------------------------------
#> (Intercept) |       20.37 |     0.62 | [18.95, 21.79] | 32.99 | < .001
#> center(X)   |       -0.02 | 7.20e-03 | [-0.04, -0.01] | -3.46 | 0.009

# Predict for X=mean(X)
pred_grid <- data.frame(X = mean(data$X))
predict(model1, newdata = pred_grid)
#>     1 
#> 20.37

Note that no matter how we centered X, its coefficient did not change - only the intercept changed!

Categorical Predictor

When adding categorical predictors (in R we call these “factors”), R converts them to a set of \(k-1\) variables (\(k\) being the number of levels in the factor) that represent a set of contrasts between the levels of the variable. For our G factor that has 3 levels, we get two such variables: G[g2] and G[g3].

model6 <- lm(Y ~ G,
             data = data)

model_parameters(model6)
#> Parameter   | Coefficient |   SE |          95% CI |  t(7) |      p
#> -------------------------------------------------------------------
#> (Intercept) |       23.33 | 0.96 | [ 21.05, 25.61] | 24.21 | < .001
#> G [g2]      |       -3.19 | 1.22 | [ -6.08, -0.31] | -2.62 | 0.034 
#> G [g3]      |       -6.83 | 1.52 | [-10.44, -3.23] | -4.49 | 0.003

By default, R uses treatment contrasts (stats::contr.treatment()). These are built up by constructing each variable to be 0 for all but 1 level of the factor, with the 1st level being left out. We can see this by playing with stats::contr.treatment():

contr.treatment(levels(data$G))
#>    g2 g3
#> g1  0  0
#> g2  1  0
#> g3  0  1

The “variable” (called a dummy variable) of g2 is only 1 for the level of G=g2, the dummy variable of g3 is only 1 for the level of G=g3, while g1 is left out - called the “reference level”.

Recall that the (Intercept) always represents the predicted values when all the predictors are set to 0. In this case, setting both g2 and g3 to 0 is the same as representing G=g1. And indeed, that is what the intercept represents:

model_parameters(model6)
#> Parameter   | Coefficient |   SE |          95% CI |  t(7) |      p
#> -------------------------------------------------------------------
#> (Intercept) |       23.33 | 0.96 | [ 21.05, 25.61] | 24.21 | < .001
#> G [g2]      |       -3.19 | 1.22 | [ -6.08, -0.31] | -2.62 | 0.034 
#> G [g3]      |       -6.83 | 1.52 | [-10.44, -3.23] | -4.49 | 0.003

# Predict for G=g1
pred_grid <- data.frame(G = "g1")
predict(model6, newdata = pred_grid)
#>        1 
#> 23.33333

Interpreting Dummy Variables

The other variables now represent the differences between each level and the reference level. It is hard to see, but we can use some linear algebra to help us here:

(mm1 <- cbind(
  `(Intercept)` = 1,
  contr.treatment(levels(data$G))
))
#>    (Intercept) g2 g3
#> g1           1  0  0
#> g2           1  1  0
#> g3           1  0  1

t(solve(mm1))
#>    (Intercept) g2 g3
#> g1           1 -1 -1
#> g2           0  1  0
#> g3           0  0  1

Each column now holds the linear contrast weights represented by that coefficient. We can confirm this:

model_parameters(model6)
#> Parameter   | Coefficient |   SE |          95% CI |  t(7) |      p
#> -------------------------------------------------------------------
#> (Intercept) |       23.33 | 0.96 | [ 21.05, 25.61] | 24.21 | < .001
#> G [g2]      |       -3.19 | 1.22 | [ -6.08, -0.31] | -2.62 | 0.034 
#> G [g3]      |       -6.83 | 1.52 | [-10.44, -3.23] | -4.49 | 0.003

# Predict for G = g1, g2, g3
pred_grid <- data.frame(G = c("g1", "g2", "g3"))
(predicted_means <- predict(model6, newdata = pred_grid))
#>        1        2        3 
#> 23.33333 20.14000 16.50000

# Treatment contrasts
predicted_means[2] - predicted_means[1]
#>         2 
#> -3.193333
predicted_means[3] - predicted_means[1]
#>         3 
#> -6.833333

Sum-to-Zero Contrasts

But say we wanted the intercept to represent the predicted value that is the grand mean - averaged across all levels of G (\(\frac{1}{k}\sum{y|G=g_i}\)) - how would we do that? By using sum-to-zero contrasts!

Sum-to-zero contrasts, like treatment coding (a.k.a dummy coding) convert a factor to a set of contrasts between the levels of the variable. Unlike dummy coding, the weights of each contrast sum to 0.

There are many coding schemes that achieve this (Helmert contrasts, polynomial contrasts, …). The most popular one is effect coding (in R, confusingly built with stats::contr.sum()).

contr.sum(levels(data$G))
#>    [,1] [,2]
#> g1    1    0
#> g2    0    1
#> g3   -1   -1

We can see that each contrast (column) sums to 0.

Side note about changing factor coding

There are 3 ways of changing the coding scheme (contrast matrix) of a factor in R:

# 1. assign to the contrasts() function
contrasts(data$G) <- contr.sum

# 2. The C() function (BIG C!)
data$G <- C(data$G, contr = contr.sum)

# 3. In the model
lm(Y ~ G,
   data = data,
   contrasts = list(G = contr.sum))

These should all produce the same results.

I will be using the 3rd option here for verbosity.

Let’s see how this affects our coefficients:

model7 <- lm(Y ~ G,
             data = data,
             contrasts = list(G = contr.sum))

model_parameters(model7)
#> Parameter   | Coefficient |   SE |         95% CI |  t(7) |      p
#> ------------------------------------------------------------------
#> (Intercept) |       19.99 | 0.57 | [18.65, 21.33] | 35.35 | < .001
#> G [1]       |        3.34 | 0.79 | [ 1.47,  5.22] |  4.21 | 0.004 
#> G [2]       |        0.15 | 0.71 | [-1.53,  1.83] |  0.21 | 0.840

# Predict for G = g1, g2, g3
pred_grid <- data.frame(G = c("g1", "g2", "g3"))
(predicted_means <- predict(model6, newdata = pred_grid))
#>        1        2        3 
#> 23.33333 20.14000 16.50000

# Mean prediction
mean(predicted_means)
#> [1] 19.99111

Success - the intercept is equal to the grand mean!

What about the other variables? It is hard to see, but we can again use some linear algebra to help us here:

(mm2 <- cbind(
  `(Intercept)` = 1,
  contr.sum(levels(data$G))
))
#>    (Intercept)      
#> g1           1  1  0
#> g2           1  0  1
#> g3           1 -1 -1

t(solve(mm2))
#>    (Intercept)                      
#> g1   0.3333333  0.6666667 -0.3333333
#> g2   0.3333333 -0.3333333  0.6666667
#> g3   0.3333333 -0.3333333 -0.3333333

Again, each column holds the linear contrast weights represented by that coefficient. It might not be obvious, but these are the deviations of each mean from the grand mean.

Again, we can confirm this:

model_parameters(model7)
#> Parameter   | Coefficient |   SE |         95% CI |  t(7) |      p
#> ------------------------------------------------------------------
#> (Intercept) |       19.99 | 0.57 | [18.65, 21.33] | 35.35 | < .001
#> G [1]       |        3.34 | 0.79 | [ 1.47,  5.22] |  4.21 | 0.004 
#> G [2]       |        0.15 | 0.71 | [-1.53,  1.83] |  0.21 | 0.840

# Predict for G = g1, g2, g3
pred_grid <- data.frame(G = c("g1", "g2", "g3"))
(predicted_means <- predict(model6, newdata = pred_grid))
#>        1        2        3 
#> 23.33333 20.14000 16.50000

# Effect contrasts
predicted_means[1] - mean(predicted_means)
#>        1 
#> 3.342222
predicted_means[2] - mean(predicted_means)
#>         2 
#> 0.1488889

I personally like the interpretability of dummy-coding better - as differences from the reference level. Thankfully, we can preserve that interpretability while also having the intercept represent the grand mean, using deviation contrasts - another type of sum-to-zero contrasts, available with datawizard::contr.deviation().

contr.deviation(levels(data$G))
#>            g2         g3
#> g1 -0.3333333 -0.3333333
#> g2  0.6666667 -0.3333333
#> g3 -0.3333333  0.6666667
model8 <- lm(Y ~ G,
             data = data,
             contrasts = list(G = contr.deviation))

model_parameters(model8)
#> Parameter   | Coefficient |   SE |          95% CI |  t(7) |      p
#> -------------------------------------------------------------------
#> (Intercept) |       19.99 | 0.57 | [ 18.65, 21.33] | 35.35 | < .001
#> G [g2]      |       -3.19 | 1.22 | [ -6.08, -0.31] | -2.62 | 0.034 
#> G [g3]      |       -6.83 | 1.52 | [-10.44, -3.23] | -4.49 | 0.003

# Predict for G = g1, g2, g3
pred_grid <- data.frame(G = c("g1", "g2", "g3"))
(predicted_means <- predict(model6, newdata = pred_grid))
#>        1        2        3 
#> 23.33333 20.14000 16.50000
mean(predicted_means)
#> [1] 19.99111

# Deviation / treatment contrasts
predicted_means[2] - predicted_means[1]
#>         2 
#> -3.193333
predicted_means[3] - predicted_means[1]
#>         3 
#> -6.833333

In Moderation Analysis

All this is great, but who cares about the intercept?

Well, this all comes into play much more substantially when moderators are involved!

When we have a moderator, the coefficient of the focal predictor represents a conditional effect. Specifically the effect of that predictor when the moderator is 0!

For two continuous predictors (\(Z\) the focal predictor and \(W\) the moderator), we can write out the regression equation as:

\[ \hat{Y} = b_0 + b_1 Z + b_2 W + b_3 W Z \] We can take out \(Z\) as a common factor:

\[ \hat{Y} = b_0 + (b_1 + b_3 W) Z + b_2 W \] In other words

\[ \text{Slope of } Z = b_1 + b_3 W \] Within this regression equation, \(b_1\) is the “intercept” - the predicted slope of \(Z\) (the conditional slope) when \(W=0\)!

If we condition on the mean of the moderator, this is often called the main effect, or the average effect.2

Let’s see this in action.

Continuous Predictor, Continuous Moderator

model9 <- lm(Y ~ X * Q,
             data = data)

model_parameters(model9)
#> Parameter   | Coefficient |    SE |            95% CI |  t(6) |     p
#> ---------------------------------------------------------------------
#> (Intercept) |        3.13 | 54.24 | [-129.58, 135.85] |  0.06 | 0.956
#> X           |        0.10 |  0.27 | [  -0.57,   0.76] |  0.35 | 0.737
#> Q           |        6.57 | 15.60 | [ -31.59,  44.73] |  0.42 | 0.688
#> X × Q       |       -0.04 |  0.08 | [  -0.24,   0.16] | -0.45 | 0.670

X’s coefficient is the conditional slope (a.k.a., the simple slope) of X when Q is 0:

# Predict for X = [0, 1] and Q = 0
pred_grid <- data.frame(X = c(0, 1), Q = 0)
(predicted_means <- predict(model9, newdata = pred_grid))
#>        1        2 
#> 3.133346 3.228691

# Conditional slope of X
predicted_means[2] - predicted_means[1]
#>          2 
#> 0.09534522

When we center Q, it is the conditional slope (a.k.a., the simple slope) of X at the mean of Q - also called the “main effect”:

model10 <- lm(Y ~ X * center(Q),
             data = data)

model_parameters(model10)
#> Parameter     | Coefficient |    SE |          95% CI |  t(6) |      p
#> ----------------------------------------------------------------------
#> (Intercept)   |       26.39 |  2.90 | [ 19.30, 33.48] |  9.11 | < .001
#> X             |       -0.03 |  0.02 | [ -0.08,  0.02] | -1.64 | 0.151 
#> center(Q)     |        6.57 | 15.60 | [-31.59, 44.73] |  0.42 | 0.688 
#> X × center(Q) |       -0.04 |  0.08 | [ -0.24,  0.16] | -0.45 | 0.670

# Predict for X = [0, 1] and Q = mean(Q)
pred_grid <- data.frame(X = c(0, 1), Q = mean(data$Q))
(predicted_means <- predict(model9, newdata = pred_grid))
#>        1        2 
#> 26.38931 26.35590

# Conditional slope of X
predicted_means[2] - predicted_means[1]
#>           2 
#> -0.03341344

Note that the coefficient of Q is not affected - it is the conditional slope of Q when X is 0. If we also wanted the main effect of Q we would need to also center X.

Continuous Predictor, Categorical Moderator

We can observe this same idea when the moderator is a factor.

# default dummy coding
model11 <- lm(Y ~ Q * A,
              data = data)

model_parameters(model11)
#> Parameter   | Coefficient |    SE |          95% CI |  t(6) |     p
#> -------------------------------------------------------------------
#> (Intercept) |       -2.40 | 11.37 | [-30.23, 25.43] | -0.21 | 0.840
#> Q           |        5.97 |  3.20 | [ -1.85, 13.80] |  1.87 | 0.111
#> A [a2]      |       14.76 | 13.58 | [-18.47, 47.99] |  1.09 | 0.319
#> Q × A [a2]  |       -3.40 |  3.81 | [-12.74,  5.93] | -0.89 | 0.406

Q’s coefficient is the conditional slope (a.k.a., the simple slope) of Q when the dummy variables (in this case just one) of A are set to 0 - when A is a1:

# Predict for Q = [0, 1] and A = a1
pred_grid <- data.frame(Q = c(0, 1), A = "a1")
(predicted_means <- predict(model11, newdata = pred_grid))
#>         1         2 
#> -2.400173  3.574452

# Conditional slope of Q
predicted_means[2] - predicted_means[1]
#>        2 
#> 5.974625

We can again use sum-to-zero contrasts to get the grand-average slope of X - the “main effect”:

model12 <- lm(Y ~ Q * A,
              data = data,
              contrasts = list(A = contr.sum))

model_parameters(model12)
#> Parameter   | Coefficient |   SE |          95% CI |  t(6) |     p
#> ------------------------------------------------------------------
#> (Intercept) |        4.98 | 6.79 | [-11.63, 21.59] |  0.73 | 0.491
#> Q           |        4.27 | 1.91 | [ -0.39,  8.94] |  2.24 | 0.066
#> A [1]       |       -7.38 | 6.79 | [-23.99,  9.23] | -1.09 | 0.319
#> Q × A [1]   |        1.70 | 1.91 | [ -2.96,  6.37] |  0.89 | 0.406

# Predict for all combinations of Q = [0, 1] and A = [a1, a2]
pred_grid <- expand.grid(Q = c(0, 1), 
                         A = c("a1", "a2"))
(predicted_means <- predict(model11, newdata = pred_grid))
#>         1         2         3         4 
#> -2.400173  3.574452 12.358596 14.929210

# Average slope of X
mean(c(predicted_means[2] - predicted_means[1], 
       predicted_means[4] - predicted_means[3]))
#> [1] 4.272619

Categorical Predictor, Continuous Moderator

This holds true when the focal predictor is a factor as well.

# default dummy coding
model13 <- lm(Y ~ A * Q,
              data = data)

model_parameters(model13)
#> Parameter   | Coefficient |    SE |          95% CI |  t(6) |     p
#> -------------------------------------------------------------------
#> (Intercept) |       -2.40 | 11.37 | [-30.23, 25.43] | -0.21 | 0.840
#> A [a2]      |       14.76 | 13.58 | [-18.47, 47.99] |  1.09 | 0.319
#> Q           |        5.97 |  3.20 | [ -1.85, 13.80] |  1.87 | 0.111
#> A [a2] × Q  |       -3.40 |  3.81 | [-12.74,  5.93] | -0.89 | 0.406

A’s coefficient (A [a2]) is the conditional difference (a.k.a., the simple effect) between the levels of A when Q is 0:

# Predict for A = [a1, a2] and Q = 0
pred_grid <- data.frame(A = c("a1", "a2"), Q = 0)
(predicted_means <- predict(model13, newdata = pred_grid))
#>         1         2 
#> -2.400173 12.358596

# Conditional treatment contrast
predicted_means[2] - predicted_means[1]
#>        2 
#> 14.75877

When we center Q, it is the conditional difference at the mean of Q - also called the “main effect”:

model14 <- lm(Y ~ A * center(Q),
              data = data)

model_parameters(model14)
#> Parameter          | Coefficient |   SE |          95% CI |  t(6) |      p
#> --------------------------------------------------------------------------
#> (Intercept)        |       18.74 | 1.15 | [ 15.92, 21.56] | 16.26 | < .001
#> A [a2]             |        2.72 | 1.49 | [ -0.93,  6.36] |  1.82 | 0.118 
#> center(Q)          |        5.97 | 3.20 | [ -1.85, 13.80] |  1.87 | 0.111 
#> A [a2] × center(Q) |       -3.40 | 3.81 | [-12.74,  5.93] | -0.89 | 0.406

# Predict for A = [a1, a2] and Q = mean(Q)
pred_grid <- data.frame(A = c("a1", "a2"), Q = mean(data$Q))
(predicted_means <- predict(model13, newdata = pred_grid))
#>        1        2 
#> 18.73805 21.45343

# Conditional treatment contrast
predicted_means[2] - predicted_means[1]
#>        2 
#> 2.715377

Note that in all cases, regardless of centering, the interaction coefficient was unchanged. This will always be true for the highest level interaction.

In ANOVA

Finally, we have Categorical-by-Categorical moderation - more often seen in ANOVA tables. As you might expect, when looking at coefficient tables, centering affects the values and meaning of the coefficient. When looking at ANOVA tables things get more complicated, as it depends on the type of sum-of-squares computed (1, 2, or 3).

You can read more about those differences in my blog post Everything You Always Wanted to Know About ANOVA.

The TL;DR is - yes, you should center your predictors / use sum-to-zero contrasts.

Summary

Centering makes coefficient tables (and ANOVA tables) more interpretable (and in Bayesian analysis can make prior elicitation easier) - a very important goal to itself, easing the interaction between the numbers and our squishy-human-flesh-brains - but that’s really all it does: it does not affect the model as a whole (doesn’t change \(R^2\) or global multicollinearity).

It is quite common to see conditional slopes interpreted as main effects even when moderators have not been properly centered. Now at least you know better!

Session Info 
sessionInfo()
#> R version 4.3.1 (2023-06-16 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 11 x64 (build 22621)
#> 
#> Matrix products: default
#> 
#> 
#> locale:
#> [1] LC_COLLATE=English_Israel.utf8  LC_CTYPE=English_Israel.utf8   
#> [3] LC_MONETARY=English_Israel.utf8 LC_NUMERIC=C                   
#> [5] LC_TIME=English_Israel.utf8    
#> 
#> time zone: Asia/Jerusalem
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] parameters_0.21.3 datawizard_0.9.1 
#> 
#> loaded via a namespace (and not attached):
#>  [1] cli_3.6.2           knitr_1.45          TH.data_1.1-2      
#>  [4] rlang_1.1.2         xfun_0.41           estimability_1.4.1 
#>  [7] xtable_1.8-4        jsonlite_1.8.8      zoo_1.8-12         
#> [10] MSBMisc_0.0.1.14    htmltools_0.5.7     rmarkdown_2.25     
#> [13] grid_4.3.1          evaluate_0.23       MASS_7.3-60        
#> [16] fastmap_1.1.1       yaml_2.3.8          mvtnorm_1.2-4      
#> [19] insight_0.19.7.4    compiler_4.3.1      multcomp_1.4-25    
#> [22] codetools_0.2-19    sandwich_3.1-0      emmeans_1.9.0      
#> [25] coda_0.19-4         htmlwidgets_1.6.4   rstudioapi_0.15.0  
#> [28] lattice_0.21-8      digest_0.6.33       splines_4.3.1      
#> [31] Matrix_1.6-4        tools_4.3.1         survival_3.5-5     
#> [34] bayestestR_0.13.1.7

References