Setup

library(tidyverse)
library(lubridate)
library(ggthemes)
library(mirt)
library(knitr)
library(dplyr)
library(ggpubr)

Load data

Load 2afc data, which has already been cleaned

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.

There were 393 items in this initial version, including CLIP and DAIVT items.

Multiafc version

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

Get rid of participants who missed more than 1 catch trial (out of three)

catch_low <- school_data_multi_afc %>%
  group_by(pid, itemType) %>%
  summarize(pc = mean(correct), num_catch = n()) %>%
  filter(itemType=='catch') %>%
  filter(pc <.5)

school_data_multi_afc_cleaned <- school_data_multi_afc %>%
  filter(!pid %in% catch_low$pid) %>% # exclude low acc kids
  filter(!is.na(as.numeric(age_group))) %>% # exclude adults
  filter(itemType=='test') %>% # only test trials
  filter(completed==TRUE) 
  
adults_multiafc_cleaned <- school_data_multi_afc %>%
  filter(!pid %in% catch_low$pid) %>% # exclude low acc kids
  filter(!pid %in% c('wm-try')) %>% # test pid
  filter(age_group == 'Adult') %>% # include adults
  filter(itemType=='test') %>% # only test trials
  filter(completed==TRUE)

We have 1428 rocketship who completed the multiafc version.

school_data_2afc_cleaned_clip_only <- school_data_2afc_cleaned %>%
  filter(word1 %in% unique(school_data_multi_afc$targetWord))

There were 393 items in this initial version, including CLIP and DAIVT items.

There were 218 items in this initial version, including only CLIP items.

Join together pids from both administrations

completed_twice_pids <- school_data_2afc_cleaned_clip_only %>%
  filter(pid %in% unique(school_data_multi_afc_cleaned$pid) ) %>%
  distinct(pid)

There are 871 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_cleaned %>%
  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_cleaned %>%
  mutate(age_group = as.numeric(age_group)) %>%
  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_cleaned %>%
  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 accuracy 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)

Same when we look at pc correlations overall, doesn’t seem to be about some kids always being at ceiling – there’s a range

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 multiafc accuracy (overall)')

Correlations between first and second administration aren’t that high, and it doesn’t matter which AFC trials we look at

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

OK, so maybe these correlations are low because there is variability in the items, and we want to see what it’s likke when we estimate thetas for each kids.

To do so, let’s go ahead and use ALL of the data to fit a model predicting thetas for all kids together

First read back in prior trial level data and munge it

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

We need “item_ids” to be consistent across the datasets, so we have so munging to do

First get meta info on distractors for 2afc multi-afc hard trials so we can merge across the same trial tyope from both experiments

multiafc_distractors = read_csv(here::here('data/preprocessed_data/exp1_all_stimuli_april_2023.csv'))  %>%
  mutate(targetWord = Word1) %>%
  filter(trial_type == 'hard')  %>%
  select(targetWord, Word2)

Load all 2afc data in prolific/school cohort and make word1_word2 trialId

all_2afc_data <- read_csv(file = here::here('data/preprocessed_data/all_preprocessed_data.csv')) %>%
  filter(task=='test_response') %>% # only test trials should be true
  filter(completed==TRUE) %>% # completed the assessment
  filter(cohort %in% c('school','prolific')) %>%
  filter(word1 %in% multiafc_distractors$targetWord) %>% # no daivt distractors
  mutate(itemId = paste(word1, word2, sep = "_")) %>%
  mutate(numAFC = 2) %>%
  rename(targetWord = word1) %>%
  select(pid, targetWord, itemId, correct, numAFC)

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

Now munge data for identical 2afc trials from the multiafc experiment to have the same itemID format (word1_word2)

all_multiafc_data_2afc <- school_data_multi_afc_cleaned %>%
  filter(itemType=='test') %>%
  filter(numAFC==2) %>%
  mutate(administration = 'multiafc') %>%
  right_join(multiafc_distractors) %>%
  mutate(itemId = paste(targetWord, Word2, sep = "_"))  %>%
  select(pid, targetWord, itemId, correct, numAFC)

Make format for multiafc (word1_numafc) for other trials and then merge

school_data_multi_afc_munged <- school_data_multi_afc_cleaned %>%
  mutate(administration = 'multiafc') %>%
  filter(numAFC > 2) %>%
  mutate(itemId = paste(targetWord, numAFC, sep = "_")) %>%
  full_join(all_multiafc_data_2afc)

Merge multiafc and 2afc trial level data with same format

all_data <- school_data_multi_afc_munged %>%
  full_join(all_2afc_data %>% mutate(administration = '2afc'))

Identify ceiling items that we can’t model – none actually are at 100%

# ceiling <- all_data %>%
  # group_by(itemId) %>%
  # summarize(pc = mean(correct)) %>%
  # filter(pc==1)

Filter these out

# all_data <- all_data %>%
  # filter(!itemId %in% unique(ceiling$itemId))

