data-explorations.Rmd

---
title: "Spanish CDI explorations"
author: "Paulina & Mike"
date: "2022-10-13"
output: html_document
---

# Intro

Goal of this project is to use the five Spanish CDI datasets in Wordbank to try and investigate dialect variation in Spanish language acquisition. 

Here are some potential questions:

* Do the sumscores across the intersecting items look the same at the population level (correcting for demographics)?
* Do individual common items have similar psychometric / developmental properties?
* What are the properties of items with the same unilemma but different item definitions?
* What are the properties of items that are not shared across dialects? (Could also compare Spanish (Mexican) between monolingual and bilingual populations)


```{r setup}
library(tidyverse)
library(wordbankr)
library(arm)
```

# Data loading

Start with summary scores. 

```{r}
eu_ws <- get_administration_data(language = "Spanish (European)", 
                                 form = "WS", 
                                 include_demographic_info = TRUE,
                                 include_language_exposure = TRUE) %>% 
  filter(language_exposures=="NULL")
mx_ws <- get_administration_data(language = "Spanish (Mexican)", 
                                 form = "WS", 
                                 include_demographic_info = TRUE,
                                 include_language_exposure = TRUE) %>% 
  filter(language_exposures=="NULL")
pr_ws <- get_administration_data(language = "Spanish (Peruvian)", 
                                 form = "WS", 
                                 include_demographic_info = TRUE,
                                 include_language_exposure = TRUE) %>% 
  filter(language_exposures=="NULL")
ar_ws <- get_administration_data(language = "Spanish (Argentinian)", 
                                 form = "WS", 
                                 include_demographic_info = TRUE,
                                 include_language_exposure = TRUE) %>% 
  filter(language_exposures=="NULL")

sp_ws <- bind_rows(eu_ws, 
                   mx_ws, 
                   pr_ws, 
                   ar_ws)
```

Make a plot!

```{r}
ggplot(sp_ws, aes(x = age, y = production)) + 
  geom_jitter(width = .2, alpha = .2) + 
  geom_smooth() + 
  facet_wrap(~language)
```


# Comparison on the intersection of items

QUESTION: Do the sumscores across the intersecting items look the same at the population level (correcting for demographics)?

```{r}
langs <- c("Spanish (European)", "Spanish (Mexican)", "Spanish (Peruvian)", "Spanish (Argentinian)")

d_ws <- map_df(langs, function(x) get_instrument_data(language = x, 
                                                      form = "WS", 
                                                      administration_info = TRUE, 
                                                      item_info = TRUE,
                                                      include_language_exposure = TRUE)) %>% 
  filter(language_exposures == "NULL")
```

Find the overlapping unilemmas.

For now, pull those unilemmas that are 1) in all languages, 2) only once in each language. 

```{r}
items <- map_df(langs, function(x) get_item_data(language = x, form = "WS"))

intersection <- items |>
  group_by(uni_lemma) |>
  summarise(n_langs = length(unique(language)), 
            n = n()) |>
  filter(n_langs == 4, n == 4) |>
  pull(uni_lemma)
# 224 common items
``` 

Filter data and replot. 

```{r}


ggplot(ms_ws, aes(x = age, y = produces)) + 
  geom_jitter(width = .2, alpha = .2) + 
  geom_smooth() + 
  facet_wrap(~language)
```



# Comparison of developmental properties

QUESTION: Do individual common items have similar psychometric / developmental properties?

```{r}
source("scripts/fit_models.R")

wb_data <- d_ws |> 
  filter(uni_lemma %in% intersection) |>
  group_by(uni_lemma, language, age) |> 
  summarise(total = n(),
            num_true = sum(produces, na.rm = TRUE))

aoas <- fit_aoas(wb_data)
```
DIF
```{r}
#language <- round(runif(nrow(df.data))) |> as.factor()

#fit.2plmulti <- multipleGroup(df.data,
                             # model = modstr,
                            #  itemtype = "2PL",
                             # group = grp,
                             # verbose = FALSE)
```


```{r}
ggplot(data = aoas) +
  geom_histogram(aes(x = aoa))+ 
  facet_wrap(~language)
```

Calculate correlations between languages

```{r}
cor_data <- aoas |> 
  pivot_wider(id_cols = uni_lemma, 
              names_from = language,
              values_from = aoa) |> 
  ungroup()

correl <- cor(cor_data |> dplyr::select(-"uni_lemma"), use = "complete.obs")
```

Plot correlogram
```{r}
library(corrgram)
corrgram(correl, type = "cor", panel = panel.cor)
```

```{r}
plot_data <- cor_data |> 
  pivot_longer(cols = c(`Spanish (Argentinian)`, `Spanish (Mexican)`, `Spanish (Peruvian)`), 
               names_to = "Dialect")
ggplot(plot_data, aes(x = `Spanish (European)`, y = value, col = Dialect)) + 
  geom_abline(slope = 1, intercept = 0, lty = "dashed", col = "#444444") +
  geom_point(alpha=.3) +
  geom_smooth(method = "lm") +
  theme_classic() + 
  theme(legend.position = "bottom") +
  labs(x = "Spanish (European) AoA",
       y = "Spanish (Latin American) AoA") +
  coord_cartesian(xlim = c(0, 40), ylim = c(0, 40))
```

