Analyses
R
Bayesian

Figure 1 illustrates the design of Experiment 2. The stages of the experiment (i.e. training, testing no-feedback, test with feedback), are identical to that of Experiment 1. The only change is that Experiment 2 participants train, and then test, on bands in the reverse order of Experiment 1 (i.e. training on the softer bands; and testing on the harder bands).

cluster Test Phase (Counterbalanced Order)data1 Varied Training 100-300350-550600-800Test1 Test  Novel Bands  800-10001000-12001200-1400data1->Test1 data2 Constant Training 600-800data2->Test1 Test3    Final Test  Novel With Feedback  800-10001000-12001200-1400Test2  Test  Varied Training Bands  100-300350-550600-800Test1->Test2 Test2->Test3
Figure 1: Experiment 2 Design. Constant and Varied participants complete different training conditions. The training and testing bands are the reverse of Experiment 1.

Results

Testing Phase - No feedback.

In the first part of the testing phase, participants are tested from each of the velocity bands, and receive no feedback after each throw.

Deviation From Target Band

Descriptive summaries testing deviation data are provided in Table 1 and Figure 2. To model differences in accuracy between groups, we used Bayesian mixed effects regression models to the trial level data from the testing phase. The primary model predicted the absolute deviation from the target velocity band (dist) as a function of training condition (condit), target velocity band (band), and their interaction, with random intercepts and slopes for each participant (id).

\[\begin{equation} dist_{ij} = \beta_0 + \beta_1 \cdot condit_{ij} + \beta_2 \cdot band_{ij} + \beta_3 \cdot condit_{ij} \cdot band_{ij} + b_{0i} + b_{1i} \cdot band_{ij} + \epsilon_{ij} \end{equation}\]

Code
result <- test_summary_table(testE2, "dist","Deviation", mfun = list(mean = mean, median = median, sd = sd))
result$constant |> kable()
result$varied |> kable()
# make kable table with smaller font size
# result$constant |> kbl(caption="Constant Testing - Deviation",booktabs=T,escape=F) |> kable_styling(font_size = 7)
Table 1: Testing Deviation - Empirical Summary
(a) Constant Testing - Deviation
Band Band Type Mean Median Sd
100-300 Extrapolation 206 48 317
350-550 Extrapolation 194 86 268
600-800 Trained 182 112 240
800-1000 Extrapolation 200 129 233
1000-1200 Extrapolation 238 190 234
1200-1400 Extrapolation 311 254 288
(b) Varied Testing - Deviation
Band Band Type Mean Median Sd
100-300 Trained 153 25 266
350-550 Trained 138 53 233
600-800 Trained 160 120 183
800-1000 Extrapolation 261 207 257
1000-1200 Extrapolation 305 258 273
1200-1400 Extrapolation 363 314 297
Code
testE2 |>  ggplot(aes(x = vb, y = dist,fill=condit)) +
    stat_summary(geom = "bar", position=position_dodge(), fun = mean) +
    stat_summary(geom = "errorbar", position=position_dodge(.9), fun.data = mean_se, width = .4, alpha = .7) + 
  labs(x="Band", y="Deviation From Target")
Figure 2: E2. Deviations from target band during testing without feedback stage.
Code
#options(brms.backend="cmdstanr",mc.cores=4)
modelFile <- paste0(here::here("data/model_cache/"), "e2_dist_Cond_Type_RF_2")
bmtd2 <- brm(dist ~ condit * bandType + (1|bandInt) + (1|id), 
    data=testE2, file=modelFile,
    iter=5000,chains=4, control = list(adapt_delta = .94, max_treedepth = 13))
                        
#bayestestR::describe_posterior(bmtd)

mted2 <- as.data.frame(describe_posterior(bmtd2, centrality = "Mean"))[, c(1,2,4,5,6)]
colnames(mted2) <- c("Term", "Estimate","95% CrI Lower", "95% CrI Upper", "pd")

mted2 |> mutate(across(where(is.numeric), \(x) round(x, 2))) |>
  tibble::remove_rownames() |> 
  mutate(Term = stringr::str_remove(Term, "b_")) |>
  kable(booktabs=TRUE) 
