Skip to contents

TDI estimation and inference

Source: vignettes/v2.TDI.Rmd
v2.TDI.Rmd

Introduction

The statistical analysis of agreement consists in the evaluation of the degree of similarity between measurements taken by different raters. Two raters would perfectly agree if we obtained the same measurements by both of them. However, in practice, the measurements obtained on the same subject or under the same conditions are generally not identical because of unavoidable errors that affect most measurement procedures. Hence, there is a need to quantify the level of agreement between raters.

For quantitative measurements a way to do this is via the total deviation index (TDI), an unscaled agreement index based on the difference between paired-measurements of both raters. Let \(X\) and \(Y\) be random variables that represent a pair of continuous measurements taken by two raters. Let their difference be \(D = X-Y\). The exact TDI is defined as the value of \(\kappa\) that yields \(\text{P}(|D|<\kappa)=p\), where \(0<p<1\), which is generally called \(\kappa_p\). This is, a boundary such that a certain proportion, \(p\), of these differences are within it.

To date, the existing methodologies for TDI estimation and inference include both parametric and non-parametric approaches. Although some methodologies exists for the estimation of the TDI in settings with more than two raters, this package focuses exclusively on the two-rater scenario.

In this vignette, we show how to use the function TDI(), which implements the estimation of the TDI and its corresponding \(100(1-\alpha)\%\) upper bound (UB), where \(\alpha\) is the significance level. We have implemented the parametric methods proposed by Choudhary (2007) and Escaramis et al. (2010), which are based on linear mixed models and assume the differences between paired-measurements to follow a normal distribution and linearity between the response and the effects; and the non-parametric methods of Choudhary (2010) and Perez-Jaume and Carrasco (2015), which are based on the estimation of quantiles of the absolute value of the paired-differences between raters.

Regarding the method of Perez-Jaume and Carrasco, in this package we have implemented an extended version of the method published in 2015 in Statistics in Medicine. We consider two non-parametric bootstrap approaches: the bootstrap of the within-subject paired-differences and the cluster bootstrap at subject level. We note that these approaches only differ when the study involves replicates, otherwise the two bootstrap approaches are equivalent. We also consider four strategies to estimate the UB: the ones based on the basic percentile and the normal distribution from Perez-Jaume and Carrasco (2015) and two additional ones based on empirical quantiles and the bias-corrected and accelerated (BC\(_a\)) method.

To illustrate the use of this function, it is applied to the same sample dataset used in the vignette for the function descriptive.plots().

We begin by loading the dataset and getting it ready for the analysis.

library(cccrm)
bp <- bpres[, c("ID", "SIS", "METODE", "NM", "EDAD", "SEXO")]
colnames(bp) <- c("id", "sys", "met", "rep", "age", "sex")
bp[, c("id", "met", "rep", "sex")] <- lapply(bp[, c("id", "met", "rep", "sex")], as.factor)

Note: In this study we identified five measurements in which the results from the two raters differ substantially. These readings originated from three specific individuals whose identification numbers were 24, 199 and 226. After discussions with clinical practitioners, these discrepancies were attributed to transcription errors (and not to the error introduced by the device itself). Since we could not recover the correct values of these measurements, we decided to discard all measurements from these patients. Consequently, the filtered dataset is used for all subsequent analyses and examples in this vignette. Eliminating these patients will leave empty levels in the id column, causing warnings to appear in subsequent function calls. While these warnings can be safely ignored, it is good practice to clean up unused levels using droplevels() before proceeding:

bp <- bp[!bp$id %in% c(24, 199, 226), ]
bp$id <- droplevels(bp$id)

The default call to the function TDI() is the following:

set.seed(17)
library(TDIagree)
x <- TDI(bp, sys, id, met, rep)

The function is called passing the name of the dataset (of class data.frame) in first instance. Afterwards, the column names of the quantitative measurements, the subject effect, the rater effect and the replicate effect (when there are no replicates the user should use rep = NA, the default value).

The object generated is of class tdi, which is a list with five components:

result: an object of class data.frame with the TDI estimates and UBs of the methods specified for every proportion (\(p\) in the definition above).