Calculate differences in age of acquisition (AoA) between all dialects for each unilemma
```{r}
clean_aoas <- aoas |>
  dplyr::select(uni_lemma, language, aoa) |>
  pivot_wider(names_from = language, values_from = aoa)

difference_argentinian_european <- data.frame(difference_argentinian_european = c(clean_aoas$`Spanish (Argentinian)`- clean_aoas$`Spanish (European)`))

difference_argentinian_mexican <- data.frame(difference_argentinian_mexican = c(clean_aoas$`Spanish (Argentinian)`- clean_aoas$`Spanish (Mexican)`))

difference_argentinian_peruvian <- data.frame(difference_argentinian_peruvian = c(clean_aoas$`Spanish (Argentinian)`- clean_aoas$`Spanish (Peruvian)`))

difference_european_mexican <- data.frame(difference_european_mexican = c(clean_aoas$`Spanish (European)`- clean_aoas$`Spanish (Mexican)`))

difference_european_peruvian <- data.frame(difference_european_peruvian = c(clean_aoas$`Spanish (European)`- clean_aoas$`Spanish (Peruvian)`))

difference_mexican_peruvian <- data.frame(difference_mexican_peruvian = c(clean_aoas$`Spanish (Mexican)`- clean_aoas$`Spanish (Peruvian)`))

difference_AoA <- data.frame(uni_lemma = c(clean_aoas$uni_lemma), difference_argentinian_european, difference_argentinian_mexican, difference_argentinian_peruvian, difference_european_mexican, difference_european_peruvian, difference_mexican_peruvian)
```
Comparing identical/nonidentical surface forms
```{r}
# Filter "items" data-set to "words" only, then filter again so that it contains only those uni_lemmas contained in the "aoas" data frame
items_filtered <- filter(items, items$item_kind == 'word') %>% 
  filter(uni_lemma %in% c((clean_aoas$uni_lemma)))
# Arrange "uni_lemma" and "language" columns to match with order of said columns in "aoas" data frame
  arrange(items_filtered, uni_lemma, language)
# Trim dataframe to include only necessary columns
items_filtered <- items_filtered %>%
  dplyr::select(language, item_definition, uni_lemma)
# Clean up discrepancies and match up with AoA data in a CSV file. Then read into a new data frame
uni_lemma_aoa <- read.csv("cleaned_items_aoa.csv")
# Pivot into a wider dataframe which includes language, uni_lemma, and cleaned surface forms
cleaned_uni_lemma_wider <- uni_lemma_aoa %>%
  dplyr::select(language, uni_lemma, aoa, cleaned_item) %>%
  pivot_wider(names_from = language, values_from = cleaned_item, uni_lemma)
# Select for just European Spanish and Mexican Spanish
Eu_Mx <- cleaned_uni_lemma_wider %>% 
  dplyr::select(uni_lemma, `Spanish (European)`, `Spanish (Mexican)`)
#Join in AoA data
clean_aoas_Eu_Mx <- clean_aoas %>%
  dplyr::select(uni_lemma, `Spanish (European)`, `Spanish (Mexican)`) %>%
  left_join(Eu_Mx, join_by(uni_lemma))
# Filter for identical surface forms and calculate correlation
Eu_Mx_Identical <- clean_aoas_Eu_Mx %>% 
  filter(`Spanish (European).y` == `Spanish (Mexican).y`)
# Calculate correlation of AoA values
correlation_Eu_Mx_Identical <- cor(Eu_Mx_Identical$`Spanish (European).x`, Eu_Mx_Identical$`Spanish (Mexican).x`)
# Filter for non-identical surface forms
Eu_Mx_NonIdentical <- clean_aoas_Eu_Mx %>%
  filter(`Spanish (European).y` != `Spanish (Mexican).y`)
# Calculate correlation
correlation_Eu_Mx_NonIdentical <- cor(Eu_Mx_NonIdentical$`Spanish (European).x`, Eu_Mx_NonIdentical$`Spanish (Mexican).x`)
```

t-test for AoA differences
```{r}
# Generate a data frame with all AoA and uni_lemma values
uni_lemma_aoa_wider <- uni_lemma_aoa %>%
  dplyr::select(language, uni_lemma, aoa, cleaned_item) %>%
  pivot_wider(names_from = language, values_from = cleaned_item, uni_lemma)
uni_lemma_aoa_master <- clean_aoas %>%
  left_join(uni_lemma_aoa_wider, join_by(uni_lemma))
# Filter to identical and nonidentical surface forms for all dialect pairings
Eu_Ar_Identical <- uni_lemma_aoa_master %>%
  filter(`Spanish (European).y` == `Spanish (Argentinian).y`)
Eu_Ar_NonIdentical <- uni_lemma_aoa_master %>%
  filter(`Spanish (European).y` != `Spanish (Argentinian).y`)

Eu_Pr_Identical <- uni_lemma_aoa_master %>% 
  filter(`Spanish (European).y` == `Spanish (Peruvian).y`)
Eu_Pr_NonIdentical <- uni_lemma_aoa_master %>% 
  filter(`Spanish (European).y` != `Spanish (Peruvian).y`)

Mx_Ar_Identical <- uni_lemma_aoa_master %>% 
  filter(`Spanish (Mexican).y` == `Spanish (Argentinian).y`)
Mx_Ar_NonIdentical <- uni_lemma_aoa_master %>% 
  filter(`Spanish (Mexican).y` != `Spanish (Argentinian).y`)

Mx_Pr_Identical <- uni_lemma_aoa_master %>% 
  filter(`Spanish (Mexican).y` == `Spanish (Peruvian).y`)
Mx_Pr_NonIdentical <- uni_lemma_aoa_master %>% 
  filter(`Spanish (Mexican).y` != `Spanish (Peruvian).y`)

Pr_Ar_Identical <- uni_lemma_aoa_master %>% 
  filter(`Spanish (Peruvian).y` == `Spanish (Argentinian).y`)
Pr_Ar_NonIdentical <- uni_lemma_aoa_master %>% 
  filter(`Spanish (Peruvian).y` != `Spanish (Argentinian).y`)
# Generate vector of differences for each possible dialect pairing
dif_Eu_Mx_Identical <- Eu_Mx_Identical$`Spanish (European).x` - Eu_Mx_Identical$`Spanish (Mexican).x`
dif_Eu_Mx_NonIdentical <- Eu_Mx_NonIdentical$`Spanish (European).x` - Eu_Mx_NonIdentical$`Spanish (Mexican).x`

dif_Eu_Ar_Identical <- Eu_Ar_Identical$`Spanish (European).x`- Eu_Ar_Identical$`Spanish (Argentinian).x`
dif_Eu_Ar_NonIdentical <- Eu_Ar_NonIdentical$`Spanish (European).x` - Eu_Ar_NonIdentical$`Spanish (Argentinian).x`

dif_Eu_Pr_Identical <- Eu_Pr_Identical$`Spanish (European).x` - Eu_Pr_Identical$`Spanish (Peruvian).x`
dif_Eu_Pr_NonIdentical <- Eu_Pr_NonIdentical$`Spanish (European).x` - Eu_Pr_NonIdentical$`Spanish (Peruvian).x`

dif_Mx_Ar_Identical <- Mx_Ar_Identical$`Spanish (Mexican).x` - Mx_Ar_Identical$`Spanish (Argentinian).x`
dif_Mx_Ar_NonIdentical <- Mx_Ar_NonIdentical$`Spanish (Mexican).x` - Mx_Ar_NonIdentical$`Spanish (Argentinian).x`

dif_Mx_Pr_Identical <- Mx_Pr_Identical$`Spanish (Mexican).x` - Mx_Pr_Identical$`Spanish (Peruvian).x`
dif_Mx_Pr_NonIdentical <- Mx_Pr_NonIdentical$`Spanish (Mexican).x` - Mx_Pr_NonIdentical$`Spanish (Peruvian).x`

dif_Pr_Ar_Identical <- Pr_Ar_Identical$`Spanish (Peruvian).x` - Pr_Ar_Identical$`Spanish (Argentinian).x`
dif_Pr_Ar_NonIdentical <- Pr_Ar_NonIdentical$`Spanish (Peruvian).x` - Pr_Ar_NonIdentical$`Spanish (Argentinian).x`

# Perform a t-test for these differences and place p-values into a data frame
Identical_NonIdentical_p_values <- data.frame(
  Eu_Mx = c(t.test(dif_Eu_Mx_Identical, dif_Eu_Mx_NonIdentical)$p.value),
  Eu_Ar = c(t.test(dif_Eu_Ar_Identical, dif_Eu_Ar_NonIdentical)$p.value),
  Eu_Pr = c(t.test(dif_Eu_Pr_Identical, dif_Eu_Pr_NonIdentical)$p.value),
  Mx_Ar = c(t.test(dif_Mx_Ar_Identical, dif_Mx_Ar_NonIdentical)$p.value),
  Mx_Pr = c(t.test(dif_Mx_Pr_Identical, dif_Mx_Pr_NonIdentical)$p.value),
  Pr_Ar = c(t.test(dif_Pr_Ar_Identical, dif_Pr_Ar_NonIdentical)$p.value)
)
```

