model_fit.Rmd

---
title: "Model fit"
author: "Piotr Król"
date: "23 08 2021"
output: html_document
---

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

Content:

* 670 items model
  + Items fit
  + Local dependence
* Cramer's V for responses
* 667 items model

```{r libraries, warning=FALSE, message=FALSE}
library(mirt)
library(DescTools) #For classical Cramer's V calculation
```

Load already prepared data (data_preparation.R) :

```{r data loading, eval=F}
responses_demo <- read.csv("Data/responses_demo.csv", encoding = "UTF-8")
responses <- as.matrix(responses_demo[,5:ncol(responses_demo)])
cdi <- read.csv("Data/cdi.csv", encoding = "UTF-8")
```

```{r ready report data, include=F}
load("Data/model_fit.RData")
```

# 670 items model

Table of log-likelihoods and quadratures (cycles were increased properly to reach convergence):

| Quadratures | LogLik      |
|-------------|-------------|
| 15          | -472591.2   | 
| 61          | -463630.5   | 
| 151         | -462896.9   | 
| 301         | -462860.8   | 
| 451         | -462860.681 | 
| 601         | -462860.656 | 

Based on 3rd report - "Quadratures" from 28.06.2021

601 quadrature points seems to be an optimal number - differece in LogLik is really marginal.

Create model:

```{r eval=F}
mod <- mirt(data = responses, model = 1, SE = TRUE, quadpts = 601, technical = list(NCYCLES = 10000))
save(mod, file = "Data/mod")
```

```{r}
mod
```

## Items fit

```{r eval=F}
itemfit <- itemfit(mod, method = "MAP") #Use default S_X2 statistic and MAP ability estimation
```

```{r}
head(itemfit)
```

Misfits for p = 0.01:

```{r}
length(which(itemfit$p.S_X2 < 0.01))
```
```{r}
itemfit[itemfit$p.S_X2 < 0.01, ]
```

```{r}
cdi[row.names(itemfit[itemfit$p.S_X2 < 0.01, ]), ]
```

Misfits for p = 0.001:

```{r}
length(which(itemfit$p.S_X2 < 0.001))
```
```{r}
itemfit[itemfit$p.S_X2 < 0.001, ]
```

```{r}
cdi[row.names(itemfit[itemfit$p.S_X2 < 0.001, ]), ]
```

```{r}
hist(itemfit$RMSEA.S_X2)
```

## Local dependence

### Residuals 

"By default, residuals() computes the local dependence (LD)
pairwise statistic between each pair of items, which is very similar to a signed 2 value (Chen
and Thissen 1997). Also, a standardized version of the LD statistic (Cramer's V) is printed
above the diagonal to aid in interpretation when items contain more than two response options
and hence more degrees of freedom" ~ mirt documentation

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

```{r}
residuals_up <- residuals
residuals_down <- residuals
residuals_up[lower.tri(residuals_up)] <- NA
residuals_down[upper.tri(residuals_down)] <- NA
```

```{r}
hist(residuals_up, main="Histogram of signed Cramer's V coefficients", xlab= "Cramer's V")
```
```{r}
summary(as.vector(residuals_up))
```
Check positions with the highest Cramer's V:

```{r}
indexes <- arrayInd(which(residuals_up > 0.3), dim(residuals))
data.frame(pos1 = cdi[indexes[, 1], "position"], pos2 = cdi[indexes[, 2], "position"], CramersV = residuals_up[indexes])
```
Check positions with the lowest Cramer's V:

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

Check LD statistic:

```{r}
hist(residuals_down, main="Histogram of LD", xlab= "LD")
```
```{r}
summary(as.vector(residuals_down))
```
Check positions with the highest LD:

```{r}
indexes <- arrayInd(which(residuals_down > 200), dim(residuals))
data.frame(pos1 = cdi[indexes[, 1], "position"], pos2 = cdi[indexes[, 2], "position"], Chi.sq = residuals_down[indexes])
```

How many items have LD > 100 with at least one item:

```{r}
indexes <- arrayInd(which(residuals_down > 100), dim(residuals))
rows <- indexes[, 1]
columns <- indexes[, 1]
unique <- unique(c(unique(rows), unique(columns)))
length(unique)
```

Check positions with the lowest LD:

```{r}
indexes <- arrayInd(which(residuals_down < 0.009), dim(residuals))
data.frame(pos1 = cdi[indexes[, 1], "position"], pos2 = cdi[indexes[, 2], "position"], Chi.sq = residuals_down[indexes])
```

# Cramer's V for responses matrix

Check classical Cramer's V values:

```{r eval=F}
cramersV <- PairApply(responses, FUN = CramerV)
cramersV[lower.tri(cramersV)] <- NA
diag(cramersV) <- NA
```

```{r}
hist(cramersV, main="Histogram of Cramer's V values", xlab= "Cramer's V")
```
How many items correlate higher than 0.8 with at least one other item:

```{r}
indexes <- arrayInd(which(cramersV > 0.8), dim(cramersV))
rows <- indexes[, 1]
columns <- indexes[, 1]
unique <- unique(c(unique(rows), unique(columns)))
length(unique)
```

What are these items:

```{r}
indexes <- arrayInd(which(cramersV > 0.8), dim(cramersV))
data.frame(pos1 = cdi[indexes[, 1], "position"], pos2 = cdi[indexes[, 2], "position"], cramersV = cramersV[indexes])
```

# 667 items model

Create model on a restricted set of items (remove items marked as misfits for p < 0.001):

```{r eval=F}
items_to_remove <- as.numeric(row.names(itemfit[itemfit$p.S_X2 < 0.001, ]))
save(items_to_remove, file="Data/items_to_remove")
```

```{r}
items_to_remove
```

```{r eval=F}
mod_r <- mirt(data = responses[, -c(392, 417, 652)], model = 1, SE = TRUE, quadpts = 601, technical = list(NCYCLES = 10000))
save(mod_r, file="Data/mod_r")
```
```{r}
mod_r
```