Sections marked (*) indicate that they contain figures and analyses either on, or directly relevant to, the poster!

Load required packages, and specify files to be read.

library(tidyverse)
library(ggthemes)
library(ggplot2)
library(ggpubr)
library(here)
library(BayesFactor)

# read in data for experiments 1 & 2
data_sum <- read_csv(here("data_tidy","bandit-S1S2-data_sum-poster.csv"))

1 Do children over-explore compared to adults?: Analyzing switching and non-maximizing choices

1.1 Switching (*)

Let’s first visualize switching behaviour - i.e., choosing an option (monster) that was different from the last option chosen.

Using the BayesFactor package, we examine if switching behaviour differ between adults and children.

switchBF = ttestBF(formula = switch ~ group, data = data_sum)
switchBF 
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 4.56554e+43 ±0%
## 
## Against denominator:
##   Null, mu1-mu2 = 0 
## ---
## Bayes factor type: BFindepSample, JZS

1.1.1 Check of assumptions

Let us the test normality assumption of this procedure.

data_sum %>%
  filter(group == "adult") %>%
  pull(switch) %>%
  shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.76418, p-value = 1.139e-13
data_sum %>%
  filter(group == "child") %>%
  pull(switch) %>%
  shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.72976, p-value = 2.307e-10
ggdensity(data_sum, 
          x = "switch",
          xlab = "switching",
          fill = "lightgray", 
          add = "mean",
          rug = TRUE,
          color = "group")

Given the way the data looks, let’s try Rubin’s Bayesian bootstrap.

library(bayesboot)

switch_c <- data_sum %>%
  filter(group == "child") %>%
  pull(switch)

switch_a <- data_sum %>%
  filter(group == "adult") %>%
  pull(switch)

# Running the Bayesian bootstrap for both datasets
b_switch_c <- bayesboot(switch_c, weighted.mean, use.weights = TRUE)
b_switch_a <- bayesboot(switch_a, weighted.mean, use.weights = TRUE)

# Calculating the posterior difference for plotting
b_switch_diff <- as.bayesboot(b_switch_c - b_switch_a)

The following plot shows us the posterior difference distribution, and associated 95% HDI, for between-groups difference in switching choices.

plot(b_switch_diff)

sum(b_switch_diff <= 0) # Was zero a value for any of the 4000 BB replications?
## [1] 0

A value of zero was not a value for any of the 4000 Bayesian bootstrap replications. So, we conclude that there is strong evidence that switching behaviour differs between adults and children.

1.2 Non-maximizing choices

Let’s visualize the non-maximizing choice proportions - another way we had operationlized exploratory behavior - for adults and children. This was the proportion of choices where participants pick an option other than the one yielding the highest reward observed so far (e.g., if Tim had chosen the 6-star monster on trial 2, but chose the 2-star monster on trial 3, Tim’s response on trial 3 would be classified as a non-maximizing choice).

We take a similar analysis approach as analyzing the proportion of switching responses.

nonMaxBF = ttestBF(formula = explore ~ group, data = data_sum)
nonMaxBF
## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 3.610737e+27 ±0%
## 
## Against denominator:
##   Null, mu1-mu2 = 0 
## ---
## Bayes factor type: BFindepSample, JZS

1.2.1 Check of assumptions

Let us test the normality assumption of the BayesFactor package.

data_sum %>%
  filter(group == "adult") %>%
  pull(explore) %>%
  shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.60512, p-value < 2.2e-16
data_sum %>%
  filter(group == "child") %>%
  pull(explore) %>%
  shapiro.test()
## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.79079, p-value = 7.003e-09
ggdensity(data_sum, 
          x = "explore",
          fill = "lightgray", 
          add = "mean",
          rug = TRUE,
          color = "group",
          xlab="non-maximizing")

Given the way the data looks, let’s try the Bayesian bootstrap.

nonMax_c <- data_sum %>%
  filter(group == "child") %>%
  pull(explore)

nonMax_a <- data_sum %>%
  filter(group == "adult") %>%
  pull(explore)

# Running the Bayesian bootstrap for both datasets
b_nonMax_c <- bayesboot(nonMax_c, weighted.mean, use.weights = TRUE)

b_nonMax_a <- bayesboot(nonMax_a, weighted.mean, use.weights = TRUE)

# Calculate the posterior difference for plotting
b_nonMax_diff <- as.bayesboot(b_nonMax_c - b_nonMax_a)

This plot shows us the posterior difference distribution, and associated 95% HDI, for between-groups difference in non-maximizing choices.

plot(b_nonMax_diff)

sum(b_nonMax_diff <= 0) # Was zero a value for any of the 4000 BB replications?
## [1] 0

We thus conclude that there is strong evidence that non-maximizing choice differs between adults and children.

Back to top

2 What leads to better change detection?: Analyzing change discovery

2.1 Discovering change in reward value in-task (*)

