DLL-evaluation_v2.Rmd

---
title: "DLL-ES IRT Analysis and Recommendations"
author: "George Kachergis"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r load-fitted-IRT-models, echo=F, message=F}
require(here)
require(mirt)
require(tidyverse)
require(gridExtra)

# use combined WS+WG data, or just WS?
load(here("data/eng_ws_wg_mod_2pl.Rds"))
coefs = list(en = coefs_2pl)
fscores = list(en = fscores_2pl)
mods = list(en = mod_2pl)

load(here("data/sp_wg_ws_mod_2pl.Rds"))
coefs$sp = coefs_2pl
fscores$sp = fscores_2pl
mods$sp = mod_2pl
rm(coefs_2pl, fscores_2pl, mod_2pl)
```

```{r load-cdi-short-forms, echo=F, message=F}
wg_short_en <- read_csv(here("DLL/eng-wg-short.csv")) # 89
ws_short_enA <- read_csv(here("DLL/eng-ws-shortA.csv")) # 100
#ws_short_enB <- read_csv(here("DLL/eng-ws-shortB.csv")) # 100 - note they did not use this one

#length(intersect(wg_short_en$word, coefs$en$definition)) # all match

#length(intersect(ws_short_enA$word, coefs$en$definition)) # all match, but:
# should "swing" (in Fenson2000) be (action) or (object) -- I made it object, since surrounded by nouns

#length(intersect(ws_short_enB$word, coefs$en$definition)) # all match
# replaced 'fish' with fish (animal) since surrounded by animals;
# 'chicken' -> chicken (food) given list context
# replaced 'feet' with 'foot' ...
# replaced 'drink' with drink (action)

# Spanish short forms
wg_short_sp <- read_csv(here("DLL/spanish-wg-short.csv"))
ws_short_sp <- read_csv(here("DLL/spanish-ws-short.csv"))
#(intersect(wg_short_sp$word, coefs$sp$definition)) # 102 - tambien not in wordbank
# 102 - wordbank errors: comer(se), adios (accent not on i), ver(se) ?, vasos (why plural?)
#setdiff(wg_short_sp$word, coefs$sp$definition)
#length(intersect(ws_short_sp$word, coefs$sp$definition)) # 100
```


## Goals

The goals are 1) to create a word list that is informative about both English and Spanish vocabulary size and 2) to ensure that there are sufficient doublets to estimate lexical overlap.
On an IRT view, we can't perfectly assess 2 (at least not without better bilingual CDI data), but we can assess criterion 1 - that is, we can look at whether the reduced word list is a good sub-test for the full CDI in each language. 
Our original concern was that the current DLL-ES test might not perform well for older or high ability kids due to the lack of abstract words, and we can test this formally.

The DLL lists are meant to be used together with the original English and Spanish MCDI short forms.

```{r load-all-DLL-lists, echo=F, message=F}
dll1short_raw <- read_csv(here("DLL/Sandy/DLL-ES Supplement Words - DLL-ES Short Form Supplements (Level 1).csv")) %>% 
  rename(english=`English Word`, 
         spanish=`Spanish Word`) 
dll2short_raw <- read_csv(here("DLL/Sandy/DLL-ES Supplement Words - DLL-ES Short Form Supplements (Level 2).csv")) %>% 
  rename(english=`English Word`, 
         spanish=`Spanish Word`) 

# Level 1 exclusions: examine to see if we find some equal difficulty, informative matches?
dll1exclusions <- read_csv(here("DLL/Sandy/DLL-ES Supplement Words - DLL-ES Long Form Excluded Matches (Level 1).csv")) %>% 
  rename(english=`English Word`, 
         spanish=`Spanish Word`) # 94
dll2exclusions <- read_csv(here("DLL/Sandy/DLL-ES Supplement Words - DLL-ES Long Form Excluded Matches (Level 2).csv")) %>% 
  rename(english=`English Word`, 
         spanish=`Spanish Word`) # 210

dll1long_raw <- read_csv(here("DLL/Sandy/DLL-ES Supplement Words - DLL-ES Long Form Supplement (Level 1).csv")) %>% 
  rename(english=`English Word`, 
         spanish=`Spanish Word`) 
dll2long_raw <- read_csv(here("DLL/Sandy/DLL-ES Supplement Words - DLL-ES Long Form Supplement (Level 2).csv")) %>% 
  rename(english=`English Word`, 
         spanish=`Spanish Word`) 
```

Now we need to match the DLL words to the wordbank items for which we have IRT parameters.

### Matching DLL short Level 1

```{r match-english-words, message=F, warning=F}
# use with EN WG short form (12-18 mos)
#dll1ENshort <- read_csv(here("DLL/DLL-ES1-short-English.csv")) # 79 items (81 after splitting defs)

dll1short <- dll1short_raw %>% # this is now 168 items...it has the WG short form included
  mutate(english = tolower(english)) %>%
  mutate(english = case_when(english == 'i' ~ 'I',
                          english == 'tv (television)' ~ 'TV',
                          english == 'water' ~ 'water (beverage)',
                          english == 'grandma (or word used in your family)' ~ 'grandma*',
                          english == 'mommy (or word used in your family)' ~ 'mommy*',
                          english == 'choo choo (train sound)' ~ 'choo choo',
                          english == 'patty cake' ~ 'pattycake',
                          english == 'bye or bye bye' ~ 'bye',
                          english == 'teddy bear' ~ 'teddybear',
                          english == 'his/her' ~ 'his', # hers ? (and should DLL say 'his/hers' ?)
                          english == 'to have' ~ 'have',
                          english == 'shh' ~ 'shh/shush/hush',
                          english == 'to sit' ~ 'sit',
                          english == 'to be' ~ 'be',
                          english == 'put, put on' ~ 'put',
                          english == 'to write' ~ 'write',
                          english == 'arms' ~ 'arm',
                          english == 'church (or word used in your family)' ~ 'church*',
                          english == 'want' ~ 'wanna/want to',
                          english == 'in' ~ 'inside/in', # into ?
                          TRUE ~ english))

dll1ENshort_num_matching = length(intersect(dll1short$english, coefs$en$definition)) 
# 75/81 (GK version)
# 160/168 (Sandy)
#setdiff(dll1short$english, coefs$en$definition) 