x$result
#>     p tdi.ch.p  ub.ch.p tdi.es.p  ub.es.p tdi.ch.np ub.ch.np tdi.pc.np
#> 1 0.9 15.29323 15.83926 15.30243 15.92305        15       16        15
#>   ub_p_db.pc.np ub_n_db.pc.np ub_e_db.pc.np ub_b_db.pc.np ub_p_cb.pc.np
#> 1      15.23333      15.93432            16            15            16
#>   ub_n_cb.pc.np ub_e_cb.pc.np ub_b_cb.pc.np
#> 1       16.1374            16            15

fitted.models: a named list with the fitted models of the parametric methods of Choudhary (2007) and Escaramis et al. (2010). For each parametric method, the list includes a character string describing the selected model structure (e.g., presence of interaction terms or heteroscedasticity) along with the corresponding fitted lme object.

x$fitted.models
#> $choudhary
#> $choudhary$structure
#> [1] "Without interaction term and with unequal residual variances"
#> 
#> $choudhary$model
#> Linear mixed-effects model fit by maximum likelihood
#>   Data: data.long 
#>   Log-likelihood: -5692.729
#>   Fixed: y ~ met - 1 
#>     met1     met2 
#> 133.1864 131.3163 
#> 
#> Random effects:
#>  Formula: ~1 | id
#>         (Intercept) Residual
#> StdDev:    19.57625 6.209965
#> 
#> Variance function:
#>  Structure: Different standard deviations per stratum
#>  Formula: ~1 | met 
#>  Parameter estimates:
#>        1        2 
#> 1.000000 1.072757 
#> Number of Observations: 1524
#> Number of Groups: 381 
#> 
#> 
#> $escaramis
#> $escaramis$structure
#> [1] "Without interaction term"
#> 
#> $escaramis$model
#> Linear mixed-effects model fit by REML
#>   Data: data.long 
#>   Log-restricted-likelihood: -5693.294
#>   Fixed: y ~ met 
#> (Intercept)        met2 
#>  133.186352   -1.870079 
#> 
#> Random effects:
#>  Formula: ~1 | id
#>         (Intercept) Residual
#> StdDev:    19.60937 6.443872
#> 
#> Number of Observations: 1524
#> Number of Groups: 381

params: a list with the values dec.est, dec.p, ub, method and alpha to be used in the method print.tdi() and in the method plot.tdi().

x$params
#> $dec.est
#> [1] 3
#> 
#> $dec.p
#> [1] 2
#> 
#> $ub
#> [1] TRUE
#> 
#> $method
#> [1] "choudhary p"              "escaramis et al."        
#> [3] "choudhary np"             "perez-jaume and carrasco"
#> 
#> $alpha
#> [1] 0.05

data.long: an object of class data.frame with columns y, id, met (and rep if it applies) with the values of the measurement, subject identifiers, rater (and replicate number if it applies) from the original data frame provided.

head(x$data.long)
#>     y id met rep
#> 1 146  1   1   1
#> 2 156  1   1   2
#> 3 166  1   2   1
#> 4 140  1   2   2
#> 5 100  2   1   1
#> 6 100  2   1   2

data.wide: an object of class data.frame with either:

  • columns id, y.met1, y.met2 (in the case of no replicates) with the measurements of each method,

  • columns id, y.met1.rep1,…, y.met1.rep\(m\), y.met2.rep1,…, y.met2.rep\(m\), with the measurements of each method,

where met1 and met2 correspond to the names of the first and second level of the column met in data, respectively. While rep1,…, rep\(m\) correspond to the names of the first,…, \(m\)-th level of the column rep in data, respectively.

head(x$data.wide)
#>    id y.1.1 y.1.2 y.2.1 y.2.2
#> 1   1   146   156   166   140
#> 5   2   100   100   107   115
#> 9   4   120   120   120   120
#> 13  5   112   114   114   113
#> 17  6   124   124   120   121
#> 21  8   120   120   122   112

By default, the arguments of the function estimate the TDI for all four methods implemented and for the proportion \(p=0.90\). The UBs at \(95\%\) are calculated with all the techniques considered, both by the bootstrap of the differences and the cluster bootstrap in the method of Perez-Jaume and Carrasco. The models of the parametric methods are automatically chosen minimizing the Akaike Information Criterion (AIC) and the tolerance considered for the method of Escaramís et al. is set to \(10^{-8}\).

General arguments