Sanity check how much data we have by items

data_by_item <- all_data %>%
  group_by(itemId) %>%
  summarize(num_responses = length(unique(pid)))

Once bunch of items that were in both versions, other clusters from the multiafc versions only

hist(data_by_item$num_responses)

Get back this itemID list for merging after running the models

itemIdforIRT <- all_data %>%
  distinct(targetWord, numAFC, itemId)

We have FOUR trial types per word – 2AFC hard 2AFC easy 3AFC 4AFC

trials_per_word <- itemIdforIRT %>%
  group_by(targetWord) %>%
  summarize(numAFC = length(unique(itemId)))

Make structure for IRT models

This structure aggregates across administrations for each participant, e.g., if they got a trial right on one take and wrong on another the value is .5

d_wide_all<- all_data %>%
  select(itemId, correct, pid) %>%
  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$pid

Construct 3PL model where guess rate can vary

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.0)' # why was this at .5?
  )
mm = mirt.model(sprintf(mm.string, 
                        start.dim,
                        start.dim, 
                        start.dim,
                        start.dim, 
                        start.dim,
                        start.dim))

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

## 
Iteration: 1, Log-Lik: -53840.030, Max-Change: 3.88448
Iteration: 2, Log-Lik: -43211.581, Max-Change: 0.61329
Iteration: 3, Log-Lik: -42703.772, Max-Change: 0.50135
Iteration: 4, Log-Lik: -42377.940, Max-Change: 0.33983
Iteration: 5, Log-Lik: -42146.205, Max-Change: 0.39752
Iteration: 6, Log-Lik: -41982.317, Max-Change: 0.41365
Iteration: 7, Log-Lik: -41867.395, Max-Change: 0.50906
Iteration: 8, Log-Lik: -41787.636, Max-Change: 0.71529
Iteration: 9, Log-Lik: -41734.913, Max-Change: 0.54049
Iteration: 10, Log-Lik: -41701.562, Max-Change: 0.38914
Iteration: 11, Log-Lik: -41680.912, Max-Change: 0.33943
Iteration: 12, Log-Lik: -41667.931, Max-Change: 0.40546
Iteration: 13, Log-Lik: -41659.446, Max-Change: 0.33233
Iteration: 14, Log-Lik: -41654.178, Max-Change: 0.10536
Iteration: 15, Log-Lik: -41650.932, Max-Change: 0.09902
Iteration: 16, Log-Lik: -41645.368, Max-Change: 0.13402
Iteration: 17, Log-Lik: -41645.279, Max-Change: 0.86698
Iteration: 18, Log-Lik: -41645.211, Max-Change: 0.02201
Iteration: 19, Log-Lik: -41645.211, Max-Change: 0.02144
Iteration: 20, Log-Lik: -41645.210, Max-Change: 0.00485
Iteration: 21, Log-Lik: -41645.210, Max-Change: 0.00329
Iteration: 22, Log-Lik: -41645.210, Max-Change: 0.00014
Iteration: 23, Log-Lik: -41645.210, Max-Change: 0.00024
Iteration: 24, Log-Lik: -41645.210, Max-Change: 0.00019

Construct a 2PL model where guess rate is specificed

Make a list for specifying guess rate in IRT instead of letting it be a paramters?

itemIdforIRTGuess <- itemIdforIRT %>%
  arrange(itemId) %>%
  mutate(guess_rate = case_when( numAFC==2 ~0.5 ,
                                 numAFC==3 ~0.33 ,
                                 numAFC==4 ~0.25 )) 

list.guessing = itemIdforIRTGuess$guess_rate

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

mod.vocab.2pl <- mirt(d_mat_all, model = mm, itemtype = '2PL', guess = list.guessing)

## 
Iteration: 1, Log-Lik: -52019.967, Max-Change: 5.42488
Iteration: 2, Log-Lik: -42361.606, Max-Change: 0.50588
Iteration: 3, Log-Lik: -42182.162, Max-Change: 0.22023
Iteration: 4, Log-Lik: -42097.826, Max-Change: 0.16255
Iteration: 5, Log-Lik: -42044.808, Max-Change: 0.13601
Iteration: 6, Log-Lik: -42010.819, Max-Change: 0.11239
Iteration: 7, Log-Lik: -41989.080, Max-Change: 0.09160
Iteration: 8, Log-Lik: -41975.227, Max-Change: 0.07404
Iteration: 9, Log-Lik: -41966.426, Max-Change: 0.05955
Iteration: 10, Log-Lik: -41960.844, Max-Change: 0.04770
Iteration: 11, Log-Lik: -41957.309, Max-Change: 0.03811
Iteration: 12, Log-Lik: -41955.074, Max-Change: 0.03040
Iteration: 13, Log-Lik: -41951.406, Max-Change: 0.00685
Iteration: 14, Log-Lik: -41951.345, Max-Change: 0.00506
Iteration: 15, Log-Lik: -41951.307, Max-Change: 0.00397
Iteration: 16, Log-Lik: -41951.246, Max-Change: 0.00090
Iteration: 17, Log-Lik: -41951.244, Max-Change: 0.00070
Iteration: 18, Log-Lik: -41951.244, Max-Change: 0.00056
Iteration: 19, Log-Lik: -41951.243, Max-Change: 0.00026
Iteration: 20, Log-Lik: -41951.243, Max-Change: 0.00005