Term Estimate 95% CrI Lower 95% CrI Upper pd
Intercept 190.91 125.03 259.31 1.00
conditVaried -20.58 -72.94 33.08 0.78
bandTypeExtrapolation 38.09 -6.94 83.63 0.95
conditVaried:bandTypeExtrapolation 82.00 41.89 121.31 1.00
Code
pe1td <- testE2 |>  ggplot(aes(x = vb, y = dist,fill=condit)) +
    stat_summary(geom = "bar", position=position_dodge(), fun = mean) +
    stat_summary(geom = "errorbar", position=position_dodge(.9), fun.data = mean_se, width = .4, alpha = .7) + 
  theme(legend.title=element_blank(),axis.text.x = element_text(angle = 45, hjust = 0.5, vjust = 0.5)) +
  labs(x="Band", y="Deviation From Target")

condEffects <- function(m,xvar){
  m |> ggplot(aes(x = {{xvar}}, y = .value, color = condit, fill = condit)) + 
  stat_dist_pointinterval() + 
  stat_halfeye(alpha=.1, height=.5) +
  theme(legend.title=element_blank(),axis.text.x = element_text(angle = 45, hjust = 0.5, vjust = 0.5)) 
  
}

pe1ce <- bmtd2 |> emmeans( ~condit + bandType) |>
  gather_emmeans_draws() |>
 condEffects(bandType) + labs(y="Absolute Deviation From Band", x="Band Type")
Loading required namespace: rstanarm
Code
(pe1td + pe1ce) + plot_annotation(tag_levels= 'A')
Figure 3: E2. Deviations from target band during testing without feedback stage.
Code
#contrasts(test$condit) 
# contrasts(testE2$vb)

modelName <- "e2_testDistBand_RF_5K"
e2_distBMM <- brm(dist ~ condit * bandInt + (1 + bandInt|id),
                      data=testE2,file=paste0(here::here("data/model_cache",modelName)),
                      iter=5000,chains=4)
mp2 <- GetModelStats(e2_distBMM) |> kable(escape=F,booktabs=T)
mp2
e2_distBMM |> 
  emmeans("condit",by="bandInt",at=list(bandInt=c(100,350,600,800,1000,1200)),
          epred = TRUE, re_formula = NA) |> 
  pairs() |> gather_emmeans_draws()  |> 
   summarize(median_qi(.value),pd=sum(.value>0)/n()) |>
   select(contrast,Band=bandInt,value=y,lower=ymin,upper=ymax,pd) |> 
   mutate(across(where(is.numeric), \(x) round(x, 2)),
          pd=ifelse(value<0,1-pd,pd)) |>
   kable(caption="Contrasts")
coef_details <- get_coef_details(e2_distBMM, "conditVaried")
Table 2: Experiment 2. Bayesian Mixed Model predicting absolute deviation as a function of condition (Constant vs. Varied) and Velocity Band
Term Estimate 95% CrI Lower 95% CrI Upper pd
Intercept 151.71 90.51 215.86 1.00
conditVaried -70.33 -156.87 16.66 0.94
Band 0.10 0.02 0.18 1.00
condit*Band 0.12 0.02 0.23 0.99
Contrasts
contrast Band value lower upper pd
Constant - Varied 100 57.57 -20.48 135.32 0.93
Constant - Varied 350 26.60 -30.93 83.84 0.83
Constant - Varied 600 -4.30 -46.73 38.52 0.58
Constant - Varied 800 -29.30 -69.38 11.29 0.92
Constant - Varied 1000 -54.62 -101.06 -5.32 0.98
Constant - Varied 1200 -79.63 -139.47 -15.45 0.99

The model predicting absolute deviation showed a modest tendency for the varied training group to have lower deviation compared to the constant training group (β = -70.33, 95% CI [-156.87, 16.66]),with 94% of the posterior distribution being less than 0. This suggests a potential benefit of training with variation, though the evidence is not definitive.

(SHOULD PROBABLY DO ALTERNATE ANALYSIS THAT ONLY CONSIDERS THE NOVEL EXTRAPOLATION BANDS)

Code
condEffects <- function(m){
  m |> ggplot(aes(x = bandInt, y = .value, color = condit, fill = condit)) + 
    stat_dist_pointinterval() + stat_halfeye(alpha=.2) +
    stat_lineribbon(alpha = .25, size = 1, .width = c(.95)) +
    theme(axis.text.x = element_text(angle = 45, hjust = 0.5, vjust = 0.5)) +
    ylab("Predicted X Velocity") + xlab("Band")
}


e2_distBMM |> emmeans( ~condit + bandInt, 
                       at = list(bandInt = c(100, 350, 600, 800, 1000, 1200))) |>
  gather_emmeans_draws() |>
 condEffects()+
  scale_x_continuous(breaks = c(100, 350, 600, 800, 1000, 1200), 
                     labels = levels(testE2$vb), 
                     limits = c(0, 1400)) 