# are these items also on the short form CDIs?
#length(intersect(dll1short$english, wg_short_en$word)) # all of the WG
#length(intersect(dll1short$english, ws_short_enA$word)) # 36 of the WS
```

```{r match-spanish-words, message=F, warning=F}
# use with SP WG short form (12-18 mos)
#dll1SPshort <- read_csv(here("DLL/DLL-ES1-short-Spanish.csv")) # 67 items

dll1short <- dll1short %>% # this is now 168 items...it has the WG short form included
  mutate(spanish = tolower(spanish)) %>%
  separate(col = spanish, into = c("spanish", NA), sep=" \\(") %>%
  mutate(spanish = case_when(spanish == 'mamá/mami' ~ 'mamá',
                             spanish == 'calcetines' ~ 'calcetín',
                             spanish == 'tomar baño / bañarse' ~ 'baño', # verb -> noun.. tomar(se) ?
                             spanish == 'espera' ~ 'esperar(se)', # close enough ?
                             spanish == 'acabar(se)' ~ 'acabar',
                             spanish == 'no hay más' ~ 'no hay', # or "más" ?
                             spanish == 'rapido' ~ 'rápido (descriptive)', # or rápido (quantifiers)
                             spanish == 'lastimado' ~ 'lastimar(se)', # close enough ?
                             spanish == 'otro' ~ 'otro/otra vez', # close enough ?
                             spanish == 'quiquiriqui' ~ 'quiquiriquí', # DLL missing acent
                             spanish == 'brazos' ~ 'brazo',
                             spanish == 'manos' ~ 'mano',
                             spanish == 'vaso' ~ 'vasos',
                             spanish == 'llaves' ~ 'llave',
                             spanish == 'adiós/byebye' ~ 'adíos/byebye', # wordbank accent is incorrect
                             spanish == 'uno, dos, tres' ~ 'uno dos tres...',
                             spanish == 'shh' ~ 'shhh',
                             spanish == 'ver' ~ 'ver(se)',
                             spanish == '¿dónde está?' ~ 'dónde', # close enough ? 
                             TRUE ~ spanish))

dll1SPshort_num_matching = length(intersect(dll1short$spanish, coefs$sp$definition)) # was 59 on GK version
# 158 / 168
#setdiff(dll1short$spanish, coefs$sp$definition)
# no match: tostada, alimentar, sonreír, algunos, también

# are these items also on the short form CDIs?
#length(intersect(dll1short$spanish, wg_short_sp$word)) # all of the WG
#length(intersect(dll1short$spanish, ws_short_sp$word)) # 63 of the WS
```

### Matching DLL short Level 2

```{r dll2-short, message=F, warning=F}
# DLL short Level 2 English / Spanish
dll2short <- dll2short_raw %>% mutate(english = tolower(english)) %>% filter(!is.na(english)) %>%
  mutate(english = case_when(english == 'church (or word used in your family)' ~ 'church*',
                             english == 'tv' ~ 'TV',
                             english == 'to be' ~ 'be',
                             english == 'to be, there is' ~ 'there', # DLL typo?
                             english == 'mommy' ~ 'mommy*',
                             english == 'in the morning' ~ 'morning',
                             english == 'swing' ~ 'swing (object)', # or swing (action) ?
                             english == 'dress' ~ 'dress (object)',
                          TRUE ~ english))
# not on CDI: snake, drum, rice, skirt, mustache, pot, matches, newspaper, bell, godmother, 
# let's go (but "go"), know, turn on (but "on"), win

#length(intersect(dll2short$english, coefs$en$definition)) # 162 / 177
#setdiff(dll2short$english, coefs$en$definition) 

dll2short <- dll2short %>% 
  mutate(spanish = tolower(spanish)) %>%
  separate(col = spanish, into = c("spanish", NA), sep=" \\(") %>%
  mutate(spanish = case_when(spanish == 'zapatos' ~ 'zapato',
                             spanish == 'adiós/bye bye' ~ 'adíos/byebye', # wordbank typo
                             spanish == 'banco' ~ 'banco (outside)', # or banco (places)
                             spanish == 'coca' ~ 'soda/refresco', # close enough?
                             spanish == 'acabar/terminar' ~ 'acabar', # or terminar
                             spanish == 'quiquiriqui' ~ 'quiquiriquí', # DLL missing accent
                             #spanish == 'debajo' ~ 'abajo', # down vs. below...close enough?
                             spanish == 'en la noche/esta noche' ~ 'en la noche',
                             spanish == 'haber' ~ 'haber (hay)',
                             spanish == 'ir de compras' ~ 'comprar', # close enough? or 'ir(se)'
                             spanish == 'llevar' ~ 'llevar(se)',
                             #spanish == 'no más' ~ 'más', # close enough? or 'no'..or 'no hay'
                             spanish == 'puede' ~ 'poder', # infinitive--close enough?
                             spanish == 'que' ~ 'que (connection word)',
                             spanish == 'rapido' ~ 'rápido (descriptive)', # or rápido (quantifiers)
                             TRUE ~ spanish))
                             
length(intersect(dll2short$spanish, coefs$sp$definition)) # 153 / 177
sort(setdiff(dll2short$spanish, coefs$sp$definition) )

#length(intersect(ws_short_enA$word, dll2short$english)) # 100
#length(intersect(ws_short_enB$word, dll2short$english)) # ..
```

### Matching DLL long supplement Level 1

```{r dll1-long, message=F}
dll1long <- dll1long_raw %>% mutate(english = tolower(english)) %>%
  mutate(english = case_when(english == 'daddy (or word used in your family)' ~ 'daddy*',
                          english == 'toy' ~ 'toy (object)',
                          english == 'swing' ~ 'swing (object)', # or swing (action) ?
                          english == 'dress' ~ 'dress (object)',
                          TRUE ~ english))


dll1long_EN_num_matching = length(intersect(dll1long$english, coefs$en$definition)) # 74 / 74 match