We can modify the default execution using several general arguments. For instance, the method argument specifies the method (or methods) to use for estimating the TDI and UB. The available options are: "Choudhary P" (Choudhary, 2007), "Escaramis et al." (Escaramis et al., 2010), "Choudhary NP" (Choudhary, 2010) and "Perez-Jaume and Carrasco" (Perez-Jaume and Carrasco, 2015). This argument is not case-sensitive and is passed to match.arg().

The package also contains a print.tdi() method for objects of class tdi, which formats and presents the results in a nice readable table.

x <- TDI(bp, sys, id, met, rep, method = c("es", "choudhary np"))
x
Parametric methods
Non-parametric methods
pp
Escaramis et al.
Choudhary NP
κ^p\hat{\kappa}_p UB95%(κp)\text{UB}_{95\%}(\kappa_p) κ^p\hat{\kappa}_p UB95%(κp)\text{UB}_{95\%}(\kappa_p)
0.90 15.302 15.923 15.000 16.000
pp, proportion; κ^p\hat{\kappa}_p, TDI estimate; UB\text{UB}, upper bound

This means, that 90% of the absolute differences between the sphygmomanometer and the OMRON 711 measurements fall within 15.302 mmHg according to the method of Escaramís et al., and within 15.000 mmHg according to the non-parametric method of Choudhary. Furthermore, with 95% confidence, we can assert that 90% of the absolute differences do not exceed 15.923 mmHg (Escaramís et al.) or 16.000 mmHg (Choudhary NP).

We can also choose the proportion (or proportions) to estimate the TDI.

x <- TDI(bp, sys, id, met, rep, method = c("es", "choudhary np"),
         p = c(0.85, 0.90, 0.95))
x
Parametric methods
Non-parametric methods
pp
Escaramis et al.
Choudhary NP
κ^p\hat{\kappa}_p UB95%(κp)\text{UB}_{95\%}(\kappa_p) κ^p\hat{\kappa}_p UB95%(κp)\text{UB}_{95\%}(\kappa_p)
0.85 13.393 13.962 12.000 13.000
0.90 15.302 15.923 15.000 16.000
0.95 18.231 18.937 19.000 20.000
pp, proportion; κ^p\hat{\kappa}_p, TDI estimate; UB\text{UB}, upper bound

If we do not want to compute the UBs we can call ub = FALSE.

x <- TDI(bp, sys, id, met, rep,
         p = c(0.85, 0.90, 0.95),
         ub = FALSE)
x
Parametric methods
Non-parametric methods
pp
Choudhary P
Escaramis et al.
Choudhary NP
Perez-Jaume and Carrasco
κ^p\hat{\kappa}_p κ^p\hat{\kappa}_p κ^p\hat{\kappa}_p κ^p\hat{\kappa}_p
0.85 13.385 13.393 12.000 12.000
0.90 15.293 15.302 15.000 15.000
0.95 18.220 18.231 19.000 19.000
pp, proportion; κ^p\hat{\kappa}_p, TDI estimate

We can modify the significance level with alpha.

x <- TDI(bp, sys, id, met, rep, method = c("es", "choudhary np"),
         p = c(0.85, 0.90, 0.95), alpha = 0.1)
x
Parametric methods
Non-parametric methods
pp
Escaramis et al.
Choudhary NP
κ^p\hat{\kappa}_p UB90%(κp)\text{UB}_{90\%}(\kappa_p) κ^p\hat{\kappa}_p UB90%(κp)\text{UB}_{90\%}(\kappa_p)
0.85 13.393 13.835 12.000 13.000
0.90 15.302 15.784 15.000 15.000
0.95 18.231 18.778 19.000 20.000
pp, proportion; κ^p\hat{\kappa}_p, TDI estimate; UB\text{UB}, upper bound

Finally, there are two arguments to modify the number of decimals displayed on the results table: dec.p for the proportion and dec.est for the TDI estimates and UBs.

x <- TDI(bp, sys, id, met, rep, method = c("es", "choudhary np"),
         p = c(0.85, 0.90, 0.95), dec.p = 3, dec.est = 2)