Figure 4: E2. Conditioinal Effect of Training Condition and Band. Ribbon indicated 95% Credible Intervals.
Discrimination between Velocity Bands

In addition to accuracy/deviation. We also assessed the ability of participants to reliably discriminate between the velocity bands (i.e. responding differently when prompted for band 600-800 than when prompted for band 150-350). Table 3 shows descriptive statistics of this measure, and Figure 1 visualizes the full distributions of throws for each combination of condition and velocity band. To quantify discrimination, we again fit Bayesian Mixed Models as above, but this time the dependent variable was the raw x velocity generated by participants.

\[\begin{equation} vx_{ij} = \beta_0 + \beta_1 \cdot condit_{ij} + \beta_2 \cdot bandInt_{ij} + \beta_3 \cdot condit_{ij} \cdot bandInt_{ij} + b_{0i} + b_{1i} \cdot bandInt_{ij} + \epsilon_{ij} \end{equation}\]

Code
testE2 %>% group_by(id,vb,condit) |> plot_distByCondit()
Figure 5: E2 testing x velocities. Translucent bands with dash lines indicate the correct range for each velocity band.
Code
result <- test_summary_table(testE2, "vx","X Velocity" ,mfun = list(mean = mean, median = median, sd = sd))
result$constant |> kable()
result$varied |> kable()
Table 3: Testing vx - Empirical Summary
(a) Constant Testing - vx
Band Band Type Mean Median Sd
100-300 Extrapolation 457 346 354
350-550 Extrapolation 597 485 368
600-800 Trained 728 673 367
800-1000 Extrapolation 953 913 375
1000-1200 Extrapolation 1064 1012 408
1200-1400 Extrapolation 1213 1139 493
(b) Varied Testing - vx
Band Band Type Mean Median Sd
100-300 Trained 410 323 297
350-550 Trained 582 530 303
600-800 Trained 696 641 316
800-1000 Extrapolation 910 848 443
1000-1200 Extrapolation 1028 962 482
1200-1400 Extrapolation 1095 1051 510
Code
e2_vxBMM <- brm(vx ~ condit * bandInt + (1 + bandInt|id),
                        data=testE2,file=paste0(here::here("data/model_cache", "e2_testVxBand_RF_5k")),
                        iter=5000,chains=4,silent=0,
                        control=list(adapt_delta=0.94, max_treedepth=13))
mt3 <-GetModelStats(e2_vxBMM ) |> kable(escape=F,booktabs=T)
mt3
cd1 <- get_coef_details(e2_vxBMM, "conditVaried")
sc1 <- get_coef_details(e2_vxBMM, "bandInt")
intCoef1 <- get_coef_details(e2_vxBMM, "conditVaried:bandInt")
Table 4: Experiment 2. Bayesian Mixed Model Predicting Vx as a function of condition (Constant vs. Varied) and Velocity Band
Term Estimate 95% CrI Lower 95% CrI Upper pd
Intercept 362.64 274.85 450.02 1.00
conditVaried -8.56 -133.97 113.98 0.55
Band 0.71 0.58 0.84 1.00
condit*Band -0.06 -0.24 0.13 0.73

See Table 4 for the full model results.

When examining discrimination ability using the model predicting raw x-velocity, the results were less clear than those of the absolute deviation analysis. The slope on Velocity Band (β = 0.71, 95% CrI [0.58, 0.84]) indicates that participants showed good discrimination between bands overall. However, the interaction term suggested this effect was not modulated by training condition (β = -0.06, 95% CrI [-0.24, 0.13]) Thus, while varied training may provide some advantage for accuracy, both training conditions seem to have similar abilities to discriminate between velocity bands.

Code
e2_vxBMM |> emmeans( ~condit + bandInt, 
                       at = list(bandInt = c(100, 350, 600, 800, 1000, 1200))) |>
  gather_emmeans_draws() |>
  ggplot(aes(x = bandInt, y = .value, color = condit, fill = condit)) + 
  stat_dist_pointinterval() +
  stat_lineribbon(alpha = .25, size = 1, .width = c(.95)) +
  ylab("Predicted X Velocity") + xlab("Band")+
  scale_x_continuous(breaks = c(100, 350, 600, 800, 1000, 1200), 
                     labels = levels(testE2$vb), 
                     limits = c(0, 1400)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 0.5, vjust = 0.5)) 
