compute_reliability_rocketship

Setup

library(tidyverse)
library(lubridate)
library(ggthemes)
# library(langcog)
library(mirt)
library(knitr)
library("tidyverse") # for data wrangling 
library(dplyr)
library(ggpubr)

Load data

2afc data

school_data_2afc_cleaned <- read_csv(file = here::here('data/preprocessed_data/all_preprocessed_data.csv')) %>%
  mutate(word1 = str_replace(word1, ' ', '.')) %>%
  mutate(word2 = str_replace(word2, ' ', '.')) %>%
  mutate(item_pair = paste0(word1,'_',word2)) %>%
  mutate(numAFC = 2) %>%
  mutate(sub_id = pid, numAFC = as.character(numAFC)) %>%
  filter(cohort=='school') %>%
  mutate(administration = 'initial')

We have 1606 rocketship children who completed the 2afc version.

Multiafc version

school_data_multi_afc <- read_csv(here::here('data/rocketship_data/multi-afc-with-meta.csv')) 

We have 1476 rocketship (and 200 adults) who completed the multiafc version.

completed_twice_pids <- school_data_multi_afc %>%
  filter(pid %in% unique(school_data_2afc_cleaned$pid) ) %>%
  distinct(pid)

There are 904 kids who completed both versions.

Compute correct by pid

correct_by_pid <- school_data_2afc_cleaned %>%
  group_by(pid) %>%
  summarize(prop_correct_2afc = mean(correct))  

Meta data for joining back

pid_by_age <- school_data_multi_afc %>%
  distinct(pid, age_group, age_numeric) 

Compute data structure with accuracy overall for repeat kids

correct_by_pid_by_multi_afc <- school_data_multi_afc %>%
  mutate(age_group = as.numeric(age_group)) %>%
  filter(itemType=='test') %>%
  filter(completed==TRUE) %>%
  group_by(pid, age_group) %>%
  summarize(prop_correct_overall = mean(correct)) 


completed_twice_overall <- correct_by_pid %>% # 2afc initial
  left_join(correct_by_pid_by_multi_afc, by=c('pid')) %>% # join overall multiafc
  filter(!is.na(prop_correct_overall))

Make data structure now with multiafc

correct_by_pid_by_afc <- school_data_multi_afc %>%
  filter(itemType=='test') %>%
  filter(completed==TRUE) %>%
  group_by(pid, numAFC, age_group) %>%
  summarize(prop_correct_nafc = mean(correct)) 

Make it wider to join

correct_by_pid_by_afc_wide <- correct_by_pid_by_afc %>%
  rowwise() %>%
  mutate(numAFC_string = paste0('num',numAFC)) %>%
  select(-numAFC) %>%
  pivot_wider(names_from = numAFC_string, values_from = 'prop_correct_nafc') 

Now join back with 2afc data, only get kids in both

completed_twice <- correct_by_pid %>%
  left_join(correct_by_pid_by_afc) %>%
  filter(!is.na(prop_correct_nafc))

Make a long version

completed_twice_long <- correct_by_pid %>%
  left_join(correct_by_pid_by_afc_wide) %>%
  filter(!is.na(num2)) 

Make plots of overall accuradcy correlations

ggplot(data=completed_twice_overall %>% mutate(age_group = as.numeric(age_group)), aes(x=prop_correct_2afc, y=prop_correct_overall, col=age_group)) +
  geom_point(alpha=.7) +
  geom_smooth() +
  theme_few(base_size=12) +
  ggtitle('Initial (2afc) vs. second (multiafc) visual vocab admin') +
  ggpubr::stat_cor(size=3) +
  xlab('Initial 2afc accuracy') +
  ylab('Second multiaf accuracy (overall)') +
  facet_wrap(~age_group)

Correlations between first and second administration aren’t that high.

ggplot(data=completed_twice, aes(x=prop_correct_2afc, y=prop_correct_nafc, col=as.numeric(age_group))) +
  geom_point(alpha=.6) +
  geom_smooth() +
  theme_few(base_size=16) +
  ggtitle('Initial vs. second adminstration by AFC') +
  facet_wrap(~numAFC) +
  xlab('Initial 2afc acc') +
  ylab('Second multiafc acc') +
  stat_cor() +
  theme(legend.position='none')