x
Parametric methods
Non-parametric methods
pp
Escaramis et al.
Choudhary NP
κ^p\hat{\kappa}_p UB95%(κp)\text{UB}_{95\%}(\kappa_p) κ^p\hat{\kappa}_p UB95%(κp)\text{UB}_{95\%}(\kappa_p)
0.850 13.39 13.96 12.00 13.00
0.900 15.30 15.92 15.00 16.00
0.950 18.23 18.94 19.00 20.00
pp, proportion; κ^p\hat{\kappa}_p, TDI estimate; UB\text{UB}, upper bound

Method-specific arguments

We can estimate the TDI and the UB with only one method and modify its specific arguments. There is no need to call only one of the methods to modify its specific arguments, but we will do so below to illustrate it more schematically.

Parametric method of Choudhary

By calling method = "choudhary p" the TDI and the UB are computed only via the parametric method of Choudhary.

x <- TDI(bp, sys, id, met, rep, method = "choudhary p")
x
Parametric methods
pp
Choudhary P
κ^p\hat{\kappa}_p UB95%(κp)\text{UB}_{95\%}(\kappa_p)
0.90 15.293 15.839
pp, proportion; κ^p\hat{\kappa}_p, TDI estimate; UB\text{UB}, upper bound

In this case, the best model according to the AIC was the one with residual homoscedasticity between raters.

x$fitted.models
#> $choudhary
#> $choudhary$structure
#> [1] "Without interaction term and with unequal residual variances"
#> 
#> $choudhary$model
#> Linear mixed-effects model fit by maximum likelihood
#>   Data: data.long 
#>   Log-likelihood: -5692.729
#>   Fixed: y ~ met - 1 
#>     met1     met2 
#> 133.1864 131.3163 
#> 
#> Random effects:
#>  Formula: ~1 | id
#>         (Intercept) Residual
#> StdDev:    19.57625 6.209965
#> 
#> Variance function:
#>  Structure: Different standard deviations per stratum
#>  Formula: ~1 | met 
#>  Parameter estimates:
#>        1        2 
#> 1.000000 1.072757 
#> Number of Observations: 1524
#> Number of Groups: 381 
#> 
#> 
#> $escaramis
#> [1] "Method of Escaramis et al. not called"

If we want to manually choose the model to be fitted, we can call choose.model.ch.p = FALSE and indicate whether there is residual homoscedasticity between raters with the argument var.equal.

x <- TDI(bp, sys, id, met, rep, method = "choudhary p", 
         choose.model.ch.p = FALSE, var.equal = FALSE)
x
Parametric methods
pp
Choudhary P
κ^p\hat{\kappa}_p UB95%(κp)\text{UB}_{95\%}(\kappa_p)
0.90 15.293 15.839
pp, proportion; κ^p\hat{\kappa}_p, TDI estimate; UB\text{UB}, upper bound

We can see now that the model fitted is the one without residual homoscedasticity between raters.

x$fitted.models
#> $choudhary
#> $choudhary$structure
#> [1] "Without interaction term and with unequal residual variances"
#> 
#> $choudhary$model
#> Linear mixed-effects model fit by maximum likelihood
#>   Data: data.long 
#>   Log-likelihood: -5692.729
#>   Fixed: y ~ met - 1 
#>     met1     met2 
#> 133.1864 131.3163 
#> 
#> Random effects:
#>  Formula: ~1 | id
#>         (Intercept) Residual
#> StdDev:    19.57625 6.209965
#> 
#> Variance function:
#>  Structure: Different standard deviations per stratum
#>  Formula: ~1 | met 
#>  Parameter estimates:
#>        1        2 
#> 1.000000 1.072757 
#> Number of Observations: 1524
#> Number of Groups: 381 
#> 
#> 
#> $escaramis
#> [1] "Method of Escaramis et al. not called"

Parametric method of Escaramís et al.

Analogously, we can estimate the TDI and UB with only the parametric method of Escaramís et al.

x <- TDI(bp, sys, id, met, rep, method = "es")
x
Parametric methods
pp
Escaramis et al.
κ^p\hat{\kappa}_p UB95%(κp)\text{UB}_{95\%}(\kappa_p)
0.90 15.302 15.923
pp, proportion; κ^p\hat{\kappa}_p, TDI estimate; UB\text{UB}, upper bound

The best model according to the AIC is the one without interaction between subjects and methods.