Figure 6: Conditional effect of training condition and Band. Ribbons indicate 95% HDI.
Code
new_data_grid=map_dfr(1, ~data.frame(unique(testE2[,c("id","condit","bandInt")]))) |> 
  dplyr::arrange(id,bandInt) |> 
  mutate(condit_dummy = ifelse(condit == "Varied", 1, 0)) 

indv_coefs <- coef(e2_vxBMM)$id |> 
  as_tibble(rownames="id") |> 
  select(id, starts_with("Est")) |>
  left_join(e2Sbjs, by=join_by(id) ) |> 
  group_by(condit) |> 
  mutate(rank = rank(desc(Estimate.bandInt)),
         intErrorRank=rank((Est.Error.Intercept)),
         bandErrorRank=rank((Est.Error.bandInt)),
         nCond = n()) |> arrange(intErrorRank)

fixed_effects <- e2_vxBMM |> 
  spread_draws(`^b_.*`,regex=TRUE) |> arrange(.chain,.draw,.iteration)


random_effects <- e2_vxBMM |> 
  gather_draws(`^r_id.*$`, regex = TRUE, ndraws = 2000) |> 
  separate(.variable, into = c("effect", "id", "term"), sep = "\\[|,|\\]") |> 
  mutate(id = factor(id,levels=levels(testE2$id))) |> 
  pivot_wider(names_from = term, values_from = .value) |> arrange(id,.chain,.draw,.iteration)

cd <- left_join(random_effects, fixed_effects, by = join_by(".chain", ".iteration", ".draw")) |> 
  rename(bandInt_RF = bandInt) |>
  mutate(Slope=bandInt_RF+b_bandInt) |> group_by(id) 

cdMed <- cd |> group_by(id) |> median_qi(Slope)  |> 
  left_join(e2Sbjs, by=join_by(id)) |> group_by(condit) |>
  mutate(rankSlope=rank(Slope)) |> arrange(rankSlope)

cdMed %>% ggplot(aes(y=rankSlope, x=Slope,fill=condit,color=condit)) + 
  geom_pointrange(aes(xmin=.lower , xmax=.upper)) + 
  labs(x="Estimated Slope", y="Participant")  + facet_wrap(~condit)  

# cdMed |>  ggplot(aes(x = condit, y = Slope,fill=condit)) +
#     stat_summary(geom = "bar", position=position_dodge(), fun = mean) +
#     stat_summary(geom = "errorbar", position=position_dodge(.9), fun.data = mean_se, width = .4, alpha = .7) + 
#   geom_jitter()
#   labs(x="Band", y="Deviation From Target")
Code
new_data_grid=map_dfr(1, ~data.frame(unique(testE2[,c("id","condit","bandInt")]))) |> 
  dplyr::arrange(id,bandInt) |> 
  mutate(condit_dummy = ifelse(condit == "Varied", 1, 0)) 

indv_coefs <- as_tibble(coef(e2_vxBMM)$id, rownames="id")|> 
  select(id, starts_with("Est")) |>
  left_join(e2Sbjs, by=join_by(id) ) 


fixed_effects <- e2_vxBMM |> 
  spread_draws(`^b_.*`,regex=TRUE) |> arrange(.chain,.draw,.iteration)


random_effects <- e2_vxBMM |> 
  gather_draws(`^r_id.*$`, regex = TRUE, ndraws = 1500) |> 
  separate(.variable, into = c("effect", "id", "term"), sep = "\\[|,|\\]") |> 
  mutate(id = factor(id,levels=levels(testE2$id))) |> 
  pivot_wider(names_from = term, values_from = .value) |> arrange(id,.chain,.draw,.iteration)
Warning: Expected 3 pieces. Additional pieces discarded in 330000 rows [1, 2, 3, 4, 5,
6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].
Code
 indvDraws <- left_join(random_effects, fixed_effects, by = join_by(".chain", ".iteration", ".draw")) |> 
  rename(bandInt_RF = bandInt,RF_Intercept=Intercept) |>
  right_join(new_data_grid, by = join_by("id")) |> 
  mutate(
    Slope = bandInt_RF+b_bandInt,
    Intercept= RF_Intercept + b_Intercept,
    estimate = (b_Intercept + RF_Intercept) + (bandInt*(b_bandInt+bandInt_RF)) + (bandInt * condit_dummy) * `b_conditVaried:bandInt`,
    SlopeInt = Slope + (`b_conditVaried:bandInt`*condit_dummy)
  ) 
