Agree or Tolerate?

Francq, Berger, and Boachie (2020) pose this question in a paper published in Statistics in Medicine. Traditionally, those working in medicine or physiology have defaulted to calculating some form of “limits of agreement” that Bland and Altman (1986) recommended in their seminal paper. The recommendation by Bland and Altman (1986) was only an approximation, and that has undergone many modifications (e.g., Bland and Altman 1999; Zou 2011) to improve the accuracy of the agreement interval and their associated confidence intervals. Meanwhile, the field of industrial statistics has focused on calculating tolerance limits. There are R packages, such as the tolerance R package available on CRAN, that are entirely dedicated to calculating tolerance limits. It is important to note that tolerance is not limited to the normal distribution and can be applied to other distributions (please see the tolerance package by Young (2010) for more details). However, in the agreement studies typically seen in medicine, tolerance limits may be a more accurate way of determining whether two measurements are adequately close to one another.

To quote Francq, Berger, and Boachie (2020):

In terms of terminology, tolerance means, in this context, that some difference between the methods is tolerated (the measurements are still comparable in practice). Furthermore, the tolerance interval is exact and therefore more appropriate than the agreement interval.

In this package, tolerance limits refer to the “tolerance” associated with the prediction interval for the difference between 2 measurements. The calculative approach (detailed much further below) involves calculating the prediction interval for the difference between two methods (i.e., an estimate of an interval in which a future observation will fall, with a certain probability, given what has already been observed) and then calculating the confidence in the interval (i.e., tolerance). Therefore, if we want a 95% prediction interval with 95% tolerance limits, we are calculating the interval in which 95% of future observations (prediction) with only a 5% probability (1-tolerance) the “true” prediction interval limits are greater/less than the upper/lower tolerance limits.

Personally, I find the use of prediction intervals and tolerance limits more attractive for 2 reasons: 1) the coverage of the prediction intervals and their tolerance limits is often better than confidence intervals for agreement limits, and 2) the interpretation of the tolerance limits is much clearer. For a greater discussion of this topic, please see the manuscript by Francq, Berger, and Boachie (2020) and check out their R package BivRegBLS.

Tolerance

In SimplyAgree we utilize a generalized least square (GLS) model to estimate the prediction interval and tolerance limits. The function uses the gls function from the nlme R package to build the model. This allows a very flexible approach for estimating prediction intervals and tolerance limits.

We can use the tolerance_limits function demonstrate the basic calculations. In this example (below), we use the temps data set to measure the differences between esophageal and rectal temperatures between different times of day (tod) and controlling for the intra-subject correlation.

tolerance_limit(
  data = temps,
  x = "trec_pre", # First measure
  y = "teso_pre", # Second measure
  id = "id", # Subject ID
  condition = "tod", # Identify condition that may affect differences
  cor_type = "sym" # Set correlation structure as Compound Symmetry
)
#> Agreement between Measures (Difference: x-y)
#> 95% Prediction Interval with 95% Tolerance Limits
#> 
#>  Condition   Bias          Bias CI Prediction Interval  Tolerance Limits
#>         AM 0.1537 [0.0595, 0.2478]   [-0.2919, 0.5993] [-0.4982, 0.8055]
#>         PM 0.2280 [0.1342, 0.3218]   [-0.3162, 0.7722] [-0.6982, 1.1542]

test1 = tolerance_limit(data = temps,
                        x = "teso_pre",
                        y = "trec_pre",
                        id = "id",
                        condition = "tod")

test1
#> Agreement between Measures (Difference: x-y)
#> 95% Prediction Interval with 95% Tolerance Limits
#> 
#>  Condition    Bias            Bias CI Prediction Interval  Tolerance Limits
#>         AM -0.1537 [-0.2478, -0.0595]   [-0.5993, 0.2919] [-0.8055, 0.4982]
#>         PM -0.2280 [-0.3218, -0.1342]   [-0.7722, 0.3162] [-1.1542, 0.6982]

Agreement

