---
title: "TDI estimation and inference"
output: 
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{TDI estimation and inference}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
<style>
body {
text-align: justify}
</style>

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

## 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()`](v1.descriptive.plots.html).

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

```{r}
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:

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

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

```{r}
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).

```{r}
x$result
```

`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.

```{r}
x$fitted.models
```

`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()`.

```{r}
x$params
```

`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.

```{r}
head(x$data.long)
```

`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.

```{r}
head(x$data.wide)
```

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.

```{r}
x <- TDI(bp, sys, id, met, rep, method = c("es", "choudhary np"))
x
```

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.

```{r}
x <- TDI(bp, sys, id, met, rep, method = c("es", "choudhary np"),
         p = c(0.85, 0.90, 0.95))
x
```

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

```{r}
x <- TDI(bp, sys, id, met, rep,
         p = c(0.85, 0.90, 0.95),
         ub = FALSE)
x
```

We can modify the significance level with `alpha`.

```{r}
x <- TDI(bp, sys, id, met, rep, method = c("es", "choudhary np"),
         p = c(0.85, 0.90, 0.95), alpha = 0.1)
x
```

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.

```{r}
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
```


## 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. 

```{r}
x <- TDI(bp, sys, id, met, rep, method = "choudhary p")
x
```

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

```{r}
x$fitted.models
```

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`.

```{r}
x <- TDI(bp, sys, id, met, rep, method = "choudhary p", 
         choose.model.ch.p = FALSE, var.equal = FALSE)
x
```

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

```{r}
x$fitted.models
```

### 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.* 

```{r}
x <- TDI(bp, sys, id, met, rep, method = "es")
x
```

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

```{r}
x$fitted.models
```

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

```{r}
x <- TDI(bp, sys, id, met, rep, method = "es",
         choose.model.es = FALSE, int = TRUE)
x
```

```{r}
x$fitted.models
```

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.

```{r}
x <- TDI(bp, sys, id, met, rep, method = "es",
         tol = 10^(-7), add.es = c(sex, age))
x
```

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.

```{r}
x$fitted.models
```

### Non-parametric method of Choudhary

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

```{r}
x <- TDI(bp, sys, id, met, rep, method = "choudhary np")
x
```

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.

```{r}
x <- TDI(bp, sys, id, met, rep, method = "pe")
x
```

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`.

```{r}
x <- TDI(bp, sys, id, met, rep, method = "pe",
         boot.type = "cl", R = 1000)
x
```


## 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.