add validation analysis

.gitignore

@@ -182,3 +182,6 @@ firestore.indexes.json
 
 # Typedoc generated docs
 docs
+
+validation/data
+validation/simulation/resp

validation/analysis/jsCAT-validation.Rmd

@@ -0,0 +1,359 @@
+---
+title: "jsCAT Validation"
+author: "Wanjing Anya Ma"
+date:  "`r Sys.Date()`"
+output:
+  bookdown::html_document2:
+    toc: true
+    toc_depth: 4
+    theme: cosmo
+    highlight: tango
+
+---
+
+# Data Preprocessing
+## Load Packages
+Let's first load the packages that we need for this chapter. 
+
+```{r, message=FALSE, warning = FALSE, echo = FALSE}
+library("knitr") # for rendering the RMarkdown file
+library("tidyverse") # for data wrangling 
+library(dplyr)
+library(mirt)
+library('ggpubr')
+library(Metrics)
+library(hash)
+library(ggplot2)
+library(GGally)
+library(gridExtra)
+library("catR")
+library(corrplot)
+library(patchwork)
+```
+
+## Settings 
+```{r echo = FALSE}
+# sets how code looks in knitted document
+opts_chunk$set(comment = "")
+
+# suppresses warning about grouping 
+options(dplyr.summarise.inform = F)
+```
+
+## Data
+
+# Theoretical Simulation
+## create resp
+```{r}
+##we'll start small just so that you can see everything, but you'll want to make this bigger downstream.
+## number of items
+ni <-200
+## number of participants for validation
+np <-100
+## number of participants for model calibration
+np.calibration <- 1000
+##now we're going to simulate data according to this model and examine some key properties
+set.seed(1)
+
+th<-rnorm(np.calibration, mean = 0, sd = 2)
+```
+
+```{r}
+################################################################
+##now we have to put this stuff together. what we want is a probability of a correct response for a person to an item
+##we're going to use what you may know from logistic regression
+
+inv_logit<-function(x) {
+  return (exp(x)/(1+exp(x)))
+}
+```
+
+randomly sample b params from a uniform distribution
+```{r}
+b.mat.random <- runif(ni, min = -4, max = 4)
+# item.bank.random <- item.bank %>%
+#    mutate(b = b.mat.random)
+```
+
+```{r}
+func.create.response <- function(b, th, np, ni) {
+  th.mat<-matrix(th,np,ni,byrow=FALSE) #these are the true abilities. we don't observe them,    which will be a real problem for us downstream. but we'll not worry about that today.
+#th.mat #abilities, on per row
+##now the item level part.
+  b.mat <- matrix(rep(b, np), nrow = np, byrow = TRUE)
+
+  pr<-inv_logit(th.mat-b.mat)
+
+  resp <- pr
+
+  for (i in 1:ncol(resp)) {
+    resp[,i]<-rbinom(nrow(resp),1,resp[,i])
+  }
+
+  return (data.frame(resp))
+}
+```
+
+```{r}
+resp.calibration<- func.create.response(b.mat.random, th, np.calibration, ni)
+```
+
+```{r}
+mod <- mirt(data.frame(resp.calibration), 1, itemtype="Rasch")
+
+params.b <- data.frame(coef(mod, IRTpars = TRUE, simplify = TRUE))$items.b
+```
+
+```{r}
+th.sample <- rnorm(np, mean = 0, sd = 2)
+resp.sample <- func.create.response(params.b, th.sample, np, ni)
+```
+```{r}
+resp.sample.with.pid <- resp.sample %>% 
+  add_column(pid = seq(1, 100, 1)) %>% 
+  mutate(pid = paste0("p", pid, "")) %>% 
+  relocate(pid)
+```
+
+```{r}
+#write.csv(resp.sample.with.pid, "../src/data/resp_sample_with_pid.csv")
+```
+
+```{r}
+th_SE.sample <- as_tibble(fscores(mod, method = "ML", response.pattern = resp.sample, max_theta = 6, min_theta = -6, theta_lim = c(-6, 6)))
+```
+
+## catR simulation
+### item bank
+```{r}
+item.bank <- data.frame(a = 1, b = params.b, c = 0, d = 1)
+```
+```{r}
+item.bank.with.item.name <- item.bank %>% 
+  mutate(b = round(b, digits = 8)) %>% 
+  add_column(item = seq(1, 200, 1)) %>% 
+  mutate(item = paste0("X", item, ""))
+```
+```{r}
+#write.csv(item.bank.with.item.name, "../src/wordlist/theoretical.item.bank.csv")
+```
+
+
+### simulation
+```{r}
+func.catSim <- function(resp, item.bank, method){
+  ni = length(resp)
+  np = length(resp[[1]])
+  list.thetas <- NULL
+  
+  list.se <- NULL
+  
+  list.pid <- NULL
+  
+  # define cat
+  start <- list(nrItems = 1, theta = 0)
+  test <- list(method = 'ML', itemSelect = method, range = c(-6, 6))
+  stop <- list(rule = 'length',thr = ni)
+  final <- list(method = 'ML',  range = c(-6, 6))
+  for (pid in 1:np){
+    res <- randomCAT(itemBank = item.bank, responses = as.numeric(resp[pid,]), start = start, test = test, final = final, stop = stop)
+    list.pid <- c(list.pid, rep(pid, ni))
+    list.thetas <- c(list.thetas, res$thetaProv)
+    list.se <- c(list.se, res$seProv)
+  }
+  list.trialNumBlock <- rep(1:ni, np)
+  length(list.trialNumBlock)
+ 
+  return(data.frame(pid = list.pid, trialNumTotal = list.trialNumBlock, F1 = list.thetas, SE_F1 = list.se))
+
+}
+```
+```{r}
+df.mfi <- func.catSim(resp.sample, item.bank, 'MFI')
+df.random <- func.catSim(resp.sample, item.bank, 'random')
+```
+
+
+## jsCAT simulation
+```{r}
+jsCAT.simulation.df <- read.csv("../data/jsCAT-simulation-data.csv")
+```
+
+### Validate theta estimate and standard error of measurement 
+```{r}
+jsCAT.simulation.df.final.estimate <- jsCAT.simulation.df %>% 
+  filter(trialNumTotal == 200, variant == "mfi") %>% 
+    mutate(F1 = thetaEstimate, SE_F1 = thetaSE) %>% 
+    select(pid, trialNumTotal, F1, SE_F1)
+
+g.1 <- ggplot(mapping = aes(x = th_SE.sample$F1,
+                     y = jsCAT.simulation.df.final.estimate$F1)) + 
+  geom_point() +
+  labs(title = "Theta estimate", 
+       x = "Theta estimate from mirt", 
+       y = "Theta estimate from jsCAT") +
+  geom_abline(color = 'gray') + 
+  xlim(-6, 6) + 
+  ylim(-6, 6) + 
+  coord_equal()
+
+g.2 <- ggplot(mapping = aes(x = th_SE.sample$SE_F1,
+                     y = jsCAT.simulation.df.final.estimate$SE_F1)) + 
+  geom_point() +
+   labs(title = "Theta standard error", 
+       x = "Standard error from mirt", 
+       y = "Standard error from jsCAT") +
+  geom_abline(color = 'gray') + 
+  coord_equal()
+
+g.1 + g.2
+
+ggsave("../plots/jsCAT_validation_1.png")
+ggsave("../plots/jsCAT_validation_1.pdf")
+```
+
+### Validate mfi effeciency 
+```{r}
+df.software.compare <- jsCAT.simulation.df %>% 
+  group_by(variant, trialNumTotal) %>% 
+  summarise(sem = mean(thetaSE)) %>% 
+  add_column(software = "jsCAT") %>% 
+  rbind(df.mfi %>% 
+  add_column(variant = "mfi") %>% 
+  rbind(df.random %>% add_column(variant = "random")) %>% 
+  group_by(variant, trialNumTotal) %>% 
+  summarise(sem = mean(SE_F1)) %>% 
+  add_column(software = "catR"))
+  
+```
+
+```{r}
+g.software.compare.1 <- ggplot(df.software.compare, mapping = aes(x = trialNumTotal,
+                     y = sem, 
+                     color = software, 
+                     shape = variant)) + 
+  geom_point() + 
+  ylim(0.2, 1) + 
+  labs(x = "Number of test items",
+       y = "Standard error of measurement",
+       title = "jsCAT vs. catR (Rasch model)") + 
+  theme(legend.position = "none")
+```
+
+# Post-hoc Simulation
+
+## Data Loading
+```{r}
+load("../data/state_c1_2007_7_responses.Rdata")
+```
+
+```{r}
+#grade 7, 
+df.mc <- df %>% 
+  filter(substring(item, 1, 2) == "mc") 
+
+# remove participants who have missing response
+df.mc.keep.pid <- df.mc %>% 
+  filter(!is.na(resp)) %>% 
+  group_by(id) %>% 
+  tally() %>% 
+  filter(n == 55)
+
+df.reading.resp.with.pid <- df.mc %>% 
+  left_join(df.mc.keep.pid, by = "id") %>% 
+  filter(!is.na(n)) %>% 
+  select(-n) %>% 
+  pivot_wider(names_from = item, values_from = resp) %>% 
+  dplyr :: rename(pid = id)
+
+df.reading.resp <- df.reading.resp.with.pid %>% 
+  select(-pid)
+```
+
+
+## Fit 3PL 
+```{r}
+mod.2 <- mirt(data.frame(df.reading.resp), 1, itemtype="3PL")
+```
+```{r}
+params.2 <- data.frame(coef(mod.2, IRTpars = TRUE, simplify = TRUE))
+```
+
+
+## Create item bank
+```{r}
+posthoc.item.bank <- params.2 %>% 
+  mutate(a = round(items.a, digits = 8), 
+         b = round(items.b, digits = 8), 
+         c = round(items.g, digits = 8), 
+         d = round(items.u, digits = 8)) %>% 
+  select(-c(means, F1, items.a, items.b, items.g, items.u)) %>% 
+  tibble::rownames_to_column(var = "item")
+```
+
+```{r}
+#write.csv(posthoc.item.bank, "../src/wordlist/posthoc.item.bank.csv")
+```
+
+## Create resp
+```{r}
+set.seed(123) # Setting seed for reproducibility
+df.reading.resp.with.pid.200 <- df.reading.resp.with.pid %>% 
+  sample_n(200) %>% 
+  mutate(pid = seq(1, 200, 1)) 
+```
+
+```{r}
+#write.csv(df.reading.resp.with.pid.200, "../src/data/resp_posthoc_with_pid.csv")
+```
+
+```{r}
+th_SE.posthoc <- as_tibble(fscores(mod.2, method = "ML", response.pattern = df.reading.resp.with.pid.200 %>% select(-pid), max_theta = 6, min_theta = -6, theta_lim = c(-6, 6)))
+```
+## catR simulation
+```{r}
+df.mfi.posthoc <- func.catSim(df.reading.resp.with.pid.200 %>% select(-pid), posthoc.item.bank %>% select(-item), 'MFI')
+df.random.posthoc <- func.catSim(df.reading.resp.with.pid.200 %>% select(-pid), posthoc.item.bank %>% select(-item), 'random')
+```
+
+## jsCAT simulation
+```{r}
+jsCAT.simulation.posthoc.df <- read.csv("../data/jsCAT-simulation-data-posthoc.csv")
+```
+
+### Validate MFI effeciency 
+```{r}
+df.software.compare.posthoc <- jsCAT.simulation.posthoc.df %>% 
+  group_by(variant, trialNumTotal) %>% 
+  summarise(sem = mean(thetaSE)) %>% 
+  add_column(software = "jsCAT") %>% 
+  rbind(df.mfi.posthoc %>% 
+  add_column(variant = "mfi") %>% 
+  rbind(df.random.posthoc %>% add_column(variant = "random")) %>% 
+  group_by(variant, trialNumTotal) %>% 
+  summarise(sem = mean(SE_F1)) %>% 
+  add_column(software = "catR"))
+  
+```
+
+```{r fig.width=4,fig.height=2}
+g.software.compare.2 <-ggplot(df.software.compare.posthoc, mapping = aes(x = trialNumTotal,
+                     y = sem, 
+                     color = software, 
+                     shape = variant)) + 
+  geom_point() + 
+  ylim(0.2, 1) + 
+  labs(x = "Number of test items",
+       y = "Standard error of measurement",
+       title = "jsCAT vs. catR (3PL model)") + 
+  theme(legend.position = "bottom", 
+        legend.title = element_blank(),
+        legend.key.size = unit(1, "lines")) 
+
+g.software.compare.1 + g.software.compare.2
+
+
+ggsave("../plots/jsCAT_validation_2.png")
+ggsave("../plots/jsCAT_validation_2.pdf")
+```
+