Notation

• \(d\) is the difference between the two measurements (x - y), and \(b\) is its estimated mean (the bias, the estimated marginal mean from emmeans) for a given row of the output (a condition and/or a value of the average when prop_bias = TRUE).
• \(\beta\) = pred_level is the content, the proportion of differences the limits should cover, and \(\alpha_1 = 1 - \beta\).
• \(\gamma\) = tol_level is the confidence of the tolerance limits, and \(\alpha_2 = 1 - \gamma\).
• \(z = z_{1 - \alpha_1/2}\), the standard normal quantile (1.96 for \(\beta = 0.95\)).
• \(\chi^2_{p, \nu}\) is the \(p\) quantile of a chi-square distribution with \(\nu\) degrees of freedom, and \(t_{p, \nu}\) the \(p\) quantile of a \(t\) distribution.
• With repeated measures, every interval refers to a single new difference from a new subject: its variance includes the between-subject and within-subject variance.

Arguments Influencing the Model

The only required arguments are x, y, and data, which give the data frame and the columns that contain the two measurements.

• id identifies the subjects (or other clusters) within which differences are correlated. With model = "lme", it can name two nested levels, outer level first (e.g., id = c("golfer", "club")).
• time gives the order of the measurements within a subject, for the autoregressive correlation structures.
• condition identifies a factor that may change the mean and the variance of the differences; it adds a separate mean and a separate residual variance for each condition.
• prop_bias = TRUE adds the average of the two measurements as a covariate, and the limits are reported at its minimum, median, and maximum.
• model and cor_type set the model for the correlation between differences from the same subject (below).
• weights and correlation let you specify an nlme variance function or correlation structure directly.

The model

With model = "gls" (default), the differences follow a marginal generalized least squares model (nlme::gls()):

