Confidence-interval methods

When the default under-covers

The boundary is not the default’s only weak spot, and the other one is easier to miss because nothing about the output looks wrong. The Monte-Carlo draw assumes the fitted parameters are approximately normally distributed around the truth. That assumption is about the estimates, not the ratings. But it degrades when the subject effects are strongly skewed or heavy-tailed, and then a nominal 95% interval covers considerably less often than 95% of the time.

A one-way simulation study measured this across four subject-effect distributions. Where the default produced an interval at all, coverage fell to 0.6725 at its worst: chi-square(1) subject effects, a true ICC of 0.6, 50 subjects and 5 raters. Those runs did not abort, warn, or widen, so the shortfall is invisible in the interval itself.

Two patterns in that study are worth carrying away, and one tempting reading of it is wrong. At 5 raters per subject, coverage falls as the subject count rises, once the true ICC is moderate or high. The largest designs are the exposed ones, not the smallest. And near-normal or uniform subject effects under-covered only in cells where many runs aborted. Wherever the default almost always returned an interval, those distributions were fine. The wrong reading is that fewer raters is safer: in every cell where both were measured, 2 raters covered worse than 5. What changes is that a larger share of the 2-rater runs abort outright, and an abort is a visible failure rather than a quiet one.

The held-out battery agrees where the geometry matches: lognormal and Laplace subject effects covered 0.825 and 0.84 at that same 50-subject, 5-rater geometry, while their 20-subject, 3-rater cells came out near nominal.

The natural reaction, switching to a closed form, does not help. In every cell where the default under-covered without also aborting often, "searle" and "burch" under-covered as well, usually by more: "burch", the one this article used to recommend for heavy tails, bottoms out at 0.6655. The remaining methods were never run on that study, so this article recommends none of them in its place.

What to do instead is ordinary statistical hygiene rather than a package setting. Look at the distribution of the subject means before trusting a narrow interval. Report the variance components alongside the coefficient. And treat an interval on visibly skewed data as optimistic about its own precision.

The transformed bootstrap-t (ci_method = "npbootstrap")

The non-parametric transformed bootstrap-t of Ukoumunne et al. (2003) serves the one-way random design (model = "oneway"), on balanced and unbalanced data alike. Of the four opt-in methods it is the only one that serves unbalanced one-way data. It resamples whole subjects with replacement rather than simulating from the fitted model, so it is the only opt-in method that takes a seed (and boot_samples). Pin both for a reproducible interval. Any conf_level in (0, 1) is accepted. unit = "average", the ICC(k), is the exact monotone Spearman-Brown image of the ICC(1) endpoints, so its coverage is inherited by construction, balanced or not. A numeric unit (a D-study projection) is restricted to balanced data. Reach for it for boundary robustness, an interval that exists where the Monte-Carlo default aborts, and for robustness to non-normal subject effects. Its endpoints are deliberately left untruncated on the estimator’s own support (Ukoumunne et al. 2003, §5.2), so a near-boundary lower limit can be negative. That is honest disclosure, not an error.

The classical closed forms (ci_method = "searle" and "burch")

Two deterministic classical intervals for the balanced one-way random design: closed forms with no resampling, so mc_samples, boot_samples, and seed do not apply and no std.error is reported. Any conf_level in (0, 1) is accepted. Both project ICC(k) through the same Spearman-Brown image as "npbootstrap", and a numeric unit with it. The exact-F interval ("searle"; Searle 1971, the McGraw & Wong 1996 Table 7 limits) is exact under normality and best-calibrated when the data are approximately normal. The Burch interval ("burch"; Burch 2011) is REML-based and kurtosis-adjusted: its width tracks the data’s tail weight, which buys it some robustness to mild non-normality. It is not, however, a remedy for heavy tails: on strongly skewed subject effects it under-covers about as badly as the default (see When the default under-covers).

Which is the tighter interval? Neither, reliably, and the margin between them is not a fixed one. Over the larger grid’s 64 cells, spanning four distribution families, "burch" is the narrower of the two in 59 of 64 cells of the larger grid. How much narrower depends on where in the design you look. Below, the median width ratio is "burch"’s width over "searle"’s, so a value under 1 means "burch" is narrower and 1 means they are the same width.

