---
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. Every number is computed live as the page builds. Any unfamiliar
term is defined in the [*Glossary*](glossary.html).

The comparison packages are `psych`, `irr`, and `irrICC`. `psych` is a general
package for psychological measurement, and its `psych::ICC` is the most
widely used ANOVA ICC among R users. `irr` is a classical inter-rater-reliability toolkit, and its
`irr::icc` computes the same ANOVA family. `irrICC` implements 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, every subject is rated by every rater. There 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. That is a two-way design: each rater is tracked across
the subjects they score.

`intraclass` estimates the coefficients from [variance
components](glossary.html#variance-component), each a share of the total
variation traced to one source. It fits them by [REML](glossary.html#reml) in
a mixed model. REML, restricted maximum likelihood, is a way to estimate
variances that corrects maximum likelihood's downward bias. `psych` and `irr`
instead derive the coefficients 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, 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)
```

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

Every coefficient agrees to within 0.00001, so no two tools differ by more
than one in the fifth decimal place. 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`.** The
match with `psych` 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 absolute-agreement coefficient (`icc2r`), for which raters
give the same score, agrees with `intraclass`'s `ICC(A,1)` to within 0.00001:

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

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

## 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, 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")
```

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, 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)
```

An ICC computed from two subjects is not usable, whatever its value.
`intraclass` instead fits the mixed model to **every observed rating**. It
reports an [effective number of
ratings](glossary.html#effective-number-of-ratings-k_eff) (`k_eff`): the
harmonic mean of the per-subject rating counts, an average that leans toward
the smaller values. That count 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)
```

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

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. 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 guide to intent, not a scorecard: each tool is excellent at what it was
designed for. Two terms in it need a gloss. Raters are
[fixed](glossary.html#fixed-vs--random-raters) when the observed raters are
the whole population of interest. An interval is
[boundary-aware](glossary.html#zero-variance-boundary) when an estimate can
land exactly at zero and the interval still behaves.

| 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 is it the error-variance framing that distinguishes
agreement from consistency, where raters agree apart from a constant offset
per rater. `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): subjects nested in
  clusters, where the subject level is reliability within a cluster, and the
  cluster level is reliability of cluster means.
- [*Interval methods*](interval-methods.html): the Monte-Carlo interval, built
  by drawing parameter values from the fitted model's uncertainty. The article
  also covers the bootstrap, which refits the model on simulated data many
  times.
- [*Estimation engines*](engines.html): the engine, the software that does the
  fitting, in its mixed-model, SEM, and Bayesian forms.

## 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 and familiar. As the
table above shows, they agree with `intraclass` there to within 0.00001. 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**. And reach for it when you want
the package to help you **choose and justify** the coefficient in the first
place.