# Conditional Effects in BRMS for Different Interactions

I ran an exploratory study to analyze the effects and interaction effects of 3 experimental conditions. In short, these 3 categorical within-subjects factors define my data: *A: 1, 2, and 3| B: 1, 2, and 3 | C: 1, and 2. In this example, I also include a control variable for the subject/participant’s gender as well as a random intercept for the subject/participant. My response variable is from a slider scale, range: [0,1]. I follow Barr, Dale J.'s recommendation in Learning statistical models through simulation in R: An interactive textbook to include only a random intercept because of the within-subjects design and no more than a single observation per cell per subject. I’m using the package ordbetareg (Ordered Beta Regression) as proposed by Kubinec, Robert in Ordered Beta Regression: A Parsimonious, Well-Fitting Model for Continuous Data with Lower and Upper Bounds because of its demonstrated efficacy over ZOIB.

I started with the model y ~ AB + AC + Participant_Gender + (1|Participant.id)

``````allData <- within(allData, A<- relevel(A, ref = "1"))
allData <- within(allData, B <- relevel(B, ref = "2"))
allData <- within(allData, C <- relevel(C, ref = "1"))

ord_fit_mean <- ordbetareg(
formula=y ~ A*B + A*C + Participant_Gender + (1|Participant.id),
data=allData,
cores=4,
chains=4,
iter=2000,
init="0",
inits=NULL,

print(summary(ord_fit_mean))
``````

And here is the output:

`````` Family: ord_beta_reg
Links: mu = identity; phi = identity; cutzero = identity; cutone = identity
Formula: y ~ A * B + A * C + Participant_Gender + (1 | Participant.id)
Data: data (Number of observations: 1830)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000

Group-Level Effects:
~Participant.id (Number of levels: 294)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.52      0.04     0.45     0.60 1.00     1096     1913

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept                       0.20      0.09     0.03     0.37 1.00     1913     2431
A2                              0.69      0.10     0.49     0.89 1.00     1965     2626
A3                              0.84      0.10     0.64     1.03 1.00     1933     2719
B1                              0.02      0.09    -0.16     0.19 1.00     1998     2725
B3                              0.20      0.09     0.03     0.38 1.00     2277     2831
C2                              0.02      0.07    -0.12     0.16 1.00     2511     2815
Participant_GenderWomen         0.16      0.08     0.01     0.32 1.00     1788     2122
Participant_GenderNonbinary    -0.21      0.31    -0.80     0.43 1.00     2208     2798
A2:B1                          -0.10      0.12    -0.32     0.14 1.00     2404     3085
A3:B1                           0.07      0.12    -0.16     0.29 1.00     2268     2730
A2:B3                          -0.17      0.12    -0.40     0.05 1.00     2368     2740
A3:B3                          -0.15      0.12    -0.37     0.09 1.00     2288     2814
A2:C2                           0.08      0.10    -0.11     0.26 1.00     2734     2857
A3:C2                           0.01      0.10    -0.17     0.21 1.00     2853     2720

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi         6.91      0.25     6.43     7.41 1.00     3878     3010
cutzero    -3.57      0.20    -3.98    -3.20 1.00     3517     2291
cutone      1.86      0.03     1.80     1.93 1.00     3267     2296

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).
``````

Finally, I plot the conditional effects.

``````plot(conditional_effects(ord_fit_mean))
``````

*Note: I’ve omitted the Participant_Gender plot as it is irrelevant to this discussion.

I’m not 100% sure of the technical details of how these plots are derived. I’ve read that they are draws from the posterior predictive distribution. I also believe that the reported main effects coefficients are simply estimates from the base category of the factors it interacts with. I believe this also extends to the conditional effects plots for those effects. For instance, the main effects plot for A1, A2, and A3 corresponds to estimates for A where B = 2 and C = 1. This is very clear for the plot for B1, B2, B3 where that distribution of values matches the distribution of A1 for B1, B2, B3 for the A:B interaction plot.

