Analyses
R
Bayesian
Published

March 16, 2024

The major manipulation adjustment of experiment 3 is for participants to receive ordinal feedback during training, in contrast to the continuous feedback of the earlier experiments. Ordinal feedback informs participants whether a throw was too soft, too hard, or fell within the target velocity range. Experiment 3 participants were randomly assigned to both a training condition (Constant vs. Varied) and a Band Order condition (original order used in Experiment 1, or the Reverse order of Experiment 2).

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. Note that these no-feedback testing trials are identical to those of Experiment 1 and 2, as the ordinal feedback only occurs during the training phase, and final testing phase, of Experiment 3.

Deviation From Target Band

Descriptive summaries testing deviation data are provided in Table 1 and Figure 1. To model differences in accuracy between groups, we fit 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).

Code
resultOrig <- test_summary_table(testE3 |> filter(bandOrder=="Original"), "dist","Deviation", mfun = list(mean = mean, median = median, sd = sd))
resultOrig$constant |> kable() 
resultOrig$varied |> kable() 
resultRev <- test_summary_table(testE3 |> filter(bandOrder=="Reverse"), "dist","Deviation", mfun = list(mean = mean, median = median, sd = sd))
resultRev$constant |> kable() 
resultRev$varied |> kable()
Table 1: Testing Deviation - Empirical Summary
(a) Constant Testing - Deviation
Band Band Type Mean Median Sd
100-300 Extrapolation 396 325 350
350-550 Extrapolation 278 176 299
600-800 Extrapolation 173 102 215
800-1000 Trained 225 126 284
1000-1200 Extrapolation 253 192 271
1200-1400 Extrapolation 277 210 262
(b) Varied Testing - Deviation
Band Band Type Mean Median Sd
100-300 Extrapolation 383 254 385
350-550 Extrapolation 287 154 318
600-800 Extrapolation 213 140 244
800-1000 Trained 199 142 209
1000-1200 Trained 222 163 221
1200-1400 Trained 281 227 246
Band Band Type Mean Median Sd
100-300 Extrapolation 403 334 383
350-550 Extrapolation 246 149 287
600-800 Trained 155 82 209
800-1000 Extrapolation 207 151 241
1000-1200 Extrapolation 248 220 222
1200-1400 Extrapolation 322 281 264
Band Band Type Mean Median Sd
100-300 Trained 153 0 307
350-550 Trained 147 55 258
600-800 Trained 159 107 192
800-1000 Extrapolation 221 160 235
1000-1200 Extrapolation 244 185 235
1200-1400 Extrapolation 324 264 291
Code
testE3 |>  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") + facet_wrap(~bandOrder)
Figure 1: e3. Deviations from target band during testing without feedback stage.
Code
#contrasts(test$condit) 

# contrasts(testE3$vb)

modelName <- "e3_testDistBand_RF_5K"
e3_distBMM <- brm(dist ~ condit * bandOrder * bandInt + (1 + bandInt|id),
                      data=testE3,file=paste0(here::here("data/model_cache",modelName)),
                      iter=5000,chains=4)


#bayestestR::describe_posterior(e3_distBMM)
m1 <- as.data.frame(describe_posterior(e3_distBMM, centrality = "Mean"))
m2 <- fixef(e3_distBMM)
mp3 <- m1[, c(1,2,4,5,6)]
colnames(mp3) <- c("Term", "Estimate","95% CrI Lower", "95% CrI Upper", "pd")
                        
mp3 |> mutate(across(where(is.numeric), \(x) round(x, 2))) |>
  tibble::remove_rownames() |> 
  mutate(Term = stringr::str_replace_all(Term, "b_bandInt", "Band")) |>
  kable(escape=F,booktabs=T)
cd1 <- get_coef_details(e3_distBMM, "conditVaried")
sc1 <- get_coef_details(e3_distBMM, "bandInt")
intCoef1 <- get_coef_details(e3_distBMM, "conditVaried:bandInt")
Table 2: Experiment 3. 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
b_Intercept 342.85 260.18 426.01 1.00
b_conditVaried 7.38 -116.96 133.20 0.54
b_bandOrderReverse -64.99 -179.19 49.75 0.86
Band -0.13 -0.22 -0.04 1.00
b_conditVaried:bandOrderReverse -185.30 -360.16 -8.89 0.98
b_conditVaried:bandInt 0.00 -0.15 0.13 0.52
b_bandOrderReverse:bandInt 0.11 -0.01 0.24 0.96
b_conditVaried:bandOrderReverse:bandInt 0.19 -0.01 0.38 0.97

