markdly/2pl_expected_score_plot_tam.md

## 2pl_expected_score_plot_tam.md

      
    Raw
  

              2pl_expected_score_plot_tam.md
            
          
    library(tidyverse)
library(TAM)
#> Loading required package: CDM
#> Loading required package: mvtnorm
#> **********************************
#> ** CDM 6.2-91 (2018-05-02 17:29:27)      
#> ** Cognitive Diagnostic Models  **
#> **********************************
#> * TAM 2.12-18 (2018-05-25 11:14:38)
prep_plot_data <- function(params) {
  # to pass tam model information to our expected score function
  params %>% 
    select(item, 
           d1 = AXsi_.Cat1, 
           d2 = AXsi_.Cat2,
           a1 = B.Cat1.Dim1,
           a2 = B.Cat2.Dim1) %>% 
    mutate(expected = pmap(list(d1, d2, a1, a2), expected_score)) %>% 
    select(-d1, -d2, -a1, -a2) %>% 
    unnest()
}

expected_score <- function(d1, d2, a1, a2, theta = seq(-3, 3, by = 0.1)) {
  # TODO: generalise using purrr::accumulate
  t0 <- 1                           # Pr(X = 0) numerator
  t1 <- exp(a1 * theta - d1)        # Pr(X = 1) numerator
  t2 <- exp(a2 * theta - (d1 + d2)) # Pr(X = 2) numerator
  D <- t0 + t1 + t2                 # Pr(X = x) denominator
  E <- 0 * t0/D + a1 * t1/D + a2 * t2/D  # Expected score
  tibble(theta = theta, exp_scr = E)
}

expected_score_plot <- function(gg) {
  ggplot(gg, aes(x = theta, y = exp_scr, colour = item,
                 xmin = -3, xmax = 3, ymin = 0)) + 
    geom_line() +
    scale_x_continuous(breaks = seq(-3, 3))
}


data(data.timssAusTwn.scored)
df <- data.timssAusTwn.scored %>% select(M032166:M032673)

### PCM
# This is our 'control case' to show our custom plot works similar to plot.tam at least sometimes
mod1 <- tam.mml(df, verbose = FALSE)
plot(mod1, items=3, export=FALSE, observed=FALSE)
#> Iteration in WLE/MLE estimation  1   | Maximal change  0.1628 
#> Iteration in WLE/MLE estimation  2   | Maximal change  0.0247 
#> Iteration in WLE/MLE estimation  3   | Maximal change  0.0024 
#> Iteration in WLE/MLE estimation  4   | Maximal change  3e-04 
#> Iteration in WLE/MLE estimation  5   | Maximal change  0 
#> ----
#>  WLE Reliability= 0.773

mod1$item %>% 
  filter(item == "M032757") %>% 
  prep_plot_data() %>% 
  expected_score_plot()

### 2PL
# plots look different when using the 2PL model though
mod2 <- tam.mml.2pl(df, verbose = FALSE)

# plot has max of 2 when plotted using plot.tam
plot(mod2, items=3, export=FALSE, observed=FALSE)
#> Iteration in WLE/MLE estimation  1   | Maximal change  2.9614 
#> Iteration in WLE/MLE estimation  2   | Maximal change  2.8004 
#> Iteration in WLE/MLE estimation  3   | Maximal change  2.04 
#> Iteration in WLE/MLE estimation  4   | Maximal change  2.4486 
#> Iteration in WLE/MLE estimation  5   | Maximal change  1.1065 
#> Iteration in WLE/MLE estimation  6   | Maximal change  0.5179 
#> Iteration in WLE/MLE estimation  7   | Maximal change  0.2052 
#> Iteration in WLE/MLE estimation  8   | Maximal change  0.0079 
#> Iteration in WLE/MLE estimation  9   | Maximal change  0.0011 
#> Iteration in WLE/MLE estimation  10   | Maximal change  2e-04 
#> Iteration in WLE/MLE estimation  11   | Maximal change  0 
#> ----
#>  WLE Reliability= 0.771

# our expected score chart has max of above 2
mod2$item %>% 
  filter(item == "M032757") %>% 
  prep_plot_data() %>% 
  expected_score_plot()

mod2$item %>% filter(item == "M032757")
#>      item    N        M   xsi.item AXsi_.Cat1 AXsi_.Cat2 B.Cat1.Dim1
#> 1 M032757 1769 1.413793 -0.8237062   1.306939  -1.647412   0.7208611
#>   B.Cat2.Dim1
#> 1    2.400068
Created on 2018-06-12 by the reprex package (v0.2.0).