I presented these graphics and coefficients to a colleague and they recommended removing the main effect term for B and C because they appear not to be interesting. The point I want to investigate further is the behavior associated with the A1:B3 estimates. As such, they recommended fitting the model y ~ A + A:B + A:C + Participant_Gender + (1|Participant.id). I am still hesitant about this suggestion because all paths on the internet talk about how the main effects term should almost always be present if the interaction is present. But I fit the model anyway and here are the results.

``````allData <- within(allData, A<- relevel(A, ref = "1"))
allData <- within(allData, B <- relevel(B, ref = "2"))
allData <- within(allData, C <- relevel(C, ref = "1"))

ord_fit_mean <- ordbetareg(
formula=y ~ A + A:B + A:C + Participant_Gender + (1|Participant.id),
data=allData,
cores=4,
chains=4,
iter=2000,
init="0",
inits=NULL,

print(summary(ord_fit_mean))
``````

And here is the output:

`````` Family: ord_beta_reg
Links: mu = identity; phi = identity; cutzero = identity; cutone = identity
Formula: Normal_Helpful_Score ~ A + A:B + A:C + Participant_gender + (1 | Participant.id)
Data: data (Number of observations: 1830)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000

Group-Level Effects:
~Participant.id (Number of levels: 294)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.52      0.04     0.45     0.60 1.00     1156     1852

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept                       0.21      0.09     0.04     0.37 1.00     1977     2642
A2                              0.68      0.10     0.48     0.88 1.00     2473     2175
A3                              0.83      0.10     0.64     1.02 1.00     2407     2679
Participant_GenderWomen         0.16      0.08     0.01     0.31 1.00     1373     2201
Participant_GenderNonbinary    -0.21      0.32    -0.84     0.40 1.00     1851     2213
A1:B1                           0.01      0.09    -0.16     0.19 1.00     3124     2872
A2:B1                          -0.08      0.07    -0.22     0.07 1.00     4026     3064
A3:B1                           0.08      0.08    -0.07     0.23 1.00     4092     3211
A1:B3                           0.20      0.09     0.02     0.38 1.00     3188     2894
A2:B3                           0.03      0.08    -0.12     0.19 1.00     3755     2961
A3:B3                           0.06      0.08    -0.10     0.21 1.00     4044     3055
A1:C2                           0.02      0.07    -0.12     0.16 1.00     3851     3214
A2:C2                           0.10      0.06    -0.03     0.21 1.00     4404     2804
A3:C2                           0.03      0.06    -0.09     0.16 1.00     4936     2755

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi         6.92      0.26     6.43     7.44 1.00     3837     3307
cutzero    -3.57      0.20    -3.97    -3.21 1.00     2734     2697
cutone      1.86      0.03     1.80     1.93 1.00     2591     2547

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).
``````

Again, I plot the conditional effects.

``````plot(conditional_effects(ord_fit_mean))
``````

*Note: I’ve omitted the Participant_Gender plot as it is irrelevant to this discussion.

A few follow-up questions:

1. Why are the conditional effects plots for both models identical when the models are defined differently?
2. I find the second model’s reported coefficient estimates easier to understand. For instance, the estimate of 0.20 for A1:B3 is its effect compared to A1:B2 but that doesn’t mean it is defined correctly. Is one model “more correct” than the other? Does it matter that this is an exploratory study to determine if these effects are relevant?
3. What is the best way to report or analyze the main effect of A? My understanding is any model that includes an interaction term with A will reduce the constrain the reported estimate to the base categories it interacts with. I’ve considered re-running the model with the same definition and releveling the interaction terms to get A’s effect across all categories but i was not able to find any description of others doing the same.
4. Any general observations or comments on glaring issues with my reported methods are welcomed. I am open to any recommendations for alternative methods/approaches.

Please reach out if anything is unclear or needs further clarification.