Total Identical/NonIdentical Items
```{r}
 AoA_Totals <- data.frame(
     Eu_Ar_I = c(nrow(Eu_Ar_Identical)),
     Eu_Ar_NI = c(nrow(Eu_Ar_NonIdentical)),
     Eu_Mx_I = c(nrow(Eu_Ar_Identical)),
     Eu_Mx_NI = c(nrow(Eu_Ar_NonIdentical)),
     Eu_Pr_I = c(nrow(Eu_Pr_Identical)),
     Eu_Pr_NI = c(nrow(Eu_Pr_NonIdentical)),
     Mx_Ar_I = c(nrow(Mx_Ar_Identical)),
     Mx_Ar_NI = c(nrow(Mx_Ar_NonIdentical)),
     Mx_Pr_I = c(nrow(Mx_Pr_Identical)),
     Mx_Pr_NI = c(nrow(Mx_Pr_NonIdentical)),
     Pr_Ar_I = c(nrow(Pr_Ar_Identical)),
     Pr_Ar_NI = c(nrow(Pr_Ar_NonIdentical))
 )
```

Permutation testing on dialect correlations
```{r} 
# Function to perform permutation test
perm_test <- function(dialect1, dialect2, n_permutations = 10000) {
  
  # Remove NA values
  valid_indices <- complete.cases(dialect1, dialect2)
  dialect1 <- dialect1[valid_indices]
  dialect2 <- dialect2[valid_indices]
  
  # Check for zero variance
  if (var(dialect1) == 0 | var(dialect2) == 0) {
    return(NA)
  }
  
  observed_correlation <- cor(dialect1, dialect2)
  permuted_correlations <- numeric(n_permutations)
  
  for (i in 1:n_permutations) {
    permuted_dialect2 <- sample(dialect2)
    permuted_correlations[i] <- cor(dialect1, permuted_dialect2)
  }
  
  p_value <- mean(permuted_correlations >= observed_correlation)
  return(p_value)
}

# Extract the data for each dialect
argentinian <- cor_data$`Spanish (Argentinian)`
european <- cor_data$`Spanish (European)`
mexican <- cor_data$`Spanish (Mexican)`
peruvian <- cor_data$`Spanish (Peruvian)`

# Create an empty data frame to store p-values
p_values <- data.frame(
  Dialect1 = character(),
  Dialect2 = character(),
  P_Value = numeric(),
  stringsAsFactors = FALSE
)

# Perform permutation tests for each pair of dialects
dialects <- list(argentinian, european, mexican, peruvian)
dialect_names <- c("Spanish (Argentinian)", "Spanish (European)", "Spanish (Mexican)", "Spanish (Peruvian)")

for (i in 1:length(dialects)) {
  for (j in i:length(dialects)) {
    if (i != j) {
      p_val <- perm_test(dialects[[i]], dialects[[j]])
      p_values <- data.frame(rbind(p_values, data.frame(Dialect1 = dialect_names[i], Dialect2 = dialect_names[j], P_Value = p_val)))
    }
  }
}
```