The agreement limit calculations in SimplyAgree are a tad more constrained than the tolerance limit. There are only 3 types of calculations that can be made: simple, replicate, and nested. The simple calculation each pair of observations (x and y) are independent; meaning that each pair represents one subject/participant. Sometimes there are multiple measurements taken within subjects when comparing two measurements tools. In some cases the true underlying value will not be expected to vary (i.e., replicates or “reps”), or multiple measurements may be taken within an individual and these values are expected to vary (i.e., nested).

The agreement_limit function, unlike the other agreement functions in the package (i.e., agree_test, agree_reps, and agree_nest), allows users to make any of the three calculations all-in-one function. Further, agreement_limit assumes you are only interested in the outermost confidence interval of the limits of agreement and therefore one-tailed confidence intervals are only calculated for the LoA (i.e., there is no TOST argument).

# Calc. LoA
a1 = agreement_limit(data = reps,
                     x = "x",
                     y = "y")
# print
a1
#> MOVER Limits of Agreement (LoA)
#> 95% LoA @ 5% Alpha-Level
#> Independent Data Points
#> 
#>    Bias           Bias CI Lower LoA Upper LoA            LoA CI
#>  0.4383 [-0.1669, 1.0436]    -1.947     2.824 [-3.0117, 3.8884]
#> 
#> SD of Differences = 1.217

a2 = agreement_limit(x = "x",
                y = "y",
                id = "id",
                data = reps,
                data_type = "reps",
                agree.level = .8) 

a2
#> MOVER Limits of Agreement (LoA)
#> 80% LoA @ 5% Alpha-Level
#> Data with Replicates
#> 
#>    Bias           Bias CI Lower LoA Upper LoA           LoA CI
#>  0.7152 [-1.5287, 2.9591]    -1.212     2.642 [-4.797, 6.2274]
#> 
#> SD of Differences = 1.5036

a3 = agreement_limit(x = "x",
                y = "y",
                id = "id",
                data = reps,
                data_type = "nest",
                loa_calc = "mover",
                agree.level = .95)
a3
#> MOVER Limits of Agreement (LoA)
#> 95% LoA @ 5% Alpha-Level
#> Nested Data
#> 
#>    Bias           Bias CI Lower LoA Upper LoA            LoA CI
#>  0.7046 [-1.5572, 2.9664]    -2.153     3.562 [-7.4979, 8.9071]
#> 
#> SD of Differences = 1.4581

How to Test a Hypothesis

Unlike the “test” functions, agreement_limits and tolerance_limits do not test a hypothesis agreement, but only estimate the limits. It is up to the user to determine if the limits are within a maximal allowable difference. However, the results a maximal allowable difference can be visualized against the tolerance/agreement limits using the delta argument for the plot method.

res1 = tolerance_limit(
  data = temps,
  x = "trec_pre", # First measure
  y = "teso_pre", # Second measure
  id = "id", # Subject ID
  condition = "tod", # Identify condition that may affect differences
  cor_type = "sym" # Set correlation structure as Compound Symmetry
)
plot(res1, delta = .25) # Set maximal allowable difference to .25 units

Checking Assumptions

The assumptions of normality, heteroscedasticity, and proportional bias can all be checked using the check method for either agreement or tolerance limits.

The function will provide 3 plots: Q-Q normality plot, standardized residuals plot, and proportional bias plot.

All 3 plots will have a statistical test in the bottom right corner². The Shapiro-Wilk test is included for the normality plot, the Bagan-Preusch test for heterogeneity, and the test for linear slope on the residuals plot. Please note that there is no formal test of proportional bias for the tolerance limits, but a plot is still included for visual checks.

test_agree = agreement_limit(x = "x",
                             y = "y",
                             data = reps)

check(test_agree)

test_tol = tolerance_limit(x = "x",
                           y = "y",
                           data = reps)

check(test_tol)

Proportional Bias

As the check plots for a1 show, proportional bias can sometimes occur. In these cases Bland and Altman (1999) recommended adjusting the bias and LoA for the proportional bias. This is simply done by include a slope for the average of both measurements (i.e, using an intercept + slope model rather than intercept only model).