x$fitted.models
#> $choudhary
#> [1] "Parametric method of Choudhary not called"
#> 
#> $escaramis
#> $escaramis$structure
#> [1] "Without interaction term"
#> 
#> $escaramis$model
#> Linear mixed-effects model fit by REML
#>   Data: data.long 
#>   Log-restricted-likelihood: -5693.294
#>   Fixed: y ~ met 
#> (Intercept)        met2 
#>  133.186352   -1.870079 
#> 
#> Random effects:
#>  Formula: ~1 | id
#>         (Intercept) Residual
#> StdDev:    19.60937 6.443872
#> 
#> Number of Observations: 1524
#> Number of Groups: 381

If we want to manually add the interaction we can do it as follows:

x <- TDI(bp, sys, id, met, rep, method = "es",
         choose.model.es = FALSE, int = TRUE)
x
Parametric methods
pp
Escaramis et al.
κ^p\hat{\kappa}_p UB95%(κp)\text{UB}_{95\%}(\kappa_p)
0.90 15.302 16.016
pp, proportion; κ^p\hat{\kappa}_p, TDI estimate; UB\text{UB}, upper bound 
x$fitted.models
#> $choudhary
#> [1] "Parametric method of Choudhary not called"
#> 
#> $escaramis
#> $escaramis$structure
#> [1] "With interaction term"
#> 
#> $escaramis$model
#> Linear mixed-effects model fit by REML
#>   Data: data.long 
#>   Log-restricted-likelihood: -5693.294
#>   Fixed: y ~ met 
#> (Intercept)        met2 
#>  133.186352   -1.870079 
#> 
#> Random effects:
#>  Composite Structure: Blocked
#> 
#>  Block 1: (Intercept)
#>  Formula: ~1 | id
#>         (Intercept)
#> StdDev:    19.60937
#> 
#>  Block 2: met1, met2
#>  Formula: ~-1 + met | id
#>  Structure: Multiple of an Identity
#>                 met1         met2 Residual
#> StdDev: 0.0005289496 0.0005289496 6.443872
#> 
#> Number of Observations: 1524
#> Number of Groups: 381

In this method, we can also change the tolerance with the argument tol and we can enter other covariates to the model with the parameter add.es, such as the sex or the age of the subjects.

x <- TDI(bp, sys, id, met, rep, method = "es",
         tol = 10^(-7), add.es = c(sex, age))
x
Parametric methods
pp
Escaramis et al.
κ^p\hat{\kappa}_p UB95%(κp)\text{UB}_{95\%}(\kappa_p)
0.90 15.302 15.923
pp, proportion; κ^p\hat{\kappa}_p, TDI estimate; UB\text{UB}, upper bound

The inclusion of covariates in the model did not modify the parameter estimates used to calculate the TDI, this is, the fixed-effect of the rater and the error variability estimates, and therefore the TDI estimates as well as their \(95\%\) UB remain the same.

x$fitted.models
#> $choudhary
#> [1] "Parametric method of Choudhary not called"
#> 
#> $escaramis
#> $escaramis$structure
#> [1] "Without interaction term"
#> 
#> $escaramis$model
#> Linear mixed-effects model fit by REML
#>   Data: data.long 
#>   Log-restricted-likelihood: -5597.299
#>   Fixed: formula_obj 
#> (Intercept)        met2        sex2         age 
#>  95.7364043  -1.8700787  -8.8516134   0.8501193 
#> 
#> Random effects:
#>  Formula: ~1 | id
#>         (Intercept) Residual
#> StdDev:    15.11408 6.443872
#> 
#> Number of Observations: 1524
#> Number of Groups: 381

Non-parametric method of Choudhary

Estimating the TDI and UB with the non-parametric method of Choudhary:

x <- TDI(bp, sys, id, met, rep, method = "choudhary np")
x
Non-parametric methods
pp
Choudhary NP
κ^p\hat{\kappa}_p UB95%(κp)\text{UB}_{95\%}(\kappa_p)
0.90 15.000 16.000
pp, proportion; κ^p\hat{\kappa}_p, TDI estimate; UB\text{UB}, upper bound

There are no specific arguments for this method.

Non-parametric method of Perez-Jaume and Carrasco

With the default setting of the non-parametric method of Perez-Jaume and Carrasco we obtain the TDI estimate and the UB calculated using eight different techniques: the combination of the two non-parametric bootstrap approaches and the four strategies mentioned in the introduction.

