produce_tables_urine.Rmd

---
title: "Affinity-Assayed Urine Biomarker Stability Analysis"
author: "Alex Spiers"
date: "09/12/2021"
output: html_document
---

## Test for stability: repeated measures ANOVA

Conduct repeated measures ANOVA and extract F-statistic and related p-value to test null hypothesis of stability:

```{r process_data, include=FALSE}
library(tidyverse)
library(brms)
library(knitr)
library(jtools)
library(kableExtra)
library(lme4)
source("./process_affinity_data.R")
```


### Initial plots of raw data

Plots of raw salivary steroid concentration against time

```{r plot_raw, echo=TRUE}
ggplot(data = affinity_urine,
       aes(x=days, y=value, col=sample_id, group=sample_id)) +
  geom_line() +
  theme_minimal() +
  theme(legend.position = "none") +
  facet_wrap(vars(biomarker), scales = "free")
```

### Initial plots of transformed (scaled as proportion of baseline) data

A simple spline model is fitted to estimate general trend

```{r plot_proportions, echo=TRUE}
ggplot(data = affinity_urine,
       aes(x = days, y = 100 * (fraction - 1))) +
  ylab("% Change") +
  geom_line(aes(col=sample_id, group=sample_id)) +
  theme_minimal() +
  theme(legend.position = "none") +
  facet_wrap(vars(biomarker), scales = "free") +
  stat_smooth(method = "gam", formula = y ~ s(x, k = 4), size = 1)
```

```{r ANOVA, results='asis', echo=FALSE, warning=FALSE, message=FALSE}
# REMOVE OUTLIER
affinity_urine <- affinity_urine %>% filter(fraction < 5)

anova_results <- data.frame()
biomarker_name <- vector()

for (bio in unique(affinity_urine$biomarker)){
    df <- affinity_urine %>% filter(biomarker == bio) %>%
      mutate(time = years)
    
    lmer_fit <- lme4::lmer(formula = value ~ time + (1 | sample_id), data = df)
    small_model <- lme4::lmer(formula = value ~  (1 | sample_id), data = df)
    biomarker_name <- c(biomarker_name, bio)
    # Small Sample Inference for Fixed Effects from Restricted Maximum Likelihood - Kenwood Roger
    kr_test <- pbkrtest::KRmodcomp(lmer_fit, small_model) 
    anova_results <- rbind(anova_results, kr_test$test[1,])
    # print(plot(lmer_fit)) # Assess residuals for model fit
}

data.frame(anova_results) %>%
    as.data.frame(row.names = 1:nrow(.)) %>%
    mutate(biomarker = biomarker_name) %>%
    select(biomarker, everything()) %>%
    kbl() %>%
    kable_styling()

```

### Post hoc estimation of rate constant in exponential decay model


```{r posthoc_degrad, results='asis', echo=FALSE, warning=FALSE, message=FALSE}

df <- affinity_urine %>% filter(biomarker == "bio")

# linear_hier <- readRDS(paste0("./models/", bio, "_linear_mod2.rds"))
sfo_hier <- readRDS(paste0("./models/", "FSH_sfo_model2.rds"))

# linear_est <- fixef(linear_hier)[2,1]
# linear_se <- fixef(linear_hier)[2,2]
# linear_lower <- fixef(linear_hier)[2,3]
# linear_upper<- fixef(linear_hier)[2,4]

sfo_est <- fixef(sfo_hier)[2,1]
sfo_se <- fixef(sfo_hier)[2,2]
sfo_lower <- fixef(sfo_hier)[2,3]
sfo_upper<- fixef(sfo_hier)[2,4]

degrad_results <- data.frame(
                    # linear_est = linear_est,
                    # linear_se = linear_se,
                    # linear_lower = linear_lower,
                    # linear_upper = linear_upper,
                    sfo_est = sfo_est,
                    sfo_se = sfo_se,
                    sfo_lower = sfo_lower,
                    sfo_upper = sfo_upper
  )

degrad_results <- degrad_results %>%
  mutate(biomarker = "FSH") %>%
  select(biomarker, everything()) %>%
  mutate(across(sfo_est:sfo_upper, ~ signif(.x, 3)))


cat("\n")
cat("### Rate constant (years^-1^) \n") 
cat("\n")

#Publish to html
degrad_results %>%
  kbl() %>%
  kable_styling()

cat("\n")
cat("### Corresponding half-lives (years) \n") 
cat("\n")

half_lives <- data.frame(
  biomarker = "FSH",
  # linear_hl = -0.5 / degrad_results$linear_est,
  # linear_hl_lower = -0.5 / degrad_results$linear_lower,
  # linear_hl_upper = -0.5 / degrad_results$linear_upper,
  sfo_hl = log(
    2) / degrad_results$sfo_est,
    sfo_hl_lower = log(2) / degrad_results$sfo_upper,
    sfo_hl_upper = log(2) / degrad_results$sfo_lower) %>%
  mutate(across(sfo_hl:sfo_hl_upper, ~ signif(.x, 3)))

half_lives %>%
  kbl() %>%
  kable_styling()

cat("\n")
cat("### Estimated annual degradation as percentage \n") 
cat("\n")

annual_degradation <- data.frame(
  biomarker = "FSH",
  # linear_hl = -0.5 / degrad_results$linear_est,
  # linear_hl_lower = -0.5 / degrad_results$linear_lower,
  # linear_hl_upper = -0.5 / degrad_results$linear_upper,
  sfo_annual = (
    1 - exp(-degrad_results$sfo_est)) * 100,
    sfo_annual_lower = (1 - exp(-degrad_results$sfo_lower)) * 100,
    sfo_annual_upper = (1 - exp(-degrad_results$sfo_upper)) * 100
    ) %>%
  mutate(across(sfo_annual:sfo_annual_upper, ~ signif(.x, 3)))

annual_degradation %>%
  kbl() %>%
  kable_styling()
```


### Post hoc estimation of rate constant - linear accumulation in ELISA-assayed steroids

For Creatinine only

```{r posthoc_accum, results='asis', echo=FALSE, warning=FALSE, message=FALSE}

df <- affinity_urine %>% filter(biomarker == bio)

linear_hier <- readRDS(paste0("./models/", "Creatinine_linear_increase.rds"))

linear_est <- fixef(linear_hier)[2, 1]
linear_se <- fixef(linear_hier)[2, 2]
linear_lower <- fixef(linear_hier)[2, 3]
linear_upper <- fixef(linear_hier)[2, 4]

accum_results <- data.frame(
    linear_est = linear_est * 100,
    # linear_se = linear_se,
    linear_lower = linear_lower * 100,
    linear_upper = linear_upper * 100
)

accum_results <- accum_results %>%
    mutate(biomarker = "Creatinine") %>%
    select(biomarker, everything()) %>%
    mutate(across(linear_est:linear_upper, ~ signif(.x, 3)))

cat("\n")
cat("### Annual accumulation \n") 
cat("\n")

accum_results %>%
    kbl() %>%
    kable_styling()
```