---
title: "Re-analysis of an Agreement Study"
output: rmarkdown::html_vignette
bibliography: refs.bib
vignette: >
  %\VignetteIndexEntry{Re-analysis of an Agreement Study}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

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

```{r setup, message=FALSE,warning=FALSE}
library(SimplyAgree)
#library(tidyverse)
library(dplyr)
library(tidyr)
library(ggplot2)
library(magrittr)
data("temps")
df_temps = temps
```

# Re-analysis of a Previous Study of Agreement

In the study by @ravanelli2020change, they attempted to estimate the effect of varying the time of day (AM or PM) on the measurement of thermoregulatory variables (e.g., rectal and esophageal temperature). In total, participants completed 6 separate trials wherein these variables were measured. While this is a robust study of these variables the analyses focused on ANOVAs and t-tests to determine whether or not the time-of-day (AM or PM). This poses a problem because 1) they were trying to test for equivalence and 2) this is a study of *agreement* not *differences* (See @Lin1989). Due to the latter point, the use of t-test or ANOVAs (F-tests) is rather inappropriate since they provide an answer to different, albeit related, question.

Instead, the authors could test their hypotheses by using tools that estimate the absolute *agreement* between the AM and PM sessions within each condition. This is rather complicated because we have multiple measurement within each participant. However, with the tools included in `SimplyAgree`^[ @cccrm is another package to check out] I believe we can get closer to the right answer.

In order to understand the underlying processes of these functions and procedures I highly recommend reading the statistical literature that documents methods within these functions. For the `cccrm` package please see the work by @carrasco2003, @carrasco2009, and @carrasco2013. The `tolerance_limit` function was inspired by the work of @francq2020tolerate which documented how to implement tolerance limits to measure agreement. 

# Concordance

An easy approach to measuring agreement between 2 conditions or measurement tools is through the concordance correlation coefficient (CCC). The CCC essentially provides a single coefficient (values between 0 and 1) that provides an estimate to how closely one measurement is to another. It is a type of intraclass correlation coefficient that takes into account the mean difference between two measurements. In other words, if we were to draw a line of identity on a graph and plot two measurements (X & Y), the closer those points are to the line of identity the higher the CCC (and vice versa). Please see the `cccrm` package for more details.

```{r, fig.cap="Example of the Line of Identity"}
qplot(1,1) + geom_abline(intercept = 0, slope = 1)
```

In the following sections, let us see how well esophageal and rectal temperature are in agreement after exercising in the heat for 1 hour at differing conditions.

## Rectal Temperature

We can visualize the concordance between the two different types of measurements and the respective time-of-day and conditions. From the plot we can see there is clear bias in the raw post exercise values (higher in the PM), but even when "correcting for baseline differences" by calculating the differences scores we can see a higher degree of disagreement between the two conditions. 

```{r pltsrec,fig.cap="Concordance Plots of Rectal Temperature",echo=FALSE,fig.width=8.5,fig.height=5}
df_rec.delta = df_temps %>%
  mutate(id_spec = paste0(id,"_",trial_condition)) %>%
  select(id,id_spec,trec_delta,tod,trial_condition) %>%
  pivot_wider(id_cols = c(id,id_spec,trial_condition),
              names_from = tod,
              values_from = trec_delta) %>%
  mutate(diff = PM - AM,
         Average = (AM + PM)/2)

df_rec.post = df_temps %>%
  mutate(id_spec = paste0(id,"_",trial_condition)) %>%
  select(id,id_spec,trec_post,tod,trial_condition) %>%
  pivot_wider(id_cols = c(id,id_spec,trial_condition),
              names_from = tod,
              values_from = trec_post) %>%
  mutate(diff = PM - AM,
         Average = (AM + PM)/2)

p_rec.delta <- ggplot(df_rec.delta, aes(x=AM, y=PM)) +
  stat_smooth(aes(color=trial_condition), 
              method = "lm",
              formula = 'y ~ x',
              fullrange = TRUE,
              geom='line', alpha=0.5, se=FALSE)+
  scale_x_continuous("AM - Trec (delta)",
                     limits = c(.15,1.1))+
  scale_y_continuous("PM - Trec (delta)",
                     limits = c(.15,1.1))+
  geom_abline(intercept = 0, slope = 1, size = .85,
              alpha = .75) +
  geom_point(aes(color=trial_condition),
             alpha = .75) +
  theme_classic() +
  theme(legend.position = "bottom") +
  scale_color_viridis_d()

p_rec.post <- ggplot(df_rec.post, aes(x=AM, y=PM)) +
  stat_smooth(aes(color=trial_condition), 
              method = "lm",
              formula = 'y ~ x',
              fullrange = TRUE,
              geom='line', alpha=0.75, se=FALSE)+
  scale_x_continuous("AM - Trec (post)",
                     limits = c(36.4,38.5))+
  scale_y_continuous("PM - Trec (post)",
                     limits = c(36.4,38.5))+
  geom_abline(intercept = 0, slope = 1, size = .85,
              alpha = .75) +
  geom_point(aes(color=trial_condition),
             alpha = .75) +
  theme_classic() +
  theme(legend.position = "bottom") +
  scale_color_viridis_d()
p_rec.post
p_rec.delta
```

# Esophageal Temperature