The effect of training condition in Experiment 3 showed a similar pattern to Experiment 2, with the varied group tending to have lower deviation than the constant group (β = 7.38, 95% CrI [-116.96, 133.2]), with 97% of the posterior distribution falling under 0.

(NEED TO CONTROL FOR BAND ORDER HERE)

Code
e3_distBMM |> emmeans( ~condit * bandInt * bandOrder, 
                       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 Deviation") + xlab("Velocity Band")+
  scale_x_continuous(breaks = c(100, 350, 600, 800, 1000, 1200), 
                     labels = levels(testE3$vb), 
                     limits = c(0, 1400)) +
  facet_wrap(~bandOrder) +
  theme(axis.text.x = element_text(angle = 45, hjust = 0.5, vjust = 0.5)) 
Loading required namespace: rstanarm
Figure 2: e3. 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
# testE3 |> filter(bandOrder=="Original")|> group_by(id,vb,condit) |> plot_distByCondit()
# testE3 |> filter(bandOrder=="Reverse")|> group_by(id,vb,condit) |> plot_distByCondit() +ggtitle("test")

testE3 |> group_by(id,vb,condit,bandOrder) |> plot_distByCondit() + 
  facet_wrap(bandOrder~condit,scale="free_x")
Figure 3: e3 testing x velocities. Translucent bands with dash lines indicate the correct range for each velocity band.
Code
resultOrig <- test_summary_table(testE3 |> filter(bandOrder=="Original"), "vx","X Velocity", mfun = list(mean = mean, median = median, sd = sd))
resultOrig$constant |> kable() 
resultOrig$varied |> kable() 
resultRev <- test_summary_table(testE3 |> filter(bandOrder=="Reverse"), "vx","X Velocity", mfun = list(mean = mean, median = median, sd = sd))
resultRev$constant |> kable() 
resultRev$varied |> kable()
Table 3: Testing vx - Empirical Summary
(a) Constant Testing - vx
Band Band Type Mean Median Sd
100-300 Extrapolation 680 625 370
350-550 Extrapolation 771 716 357
600-800 Extrapolation 832 786 318
800-1000 Trained 1006 916 417
1000-1200 Extrapolation 1149 1105 441
1200-1400 Extrapolation 1180 1112 443
(b) Varied Testing - vx
Band Band Type Mean Median Sd
100-300 Extrapolation 667 554 403
350-550 Extrapolation 770 688 383
600-800 Extrapolation 869 814 358
800-1000 Trained 953 928 359
1000-1200 Trained 1072 1066 388
1200-1400 Trained 1144 1093 426
Band Band Type Mean Median Sd
100-300 Extrapolation 684 634 406
350-550 Extrapolation 729 679 350
600-800 Trained 776 721 318
800-1000 Extrapolation 941 883 387
1000-1200 Extrapolation 1014 956 403
1200-1400 Extrapolation 1072 1014 442
Band Band Type Mean Median Sd
100-300 Trained 392 270 343
350-550 Trained 540 442 343
600-800 Trained 642 588 315
800-1000 Extrapolation 943 899 394
1000-1200 Extrapolation 1081 1048 415
1200-1400 Extrapolation 1185 1129 500
Code
e3_vxBMM <- brm(vx ~ condit * bandOrder * bandInt + (1 + bandInt|id),
                        data=testE3,file=paste0(here::here("data/model_cache", "e3_testVxBand_RF_5k")),
                        iter=5000,chains=4,silent=0,
                        control=list(adapt_delta=0.94, max_treedepth=13))

# mt4 <-GetModelStats(e3_vxBMM ) |> kable(escape=F,booktabs=T)
# mt4

#bayestestR::describe_posterior(e3_vxBMM)
m1 <- as.data.frame(describe_posterior(e3_vxBMM, centrality = "Mean"))
m2 <- fixef(e3_vxBMM)
mp3 <- m1[, c(1,2,4,5,6)]
colnames(mp3) <- c("Term", "Estimate","95% CrI Lower", "95% CrI Upper", "pd")
                        
mp3 |> mutate(across(where(is.numeric), \(x) round(x, 2))) |>
  tibble::remove_rownames() |> 
  mutate(Term = stringr::str_replace_all(Term, "b_bandInt", "Band")) |>
  kable(escape=F,booktabs=T)
cd1 <- get_coef_details(e3_vxBMM, "conditVaried")
sc1 <- get_coef_details(e3_vxBMM, "bandInt")
intCoef1 <- get_coef_details(e3_vxBMM, "conditVaried:bandInt")
Table 4: Experiment 3. 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
b_Intercept 601.83 504.75 699.42 1.00
b_conditVaried 12.18 -134.94 162.78 0.56
b_bandOrderReverse 13.03 -123.89 144.67 0.58
Band 0.49 0.36 0.62 1.00
b_conditVaried:bandOrderReverse -338.15 -541.44 -132.58 1.00
b_conditVaried:bandInt -0.04 -0.23 0.15 0.67
b_bandOrderReverse:bandInt -0.10 -0.27 0.08 0.86
b_conditVaried:bandOrderReverse:bandInt 0.42 0.17 0.70 1.00

See Table 4 for the full model results.

Slope estimates for experiment 3 suggest that participants were capable of distinguishing between velocity bands even when provided only ordinal feedback during training (β = 0.49, 95% CrI [0.36, 0.62]). Unlike the previous two experiments, the posterior distribution for the interaction between condition and band was consistently positive, suggestive of superior discrimination for the varied participants β = -0.04, 95% CrI [-0.23, 0.15].

Code
e3_vxBMM |> emmeans( ~condit* bandOrder* 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(testE3$vb), 
                     limits = c(0, 1400)) +
  facet_wrap(~bandOrder) +
  theme(axis.text.x = element_text(angle = 45, hjust = 0.5, vjust = 0.5))
Figure 4: Conditional effect of training condition and Band. Ribbons indicate 95% HDI.
Code
new_data_grid=map_dfr(1, ~data.frame(unique(testE3[,c("id","condit","bandInt")]))) |> 
  dplyr::arrange(id,bandInt) |> 
  mutate(condit_dummy = ifelse(condit == "Varied", 1, 0)) 

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


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


random_effects <- e3_vxBMM |> 
  gather_draws(`^r_id.*$`, regex = TRUE, ndraws = 1500) |> 
  separate(.variable, into = c("effect", "id", "term"), sep = "\\[|,|\\]") |> 
  mutate(id = factor(id,levels=levels(testE3$id))) |> 
  pivot_wider(names_from = term, values_from = .value) |> arrange(id,.chain,.draw,.iteration)
Warning: Expected 3 pieces. Additional pieces discarded in 585000 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(e3Sbjs, by=join_by(id)) |> group_by(condit) |>
    select(id,condit,bandOrder,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()
  indvSlopes |> mutate(Condition=condit) |>  group_by(Condition, bandOrder) |> 
    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(bandOrder,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.3424400 0.1739116 0.4473200 0.6960787 1.934017
Varied -0.3963599 0.0993045 0.4306245 0.7401616 1.436481
bandOrder Condition Q_0%_mean Q_25%_mean Q_50%_mean Q_75%_mean Q_100%_mean
Original Constant -0.3424400 0.2684095 0.4441671 0.6619190 1.934017
Original Varied -0.3284811 0.1605767 0.4195062 0.7098967 1.339024
Reverse Constant -0.2164930 0.1596895 0.4577490 0.7085250 1.893145
Reverse Varied -0.3963599 0.0072808 0.4322739 0.7882834 1.436481

Figure 5 shows 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") 



indvSlopes |> 
  #left_join(e3Sbjs, by=c("id","condit")) |>
  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(~bandOrder+condit) + 
  {
  ggplot(indvSlopes,
    #left_join(e3Sbjs, by=c("id","condit")), 
    aes(x = SlopeInt, color = condit)) + 
    geom_density() + 
        facet_wrap(~bandOrder) +
    labs(x="Slope Coefficient",y="Density")  }
(a) Slope estimates by participant - ordered from lowest to highest within each condition.
(b) Destiny of slope coefficients by training group
Figure 5: Slope distributions between condition
Code
nSbj <- 3
indvDraws  |> indv_model_plot(indvSlopes, testE3Avg, SlopeInt,rank_variable=Slope,n_sbj=nSbj,"max")
indvDraws |> indv_model_plot(indvSlopes, testE3Avg,SlopeInt, rank_variable=Slope,n_sbj=nSbj,"min")
(a) subset with largest slopes
(b) subset with smallest slopes
Figure 6: 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.