## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  fig.width = 6,
  fig.height = 4,
  fig.align = "center"
)
library(csemGT)

## ----ex1-subsample------------------------------------------------------------
set.seed(1998)
idx <- sample(nrow(iowa_like), size = 600)
sub <- iowa_like[idx, ]
n_i <- ncol(sub)

fit <- csem_gt(sub, error_type = "absolute", method = "full")

## ----ex1-lord-check-----------------------------------------------------------
score_tbl <- fit$by_score[, c("observed_score", "csem.absolute")]
score_tbl$csem_lord <- with(
  score_tbl,
  sqrt(observed_score * (n_i - observed_score) / (n_i - 1)) / n_i
)
score_tbl$diff <- score_tbl$csem.absolute - score_tbl$csem_lord
head(score_tbl, 8)

max(abs(score_tbl$diff))

## ----ex1-plot, fig.cap = "Per-person absolute CSEM for the iowa_like subsample (n = 600, I = 40), full estimator."----
plot(fit, main = "iowa_like: absolute CSEM by observed score")

## ----ex1-bootstrap------------------------------------------------------------
fit_boot <- csem_gt(
  sub,
  error_type        = "absolute",
  method            = "full",
  bootstrap         = TRUE,
  return_analytical = TRUE,
  R                 = 200L,
  bootstrap_type    = "item",
  ci_method         = "percentile",
  seed              = 1998
)

## ----ex1-bootstrap-table------------------------------------------------------
cols <- c("observed_score", "csem.absolute",
          "se.analytic.absolute",     "se.boot.absolute",
          "ci_low.analytic.absolute", "ci_up.analytic.absolute",
          "ci_low.boot.absolute",     "ci_up.boot.absolute")
fit_boot$estimates[1:5, cols]

## ----ex2-deps-----------------------------------------------------------------
has_mirt <- requireNamespace("mirt", quietly = TRUE)

## ----ex2-key, eval = has_mirt-------------------------------------------------
key <- c(1, 4, 5, 2, 3, 1, 2, 1,
         3, 1, 2, 4, 2, 1, 5, 3,
         4, 4, 1, 4, 3, 3, 4, 1,
         3, 5, 1, 3, 1, 5, 4, 3)

scored <- mirt::key2binary(mirt::SAT12, key)
scored <- scored[stats::complete.cases(scored), ]
dim(scored)

## ----ex2-fit, eval = has_mirt-------------------------------------------------
fit_sat <- csem_gt(
  scored,
  error_type = "relative",
  method     = c("full", "large_a", "uncorrelated"),
  smoother   = "none"
)

## ----ex2-plot, eval = has_mirt, fig.cap = "Relative CSEM under the three estimators implemented by csem_gt() for the SAT12 science test."----
plot(
  fit_sat,
  compare_methods = TRUE,
  compare_points  = TRUE,
  show_smooth     = FALSE
)

## ----ex3-fit------------------------------------------------------------------
fit_ip <- csem_gt(
  ipip_like,
  error_type = "relative",
  method     = "full",
  smoother   = "polynomial"
)

## ----ex3-components-----------------------------------------------------------
fit_ip$variance_components$person
fit_ip$variance_components$item
fit_ip$variance_components$residual
fit_ip$variance_components$reliability_coefficients$erho2

## ----ex3-as-df----------------------------------------------------------------
as.data.frame(fit_ip)[1:5, ]

## ----ex3-plot, fig.cap = "Relative CSEM against the observed score for ipip_like, with the 95% CI band around the quadratic fit."----
plot(fit_ip, plot_type = "both", cibands = "model")

## ----ex3-dstudy---------------------------------------------------------------
fit_05 <- csem_gt(
  ipip_like,
  error_type = "relative", method = "full",
  smoother   = "none",
  n_items_D  = 5L
)
fit_20 <- csem_gt(
  ipip_like,
  error_type = "relative", method = "full",
  smoother   = "none",
  n_items_D  = 20L
)

rho_ref <- fit_ip$variance_components$reliability_coefficients$erho2
spearman_brown <- function(rho, K) (K * rho) / (1 + (K - 1) * rho)

dstudy <- data.frame(
  test_length = c(5L, 10L, 20L),
  erho2_gt    = c(
    fit_05$variance_components$reliability_coefficients$erho2,
    rho_ref,
    fit_20$variance_components$reliability_coefficients$erho2
  ),
  erho2_sb    = c(
    spearman_brown(rho_ref, 0.5),
    rho_ref,
    spearman_brown(rho_ref, 2.0)
  )
)
dstudy$diff <- dstudy$erho2_gt - dstudy$erho2_sb
dstudy