Let us use survival analysis to examine factors influencing discovery of the 8-star option after this change occurs during the task.

Load the required package.

library(survival)
library(survminer)

We first visualise the data using Kaplan-Meier curves.

survfit <- survfit(
  Surv(time, status) ~ group, 
  conf.int=.95,
  type=c("kaplan-meier"),
  data = data_sum)

2.1.1 Interpreting Kaplan-Meier curves (*)

The y-axis denotes the survival function \(S(t)\); the probability that a participant will ‘survive’ (i.e., that a participant will still not have chosen the 8-star option) past trial \(t\). As expected, at \(0 \le t \le 40\) , \(S(t) = 1\) since the 8-star option has not yet become available.

Here, we observe that there is a notable difference between groups for survival probability over time. The dashed vertical line indicates the median survival; here, only the median survival for the child group lies within the plotted space, as the median survival for the adult group is censored.

2.1.2 Predictors of change discovery

We use Cox proportional-hazards regression to ascertain the best predictors of change discovery. Let us first list the models corresponding to our questions of interest.

2.1.2.1 1. Are children better at detecting change than adults? (*)

discovery.cox.1 <- 
  coxph(
    Surv(time, status) ~ group, 
    data = data_sum
  )
summary(discovery.cox.1) 
## Call:
## coxph(formula = Surv(time, status) ~ group, data = data_sum)
## 
##   n= 107, number of events= 62 
##    (106 observations deleted due to missingness)
## 
##              coef exp(coef) se(coef)     z Pr(>|z|)    
## groupchild 1.4523    4.2727   0.2772 5.239 1.62e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## groupchild     4.273      0.234     2.482     7.357
## 
## Concordance= 0.668  (se = 0.034 )
## Rsquare= 0.244   (max possible= 0.993 )
## Likelihood ratio test= 29.97  on 1 df,   p=4e-08
## Wald test            = 27.44  on 1 df,   p=2e-07
## Score (logrank) test = 31.98  on 1 df,   p=2e-08

From the Wald statistic value, we can conclude that the variable group has statistically significant coefficients. The beta coefficient for groupchild indicates that children have a higher ‘risk’ of discovering change. exp(coef) - the hazard ratio - gives the effect size of the covariate. In this case, being a child increases the ‘hazard’ (i.e., discovery of 8-star option) by a factor of approximately 4 times. The upper and lower 95% confidence intervals for the hazard ratio (exp(coef)) is also given.

Taken together, we conclude that children are more likely to detect change than adults.

2.1.2.2 2. Group, switching, and ‘explore’ choices as covariates.

discovery.cox.2 <- 
  coxph(
    Surv(time, status) ~ group + explore + switch, 
    data = data_sum
  )
summary(discovery.cox.2) 
## Call:
## coxph(formula = Surv(time, status) ~ group + explore + switch, 
##     data = data_sum)
## 
##   n= 107, number of events= 62 
##    (106 observations deleted due to missingness)
## 
##               coef exp(coef) se(coef)      z Pr(>|z|)   
## groupchild -1.3960    0.2476   0.5239 -2.665  0.00770 **
## explore     2.4192   11.2370   0.9279  2.607  0.00913 **
## switch      3.7376   41.9977   1.1861  3.151  0.00163 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## groupchild    0.2476    4.03897   0.08868    0.6913
## explore      11.2370    0.08899   1.82308   69.2614
## switch       41.9977    0.02381   4.10758  429.4028
## 
## Concordance= 0.819  (se = 0.041 )
## Rsquare= 0.636   (max possible= 0.993 )
## Likelihood ratio test= 108.2  on 3 df,   p=<2e-16
## Wald test            = 70.14  on 3 df,   p=4e-15
## Score (logrank) test = 107.8  on 3 df,   p=<2e-16

In this multivariate analysis, the covariate switching remains significant, but neither group nor explore (non-maximizing choices) are significant. Switching between options is the best predictor of change discovery, and drives children’s increased likelihood to discover changes in-task.

2.1.2.2.1 Testing the proportional hazards assumption

We then test the assumptions of Cox regression for each of our models. We check the proportional hazards assumption for each covariate included in the model fit by testing for independence between corresponding sets of scaled Schoenfield residuals with time.

cox.zph(discovery.cox.1)
##              rho chisq     p
## groupchild 0.106 0.713 0.399
cox.zph(discovery.cox.2)
##                rho   chisq     p
## groupchild -0.0197 0.03315 0.856
## explore     0.0107 0.00741 0.931
## switch      0.0442 0.17299 0.677
## GLOBAL          NA 0.80812 0.848

As none of the tests for each of the covariates nor the global test of each model returned statistically significant, we assume proportional hazards.

Back to top

2.2 Performance on the posttest recall task: What predicts better recall of choice-reward associations?

2.2.1 Overall post-test response accuracy (*)

Let us first visualise the data.

We use ordinal logistic regression to examine differences in overall accuracy for post-test questions between adults and children.

library(ordinal)