x <- TDI(bp, sys, id, met, rep, method = "pe")
x
Non-parametric methods
Perez-Jaume and Carrasco
pp κ^p\hat{\kappa}_p
Differences bootstrap
Cluster bootstrap
UB95%p(κp)\text{UB}_{95\%}^{\text{p}}(\kappa_p) UB95%n(κp)\text{UB}_{95\%}^{\text{n}}(\kappa_p) UB95%e(κp)\text{UB}_{95\%}^{\text{e}}(\kappa_p) UB95%BCa(κp)\text{UB}_{95\%}^{\text{BC}_a}(\kappa_p) UB95%p(κp)\text{UB}_{95\%}^{\text{p}}(\kappa_p) UB95%n(κp)\text{UB}_{95\%}^{\text{n}}(\kappa_p) UB95%e(κp)\text{UB}_{95\%}^{\text{e}}(\kappa_p) UB95%BCa(κp)\text{UB}_{95\%}^{\text{BC}_a}(\kappa_p)
0.90 15.000 15.233 15.932 16.000 15.000 16.000 16.120 16.000 15.000
pp, proportion; κ^p\hat{\kappa}_p, TDI estimate; UB\text{UB}, upper bound; p\text{p}, strategy based on percentiles; n\text{n}, strategy based on the normal distribution; e\text{e}, strategy based on the empirical quantiles; BCa\text{BC}_a, strategy based on the bias-corrected and accelerated method

We can choose the bootstrap type with the argument boot.type; the algorithm for quantile computations is chosen with the argument type; and the number of bootstrap replicates with the argument R.

x <- TDI(bp, sys, id, met, rep, method = "pe",
         boot.type = "cl", R = 1000)
x
Non-parametric methods
Perez-Jaume and Carrasco
pp κ^p\hat{\kappa}_p
Cluster bootstrap
UB95%p(κp)\text{UB}_{95\%}^{\text{p}}(\kappa_p) UB95%n(κp)\text{UB}_{95\%}^{\text{n}}(\kappa_p) UB95%e(κp)\text{UB}_{95\%}^{\text{e}}(\kappa_p) UB95%BCa(κp)\text{UB}_{95\%}^{\text{BC}_a}(\kappa_p)
0.90 15.000 16.000 16.152 16.000 15.000
pp, proportion; κ^p\hat{\kappa}_p, TDI estimate; UB\text{UB}, upper bound; p\text{p}, strategy based on percentiles; n\text{n}, strategy based on the normal distribution; e\text{e}, strategy based on the empirical quantiles; BCa\text{BC}_a, strategy based on the bias-corrected and accelerated method

Considerations

Choosing the appropriate method for estimating the TDI depends heavily on the distributional characteristics of the data. Performing an initial exploratory data analysis is essential to inspect normality, skewness, and potential outliers before deciding on the estimation approach. In this package, we have implemented the function descriptive.plots() function to facilitate a visual exploration of the data.

Under normality, parametric methods generally outperform non-parametric ones by providing better accuracy, consistency and efficiency than the non-parametric methods. Therefore, in situations where the assumptions of the parametric methods hold, we recommend assessing agreement using either the parametric method of Choudhary or the method of Escaramís et al.

In contrast, when dealing with skewed data, parametric methods are noticeably affected by the presence of outliers in the paired-differences, whereas the non-parametric proposals are thought to be more robust in this regard. Thus, when working with skewed data, we recommend assessing agreement via the non-parametric method of Choudhary or the method of Perez-Jaume and Carrasco.

References

Lin LI. Total deviation index for measuring individual agreement with applications in laboratory performance and bioequivalence. Statistics in Medicine. 2000;19(2):255-270.

Choudhary PK. A tolerance interval approach for assessment of agreement with left censored data. Journal of Biopharmaceutical Statistics. 2007;17(4):583-94.

Escaramís G, Ascaso C, Carrasco JL. The total deviation index estimated by tolerance intervals to evaluate the concordance of measurement devices. BMC Medical Research Methodology. 2010;10:1-12.

Choudhary PK. A unified approach for nonparametric evaluation of agreement in method comparison studies. The International Journal of Biostatistics. 2010;6(1).

Perez-Jaume S, Carrasco JL. A non-parametric approach to estimate the total deviation index for non-normal data. Statistics in Medicine. 2015;34(25):3318-35.