Is 2pl vs 3pl better?

anova(mod.vocab.2pl, mod.vocab.3pl)

##                    AIC    SABIC       HQ      BIC    logLik   logPost  df
## mod.vocab.2pl 83356.76 85511.79 85149.55 88256.82 -40814.38 -41951.24  NA
## mod.vocab.3pl 83802.70 87035.24 86491.89 91152.79 -40605.35 -41645.21 432

Get some coefficients for these models – chosing 2pl since lower AIC/BIC

co <- coef(mod.vocab.2pl,simplify=TRUE, IRTpars = FALSE) # Get coeeficients
coefs <- tibble::rownames_to_column(as.data.frame(co$items),'itemId')

Very good marginal_rxx?

mirt::marginal_rxx(mod.vocab.2pl)

## [1] 0.9694614

item_fits = itemfit(mod.vocab.2pl, fit_stats='infit')

It looks like we have really low infit statistics – I’m not sure what is going on here.

hist(item_fits$infit)

Outfits mostly look ok?

hist(item_fits$outfit)

plot(mod.vocab.2pl)

Trim off items with low slopes & rerunning 2pl

coefs_trimmed <- coefs %>%
  rename(slope = a1) %>%
  filter(slope>.8) %>%
  left_join(itemIdforIRT, by=c('itemId'))  %>%
  group_by(targetWord) %>%
  mutate(average_slope = mean(slope), both_items = length(unique(itemId)))

We got rid of 71 of the items with lowest slopes.

trimmed_items <- coefs %>%
  rename(slope = a1) %>%
  filter(slope<=.8) %>%
  left_join(itemIdforIRT, by=c('itemId'))  %>%
  arrange(targetWord)

Looks like these items come from 45 unqique words, which is nice.

kable(trimmed_items %>% arrange(slope))

