---
title: "Comparison with other packages"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Comparison with other packages}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

```{r setup}
library(intraclass)
```

If you already use another R package for intraclass correlations, two questions
matter before switching: **does `intraclass` agree with the tool I trust on the
problems that tool handles?** and **what does it do that my current tool cannot?**
This article answers both, on the package's own shipped datasets, with every
number computed live as the page builds. Any unfamiliar term is defined in the
[*Glossary*](glossary.html).

The comparison packages are `psych` (Revelle's `psych::ICC`, the most widely used
ANOVA ICC in R), `irr` (`irr::icc`, a classical inter-rater-reliability toolkit),
and `irrICC` (Gwet's model-based ICCs). All three are optional: the code chunks
below only run when the package is installed.

```{r 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])
}
```

## Does it agree? (validation)

On a **balanced** design, one where every subject is rated by every rater, the
whole ICC family is defined for all of these tools. So we can line them up
coefficient by coefficient. The `ratings` dataset is six subjects each scored by
the same four raters.

`intraclass` estimates the coefficients from [variance
components](glossary.html#variance-component) fitted by [REML](glossary.html#reml)
(a mixed model), whereas `psych` and `irr` derive them from classical **ANOVA
mean squares**. Those are different computational routes to the same population
quantity, and they are known to converge to each other. The table shows how close
they land here:

```{r 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)
```

```{r validation-gap, eval = have_psych && have_irr}
max_gap <- max(abs(comparison$intraclass - comparison$psych),
               abs(comparison$intraclass - comparison$irr))
```

Every coefficient matches to five decimal places. The largest disagreement
anywhere in the table is `r if (exists("max_gap")) format(max_gap, digits = 2) else "<1e-5"`.
That residual is not error in either tool. It is the small-sample gap between a
REML fit and ANOVA mean squares, which vanishes as the sample grows. **On the
designs classical tools handle, you lose nothing by using `intraclass`**, and the
`psych` agreement is in fact checked on every test run of this package.

A model-based tool from a different lineage agrees too. `irrICC` implements
Gwet's ICCs, estimated by a moment method rather than either REML or ANOVA. Its
two-way random agreement coefficient (`icc2r`) reproduces `intraclass`'s
`ICC(A,1)`:

```{r 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)
)
```

## What does it add? (differentiation)

The classical tools were built for the balanced, complete case. Real rating data
are rarely so tidy, and that is where the packages diverge.

### Incomplete and unbalanced data

The `ratings_incomplete` dataset is the same study with four ratings missing. In
particular, the second rater scored only two of the six subjects:

```{r 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)
```

A classical ANOVA ICC needs a complete rectangle, so `psych` and `irr`
**listwise-delete** any subject with a missing cell. Here that discards the four
subjects rater 2 skipped, leaving only two:

```{r 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)
```

An ICC computed from two subjects is not usable, whatever its value. `intraclass`
instead fits the mixed model to **every observed rating** and reports an
[effective number of ratings](glossary.html#effective-number-of-ratings-k_eff)
(`k_eff`) that accounts for the imbalance:

```{r 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)
```

All six subjects and all twenty observed ratings contribute, and nothing is thrown
away. `irrICC` can also fit incomplete data with its own model, as the capability
matrix below shows, but the mean-squares tools cannot.

### The bigger picture

Agreement on balanced data and graceful handling of missing data are two entries
in a wider gap. The table below summarizes what each package computes. It is a map
of intent, not a scorecard: each tool is excellent at what it was designed for.

| Capability | `psych` | `irr` | `irrICC` | `intraclass` |
|---|:--:|:--:|:--:|:--:|
| Balanced ANOVA ICC family | ✅ | ✅ | ✅ | ✅ |
| Incomplete / unbalanced data | no | no | ✅ | ✅ |
| Multilevel (subject **and** cluster) IRR | no | no | no | ✅ |
| Boundary-aware interval | no | no | partial | ✅ |
| [Fixed vs. random](glossary.html#fixed-vs--random-raters) rater framing | partial | partial | no | ✅ |
| Guidance on *which* ICC to report | no | no | no | ✅ |

Two rows deserve a word. Model-based extractors such as `performance::icc`
return **variance components** or a variance-partition coefficient. That is the
raw material of an ICC, but not the inter-rater-reliability coefficient family
itself, nor the error-variance framing that distinguishes agreement from
consistency. `intraclass`'s own generalizability coefficients were validated
against `gtheory`, agreeing to within 0.001. `gtheory` is a
generalizability-theory package archived from CRAN in March 2025, and is not a
dependency here. Those committed reference values live in the package's reference
notes. And an **interval** that is
[boundary-aware](glossary.html#monte-carlo-interval) is something none of the
classical tools provide. Such an interval behaves correctly when a variance
component is estimated at its [zero
boundary](glossary.html#zero-variance-boundary), where a normal-approximation
interval silently misbehaves.

`intraclass` earns its extra machinery on exactly these cases. For the details of
each, see the companion articles:

- [*Choosing an ICC*](choosing-an-icc.html): the selection framework the last
  matrix row points to.
- [*Multilevel designs*](multilevel-designs.html): subject- and cluster-level
  reliability when raters are nested.
- [*Interval methods*](interval-methods.html): the boundary-aware Monte-Carlo
  and bootstrap intervals.
- [*Estimation engines*](engines.html): the mixed-model, SEM, and Bayesian
  backends behind these numbers.

## When to use which

If your design is **balanced and complete** and you only need the classic
McGraw–Wong coefficients, `psych` and `irr` are mature, familiar, and, as the
table above shows, numerically identical to `intraclass`. Reach for `intraclass`
when your data are **incomplete or unbalanced**, or when raters are **nested in
clusters**. Reach for it too when you need an **interval you can trust near the
boundary**, or when you want the package to help you **choose and justify** the
coefficient in the first place.