run_affinity_urine_stability.Rmd

---
title: "Affinity Urine Stability Analysis"
author: "Alex Spiers"
date: "12/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")
if (getwd() == "C:/Users/as819/OneDrive - Imperial College London/PhD/SCAMP/BASS/Stability Analysis/stability-analysis"){
  cmdstanr::set_cmdstan_path("C:/Users/as819/Documents/cmdstan-2.28.1")
}
```


### 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
    print(paste("summary for", bio))
    print(summ(lmer_fit))
    print(confint(lmer_fit))
}

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

```

## Stability of Creatinine
## Fitting mixed effects linear model 
as Creatinine showed moderate increase

```{r creatinine_fit, echo=FALSE, debug=TRUE}

REFIT_SAVED_MODEL <- FALSE

b_class_prior <- paste0("normal(0, 0.5)")
intercept_prior <- paste0("normal(1, 0.02)")
sigma_prior <- paste0("student_t(4, 0, 0.2)")
group_variance_prior <- paste0("student_t(4, 0, 0.1)")
corr_matrix_prior <- "lkj_corr_cholesky(2)"

prior_linear_mod <- c(
  prior_string(b_class_prior, class = "b"),
  prior_string(intercept_prior, class = "Intercept"),
  prior_string(group_variance_prior, class = "sd"),
  prior_string(sigma_prior, class = "sigma"),
  prior_string(corr_matrix_prior, class = "L")
)

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


if (!file.exists("./models/Creatinine_linear_increase.rds") | REFIT_SAVED_MODEL){
  linear_mod <- brm(
      formula = fraction ~ years + (years | sample_id),
      data = df,
      prior = prior_linear_mod,
      #sample_prior = TRUE,
      family = gaussian(),
      chains = 6,
      cores = 12,
      iter = 3000,
      warmup = 1000,
      backend = "cmdstanr"#,
      #control=list(adapt_delta=0.99, max_treedepth=15)
  )

  saveRDS(linear_mod, "./models/Creatinine_linear_increase.rds")
}

linear_mod1 <- readRDS("./models/Creatinine_linear_increase.rds")

summary(linear_mod1)

```


## Stability of FSH
## Fitting mixed effects linear model 
as Creatinine showed moderate increase