itemId	slope	d	g	u	targetWord	numAFC
gondola_trolley	0.1067988	-0.6478280	0.50	1	gondola	2
candlestick_4	0.1701950	-1.3805171	0.25	1	candlestick	4
gondola_3	0.2209317	-0.3495146	0.33	1	gondola	3
gondola_4	0.3728742	-0.1408634	0.25	1	gondola	4
tuxedo_4	0.4283888	-0.8250671	0.25	1	tuxedo	4
freezer_4	0.4303087	-0.3696897	0.25	1	freezer	4
tuxedo_3	0.4536586	-1.3099747	0.33	1	tuxedo	3
taillight_3	0.4604713	0.4208304	0.33	1	taillight	3
koala_bee	0.4647243	3.2771571	0.50	1	koala	2
grate_3	0.4695413	-1.8370996	0.33	1	grate	3
footbath_trough	0.4781102	3.0808308	0.50	1	footbath	2
sauerkraut_dashboard	0.4955024	0.8802140	0.50	1	sauerkraut	2
turtle_frog	0.5140553	1.5912830	0.50	1	turtle	2
teapot_3	0.5151518	0.7788158	0.33	1	teapot	3
sorbet_tamale	0.5328676	0.4471138	0.50	1	sorbet	2
pajamas_pants	0.5383338	-0.6455522	0.50	1	pajamas	2
fox_3	0.5472977	2.1476012	0.33	1	fox	3
taillight_bumper	0.5549854	-0.0914494	0.50	1	taillight	2
spatula_paddle	0.5632006	0.2050834	0.50	1	spatula	2
freezer_3	0.5662422	-0.6941667	0.33	1	freezer	3
kimono_turntable	0.5736853	-0.1410427	0.50	1	kimono	2
cymbal_kaleidoscope	0.5750447	-1.0114297	0.50	1	cymbal	2
oil_3	0.5805746	1.1164655	0.33	1	oil	3
koala_3	0.5813491	1.8556671	0.33	1	koala	3
fox_4	0.5910830	2.4237912	0.25	1	fox	4
headdress_3	0.6070051	0.0832789	0.33	1	headdress	3
teapot_4	0.6118214	1.0306691	0.25	1	teapot	4
lollipop_doorbell	0.6125777	2.8564999	0.50	1	lollipop	2
pajamas_cage	0.6145629	2.1364722	0.50	1	pajamas	2
stretcher_thermostat	0.6235046	1.5316646	0.50	1	stretcher	2
footbath_vest	0.6444995	3.1257709	0.50	1	footbath	2
tapestry_3	0.6467081	0.5350411	0.33	1	tapestry	3
sauerkraut_4	0.6477774	-1.4705269	0.25	1	sauerkraut	4
artichoke_sunroof	0.6485020	1.1577057	0.50	1	artichoke	2
stew_4	0.6485624	-0.2816407	0.25	1	stew	4
omelet_pancake	0.6518754	2.5997755	0.50	1	omelet	2
teabag_4	0.6750026	1.2623273	0.25	1	teabag	4
pump_coat	0.6760156	2.1812969	0.50	1	pump	2
grate_4	0.6893575	-1.2594698	0.25	1	grate	4
fox_wolf	0.7038353	1.8195733	0.50	1	fox	2
headdress_poster	0.7076194	0.9659507	0.50	1	headdress	2
mandolin_gourd	0.7095347	0.2027028	0.50	1	mandolin	2
headdress_crown	0.7148810	0.5808639	0.50	1	headdress	2
elbow_4	0.7149578	1.6821833	0.25	1	elbow	4
horn_bone	0.7161944	0.3539825	0.50	1	horn	2
carrot_vegetable	0.7167608	3.4747239	0.50	1	carrot	2
candlestick_3	0.7197428	-2.0724587	0.33	1	candlestick	3
cornbread_wreath	0.7234291	3.3167238	0.50	1	cornbread	2
map_marker	0.7241045	2.7809587	0.50	1	map	2
gondola_clipboard	0.7254850	1.6753324	0.50	1	gondola	2
telescope_3	0.7269019	2.8295066	0.33	1	telescope	3
mulch_compost	0.7364623	-5.3051334	0.50	1	mulch	2
mandolin_4	0.7372156	-0.1763239	0.25	1	mandolin	4
marshmallow_snowball	0.7401504	3.0968672	0.50	1	marshmallow	2
grate_crate	0.7452091	-4.2467257	0.50	1	grate	2
turkey_swan	0.7480293	3.1128110	0.50	1	turkey	2
hoe_peg	0.7481080	1.3235930	0.50	1	hoe	2
koala_panda	0.7487401	2.1059978	0.50	1	koala	2
spatula_4	0.7526373	0.5589374	0.25	1	spatula	4
sloth_3	0.7562491	1.4303120	0.33	1	sloth	3
cake_dessert	0.7599294	2.9353615	0.50	1	cake	2
sunflower_rim	0.7607375	3.4876570	0.50	1	sunflower	2
sandbag_3	0.7636237	0.8599957	0.33	1	sandbag	3
cymbal_4	0.7717604	-0.2918655	0.25	1	cymbal	4
sauerkraut_3	0.7728945	-1.0907018	0.33	1	sauerkraut	3
candlestick_candle	0.7785190	-4.6791840	0.50	1	candlestick	2
waterwheel_4	0.7794844	1.0899263	0.25	1	waterwheel	4
headdress_4	0.7914213	0.6841997	0.25	1	headdress	4
cymbal_mallet	0.7915190	1.3359524	0.50	1	cymbal	2
caramel_amber	0.7961236	1.0165539	0.50	1	caramel	2
saddle_figurine	0.7971800	0.6279682	0.50	1	saddle	2

Construct a 2PL with trimmed items

Construct trimmed dataframe for refitting 2PL model

d_wide_trimmed<- all_data %>%
  filter(itemId %in% unique(coefs_trimmed$itemId)) %>%
  select(itemId, correct, pid) %>%
  ungroup() %>%
  arrange(itemId) %>%
  ungroup() %>%
  pivot_wider(names_from=itemId, values_from=correct, values_fn = ~mean(.x)) %>%
  ungroup()

d_mat_all_trimmed <- d_wide_trimmed %>%
  select(-pid) %>%
  data.frame %>%
  data.matrix 

rownames(d_mat_all_trimmed) <- d_wide_trimmed$pid

Set model parameters

start.dim <- dim(d_mat_all_trimmed)[2]
mm.trimmed.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.0)' # why was this at .5?
  )
mm.trimmed = mirt.model(sprintf(mm.trimmed.string, 
                        start.dim,
                        start.dim, 
                        start.dim,
                        start.dim, 
                        start.dim,
                        start.dim))

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

First get guessing list

itemIdforIRTTrimmed <- all_data %>%
  filter(itemId %in% unique(coefs_trimmed$itemId)) %>%
  distinct(itemId, targetWord,numAFC) %>%
  arrange(itemId) %>%
  mutate(guess_rate = case_when( numAFC==2 ~0.5 ,
                                 numAFC==3 ~0.33 ,
                                 numAFC==4 ~0.25 )) 

list.guessing.trimmed = itemIdforIRTTrimmed$guess_rate

mod.vocab.2pl.trimmed <- mirt(d_mat_all_trimmed, model = mm.trimmed, itemtype = '2PL', guess = list.guessing.trimmed)

