SiZ.Rmd

---
title: "IRMiK: Słowa i Zdania"
author: "Piotr Król"
date: "26 01 2022"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.align = "center")
```

Load libraries and functions:

```{r libraries, warning=FALSE, message=FALSE}
library(mirt)
library(mirtCAT)
library(parallel)
library(plyr)
library(ggpubr)
source(paste0(getwd(),"/Functions/get_sim_results.R"))
source(paste0(getwd(),"/Functions/get_first_item_responses.R"))
source(paste0(getwd(),"/Functions/misfits_removal.R"))
source(paste0(getwd(),"/Functions/prepare_data.R"))
source(paste0(getwd(),"/Functions/prepare_input_files.R"))
source(paste0(getwd(),"/Functions/get_start_thetas.R"))
```

# Data preparation

```{r include=F}
load("Data/SiZ.RData")
```

Load data downloaded from https://drive.google.com/file/d/1i4YupJCfE6qt1iZ9A6iqZaxXv71gZ_JV/view and prepare it with ready function: 

```{r eval=F}
data <- read.csv("Data/irmik_lex.csv", encoding = "UTF-8")
names(data)[names(data) == 'Płeć'] <- 'Plec' #Polish symbols causes some problems later when preparing data
prepare_data(data, "Z")
rm(d)
```

Prepare responses matrix

```{r}
responses <- as.matrix(responses_demo[,5:ncol(responses_demo)])
print(paste0("There are ", nrow(responses), " respondents that answered to ", ncol(responses), " items and there are ", length(which(is.na(responses))), " missing responses."))
```

# Model creation

Load model sent by Grzegorz

```{r eval=F}
load("Data/output_ws.Rds")
mod <- output[[1]]
items_removed <- output[[2]]
```

```{r}
mod
```

## Items fit

```{r}
items_removed
```

```{r}
items_to_remove_ids <- as.numeric(laply(unlist(items_removed), function(x) substr(x, 5, nchar(x))))
cdi[which(cdi$number.ws %in% items_to_remove_ids), ]
```

Create new responses matrix and new cdi df without this item:

```{r}
responses_r <- responses[, -which((colnames(responses) %in% unlist(items_removed)))]
cdi_r <- cdi[-which(cdi$number.ws %in% items_to_remove_ids), ]
if (nrow(cdi_r) == ncol(responses_r)) print(paste0(ncol(responses_r), " items left after removal."))
```

## Local dependence

```{r eval=F}
residuals <- residuals(mod)
```

```{r}
residuals_up <- residuals
residuals_up[lower.tri(residuals_up)] <- NA
hist(residuals_up, main="Histogram of signed Cramer's V coefficients", xlab= "Cramer's V")
```

```{r}
summary(as.vector(residuals_up))
```

The most highly dependent positions:

```{r}
indexes <- arrayInd(which(residuals_up > 0.3), dim(residuals))
data.frame(pos1 = cdi_r[indexes[, 1], "position"], pos2 = cdi_r[indexes[, 2], "position"], CramersV = residuals_up[indexes])
```

But no pair reach the threshold of 0.5 so I don't remove anything.

# Model exploration

## IRT items parameteres

```{r eval=F}
params_IRT <- as.data.frame(coef(mod, simplify = T, IRTpars = T)$items)[1:2]
params_cdi <- cbind(params_IRT, cdi_r)
```

```{r fig.width=12, fig.height=8}
ggplot(params_cdi, aes(b, a, colour=category, label=position)) +
  geom_text(check_overlap = T) +
  geom_point(alpha = 0.3) +
  xlab("Difficulty") +
  ylab("Discrimination") +
  labs(colour = "Category") +
  theme(
  panel.background = element_rect(fill = NA),
  panel.grid.major = element_line(colour = "grey"),
  legend.position = "top"
)
```

## Ability vs. SE

```{r eval=F}
fscores <- as.data.frame(fscores(mod, method = "MAP", full.scores.SE = T))
```

```{r fig.width=12, fig.height=8}
ggplot(fscores, aes(F1, SE_F1)) +
  geom_point(alpha = 0.3) +
  xlab("Ability level (Θ)") +
  ylab("Standard error (SE)") +
  labs(title = "CDI: WS (all children)") +
  theme(
  panel.background = element_rect(fill = NA),
  panel.grid.major = element_line(colour = "lightgrey"),
  text = element_text(size=16)
)
```

```{r}
summary(fscores$SE_F1)
```

## SE vs. age

```{r fig.width=12, fig.height=8 }
ggplot(data.frame(age = as.factor(responses_demo[responses_demo$participant %in% responses_demo$participant, ]$months), SE = fscores$SE_F1), aes(age, SE)) +
  geom_boxplot() +
  xlab("Age (months)") +
  ylab("Standard error of ability estimations") +
  labs(title = "CDI: WS") +
  theme(
  panel.background = element_rect(fill = NA),
  panel.grid.major = element_line(colour = "lightgrey"),
  text = element_text(size=16)
  )
