Supplementary Material
This code contains the process, models, figures, and results that were used to build and analyze the WICE model. Section 3 details the implementation of the WICE model along with the model summary. Section 4 analyzes the estimated maximum endurance time (ET) from each model using leave-one-out cross validation and Pearson correlation to compare the performance of different models, with posterior prediction plots provided as visual evidence. Section 5 examines the estimated consumed endurance (CE) at 60 seconds using Pearson correlation, and includes posterior prediction plots to further validate the model performance. Section 6 demonstrates the application of the WICE model to compare endurance differences between using hands-only versus controllers, as well as other devices (e.g., haptic devices), as discussed in the paper Discussion section. It also presents the application of the WICE model to investigate the effects of weight and location on endurance time.
The key libraries used for model generation and plotting, along with their respective versions, are listed below.
cat("tidyverse:", as.character(packageVersion("tidyverse")),
"dplyr:", as.character(packageVersion("dplyr")),
"brms:", as.character(packageVersion("brms")),
"mgcv:", as.character(packageVersion("mgcv")),
"loo:", as.character(packageVersion("loo")),
"ggplot2:", as.character(packageVersion("ggplot2")))tidyverse: 2.0.0 dplyr: 1.1.4 brms: 2.22.0 mgcv: 1.9.1 loo: 2.8.0 ggplot2: 3.5.1This section presents the implementation of WICE model and the corresponding model summary. The WICE model includes both the new calculation of T_shoulder and the individual arm mass information.
formula <- bf(
ET ~ a / (((((T_shoulder_CoM + C_theta) / T_max_NICER) * 100)^b) * c),
a ~ 1,
b ~ 1 + (1 | participant),
c ~ 1,
nl = TRUE
)
priors <- c(
prior(normal(14.86, 10), nlpar = "a", lb = 0),
prior(normal(1.83, 2), nlpar = "b", lb = 0),
prior(normal(0.000218, 0.1), nlpar = "c", lb = 0)
)
ET_WICE <- brm(
formula,
data = data,
family = shifted_lognormal(),
prior = priors,
iter = 4000,
warmup = 2000,
save_pars = save_pars(all = TRUE)
)
saveRDS(ET_WICE, file = "brm_models/ET_WICE.rds")ET_WICE <- readRDS("brm_models/ET_WICE.rds")
summary(ET_WICE) Family: shifted_lognormal
Links: mu = identity; sigma = identity; ndt = identity
Formula: ET ~ a/(((((T_shoulder_BodyM + C_theta)/T_max_NICER) * 100)^b) * c)
a ~ 1
b ~ 1 + (1 | participant)
c ~ 1
Data: data (Number of observations: 190)
Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
total post-warmup draws = 16000
Multilevel Hyperparameters:
~participant (Number of levels: 38)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(b_Intercept) 0.04 0.01 0.03 0.05 1.00 1981 3180
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
a_Intercept 2.18 1.15 0.38 4.75 1.00 5076 3848
b_Intercept 0.39 0.03 0.33 0.45 1.00 4552 7444
c_Intercept 0.13 0.07 0.02 0.29 1.00 5336 3834
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.32 0.02 0.28 0.36 1.00 7354 10182
ndt 5.35 3.19 0.34 11.85 1.00 6116 5268
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).This section contains the estimated maximum endurance time (ET) from each model (CE, NICE, NICER, WICE) using leave-one-out cross validation and Pearson correlation to compare the performance of different models, with posterior prediction plots provided as visual evidence.
This section presents the overall performance of models on estimating the maximum endurance time using leave-one-out method.
loo_ET_CE <- loo(ET_CE, newdata = data)
loo_ET_NICE <- loo(ET_NICE, newdata = data)
loo_ET_NICER <- loo(ET_NICER, newdata = data)
loo_WICE <- loo(ET_WICE, newdata = data)
# Overall Performance
loo_compare(loo_ET_CE, loo_ET_NICE, loo_ET_NICER, loo_WICE) elpd_diff se_diff
ET_WICE 0.0 0.0
ET_NICER -366.3 15.1
ET_CE -616.4 17.1
ET_NICE -863.8 17.1 This section presents the models’ performance on the female dataset for estimating maximum endurance time, using a leave-one-out approach.
data_female <- data %>% filter(gender == "Female")
loo_ET_CE <- loo(ET_CE, newdata = data_female)
loo_ET_NICE <- loo(ET_NICE, newdata = data_female)
loo_ET_NICER <- loo(ET_NICER, newdata = data_female)
loo_WICE <- loo(ET_WICE, newdata = data_female)
# Female Performance
loo_compare(loo_ET_CE, loo_ET_NICE, loo_ET_NICER, loo_WICE) elpd_diff se_diff
ET_WICE 0.0 0.0
ET_NICER -173.5 12.9
ET_CE -320.4 14.5
ET_NICE -423.2 14.9 This section presents the models’ performance on the male dataset for estimating maximum endurance time, using a leave-one-out approach.
data_male <- data %>% filter(gender == "Male")
loo_ET_CE <- loo(ET_CE, newdata = data_male)
loo_ET_NICE <- loo(ET_NICE, newdata = data_male)
loo_ET_NICER <- loo(ET_NICER, newdata = data_male)
loo_WICE <- loo(ET_WICE, newdata = data_male)
# Female Performance
loo_compare(loo_ET_CE, loo_ET_NICE, loo_ET_NICER, loo_WICE) elpd_diff se_diff
ET_WICE 0.0 0.0
ET_NICER -192.9 7.8
ET_CE -296.0 8.9
ET_NICE -440.6 8.2 This section presents the models’ performance on the bare hand dataset for estimating maximum endurance time, using a leave-one-out approach.
data_noweight <- subset(data, M_object == 0)
loo_ET_CE <- loo(ET_CE, newdata = data_noweight)
loo_ET_NICE <- loo(ET_NICE, newdata = data_noweight)
loo_ET_NICER <- loo(ET_NICER, newdata = data_noweight)
loo_WICE <- loo(ET_WICE, newdata = data_noweight)
# Female Performance
loo_compare(loo_ET_CE, loo_ET_NICE, loo_ET_NICER, loo_WICE) elpd_diff se_diff
ET_WICE 0.0 0.0
ET_NICER -76.8 5.1
ET_CE -113.6 6.6
ET_NICE -171.2 5.6 This section presents the models’ performance on the weighted object dataset for estimating maximum endurance time, using a leave-one-out approach.
data_weight <- subset(data, M_object != 0)
loo_ET_CE <- loo(ET_CE, newdata = data_weight)
loo_ET_NICE <- loo(ET_NICE, newdata = data_weight)
loo_ET_NICER <- loo(ET_NICER, newdata = data_weight)
loo_WICE <- loo(ET_WICE, newdata = data_weight)
# Female Performance
loo_compare(loo_ET_CE, loo_ET_NICE, loo_ET_NICER, loo_WICE) elpd_diff se_diff
ET_WICE 0.0 0.0
ET_NICER -289.5 14.2
ET_CE -502.8 15.6
ET_NICE -692.6 16.1 This section presents the calculation of Pearson correlation for the posterior predictions of maximum endurance time.
epred_WICE <- posterior_predict(ET_WICE, newdata = data)
data <- data %>% mutate(ET_WICE_pred = colMeans(epred_WICE))
calc_cor_stats <- function(pred_matrix, y_true) {
cor_draws <- apply(pred_matrix, 1, function(pred_row) {
cor(pred_row, y_true, method = "pearson")
})
tibble(
cor_median = median(cor_draws),
cor_lower = quantile(cor_draws, 0.025),
cor_upper = quantile(cor_draws, 0.975)
)
}
pred_cols <- c("ET_CE", "ET_NICE", "ET_NICER", "ET_WICE_pred")This section presents the Pearson correlations of models on the actual maximum endurance time.
overall_cor <- map_dfr(pred_cols, function(col) {
test <- cor.test(data[[col]], data$ET, method = "pearson")
tibble(
Model = col,
correlation = test$estimate,
lower_CI = test$conf.int[1],
upper_CI = test$conf.int[2]
)
})
print(overall_cor)# A tibble: 4 × 4
Model correlation lower_CI upper_CI
<chr> <dbl> <dbl> <dbl>
1 ET_CE -0.149 -0.285 -0.00694
2 ET_NICE -0.152 -0.288 -0.00952
3 ET_NICER 0.138 -0.00462 0.275
4 ET_WICE_pred 0.841 0.794 0.879 This section presents the Pearson correlations of the models estimated maximum endurance time on the female dataset with respect to the actual maximum endurance time.
data_female <- data %>% filter(gender == "Female")
female_cor <- data_female %>%
group_modify(~ {
map_dfr(pred_cols, function(col) {
test <- cor.test(.x[[col]], .x$ET, method = "pearson")
tibble(
Model = col,
correlation = test$estimate,
lower_CI = test$conf.int[1],
upper_CI = test$conf.int[2]
)
})
})
print(female_cor)# A tibble: 4 × 4
Model correlation lower_CI upper_CI
<chr> <dbl> <dbl> <dbl>
1 ET_CE -0.0926 -0.289 0.111
2 ET_NICE -0.0898 -0.286 0.114
3 ET_NICER -0.107 -0.302 0.0961
4 ET_WICE_pred 0.841 0.770 0.892 This section presents the Pearson correlations of the models estimated maximum endurance time on the male dataset with respect to the actual maximum endurance time.
data_male <- data %>% filter(gender == "Male")
male_cor <- data_male %>%
group_modify(~ {
map_dfr(pred_cols, function(col) {
test <- cor.test(.x[[col]], .x$ET, method = "pearson")
tibble(
Model = col,
correlation = test$estimate,
lower_CI = test$conf.int[1],
upper_CI = test$conf.int[2]
)
})
})
print(male_cor)# A tibble: 4 × 4
Model correlation lower_CI upper_CI
<chr> <dbl> <dbl> <dbl>
1 ET_CE -0.0201 -0.221 0.182
2 ET_NICE -0.0282 -0.228 0.174
3 ET_NICER -0.0835 -0.280 0.120
4 ET_WICE_pred 0.832 0.757 0.885This section presents the Pearson correlations of the models estimated maximum endurance time on the bare hand dataset with respect to the actual maximum endurance time.
data_noweight <- subset(data, M_object == 0)
noweight_cor <- data_noweight %>%
group_modify(~ {
map_dfr(pred_cols, function(col) {
test <- cor.test(.x[[col]], .x$ET, method = "pearson")
tibble(
Model = col,
correlation = test$estimate,
lower_CI = test$conf.int[1],
upper_CI = test$conf.int[2]
)
})
})
print(noweight_cor)# A tibble: 4 × 4
Model correlation lower_CI upper_CI
<chr> <dbl> <dbl> <dbl>
1 ET_CE -0.0973 -0.401 0.225
2 ET_NICE -0.121 -0.420 0.202
3 ET_NICER 0.0938 -0.228 0.398
4 ET_WICE_pred 0.778 0.613 0.878This section presents the Pearson correlations of the models estimated maximum endurance time on the weighted object dataset with respect to the actual maximum endurance time.
data_weight <- subset(data, M_object != 0)
weight_cor <- data_weight %>%
group_modify(~ {
map_dfr(pred_cols, function(col) {
test <- cor.test(.x[[col]], .x$ET, method = "pearson")
tibble(
Model = col,
correlation = test$estimate,
lower_CI = test$conf.int[1],
upper_CI = test$conf.int[2]
)
})
})
print(weight_cor)# A tibble: 4 × 4
Model correlation lower_CI upper_CI
<chr> <dbl> <dbl> <dbl>
1 ET_CE -0.211 -0.358 -0.0527
2 ET_NICE -0.206 -0.354 -0.0476
3 ET_NICER 0.233 0.0761 0.379
4 ET_WICE_pred 0.846 0.794 0.886 This section presents the visualisation of mean posterior predictions.
ET_mean
This section contains the estimated consumed endurance (CE) at 60 seconds using Pearson correlation, and includes posterior prediction plots to further validate the model performance.
data$WICE_pred <- 100 * (60 / data$ET_WICE_pred)
pred_cols <- c("CE", "NICE", "NICER", "WICE_pred")This section presents the overall performance of models in estimating consumed endurance relative to the ground truth.
overall_cor <- map_dfr(pred_cols, function(col) {
test <- cor.test(data[[col]], data$actual_CE, method = "pearson")
tibble(
Model = col,
correlation = test$estimate,
lower_CI = test$conf.int[1],
upper_CI = test$conf.int[2]
)
})
print(overall_cor)# A tibble: 4 × 4
Model correlation lower_CI upper_CI
<chr> <dbl> <dbl> <dbl>
1 CE -0.236 -0.366 -0.0973
2 NICE -0.236 -0.366 -0.0965
3 NICER 0.183 0.0419 0.317
4 WICE_pred 0.889 0.855 0.915 This section presents the performance of the models on the female dataset in estimating consumed endurance relative to the ground truth.
data_female <- data %>% filter(gender == "Female")
female_cor <- data_female %>%
group_modify(~ {
map_dfr(pred_cols, function(col) {
test <- cor.test(.x[[col]], .x$actual_CE, method = "pearson")
tibble(
Model = col,
correlation = test$estimate,
lower_CI = test$conf.int[1],
upper_CI = test$conf.int[2]
)
})
})
print(female_cor)# A tibble: 4 × 4
Model correlation lower_CI upper_CI
<chr> <dbl> <dbl> <dbl>
1 CE -0.160 -0.350 0.0433
2 NICE -0.151 -0.342 0.0518
3 NICER -0.157 -0.347 0.0462
4 WICE_pred 0.897 0.849 0.930 This section presents the performance of the models on the male dataset in estimating consumed endurance relative to the ground truth.
data_male <- data %>% filter(gender == "Male")
male_cor <- data_male %>%
group_modify(~ {
map_dfr(pred_cols, function(col) {
test <- cor.test(.x[[col]], .x$actual_CE, method = "pearson")
tibble(
Model = col,
correlation = test$estimate,
lower_CI = test$conf.int[1],
upper_CI = test$conf.int[2]
)
})
})
print(male_cor)# A tibble: 4 × 4
Model correlation lower_CI upper_CI
<chr> <dbl> <dbl> <dbl>
1 CE -0.0353 -0.235 0.167
2 NICE -0.0393 -0.239 0.164
3 NICER -0.0737 -0.271 0.130
4 WICE_pred 0.823 0.745 0.879This section presents the performance of the models on the bare hand dataset in estimating consumed endurance relative to the ground truth.
data_noweight <- subset(data, M_object == 0)
noweight_cor <- data_noweight %>%
group_modify(~ {
map_dfr(pred_cols, function(col) {
test <- cor.test(.x[[col]], .x$actual_CE, method = "pearson")
tibble(
Model = col,
correlation = test$estimate,
lower_CI = test$conf.int[1],
upper_CI = test$conf.int[2]
)
})
})
print(noweight_cor)# A tibble: 4 × 4
Model correlation lower_CI upper_CI
<chr> <dbl> <dbl> <dbl>
1 CE -0.0603 -0.369 0.260
2 NICE -0.0678 -0.375 0.253
3 NICER 0.138 -0.185 0.435
4 WICE_pred 0.760 0.585 0.868This section presents the performance of the models on the weighted object dataset in estimating consumed endurance relative to the ground truth.
data_weight <- subset(data, M_object != 0)
weight_cor <- data_weight %>%
group_modify(~ {
map_dfr(pred_cols, function(col) {
test <- cor.test(.x[[col]], .x$actual_CE, method = "pearson")
tibble(
Model = col,
correlation = test$estimate,
lower_CI = test$conf.int[1],
upper_CI = test$conf.int[2]
)
})
})
print(weight_cor)# A tibble: 4 × 4
Model correlation lower_CI upper_CI
<chr> <dbl> <dbl> <dbl>
1 CE -0.267 -0.409 -0.112
2 NICE -0.260 -0.403 -0.104
3 NICER 0.225 0.0674 0.371
4 WICE_pred 0.876 0.833 0.909data$WICE_pred <- 100 * (60 / data$ET_WICE_pred)
pred_cols <- c("CE", "NICE", "NICER", "WICE_pred")This section presents the Pearson Correlations of models in estimating consumed endurance relative to the Borg CR10 ratings.
overall_cor <- map_dfr(pred_cols, function(col) {
test <- cor.test(data[[col]], data$Borg, method = "pearson")
tibble(
Model = col,
correlation = test$estimate,
lower_CI = test$conf.int[1],
upper_CI = test$conf.int[2]
)
})
print(overall_cor)# A tibble: 4 × 4
Model correlation lower_CI upper_CI
<chr> <dbl> <dbl> <dbl>
1 CE -0.229 -0.360 -0.0900
2 NICE -0.231 -0.361 -0.0914
3 NICER 0.141 -0.000880 0.278
4 WICE_pred 0.787 0.726 0.836 This section presents the performance of the models on the female dataset in estimating consumed endurance relative to the Borg CR10 ratings.
data_female <- data %>% filter(gender == "Female")
female_cor <- data_female %>%
group_modify(~ {
map_dfr(pred_cols, function(col) {
test <- cor.test(.x[[col]], .x$Borg, method = "pearson")
tibble(
Model = col,
correlation = test$estimate,
lower_CI = test$conf.int[1],
upper_CI = test$conf.int[2]
)
})
})
print(female_cor)# A tibble: 4 × 4
Model correlation lower_CI upper_CI
<chr> <dbl> <dbl> <dbl>
1 CE -0.182 -0.370 0.0200
2 NICE -0.178 -0.366 0.0245
3 NICER -0.187 -0.375 0.0149
4 WICE_pred 0.772 0.675 0.842 This section presents the performance of the models on the male dataset in estimating consumed endurance relative to the Borg CR10 ratings.
data_male <- data %>% filter(gender == "Male")
male_cor <- data_male %>%
group_modify(~ {
map_dfr(pred_cols, function(col) {
test <- cor.test(.x[[col]], .x$Borg, method = "pearson")
tibble(
Model = col,
correlation = test$estimate,
lower_CI = test$conf.int[1],
upper_CI = test$conf.int[2]
)
})
})
print(male_cor)# A tibble: 4 × 4
Model correlation lower_CI upper_CI
<chr> <dbl> <dbl> <dbl>
1 CE -0.0625 -0.261 0.141
2 NICE -0.0660 -0.264 0.137
3 NICER -0.0922 -0.288 0.111
4 WICE_pred 0.817 0.737 0.874This section presents the performance of the models on the bare hand dataset in estimating consumed endurance relative to the Borg CR10 ratings.
data_noweight <- subset(data, M_object == 0)
noweight_cor <- data_noweight %>%
group_modify(~ {
map_dfr(pred_cols, function(col) {
test <- cor.test(.x[[col]], .x$Borg, method = "pearson")
tibble(
Model = col,
correlation = test$estimate,
lower_CI = test$conf.int[1],
upper_CI = test$conf.int[2]
)
})
})
print(noweight_cor)# A tibble: 4 × 4
Model correlation lower_CI upper_CI
<chr> <dbl> <dbl> <dbl>
1 CE -0.0975 -0.401 0.225
2 NICE -0.101 -0.404 0.221
3 NICER 0.0950 -0.227 0.399
4 WICE_pred 0.653 0.425 0.803This section presents the performance of the models on the weighted object dataset in estimating consumed endurance relative to the Borg CR10 ratings.
data_weight <- subset(data, M_object != 0)
weight_cor <- data_weight %>%
group_modify(~ {
map_dfr(pred_cols, function(col) {
test <- cor.test(.x[[col]], .x$Borg, method = "pearson")
tibble(
Model = col,
correlation = test$estimate,
lower_CI = test$conf.int[1],
upper_CI = test$conf.int[2]
)
})
})
print(weight_cor)# A tibble: 4 × 4
Model correlation lower_CI upper_CI
<chr> <dbl> <dbl> <dbl>
1 CE -0.287 -0.427 -0.134
2 NICE -0.279 -0.420 -0.125
3 NICER 0.207 0.0487 0.355
4 WICE_pred 0.755 0.677 0.816This section presents the visualisation of mean posterior predictions.
CE_mean
This section contains the application of the WICE model to compare endurance differences between using hands-only versus controllers, as well as other devices (e.g., haptic devices), as discussed in the paper Discussion section. It also presents the application of the WICE model to investigate the effects of weight and location on endurance time.
predictions <- predict(ET_WICE, newdata = newdata, allow_new_levels = TRUE)
newdata$Estimate <- round(predictions[, 1], 0)
newdata$Est.Error <- round(predictions[, 2], 0)
newdata$Q2.5 <- round(predictions[, 3], 0)
newdata$Q97.5 <- round(predictions[, 4], 0)
grouped_means <- newdata %>%
group_by(product, M_object, gender) %>%
summarise(Estimate = round(mean(Estimate, na.rm = TRUE), 0),
Error = round(mean(Est.Error, na.rm = TRUE), 0),
Q2.5 = round(mean(Q2.5, na.rm = TRUE), 0),
Q97.5 = round(mean(Q97.5, na.rm = TRUE), 0),
.groups = "drop")
print(as.data.frame(grouped_means)) product M_object gender Estimate Error Q2.5 Q97.5
1 Apple Vision Pro 0.000 Female 227 80 111 419
2 Apple Vision Pro 0.000 Male 281 99 137 519
3 Cybergrasp 0.450 Female 189 66 93 347
4 Cybergrasp 0.450 Male 244 85 119 448
5 Dexmo 0.300 Female 200 70 98 368
6 Dexmo 0.300 Male 255 89 124 470
7 ELAXO 0.540 Female 183 63 90 335
8 ELAXO 0.540 Male 237 83 116 436
9 HandMorph 0.171 Female 211 74 103 388
10 HandMorph 0.171 Male 265 93 129 489
11 HapThimble 0.100 Female 218 77 106 401
12 HapThimble 0.100 Male 272 96 132 502
13 Haptic Glove 0.640 Female 176 61 87 323
14 Haptic Glove 0.640 Male 230 80 113 423
15 HTC Vive Focus 3 0.142 Female 214 75 104 394
16 HTC Vive Focus 3 0.142 Male 268 94 130 494
17 Meta Quest 3 0.126 Female 215 75 105 396
18 Meta Quest 3 0.126 Male 270 95 131 497
19 Meta Quest 3S 0.103 Female 218 76 106 401
20 Meta Quest 3S 0.103 Male 271 95 132 500
21 Meta Quest Pro 0.164 Female 212 74 104 390
22 Meta Quest Pro 0.164 Male 266 94 130 491
23 Pico 4 0.184 Female 210 74 103 386
24 Pico 4 0.184 Male 265 93 129 487
25 PlayStation VR 0.162 Female 212 74 104 390
26 PlayStation VR 0.162 Male 266 94 130 491gender
location
weight