dll1long <- dll1long %>% mutate(spanish = tolower(spanish)) %>%
  separate(col = spanish, into = c("spanish", NA), sep=" \\(") %>%
  mutate(spanish = case_when(spanish == 'pipi' ~ 'pipí', # is this a match?
                             spanish == 'orejas' ~ 'oreja',
                             spanish == 'dedos' ~ 'dedo',
                             spanish == 'escalera' ~ 'escaleras', 
                             spanish == 'bolsa' ~ 'bolsa (clothing)', # or bolsa (household) ?
                             spanish == 'papá/papi' ~ 'papá', 
                             spanish == 'cosquillita' ~ 'cosquillitas', 
                             spanish == 'hacer la meme' ~ 'siesta', # close enough ?  or hacer?
                             TRUE ~ spanish))

dll1long_SP_num_matching = length(intersect(dll1long$spanish, coefs$sp$definition)) # 74 / 74 match
```

### Matching DLL long supplement Level 2

```{r dll2-long, message=F}
dll2long <- dll2long_raw %>% mutate(english = tolower(english)) %>%
  mutate(english = case_when(english == 'daddy (or name/word used in your family)' ~ 'daddy*',
                          english == 'teddy bear' ~ 'teddybear',
                          english == 'patty cake' ~ 'pattycake', # or swing (action) ?
                          english == 'dress' ~ 'dress (object)',
                          english == 'i' ~ 'I',
                          english == 'penis (or word used in your family)' ~ 'penis*',
                          english == 'water' ~ 'water (beverage)',
                          english == 'orange' ~ 'orange (food)',
                          english == 'clock/watch' ~ 'clock', # or watch (object)
                          english == 'drink' ~ 'drink (action)',
                          english == 'feet' ~ 'foot',
                          english == 'picture (\"or photo\")' ~ 'picture',
                          english == 'buttocks/bottom (or word used in your family)' ~ 'buttocks/bottom*',
                          TRUE ~ english))

dll2long <- dll2long %>% mutate(spanish = tolower(spanish)) %>%
  separate(col = spanish, into = c("spanish", NA), sep=" \\(") %>%
  mutate(spanish = case_when(spanish == 'pipi' ~ 'pipí',
                          spanish == 'shh' ~ 'shhh',
                          spanish == 'bolsa' ~ 'bolsa (clothing)', # or bolsa (household) ?
                          spanish == 'qué' ~ 'qué (question_words)',
                          TRUE ~ spanish))

#length(intersect(dll2long$spanish, coefs$sp$definition)) 
#setdiff(dll2long$spanish, coefs$sp$definition)
```


```{r load-wordbank-data, echo=F, message=F}
load(here("data/wordbank_eng_ws_wg_webcdi31-36mos.Rds"))

too_young <- which(d_demo$age < 12) # 378 children can't be producing any words yet

d_demo_en = d_demo[-too_young,] 
d_mat_en = d_mat[-too_young,]

# DLL items not named the same
#dll1ENshort[which(!is.element(dll1ENshort$word, colnames(d_mat))),]
#dll1SPshort[which(!is.element(dll1SPshort$word, coefs$sp$definition)),] # tostada, 
```


English DLL items not in our wordbank IRT model: *one, two, three*, *family*, *drum*, *good morning*, *also*, and *many*.
Spanish DLL items not in our wordbank IRT model: *tostada*, *algunos*, *alimentar*, *sonreír*, *no hay más* (although we have *no* and *no hay*, as well as *más*), and lastimado (but we have *lastimar(se)*).

DLL2-Spanish has "esta noche" translated as "at night", but should be "tonight".
DLL also has "puede", a conjugation of infinitive verb *poder* (which is in Wordbank).
Wordbank and DLL also often disagree about plural, e.g. *escalera*, and *pierna* on the DLL correspond to *escaleras* and *piernas* in Wordbank. (See also *brazos*, *manos*, *vasos*, and *llaves*.)

## Does the DLL short form recover full form scores?

### English DLL Level 1: Production

Using data from `r nrow(subset(d_demo_en, age<19))` English-speaking children 12-18 month of age from Wordbank, we test how well sumscores from the DLL-ES1 short English form + CDI:WG short form predict children's  production scores from the full CDI (WG/WS).
The left panel shows full CDI scores vs. the DLL-ES1 short + CDI:WG short score, and the right panel shows the full CDI scores vs. just the CDI:WG short form score.

```{r cdi-short-DLL1-vs-full-cdi-English-prod, echo=F, message=F, fig.width=8, fig.height=4}
en_wg_short_cols = na.omit(match(wg_short_en$word, colnames(d_mat_en)))
en_wg_short_dll_cols = na.omit(match(dll1short$english, colnames(d_mat_en))) # 160 items
d_demo_en$production = rowSums(d_mat_en, na.rm=T)
d_demo_en$DLLsum = rowSums(d_mat_en[,en_wg_short_dll_cols], na.rm=T)
d_demo_en$WGshort = rowSums(d_mat_en[,en_wg_short_cols], na.rm=T)

young_en = d_demo_en %>% filter(age <= 18, production < 396) # want only WG kids

dll_en_cor = cor(young_en$production, young_en$DLLsum) 

p1 <- young_en %>%
  ggplot(aes(x=production, y=DLLsum, group=sex, color=sex)) + 
  geom_point(alpha=.3) + theme_classic() + 
  xlab("Full CDI:WG Production Score") + ylab("CDI WG Short + DLL Score") + 
  geom_abline(slope=length(en_wg_short_dll_cols) / 395, intercept=0, linetype = 'dashed') +
  geom_smooth()

p2 <- young_en %>%
  ggplot(aes(x=production, y=WGshort, group=sex, color=sex)) + 
  geom_point(alpha=.3) + theme_classic() + 
  xlab("Full CDI:WG Production Score") + ylab("CDI WG Short Score") + 
  geom_abline(slope=length(en_wg_short_cols) / 395, intercept=0, linetype = 'dashed') +
  geom_smooth() + theme(legend.position = "none")
  
grid.arrange(p1, p2, nrow = 1, widths=c(4.5,3.9), heights=3.5)
```

Overall, the correlation of children's CDI:WG short + DLL scores and their full CDI production scores is quite high ($r=`r round(dll_en_cor,2)`$), but as shown above, for small vocabulary sizes the DLL score overestimates full CDI:WG production scores, while for higher full CDI:WG scores the DLL underestimates vocab size (dotted line has slope $=`r length(en_wg_short_dll_cols)` / `r 395`$).
However, the CDI:WG short form alone (right panel) shows a similar (and more extreme) overestimation for small vocabulary sizes.

### English DLL Level 1: Comprehension

Now we do the same for comprehension (receptive vocabulary) using Wordbank's CDI:WG English data.

```{r cdi-short-DLL1-vs-full-cdi-English-comp, echo=F, message=F, fig.width=8, fig.height=4}
load(here("data/comp/wordbank_eng_wg_webcdi.Rds"))