Regression for all 4 dialects
```{r}
#Load gamlss library
library(gamlss)

#Clean sp_ws data frame to include only necessary rows
 sp_ws_minimal <- dplyr::select(sp_ws, age, production, language)

#Find prop_produced for each dialect
items_filtered <- filter(items, items$item_kind == 'word')

sp_ws_left_join <- left_join(sp_ws_minimal,(items_filtered |> group_by(language) |> summarise(n = n())))

sp_ws_prop_produced <- mutate(sp_ws_left_join, prop_produced = sp_ws_left_join$production / sp_ws_left_join$n)

#Run regression
gam_mod <- gamlss(prop_produced ~ pbm(age, lambda = 10000) * as.factor(language),
                  sigma.formula = ~ pbm(age, lambda = 10000) * as.factor(language),
                  data = sp_ws_prop_produced)
summary(gam_mod)
```
Single dialect regression/visualization
```{r}
# Filter data for each dialect and add a Dialect column
sp_ws_prop_produced_Eu <- sp_ws_prop_produced %>% 
  filter(language == "Spanish (European)") 

sp_ws_prop_produced_Ar <- sp_ws_prop_produced %>% 
  filter(language == "Spanish (Argentinian)") 

sp_ws_prop_produced_Mx <- sp_ws_prop_produced %>% 
  filter(language == "Spanish (Mexican)")

sp_ws_prop_produced_Pr <- sp_ws_prop_produced %>% 
  filter(language == "Spanish (Peruvian)")

# Fit GAMLSS models for each dialect
gam_mod_Eu <- gamlss(prop_produced ~ pb(age, lambda = 10000), 
                     sigma.formula = ~ pb(age, lambda = 10000),
                     data = sp_ws_prop_produced_Eu)

gam_mod_Ar <- gamlss(prop_produced ~ pb(age, lambda = 10000), 
                     sigma.formula = ~ pb(age, lambda = 10000),
                     data = sp_ws_prop_produced_Ar)

gam_mod_Mx <- gamlss(prop_produced ~ pb(age, lambda = 10000), 
                     sigma.formula = ~ pb(age, lambda = 10000),
                     data = sp_ws_prop_produced_Mx)

gam_mod_Pr <- gamlss(prop_produced ~ pb(age, lambda = 10000), 
                     sigma.formula = ~ pb(age, lambda = 10000),
                     data = sp_ws_prop_produced_Pr)

# Create age range and predictions for each dialect
age_range_Eu <- seq(min(sp_ws_prop_produced_Eu$age), max(sp_ws_prop_produced_Eu$age), by = 0.1)
pred_data_Eu <- data.frame(age = age_range_Eu)
pred_data_Eu$predicted <- predict(gam_mod_Eu, newdata = pred_data_Eu, type = "response")
pred_data_Eu$language <- "Spanish (European)"

age_range_Ar <- seq(min(sp_ws_prop_produced_Ar$age), max(sp_ws_prop_produced_Ar$age), by = 0.1)
pred_data_Ar <- data.frame(age = age_range_Ar)
pred_data_Ar$predicted <- predict(gam_mod_Ar, newdata = pred_data_Ar, type = "response")
pred_data_Ar$language <- "Spanish (Argentinian)"

age_range_Mx <- seq(min(sp_ws_prop_produced_Mx$age), max(sp_ws_prop_produced_Mx$age), by = 0.1)
pred_data_Mx <- data.frame(age = age_range_Mx)
pred_data_Mx$predicted <- predict(gam_mod_Mx, newdata = pred_data_Mx, type = "response")
pred_data_Mx$language <- "Spanish (Mexican)"

age_range_Pr <- seq(min(sp_ws_prop_produced_Pr$age), max(sp_ws_prop_produced_Pr$age), by = 0.1)
pred_data_Pr <- data.frame(age = age_range_Pr)
pred_data_Pr$predicted <- predict(gam_mod_Pr, newdata = pred_data_Pr, type = "response")
pred_data_Pr$language <- "Spanish (Peruvian)"

# Combine the observed data and prediction data
combined_data <- bind_rows(
  sp_ws_prop_produced_Eu,
  sp_ws_prop_produced_Ar,
  sp_ws_prop_produced_Mx,
  sp_ws_prop_produced_Pr
)

combined_pred_data <- bind_rows(
  pred_data_Eu,
  pred_data_Ar,
  pred_data_Mx,
  pred_data_Pr
)

# Plot using ggplot2 and facet_wrap
ggplot(combined_data, aes(x = age, y = prop_produced, color = language)) +
  geom_point(alpha = 0.5) +
  geom_jitter(width = .5, alpha = .1) +
  geom_line(data = combined_pred_data, aes(x = age, y = predicted, group = language), size = 1) +
  facet_wrap(~ language) +
  labs(
    title = "Logistic Regression: Age vs. Proportion Produced by Dialect",
    x = "Age",
    y = "Proportion Produced"
  ) +
  theme_minimal() +
  scale_color_manual(values = c("blue", "green", "red", "turquoise"))

# Plot all dialects on the same chart
ggplot(combined_data, aes(x = age, y = prop_produced, color = language)) +
  geom_point(alpha = 0.1) +  # Observed data points
  geom_line(data = combined_pred_data, aes(x = age, y = predicted, group = language, color = language), size = 1.5) +  # Predicted curves
  labs(
    title = "Beta Regression: Age vs. Proportion Produced by Dialect",
    x = "Age",
    y = "Proportion Produced"
  ) +
  theme_minimal() +
  scale_color_manual(values = c("blue", "green", "red", "turquoise")) +
  theme(legend.position = "bottom")  # Place legend at the bottom for better readability
```
IRT 
```{r}
# Download necessary libraries
library(psych) # for general functions
library(glue) # for string gluing
library(mirt) # for IRT models
library(ggrepel) # for plot labels

# devtools::install_github("masurp/ggmirt")
library(ggmirt) # extension for 'mirt' 

knitr::knit_hooks$set(inline = function(x) {
  x <- sprintf("%1.2f", x)
  paste(x, collapse = ", ")
})

# Trim and pivot "d_ws" data frame to fit the necessary specifications for the IRT model. Do this for each of the four dialects."
d_ws_wider_Eu <- d_ws %>%
  filter(item_kind == "word") %>%
  filter(language == "Spanish (European)") %>% 
  dplyr::select(data_id, item_id, produces) %>% 
  pivot_wider(names_from = item_id, values_from = produces) %>% 
  dplyr::select(-data_id) %>% 
  mutate_all(as.integer)

d_ws_wider_Ar <- d_ws %>%
  filter(item_kind == "word") %>%
  filter(language == "Spanish (Argentinian)") %>% 
  dplyr::select(data_id, item_id, produces) %>% 
  pivot_wider(names_from = item_id, values_from = produces) %>% 
  dplyr::select(-data_id) %>% 
  mutate_all(as.integer)

d_ws_wider_Mx <- d_ws %>%
  filter(item_kind == "word") %>%
  filter(language == "Spanish (Mexican)") %>% 
  dplyr::select(data_id, item_id, produces) %>% 
  pivot_wider(names_from = item_id, values_from = produces) %>% 
  dplyr::select(-data_id) %>% 
  mutate_all(as.integer)

d_ws_wider_Pr <- d_ws %>%
  filter(item_kind == "word") %>%
  filter(language == "Spanish (Peruvian)") %>% 
  dplyr::select(data_id, item_id, produces) %>% 
  pivot_wider(names_from = item_id, values_from = produces) %>% 
  dplyr::select(-data_id) %>% 
  mutate_all(as.integer)

# Fit the 2PL models for each dialect
modstr_Eu <- glue("F = 1 - {ncol(d_ws_wider_Eu)}")
fit.2pl.Eu <- mirt(d_ws_wider_Eu,
                model = modstr_Eu,
                itemtype = "2PL",
                verbose = FALSE)

modstr_Ar <- glue("F = 1 - {ncol(d_ws_wider_Ar)}")
fit.2pl.Ar <- mirt(d_ws_wider_Ar,
                model = modstr_Ar,
                itemtype = "2PL",
                verbose = FALSE)

modstr_Mx <- glue("F = 1 - {ncol(d_ws_wider_Mx)}")
fit.2pl.Mx <- mirt(d_ws_wider_Mx,
                model = modstr_Mx,
                itemtype = "2PL",
                verbose = FALSE)

modstr_Pr <- glue("F = 1 - {ncol(d_ws_wider_Pr)}")
fit.2pl.Pr <- mirt(d_ws_wider_Pr,
                model = modstr_Pr,
                itemtype = "2PL",
                verbose = FALSE)

# model coefficients for each dialect
coef.2pl.Eu <- as_tibble(coef(fit.2pl.Eu, simplify = TRUE)$items,
                      rownames = "definition")
coef.2pl.Ar <- as_tibble(coef(fit.2pl.Ar, simplify = TRUE)$items,
                      rownames = "definition")
coef.2pl.Mx <- as_tibble(coef(fit.2pl.Mx, simplify = TRUE)$items,
                      rownames = "definition")
coef.2pl.Pr <- as_tibble(coef(fit.2pl.Pr, simplify = TRUE)$items,
                      rownames = "definition")
```
Comparing identical/nonidentical surface forms
```{r}
# Filter "items" data-set to "words" only, then filter again so that it contains only those uni_lemmas contained in the "aoas" data frame
items_filtered <- filter(items, items$item_kind == 'word') %>% 
  filter(uni_lemma %in% c((clean_aoas$uni_lemma)))
# Arrange "uni_lemma" and "language" columns to match with order of said columns in "aoas" data frame
  arrange(items_filtered, uni_lemma, language)
# Trim dataframe to include only necessary columns
items_filtered <- items_filtered %>%
  dplyr::select(language, item_definition, uni_lemma)
# Clean up discrepancies and match up with AoA data in a CSV file. Then read into a new data frame
uni_lemma_aoa <- read.csv("cleaned_items_aoa.csv")
# Pivot into a wider dataframe which includes language, uni_lemma, and cleaned surface forms
cleaned_uni_lemma_wider <- uni_lemma_aoa %>%
  dplyr::select(language, uni_lemma, aoa, cleaned_item) %>%
  pivot_wider(names_from = language, values_from = cleaned_item, uni_lemma)
# Select for just European Spanish and Mexican Spanish
Eu_Mx <- cleaned_uni_lemma_wider %>% 
  dplyr::select(uni_lemma, `Spanish (European)`, `Spanish (Mexican)`)
#Join in AoA data
clean_aoas_Eu_Mx <- clean_aoas %>%
  dplyr::select(uni_lemma, `Spanish (European)`, `Spanish (Mexican)`) %>%
  left_join(Eu_Mx, join_by(uni_lemma))
# Filter for identical surface forms and calculate correlation
Eu_Mx_Identical <- clean_aoas_Eu_Mx %>% 
  filter(`Spanish (European).y` == `Spanish (Mexican).y`)
# Calculate correlation of AoA values
correlation_Eu_Mx_Identical <- cor(Eu_Mx_Identical$`Spanish (European).x`, Eu_Mx_Identical$`Spanish (Mexican).x`)
# Filter for non-identical surface forms
Eu_Mx_NonIdentical <- clean_aoas_Eu_Mx %>%
  filter(`Spanish (European).y` != `Spanish (Mexican).y`)
# Calculate correlation
correlation_Eu_Mx_NonIdentical <- cor(Eu_Mx_NonIdentical$`Spanish (European).x`, Eu_Mx_NonIdentical$`Spanish (Mexican).x`)
```