Estimate thetas for each kid and correlate

Read in prior trial level data

Use preprocessed data from other children/adults who have taken any version to estimate thetas for each child.

Get all 2afc data

all_2afc_data <- read_csv(file = here::here('data/preprocessed_data/all_preprocessed_data.csv')) %>%
  mutate(itemId = paste(word1, 2, sep = "_")) %>%
  filter(cohort %in% c('school','prolific')) %>%
  select(pid, itemId, correct)

All multiafc data - this includes the prolific data as well, they are marked as “adults”

all_multiafc_data <- school_data_multi_afc %>%
  mutate(itemId = paste(targetWord, numAFC, sep = "_"))  %>%
  select(pid, itemId, correct)

Merge it back together

all_data <- all_multiafc_data %>%
  full_join(all_2afc_data)

Make structure for IRT models

d_wide_all<- all_data %>%
  ungroup() %>%
  arrange(itemId) %>%
  ungroup() %>%
  pivot_wider(names_from=itemId, values_from=correct, values_fn = ~mean(.x)) %>%
  ungroup()

d_mat_all <- d_wide_all %>%
  select(-pid) %>%
  data.frame %>%
  data.matrix 

rownames(d_mat_all) <- d_wide_all$sub_id

Construct 3PL

# d_mat_all <- processed_data_resp %>% 
  # select(-...1)

start.dim <- dim(d_mat_all)[2]
mm.string = (
  'F = 1-%d,
  PRIOR = (1-%d, a1, norm, 1, 0.5),
  PRIOR = (1-%d, d, norm, 0, 2), 
  LBOUND = (1-%g, g, 0.25), 
  UBOUND = (1-%g, g, 1), 
  LBOUND = (1-%d, a1, 0.5)'
  )
mm = mirt.model(sprintf(mm.string, 
                        start.dim,
                        start.dim, 
                        start.dim,
                        start.dim, 
                        start.dim,
                        start.dim))

Run basic 3PL model with priors specified above and guessing rate per item

mod.vocab <- mirt(d_mat_all, model = mm, itemtype = '3PL')

## 
Iteration: 1, Log-Lik: -86297.621, Max-Change: 4.05003
Iteration: 2, Log-Lik: -71565.441, Max-Change: 0.42313
Iteration: 3, Log-Lik: -71055.576, Max-Change: 0.46055
Iteration: 4, Log-Lik: -70646.263, Max-Change: 0.29509
Iteration: 5, Log-Lik: -70319.275, Max-Change: 0.39677
Iteration: 6, Log-Lik: -70028.272, Max-Change: 0.34248
Iteration: 7, Log-Lik: -69804.968, Max-Change: 0.22445
Iteration: 8, Log-Lik: -69636.958, Max-Change: 0.21053
Iteration: 9, Log-Lik: -69508.048, Max-Change: 0.21068
Iteration: 10, Log-Lik: -69406.137, Max-Change: 0.17941
Iteration: 11, Log-Lik: -69327.243, Max-Change: 0.26429
Iteration: 12, Log-Lik: -69266.051, Max-Change: 0.31661
Iteration: 13, Log-Lik: -69219.569, Max-Change: 0.29294
Iteration: 14, Log-Lik: -69184.245, Max-Change: 0.21879
Iteration: 15, Log-Lik: -69157.528, Max-Change: 0.09157
Iteration: 16, Log-Lik: -69137.755, Max-Change: 0.08051
Iteration: 17, Log-Lik: -69122.970, Max-Change: 0.08404
Iteration: 18, Log-Lik: -69111.731, Max-Change: 0.06211
Iteration: 19, Log-Lik: -69103.343, Max-Change: 0.05318
Iteration: 20, Log-Lik: -69097.010, Max-Change: 0.04713
Iteration: 21, Log-Lik: -69092.231, Max-Change: 0.04067
Iteration: 22, Log-Lik: -69088.641, Max-Change: 0.04252
Iteration: 23, Log-Lik: -69085.895, Max-Change: 0.03457
Iteration: 24, Log-Lik: -69083.874, Max-Change: 0.02638
Iteration: 25, Log-Lik: -69078.630, Max-Change: 0.01786
Iteration: 26, Log-Lik: -69078.414, Max-Change: 0.05392
Iteration: 27, Log-Lik: -69078.252, Max-Change: 0.03492
Iteration: 28, Log-Lik: -69078.033, Max-Change: 0.00548
Iteration: 29, Log-Lik: -69077.969, Max-Change: 0.00835
Iteration: 30, Log-Lik: -69077.921, Max-Change: 0.00860
Iteration: 31, Log-Lik: -69077.790, Max-Change: 0.00714
Iteration: 32, Log-Lik: -69077.787, Max-Change: 0.00923
Iteration: 33, Log-Lik: -69077.785, Max-Change: 0.00089
Iteration: 34, Log-Lik: -69077.784, Max-Change: 0.00094
Iteration: 35, Log-Lik: -69077.782, Max-Change: 0.00071
Iteration: 36, Log-Lik: -69077.781, Max-Change: 0.00082
Iteration: 37, Log-Lik: -69077.778, Max-Change: 0.00019
Iteration: 38, Log-Lik: -69077.778, Max-Change: 0.00016
Iteration: 39, Log-Lik: -69077.778, Max-Change: 0.00017
Iteration: 40, Log-Lik: -69077.778, Max-Change: 0.00021
Iteration: 41, Log-Lik: -69077.778, Max-Change: 0.00013

