## ----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, eval = have_psych && have_irr--------------------------------
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
  )
}))

knitr::kable(comparison, digits = 5, row.names = FALSE)

## ----validation-gap, eval = have_psych && have_irr----------------------------
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)"])

data.frame(
  source   = c("intraclass ICC(A,1)", "irrICC icc2r (Gwet)"),
  estimate = c(intraclass_a1, gwet_agree)
)

## ----incomplete-show----------------------------------------------------------
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))
knitr::kable(wide_incomplete, row.names = FALSE)

## ----incomplete-classical-----------------------------------------------------
wm_inc <- to_wide(ratings_incomplete)
surviving <- sum(stats::complete.cases(wm_inc))
c(observed_cells = nrow(ratings_incomplete),
  possible_cells = nrow(ratings),
  subjects_after_listwise_deletion = surviving)

## ----incomplete-intraclass----------------------------------------------------
fit_inc <- icc(ratings_incomplete, subject = subject, rater = rater, score = score,
               model = "twoway", type = "agreement", unit = "average")
gl_inc <- glance(fit_inc)
c(estimate = with(tidy(fit_inc), estimate[term == "ICC(A,k)"]),
  subjects_used = gl_inc$n_subjects,
  ratings_used = gl_inc$n_obs,
  k_eff = gl_inc$k_eff)