too_young <- which(d_demo$age < 12) # 378 children can't be producing any words yet

d_demo_en_comp = d_demo[-too_young,] 
d_mat_en_comp = d_mat_wg[-too_young,]

en_wg_short_colsC = na.omit(match(wg_short_en$word, colnames(d_mat_en_comp)))
en_wg_short_dll_colsC = na.omit(match(dll1short$english, colnames(d_mat_en_comp))) # 160 items
d_demo_en_comp$comprehension = rowSums(d_mat_en_comp, na.rm=T)
d_demo_en_comp$DLLsum = rowSums(d_mat_en_comp[,en_wg_short_dll_colsC], na.rm=T)
d_demo_en_comp$WGshort = rowSums(d_mat_en_comp[,en_wg_short_colsC], na.rm=T)

dll_en_cor_comp = cor(d_demo_en_comp$comprehension, d_demo_en_comp$DLLsum) 

p1 <- d_demo_en_comp %>%
  ggplot(aes(x=comprehension, y=DLLsum, group=sex, color=sex)) + 
  geom_point(alpha=.3) + theme_classic() + 
  xlab("Full CDI:WG Comprehension Score") + ylab("CDI WG Short + DLL Score") + 
  geom_abline(slope=length(en_wg_short_dll_colsC) / 395, intercept=0, linetype = 'dashed') +
  geom_smooth()

p2 <- d_demo_en_comp %>%
  ggplot(aes(x=comprehension, y=WGshort, group=sex, color=sex)) + 
  geom_point(alpha=.3) + theme_classic() + 
  xlab("Full CDI:WG Comprehension Score") + ylab("CDI WG Short Score") + 
  geom_abline(slope=length(en_wg_short_colsC) / 395, intercept=0, linetype = 'dashed') +
  geom_smooth() + theme(legend.position = "none")
  
grid.arrange(p1, p2, nrow = 1, widths=c(4.5,3.9), heights=3.5)
```

The correlation of children's CDI:WG short + DLL scores and their full CDI comprehension scores is quite high ($r=`r round(dll_en_cor_comp,2)`$), and extrapolating from DLL to full CDI:WG scores shows very little overestimation (less than when using the CDI:WG short form alone).

### English DLL Level 2 (Production)

Using data from `r nrow(subset(d_demo_en, age>15 & age<31))` English-speaking children 16-30 month of age from Wordbank, we test how well sumscores from the DLL-ES1 short English form + CDI:WG short form predict children's full production scores from the CDI:WS.
The left panel shows full CDI scores vs. the DLL-ES2 short + CDI:WS short form (A) score, and the right panel shows the full CDI scores vs. just the CDI:WG short form (A) score.

```{r cdi-short-DLL2-vs-full-cdi-English, echo=F, message=F, fig.width=8, fig.height=4}
en_ws_short_cols = na.omit(match(ws_short_enA$word, colnames(d_mat_en))) # 
en_ws_short_dll_cols = na.omit(match(dll2short$english, colnames(d_mat_en))) # 162 items
d_demo_en$production = rowSums(d_mat_en, na.rm=T)
d_demo_en$DLL2sum = rowSums(d_mat_en[,en_ws_short_dll_cols], na.rm=T)
d_demo_en$WSshort = rowSums(d_mat_en[,en_ws_short_cols], na.rm=T)

old_en = d_demo_en %>% filter(age > 15, age < 31)

dll2_en_cor = cor(old_en$production, old_en$DLLsum) 


p1 <- old_en %>% 
  ggplot(aes(x=production, y=DLL2sum, group=sex, color=sex)) + 
  geom_point(alpha=.1) + theme_classic() + 
  xlab("Full CDI Production Score") + ylab("CDI WS Short + DLL2 Score") + 
  geom_abline(slope=length(en_ws_short_dll_cols) / 680, intercept=0, linetype = 'dashed') +
  geom_smooth()

p2 <- old_en %>% 
  ggplot(aes(x=production, y=WSshort, group=sex, color=sex)) + 
  geom_point(alpha=.1) + theme_classic() + 
  xlab("Full CDI Production Score") + ylab("CDI WS Short Score") + 
  geom_abline(slope=length(en_ws_short_cols) / 680, intercept=0, linetype = 'dashed') +
  geom_smooth() + theme(legend.position = "none")

grid.arrange(p1, p2, nrow = 1, widths=c(4.5,3.9), heights=3.5)
```

Overall, the correlation of children's CDI:WS short + DLL2 scores and their full CDI production scores is quite high ($r=`r round(dll2_en_cor,2)`$), but as shown above, the DLL2 again mostly overestimates  production scores on the full CDI (dotted line has slope $=`r length(en_ws_short_dll_cols)` / `r 680`$).
In comparison, the CDI:WS short form (A) score only overestimates full CDI scores for smaller vocabulary sizes (<400).


### Spanish DLL Level 1

Now we look at overestimation for Spanish DLLs + CDI short forms.

```{r DLL-vs-full-score-Spanish, echo=F, message=F, fig.width=8, fig.height=4}
load(here("data/wordbank_sp_ws_wg_webcdi12-30mos.Rds"))
d_demo_sp = d_demo
d_mat_sp = d_mat

sp_wg_short_cols = na.omit(match(wg_short_sp$word, colnames(d_mat_sp)))
sp_ws_short_cols = na.omit(match(ws_short_sp$word, colnames(d_mat_sp)))


sp_wg_short_dll_cols = na.omit(match(dll1short$spanish, colnames(d_mat_sp)))
sp_ws_short_dll_cols = na.omit(match(dll2short$spanish, colnames(d_mat_sp)))