data_sum$correct <- as.factor(data_sum$correct)

fit_group <- clm2(correct~group,data=data_sum)  
summary(fit_group)
## Call:
## clm2(location = correct ~ group, data = data_sum)
## 
## Location coefficients:
##            Estimate Std. Error z value Pr(>|z|) 
## groupchild  0.8357   0.2853     2.9287 0.0034041
## 
## No scale coefficients
## 
## Threshold coefficients:
##         Estimate Std. Error z value
## 0|0.2   -4.0292   0.5855    -6.8822
## 0.2|0.4 -3.3198   0.4198    -7.9076
## 0.4|0.6 -1.5323   0.2085    -7.3505
## 0.6|0.8 -0.7264   0.1748    -4.1549
## 0.8|1    0.2029   0.1672     1.2132
## 
## log-likelihood: -272.1656 
## AIC: 556.3312 
## Condition number of Hessian: 70.22688

We see here that age group is a significant predictor of posttest performance - children perform better than adults in recall of choice-reward associations.

fit2 <- clm2(correct~ group + switch ,data=data_sum)
summary(fit2)
## Call:
## clm2(location = correct ~ group + switch, data = data_sum)
## 
## Location coefficients:
##            Estimate Std. Error z value Pr(>|z|) 
## groupchild  0.2315   0.4375     0.5291 0.5967172
## switch      1.1036   0.6261     1.7625 0.0779812
## 
## No scale coefficients
## 
## Threshold coefficients:
##         Estimate Std. Error z value
## 0|0.2   -3.8500   0.5936    -6.4862
## 0.2|0.4 -3.1426   0.4309    -7.2935
## 0.4|0.6 -1.3571   0.2300    -5.9008
## 0.6|0.8 -0.5472   0.2016    -2.7140
## 0.8|1    0.3968   0.2005     1.9792
## 
## log-likelihood: -270.6316 
## AIC: 555.2632 
## Condition number of Hessian: 91.23175
fit2_nogroup <- clm2(correct~ switch ,data=data_sum)
summary(fit2_nogroup)
## Call:
## clm2(location = correct ~ switch, data = data_sum)
## 
## Location coefficients:
##        Estimate Std. Error z value Pr(>|z|)  
## switch  1.3542   0.4092     3.3095 0.00093472
## 
## No scale coefficients
## 
## Threshold coefficients:
##         Estimate Std. Error z value
## 0|0.2   -3.8358   0.5927    -6.4715
## 0.2|0.4 -3.1289   0.4298    -7.2800
## 0.4|0.6 -1.3450   0.2284    -5.8879
## 0.6|0.8 -0.5360   0.2001    -2.6782
## 0.8|1    0.4085   0.1990     2.0526
## 
## log-likelihood: -270.7728 
## AIC: 553.5457 
## Condition number of Hessian: 75.00258
anova(fit2,fit2_nogroup) 
## Likelihood ratio tests of cumulative link models
## 
## Response: correct
##                  Model Resid. df -2logLik   Test    Df  LR stat.   Pr(Chi)
## 1         switch |  |        207 541.5457                                 
## 2 group + switch |  |        206 541.2632 1 vs 2     1 0.2824349 0.5951097

When controlling for switching choices, group is no longer a significant predictor of posttest performance. Hence, this suggests that it is children’s tendency to explore (switch between options) that explains their better recall of choice-reward associations.

2.2.2 Correct identification of 8-star option

Let’s take a preliminary glimpse at the difference in proportion of participants correctly identifying the 8-star option, between adult and child groups.

data_sum %>% 
  filter(condition=="dynamic") %>% 
  group_by(group) %>% 
  summarise(mean(correct_8))
## # A tibble: 2 x 2
##   group `mean(correct_8)`
##   <chr>             <dbl>
## 1 adult             0.345
## 2 child             0.82

We find that in the dynamic task, many more children correctly identify the 8-star monster compared to adults.

Using logistic regression, we find that age group is a significant predictor of the correct identification of the 8-star option.

dataDynamic <- data_sum %>% filter(condition == "dynamic")
glm.fit <- glm(correct_8 ~ group, data = dataDynamic, family= "binomial")

summary(glm.fit)
## 
## Call:
## glm(formula = correct_8 ~ group, family = "binomial", data = dataDynamic)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8519  -0.9196   0.6300   0.6300   1.4592  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.6419     0.2763  -2.323   0.0202 *  
## groupchild    2.1582     0.4602   4.689 2.74e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 147.90  on 107  degrees of freedom
## Residual deviance: 121.86  on 106  degrees of freedom
## AIC: 125.86
## 
## Number of Fisher Scoring iterations: 4

We use the confint function to obtain confidence intervals for the coefficient estimates; CIs using profiled log-likelihood.

confint(glm.fit) 
##                 2.5 %     97.5 %
## (Intercept) -1.202014 -0.1122694
## groupchild   1.291556  3.1080330

Back to top