```{r pltseso,fig.cap="Concordance Plots of Esophageal Temperature",echo=FALSE,fig.width=8.5,fig.height=5}
df_eso.delta = df_temps %>%
  mutate(id_spec = paste0(id,"_",trial_condition)) %>%
  select(id,id_spec,teso_delta,tod,trial_condition) %>%
  pivot_wider(id_cols = c(id,id_spec,trial_condition),
              names_from = tod,
              values_from = teso_delta) %>%
  mutate(diff = PM - AM,
         Average = (AM + PM)/2)

df_eso.post = df_temps %>%
  mutate(id_spec = paste0(id,"_",trial_condition)) %>%
  select(id,id_spec,teso_post,tod,trial_condition) %>%
  pivot_wider(id_cols = c(id,id_spec,trial_condition),
              names_from = tod,
              values_from = teso_post) %>%
  mutate(diff = PM - AM,
         Average = (AM + PM)/2)

p_eso.delta <- ggplot(df_eso.delta, aes(x=AM, y=PM)) +
  stat_smooth(aes(color=trial_condition), 
              method = "lm",
              formula = 'y ~ x',
              fullrange = TRUE,
              geom='line', alpha=0.5, se=FALSE)+
  scale_x_continuous("AM - Teso (delta)",
                     limits = c(.15,1.1))+
  scale_y_continuous("PM - Teso (delta)",
                     limits = c(.15,1.1))+
  geom_abline(intercept = 0, slope = 1, size = .85,
              alpha = .75) +
  geom_point(aes(color=trial_condition),
             alpha = .75) +
  theme_classic() +
  theme(legend.position = "bottom") +
  scale_color_viridis_d()

p_eso.post <- ggplot(df_eso.post, aes(x=AM, y=PM)) +
  stat_smooth(aes(color=trial_condition), 
              method = "lm",
              formula = 'y ~ x',
              fullrange = TRUE,
              geom='line', alpha=0.75, se=FALSE)+
  scale_x_continuous("AM - Teso (post)",
                     limits = c(36.4,38.5))+
  scale_y_continuous("PM - Teso (post)",
                     limits = c(36.4,38.5))+
  geom_abline(intercept = 0, slope = 1, size = .85,
              alpha = .75) +
  geom_point(aes(color=trial_condition),
             alpha = .75) +
  theme_classic() +
  theme(legend.position = "bottom") +
  scale_color_viridis_d()

p_eso.post

p_eso.delta
```

# Tolerance Limits to Assess Agreement

The `tolerance_limit` function can be used to calculate the "tolerance limits". Typically a 95% prediction interval is calculated which provides the predicted difference between two measuring systems for 95% of future measurements pairs. However, we have to account for sampling error so we also need to estimate the confidence in the prediction intervals, and therefore we calculate the tolerance limits. So, when we have 95% tolerance limits for a 95% prediction interval, we can conclude that there is only a 5% probability (1-tolerance) that our tolerance limits do not contain the *true* prediction interval/limit.

## Rectal Temperature

So we will calculate the tolerance using the `loa_lme` function. We will need to identify the columns with the right information using the `diff`, `avg`, `condition`, and `id` arguments. We then select the right data set using the `data` argument. Lastly, we specify the specifics of the conditions for how the limits are calculated. For this specific analysis I decided to calculate 95% prediction intervals with 95% tolerance limits, and I will use percentile bootstrap confidence intervals.

```{r recpost}
# note: more accurate tolerance limits are given by tol_method = "perc"
rec.post_tol = tolerance_limit(
  data = df_rec.post,
  x = "PM",
  y = "AM",
  id = "id",
  condition = "trial_condition"
)


```

When we print a table of the tolerance limits, at least for Trec post exercise, are providing the same conclusion (poor agreement).

```{r}
print(rec.post_tol)
```

Furthermore, we can visualize the results with a Bland-Altman style plot of the tolerance.
Notice, despite the tighter cluster in the 3rd condition, the prediction/tolerance limits are wider.
This is a curious result, but if we inspect the results further (`rec.post_tol$emmeans`) we can see the degrees of freedom for this condition are horribly low.

```{r, fig.cap="Tolerance Limits for Trec Post Exercise",fig.width=7,fig.height=5}
plot(rec.post_tol)
```

Now, when we look at the Delta values for Trec we find that there is much closer agreement (maybe even acceptable agreement) when we look at tolerance limits. However, we cannot say that the average difference would be less than 0.25 which may not be acceptable for some researchers.

```{r}

rec.delta_tol = tolerance_limit(
  x = "PM",
  y = "AM",
  condition = "trial_condition",
  id = "id",
  data = df_rec.delta
)
  
rec.delta_tol
```

```{r, fig.cap="Tolerance Limits for Delta Trec",fig.width=7,fig.height=5}

# Plot Maximal Allowable Difference with delta argument
plot(rec.delta_tol,
     delta = .25)
```

## Esophageal Temperature

We can repeat the process for esophageal temperature. Overall, the results are fairly similar, and while there is better agreement on the delta (change scores), it is still fairly difficult to determine that there is "good" agreement between the AM and PM measurements.

```{r}

eso.post_tol = tolerance_limit(
  x = "AM",
  y = "PM",
  condition = "trial_condition",
  id = "id",
  data = df_eso.post
)

eso.delta_tol = tolerance_limit(
  x = "AM",
  y = "PM",
  condition = "trial_condition",
  id = "id",
  data = df_eso.delta
)
```

```{r}
eso.post_tol
```


```{r, fig.cap="Limits of Agreement for Teso Post Exercise",fig.width=7,fig.height=5}
plot(eso.post_tol)
```


```{r}
eso.delta_tol
```

```{r, fig.cap="Limits of Agreement for Delta Teso",fig.width=7,fig.height=5}
plot(eso.delta_tol,
     delta = .25)
```


# References