d_demo_sp$production = rowSums(d_mat, na.rm=T)
d_demo_sp$DLL1sum = rowSums(d_mat_sp[,sp_wg_short_dll_cols], na.rm=T)
d_demo_sp$DLL2sum = rowSums(d_mat_sp[,sp_ws_short_dll_cols], na.rm=T)
d_demo_sp$WGshort = rowSums(d_mat_sp[,sp_wg_short_cols], na.rm=T)
d_demo_sp$WSshort = rowSums(d_mat_sp[,sp_ws_short_cols], na.rm=T)

young_sp <- d_demo_sp %>% filter(age<19, production<=428)
old_sp <- d_demo_sp %>% filter(age>15)

dll1_sp_cor = cor(young_sp$production, young_sp$DLL1sum) 
dll2_sp_cor = cor(old_sp$production, old_sp$DLL2sum) 

p1 <- young_sp %>% ggplot(aes(x=production, y=DLL1sum, group=sex, color=sex)) + 
  geom_point(alpha=.3) + theme_classic() + 
  xlab("Full CDI Production Score") + ylab("DLL1 Score") + 
  geom_abline(slope=length(sp_wg_short_dll_cols) / 428, intercept=0, linetype = 'dashed') +
  geom_smooth()

p2 <- young_sp %>% ggplot(aes(x=production, y=WGshort, group=sex, color=sex)) + 
  geom_point(alpha=.3) + theme_classic() + 
  xlab("Full CDI Production Score") + ylab("CDI WG Short Score") + 
  geom_abline(slope=length(sp_wg_short_cols) / 428, intercept=0, linetype = 'dashed') +
  geom_smooth() + theme(legend.position = "none")

grid.arrange(p1, p2, nrow = 1, widths=c(4.5,3.9), heights=3.5)
```


Using Wordbank data from `r nrow(young_sp)` Spanish-speaking children aged 12-18 months, we test how well sumscores from the DLL-ES1 short Spanish form correlate with children's full CDI:WG production scores.

As for English, the correlation of Spanish-speaking children's DLL scores and their full CDI:WG production scores is quite high ($r=`r round(dll1_sp_cor,2)`$), but as shown above, their DLL score overestimates the production score on the full CDI at smaller vocabulary sizes (dotted line has slope $=`r length(sp_wg_short_dll_cols)` / `r 428`$).
Do note that few children in this dataset have large productive vocabularies.

### Spanish DLL Level 1: Comprehension

Now we do the same for comprehension (receptive vocabulary) using Wordbank's CDI:WG Spanish data.

```{r cdi-short-DLL1-vs-full-cdi-Spanish-comp, echo=F, message=F, fig.width=8, fig.height=4}
load(here("data/comp/wordbank_sp_wg_webcdi.Rds"))

too_young <- which(d_demo$age < 12) # 242 children can't be producing any words yet...leave in for comp?

d_demo_sp_comp = d_demo[-too_young,] 
d_mat_sp_comp = d_mat_wg[-too_young,]

sp_wg_short_colsC = na.omit(match(wg_short_sp$word, colnames(d_mat_sp_comp)))
sp_wg_short_dll_colsC = na.omit(match(dll1short$spanish, colnames(d_mat_sp_comp))) # 160 items
d_demo_sp_comp$comprehension = rowSums(d_mat_sp_comp, na.rm=T)
d_demo_sp_comp$DLLsum = rowSums(d_mat_sp_comp[,sp_wg_short_dll_colsC], na.rm=T)
d_demo_sp_comp$WGshort = rowSums(d_mat_sp_comp[,sp_wg_short_colsC], na.rm=T)

dll_sp_cor_comp = cor(d_demo_sp_comp$comprehension, d_demo_sp_comp$DLLsum) 

p1 <- d_demo_sp_comp %>%
  ggplot(aes(x=comprehension, y=DLLsum, group=sex, color=sex)) + 
  geom_point(alpha=.3) + theme_classic() + 
  xlab("Full CDI:WG Comprehension Score") + ylab("CDI WG Short + DLL Score") + 
  geom_abline(slope=length(sp_wg_short_dll_colsC) / 428, intercept=0, linetype = 'dashed') +
  geom_smooth()

p2 <- d_demo_sp_comp %>%
  ggplot(aes(x=comprehension, y=WGshort, group=sex, color=sex)) + 
  geom_point(alpha=.3) + theme_classic() + 
  xlab("Full CDI:WG Comprehension Score") + ylab("CDI WG Short Score") + 
  geom_abline(slope=length(sp_wg_short_colsC) / 428, intercept=0, linetype = 'dashed') +
  geom_smooth() + theme(legend.position = "none")
  
grid.arrange(p1, p2, nrow = 1, widths=c(4.5,3.9), heights=3.5)
```

The correlation of children's CDI:WG short + DLL scores and their full CDI comprehension scores is quite high ($r=`r round(dll_sp_cor_comp,2)`$), and extrapolating from DLL to full CDI:WG scores shows very little overestimation--similar to the level shown by using the CDI:WG alone.

### Spanish DLL Level 2

```{r spanish-dll2, echo=F, message=F, fig.width=8, fig.height=4}

p1 <- old_sp %>% ggplot(aes(x=production, y=DLL1sum, group=sex, color=sex)) + 
  geom_point(alpha=.3) + theme_classic() + 
  xlab("Full CDI Production Score") + ylab("DLL2 Score") + 
  geom_abline(slope=length(sp_ws_short_dll_cols) / 680, intercept=0, linetype = 'dashed') +
  geom_smooth()

p2 <- old_sp %>% ggplot(aes(x=production, y=WGshort, group=sex, color=sex)) + 
  geom_point(alpha=.3) + theme_classic() + 
  xlab("Full CDI Production Score") + ylab("CDI WS Short Score") + 
  geom_abline(slope=length(sp_ws_short_cols) / 680, intercept=0, linetype = 'dashed') +
  geom_smooth() + theme(legend.position = "none")

grid.arrange(p1, p2, nrow = 1, widths=c(4.5,3.9), heights=3.5)

```


## Recommendations

```{r dll-item-ease, echo=F}
en_dll1_ease <- coefs$en %>% mutate(onDLL = 
                  ifelse(is.element(definition, dll1short$english), 1, 0)) %>%
  group_by(onDLL) %>% 
  summarise(easiness=mean(d), n=n())

en_dll2_ease <- coefs$en %>% mutate(onDLL = 
                  ifelse(is.element(definition, dll2short$english), 1, 0)) %>%
  group_by(onDLL) %>% 
  summarise(easiness=mean(d), n=n())