```

# Simulations

```{r eval=F}
params <- as.data.frame(coef(mod, simplify = T)$items)[1:2]
mo <- generate.mirt_object(params, '2PL')
cl <- makeCluster(detectCores())
```

## Whole pool for everyone

```{r eval=F}
results_all <- mirtCAT(mo = mo, method = "MAP", criteria = "MI", start_item = "MI", local_pattern = responses_r, cl = cl, design = list(min_items = ncol(responses_r)))
```

```{r}
get_sim_results(results_all, fscores, params)
```

```{r fig.width=12, fig.height=8}
SE_history <- laply(results_all, function(x) x$thetas_SE_history)
meanSE_item <- as.data.frame(colMeans(SE_history))
meanSE_item$item <- row.names(meanSE_item)
colnames(meanSE_item) <- c("mean_SE", "item")
meanSE_item$item <- as.numeric(meanSE_item$item)

ggline(meanSE_item, x = "item", y = "mean_SE", numeric.x.axis = TRUE) +
  xlab ("Test length") +
  ylab ("Mean SE") +
  labs(title = "CDI: WS") +
  theme(
    panel.background = element_rect(fill = NA),
    panel.grid.major = element_line(colour = "lightgrey"),
    text = element_text(size=16)
  ) + 
  scale_y_continuous(breaks = c(0.1, 0.15, 0.2, 0.25, 0.5, 0.75, 1)) + 
  scale_x_continuous(breaks = seq(0, nrow(params), by=10), limits = c(1, 100))
```

## SE = 0.1 stop criterion

```{r eval=F}
results_SE1 <- mirtCAT(mo = mo, method = "MAP", criteria = "MI", start_item = "MI", local_pattern = responses_r, cl = cl, design = list(min_SEM = 0.1))
```

```{r}
get_sim_results(results_SE1, fscores, params)
```

```{r fig.width=12, fig.height=8}
#Prepare cuts
tests_lengths <- laply(results_SE1, function(x) length(x$items_answered))
cuts <- cut(tests_lengths, breaks = c(0, 10, 15, 25, 35, 50, 75, 100, 665, 666), labels = c("1 - 10", "11 - 15", "16 - 25", "26 - 35", "36 - 50", "51 - 75", "76 - 100", "101 - 665", "666 (all)"))

#Plot
ggplot(data.frame(round(table(cuts) / nrow(responses_r) * 100, 1)), aes(x = cuts, y = Freq)) +
  xlab("Number of items administered") +
  ylab("Percent of respondents (%)") +
  geom_bar(stat = "identity") +
  geom_text(aes(label = Freq), vjust = -0.3, size=6) + 
  labs(title = "CDI: WS. SE = 0.1 stop criterion") +
  theme_pubclean() +
  theme(text = element_text(size=16))
```

## SE = 0.15 stop criterion

```{r eval=F}
results_SE15 <- mirtCAT(mo = mo, method = "MAP", criteria = "MI", start_item = "MI", local_pattern = responses_r, cl = cl, design = list(min_SEM = 0.15))
```

```{r}
get_sim_results(results_SE15, fscores, params)
```

```{r fig.width=12, fig.height=8}
#Prepare cuts
tests_lengths <- laply(results_SE15, function(x) length(x$items_answered))
cuts <- cut(tests_lengths, breaks = c(0, 10, 15, 25, 35, 50, 75, 100, 665, 666), labels = c("1 - 10", "11 - 15", "16 - 25", "26 - 35", "36 - 50", "51 - 75", "76 - 100", "101 - 665", "666 (all)"))

#Plot
ggplot(data.frame(round(table(cuts) / nrow(responses_r) * 100, 1)), aes(x = cuts, y = Freq)) +
  xlab("Number of items administered") +
  ylab("Percent of respondents (%)") +
  geom_bar(stat = "identity") +
  geom_text(aes(label = Freq), vjust = -0.3, size=6) + 
  labs(title = "CDI: WS. SE = 0.15 stop criterion") +
  theme_pubclean() +
  theme(text = element_text(size=16))
```

## Different start thetas

```{r eval=F}
output2 <- get_start_thetas(responses_demo[which(responses_demo$months >= 16), ], fscores)
fscores_aggr <- output2[[1]]
fscores_aggr$F1 <- fscores_aggr$F1 - 2
start_thetas <- output2[[2]]
start_thetas <- start_thetas - 2
```

TODO: Simulations with different start thetas when stop criterion would be specified

# Input files preparation

```{r eval=F}
prepare_input_files(params, cdi_r, fscores_aggr, question = "Czy Twoje dziecko mówi", group = "mowa", file1 = "items-WS.csv", file2 = "startThetas-WS.csv")
```