Two things to read off. The first is the pattern along the true ICC. "burch"’s width advantage holds much the same up to a true ICC of 0.3 rather than shrinking as the true ICC rises (on the larger grid; the smaller grid’s margin does shrink across its levels). On the larger grid, measured by level medians the largest margin is at a true ICC of 0.1, not at the bottom of the range. Cell by cell, though, "burch"’s width margin is the larger one at the bottom level in 11 of 16 paired cells. The advantage then collapses to near parity at a true ICC of 0.6, on the one grid reaching that value, where every cell favouring "searle" sits. The second is the pattern along the subject count. There "burch"’s width margin shrinks steadily as the subject count grows, measured at 5 raters, which is the cut the second table takes. That cut is at 5 raters because 10 subjects is the only subject count in either grid where the rater count varies, so an unstratified row there would be confounded with it. The smaller grid’s 16 cells carry only the two lowest true-ICC values, so that sweep cannot show the first pattern. It shows the second in the same direction, "burch" running narrower by a median width ratio of 0.9017 at 10 subjects, 0.9611 at 30 subjects and 0.9775 at 50 subjects.

A pooled figure over both grids that vary only the subject effect would hide all of this, and it would invite a second misreading. The smaller grid’s design points are a subset of the larger one’s. So much of the gap between their pooled medians is which true-ICC values each grid covers, rather than a disagreement between two bodies of evidence. But only much of it. Restricting the larger grid to the smaller one’s design points closes most of that gap and leaves a remainder. The two are separate simulations that mostly disagree at the design points they share, agreeing closely at only a couple of them. A pooled between-grid comparison is not a clean contrast in either direction, which is why nothing above is one.

Both of the grids above draw the subject effects alone from the non-normal family and always draw the errors from a normal, and that is not an incidental detail. Burch’s own expected-length comparison is against this very exact-F interval, and it is kurtosis-conditional. He finds his interval shorter for light-tailed data but wider for symmetric heavy-tailed data, measuring with the subject effects and the errors alike drawn from the studied family. A third grid now measures that residual case here. What "burch" does against "searle" depends on what the residual is drawn from, and the three grids now measure that: the two grids that vary only the subject effect put it narrower nearly everywhere, while the third, which draws the residual from the same family as the subject effect, puts it wider at every symmetric heavy-tailed family measured (a median width ratio of 1.2963 at t(5) with 100 subjects) and narrower at every lighter-tailed one, the normal included. So the honest summary is that the ordering depends on the data, on what the residual is drawn from most of all. You should not pick between them on width.

Prefer "searle": across every distribution family in that skew study it landed closer to nominal coverage in most cells, heavy-tailed ones included. What "burch" buys is dipping below the nominal level in fewer cells overall, a more limited kind of steadiness than its kurtosis adjustment suggests. Their value over the default is a finite, well-calibrated interval at the near-zero-ICC boundary where the Monte-Carlo default aborts. One asymmetry between the siblings: on data with no between-subject variance at all, "burch" aborts (its kurtosis standardization divides by zero there) while "searle" still returns an interval. Read that interval carefully: the single-rater coefficient gets the attained minimum, and the averaged projection carries that minimum through the Spearman-Brown pole to negative infinity, which a default call prints beside it.

mc <- tidy(icc(ratings, score, subject, rater, model = "oneway", seed = 1))
se <- tidy(icc(ratings, score, subject, rater,
  model = "oneway", ci_method = "searle"
))
bu <- tidy(icc(ratings, score, subject, rater,
  model = "oneway", ci_method = "burch"
))
np <- tidy(icc(ratings, score, subject, rater,
  model = "oneway", ci_method = "npbootstrap", boot_samples = 199, seed = 1
))
data.frame(
  term = mc$term,
  estimate = round(mc$estimate, 3),
  montecarlo = sprintf("[%.2f, %.2f]", mc$conf.low, mc$conf.high),
  searle = sprintf("[%.2f, %.2f]", se$conf.low, se$conf.high),
  burch = sprintf("[%.2f, %.2f]", bu$conf.low, bu$conf.high),
  npbootstrap = sprintf("[%.2f, %.2f]", np$conf.low, np$conf.high)
)
#>     term estimate   montecarlo        searle         burch   npbootstrap
#> 1 ICC(1)    0.166 [0.01, 0.83] [-0.13, 0.72] [-0.13, 0.56] [-0.05, 0.89]
#> 2 ICC(k)    0.443 [0.03, 0.95] [-0.88, 0.91] [-0.90, 0.84] [-0.26, 0.97]

