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,
control = list(adapt_delta = 0.8))
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,
control = list(adapt_delta = 0.8))
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:
- Why are the conditional effects plots for both models identical when the models are defined differently?
- 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?
- 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.
- 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.