# mirt::marginal_rxx(mod.vocab)
# Error in integrate(fn, lower = -Inf, upper = Inf, mod = mod, den = density, :
# non-finite function value

Pre-processing individual data for theta estiamtes

Create sandbox participant who has seen all trials

list.items <- colnames(d_wide_all %>% select(-pid))
df.sandbox.participant <- 
  tibble(pid = "sandbox", itemId = list.items, correct = 1)

Create trial ID structure for matrix

duplicate_ids <- school_data_multi_afc %>%
  filter(completed==TRUE) %>%
  group_by(pid) %>%
  summarize(duplicates = length(unique(runId))) %>%
  filter(duplicates>1) 

df.rocketship.roar.trials.clean <- school_data_multi_afc %>% 
  filter(pid %in% completed_twice$pid) %>%
  filter(!pid %in% duplicate_ids$pid) %>%
  filter(completed==TRUE) %>%
  filter(studyId %in%  c("school-multiAFC"), itemType == "test")  %>%
  mutate(itemId = paste(targetWord, numAFC, sep = "_")) %>%
  dplyr :: distinct(pid, runId, studyId, start_time, numAFC, targetWord, itemId, correct, rt) 

df.rocketship.2afc.roar.trials.clean <- school_data_2afc_cleaned %>% 
  mutate(itemId = paste(word1, numAFC, sep = "_")) %>%
  dplyr :: distinct(pid, runId, studyId, start_time, numAFC, word1, word2, itemId, correct, rt) 

There were 176 trials in the 2AFC version.

Make each of the data structures for getting thetas

df.multiAFC.resp <- df.rocketship.roar.trials.clean %>% 
  select(pid, itemId, correct) %>% 
  rbind(df.sandbox.participant) %>%
  unique() %>% 
  pivot_wider(names_from = itemId, values_from = correct, values_fill = NA)

df.2AFC.resp <- df.rocketship.roar.trials.clean %>% 
  filter(numAFC == 2) %>% 
  distinct(pid, itemId, correct) %>% 
  rbind(df.sandbox.participant) %>%
  unique() %>% 
  pivot_wider(names_from = itemId, values_from = correct, values_fill = NA)

df.3AFC.resp <- df.rocketship.roar.trials.clean %>% 
  filter(numAFC == 3) %>% 
  distinct(pid, itemId, correct) %>% 
  rbind(df.sandbox.participant) %>%
  unique() %>% 
  pivot_wider(names_from = itemId, values_from = correct, values_fill = NA)

df.4AFC.resp <- df.rocketship.roar.trials.clean %>% 
  filter(numAFC == 4) %>% 
  distinct(pid, itemId, correct) %>% 
  rbind(df.sandbox.participant) %>%
  unique() %>% 
  pivot_wider(names_from = itemId, values_from = correct, values_fill = NA)

Make the data structures for getting thetas from 2afc