All four columns share the same point estimate: ci_method selects the interval, never the estimator. The visible difference is at the lower end. The three opt-in methods’ lower limits dip below zero, because their endpoints are left untruncated on the estimator’s own support, while the Monte-Carlo interval stays inside the range. And the Burch adjustment is empirical, not a one-way widening. Here it comes out narrower than the exact-F interval, because its width tracks the tail weight these particular data actually show, which is the same direction the grids above measure.

The modified profile likelihood (ci_method = "mpl")

The modified profile-likelihood interval of Xiao & Liu (2013) is the two-way counterpart. It serves the balanced, complete two-way random absolute-agreement ICC(A,1), with ICC(A,k) and any numeric-unit projection its pole-safe Spearman-Brown image. It aborts on any other design, on consistency or fixed raters, and on unbalanced or incomplete data. It is a deterministic closed form: no resampling, no seed. Its calibration fixes two fences: conf_level must be 0.90, 0.95, or 0.99 (each level carries its own calibrated correction constant, never interpolated between levels), and the calibration grid spans 2–10 raters and 10–100 subjects. Like "npbootstrap", it returns an interval at the near-zero-ICC boundary where the two-way Monte-Carlo default aborts. It is deliberately conservative: it over-covers, is wider than the Monte-Carlo interval at interior cells, and so it is an opt-in and not the default. Two reporting caveats come from ?icc. The two-sided interval is not equal-tailed, so a limit must not be read as a one-sided bound at half the complementary level, and separately, at conf_level = 0.99 with two raters the interval can be near-vacuous.

The shipped ratings data are too small for the calibration grid (six subjects), so the demonstration simulates a balanced two-way design inside it:

set.seed(88)
n_s <- 20
n_r <- 4
subj_eff <- rnorm(n_s, sd = sqrt(0.6))
rater_eff <- rnorm(n_r, sd = sqrt(0.1))
noise <- matrix(rnorm(n_s * n_r, sd = sqrt(0.2)), n_s, n_r)
sim <- data.frame(
  subject = factor(rep(seq_len(n_s), times = n_r)),
  rater = factor(rep(seq_len(n_r), each = n_s)),
  score = as.numeric(outer(subj_eff, rep(1, n_r)) +
    outer(rep(1, n_s), rater_eff) + noise)
)
mc2 <- tidy(icc(sim, score, subject, rater, type = "agreement", seed = 1))
ml <- tidy(icc(sim, score, subject, rater, type = "agreement", ci_method = "mpl"))
data.frame(
  term = mc2$term,
  estimate = round(mc2$estimate, 3),
  montecarlo = sprintf("[%.2f, %.2f]", mc2$conf.low, mc2$conf.high),
  mpl = sprintf("[%.2f, %.2f]", ml$conf.low, ml$conf.high)
)
#>       term estimate   montecarlo          mpl
#> 1 ICC(A,1)    0.709 [0.47, 0.84] [0.42, 0.87]
#> 2 ICC(A,k)    0.907 [0.78, 0.95] [0.75, 0.96]

The two point estimates agree, from the same REML fit. The "mpl" interval is the wider of the pair at this comfortably interior cell, the conservatism described above, visible on ordinary data.

Highest-posterior-density intervals

For comparison you can ask for a highest-posterior-density interval (HPDI), the narrowest interval containing 95% of the posterior mass. Request it with posterior_summary = "hpdi":

icc(ratings, score, subject, rater, engine = "brms",
  type = "agreement", posterior_summary = "hpdi", seed = 1)

#> ── Intraclass correlation: two-way random, absolute agreement ──────────────────
#> Subjects: 6 | Raters: 4 (random) | Observations: 24 of 24 cells (complete)
#> Engine: brms (MCMC) | CI: 95% posterior credible (HPDI) (4000 draws)
#>
#>   index     estimate   95% CI
#>   ICC(A,1)     0.241   [0.040, 0.601]
#>   ICC(A,k)     0.679   [0.256, 0.904]
#>
#> Variance components: subject 1.522, rater 2.653, residual 0.962
#> Shrout & Fleiss equivalent: ICC(A,1) = ICC(2,1), ICC(A,k) = ICC(2,k)

The header now flags (HPDI), and on the same draws the interval is no wider than the percentile one. That is what “narrowest” means. Here ICC(A,1) is [0.04, 0.60] against the percentile [0.07, 0.65], and the point estimate (the MAP) is unchanged. Percentile stays the default, because HPDI is not transform-invariant and can behave less well at the variance boundary. It is offered for comparison, not as an upgrade.