en_wg_short_ease <- coefs$en %>% mutate(onWGshort = 
                          ifelse(is.element(definition, wg_short_en$word), 1, 0)) %>%
  group_by(onWGshort) %>% 
  summarise(easiness=mean(d), n=n())

en_ws_short_ease <- coefs$en %>% mutate(onWSshort = 
                          ifelse(is.element(definition, ws_short_enA$word), 1, 0)) %>%
  group_by(onWSshort) %>% 
  summarise(easiness=mean(d), n=n())

sp_dll1_ease <- coefs$sp %>% mutate(onDLL = 
                  ifelse(is.element(definition, dll1short$spanish), 1, 0)) %>%
  group_by(onDLL) %>% 
  summarise(easiness=mean(d), n=n())

sp_dll2_ease <- coefs$sp %>% mutate(onDLL = 
                  ifelse(is.element(definition, dll2short$spanish), 1, 0)) %>%
  group_by(onDLL) %>% 
  summarise(easiness=mean(d), n=n())
```


Overall, it seems that many of the items on the DLL are somewhat easier than average, and thus these forms tend to overestimate children's full CDI scores (indeed, for items on the DLL1 English short form, the average easiness is `r round(en_dll1_ease[2,]$easiness,2)`, while the mean easiness of items not on the DLL is `r round(en_dll1_ease[1,]$easiness,2)`).
This is also true of the CDI:WG short English form: the average easiness is `r round(en_wg_short_ease[2,]$easiness,2)` and the average ease of items not on the WG short form is `r round(en_wg_short_ease[1,]$easiness,2)`.
The CDI:WS short English form (A) is less biased towards easy items: average easiness is `r round(en_ws_short_ease[2,]$easiness,2)` vs. `r round(en_ws_short_ease[1,]$easiness,2)` for items not on the short WS.
The histograms below show the distribution of easiness parameters for English (left) and Spanish (right) CDI words. 
Solid lines show the average ease of DLL items (DLL 1 = red, DLL 2 = orange), and dashed lines show the average of non-DLL items.

Spanish DLL1 items have an average ease of `r round(sp_dll1_ease[2,]$easiness,2)`, while other items on the full CDI have a mean ease of `r round(sp_dll1_ease[1,]$easiness,2)`.

Spanish DLL2 shows the least bias: items on it have an average ease of `r round(sp_dll2_ease[2,]$easiness,2)`, while other CDI items have a mean of `r round(en_dll2_ease[1,]$easiness,2)`.
English DLL2 items have an average ease of `r round(en_dll2_ease[2,]$easiness,2)`, while other items on the full CDI have a mean ease of `r round(en_dll2_ease[1,]$easiness,2)`.


```{r, echo=F, message=F, fig.width=8, fig.height=4}
p1 <- ggplot(coefs$en, aes(x=d)) + geom_histogram() + theme_bw() + 
  xlab("English Item Easiness") + 
  geom_vline(xintercept=en_dll1_ease$easiness, color="red", linetype=c("dashed","solid")) +
  geom_vline(xintercept=en_dll2_ease$easiness, color="orange", linetype=c("dashed","solid"))
p2 <- ggplot(coefs$sp, aes(x=d)) + geom_histogram() + theme_bw() + 
  xlab("Spanish Item Easiness") + 
  geom_vline(xintercept=sp_dll1_ease$easiness, color="red", linetype=c("dashed","solid")) +
  geom_vline(xintercept=sp_dll2_ease$easiness, color="orange", linetype=c("dashed","solid"))

grid.arrange(p1, p2, nrow = 1, widths=c(4,4), heights=3)
```


We recommend bringing the overall mean estimated IRT difficulty of the words selected for the DLLs closer to the mean difficulty of the words on the rest of the CDI.

To start, we examine IRT easiness parameters for the doublets on the existing DLL lists, looking for items with large mismatch between their English and Spanish ease.

## Do doublets have similar difficulties?

We want to whether assess doublet items have similar difficulty (operationalized by their IRT parameters) in English and in Spanish. 
For example, consider if "perro" was for some reason much more difficult than "dog", then you wouldn't want to include it because it wouldn't be a good item for estimating vocabulary overlap!

Below are shown the parameters for items from the DLL Level 1 short form: en_d = English easiness, sp_d = Spanish easiness), ordered by the most to least discrepant (difficulty difference squared); (_a1 columns show item discriminations (slopes), and sp_en_d_diff simply shows difference in Spanish and English easiness parameter).

### DLL Level 1 Short Form

For the short form, we merely want to identify items that have very different difficulty in English and Spanish, and recommend that researchers interested in estimating conceptual overlap in bilinguals not include these items in their calculations due to the bias.
Below is shown the distribution of the squared difference in difficulty for doublets from the DLL Level 1 short form.

```{r compare-doublet-difficulties-dll1-short, echo=F}
# a1 = discrimination, d = difficulty, g = lower-bound, u = upper-bound

pars_dll1short <- dll1short %>% 
  left_join(coefs$sp, by=c("spanish"="definition")) %>% 
  select(english, spanish, a1, d) %>% 
  rename(sp_a1 = a1, sp_d = d) %>% 
  left_join(coefs$en, by=c("english"="definition")) %>%
  select(english, spanish, sp_a1, sp_d, a1, d) %>% 
  rename(en_a1 = a1, en_d = d)

pars_dll1short <- pars_dll1short %>% 
  mutate(sp_en_d_diff = sp_d - en_d,
         d_diff_sq = (sp_d - en_d)^2) %>% 
  filter(!is.na(d_diff_sq)) %>%
  arrange(desc(d_diff_sq))

threshold = 7.17 # mean = 2.17, sd = 3.3... 2.17 + 1.5*3.3 = 7.17
bad_dll1 <- subset(pars_dll1short, d_diff_sq >= threshold)