duplicates_first_admin <- df.rocketship.2afc.roar.trials.clean %>%
  group_by(pid) %>%
  summarize(duplicates = length(unique(runId))>1) %>%
  filter(duplicates==TRUE)

# kids_enough_items <- df.rocketship.2afc.roar.trials.clean %>% 
#   filter(!pid %in% duplicates_first_admin$pid) %>%
#   distinct(pid, itemId, correct)  %>%
#   filter(itemId %in% list.items) %>%
#   group_by(pid) %>%
#   summarize(items_per_kid = length(unique(itemId))) %>%
#   filter(items_per_kid>20)

df.2AFC.first_admin <- df.rocketship.2afc.roar.trials.clean %>% 
  filter(!pid %in% duplicates_first_admin$pid) %>%
  # filter(pid %in% kids_enough_items$pid) %>%
  distinct(pid, itemId, correct) %>% 
  filter(itemId %in% list.items) %>%
  rbind(df.sandbox.participant) %>%
  unique() %>% 
  pivot_wider(names_from = itemId, values_from = correct, values_fill = NA)

Get theta score based on these reseponse pattenrs

getTheta <- function(df){
  thetas <- fscores(mod.vocab, method = "ML", response.pattern = df %>% select(-pid), max_theta = 6, min_theta = -8, theta_lim = c(-8, 6))
  return (data.frame(thetas)$F)
}

Overall and based on 234afc separately

df.thetas.students <- df.multiAFC.resp %>% 
    dplyr :: select(pid) %>% 
    add_column(theta.multiAFC = getTheta(df.multiAFC.resp), 
               theta.2AFC = getTheta(df.2AFC.resp),
               theta.3AFC = getTheta(df.3AFC.resp),
               theta.4AFC = getTheta(df.4AFC.resp)
               )

Overall and based on 234afc separately

df.thetas.students.first_admin <- df.2AFC.first_admin %>% 
    dplyr :: select(pid) %>% 
    add_column(theta.2AFC = getTheta(df.2AFC.first_admin))

Now get thetas for 2AFC kids with same model

dfs_all <- df.thetas.students.first_admin %>%
  rename(theta.2AFC.initial = theta.2AFC) %>%
  right_join(df.thetas.students, by=c('pid')) %>%
  left_join(pid_by_age) %>%
  mutate(age_group = as.numeric(age_group)) %>%
  filter(!is.na(age_group))

Students who did both multiAFC and 2AFC assessments.

ggplot(data=dfs_all, aes(x=theta.2AFC.initial, y=theta.multiAFC, col=age_group)) +
  geom_point(alpha=.6) +
  stat_cor(cor.coef.name = 'r', aes(label = ..r.label..))  +
  geom_smooth(span=20) +
  facet_grid(~age_group) +
  theme(legend.position='none')

ggplot(data=dfs_all, aes(x=age_numeric, y=theta.2AFC, col=age_group)) +
  geom_point(alpha=.6) +
  stat_cor(cor.coef.name = 'r', aes(label = ..r.label..))  +
  geom_smooth(span=20) +
  # facet_grid(~age_group) +
  theme(legend.position='none')

ggplot(data=dfs_all, aes(x=age_numeric, y=theta.multiAFC, col=age_group)) +
  geom_point(alpha=.6) +
  stat_cor(cor.coef.name = 'r', aes(label = ..r.label..))  +
  geom_smooth(span=20) +
  # facet_grid(~age_group) +
  theme(legend.position='none')

Sanity check – who are these kids at floor for multiafc?

TBH I don’t really understand why some kids are getting such a bad theta estimate despite have

low_theta <- dfs_all %>%
  filter(theta.multiAFC< -4) %>%
  mutate(age_group = as.numeric(age_group)) %>%
  left_join(correct_by_pid_by_afc, by=c('pid'))

ggplot(data=low_theta, aes(x=theta.multiAFC, y=prop_correct_nafc, col=numAFC)) +
  geom_point(alpha=.6) +
  stat_cor(cor.coef.name = 'r', aes(label = ..r.label..))  +
  geom_smooth(span=20) +
  facet_grid(~numAFC) +
  theme(legend.position='none') +
  ggtitle('low theta kids')