Warning in right_join(rename(left_join(random_effects, fixed_effects, by = join_by(".chain", : Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 1 of `x` matches multiple rows in `y`.
ℹ Row 1 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship =
  "many-to-many"` to silence this warning.
Code
  indvSlopes <- indvDraws |> group_by(id) |> median_qi(Slope,SlopeInt, Intercept,b_Intercept,b_bandInt) |>
  left_join(e2Sbjs, by=join_by(id)) |> group_by(condit) |>
    select(id,condit,Intercept,b_Intercept,starts_with("Slope"),b_bandInt, n) |>
  mutate(rankSlope=rank(Slope)) |> arrange(rankSlope)   |> ungroup()
 
  
  indvSlopes |> mutate(Condition=condit) |>  group_by(Condition) |> 
    reframe(enframe(quantile(SlopeInt, c(0.0,0.25, 0.5, 0.75,1)), "quantile", "SlopeInt")) |> 
  pivot_wider(names_from=quantile,values_from=SlopeInt,names_prefix="Q_") |>
  group_by(Condition) |>
  summarise(across(starts_with("Q"), list(mean = mean))) |> kable()
Table 5: Slope coefficients by quartile, per condition
Condition Q_0%_mean Q_25%_mean Q_50%_mean Q_75%_mean Q_100%_mean
Constant -0.2788310 0.3898592 0.697289 1.0778097 1.616840
Varied -0.2454246 0.3040574 0.679329 0.9470274 1.810293

Figure 7 visually represents the distributions of estimated slopes relating velocity band to x velocity for each participant, ordered from lowest to highest within condition. Slope values are lower overall for varied training compared to constant training. Figure Xb plots the density of these slopes for each condition. The distribution for varied training has more mass at lower values than the constant training distribution. Both figures illustrate the model’s estimate that varied training resulted in less discrimination between velocity bands, evidenced by lower slopes on average.

Code
  indvSlopes |> ggplot(aes(y=rankSlope, x=SlopeInt,fill=condit,color=condit)) + 
  geom_pointrange(aes(xmin=SlopeInt.lower , xmax=SlopeInt.upper)) + 
  labs(x="Estimated Slope", y="Participant")  + facet_wrap(~condit) +
   ggplot(indvSlopes, aes(x = SlopeInt, color = condit)) + 
  geom_density() + labs(x="Slope Coefficient",y="Density")
(a) Slope estimates by participant - ordered from lowest to highest within each condition.
Figure 7: Slope distributions between condition
Code
nSbj <- 3
indvDraws  |> indv_model_plot(indvSlopes, testE2Avg, SlopeInt,rank_variable=Slope,n_sbj=nSbj,"max")
indvDraws |> indv_model_plot(indvSlopes, testE2Avg,SlopeInt, rank_variable=Slope,n_sbj=nSbj,"min")
(a) subset with largest slopes
(b) subset with smallest slopes
Figure 8: Subset of Varied and Constant Participants with the smallest and largest estimated slope values. Red lines represent the best fitting line for each participant, gray lines are 200 random samples from the posterior distribution. Colored points and intervals at each band represent the empirical median and 95% HDI.

control for training end performance

Code
testE2 |> group_by(id,condit) |> pivot_longer(c("dist","train_end"),names_to="var",values_to="value") |> 
  ggplot(aes(x=var,y=value, fill=condit)) + stat_bar + facet_wrap(~var)

Code
testE2 |> ggplot(aes(x=train_end,y=dist,fill=condit)) + 
  stat_summary(geom = "line", position=position_dodge(), fun = mean) +
    stat_summary(geom = "errorbar", position=position_dodge(.9), fun.data = mean_se, width = .4, alpha = .7) + 
    facet_wrap(~vb) +
  labs(x="Band", y="Deviation From Target")
Warning: Width not defined
ℹ Set with `position_dodge(width = ...)`
Warning: `position_dodge()` requires non-overlapping x intervals.
`position_dodge()` requires non-overlapping x intervals.
`position_dodge()` requires non-overlapping x intervals.
`position_dodge()` requires non-overlapping x intervals.
`position_dodge()` requires non-overlapping x intervals.
`position_dodge()` requires non-overlapping x intervals.

Code
testE2 |> ggplot(aes(x=train_end,y=dist,fill=condit,col=condit)) + 
  #geom_point() +
  geom_smooth(method="loess") +
    facet_wrap(~vb) +
  labs(x="Band", y="Deviation From Target")
`geom_smooth()` using formula = 'y ~ x'

Code
# create quartiles for train_end
testE2 |> group_by(condit,vb) |> 
  mutate(train_end_q = ntile(train_end,4)) |> 
  ggplot(aes(x=train_end_q,y=dist,fill=condit)) + 
 stat_bar + 
    facet_wrap(~vb) +
  labs(x="quartile", y="Deviation From Target")

Code
bmtd3 <- brm(dist ~ condit * bandType * train_end + (1|bandInt) + (1|id), 
    data=testE2, 
    file=paste0(here::here("data/model_cache","e2_trainEnd_BT_RF2")),
    iter=1000,chains=2, control = list(adapt_delta = .92, max_treedepth = 11))
summary(bmtd3)
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: dist ~ condit * bandType * train_end + (1 | bandInt) + (1 | id) 
   Data: testE2 (Number of observations: 6709) 
  Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
         total post-warmup draws = 1000

Group-Level Effects: 
~bandInt (Number of levels: 6) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    56.29     23.16    27.44   117.13 1.01      368      453

~id (Number of levels: 110) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    96.88      7.59    83.12   112.62 1.00      238      391

Population-Level Effects: 
                                             Estimate Est.Error l-95% CI
Intercept                                       85.19     40.60     2.29
conditVaried                                    46.19     40.72   -31.33
bandTypeExtrapolation                          101.67     27.77    49.32
train_end                                        1.10      0.23     0.67
conditVaried:bandTypeExtrapolation             -11.41     31.46   -74.92
conditVaried:train_end                          -0.85      0.32    -1.51
bandTypeExtrapolation:train_end                 -0.69      0.19    -1.06
conditVaried:bandTypeExtrapolation:train_end     0.93      0.23     0.46
                                             u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept                                      163.45 1.01      222      298
conditVaried                                   127.49 1.01      237      363
bandTypeExtrapolation                          159.00 1.00      490      727
train_end                                        1.59 1.00      289      372
conditVaried:bandTypeExtrapolation              48.65 1.00      485      722
conditVaried:train_end                          -0.23 1.01      242      417
bandTypeExtrapolation:train_end                 -0.34 1.00      402      646
conditVaried:bandTypeExtrapolation:train_end     1.35 1.00      385      602

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma   241.83      2.08   237.75   245.98 1.00     1779      602

Draws were sampled using sample(hmc). 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).
Code
bayestestR::describe_posterior(bmtd3)
Summary of Posterior Distribution

Parameter                                    | Median |           95% CI |     pd |            ROPE | % in ROPE |  Rhat |    ESS
--------------------------------------------------------------------------------------------------------------------------------
(Intercept)                                  |  85.92 | [  2.29, 163.45] | 97.60% | [-27.00, 27.00] |     5.05% | 1.015 | 207.00
conditVaried                                 |  44.84 | [-31.33, 127.49] | 87.20% | [-27.00, 27.00] |    30.21% | 1.014 | 206.00
bandTypeExtrapolation                        | 100.74 | [ 49.32, 159.00] |   100% | [-27.00, 27.00] |        0% | 0.999 | 478.00
train_end                                    |   1.10 | [  0.67,   1.59] |   100% | [-27.00, 27.00] |      100% | 1.003 | 280.00
conditVaried:bandTypeExtrapolation           | -10.78 | [-74.92,  48.65] | 64.00% | [-27.00, 27.00] |    61.68% | 1.000 | 473.00
conditVaried:train_end                       |  -0.83 | [ -1.51,  -0.23] | 99.70% | [-27.00, 27.00] |      100% | 1.006 | 236.00
bandTypeExtrapolation:train_end              |  -0.70 | [ -1.06,  -0.34] |   100% | [-27.00, 27.00] |      100% | 1.000 | 385.00
conditVaried:bandTypeExtrapolation:train_end |   0.93 | [  0.46,   1.35] |   100% | [-27.00, 27.00] |      100% | 1.001 | 374.00
Code
condEffects <- function(m,xvar){
  m |> ggplot(aes(x = {{xvar}}, y = .value, color = condit, fill = condit)) + 
  stat_dist_pointinterval() + 
  stat_halfeye(alpha=.1, height=.5) +
  theme(legend.title=element_blank(),axis.text.x = element_text(angle = 45, hjust = 0.5, vjust = 0.5)) 
  
}

 bmtd3 |> emmeans( ~condit * bandType * train_end) |>
  gather_emmeans_draws() |>
 condEffects(bandType) + labs(y="Absolute Deviation From Band", x="Band Type")

Code
ce_bmtd3 <- plot(conditional_effects(bmtd3),points=FALSE,plot=FALSE)
wrap_plots(ce_bmtd3)