ANOVA Tests for Lexical Categories, Pairwise t.tests, Density Curves
```{r}
source("scripts/fit_models.R")

# Make a new "wb_data" that includes all uni_lemmas and fit new AoA data
wb_data_full <- d_ws |> 
  group_by(uni_lemma, language, age) |> 
  summarise(total = n(),
            num_true = sum(produces, na.rm = TRUE))

aoas_full <- fit_aoas(wb_data_full)

# Wrangle data to fit proper form for ANOVA test
lexical_categories <- d_ws %>% 
  distinct(uni_lemma, .keep_all=TRUE) %>% 
  dplyr::select(uni_lemma, lexical_category) %>% 
  arrange(uni_lemma)

aoas_full_wider <- aoas_full %>%
  dplyr::select(uni_lemma, language, aoa) %>% 
  pivot_wider(names_from = language, values_from = aoa) %>% 
  left_join(lexical_categories, by = "uni_lemma")

# Perform ANOVA test for 6 possible dialect combinations

aoas_Eu_Mx <- aoas_full_wider %>% 
  dplyr::select(uni_lemma,"Spanish (European)", "Spanish (Mexican)", lexical_category) %>% 
  na.omit() %>%
  mutate(diff_aoa = `Spanish (European)` - `Spanish (Mexican)`)

aov_aoas_Eu_Mx <- aov(diff_aoa ~ lexical_category, data = aoas_Eu_Mx)
summary(aov_aoas_Eu_Mx)

aoas_Eu_Ar <- aoas_full_wider %>% 
  dplyr::select(uni_lemma,"Spanish (European)", "Spanish (Argentinian)", lexical_category) %>% 
  na.omit() %>%
  mutate(diff_aoa = `Spanish (European)` - `Spanish (Argentinian)`)

aov_aoas_Eu_Ar <- aov(diff_aoa ~ lexical_category, data = aoas_Eu_Ar)
summary(aov_aoas_Eu_Ar)

aoas_Eu_Pr <- aoas_full_wider %>% 
  dplyr::select(uni_lemma,"Spanish (European)", "Spanish (Peruvian)", lexical_category) %>% 
  na.omit() %>%
  mutate(diff_aoa = `Spanish (European)` - `Spanish (Peruvian)`)

aov_aoas_Eu_Pr <- aov(diff_aoa ~ lexical_category, data = aoas_Eu_Pr)
summary(aov_aoas_Eu_Pr)

aoas_Ar_Mx <- aoas_full_wider %>% 
  dplyr::select(uni_lemma,"Spanish (Argentinian)", "Spanish (Mexican)", lexical_category) %>% 
  na.omit() %>%
  mutate(diff_aoa = `Spanish (Argentinian)` - `Spanish (Mexican)`)

aov_aoas_Ar_Mx <- aov(diff_aoa ~ lexical_category, data = aoas_Ar_Mx)
summary(aov_aoas_Ar_Mx)

aoas_Ar_Pr <- aoas_full_wider %>% 
  dplyr::select(uni_lemma,"Spanish (Argentinian)", "Spanish (Peruvian)", lexical_category) %>% 
  na.omit() %>%
  mutate(diff_aoa = `Spanish (Argentinian)` - `Spanish (Peruvian)`)

aov_aoas_Ar_Pr <- aov(diff_aoa ~ lexical_category, data = aoas_Ar_Pr)
summary(aov_aoas_Ar_Pr)

aoas_Mx_Pr <- aoas_full_wider %>% 
  dplyr::select(uni_lemma,"Spanish (Mexican)", "Spanish (Peruvian)", lexical_category) %>% 
  na.omit() %>%
  mutate(diff_aoa = `Spanish (Mexican)` - `Spanish (Peruvian)`)

aov_aoas_Mx_Pr <- aov(diff_aoa ~ lexical_category, data = aoas_Mx_Pr)
summary(aov_aoas_Mx_Pr)

# Run pairwise t.tests for each of the possible dialect combinations

pairwise_ttest_Eu_Ar <- pairwise.t.test(aoas_Eu_Ar$diff_aoa, aoas_Eu_Ar$lexical_category, p.adjust.method = "bonferroni")
print(pairwise_ttest_Eu_Ar)

pairwise_ttest_Ar_Mx <- pairwise.t.test(aoas_Ar_Mx$diff_aoa, aoas_Ar_Mx$lexical_category, p.adjust.method = "bonferroni")
print(pairwise_ttest_Ar_Mx)

pairwise_ttest_Ar_Pr <- pairwise.t.test(aoas_Ar_Pr$diff_aoa, aoas_Ar_Pr$lexical_category, p.adjust.method = "bonferroni")
print(pairwise_ttest_Ar_Pr)

pairwise_ttest_Eu_Mx <- pairwise.t.test(aoas_Eu_Mx$diff_aoa, aoas_Eu_Mx$lexical_category, p.adjust.method = "bonferroni")
print(pairwise_ttest_Eu_Mx)

pairwise_ttest_Eu_Pr <- pairwise.t.test(aoas_Eu_Pr$diff_aoa, aoas_Eu_Pr$lexical_category, p.adjust.method = "bonferroni")
print(pairwise_ttest_Eu_Pr)

pairwise_ttest_Mx_Pr <- pairwise.t.test(aoas_Mx_Pr$diff_aoa, aoas_Mx_Pr$lexical_category, p.adjust.method = "bonferroni")
print(pairwise_ttest_Mx_Pr)


# Visualize AoA as density curves
library(reshape2)

# Reshape the data into a long format
aoa_density_curves <- clean_aoas %>% 
  reshape2::melt(id.vars = NULL, 
                measure.vars = c("Spanish (Argentinian)", "Spanish (Mexican)", "Spanish (European)", "Spanish (Peruvian)"),
                variable.name = "Dialect", 
                value.name = "Value")

# Create the density plot
ggplot(aoa_density_curves, aes(x = Value, fill = Dialect, color = Dialect)) +
  geom_density(alpha = 0.5) +
  labs(title = "Age of Acquisition Density Curves for Different Dialects", 
       x = "Age of Acquisition (Months)", 
       y = "Density") +
  scale_x_continuous(limits = c(8, 42)) +  # Adjust x-axis limits
  scale_y_continuous(expand = expansion(mult = c(0, .2))) +  # Adjust y-axis scaling for more vertical space
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

# Create boxplots for each dialect combination
# Create a new column for each data frame containing absolute values
aoas_Ar_Mx$abs_diff_aoa <- abs(aoas_Ar_Mx$diff_aoa)
aoas_Ar_Pr$abs_diff_aoa <- abs(aoas_Ar_Pr$diff_aoa)
aoas_Eu_Ar$abs_diff_aoa <- abs(aoas_Eu_Ar$diff_aoa)
aoas_Eu_Mx$abs_diff_aoa <- abs(aoas_Eu_Mx$diff_aoa)
aoas_Eu_Pr$abs_diff_aoa <- abs(aoas_Eu_Pr$diff_aoa)
aoas_Mx_Pr$abs_diff_aoa <- abs(aoas_Mx_Pr$diff_aoa)

# Create boxplots
ggplot(aoas_Ar_Mx, aes(x = lexical_category, y = abs_diff_aoa, fill = lexical_category)) +
  geom_boxplot() +
  labs(title = "Lexical Categories: Spanish (Argentinian) vs Spanish (Mexican)",
       x = "Lexical Category",
       y = "Diff AoA Ar_Mx") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

ggplot(aoas_Ar_Pr, aes(x = lexical_category, y = abs_diff_aoa, fill = lexical_category)) +
  geom_boxplot() +
  labs(title = "Lexical Categories: Spanish (Argentinian) vs Spanish (Peruvian)",
       x = "Lexical Category",
       y = "Diff AoA Ar_Pr") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

ggplot(aoas_Eu_Ar, aes(x = lexical_category, y = abs_diff_aoa, fill = lexical_category)) +
  geom_boxplot() +
  labs(title = "Lexical Categories: Spanish (European) vs Spanish (Argentinian)",
       x = "Lexical Category",
       y = "Diff AoA Eu_Ar") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

ggplot(aoas_Eu_Mx, aes(x = lexical_category, y = abs_diff_aoa, fill = lexical_category)) +
  geom_boxplot() +
  labs(title = "Lexical Categories: Spanish (European) vs Spanish (Mexican)",
       x = "Lexical Category",
       y = "Diff AoA Eu_Mx") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

ggplot(aoas_Eu_Pr, aes(x = lexical_category, y = abs_diff_aoa, fill = lexical_category)) +
  geom_boxplot() +
  labs(title = "Lexical Categories: Spanish (European) vs Spanish (Peruvian)",
       x = "Lexical Category",
       y = "Diff AoA Eu_Pr") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

ggplot(aoas_Mx_Pr, aes(x = lexical_category, y = abs_diff_aoa, fill = lexical_category)) +
  geom_boxplot() +
  labs(title = "Lexical Categories: Spanish (Mexican) vs Spanish (Peruvian)",
       x = "Lexical Category",
       y = "Diff AoA Mx_Pr") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))
```