hist(pars_dll1short$d_diff_sq, main="", xlab="(EN diff - SP diff)^2")
abline(v=threshold,col="red")
```

The mean squared difference in difficulty of doublets on the DLL Level 1 short form is `r round(mean(pars_dll1short$d_diff_sq), 2)` (SD=`r round(sd(pars_dll1short$d_diff_sq), 2)`), and as shown above this distribution is highly skewed: most doublets are fairly well-matched (median d_diff_sq=0.99), but a few items are extremely mismatched.
We recommend that the `r nrow(bad_dll1)` items (shown below) with a squared difficulty difference of 7.17 (mean + 1.5 * SD) or more be excluded from calculations of conceptual overlap.

```{r, echo=F}
knitr::kable(bad_dll1, digits=2)
```

### DLL Level 2 Short Form

```{r compare-doublet-difficulties-dll2-short, echo=F}
pars_dll2short <- dll2short %>% 
  left_join(coefs$sp, by=c("spanish"="definition")) %>% 
  select(english, spanish, a1, d) %>% 
  rename(sp_a1 = a1, sp_d = d) %>% 
  left_join(coefs$en, by=c("english"="definition")) %>%
  select(english, spanish, sp_a1, sp_d, a1, d) %>% 
  rename(en_a1 = a1, en_d = d)

pars_dll2short <- pars_dll2short %>% 
  mutate(sp_en_d_diff = sp_d - en_d,
         d_diff_sq = (sp_d - en_d)^2) %>% 
  filter(!is.na(d_diff_sq)) %>%
  arrange(desc(d_diff_sq))

threshold = 8.22 # M = 2.53, SD = 3.79  M + 1.5*SD = 8.22
bad_dll2 <- subset(pars_dll2short, d_diff_sq >= threshold)

hist(pars_dll2short$d_diff_sq, main="", xlab="(EN diff - SP diff)^2")
abline(v=threshold,col="red")
```

The mean squared difference in difficulty of doublets on the DLL Level 2 short form is `r round(mean(pars_dll2short$d_diff_sq), 2)` (SD=`r round(sd(pars_dll2short$d_diff_sq), 2)`), and once again while most doublets are fairly well-matched, several items are extremely mismatched.
We recommend that the `r nrow(bad_dll2)` items (shown below) with a squared difficulty difference of 8.22 (mean + 1.5 * SD) or more be excluded from calculations of conceptual overlap.

```{r, echo=F}
knitr::kable(bad_dll2, digits=2)
```

Now we examine the doublets on the long supplemental DLL forms, and make recommendations of mismatched items to swap out.

### DLL Level 1 Long Supplement

Below are shown the parameters for items from the DLL Level 1 long supplement.
Clearly some of these items have quite different difficulty levels in English and Spanish, and we should try to find items that are more equivalent and swap them onto the long supplement.

```{r compare-doublet-difficulties-dll1-long, echo=F}
# a1 = discrimination, d = difficulty, g = lower-bound, u = upper-bound

pars_dll <- dll1long %>% 
  left_join(coefs$sp, by=c("spanish"="definition")) %>% 
  select(english, spanish, a1, d) %>% 
  rename(sp_a1 = a1, sp_d = d) %>% 
  left_join(coefs$en, by=c("english"="definition")) %>%
  select(english, spanish, sp_a1, sp_d, a1, d) %>% 
  rename(en_a1 = a1, en_d = d)

pars_dll <- pars_dll %>% 
  mutate(sp_en_d_diff = sp_d - en_d,
         d_diff_sq = (sp_d - en_d)^2) %>% 
  filter(!is.na(d_diff_sq)) %>%
  arrange(desc(d_diff_sq))

knitr::kable(pars_dll, digits=2)

#pars_dll <- pars_dll %>% filter(!is.element(english, wg_short_en$word),
#                                !is.element(spanish, wg_short_sp$word))
# that gets rid of all of them

#pars_dll %>% ggplot(aes(x=en_d, y=sp_d)) +
#  geom_point() + theme_minimal() # geom_text_repel() # but can only use 1lang..

# mean(diff): Spanish = -.75,  English = -.25
#colMeans(pars_dll[,c("sp_d","en_d")])
```

### DLL Level 2 Long Supplement

```{r compare-doublet-difficulties-dll2, echo=F}
# a1 = discrimination, d = difficulty, g = lower-bound, u = upper-bound

pars_dll2 <- dll2long %>% 
  left_join(coefs$sp, by=c("spanish"="definition")) %>% 
  select(english, spanish, a1, d) %>% 
  rename(sp_a1 = a1, sp_d = d) %>% 
  left_join(coefs$en, by=c("english"="definition")) %>%
  select(english, spanish, sp_a1, sp_d, a1, d) %>% 
  rename(en_a1 = a1, en_d = d)

pars_dll2 <- pars_dll2 %>% 
  mutate(sp_en_d_diff = sp_d - en_d,
         d_diff_sq = (sp_d - en_d)^2) %>% 
  filter(!is.na(d_diff_sq)) %>%
  arrange(desc(d_diff_sq))

knitr::kable(pars_dll2, digits=2)

# mean(diff): Spanish = -.72,  English = -.35
#colMeans(pars_dll2[,c("sp_d","en_d")])
```

## Recommended Item Swaps

We will use Wordbank's unilemmas to find translation-equivalent pairs that have smaller d_diff_sq values than current DLL extended items.
We first get the English / Spanish unilemmas from wordbank (both WS and WG), and below simply show the Spanish vs. English easiness parameters.

```{r unilemmas, echo=F, message=F}
# created in get_L1-L2_unilemmas.R
load("data/English-Spanish-dict.Rdata")

# should we filter out the words from the 

# add unilemmas to DLL dfs
dll1long <- dll1long %>% 
  left_join(dict %>% select(english, uni_lemma, d_diff_sq, en_d, sp_d, lexical_class), by=c("english"="english")) %>% 
  select(-`English Word Origin`, -`Spanish Word Origin`, -`Spanish Alternative Translations (experimenter or MCDI)`) %>%
  group_by(english) %>% arrange(d_diff_sq) %>% slice(1) %>% ungroup() # just take item with minimal easiness difference

dll2long <- dll2long %>% 
  left_join(dict %>% select(english, uni_lemma, d_diff_sq, en_d, sp_d, lexical_class), by=c("english"="english")) %>%
  select(-`English Word Origin`, -`Spanish Word Origin`, -`Spanish Alternative Translations (experimenter or MCDI)`) %>%
  group_by(english) %>% arrange(d_diff_sq) %>% slice(1) %>% ungroup()

