---
title: "D-studies and within-cell replicates"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{D-studies and within-cell replicates}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

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

A fitted `icc()` carries its [variance
components](glossary.html#variance-component), each a share of the total
variation traced to one source. Two things you can do with that fit. The first
is a [**decision study**](glossary.html#d-study-decision-study), which
projects the fitted variance components to other rater or
[occasion](glossary.html#occasion-within-cell-replicate) counts.
The second is
to separate the
subject-by-rater interaction from pure error. That needs **within-cell
replicates**, each one of several ratings by the same rater of the same
subject. This article also shows the `autoplot()` methods that visualize a
fit. Unfamiliar terms are defined in the [*Glossary*](glossary.html).

## How many raters do I need? A D-study

`ICC(*,1)` is the reliability of a *single* rater: the reliability of one
rater's score. `ICC(*,k)` is the reliability of the mean of the `k` raters
you actually used. A **decision (D-) study** asks a forward-looking question:
*how reliable would the mean of some other number of raters `m` be?* In
generalizability theory the absolute-agreement ICC (raters give the same
score) has its own name. It is the [dependability
coefficient](glossary.html#dependability-coefficient), the
generalizability-theory name for the agreement ICC. Projecting it to `m`
raters is just a change of the averaging divisor,

$$\Phi(m) = \frac{\sigma^2_s}{\sigma^2_s + (\sigma^2_r + \sigma^2_{res}) / m},$$

so `d_study()` reuses the fit you already have, with no refitting.

```{r dstudy, eval = requireNamespace("glmmTMB", quietly = TRUE)}
fit <- icc(ratings, score, subject, rater, type = "agreement", seed = 1)
proj <- d_study(fit, m = 1:8, seed = 1)
proj
```

Reliability climbs with more raters, but with diminishing returns. The
projection is anchored to what you observed. At `m = 4`, the number of raters
in `ratings`, `Φ(m)` is exactly the `ICC(A,k)` you would get from `icc()`
directly.

The projection has its own interval settings: `conf_level`, `mc_samples` and
`seed`. Each defaults to the fit's own whenever the fit carries them. A default
[Monte-Carlo](glossary.html#monte-carlo-interval) fit, whose interval is built
by drawing parameter values from the fitted model's uncertainty, does carry them. Pass any of them to
`d_study()` to override it for this call alone. The `seed = 1` above is the
fit's own seed passed explicitly. So the `m = 4` row reproduces the fit's
`ICC(A,k)` *interval* exactly, and not merely its estimate.

`tidy()` turns the projection into one row per projected point, and `glance()`
summarizes the projection as a whole:

```{r dstudy-tidy, eval = requireNamespace("glmmTMB", quietly = TRUE)}
tidy(proj)

glance(proj)
```

The one column here that `tidy()` on the *fit* does not also give you is
**`m`**, the rater-count column. That column is the whole point of a D-study:
the same coefficient, indexed by the rater count it is projected to.
`glance()` describes the projection rather than the model. It reports the
swept range (`n_m`, `m_min`, `m_max`) alongside the design and the interval
settings.

We fit with `type = "agreement"` here because the dependability coefficient
\eqn{\Phi(m)} above *is* the absolute-agreement projection. A default `icc()`
reports both agreement and [consistency](glossary.html#consistency), where
raters agree apart from a constant offset per rater. `d_study()` then projects
**one curve per error definition**, told apart by the `type` column. That is
handy when you want both, but for this walkthrough one curve keeps the picture
clear.

Read as a curve, this is the classic "how many raters?" picture. Plot it with
`autoplot()`, which needs **ggplot2**:

```{r dstudy-plot, eval = requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("glmmTMB", quietly = TRUE), fig.alt = "Projected reliability rising with the number of raters, with a Monte-Carlo interval band."}
library(ggplot2)
autoplot(d_study(fit, m = 1:12))
```

There is also a `plot()` method, for the console habit of drawing a picture as
a side effect. `plot()` on a projection prints the very same ggplot that
`autoplot()` builds from it, then returns the projection invisibly. We show
the call rather than run it. Running it would render another reliability
curve just like the one above:

```{r dstudy-plot-wrapper, eval = FALSE}
plot(proj)
```

### One projected value, without a projection object

When you want a single projected rater count rather than a whole curve,
`icc()` will report it inline. `unit` takes the keywords `"single"` and
`"average"`, which give the `ICC(*,1)` and `ICC(*,k)` you always get.
Alongside them it also takes any number `m >= 1`:

```{r dstudy-unit, eval = requireNamespace("glmmTMB", quietly = TRUE)}
icc(ratings, score, subject, rater,
  type = "agreement", unit = c("single", "average", 6), seed = 1
)
```

The extra row is labeled `ICC(A,6)`, and it is the same quantity `d_study()`
projects: the `m = 6` row of `proj` above carries the same estimate. The
*interval* matches too, but only because both calls pass the same `seed`. The
interval is Monte-Carlo, so two unseeded runs of the same projection agree on
the estimate and differ in the last digits of the endpoints. Which one you
reach for is a matter of shape, not of arithmetic. `unit` adds a row to a
coefficient table. `d_study()` returns a projection object you can `tidy()`,
`glance()` and plot.

**Fixed raters, where the observed raters are the whole population of
interest, cannot be projected to absolute agreement.** The rater term for
[*fixed* raters](glossary.html#fixed-vs--random-raters) is the variance of
exactly the raters you observed. There is no wider pool for a hypothetical
sixth rater to be drawn from, so `icc()` refuses rather than returning a
number:

```{r dstudy-unit-fixed, error = TRUE, eval = requireNamespace("glmmTMB", quietly = TRUE)}
icc(ratings, score, subject, rater,
  type = "agreement", raters = "fixed", unit = c("single", "average", 6),
  seed = 1
)
```

Consistency has no such problem, because the rater term drops out of it. So
the two remedies the message names are to project with `raters = "random"`,
or to ask for a fixed-rater `type = "consistency"` projection instead. The
refusal is raised only when absolute agreement is *all* you asked for. A
default `icc()` call reports both error definitions. There the consistency
projection is kept and you are told the agreement one was dropped.

**Projection is extrapolation.** The rater variance $\sigma^2_r$ is estimated
from only as many raters as you observed, so projecting far beyond that design
leans hard on that estimate. The Monte-Carlo interval widens honestly to
reflect this, rather than pretending to a precision it lacks. And projecting
absolute agreement is refused for [*fixed*
raters](glossary.html#fixed-vs--random-raters), where there is no wider rater
universe to generalize to (use `raters = "random"`). The [D-study also works
on a multilevel fit](multilevel-designs.html#how-many-raters-a-multilevel-d-study),
projecting the rater count at each level.

## Within-cell replicates: interaction vs. pure error

So far every subject-by-rater cell holds a single rating. When each rater
rates each subject **more than once**, the design has *within-cell
replicates*. Then you can separate two things that a single rating confounds.
The first is the **subject-by-rater interaction**: does a rater systematically
score a particular subject high or low, a *stable* disagreement? The second is
**pure error**, meaning how much a rater's repeat ratings of the same subject
wobble. The shipped `ratings_replicates` data are simulated for this case:
20 subjects, 4 raters, and 3 ratings in every subject-by-rater cell (see
`?ratings_replicates`). `icc()` detects the replicates and fits the
interaction model automatically:

```{r replicates-data, eval = requireNamespace("glmmTMB", quietly = TRUE)}
icc(ratings_replicates, score, subject, rater, type = "agreement", occasions = c("single", "average"))
```

The variance-components line now shows `subject:rater` (the interaction) and
`residual` (pure error) as separate terms. The **single-occasion** rows, whose
`occasions` column reads 1, are the ordinary ICCs, since a single rating's
error still includes the interaction. But they are now fit correctly, rather
than folding the interaction into the residual. The **occasion-averaged**
rows, whose `occasions` column reads 3 here, divide pure error by three
occasions per rater. `ICC(A,1)` is then the reliability of one rater's *mean
of three ratings*. `ICC(A,k)` is the reliability of a mean over all four
raters, each at that same three. That 3 is the fitted per-cell replicate
count. Averaging cuts pure error but not the interaction, so those
coefficients are higher.

Within-cell replicates extend beyond this balanced two-way random example,
where each rater is tracked across the subjects they score. **Fixed** raters (balanced) are
supported. So are **multilevel** designs (crossed Design 1 and nested Design
2, balanced), and **ragged** replicates at a single occasion. What remains
open is the occasion-*averaged* coefficient on ragged replicates. There
unequal per-cell counts leave no single effective-occasion divisor with a
validated oracle. The compound fixed-by-ragged and multilevel-by-ragged
corners are open too.

## How many occasions do I need? A D-study on the occasion facet

Just as `d_study(m = ...)` projects the number of *raters*, `d_study(n_o = ...)`
projects the number of *occasions* off a replicate fit. It holds the raters
fixed and asks "how reliable would each rater's mean of `n_o` ratings be?".
Supply exactly one axis per call (`m` **or** `n_o`).

```{r occasion-dstudy, eval = requireNamespace("glmmTMB", quietly = TRUE)}
fit_rep <- icc(ratings_replicates, score, subject, rater, type = "agreement", occasions = "average")
d_study(fit_rep, n_o = 1:6)
```

Notice the curve **flattens**. Averaging more occasions only cancels *pure
error* (`residual`). It never touches the rater or `subject:rater` variance.
So the occasion curve climbs to a **ceiling below 1**, rather than approaching
1 the way a rater projection does. That ceiling is the reliability you would
reach with perfectly repeatable ratings but the same raters. Read it as "how
much does re-rating help?", which saturates.

Because occasions are a **random** facet however the raters are treated, the
occasion projection is defined even where a rater projection is not. In
particular, **fixed-rater absolute agreement projects on the occasion axis**.
It is only the *rater* axis that is undefined for fixed absolute agreement,
having no "freshly sampled rater" to add. On a **multilevel** replicate fit
the subject level is reliability within a cluster, and the cluster level is
reliability of cluster means. There the subject-level curve rises with `n_o`,
while the [cluster-level](glossary.html#subject-level-vs--cluster-level) curve
is **flat**. The cluster-level
error set has no pure-error term, so occasions cannot change it, and
`d_study()` says so with a note.

## Visualizing a fit

Every `icc()` fit carries an `autoplot()` method, with a `plot()` wrapper. So
you can *see* the coefficients and the variance components behind them without
building a plot by hand. Both read straight off the fitted object, so the
picture can never disagree with the printed table. They need **ggplot2**, an
optional dependency.

The default, `what = "coefficients"`, is a **forest plot**: each ICC index as
a point estimate with its Monte-Carlo interval. Reusing the two-way `ratings`
fit from the D-study section above:

```{r plot-coef, eval = requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("glmmTMB", quietly = TRUE), fig.alt = "Forest plot of ICC(A,1) and ICC(A,k) for the ratings data, each a point estimate with a horizontal Monte-Carlo interval, ICC(A,k) higher than ICC(A,1), and its interval slightly wider."}
library(ggplot2)
autoplot(fit) # `fit <- icc(ratings, score, subject, rater, type = "agreement", seed = 1)`
```

`what = "components"` shows the other half of the story: the estimated
**variance components** the ratio is built from. It makes plain *why* absolute
agreement is so much lower than the averaged coefficient on `ratings`. The
**rater** component is large, and only absolute agreement counts between-rater
differences as error:

```{r plot-comp, eval = requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("glmmTMB", quietly = TRUE), fig.alt = "Bar chart of the estimated variance components for the ratings fit: subject, rater, and residual, with the rater component the largest."}
autoplot(fit, what = "components")
```

For a **multilevel** fit the forest plot facets by level. See [*Multilevel
designs*](multilevel-designs.html#visualizing-the-levels) for that example.