IRT 

```{r}
# Download necessary libraries
library(psych) # for general functions
library(glue) # for string gluing
library(mirt) # for IRT models
library(ggrepel) # for plot labels

# devtools::install_github("masurp/ggmirt")
library(ggmirt) # extension for 'mirt' 

knitr::knit_hooks$set(inline = function(x) {
  x <- sprintf("%1.2f", x)
  paste(x, collapse = ", ")
})

# Trim and pivot "d_ws" data frame to fit the necessary specifications for the IRT model. Do this for each of the four dialects."
d_ws_wider_Eu <- d_ws %>%
  filter(item_kind == "word") %>%
  filter(language == "Spanish (European)") %>% 
  dplyr::select(data_id, item_id, produces) %>% 
  pivot_wider(names_from = item_id, values_from = produces) %>% 
  dplyr::select(-data_id) %>% 
  mutate_all(as.integer)

d_ws_wider_Ar <- d_ws %>%
  filter(item_kind == "word") %>%
  filter(language == "Spanish (Argentinian)") %>% 
  dplyr::select(data_id, item_id, produces) %>% 
  pivot_wider(names_from = item_id, values_from = produces) %>% 
  dplyr::select(-data_id) %>% 
  mutate_all(as.integer)

d_ws_wider_Mx <- d_ws %>%
  filter(item_kind == "word") %>%
  filter(language == "Spanish (Mexican)") %>% 
  dplyr::select(data_id, item_id, produces) %>% 
  pivot_wider(names_from = item_id, values_from = produces) %>% 
  dplyr::select(-data_id) %>% 
  mutate_all(as.integer)

d_ws_wider_Pr <- d_ws %>%
  filter(item_kind == "word") %>%
  filter(language == "Spanish (Peruvian)") %>% 
  dplyr::select(data_id, item_id, produces) %>% 
  pivot_wider(names_from = item_id, values_from = produces) %>% 
  dplyr::select(-data_id) %>% 
  mutate_all(as.integer)

# Fit the 2PL models for each dialect
modstr_Eu <- glue("F = 1 - {ncol(d_ws_wider_Eu)}")
fit.2pl.Eu <- mirt(d_ws_wider_Eu,
                model = modstr_Eu,
                itemtype = "2PL",
                verbose = FALSE)

modstr_Ar <- glue("F = 1 - {ncol(d_ws_wider_Ar)}")
fit.2pl.Ar <- mirt(d_ws_wider_Ar,
                model = modstr_Ar,
                itemtype = "2PL",
                verbose = FALSE)

modstr_Mx <- glue("F = 1 - {ncol(d_ws_wider_Mx)}")
fit.2pl.Mx <- mirt(d_ws_wider_Mx,
                model = modstr_Mx,
                itemtype = "2PL",
                verbose = FALSE)

modstr_Pr <- glue("F = 1 - {ncol(d_ws_wider_Pr)}")
fit.2pl.Pr <- mirt(d_ws_wider_Pr,
                model = modstr_Pr,
                itemtype = "2PL",
                verbose = FALSE)

# model coefficients for each dialect
coef.2pl.Eu <- as_tibble(coef(fit.2pl.Eu, simplify = TRUE)$items,
                      rownames = "definition")
coef.2pl.Ar <- as_tibble(coef(fit.2pl.Ar, simplify = TRUE)$items,
                      rownames = "definition")
coef.2pl.Mx <- as_tibble(coef(fit.2pl.Mx, simplify = TRUE)$items,
                      rownames = "definition")
coef.2pl.Pr <- as_tibble(coef(fit.2pl.Pr, simplify = TRUE)$items,
                      rownames = "definition")


```

