## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")

## ----setup--------------------------------------------------------------------
library(intraclass)

## ----have, include = FALSE----------------------------------------------------
have_psych  <- requireNamespace("psych",  quietly = TRUE)
have_irr    <- requireNamespace("irr",    quietly = TRUE)
have_irr_icc <- requireNamespace("irrICC", quietly = TRUE)

# psych and irr consume a wide subjects-by-raters matrix; the shipped `ratings`
# data are long, so reshape once here.
to_wide <- function(d) {
  w <- reshape(d, idvar = "subject", timevar = "rater", direction = "wide")
  w <- w[order(as.integer(as.character(w$subject))), ]
  as.matrix(w[, -1])
}

## ----validation, echo = FALSE, eval = have_psych && have_irr && requireNamespace("gt", quietly = TRUE)----
wm <- to_wide(ratings)
# Scalar `type` and `unit` make this a one-row fit, so `[1]` is that row --
# the coefficient asked for. Elsewhere, select by `term`.
ic <- function(model, type, unit) {
  tidy(icc(ratings, subject = subject, rater = rater, score = score,
           model = model, type = type, unit = unit))$estimate[1]
}
ps <- psych::ICC(wm)$results
psv <- stats::setNames(ps$ICC, ps$type)

rows <- list(
  c("ICC(1)",   "oneway", "agreement",   "single",  "ICC1"),
  c("ICC(1,k)", "oneway", "agreement",   "average", "ICC1k"),
  c("ICC(A,1)", "twoway", "agreement",   "single",  "ICC2"),
  c("ICC(A,k)", "twoway", "agreement",   "average", "ICC2k"),
  c("ICC(C,1)", "twoway", "consistency", "single",  "ICC3"),
  c("ICC(C,k)", "twoway", "consistency", "average", "ICC3k")
)

comparison <- do.call(rbind, lapply(rows, function(r) {
  data.frame(
    coefficient = r[1],
    intraclass  = ic(r[2], r[3], r[4]),
    psych       = unname(psv[r[5]]),
    irr         = irr::icc(wm, model = r[2], type = r[3], unit = r[4])$value
  )
}))

gt::gt(comparison) |>
  gt::tab_header(title = "The ICC family on the balanced ratings data") |>
  gt::cols_label(
    coefficient = "Coefficient",
    intraclass = "intraclass (REML)",
    psych = "psych (ANOVA)",
    irr = "irr (ANOVA)"
  ) |>
  gt::fmt_number(columns = c(intraclass, psych, irr), decimals = 5)

## ----validation-gap, include = FALSE, eval = exists("comparison")-------------
max_gap <- max(abs(comparison$intraclass - comparison$psych),
               abs(comparison$intraclass - comparison$irr))

## ----irricc, eval = have_irr_icc----------------------------------------------
w <- reshape(ratings, idvar = "subject", timevar = "rater", direction = "wide")
w <- w[order(as.integer(as.character(w$subject))), ]
gwet_frame <- data.frame(
  Target = as.integer(as.character(w$subject)),
  J1 = w$score.1, J2 = w$score.2, J3 = w$score.3, J4 = w$score.4
)
gwet_agree <- irrICC::icc2.inter.fn(gwet_frame)$icc2r
intraclass_a1 <- with(tidy(icc(ratings, subject = subject, rater = rater, score = score,
                              model = "twoway", type = "agreement",
                              unit = "single")), estimate[term == "ICC(A,1)"])

## ----irricc-table, echo = FALSE, eval = have_irr_icc && requireNamespace("gt", quietly = TRUE)----
data.frame(
  source   = c("intraclass ICC(A,1)", "irrICC icc2r (Gwet)"),
  estimate = c(intraclass_a1, gwet_agree)
) |>
  gt::gt() |>
  gt::tab_header(title = "Single-rater absolute agreement on the ratings data") |>
  gt::cols_label(source = "Package and coefficient", estimate = "Estimate") |>
  gt::fmt_number(columns = estimate, decimals = 5)

## ----incomplete-show, echo = FALSE, eval = requireNamespace("gt", quietly = TRUE)----
wide_incomplete <- reshape(ratings_incomplete, idvar = "subject",
                           timevar = "rater", direction = "wide")
wide_incomplete <- wide_incomplete[order(as.integer(as.character(wide_incomplete$subject))), ]
colnames(wide_incomplete) <- c("subject", paste0("rater", 1:4))
gt::gt(wide_incomplete) |>
  gt::tab_header(title = "Scores in ratings_incomplete, one row per subject") |>
  gt::cols_label(
    subject = "Subject",
    rater1 = "Rater 1",
    rater2 = "Rater 2",
    rater3 = "Rater 3",
    rater4 = "Rater 4"
  ) |>
  gt::fmt_number(columns = rater1:rater4, decimals = 0) |>
  gt::sub_missing(missing_text = "not rated")

## ----incomplete-classical, echo = FALSE, eval = requireNamespace("gt", quietly = TRUE)----
wm_inc <- to_wide(ratings_incomplete)
surviving <- sum(stats::complete.cases(wm_inc))
data.frame(
  quantity = c("Observed ratings", "Possible ratings",
               "Subjects left after listwise deletion"),
  value = c(nrow(ratings_incomplete), nrow(ratings), surviving)
) |>
  gt::gt() |>
  gt::tab_header(title = "What listwise deletion leaves of ratings_incomplete") |>
  gt::cols_label(quantity = "Quantity", value = "Count") |>
  gt::fmt_number(columns = value, decimals = 0)

## ----incomplete-intraclass----------------------------------------------------
fit_inc <- icc(ratings_incomplete, subject = subject, rater = rater, score = score,
               model = "twoway", type = "agreement", unit = "average")
gl_inc <- glance(fit_inc)

## ----incomplete-intraclass-table, echo = FALSE, eval = requireNamespace("gt", quietly = TRUE)----
data.frame(
  quantity = c("ICC(A,k) estimate", "Subjects used", "Ratings used",
               "Effective number of ratings (k_eff)"),
  value = c(with(tidy(fit_inc), estimate[term == "ICC(A,k)"]),
            gl_inc$n_subjects, gl_inc$n_obs, gl_inc$k_eff)
) |>
  gt::gt() |>
  gt::tab_header(title = "intraclass on ratings_incomplete") |>
  gt::cols_label(quantity = "Quantity", value = "Value") |>
  gt::fmt_number(columns = value, rows = c(1, 4), decimals = 3) |>
  gt::fmt_number(columns = value, rows = c(2, 3), decimals = 0)