# do the same for short forms? DLL group doesn't want to, but maybe didn't understand..
# dll1short <- dll1short %>% 
#   left_join(dict %>% select(spanish, uni_lemma, d_diff_sq, en_d, sp_d, lexical_class), by=c("spanish"="spanish")) %>% select(-notes) %>%
#   group_by(word) %>% arrange(d_diff_sq) %>% slice(1) %>% ungroup()
# 
# dll2short <- dll2short %>% 
#   left_join(dict %>% select(spanish, uni_lemma, d_diff_sq, en_d, sp_d, lexical_class), by=c("spanish"="spanish")) %>% select(-notes) %>%
#   group_by(word) %>% arrange(d_diff_sq) %>% slice(1) %>% ungroup()

dict %>% ggplot(aes(x=en_d, y=sp_d, color=lexical_class)) + 
  geom_abline(slope=1, intercept=0, linetype = "dashed") +
  geom_point(alpha=.7) + geom_smooth(method='lm') + theme_bw() +
  xlab("English word easiness") + ylab("Spanish word easiness")

dict <- dict %>% filter(!is.element(english, dll1short$english),
                        !is.element(spanish, dll1short$spanish), # was 380, now only 218
                        !is.element(english, dll2short$english),
                        !is.element(spanish, dll2short$spanish)) # 183 

Nswaps = 25
```

Working with `r nrow(dict)` unilemmas that match both our English and Spanish IRT parameters and that are *not* already on the DLL short lists, for each DLL list we will simply swap the N=`r Nswaps` items with the largest easiness difference for items with minimal easiness difference.
We also attempt to find replacement items of the same lexical class.
We report the original DLL list's easiness SSE, as well as the improvement in (easiness SSE) after each swap is made.


```{r optimization-functions, echo=F}
# penalize for 1. overestimation of vocab size relative to full CDI
# 2. sum of squared difference in item difficulties

# find item with biggest squared difference and select a new pair with a small difference
# with smaller (est?) difference (and same lexical category?)
# with high k, almost always chooses argmax
propose_swap <- function(cur_list, dictn, k=5, match_lexical_class=T) {
  set.seed(42)
  item_to_replace = sample(1:nrow(cur_list), 1, prob=cur_list$d_diff_sq^k)
  if(match_lexical_class) {
    tmp = subset(dictn, lexical_class == cur_list[item_to_replace,]$lexical_class)
    if(nrow(tmp)>0) dictn = tmp 
  }
  replacement = sample(1:nrow(dictn), 1, prob = 1/dictn$d_diff_sq^k)
  improvement = cur_list[item_to_replace,]$d_diff_sq - dictn[replacement,]$d_diff_sq
  if(improvement <= 0) {
    print("No better word found!")
    return(list())
  }
  print(paste0("Replacing '", cur_list[item_to_replace,]$uni_lemma,"' with '", 
               dictn[replacement,]$uni_lemma,"' (SSE improvement = ",round(improvement,2),")"))
  return(list(remove = cur_list[item_to_replace,], 
              add = dictn[replacement,]))
}

# for each DLL list, propose a number of swaps
improve_DLL_list <- function(cur_list, dict, Nswaps = 10) {
  unilingual_items <- subset(cur_list, is.na(uni_lemma))
  dict_novel <- subset(dict, !is.element(uni_lemma, cur_list$uni_lemma))
  cur_list <- na.omit(cur_list) 
  print(paste("Original list item easiness SSE:",round(sum(cur_list$d_diff_sq), 2)))
  print(paste("Mean Spanish item easiness:",round(mean(cur_list$sp_d), 2)))
  print(paste("Mean English item easiness:",round(mean(cur_list$en_d), 2)))
  # remove items from dict that are already on DLL list
  print(paste("Selecting from",nrow(dict_novel),"words on both Eng/Sp CDIs that are not on the DLL."))
  removed = c()
  added = c()
  for(i in 1:Nswaps) {
     swap = propose_swap(cur_list, dict_novel, k=5)
     if(length(swap)>0) {
      removed = rbind(removed, swap$remove)
      added = rbind(added, swap$add)
      # remove added word from dict
      dict_novel = subset(dict_novel, uni_lemma!=swap$add$uni_lemma)
      # remove removed word from DLL
      cur_list = subset(cur_list, uni_lemma!=swap$remove$uni_lemma)
      swap$remove
      #swap$add$word = unlist(swap$add[,language])
      tmp = swap$add %>%
        select(english, spanish, uni_lemma, lexical_class, d_diff_sq, en_d, sp_d)
      cur_list = rbind(cur_list, tmp)
     }
  }
  print(paste("New list item easiness SSE:",round(sum(cur_list$d_diff_sq), 2)))
  print(paste("Mean Spanish item easiness:",round(mean(cur_list$sp_d), 2)))
  print(paste("Mean English item easiness:",round(mean(cur_list$en_d), 2)))
  
  return(list(new_list = rbind(cur_list, unilingual_items), added = added, removed = removed))
}

```

## Level 1 DLL Long Supplement

```{r}
new_dll1long <- improve_DLL_list(dll1long, dict, Nswaps=Nswaps)$new_list
write.csv(new_dll1long, file="DLL/new_DLL-ES1-long.csv")
```

Even swapping only 25 items reduced the total SSE by more than 50%.
Although it was not optimized for, the average ease of the items also decreased for both languages, coming closer to the mean.
We now do the same for the other DLL forms before determining whether the DLL overestimation has decreased.

## Level 2 DLL Long Supplement

```{r}
new_dll2long <- improve_DLL_list(dll2long, dict, Nswaps=Nswaps)$new_list
write.csv(new_dll2long, file="DLL/new_DLL-ES2-long.csv")
```

Once again, the total SSE decreased significantly, and the average ease of the items in both languages went down, and got closer to each other.


## Summary

For each of the DLL lists, swapping the `r Nswaps` items with the largest discrepancy between English and Spanish easiness for items of minimal discrepancy within the same lexical class resulted in substantially reducing the total easiness SSE, and also resulted in mean item easiness (in both languages) that are more equal and closer to the means of each language, and thus generally better for extrapolating to children's full CDI score. 
Recognizing that these swaps are chosen algorithmically, we recommend the DLL team to consider each of the above swaps and determine whether any key words have been removed, or whether any undesirable words have been added.