Plot
```{r}
testInfoPlot(fit.2pl.Ar) + scale_color_manual(values=rep("black", 1000))
```


``` {r}
# Trim and pivot "d_ws" data frame to fit the necessary specifications for the IRT model"
d_ws_wider <- d_ws %>% 
  dplyr::select(data_id, item_id, produces) %>% 
  pivot_wider(names_from = item_id, values_from = produces) %>% 
  dplyr::select(-data_id) %>% 
  mutate_all(as.integer)
```

DIF
``` {r}
# Load necessary 
library(psych) # for general functions
library(glue) # for string gluing
library(mirt) # for IRT models
library(ggrepel) # for plot labels

d_ws_DIF_Eu_Mx <- d_ws %>%
  filter(uni_lemma %in% intersection) %>% 
  dplyr::select(data_id, uni_lemma, language, produces) %>%
  pivot_wider(names_from = uni_lemma, values_from = produces) %>% 
  filter(language == c("Spanish (European)", "Spanish (Mexican)")) %>% 
  dplyr::select(-data_id) %>%  
  mutate_at(vars(-language), ~ as.numeric(.))

DIF_language <- d_ws_DIF_Eu_Mx %>% 
  dplyr::pull(language) %>% 
  as.factor()

d_ws_DIF_Eu_Mx <- d_ws_DIF_Eu_Mx %>% 
  dplyr::select(-language)
  
# Perform DIF Analysis
modstr <- glue("F = 1 - {ncol(d_ws_DIF_Eu_Mx)}")
fit_2plmulti <- multipleGroup(d_ws_DIF_Eu_Mx,
                              model = modstr,
                              method = "MHRM",
                              itemtype = "2PL",
                              group = DIF_language,
                              verbose = FALSE)
  
# Likelihood ratio test
dif_2plmulti <- DIF(fit_2plmulti,
                    which.par = c("a1", "d"),
                    method = "MHRM",
                    technical = list("NCYCLES" = 1000),
                    items2test = c(1:ncol(d_ws_DIF_Eu_Mx)),
                    p.adjust = "BH")
```