## 
Iteration: 1, Log-Lik: -42078.310, Max-Change: 4.73729
Iteration: 2, Log-Lik: -34345.330, Max-Change: 0.50435
Iteration: 3, Log-Lik: -34166.837, Max-Change: 0.22793
Iteration: 4, Log-Lik: -34077.001, Max-Change: 0.17967
Iteration: 5, Log-Lik: -34021.571, Max-Change: 0.14525
Iteration: 6, Log-Lik: -33987.152, Max-Change: 0.11620
Iteration: 7, Log-Lik: -33965.848, Max-Change: 0.09234
Iteration: 8, Log-Lik: -33952.698, Max-Change: 0.07318
Iteration: 9, Log-Lik: -33944.597, Max-Change: 0.05774
Iteration: 10, Log-Lik: -33939.614, Max-Change: 0.04550
Iteration: 11, Log-Lik: -33936.552, Max-Change: 0.03575
Iteration: 12, Log-Lik: -33934.673, Max-Change: 0.02808
Iteration: 13, Log-Lik: -33931.837, Max-Change: 0.00652
Iteration: 14, Log-Lik: -33931.782, Max-Change: 0.00482
Iteration: 15, Log-Lik: -33931.748, Max-Change: 0.00375
Iteration: 16, Log-Lik: -33931.698, Max-Change: 0.00083
Iteration: 17, Log-Lik: -33931.697, Max-Change: 0.00065
Iteration: 18, Log-Lik: -33931.697, Max-Change: 0.00052
Iteration: 19, Log-Lik: -33931.696, Max-Change: 0.00032
Iteration: 20, Log-Lik: -33931.696, Max-Change: 0.00012
Iteration: 21, Log-Lik: -33931.696, Max-Change: 0.00006

co_trimmed_model <- coef(mod.vocab.2pl.trimmed,simplify=TRUE, IRTpars = FALSE) # Get coeeficients
coefs_trimmed_model <- tibble::rownames_to_column(as.data.frame(co_trimmed_model$items),'itemId')

# coefs_trimmed_model %>%
  # arrange(a1) %>%
  # kable()

 hard_but_good <- coefs_trimmed_model %>%
  filter(d>quantile(d,.75)) %>%
  arrange(a1) %>%
  kable()

marginal_rxx(mod.vocab.2pl)

## [1] 0.9694614

Do our items load onto one factor mostly?

Read in aoa and other meta

item_meta <- read_csv(file = here::here('data/preprocessed_data/exp1_all_trials2023-04-11.csv'))

word1_aoa <- item_meta %>%
  distinct(Word1, AoA_Est_Word1)  %>%
  rename(targetWord = Word1)

item_check = as.data.frame(summary(mod.vocab.2pl.trimmed))