For either agreement_limit or tolerance_limit functions, this can be accomplished with the prop_bias argument. When this is set to TRUE, then the proportional bias adjusted model is utilized. Plots and checks of the data should always be inspected.

test_tol = tolerance_limit(x = "x",
                           y = "y",
                           data = reps,
                           prop_bias = TRUE)
print(test_tol)
#> Agreement between Measures (Difference: x-y)
#> 95% Prediction Interval with 95% Tolerance Limits
#> 
#>  Average of Both Methods    Bias           Bias CI Prediction Interval
#>                    3.905 -0.4670 [-1.3842, 0.4502]   [-2.8876, 1.9537]
#>                    5.240  0.3513 [-0.1816, 0.8842]    [-1.9513, 2.654]
#>                    7.395  1.6723  [0.5218, 2.8228]    [-0.846, 4.1906]
#>   Tolerance Limits
#>  [-3.6396, 2.7057]
#>  [-2.6667, 3.3694]
#>  [-1.6284, 4.9729]

# See effect of proportional bias on limits
plot(test_tol)

# Confirm its effects in proportional bias check plot (should be horizontal now)
check(test_tol)

Log transformation

Sometimes a log transformation may be a useful way of “normalizing” the data. Most often this is done when the error is proportional to the mean. The interpretation is also easy because the differences (when back transformed) can be interpreted as ratios. The log transformation (natural base) can be accomplished with the log_tf argument.

tolerance_limit(
  data = temps,
  log_tf = TRUE, # natural log transformation of responses
  x = "trec_pre", # First measure
  y = "teso_pre", # Second measure
  id = "id", # Subject ID
  condition = "tod", # Identify condition that may affect differences
  cor_type = "sym" # Set correlation structure as Compound Symmetry
)
#> Agreement between Measures (Ratio: x/y)
#> 95% Prediction Interval with 95% Tolerance Limits
#> 
#>  Condition  Bias          Bias CI Prediction Interval Tolerance Limits
#>         AM 1.004 [1.0016, 1.0068]    [0.9921, 1.0165] [0.9865, 1.0222]
#>         PM 1.006 [1.0036, 1.0088]    [0.9913, 1.0213]  [0.981, 1.0321]

If you prefer to interpret the differences as a percentage difference, you can do this by setting the log_tf_display argument to “sympercent”. This stands for the “symmetric percentage difference” which is the log transformed differences between the measures, \(s\% = (log(x)-log(y) ) \cdot 100\% = log(x/y) \cdot 100\%\), and can be interpreted as a percentage difference between the two paired measures.

tolerance_limit(
  data = temps,
  log_tf = TRUE, # natural log transformation of responses
  log_tf_display = "sympercent", # display results as sympercent
  x = "trec_pre", # First measure
  y = "teso_pre", # Second measure
  id = "id", # Subject ID
  condition = "tod", # Identify condition that may affect differences
  cor_type = "sym" # Set correlation structure as Compound Symmetry
)
#> Sympercent Difference between Methods (s%)
#> 95% Prediction Interval with 95% Tolerance Limits
#> 
#>  Condition   Bias          Bias CI Prediction Interval  Tolerance Limits
#>         AM 0.4188 [0.1615, 0.6761]   [-0.7959, 1.6335] [-1.3562, 2.1938]
#>         PM 0.6184 [0.3622, 0.8747]   [-0.8705, 2.1074] [-1.9196, 3.1564]

Visualizing “Big” Data

Sometimes there may be a lot of data and individual points of data on mean difference visualization may be less than ideal. In order to change the plots from showing the individual data points we can modify the geom argument.

set.seed(81346)
x = rnorm(750, 100, 10)
diff = rnorm(750, 0, 1)
y = x + diff

df = data.frame(x = x,
                y = y)

a1 = agreement_limit(data = df,
                     x = "x",
                     y = "y",
                     agree.level = .95)

plot(a1,
     geom = "geom_point")

plot(a1,
     geom = "geom_bin2d")

plot(a1,
     geom = "geom_density_2d")

plot(a1,
     geom = "geom_density_2d_filled")

plot(a1,
     geom = "stat_density_2d")