t-test for AoA differences
```{r}
# Generate a data frame with all AoA and uni_lemma values
uni_lemma_aoa_wider <- uni_lemma_aoa %>%
  dplyr::select(language, uni_lemma, aoa, cleaned_item) %>%
  pivot_wider(names_from = language, values_from = cleaned_item, uni_lemma)
uni_lemma_aoa_master <- clean_aoas %>%
  left_join(uni_lemma_aoa_wider, join_by(uni_lemma))
# Filter to identical and nonidentical surface forms for all dialect pairings
Eu_Ar_Identical <- uni_lemma_aoa_master %>%
  filter(`Spanish (European).y` == `Spanish (Argentinian).y`)
Eu_Ar_NonIdentical <- uni_lemma_aoa_master %>%
  filter(`Spanish (European).y` != `Spanish (Argentinian).y`)

Eu_Pr_Identical <- uni_lemma_aoa_master %>% 
  filter(`Spanish (European).y` == `Spanish (Peruvian).y`)
Eu_Pr_NonIdentical <- uni_lemma_aoa_master %>% 
  filter(`Spanish (European).y` != `Spanish (Peruvian).y`)

Mx_Ar_Identical <- uni_lemma_aoa_master %>% 
  filter(`Spanish (Mexican).y` == `Spanish (Argentinian).y`)
Mx_Ar_NonIdentical <- uni_lemma_aoa_master %>% 
  filter(`Spanish (Mexican).y` != `Spanish (Argentinian).y`)

Mx_Pr_Identical <- uni_lemma_aoa_master %>% 
  filter(`Spanish (Mexican).y` == `Spanish (Peruvian).y`)
Mx_Pr_NonIdentical <- uni_lemma_aoa_master %>% 
  filter(`Spanish (Mexican).y` != `Spanish (Peruvian).y`)

Pr_Ar_Identical <- uni_lemma_aoa_master %>% 
  filter(`Spanish (Peruvian).y` == `Spanish (Argentinian).y`)
Pr_Ar_NonIdentical <- uni_lemma_aoa_master %>% 
  filter(`Spanish (Peruvian).y` != `Spanish (Argentinian).y`)
# Generate vector of differences for each possible dialect pairing
dif_Eu_Mx_Identical <- Eu_Mx_Identical$`Spanish (European).x` - Eu_Mx_Identical$`Spanish (Mexican).x`
dif_Eu_Mx_NonIdentical <- Eu_Mx_NonIdentical$`Spanish (European).x` - Eu_Mx_NonIdentical$`Spanish (Mexican).x`

dif_Eu_Ar_Identical <- Eu_Ar_Identical$`Spanish (European).x`- Eu_Ar_Identical$`Spanish (Argentinian).x`
dif_Eu_Ar_NonIdentical <- Eu_Ar_NonIdentical$`Spanish (European).x` - Eu_Ar_NonIdentical$`Spanish (Argentinian).x`

dif_Eu_Pr_Identical <- Eu_Pr_Identical$`Spanish (European).x` - Eu_Pr_Identical$`Spanish (Peruvian).x`
dif_Eu_Pr_NonIdentical <- Eu_Pr_NonIdentical$`Spanish (European).x` - Eu_Pr_NonIdentical$`Spanish (Peruvian).x`

dif_Mx_Ar_Identical <- Mx_Ar_Identical$`Spanish (Mexican).x` - Mx_Ar_Identical$`Spanish (Argentinian).x`
dif_Mx_Ar_NonIdentical <- Mx_Ar_NonIdentical$`Spanish (Mexican).x` - Mx_Ar_NonIdentical$`Spanish (Argentinian).x`

dif_Mx_Pr_Identical <- Mx_Pr_Identical$`Spanish (Mexican).x` - Mx_Pr_Identical$`Spanish (Peruvian).x`
dif_Mx_Pr_NonIdentical <- Mx_Pr_NonIdentical$`Spanish (Mexican).x` - Mx_Pr_NonIdentical$`Spanish (Peruvian).x`

dif_Pr_Ar_Identical <- Pr_Ar_Identical$`Spanish (Peruvian).x` - Pr_Ar_Identical$`Spanish (Argentinian).x`
dif_Pr_Ar_NonIdentical <- Pr_Ar_NonIdentical$`Spanish (Peruvian).x` - Pr_Ar_NonIdentical$`Spanish (Argentinian).x`

# Combine the data into a single data frame for easier plotting
differences <- data.frame(
  Group = c(rep("Identical", length(dif_Eu_Mx_Identical)), rep("NonIdentical", length(dif_Eu_Mx_NonIdentical))),
  Difference = c(dif_Eu_Mx_Identical, dif_Eu_Mx_NonIdentical)
)

# Create the boxplots
library(ggplot2)

ggplot(differences, aes(x = Group, y = Difference, fill = Group)) +
  geom_boxplot() +
  labs(title = "Comparing Identical vs Non-Identical Surface Forms (European vs Mexican)",
       x = "Group",
       y = "Diff AoA") +
  theme_minimal()

# Perform a t-test for these differences and place p-values into a data frame
Identical_NonIdentical_p_values <- data.frame(
  Eu_Mx = c(t.test(dif_Eu_Mx_Identical, dif_Eu_Mx_NonIdentical)$p.value),
  Eu_Ar = c(t.test(dif_Eu_Ar_Identical, dif_Eu_Ar_NonIdentical)$p.value),
  Eu_Pr = c(t.test(dif_Eu_Pr_Identical, dif_Eu_Pr_NonIdentical)$p.value),
  Mx_Ar = c(t.test(dif_Mx_Ar_Identical, dif_Mx_Ar_NonIdentical)$p.value),
  Mx_Pr = c(t.test(dif_Mx_Pr_Identical, dif_Mx_Pr_NonIdentical)$p.value),
  Pr_Ar = c(t.test(dif_Pr_Ar_Identical, dif_Pr_Ar_NonIdentical)$p.value)
)
```
Dendrogram
```{r}
correlation_dendrogram <- correl %>% 
  dist() %>% 
  hclust()

  plot(correlation_dendrogram, ylim = c(0, max(correlation_dendrogram$height) * 1.2))
```

## Step 1: Fit the 2PL model
We fit the 2PL model using `itemtype = "2PL"`.
```{r}
modstr <- glue("F = 1 - {ncol(eu_ws)}")
eu_ws_cleaned <- dplyr::select(eu_ws, comprehension, production)
Spanish_European_2PL <- mirt(eu_ws_cleaned,
                model = modstr,
                itemtype = "2PL",
                verbose = FALSE)

# model coefficients
coef.2pl_Spanish_European_2PL <- as_tibble(coef(Spanish_European_2PL, simplify = TRUE)$items,
                      rownames = "definition")
```

# model coefficients
coef.2pl <- as_tibble(coef(fit.2pl, simplify = TRUE)$items,
                      rownames = "definition")
coef.2pl
```
Now we see that discrimination (`a1`) varies across items.
Higher discriminability estimates indicate that the item has a better ability to tell the difference between different levels of latent ability, as will be made clearer in the ICC plots.

## Step 2: Plot the Item Characteristic Curves 
```{r}
tracePlot(fit.2pl, facet = FALSE, legend = TRUE) + scale_color_brewer(palette = "Set3")
```
  
Unlike the ICCs for the 1PL model, the ICCs for the 2PL model do not all have the same shape. 
Item curves which are more "spread out" indicate lower discriminability (i.e., that individuals of a range of ability levels have some probability of getting the item correct). 
Compare this to an item with high discriminability (steep slope): for this item, we have a better estimate of the individual's latent ability based on whether they got the question right or wrong.

*A note about difficulty*: Because of the differing slopes, the rank-order of item difficulty changes across different latent ability levels. We can see that item 9 is generally still the most difficult item (i.e. lowest probability of getting correct for most latent trait values, up until about $\theta=3.3$). 

## Step 3: Plot the Item Information Curves for Items & Test
```{r}
itemInfoPlot(fit.2pl, facet = FALSE, legend = TRUE) + scale_color_brewer(palette = "Set3")
```
  
The item IICs demonstrate that some items provide more information about latent ability for different ability levels. 
The higher the item discriminability estimate, the more information an item provides about ability levels around the point where there is a .5 probability of getting the item right (i.e. the steepest point in the ICC slope). 
For example, item 9 clearly provides the most information at high ability levels, around $\theta=2.3$, but almost no information about low ability levels (< -1) because the item is already too hard for those participants. 
In contrast, item 8, which has low discriminability, doesn't give very much information overall, but covers a wide range of ability levels.

Next, we plot the item information curve for the whole test. This is the sum of all the item IICs above.
```{r}
testInfoPlot(fit.2pl)
```
  
The IIC for the whole test shows that the test provides the most information for slightly-higher-than average ability levels (about $\theta=1$), but does not provide much information about extremely high or low ability levels.

## Step 4: Examine fit of the 2PL model
Next, we test how well the 2PL model fits the data.
```{r}
itemfit(fit.2pl, na.rm = TRUE)
```

Item 9 is still problematic, as before.

## Step 5: Estimate & Plot Ability Scores 

Estimate the individual latent ability scores using the `factor.scores()` function in the `ltm` library.
```{r}
theta.2pl <- fscores(fit.2pl)

head(theta.2pl, 10)
```

Plot the density curve of the estimated ability scores
```{r}
personDist(fit.2pl)
```