##                            F    h2
## acorn_3                0.595 0.353
## acorn_4                0.613 0.375
## acorn_coconut          0.623 0.388
## acorn_key              0.534 0.286
## aloe_3                 0.604 0.365
## aloe_4                 0.709 0.503
## aloe_bracket           0.641 0.410
## aloe_cactus            0.633 0.400
## antenna_3              0.541 0.293
## antenna_4              0.594 0.353
## antenna_pinecone       0.524 0.274
## antenna_stem           0.583 0.340
## artichoke_3            0.522 0.273
## artichoke_4            0.627 0.393
## artichoke_leek         0.581 0.337
## bamboo_3               0.607 0.369
## bamboo_4               0.622 0.387
## bamboo_lumber          0.583 0.340
## bamboo_yacht           0.440 0.193
## barrel_3               0.604 0.365
## barrel_4               0.544 0.296
## barrel_bucket          0.578 0.334
## barrel_shelf           0.517 0.267
## blender_3              0.621 0.385
## blender_4              0.594 0.353
## blender_blind          0.457 0.209
## blender_mailbox        0.530 0.281
## blower_3               0.659 0.434
## blower_4               0.439 0.193
## blower_buggy           0.484 0.234
## blower_burner          0.547 0.300
## bobsled_3              0.637 0.406
## bobsled_4              0.411 0.169
## bobsled_dip            0.503 0.253
## bobsled_sidecar        0.469 0.220
## bouquet_3              0.693 0.480
## bouquet_4              0.723 0.523
## bouquet_bandanna       0.707 0.499
## bouquet_centerpiece    0.731 0.534
## buffet_3               0.550 0.303
## buffet_4               0.649 0.422
## buffet_counter         0.635 0.403
## buffet_crib            0.525 0.276
## bulldozer_3            0.507 0.257
## bulldozer_4            0.631 0.398
## bulldozer_orange       0.547 0.299
## bulldozer_tractor      0.634 0.403
## cake_3                 0.587 0.344
## cake_4                 0.597 0.356
## cake_refrigerator      0.454 0.206
## candlestick_egg        0.531 0.281
## caramel_3              0.500 0.250
## caramel_4              0.640 0.410
## caramel_olive          0.604 0.365
## carousel_3             0.512 0.262
## carousel_4             0.650 0.422
## carousel_carriage      0.548 0.300
## carousel_dish          0.552 0.305
## carrot_3               0.531 0.282
## carrot_4               0.595 0.354
## carrot_phone           0.581 0.337
## cassette_3             0.561 0.315
## cassette_4             0.591 0.349
## cassette_book          0.522 0.272
## cassette_tape          0.522 0.272
## cheese_3               0.668 0.447
## cheese_4               0.709 0.502
## cheese_butter          0.735 0.540
## cheese_mud             0.467 0.218
## cloak_3                0.634 0.403
## cloak_4                0.714 0.509
## cloak_cinnamon         0.594 0.353
## cloak_hood             0.605 0.365
## clothespin_3           0.498 0.248
## clothespin_4           0.459 0.211
## clothespin_toothpick   0.586 0.344
## clothespin_volleyball  0.582 0.339
## coaster_3              0.616 0.380
## coaster_4              0.587 0.344
## coaster_card           0.496 0.246
## coaster_painting       0.499 0.249
## cork_3                 0.670 0.449
## cork_4                 0.688 0.473
## cork_jersey            0.593 0.352
## cork_skateboard        0.622 0.387
## cornbread_3            0.574 0.330
## cornbread_4            0.531 0.282
## cornbread_corn         0.620 0.384
## corset_3               0.456 0.208
## corset_4               0.417 0.174
## corset_harness         0.700 0.491
## corset_mantle          0.494 0.244
## cymbal_3               0.453 0.205
## dumpling_3             0.584 0.341
## dumpling_4             0.544 0.296
## dumpling_meatball      0.472 0.223
## dumpling_rocket        0.647 0.418
## elbow_3                0.528 0.278
## elbow_arm              0.460 0.212
## elbow_bed              0.517 0.267
## fan_3                  0.565 0.319
## fan_4                  0.579 0.335
## fan_album              0.616 0.379
## fan_boulder            0.536 0.287
## flan_3                 0.549 0.302
## flan_4                 0.501 0.251
## flan_amplifier         0.576 0.332
## flan_fuse              0.580 0.337
## foam_3                 0.589 0.346
## foam_4                 0.640 0.410
## foam_float             0.725 0.526
## foam_quilt             0.543 0.295
## footbath_3             0.532 0.283
## footbath_4             0.644 0.415
## fox_reindeer           0.533 0.284
## freezer_cooler         0.499 0.249
## freezer_nest           0.499 0.249
## fruitcake_3            0.493 0.243
## fruitcake_4            0.663 0.439
## fruitcake_dough        0.431 0.186
## fruitcake_fudge        0.506 0.256
## grate_reel             0.515 0.265
## gutter_3               0.493 0.243
## gutter_4               0.571 0.326
## gutter_almond          0.581 0.337
## gutter_filter          0.551 0.304
## hamster_3              0.662 0.439
## hamster_4              0.623 0.388
## hamster_rabbit         0.580 0.336
## hamster_tadpole        0.634 0.402
## hedgehog_3             0.584 0.341
## hedgehog_4             0.649 0.421
## hedgehog_iguana        0.466 0.217
## hedgehog_owl           0.643 0.413
## hoe_3                  0.440 0.193
## hoe_4                  0.552 0.305
## hoe_dustpan            0.442 0.195
## hopscotch_3            0.595 0.355
## hopscotch_4            0.605 0.367
## hopscotch_chalk        0.618 0.381
## hopscotch_funnel       0.654 0.428
## horn_3                 0.452 0.205
## horn_4                 0.500 0.250
## horn_chin              0.437 0.191
## kimono_3               0.652 0.425
## kimono_4               0.542 0.294
## kimono_cardigan        0.628 0.394
## koala_4                0.595 0.354
## latch_3                0.612 0.375
## latch_4                0.605 0.366
## latch_lock             0.550 0.302
## latch_microwave        0.573 0.328
## locker_3               0.468 0.219
## locker_4               0.673 0.453
## locker_basket          0.610 0.372
## locker_cabinet         0.468 0.219
## lollipop_3             0.569 0.323
## lollipop_4             0.576 0.332
## lollipop_candy         0.612 0.375
## mandolin_3             0.504 0.254
## mandolin_mast          0.511 0.261
## map_3                  0.517 0.267
## map_4                  0.568 0.322
## map_ruby               0.495 0.245
## marshmallow_3          0.521 0.272
## marshmallow_4          0.522 0.273
## marshmallow_dryer      0.486 0.236
## mulch_3                0.634 0.402
## mulch_4                0.604 0.365
## mulch_clarinet         0.497 0.247
## net_3                  0.680 0.462
## net_4                  0.707 0.500
## net_domino             0.620 0.384
## net_tee                0.681 0.463
## oatmeal_3              0.507 0.258
## oatmeal_4              0.567 0.322
## oatmeal_cereal         0.476 0.226
## oatmeal_collar         0.474 0.225
## oil_4                  0.435 0.189
## oil_gel                0.554 0.307
## oil_grain              0.636 0.405
## omelet_3               0.692 0.479
## omelet_4               0.659 0.434
## omelet_clay            0.631 0.399
## otter_3                0.581 0.337
## otter_4                0.439 0.192
## otter_beaver           0.660 0.436
## otter_goose            0.442 0.195
## pajamas_3              0.592 0.351
## pajamas_4              0.436 0.190
## parsley_3              0.627 0.393
## parsley_4              0.651 0.424
## parsley_basil          0.516 0.267
## parsley_tiara          0.577 0.333
## pie_3                  0.508 0.258
## pie_4                  0.584 0.341
## pie_flashlight         0.484 0.235
## pie_pizza              0.588 0.346
## pistachio_3            0.555 0.307
## pistachio_4            0.687 0.472
## pistachio_avocado      0.448 0.200
## pistachio_taco         0.471 0.222
## pitcher_3              0.714 0.510
## pitcher_4              0.789 0.623
## pitcher_batter         0.734 0.538
## pitcher_tumbleweed     0.582 0.339
## potato_3               0.554 0.307
## potato_4               0.639 0.408
## potato_glasses         0.627 0.394
## potato_pot             0.610 0.372
## prism_3                0.686 0.471
## prism_4                0.717 0.515
## prism_radar            0.686 0.470
## prism_subway           0.496 0.246
## prune_3                0.531 0.282
## prune_4                0.614 0.377
## prune_headrest         0.608 0.370
## prune_raisin           0.467 0.218
## puddle_3               0.596 0.355
## puddle_4               0.660 0.436
## puddle_baseball        0.582 0.339
## puddle_bubble          0.494 0.244
## pump_3                 0.549 0.301
## pump_4                 0.631 0.399
## pump_plug              0.640 0.410
## rice_3                 0.597 0.356
## rice_4                 0.541 0.293
## rice_box               0.537 0.288
## rice_dice              0.571 0.326
## saddle_3               0.536 0.287
## saddle_4               0.665 0.442
## saddle_handle          0.623 0.388
## sandbag_4              0.439 0.193
## sandbag_sandpaper      0.451 0.203
## sandbag_valve          0.546 0.298
## sauerkraut_potpourri   0.492 0.242
## scaffolding_3          0.498 0.248
## scaffolding_4          0.603 0.364
## scaffolding_satellite  0.647 0.418
## scaffolding_veil       0.513 0.263
## scoop_3                0.518 0.269
## scoop_4                0.547 0.299
## scoop_sauce            0.493 0.244
## scoop_statue           0.474 0.224
## seagull_3              0.575 0.331
## seagull_4              0.496 0.246
## seagull_lobster        0.552 0.304
## seagull_pigeon         0.637 0.406
## ship_3                 0.553 0.305
## ship_4                 0.519 0.269
## ship_boat              0.625 0.391
## ship_nose              0.488 0.238
## shower_3               0.601 0.361
## shower_4               0.605 0.366
## shower_soap            0.513 0.263
## shower_toe             0.552 0.305
## silverware_3           0.594 0.353
## silverware_4           0.659 0.435
## silverware_garlic      0.560 0.314
## silverware_trophy      0.644 0.415
## sink_3                 0.661 0.437
## sink_4                 0.583 0.340
## sink_stair             0.661 0.436
## sink_toilet            0.580 0.336
## ski_3                  0.499 0.249
## ski_4                  0.629 0.396
## ski_ice                0.596 0.355
## ski_suitcase           0.633 0.401
## sloth_4                0.474 0.225
## sloth_clam             0.514 0.264
## sloth_donkey           0.433 0.187
## snail_3                0.592 0.351
## snail_4                0.588 0.346
## snail_cow              0.532 0.283
## snail_worm             0.619 0.383
## sorbet_3               0.470 0.221
## sorbet_4               0.550 0.302
## sorbet_palette         0.469 0.220
## spatula_3              0.489 0.239
## spatula_lava           0.593 0.352
## sprinkler_3            0.590 0.348
## sprinkler_4            0.612 0.375
## sprinkler_flower       0.532 0.283
## sprinkler_sparkler     0.486 0.236
## squash_3               0.587 0.344
## squash_4               0.481 0.232
## squash_branch          0.528 0.279
## squash_pumpkin         0.571 0.326
## squirrel_3             0.603 0.363
## squirrel_4             0.478 0.228
## squirrel_eagle         0.591 0.350
## squirrel_monkey        0.477 0.228
## stew_3                 0.628 0.394
## stew_mixer             0.554 0.307
## stew_soup              0.644 0.415
## stretcher_3            0.681 0.463
## stretcher_4            0.651 0.424
## stretcher_drill        0.599 0.359
## stump_3                0.660 0.436
## stump_4                0.589 0.347
## stump_bookshelf        0.715 0.511
## stump_log              0.707 0.500
## sunflower_3            0.469 0.220
## sunflower_4            0.538 0.290
## sunflower_pineapple    0.608 0.370
## swordfish_3            0.566 0.320
## swordfish_4            0.608 0.370
## swordfish_catfish      0.737 0.543
## swordfish_leopard      0.567 0.321
## taillight_4            0.523 0.274
## taillight_luggage      0.670 0.449
## tapestry_4             0.571 0.326
## tapestry_backdrop      0.575 0.330
## tapestry_leggings      0.506 0.256
## teabag_3               0.482 0.232
## teabag_fingerprint     0.606 0.367
## teabag_tea             0.526 0.277
## teapot_brake           0.590 0.348
## teapot_teacup          0.509 0.259
## telescope_4            0.532 0.283
## telescope_tripod       0.435 0.189
## telescope_windshield   0.594 0.353
## thermos_3              0.666 0.443
## thermos_4              0.581 0.337
## thermos_firewood       0.651 0.424
## thermos_oilcan         0.669 0.447
## treasure_3             0.537 0.289
## treasure_4             0.443 0.196
## treasure_gift          0.480 0.230
## treasure_rope          0.559 0.313
## trumpet_3              0.478 0.228
## trumpet_4              0.660 0.435
## trumpet_flute          0.546 0.298
## trumpet_violin         0.445 0.198
## tulip_3                0.499 0.249
## tulip_4                0.565 0.319
## tulip_glove            0.663 0.440
## tulip_rose             0.574 0.330
## turbine_3              0.606 0.367
## turbine_4              0.560 0.313
## turbine_generator      0.634 0.402
## turbine_tab            0.681 0.464
## turkey_3               0.491 0.241
## turkey_4               0.604 0.365
## turkey_goat            0.660 0.436
## turtle_3               0.514 0.264
## turtle_4               0.509 0.259
## turtle_horse           0.650 0.423
## tuxedo_stroller        0.504 0.254
## tuxedo_suit            0.527 0.277
## typewriter_3           0.494 0.244
## typewriter_4           0.678 0.460
## typewriter_printer     0.531 0.282
## typewriter_sunglasses  0.635 0.403
## watermelon_3           0.499 0.249
## watermelon_4           0.478 0.228
## watermelon_screwdriver 0.519 0.269
## watermelon_strawberry  0.567 0.321
## waterwheel_3           0.558 0.311
## waterwheel_pinwheel    0.641 0.411
## waterwheel_target      0.454 0.206
## 
## SS loadings:  119.323 
## Proportion Var:  0.331 
## 
## Factor correlations: 
## 
##   F
## F 1