\[ d_{ij} = \mathbf{x}_{ij}^\top \beta + e_{ij}, \qquad \text{Var}(e_{ij}) = \sigma^2 g_{ij}^2, \qquad \text{Cor}(e_{ij}, e_{ij'}) = R_{i, jj'} \]

where \(g_{ij}\) is the variance function (1 without one) and \(R_i\) is the correlation matrix for subject \(i\): compound symmetry (cor_type = "sym", a common correlation \(\rho\)), AR(1) or continuous AR(1) ("ar1", "car1"), or none. The variance of a single difference for a given row of the output is \(S^2 = \hat\sigma^2 g^2\).

With model = "lme", the differences follow a linear mixed model (nlme::lme()) with a random intercept for each subject and, optionally, for each setting (inner level) within subject:

\[ d_{ijk} = \mathbf{x}_{ijk}^\top \\beta + u_i + v_{ij} + e_{ijk}, \qquad u_i \sim N(0, \sigma_1^2), \quad v_{ij} \sim N(0, \sigma_2^2), \quad \text{Var}(e_{ijk}) = \sigma^2 g_{ijk}^2 \]

where \(v_{ij}\) is present only with two id columns, the variance function applies to the residuals only, and cor_type = "ar1" or "car1" adds serial correlation of the residuals within the innermost level. The variance of a single difference for a given row is

\[ S^2 = \hat\sigma_1^2 + \hat\sigma_2^2 + \hat\sigma^2 g^2 \]

The output reports the components as SD.between (\(\hat\sigma_1\)), SD.nested (\(\hat\sigma_2\)), and SD.within (\(\hat\sigma g\)). For gls with compound symmetry, SD.between \(= \sqrt{\hat\rho}\, S\) and SD.within \(= \sqrt{1 - \hat\rho}\, S\).

The bias and its standard error

The bias \(b\) for each row and its standard error (SEM) are the estimated marginal mean and its standard error from emmeans. The confidence interval for the bias is \(b \pm t_{1 - (1 - c)/2,\, df_b} \cdot SEM\), where \(c\) is conf_level. The degrees of freedom \(df_b\) are:

• model = "gls": the approximation of Satterthwaite (1946) (see Kuznetsova et al. (2017) for its implementation in R).
• model = "lme": the “containment” degrees of freedom, which are the number of subjects minus 1 for the models fit here. The Satterthwaite approximation for lme models in emmeans can fail (for a model with only an intercept) or not terminate (with a variance function), so it is not used.

Prediction interval

The standard error of prediction (SEP) combines the uncertainty in the bias with the variance of a single difference:

\[ SEP = \sqrt{SEM^2 + S^2} \]

and the prediction interval is

\[ PI = b \pm t_{1-\alpha_1/2,\, df_b} \cdot SEP \]

This is a “beta-expectation” tolerance interval: on average over repeated studies, it contains a proportion \(\beta\) of the differences. It is not a confidence statement.

Upper confidence bound for the SD

Both kinds of tolerance limits need a one-sided upper confidence bound \(S_U\), at level \(\gamma\), for the SD of a single difference. The difficulty with repeated measures is that \(S^2\) is a sum of variance components that are estimated with very different amounts of information: the between-subject variance with only (subjects \(- 1\)) degrees of freedom, and the within-subject variance with many more. A single chi-square distribution with a single “effective” degrees of freedom misses the skewness of the between-subject part and gives limits that are too narrow (coverage of about 91% for a 95% target in our simulations). Instead, \(S^2\) is split into pieces, each with its own degrees of freedom, and the pieces are combined with the MOVER method (Zou 2011; Donner and Zou 2012), as for the nested limits of agreement further below.

Independent data (gls without id, or with cor_type = "none"): with \(\nu\) the residual degrees of freedom (\(N - p\) without a variance function),

\[ S_U = S \sqrt{\frac{\nu}{\chi^2_{\alpha_2, \nu}}} \]

Variance components (gls with compound symmetry, or lme): write \(S^2 = \sum_j a_j\), where piece \(j\) is estimated with \(\nu_j\) degrees of freedom. Then

\[ S_U^2 = S^2 + \sqrt{\sum_j \left[ a_j \left( \frac{\nu_j}{\chi^2_{\alpha_2, \nu_j}} - 1 \right) \right]^2 } \]

Pieces with zero degrees of freedom are dropped. The pieces are the balanced nested analysis of variance decomposition of \(S^2\) into the variance of a subject mean and the remaining within-subject variation, generalized to unbalanced data and to residual correlation:

where:

• \(n\) is the number of subjects, \(n_c\) the number of settings (the inner level; \(n_c = n\) without nesting, so the middle piece is empty), and \(N\) the number of measurements.
• \(\hat\sigma_e^2 = \hat\sigma^2 g^2\) is the residual variance for the row, and \(\hat\rho\) is the compound symmetry correlation (truncated at 0).
• \(m_h = n / \sum_i m_i^{-1}\) is the harmonic mean number of measurements per subject (\(m_i\) for subject \(i\)).
• \(\kappa_s = \frac{1}{n} \sum_i \mathbf{1}^\top R_i \mathbf{1} / m_i^2\) is the residual variance of a subject mean relative to that of one measurement, averaged over subjects, and \(\kappa_c\) is the same for setting means. Without residual correlation, \(\kappa_s\) is the mean of \(1/m_i\) (so \(\kappa_s = 1/m_h\)) and \(\kappa_c\) the mean of \(1/m_c\); with AR(1) residuals, \(\mathbf{1}^\top R_i \mathbf{1}\) adds the correlations between measurements.
• \(\tau = \frac{1}{n} \sum_i \sum_c m_{ic}^2 / m_i^2\) is the share of the setting variance in the variance of a subject mean (\(1/s\) with \(s\) settings of equal size).

For balanced data without residual correlation, these are exactly the pieces of the classical decomposition \(S^2 = MS_A / (sk) + MS_B (1/k - 1/(sk)) + MS_W (1 - 1/k)\) for \(n\) subjects, \(s\) settings per subject, and \(k\) measurements per setting.

When condition is supplied, each condition has its own residual variance, and each condition’s variance components are estimated from the measurements in that condition. The counts (\(n\), \(n_c\), \(N\), \(m_i\), \(m_h\)) and \(\kappa\), \(\tau\) are then computed from the rows of that condition only. (Using the counts of the whole data set gave limits that were too narrow: about 93% coverage for a 95% target.)

Other correlation structures (gls with AR(1), continuous AR(1), or a user-specified structure): \(S_U = S \sqrt{\nu / \chi^2_{\alpha_2, \nu}}\) with an effective degrees of freedom

\[ \nu = \frac{2 \cdot S^4}{\widehat{\text{Var}}(S^2)} \]

where \(\widehat{\text{Var}}(S^2)\) comes from the delta method applied to the approximate covariance matrix of the variance parameters (\(\log \hat\sigma\) and the variance function parameters) that nlme computes (apVar). If that matrix is not available, \(df_b\) is used instead, with a warning.

The output reports \(S_U\) as SD.upper and the degrees of freedom of \(S^2\) as SD.df (for the variance components, the Satterthwaite value \(S^4 / \sum_j a_j^2 / \nu_j\), reported for information only).

Tolerance limits: bound_type = "joint"

See Joint versus Intersection-Union (IU) Bounds for when to use the joint or the IU limits.

The joint tolerance limits are a “beta-content”, “gamma-confidence” tolerance interval: with confidence \(\gamma\), at least a proportion \(\beta\) of all differences lie within the limits. They use the approximation of Howe (1969), as described by Francq et al. (2020):

\[ TI = b \pm z \cdot SEP \cdot \frac{S_U}{S} \]

For independent data with only an intercept, \(SEP = S \sqrt{1 + 1/n}\) and \(S_U / S = \sqrt{(n-1)/\chi^2_{\alpha_2, n-1}}\), and this is exactly Howe’s tolerance factor. If the limits lie within a maximal allowable difference \(\pm\Delta\), one can conclude, at level \(\alpha_2\), that at least a proportion \(\beta\) of the differences lie within \(\pm\Delta\). The tails do not need to be split equally.

Tolerance limits: bound_type = "iu"

The “iu” limits are one-sided \(\gamma\) confidence bounds on each of the two limits of agreement, \(\theta_L = \mu - z\sigma\) and \(\theta_U = \mu + z\sigma\) (the \((1 \mp \beta)/2\) percentiles of the differences). If both bounds lie within \(\pm\Delta\), one can reject, at level \(\alpha_2\), that either limit of agreement lies outside \(\pm\Delta\) (an intersection-union test). The two bounds hold together only about \(1 - 2\alpha_2\) of the time, so they are not a joint \(\gamma\) interval.

Independent data (gls without a correlation structure): the bounds are exact, from the noncentral \(t\) distribution,

\[ b \mp t'_{\gamma,\, \nu}\!\left(\frac{zS}{SEM}\right) SEM \]

where \(t'_{\gamma, \nu}(\lambda)\) is the \(\gamma\) quantile of a noncentral \(t\) distribution with \(\nu\) degrees of freedom and noncentrality \(\lambda\). (When \(\lambda > 37\), where the noncentral \(t\) quantile is numerically unreliable, the MOVER bound below is used.)

Otherwise, the MOVER bound (Zou 2011) combines the one-sided bounds for the bias and for \(S\):

\[ H = \sqrt{\left(t_{\gamma,\, df_b}\, SEM\right)^2 + z^2 \left(S_U - S\right)^2} \]

\[ L = (b - zS) - H, \qquad U = (b + zS) + H \]

Bootstrap calibration: tol_method = "boot_cal" (experimental)

The bootstrap calibrates the nominal confidence level of the analytic limits (Loh 1987; Beran 1987):

1. Simulate \(B\) (replicates) new data sets from the fitted model, using the point estimates: \(\mathbf{d}^* = \mathbf{X}\hat{\beta} + \mathbf{e}^*\), where \(\mathbf{e}^*\) has the fitted marginal covariance. For subject \(i\) this is \(D_i R_i D_i\) for gls and \(\hat\sigma_1^2 J + \hat\sigma_2^2 B_i + D_i R_i D_i\) for lme, where \(D_i\) holds the residual SDs, \(J\) is a matrix of ones, and \(B_i\) has ones for pairs of measurements in the same setting.
2. Refit the model to each data set, and compute the replicate bias \(b^*\), SD \(S^*\), SEM\(^*\), and variance components (or effective degrees of freedom). Replicates whose model fails to converge are dropped, with a warning if more than 5% fail.
3. Calibrate. For a nominal level \(\lambda\), compute the analytic limits \(L^*(\lambda)\) and \(U^*(\lambda)\) in each replicate (reusing \(df_b\) from the original fit). The fitted model plays the role of the truth, with differences distributed as \(N(b, S^2)\).
   • "joint": the coverage at level \(\lambda\) is \(C(\lambda) = \frac{1}{B} \sum \mathbb{1}\left\{ \Phi\!\left(\frac{U^* - b}{S}\right) - \Phi\!\left(\frac{L^* - b}{S}\right) \ge \beta \right\}\).
   • "iu": separately for each bound, \(C_L(\lambda) = \frac{1}{B} \sum \mathbb{1}\{L^*(\lambda) \le b - zS\}\) and \(C_U(\lambda) = \frac{1}{B} \sum \mathbb{1}\{U^*(\lambda) \ge b + zS\}\).
   The calibrated level \(\hat\lambda\) is the smallest \(\lambda \in [0.5, 1)\) with coverage of at least \(\gamma\) (found by bisection).
4. Report the analytic limits for the observed data at the calibrated level(s), which are returned as lower.TL.level and upper.TL.level.

Because the replicates are simulated from the estimated variance components as if they were known, the calibration has its own error when there are few subjects: in our simulations with 10 to 15 subjects, its coverage ranged from about 93% to 96% for a 95% target, while the analytic limits were 95% to 96.5% in the same settings. The analytic limits are therefore recommended, and the bootstrap may be useful as a check for models whose analytic limits have not been evaluated (e.g., other variance functions or proportional bias).

Summary of the degrees of freedom

Coverage in simulations

The table summarizes the coverage of the analytic limits in our simulations (95% content and 95% confidence; 1,000 simulated data sets each, so the Monte Carlo standard error is about 0.7 percentage points). The full record, with the scripts, is kept in the package’s source repository (references/general/tolerance_simulations).

Model assumptions

A few modelling choices deserve attention (see the “Model assumptions” section of ?tolerance_limit for details):

• With the default marginal model (model = "gls"), the autoregressive correlation structures (cor_type = "ar1" or "car1") have no persistent subject effect and the correlation between measurements from the same subject decays towards zero over time. If subjects have a persistent bias of their own, these structures understate the uncertainty in the bias and give limits that are too narrow. Setting model = "lme" fits a random intercept for each subject and adds the serial correlation to the residuals, which keeps the persistent subject effect (see the simulations above).
• With model = "gls", compound symmetry, and a variance function (e.g., from condition), the between-subject variance scales with the residual standard deviation of each condition. With model = "lme", the variance function applies to the residuals only, and the between-subject variance is shared across conditions.
• With model = "lme", id can name two nested levels, outer level first (e.g., id = c("golfer", "club") for shots with several clubs per golfer). Avoid setting id to the inner level alone (e.g., a golfer-by-club identifier): that drops the correlation across clubs within a golfer and gives limits that are too narrow. More than two levels are not supported.
• A slope of the differences on the average (prop_bias = TRUE) can appear without any true proportional bias when the two methods have unequal measurement error variances (Bland and Altman 1999).

For example, a random intercept model with serial correlation of the residuals, and one with nested random intercepts:

# random intercept for each subject, plus AR(1) correlation of the residuals
tolerance_limit(
  data = temps,
  x = "trec_pre",
  y = "teso_pre",
  id = "id",
  time = "trial_num",
  model = "lme",
  cor_type = "ar1"
)

# nested random intercepts: clubs within golfers
tolerance_limit(
  data = golf, # a data set with columns golfer, club, x, and y
  x = "x",
  y = "y",
  id = c("golfer", "club"),
  model = "lme"
)

Calculation Steps

The reported limits of agreement are derived from the work of Bland and Altman (1986) and Bland and Altman (1999). Throughout, \(z_A = z_{1-(1-agree)/2}\) is the normal quantile for the agreement level (1.96 for the default of 95%), and \(z_{1-\alpha}\) and \(t_{1-\alpha, df}\) use the one-sided level \(\alpha\) for each confidence bound (see the note above).

LoA

\[ LoA = \bar d \pm z_A \cdot S_d \]

wherein \(\bar d\) is the mean of the \(N\) differences and \(S_d\) is their standard deviation.

Confidence Interval

1. Calculate the standard error of the LoA

\[ S_{LoA} = S_d \cdot \sqrt{\frac{1}{N}+ \frac{z^2_A}{2 \cdot(N-1)} } \]

2. Calculate the limit of agreement margin of error (LME)

Bland-Altman Method

\[ LME = t_{1 - \alpha, \, N-1} \cdot S_{LoA} \]

MOVER Method

\[ LME = S_d \cdot \sqrt{\frac{z_{1-\alpha}^2}{N} + z^2_A \cdot \left(\sqrt{\frac{N-1}{\chi^2_{\alpha, N-1}}}-1\right)^2} \]

where \(\chi^2_{\alpha, N-1}\) is the lower \(\alpha\) quantile of the chi-square distribution.

3. Calculate Confidence Interval

\[ Lower \space LoA \space C.I. = Lower \space LoA - LME \]

\[ Upper \space LoA \space C.I. = Upper \space LoA + LME \]

Calculative Steps

With replicates, x and y need not be measured the same number of times, so their within-subject variances are estimated separately (Zou 2011).

1. Compute the subject means and the within-subject variances

For subject \(i\) (\(i = 1, \dots, n\)), with \(m_{xi}\) measurements of x and \(m_{yi}\) measurements of y:

\[ \bar d_i = \bar x_i - \bar y_i, \qquad \bar d = \frac{1}{n} \Sigma_{i=1}^{n} \bar d_i \]

and \(s_{xi}^2\) and \(s_{yi}^2\) are the sample variances of the replicates of x and of y within subject \(i\).

2. Compute pooled estimates of the within-subject variances

\[ s_{xw}^2 = \Sigma_{i=1}^{n} \frac{m_{xi}-1}{N_x-n} \cdot s_{xi}^2, \qquad s_{yw}^2 = \Sigma_{i=1}^{n} \frac{m_{yi}-1}{N_y-n} \cdot s_{yi}^2 \]

where \(N_x = \Sigma_i m_{xi}\) and \(N_y = \Sigma_i m_{yi}\).

3. Compute the variance of the subject mean differences (between-subject variance)

\[ s^2_b = \Sigma_{i=1}^n \frac{ (\bar d_i - \bar d)^2}{n-1} \]

4. Compute the harmonic means of the replicate sizes

\[ m_{xh} = \frac{n}{\Sigma_{i=1}^n m_{xi}^{-1}}, \qquad m_{yh} = \frac{n}{\Sigma_{i=1}^n m_{yi}^{-1}} \]

5. Compute the variance of a single difference

\[ s_d^2 = s^2_b + \left(1 - \frac{1}{m_{xh}}\right) s_{xw}^2 + \left(1 - \frac{1}{m_{yh}}\right) s_{yw}^2 \]

6. Calculate the LME

MOVER Method

\[ u_1 = \left[s_b^2 \left(\frac{n-1}{\chi^2_{\alpha, n-1}} - 1\right)\right]^2 \]

\[ u_2 = \left[\left(1-\frac{1}{m_{xh}}\right) s_{xw}^2 \left(\frac{N_x-n}{\chi^2_{\alpha, N_x-n}} - 1\right)\right]^2 \]

\[ u_3 = \left[\left(1-\frac{1}{m_{yh}}\right) s_{yw}^2 \left(\frac{N_y-n}{\chi^2_{\alpha, N_y-n}} - 1\right)\right]^2 \]

\[ u = s_d^2 + \sqrt{u_1 + u_2 + u_3} \]

\[ LME = \sqrt{\frac{z_{1-\alpha}^2 \cdot s_b^2}{n} + z_A^2 \cdot \left(\sqrt{u}-\sqrt{s^2_d}\right)^2} \]

Bland-Altman Method

\[ S_{LoA}^2 = \frac{s^2_b}{n} + \frac{z_A^2}{2 \cdot s_d^2} \cdot \left(\frac{s_b^4}{n-1} + \left(1-\frac{1}{m_{xh}}\right)^2 \frac{s_{xw}^4}{N_x - n} + \left(1-\frac{1}{m_{yh}}\right)^2 \frac{s_{yw}^4}{N_y - n}\right) \]

\[ LME = z_{1-\alpha} \cdot S_{LoA} \]

7. Calculate LoA

\[ LoA_{lower} = \bar d - z_A \cdot s_d \]

\[ LoA_{upper} = \bar d + z_A \cdot s_d \]

8. Calculate LoA CI

\[ Lower \space CI = LoA_{lower} - LME \]

\[ Upper \space CI = LoA_{upper} + LME \]

Calculation Steps

1. Model

A linear mixed model with a random intercept for each subject (fit with lme4::lmer, REML) is used to estimate the bias (mean difference) and the variance components:

\[ \begin{aligned} d_{ij} &\sim N \left(\mu + u_{i}, \sigma^2_w \right) \\ u_{i} &\sim N \left(0, \sigma^2_b \right) \text{, for subject } i = 1, \dots, n \end{aligned} \]

2. Extract Components

The between-subject variance \(s_b^2\) (the estimate of \(\sigma^2_b\)) and the within-subject variance \(s_w^2\) (the estimate of \(\sigma^2_w\)) are estimated from the model. Their sum is the variance of a single difference, \(s_d^2 = s_b^2 + s_w^2\), and the intercept is the bias \(\bar d\).

The harmonic mean of the number of measurements per subject (\(m_i\) for subject \(i\)) is also calculated:

\[ m_h = \frac{n}{\Sigma_{i=1}^n m_i^{-1}} \]

and \(N = \Sigma_i m_i\) is the total number of measurements.

3. Compute LME

MOVER Method

\[ u_1 = \left[s_b^2 \left(\frac{n-1}{\chi^2_{\alpha,n-1}}-1\right)\right]^2 \]

\[ u_2 = \left[\left(1-\frac{1}{m_h}\right) s_w^2 \left(\frac{N-n}{\chi^2_{\alpha,N-n}} -1\right)\right]^2 \]

\[ u = s_d^2 + \sqrt{u_1 + u_2} \]

\[ LME = \sqrt{\frac{z^2_{1-\alpha} \cdot s_b^2}{n} + z^2_A \cdot \left(\sqrt{u} - \sqrt{s^2_d}\right)^2} \]

Bland-Altman Method

\[ S_{LoA}^2 = \frac{s_b^2}{n} + \frac{z_A^2}{2 \cdot s_d^2} \cdot \left(\frac{s_b^4}{n-1} + \left(1 - \frac{1}{m_h}\right)^2 \frac{s_{w}^4}{N-n}\right) \]

\[ LME = z_{1-\alpha} \cdot S_{LoA} \]

4. Calculate LoA

\[ LoA_{lower} = \bar d - z_A \cdot s_d \]

\[ LoA_{upper} = \bar d + z_A \cdot s_d \]

5. Calculate LoA CI

\[ Lower \space CI = LoA_{lower} - LME \]

\[ Upper \space CI = LoA_{upper} + LME \]