```{r fitting_models, echo=FALSE, debug=TRUE}

REFIT_SAVED_MODELS = FALSE

SFO_k_fixed_formula <- brmsformula(
  fraction ~ C0 * exp(-k * years),
  C0~1+(1|sample_id), 
  k~1,
  nl=TRUE)

SFO_k_random_formula <- brmsformula(
  fraction ~ C0 * exp(-k * years),
  C0~1+(1|sample_id),
  k~1+(1|sample_id),
  nl=TRUE)


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

#SET PRIORS FOR linear models
location_b <- mean(df$fraction) %>% signif(., 2)
scale_b <- sd(df$fraction)

b_class_prior <- paste0("normal(0, 1)")
intercept_prior <- paste0("normal(1, 0.02)")
sigma_prior <- paste0("student_t(4, 0, 0.2)")
group_variance_prior <- paste0("student_t(4, 0, 0.1)")
corr_matrix_prior <- "lkj_corr_cholesky(2)"

prior_linear_mod1 <- c(
  prior_string(b_class_prior, class="b"),
  prior_string(intercept_prior, class="Intercept"),
  prior_string(group_variance_prior, class="sd"),
  prior_string(sigma_prior, class="sigma")
  )


if (!file.exists(paste0("./models/", "FSH_linear_mod1.rds")) | REFIT_SAVED_MODELS){
    print("Fitting linear model (fixed k)")
    linear_mod1 <- brm(
        formula = fraction ~ years + (1|sample_id),
        data = df,
        prior = prior_linear_mod1,
        #sample_prior = TRUE,
        family = gaussian(),
        chains = 6,
        cores = 12,
        iter = 3000,
        warmup = 1000,
        backend = "cmdstanr",
        #save_pars = save_pars(all = TRUE),
        control=list(adapt_delta=0.99, max_treedepth=15)
        )
    saveRDS(
        linear_mod1, 
        paste0("./models/", "FSH_linear_mod1.rds")
        )
    Sys.sleep(5)
} else {
    print("Loading from file linear model (fixed k)")
    linear_mod1 <- readRDS(paste0("./models/", "FSH_linear_mod1.rds"))
}

print(summary(linear_mod1))

prior_linear_mod2 <- c(
    prior_string(b_class_prior, class = "b"),
    prior_string(intercept_prior, class = "Intercept"),
    prior_string(group_variance_prior, class = "sd"),
    prior_string(sigma_prior, class = "sigma"),
     prior_string(corr_matrix_prior, class = "L")
)

print("Fitting linear model (k varies between subject)")

if (!file.exists(paste0("./models/", "FSH_linear_mod2.rds")) | REFIT_SAVED_MODELS){
linear_mod2 <- brm(
    formula = fraction ~ years + (years | sample_id),
    data = df,
    prior = prior_linear_mod2,
    #sample_prior = TRUE,
    family = gaussian(),
    chains = 6,
    cores = 12,
    iter = 3000,
    warmup = 1000,
    backend = "cmdstanr",
    control=list(adapt_delta=0.99, max_treedepth=15)
)
    saveRDS(
        linear_mod2, 
        paste0("./models/", "FSH_linear_mod2.rds"))
Sys.sleep(5)
} else {
    print("Loading from file linear model (between-subject k varies)")
    linear_mod2 <- readRDS(paste0("./models/", "FSH_linear_mod2.rds"))   
}

print(summary(linear_mod2))
#Set priors for SFO models
prior_SFO1 <- c(
    prior_string(intercept_prior, class="b", nlpar="C0"),
    prior_string("normal(0, 1)", lb=0, nlpar="k"),
    prior_string(group_variance_prior, class="sd", nlpar="C0"),
    prior_string(sigma_prior, class="sigma")
    )

if (!file.exists(paste0("./models/", "FSH_SFO_model1.rds")) | REFIT_SAVED_MODELS){
    print("Fitting SFO model (fixed k)")
    SFO_model1 <- brm(
        formula = SFO_k_fixed_formula,
        #sample_prior = TRUE,
        data = df,
        family=gaussian(),
        prior = prior_SFO1,
        chains = 6,
        cores = 12,
        iter = 3000,
        warmup = 1000,
        backend = "cmdstanr",
        control=list(adapt_delta=0.99, max_treedepth=15)
        )
    saveRDS(
        SFO_model1,
        paste0("./models/", "FSH_SFO_model1.rds"))

    Sys.sleep(5)
} else {
    print("Loading from file SFO model (fixed k)")
    SFO_model1 <- readRDS(paste0("./models/", "FSH_SFO_model1.rds"))   
}

print(summary(SFO_model1))

prior_SFO2 <- c(
    prior_string(intercept_prior, class="b", nlpar="C0"),
    prior_string("normal(0, 1)", lb=0, nlpar="k"),
    prior_string(group_variance_prior, class="sd", nlpar="k"),
    prior_string(group_variance_prior, class="sd", nlpar="C0"),
    prior_string(sigma_prior, class="sigma")
)
if (!file.exists(paste0("./models/", "FSH_SFO_model2.rds")) | REFIT_SAVED_MODELS){
    print("Fitting SFO model (k varies between subject)")
    SFO_model2 <- brm(
        formula = SFO_k_random_formula,
        #sample_prior = TRUE,
        data = df,
        family=gaussian(),
        prior = prior_SFO2,
        chains = 6,
        cores = 12,
        iter = 3000,
        warmup = 1000,
        backend = "cmdstanr",
        #save_pars = save_pars(all = TRUE),
        control=list(adapt_delta=0.99, max_treedepth=15)
    )
    #SFO_model2 <- add_criterion(SFO_model2, "loo", moment_match=TRUE)
    saveRDS(
        SFO_model2,
        paste0("./models/", "FSH_SFO_model2.rds")
    )
    Sys.sleep(5)
} else {
    print("Loading from file SFO model (k varies between subject)")
    SFO_model2 <- readRDS(paste0("./models/", "FSH_SFO_model2.rds"))   
}

print(summary(SFO_model2))

# print("comparing the LOO for models for")
# print(bio)
# print(loo_compare(linear_mod1, linear_mod2, SFO_model1, SFO_model2))
# print("completed analysis for")
# print(bio)
```