item_check_2 = rownames_to_column(item_check)

item_info <- item_check_2 %>%
  as_tibble() %>%
  rename(itemId = rowname) %>%
  mutate(targetWord = str_split_fixed(itemId, '_',2)[,1]) %>%
    mutate(numAFC = str_split_fixed(itemId, '_',2)[,2]) %>%
  left_join(word1_aoa) %>%
  group_by(targetWord) %>%
  mutate(avg_h2 = mean(h2)) %>%
  ungroup() %>%
  mutate(targetWordOrdered = fct_reorder(targetWord, avg_h2))

Here, each dot is a trial type, and the y-axis the factor loading for that trial – looks like we have a pretty even spead from .4 to .8

ggplot(data = item_info, aes(x=targetWordOrdered, y=F, col=AoA_Est_Word1)) +
  theme_few(base_size = 12) +
  geom_point()  +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=4), aspect.ratio = .4)

  # facet_wrap(~numAFC, nrow=1)

Now do some pre-processing to get individual theta estiamtes

Create sandbox participant who has seen all trials

list.items <- itemIdforIRTTrimmed$itemId
df.sandbox.participant <- 
  tibble(pid = "sandbox", itemId = list.items, correct = 1)

Create trial ID structure for matrix

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

df.rocketship.roar.trials.clean <- school_data_multi_afc_cleaned %>% 
  filter(pid %in% completed_twice$pid) %>%
  filter(!pid %in% duplicate_ids$pid) %>%
  filter(studyId %in%  c("school-multiAFC"), itemType == "test")  %>%
  left_join(itemIdforIRT, by=c('targetWord','numAFC')) %>%
  filter(itemId %in% list.items) %>%
  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, word2, sep = "_")) %>%
  filter(itemId %in% list.items) %>%
  dplyr :: distinct(pid, runId, studyId, start_time, numAFC, word1, word2, itemId, correct, rt)

There were 177 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_fn = mean, 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_fn = mean, 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_fn = mean, 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_fn = mean, 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_not_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_not_enough_items$pid) %>%
  filter(pid %in% unique(df.rocketship.roar.trials.clean$pid)) %>%
  distinct(pid, itemId, correct) %>% 
  rbind(df.sandbox.participant) %>%
  unique() %>% 
  pivot_wider(names_from = itemId, values_from = correct, values_fn = mean, values_fill = NA)

Get theta score based on these reseponse pattenrs

getTheta <- function(df){
  thetas <- fscores(mod.vocab.2pl.trimmed, method = "ML", response.pattern = df %>% select(-pid), max_theta = 8, min_theta = -8, theta_lim = c(-8, 8))
  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)
               )