```{r setup, echo = FALSE, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(dplyr)
library(tidyr)
library(lme4)
library(ggplot2)
library(broom.mixed)
library(tibble)
library(stringr)
library(car)
library(multcomp)
```

# Introduction {.tabset}
This script is to demonstrate some of the code that happens in the background of the application, using the simulated dataset and a simple linear model.

## Data Formatting

### Upload data in wide format and convert to long format
```{r}
# Load in wide format dataframe
# Demo data located here: https://github.com/TIDAL-modelling/TIDAL/tree/main/data

dataWide <- read.csv("data/emot_reg_emot_simulated.csv")

# Gather sdq_t columns into long format (this is a questionnaire score at each time point)
dataScore <- dataWide %>%
                gather(score_cat_test_col, score, all_of(c("sdq_t1", "sdq_t2", "sdq_t3", "sdq_t4", "sdq_t5"))) 

# Gather age_t columns into long format (this is age at each time point)
dataLong <- dataWide %>%
              gather(time_point, age, all_of(c("age_t1", "age_t2", "age_t3", "age_t4", "age_t5")))

# Merge them into one dataframe
dataLong[,"score"] <- dataScore[,"score"]

# Have a look at the dataframe
head(dataLong)

```
## Data Exploration {.tabset}

### Choose variables for modelling
```{r}
# This is similar to the drop down selections on the Data Exploration page.
ID <- "subject"
traj <- "score"
age <- "age"
modelType <- "Linear"
timePoint <- "time_point"
```

```{r}
# Mean center age (this helps the model fit)
dataLong <- dataLong %>%
          mutate(age_original = as.numeric(!!sym(age)) ) %>%
          mutate(!!sym(age) := as.numeric(!!sym(age)) - mean( as.numeric(!!sym(age)), na.rm = T ))
```


### Descriptive Statistics
```{r}
# Table:
dataLong %>%
          group_by(across( !!timePoint )) %>%
          summarise(N = sum(!is.na( !!sym(traj) )),
                    mean = mean(!!sym(traj), na.rm = T),
                    SD = sd(!!sym(traj), na.rm = T),
                    median = median(!!sym(traj), na.rm = T),
                    IQR = IQR(!!sym(traj), na.rm = T)
          )

# Plot
df.plot <- dataLong %>%
          group_by(across( !!timePoint )) %>%
          summarise(Age = mean(age_original, na.rm = T),
                    Phenotype = mean(!!sym(traj), na.rm = T),
                    SD = sd(!!sym(traj), na.rm = T),
                    n = sum(!is.na( !!sym(traj) ))) %>%
          mutate(upper = Phenotype + ( qnorm(0.975)*SD/sqrt(n) ),
                 lower = Phenotype - ( qnorm(0.975)*SD/sqrt(n) ))

ggplot(df.plot,aes(x=Age, y=Phenotype)) +
  geom_point()+
  geom_line() +
  geom_errorbar(aes(ymin = lower, ymax = upper))+
  ylab(paste0("Score (", traj, ")")) +
  xlab("Age")
```

### Model Results

#### Run the model using lme4
```{r}
# Run the model using lmer
fit <- lmer(formula = "score ~ age + (1 + age | subject)",
                          REML=F ,
                          data = dataLong,
                          control=lmerControl(optimizer="bobyqa",
                                              optCtrl=list(maxfun=2e5)))
# Summary of the model:
summary(fit)
```

#### N observations and groups
```{r}
paste0("The number of observations (measurements) is ",
               format(summary(fit)$devcomp$dims[[1]], big.mark=",", scientific=FALSE),
               " and the number of groups (people) is ",
               format(summary(fit)$ngrps[[1]], big.mark=",", scientific=FALSE) ,
               ".")
```

#### Fixed effects
```{r}
cbind(
  tidy(fit, "fixed"),
  confint(fit, "beta_", method = "Wald")) %>%
  mutate(p.z = 2 * (1 - pnorm(abs(statistic)))) %>%
  mutate(p.z = ifelse(p.z <= 0, "p < 0.001", round(p.z, 3)))
```

#### Random effects
```{r}
randomDF <- as.data.frame(VarCorr(fit),
              order = "lower.tri") %>%
              mutate(across(where(is.numeric), round, 3))
colnames(randomDF) <- c("Level", "Variable1", "Variable2", "Variance/Covariance", "SD Variance/Covariance")
randomDF
```

#### Deviance
```{r}
# This is the model fit, you can compare this to similar models.
deviance(fit)
```

### Plot
```{r}
# Add the "prediction"/model col to dataframe (adjustedScore)
# The alternative is to use the "predict" function, but this didn't work when users wanted to add categorical covariates.
# We had to calculate it manually instead. Below is just an example of how you do this for a linear model with no covariates.

# Pull out the mean centered age column
ageVec <- dataLong %>% pull(!!age)

# Get the coefficients from the model
coef <- summary(fit)$coefficients

# Calculate the model mean/prediciton. This is age multiplied by the slope plus the intercept.
adjustedScore <- ageVec * coef[2,1] + coef[1,1]

# Add this as a new column to the dataframe
modelDataEdit <- dataLong %>%
                  mutate(pred = adjustedScore)

# Add 95% CIs for these estimates using glht
ageOrig <- dataLong %>% pull(age_original)
ageOrig <- ageOrig[!is.na(ageOrig)]

# Get a vector of the ages we want to calculate 95% CIs for (instead of doing all the ages in the column)
ageCalcs <- c(min(ageOrig), seq(round(min(ageOrig), 1), round(max(ageOrig), 1), 0.5), max(ageOrig) )

# Go through each age and calculate the 95% CI
score_glht <- lapply(ageCalcs, function(x){

    # Get the row names of the coefficients to put in the glht function (this is Intercept and age)
    rowNames <- rownames(summary(fit)$coefficients)
    
    # Round age to 3 decimal places
    ageInput <- round(x - mean(ageOrig), 3)

    # Use glht to get 95% CIs
    # Essentially we are doing the intercept + slope * age == 0 as it is a linear model
    res <- multcomp::glht(fit, linfct = c( paste0(rowNames[1], " + ", rowNames[2], "*", ageInput, " == 0")) )

    # Tidy up the dataframe
    res <- tidy(confint(res))
    rowname <- paste0("Score (", traj, ")")
    res <-  res %>%
              mutate(contrast = paste0(rowname, " (95% CIs)")) %>%
              column_to_rownames(var = "contrast") %>%
              mutate(across(where(is.numeric), round, 2))
          
    return(res)
})

# Put everything into one dataframe
estimate <- do.call(rbind, score_glht) %>%
                      mutate(age = ageCalcs)

# Plot the mean trajectory and it's 95% CIs and the underlying descriptive data:
ggplot() +
          geom_line(data = modelDataEdit, aes(x= age_original ,  y = pred), color = "#1D86C7", linewidth = 1.5, na.rm=T) +
          geom_ribbon(data = estimate, aes(x= age , ymin = conf.low, ymax = conf.high), fill = "#1D86C7", alpha = 0.2, na.rm = T) +
          geom_point(data = df.plot,aes(x=Age, y=Phenotype))+
          geom_line(data = df.plot,aes(x=Age, y=Phenotype)) +
          geom_errorbar(data = df.plot, aes(x=Age, y=Phenotype, ymin = lower, ymax = upper)) +
          ylab(paste0("Score (", traj, ")")) +
          xlab("Age")
```

### Scores at ages
```{r}
# Choose some ages to calculate (these are equivalent to the checkboxes on the app)
ageInput <- c(4, 7, 12)

# Go through each age and let's calculate the score at that age using glht (similar to getting 95% CIs as above)
score <- lapply(as.numeric(ageInput), function(x){

          # Get the row names that correspond to the coefficients, ie. Intercept, Age (slope)
          rowNames <- rownames(summary(fit)$coefficients)

          # Round age and mean center it so it's on the correct scale
          ageInputRound <- round(x - mean(ageOrig), 3)
        
          # Essentially we are doing the intercept + slope * age == 0 as it is a linear model
          res <- multcomp::glht(fit, linfct = c( paste0(rowNames[1], " + ", rowNames[2], "*", ageInputRound, " == 0")) )

          # Tidy up the data frame by reformatting it and labelling it sensibly
          res <- tidy(confint(res))

          rowname <- paste0("Score (", traj, ")")

          res <-  res %>%
            mutate(contrast = paste0(rowname, " (95% CIs)")) %>%
            column_to_rownames(var = "contrast")
})

# Tidy it up into one dataframe
estimate <- lapply(score, function(df) {
   df %>%
     dplyr::select(estimate)
 })  %>% do.call(cbind, .)

colnames(estimate) <- ageInput
estimate <- estimate %>%
          gather(age, score, 1:ncol(estimate)) %>%
          mutate(age = as.numeric(age))

# Put the CIs into separate dataframes ready for plotting
conf.low <- lapply(score, function(df) {df %>% dplyr::select(conf.low)}) %>% do.call(cbind, .)
colnames(conf.low) <- ageInput
conf.low  <- conf.low %>%
          gather(age, conf.low, 1:ncol(conf.low)) %>%
          mutate(age = as.numeric(age))

conf.high <- lapply(score, function(df) {df %>% dplyr::select(conf.high)}) %>% do.call(cbind, .)
colnames(conf.high) <- ageInput
conf.high  <- conf.high %>%
          gather(age, conf.high, 1:ncol(conf.high)) %>%
          mutate(age = as.numeric(age))%>%
          dplyr::select(-age)

conf <- cbind(conf.low, conf.high, "age")

ggplot() +
  geom_line(data = modelDataEdit, aes(x= age_original ,  y = pred  ) , na.rm=T) +
  theme(legend.text = element_text(color = "black")) +
  geom_errorbar(data = conf, aes(x = age, ymin = conf.low, ymax = conf.high), width = 0.5) +
  geom_point(data = estimate, aes(x = age, y = score), col = "#1D86C7", size = 5) +
  ylab(paste0("Score (", traj, ")")) +
  xlab("Age")

# Tidy up the estimate and 95% CIs to display nicely to the user
estimateCI <- lapply(score, function(df) {
  df %>%
    mutate(estimateCI = paste0(round(estimate,3) , " (", round(conf.low,3) , " - ", round(conf.high,3) , ")")) %>%
    dplyr::select(estimateCI)
})  %>% do.call(cbind, .)
colnames(estimateCI) <- ageInput
estimateCI        
```

### Area under the curve
```{r}
# Choose 2 ages to calculate AUC for (this is the slider on the app)
AUCages <- c(4,10)

# Calculate the AUC for the ages the user has chosen
coef <- summary(fit)$coefficients

# Pull out the age column from the dataframe
ageOrig <- modelDataEdit %>%
              pull(age_original)
ageOrig <- ageOrig[!is.na(ageOrig)]

# Minus mean age so the age is on the same scale as that used in the model
age1 <- AUCages[1] - mean(ageOrig)
age2 <- AUCages[2] - mean(ageOrig)

# Remove any characters that means the deltaMethod function wont work
rowNames <- rownames(coef) %>%
              str_remove_all("I|\\(|\\^|\\)")

# Use the delta method to calculate AUC 
# AUC is calculate by doing age2 (the larger age) * intercept + slope * age2^2/2 minus age1 (the younger age) * intercept + slope * age2^2/2. This is for a linear model.
AUC <- car::deltaMethod(fit, c( paste0("(((", age2, ")*(", rowNames[1], ")) + ((", rowNames[2], ")*(", age2, ")^2/2)) - (((", age1,")*(", rowNames[1], ")) + ((", rowNames[2], ")*(", age1, ")^2/2))") ) , parameterNames = rowNames)

# Format the results
AUC <- paste0( round(AUC$Estimate, 2), " (", round(AUC$`2.5 %`,2), " - ", round(AUC$`97.5 %`,2), ")")

# Plot the results
ggplot() +
          geom_ribbon(data = modelDataEdit,
                      aes(x = age_original, ymax = pred, ymin = 0),
                      alpha = 0.1, show.legend = FALSE, fill = "#1D86C7") +
          geom_line(data = modelDataEdit, aes(x= age_original ,  y = pred ) , na.rm=T)+
          coord_cartesian(xlim = c(AUCages[1], AUCages[2])) +
          scale_colour_discrete(na.translate = F) +
          theme(legend.text = element_text(color = "black")) +
          ylab(paste0("Score (", traj, ")")) +
          xlab("Age") +
          scale_x_continuous(breaks = seq(round(min(modelDataEdit$age_original, na.rm =T)), round(max(modelDataEdit$age_original, na.rm =T)), by = 1),
                             expand = c(0, 0)) +
          scale_y_continuous(expand = c(0, 0))

# Format the results to display to the user
df <- t(data.frame(paste0(AUCages[1], " - ", AUCages[2]),
                     AUC))

rowname <- paste0("AUC (", traj, ") (95% CIs)")
rownames(df) <- c("Age Range", rowname)
df
```
## Interaction Variable (factor with two levels) {.tabset}

### Choose variable
```{r}
# Demo shown for categorical variable
# Split the trajectory by the sex variable
condition <- "female"

# Get the index for the column we want to split the trajectory on
colSplit <- which(colnames(dataLong) %in% condition)

# Make sure this column is a factor and replace some characters that may cause a problem with an underscore
modelData <- dataLong %>%
            filter(!is.na(!!sym(condition))) %>%
            mutate(!!condition := factor( str_replace_all(.[[colSplit]], " |-|/", "_") ) ) 
  

```

### Model Results
```{r}
# Run the model using lmer
fit <- lmer(formula = "score ~ age + female + age*female + (1 + age | subject)",
                          REML=F ,
                          data = modelData,
                          control=lmerControl(optimizer="bobyqa",
                                              optCtrl=list(maxfun=2e5)))
# Summary of the model:
summary(fit)

# Fixed effects
cbind(
            tidy(fit, "fixed"),
            confint(fit, "beta_", method = "Wald")) %>%
            mutate(p.z = 2 * (1 - pnorm(abs(statistic)))) %>%
            mutate(p.z = ifelse(p.z < 0.001, "p < 0.001", round(p.z, 3) ))

# Random effects
randomDF <- as.data.frame(VarCorr(fit),
                        order = "lower.tri") %>%
            mutate(across(where(is.numeric), round, 3))
          colnames(randomDF) <- c("Level", "Variable1", "Variable2", "Variance/Covariance", "SD Variance/Covariance")
          randomDF

```


```{r}
# add the "predicted" column to this dataset (it's not really a prediction because its the same dataset, it just shows the model)

# Pull out the age column
ageVec <- modelData %>% pull(!!age)

# Get the coefficients of the model
coef <- summary(fit)$coefficients

# Get the predicted column (note we did this manually in case covariates were added). If covariates were added then essentially this predicted "zero" variable is when the covariates are set to zero.

zero <- ageVec * coef[2,1] + coef[1,1]

# Get the row index from the coefficients table that correspond to the condition we are splitting on
rowIndex <- which(str_detect(string = row.names(coef),
                             pattern = condition) &
                  str_detect(string = row.names(coef),
                             pattern = ":", negate = T))

# Same as above but for the interaction term
rowIndexInteract1 <- which(str_detect(string = row.names(coef),
                                              pattern = condition) &
                                     str_detect(string = row.names(coef),
                                                pattern = ":") &
                                     str_detect(string = row.names(coef),
                                                pattern = "\\^", negate = T))

# How many levels does out condition have, eg. 2 if options are only male or female in this example
 n <- length(unique(pull(modelData, !!sym(condition))))
 
 # Calculate the prediction for the other level, ie. when the covariate is not set to zero
 predCovs <- lapply(1:(n-1), function(i){

              coef[1,1] +
              (ageVec * coef[2,1]) +
              coef[rowIndex[i],1] +
              (ageVec * coef[rowIndexInteract1[i],1])
 })

# Name of level, in this case female is coded as 0 or 1
num <- str_subset(row.names(coef), condition) %>%
        sub(paste0(".*", condition), "", .) %>%
        unique()

# Name the list by the the level it corresponds to
names(predCovs) <- paste0(condition, "_", num)

# Add these as new columns to the dataset and make a pred column which gives the correct predicted value for the corresponding condition level, ie. the female predicted value if the participant is female.
modelDataEdit <- cbind(modelData, do.call(cbind, predCovs)) %>%
            mutate(zero = zero) %>%
            mutate(!!condition := as.factor(.[[colSplit]]) ) %>%
            mutate(pred =  eval(parse(text =
                                        paste0(paste0("ifelse(", condition, " == '", num, "', ", condition, "_",num,",", collapse = " "), "zero",
                                               paste0(rep(")", length(num)), collapse = ""), collapse = "")
            )))
          
          

```

### Plot

```{r}
# use glht to calculate 95% CIs

# Pull out the age column      
ageOrig <- modelDataEdit %>% pull(age_original)
ageOrig <- ageOrig[!is.na(ageOrig)]

# Determine which ages we want to calculate 95% CI for
# Don't need to do it for every participant as age when age is fairly similar, instead do it as a sequence of 0.5 for every unit in age.
# This speeds up the calculation time
ageCalcs <-  c(min(ageOrig), seq(round(min(ageOrig), 1), round(max(ageOrig), 1), 0.5), max(ageOrig) )

# Go through each age and calculate 95% CI
CI_glht  <- lapply(ageCalcs, function(x){

          # Get number of levels in condition
          n <- length(unique(pull(modelDataEdit, !!sym(condition))))

          # Get corresponding rows for condition in coefficients table
          rowIndex <- which(str_detect(string = row.names(summary(fit)$coefficients),
                                       pattern = condition) &
                              str_detect(string = row.names(summary(fit)$coefficients),
                                         pattern = ":", negate = T))

          rowNames <- rownames(summary(fit)$coefficients)

          ageInput <- round(x - mean(ageOrig), 3)

          equations <- sapply(1:(n-1), function(i){
                paste0(rowNames[1], " + \`", rowNames[2], "\`*", ageInput, " + \`", rowNames[rowIndex[i]] , "\` + \`", rowNames[2], ":", rowNames[rowIndex[i]],"\`*",ageInput , " == 0")
              })
          res <- multcomp::glht(fit, linfct = c( paste0(rowNames[1], " + \`", rowNames[2], "\`*", ageInput, " == 0"), equations) )
          
          # Tidy results dataframe and rename rows/columns
          res <- tidy(confint(res))

          rowname <- paste0("Score (", traj, ")")
          
          levelNames <- as.character(levels(as.factor(pull(modelDataEdit, !!sym(condition)))))

            res <-  res %>%
              mutate(condition = levelNames)
          
            })

ageOrig <- modelDataEdit %>% pull(age_original)
ageOrig <- ageOrig[!is.na(ageOrig)]
ageCalcs <-  c(min(ageOrig), seq(round(min(ageOrig), 1), round(max(ageOrig), 1), 0.5), max(ageOrig) )
levelNames <- as.character(levels(as.factor(pull(modelDataEdit, !!sym(condition)))))

estimate <- do.call(rbind, CI_glht) %>%
          mutate(age = rep(ageCalcs, each = length(levelNames)) ) 

```

```{r}
ggplot() +
              geom_line(data = modelDataEdit, aes(x= age_original ,  y = pred, color = !!sym(condition) ) , na.rm=T) +
              geom_ribbon(data = estimate ,
                          aes(x= age , ymin = conf.low, ymax = conf.high, fill = condition), alpha = 0.2, na.rm = T, show.legend = FALSE) +
              scale_color_manual(values = c("orange", "blue")) +
              scale_fill_manual(values = c("orange", "blue")) +
              geom_point(data = df.plot,aes(x=Age, y=Phenotype))+
              geom_line(data = df.plot,aes(x=Age, y=Phenotype)) +
              geom_errorbar(data = df.plot, aes(x=Age, y=Phenotype, ymin = lower, ymax = upper)) +
              theme(legend.text = element_text(color = "black"))+
              ylab(paste0("Score (", traj, ")")) +
              xlab("Age")
```


### Scores at ages
```{r}
# Select ages to calculate scores for (same as the checkboxes on the app)
ageInputScore <- c(4, 9, 14)

# use glht to calculate scores at ages

ageOrig <- modelDataEdit %>% pull(age_original)
ageOrig <- ageOrig[!is.na(ageOrig)]

score_glht <- lapply(as.numeric(ageInputScore), function(x){


          n <- length(unique(pull(modelDataEdit, !!sym(condition))))


          rowIndex <- which(str_detect(string = row.names(summary(fit)$coefficients),
                                       pattern = condition) &
                              str_detect(string = row.names(summary(fit)$coefficients),
                                         pattern = ":", negate = T))

          rowNames <- rownames(summary(fit)$coefficients)

          ageInput <- round(x - mean(ageOrig), 3)

          equations <- sapply(1:(n-1), function(i){
                paste0(rowNames[1], " + ", rowNames[2], "*", ageInput, " + ", rowNames[rowIndex[i]] , " + ", rowNames[2], ":", rowNames[rowIndex[i]],"*",ageInput , " == 0")
              })
            res <- multcomp::glht(fit, linfct = c( paste0(rowNames[1], " + ", rowNames[2], "*", ageInput, " == 0"), equations) )
 
          # Tidy results dataframe and rename rows/columns
          res <- tidy(confint(res))

          rowname <- paste0("Score (", traj, ")")

          levelNames <- as.character(levels(as.factor(pull(modelDataEdit, !!sym(condition)))))

          res <-  res %>%
            mutate(contrast = paste0(rowname, " [", condition, ", level = ", levelNames, " ] (95% CIs)")) %>%
            column_to_rownames(var = "contrast")
})

# Plot the scores at the given age        
estimate <- lapply(score_glht, function(df) {
            df %>%
              dplyr::select(estimate)
          })  %>% do.call(cbind, .)

colnames(estimate) <- ageInputScore

estimate <- estimate %>%
              gather(age, score, 1:ncol(estimate)) %>%
              mutate(age = as.numeric(age))

conf.low <- lapply(score_glht, function(df) {df %>% dplyr::select(conf.low)}) %>% do.call(cbind, .)
colnames(conf.low) <- ageInputScore
conf.low  <- conf.low %>%
              gather(age, conf.low, 1:ncol(conf.low)) %>%
              mutate(age = as.numeric(age))

conf.high <- lapply(score_glht, function(df) {df %>% dplyr::select(conf.high)}) %>% do.call(cbind, .)
colnames(conf.high) <- ageInputScore
conf.high  <- conf.high %>%
              gather(age, conf.high, 1:ncol(conf.high)) %>%
              mutate(age = as.numeric(age))%>%
              dplyr::select(-age)

conf <- cbind(conf.low, conf.high, "age")

# Plot
ggplot() +
   geom_line(data = modelDataEdit, aes(x= age_original ,  y = pred, color = !!sym(condition) ) , na.rm=T) +
   theme(legend.text = element_text(color = "black")) +
   geom_errorbar(data = conf, aes(x = age, ymin = conf.low, ymax = conf.high), width = 0.5) +
   geom_point(data = estimate, aes(x = age, y = score), col = "#1D86C7", size = 5) +
   ylab(paste0("Score (", traj, ")")) +
   xlab("Age")


# Table
estimateCI <- lapply(score_glht, function(df) {
            df %>%
              mutate(estimateCI = paste0(round(estimate,3) , " (", round(conf.low,3) , " - ", round(conf.high,3) , ")")) %>%
              dplyr::select(estimateCI)
          })  %>% do.call(cbind, .)
colnames(estimateCI) <- ageInputScore
estimateCI


# Difference in scores at ages (with 95% CIs)
levelsScores <- c(0,1)
differenceScores <- lapply(as.numeric(ageInputScore), function(x){
            ageInput <- round(x - mean(ageOrig), 3)

            coef <- summary(fit)$coefficients

            rowIndex <- which(str_detect(string = row.names(coef),
                                         pattern = condition) &
                                str_detect(string = row.names(coef),
                                           pattern = ":", negate = T))

            rowNames <- rownames(coef) %>%
              str_remove_all("I|\\(|\\^|\\)|\\:")

            levelNames <- paste0(condition,levelsScores) %>%
              str_remove_all("I|\\(|\\^|\\)|\\:")

            levelNames1 <- rowNames[str_detect(rowNames, paste0(levelNames, collapse = "|"))]

            res <- deltaMethod(fit, c(paste0(
                  "(", rowNames[1], " + ", rowNames[2], "*", ageInput, " + ", levelNames1[1], " + ", levelNames1[2], "*", ageInput, ") - (",rowNames[1], " + ", rowNames[2], "*", ageInput, ")"
                )), parameterNames = rowNames )

            dif <- paste0( round(res$Estimate, 3), " (", round(res$`2.5 %`,3), " - ", round(res$`97.5 %`,3), ")")
            res <- data.frame(age = dif)
return(res)
})   

difTab <- do.call(cbind, differenceScores)
colnames(difTab) <- ageInputScore
rownames(difTab) <- paste0("Difference between ",levelsScores[1]," and ", levelsScores[2]," (95% CI)")
difTab
                
```

### Area under the curve
```{r}
# Select ages to calculate AUC for (same as the slider bar on the app)
AUCages <- c(4,10)

coef <- summary(fit)$coefficients

ageOrig <- modelDataEdit %>%
          pull(age_original)
ageOrig <- ageOrig[!is.na(ageOrig)]
age1 <- AUCages[1] - mean(ageOrig)
age2 <- AUCages[2] - mean(ageOrig)

rowNames <- rownames(coef) %>%
          str_remove_all("I|\\(|\\^|\\)|\\:")

AUC <- deltaMethod(fit, c( paste0("(((", age2, ")*(", rowNames[1], ")) + ((", rowNames[2], ")*(", age2, ")^2/2)) - (((", age1,")*(", rowNames[1], ")) + ((", rowNames[2], ")*(", age1, ")^2/2))") ) , parameterNames = rowNames)

AUC <- paste0( round(AUC$Estimate, 2), " (", round(AUC$`2.5 %`,2), " - ", round(AUC$`97.5 %`,2), ")")
        
rowIndex <- which(str_detect(string = row.names(summary(fit)$coefficients),
                                   pattern = condition) &
                          str_detect(string = row.names(summary(fit)$coefficients),
                                     pattern = ":", negate = T))

rowIndexInteract1 <- which(str_detect(string = row.names(summary(fit)$coefficients),
                                            pattern = condition) &
                                   str_detect(string = row.names(summary(fit)$coefficients),
                                              pattern = ":") &
                                   str_detect(string = row.names(summary(fit)$coefficients),
                                              pattern = "\\^", negate = T))
      
n <- length(unique(pull(modelDataEdit, !!sym(condition))))

AUCCovs <- lapply(1:(n-1), function(i){
         
            AUC <- deltaMethod(fit,
                               c( paste0(
                                 "(((",age2,")*(",rowNames[1],")) + ((",rowNames[2],")*(",age2,")^2/2) + ((",rowNames[rowIndex[i]],")*(",age2,")) + ((",rowNames[rowIndexInteract1[i]],")*(",age2,")^2/2) ) -
                                   (((",age1,")*(",rowNames[1],")) + ((",rowNames[2],")*(",age1,")^2/2) + ((",rowNames[rowIndex[i]],")*(",age1,")) + ((",rowNames[rowIndexInteract1[i]],")*(",age1,")^2/2) )"
                               ) )
                               , parameterNames = rowNames )

            AUC <- paste0( round(AUC$Estimate, 2), " (", round(AUC$`2.5 %`,2), " - ", round(AUC$`97.5 %`,2), ")")
            AUC
        })
        
AUC_delta <- list(AUC = AUC, AUCCovs = AUCCovs)
       
# Plot
ggplot(data = modelDataEdit) +
            geom_ribbon(data = modelDataEdit,
                        aes(x = age_original, ymax = pred, ymin = 0, fill = !!sym(condition)),
                        alpha = 0.1, show.legend = FALSE) +
            coord_cartesian(xlim = c(AUCages[1], AUCages[2])) +
            geom_line(aes(x = age_original, y = pred, color = !!sym(condition)), na.rm = TRUE) +
            theme(legend.text = element_text(color = "black")) +
            ylab(paste0("Score (", traj, ")")) +
            xlab("Age") +
            scale_x_continuous(breaks = seq(round(min(modelDataEdit$age_original, na.rm =T)), round(max(modelDataEdit$age_original, na.rm =T)), by = 1),
                               expand = c(0, 0))+
            scale_y_continuous(expand = c(0, 0))


# Table
df <- t(
            cbind(
              data.frame(paste0(AUCages[1], " - ", AUCages[2]),
                         AUC_delta$AUC),
              do.call(cbind, AUC_delta$AUCCovs)
            )
)

rowname <- paste0("AUC (", traj, ")")
levelNames <- as.character(levels(as.factor(pull(modelDataEdit, condition))))
rownames(df) <- c("Age Range", paste0(rowname, " [", condition, ", level = ", levelNames, " ] (95% CIs)") )
df
```


```{r}
# -------
# Difference in AUC
levelsAUC <- c(0,1)

levelNames <- paste0(condition,levelsAUC) %>%
          str_remove_all("I|\\(|\\^|\\)|\\:")

coef <- summary(fit)$coefficients

ageOrig <- modelDataEdit %>%
          pull(age_original)
ageOrig <- ageOrig[!is.na(ageOrig)]
age1 <- AUCages[1] - mean(ageOrig)
age2 <- AUCages[2] - mean(ageOrig)

rowNames <- rownames(coef) %>%
          str_remove_all("I|\\(|\\^|\\)|\\:")

levelNames1 <- rowNames[str_detect(rowNames, paste0(levelNames, collapse = "|"))]

res <- deltaMethod(fit, c( paste0(" (( ((",age2,")*(",rowNames[1],")) + ((",rowNames[2],")*(",age2,")^2/2)   + ( (",levelNames1[1],")*(",age2,") ) + ((",levelNames1[2],")*(",age2,")^2/2) ) - ( ((",age1,")*(",rowNames[1],")) + ((",rowNames[2],")*(",age1,")^2/2)   + ( (",levelNames1[1],")*(",age1,") ) + ((",levelNames1[2],")*(",age1,")^2/2) )) - (( ((",age2,")*(",rowNames[1],")) + ((",rowNames[2],")*(",age2,")^2/2) ) - ( ((",age1,")*(",rowNames[1],")) + ((",rowNames[2],")*(",age1,")^2/2)     )) ") ) ,parameterNames = rowNames)
    
dif <- paste0( round(res$Estimate, 2), " 95% CI: (", round(res$`2.5 %`,2), " - ", round(res$`97.5 %`,2), ")")
statement <- paste0("The difference between the two factor levels (for the age ranges ", AUCages[1], " - ", AUCages[2] ,") is ", dif ,".")
statement

```

## Interaction Variable (factor with four levels) {.tabset}

### Choose variable
```{r}
# Demo shown for categorical variable
# Split the trajectory by the er_bin variable
condition <- "er_bin"

# Get the index for the column we want to split the trajectory on
colSplit <- which(colnames(dataLong) %in% condition)

# Make sure this column is a factor and replace some characters that may cause a problem with an underscore
modelData <- dataLong %>%
            filter(!is.na(!!sym(condition))) %>%
            mutate(!!condition := factor( str_replace_all(.[[colSplit]], " |-|/", "_") ) ) 
  

```

### Model Results
```{r}
# Run the model using lmer
fit <- lmer(formula = "score ~ age + er_bin + age*er_bin + (1 + age | subject)",
                          REML=F ,
                          data = modelData,
                          control=lmerControl(optimizer="bobyqa",
                                              optCtrl=list(maxfun=2e5)))
# Summary of the model:
summary(fit)

# Fixed effects
cbind(
            tidy(fit, "fixed"),
            confint(fit, "beta_", method = "Wald")) %>%
            mutate(p.z = 2 * (1 - pnorm(abs(statistic)))) %>%
            mutate(p.z = ifelse(p.z < 0.001, "p < 0.001", round(p.z, 3) ))

# Random effects
randomDF <- as.data.frame(VarCorr(fit),
                        order = "lower.tri") %>%
            mutate(across(where(is.numeric), round, 3))
colnames(randomDF) <- c("Level", "Variable1", "Variable2", "Variance/Covariance", "SD Variance/Covariance")
randomDF

```


```{r}
# add the "predicted" column to this dataset (it's not really a prediction because its the same dataset, it just shows the model)

# Pull out the age column
ageVec <- modelData %>% pull(!!age)

# Get the coefficients of the model
coef <- summary(fit)$coefficients

# Get the predicted column (note we did this manually in case covariates were added). If covariates were added then essentially this predicted "zero" variable is when the covariates are set to zero.

zero <- ageVec * coef[2,1] + coef[1,1]

# Get the row index from the coefficients table that correspond to the condition we are splitting on
rowIndex <- which(str_detect(string = row.names(coef),
                             pattern = condition) &
                  str_detect(string = row.names(coef),
                             pattern = ":", negate = T))

# Same as above but for the interaction term
rowIndexInteract1 <- which(str_detect(string = row.names(coef),
                                              pattern = condition) &
                                     str_detect(string = row.names(coef),
                                                pattern = ":") &
                                     str_detect(string = row.names(coef),
                                                pattern = "\\^", negate = T))

# How many levels does out condition have, eg. 2 if options are only male or female in this example
 n <- length(unique(pull(modelData, !!sym(condition))))
 
 # Calculate the prediction for the other level, ie. when the covariate is not set to zero
 predCovs <- lapply(1:(n-1), function(i){

              coef[1,1] +
              (ageVec * coef[2,1]) +
              coef[rowIndex[i],1] +
              (ageVec * coef[rowIndexInteract1[i],1])
 })

# Name of level, in this case female is coded as 0 or 1
num <- str_subset(row.names(coef), condition) %>%
        sub(paste0(".*", condition), "", .) %>%
        unique()

# Name the list by the the level it corresponds to
names(predCovs) <- paste0(condition, "_", num)

# Add these as new columns to the dataset and make a pred column which gives the correct predicted value for the corresponding condition level, ie. the female predicted value if the participant is female.
modelDataEdit <- cbind(modelData, do.call(cbind, predCovs)) %>%
            mutate(zero = zero) %>%
            mutate(!!condition := as.factor(.[[colSplit]]) ) %>%
            mutate(pred =  eval(parse(text =
                                        paste0(paste0("ifelse(", condition, " == '", num, "', ", condition, "_",num,",", collapse = " "), "zero",
                                               paste0(rep(")", length(num)), collapse = ""), collapse = "")
            )))
          
          

```

### Plot

```{r}
# use glht to calculate 95% CIs

# Pull out the age column      
ageOrig <- modelDataEdit %>% pull(age_original)
ageOrig <- ageOrig[!is.na(ageOrig)]

# Determine which ages we want to calculate 95% CI for
# Don't need to do it for every participant as age when age is fairly similar, instead do it as a sequence of 0.5 for every unit in age.
# This speeds up the calculation time
ageCalcs <-  c(min(ageOrig), seq(round(min(ageOrig), 1), round(max(ageOrig), 1), 0.5), max(ageOrig) )

# Go through each age and calculate 95% CI
CI_glht  <- lapply(ageCalcs, function(x){

          # Get number of levels in condition
          n <- length(unique(pull(modelDataEdit, !!sym(condition))))

          # Get corresponding rows for condition in coefficients table
          rowIndex <- which(str_detect(string = row.names(summary(fit)$coefficients),
                                       pattern = condition) &
                              str_detect(string = row.names(summary(fit)$coefficients),
                                         pattern = ":", negate = T))

          rowNames <- rownames(summary(fit)$coefficients)

          ageInput <- round(x - mean(ageOrig), 3)

          equations <- sapply(1:(n-1), function(i){
                paste0(rowNames[1], " + \`", rowNames[2], "\`*", ageInput, " + \`", rowNames[rowIndex[i]] , "\` + \`", rowNames[2], ":", rowNames[rowIndex[i]],"\`*",ageInput , " == 0")
              })
          res <- multcomp::glht(fit, linfct = c( paste0(rowNames[1], " + \`", rowNames[2], "\`*", ageInput, " == 0"), equations) )
          
          # Tidy results dataframe and rename rows/columns
          res <- tidy(confint(res))

          rowname <- paste0("Score (", traj, ")")
          
          levelNames <- as.character(levels(as.factor(pull(modelDataEdit, !!sym(condition)))))

            res <-  res %>%
              mutate(condition = levelNames)
          
            })

ageOrig <- modelDataEdit %>% pull(age_original)
ageOrig <- ageOrig[!is.na(ageOrig)]
ageCalcs <-  c(min(ageOrig), seq(round(min(ageOrig), 1), round(max(ageOrig), 1), 0.5), max(ageOrig) )
levelNames <- as.character(levels(as.factor(pull(modelDataEdit, !!sym(condition)))))

estimate <- do.call(rbind, CI_glht) %>%
          mutate(age = rep(ageCalcs, each = length(levelNames)) ) 

```

```{r}
ggplot() +
              geom_line(data = modelDataEdit, aes(x= age_original ,  y = pred, color = !!sym(condition) ) , na.rm=T) +
              geom_ribbon(data = estimate ,
                          aes(x= age , ymin = conf.low, ymax = conf.high, fill = condition), alpha = 0.2, na.rm = T, show.legend = FALSE) +
              scale_color_manual(values = c("#E69F00","#56B4E9","#009E73","#F5C710")) +
              scale_fill_manual(values = c("#E69F00","#56B4E9","#009E73","#F5C710")) +
              geom_point(data = df.plot,aes(x=Age, y=Phenotype))+
              geom_line(data = df.plot,aes(x=Age, y=Phenotype)) +
              geom_errorbar(data = df.plot, aes(x=Age, y=Phenotype, ymin = lower, ymax = upper)) +
              theme(legend.text = element_text(color = "black"))+
              ylab(paste0("Score (", traj, ")")) +
              xlab("Age")
```


### Scores at ages
```{r}
# Select ages to calculate scores for (same as the checkboxes on the app)
ageInputScore <- c(4, 9, 14)

# use glht to calculate scores at ages

ageOrig <- modelDataEdit %>% pull(age_original)
ageOrig <- ageOrig[!is.na(ageOrig)]

score_glht <- lapply(as.numeric(ageInputScore), function(x){


          n <- length(unique(pull(modelDataEdit, !!sym(condition))))


          rowIndex <- which(str_detect(string = row.names(summary(fit)$coefficients),
                                       pattern = condition) &
                              str_detect(string = row.names(summary(fit)$coefficients),
                                         pattern = ":", negate = T))

          rowNames <- rownames(summary(fit)$coefficients)

          ageInput <- round(x - mean(ageOrig), 3)

          equations <- sapply(1:(n-1), function(i){
                paste0(rowNames[1], " + ", rowNames[2], "*", ageInput, " + ", rowNames[rowIndex[i]] , " + ", rowNames[2], ":", rowNames[rowIndex[i]],"*",ageInput , " == 0")
              })
            res <- multcomp::glht(fit, linfct = c( paste0(rowNames[1], " + ", rowNames[2], "*", ageInput, " == 0"), equations) )
 
          # Tidy results dataframe and rename rows/columns
          res <- tidy(confint(res))

          rowname <- paste0("Score (", traj, ")")

          levelNames <- as.character(levels(as.factor(pull(modelDataEdit, !!sym(condition)))))

          res <-  res %>%
            mutate(contrast = paste0(rowname, " [", condition, ", level = ", levelNames, " ] (95% CIs)")) %>%
            column_to_rownames(var = "contrast")
})

# Plot the scores at the given age        
estimate <- lapply(score_glht, function(df) {
            df %>%
              dplyr::select(estimate)
          })  %>% do.call(cbind, .)

colnames(estimate) <- ageInputScore

estimate <- estimate %>%
              gather(age, score, 1:ncol(estimate)) %>%
              mutate(age = as.numeric(age))

conf.low <- lapply(score_glht, function(df) {df %>% dplyr::select(conf.low)}) %>% do.call(cbind, .)
colnames(conf.low) <- ageInputScore
conf.low  <- conf.low %>%
              gather(age, conf.low, 1:ncol(conf.low)) %>%
              mutate(age = as.numeric(age))

conf.high <- lapply(score_glht, function(df) {df %>% dplyr::select(conf.high)}) %>% do.call(cbind, .)
colnames(conf.high) <- ageInputScore
conf.high  <- conf.high %>%
              gather(age, conf.high, 1:ncol(conf.high)) %>%
              mutate(age = as.numeric(age))%>%
              dplyr::select(-age)

conf <- cbind(conf.low, conf.high, "age")

# Plot
ggplot() +
   geom_line(data = modelDataEdit, aes(x= age_original ,  y = pred, color = !!sym(condition) ) , na.rm=T) +
   theme(legend.text = element_text(color = "black")) +
   geom_errorbar(data = conf, aes(x = age, ymin = conf.low, ymax = conf.high), width = 0.5) +
   geom_point(data = estimate, aes(x = age, y = score), col = "#1D86C7", size = 5) +
   ylab(paste0("Score (", traj, ")")) +
   xlab("Age")


# Table
estimateCI <- lapply(score_glht, function(df) {
            df %>%
              mutate(estimateCI = paste0(round(estimate,3) , " (", round(conf.low,3) , " - ", round(conf.high,3) , ")")) %>%
              dplyr::select(estimateCI)
          })  %>% do.call(cbind, .)
colnames(estimateCI) <- ageInputScore
estimateCI


# Difference in scores at ages (with 95% CIs)
# "Select two levels from your factor to calculate the difference in scores:"
levelsScores <- c(0,3)
differenceScores <- lapply(as.numeric(ageInputScore), function(x){
            ageInput <- round(x - mean(ageOrig), 3)

            coef <- summary(fit)$coefficients

            rowIndex <- which(str_detect(string = row.names(coef),
                                         pattern = condition) &
                                str_detect(string = row.names(coef),
                                           pattern = ":", negate = T))

            rowNames <- rownames(coef) %>%
              str_remove_all("I|\\(|\\^|\\)|\\:")

            levelNames <- paste0(condition,levelsScores) %>%
              str_remove_all("I|\\(|\\^|\\)|\\:")

            levelNames1 <- rowNames[str_detect(rowNames, paste0(levelNames, collapse = "|"))]

            res <- deltaMethod(fit, c(paste0(
                  "(", rowNames[1], " + ", rowNames[2], "*", ageInput, " + ", levelNames1[1], " + ", levelNames1[2], "*", ageInput, ") - (",rowNames[1], " + ", rowNames[2], "*", ageInput, ")"
                )), parameterNames = rowNames )

            dif <- paste0( round(res$Estimate, 3), " (", round(res$`2.5 %`,3), " - ", round(res$`97.5 %`,3), ")")
            res <- data.frame(age = dif)
return(res)
})   

difTab <- do.call(cbind, differenceScores)
colnames(difTab) <- ageInputScore
rownames(difTab) <- paste0("Difference between ",levelsScores[1]," and ", levelsScores[2]," (95% CI)")
difTab
                
```

### Area under the curve
```{r}
# Select ages to calculate AUC for (same as the slider bar on the app)
AUCages <- c(4,10)

coef <- summary(fit)$coefficients

ageOrig <- modelDataEdit %>%
          pull(age_original)
ageOrig <- ageOrig[!is.na(ageOrig)]
age1 <- AUCages[1] - mean(ageOrig)
age2 <- AUCages[2] - mean(ageOrig)

rowNames <- rownames(coef) %>%
          str_remove_all("I|\\(|\\^|\\)|\\:")

AUC <- deltaMethod(fit, c( paste0("(((", age2, ")*(", rowNames[1], ")) + ((", rowNames[2], ")*(", age2, ")^2/2)) - (((", age1,")*(", rowNames[1], ")) + ((", rowNames[2], ")*(", age1, ")^2/2))") ) , parameterNames = rowNames)

AUC <- paste0( round(AUC$Estimate, 2), " (", round(AUC$`2.5 %`,2), " - ", round(AUC$`97.5 %`,2), ")")
        
rowIndex <- which(str_detect(string = row.names(summary(fit)$coefficients),
                                   pattern = condition) &
                          str_detect(string = row.names(summary(fit)$coefficients),
                                     pattern = ":", negate = T))

rowIndexInteract1 <- which(str_detect(string = row.names(summary(fit)$coefficients),
                                            pattern = condition) &
                                   str_detect(string = row.names(summary(fit)$coefficients),
                                              pattern = ":") &
                                   str_detect(string = row.names(summary(fit)$coefficients),
                                              pattern = "\\^", negate = T))
      
n <- length(unique(pull(modelDataEdit, !!sym(condition))))

AUCCovs <- lapply(1:(n-1), function(i){
         
            AUC <- deltaMethod(fit,
                               c( paste0(
                                 "(((",age2,")*(",rowNames[1],")) + ((",rowNames[2],")*(",age2,")^2/2) + ((",rowNames[rowIndex[i]],")*(",age2,")) + ((",rowNames[rowIndexInteract1[i]],")*(",age2,")^2/2) ) -
                                   (((",age1,")*(",rowNames[1],")) + ((",rowNames[2],")*(",age1,")^2/2) + ((",rowNames[rowIndex[i]],")*(",age1,")) + ((",rowNames[rowIndexInteract1[i]],")*(",age1,")^2/2) )"
                               ) )
                               , parameterNames = rowNames )

            AUC <- paste0( round(AUC$Estimate, 2), " (", round(AUC$`2.5 %`,2), " - ", round(AUC$`97.5 %`,2), ")")
            AUC
        })
        
AUC_delta <- list(AUC = AUC, AUCCovs = AUCCovs)
       
# Plot
ggplot(data = modelDataEdit) +
            geom_ribbon(data = modelDataEdit,
                        aes(x = age_original, ymax = pred, ymin = 0, fill = !!sym(condition)),
                        alpha = 0.1, show.legend = FALSE) +
            coord_cartesian(xlim = c(AUCages[1], AUCages[2])) +
            geom_line(aes(x = age_original, y = pred, color = !!sym(condition)), na.rm = TRUE) +
            theme(legend.text = element_text(color = "black")) +
            ylab(paste0("Score (", traj, ")")) +
            xlab("Age") +
            scale_x_continuous(breaks = seq(round(min(modelDataEdit$age_original, na.rm =T)), round(max(modelDataEdit$age_original, na.rm =T)), by = 1),
                               expand = c(0, 0))+
            scale_y_continuous(expand = c(0, 0))


# Table
df <- t(
            cbind(
              data.frame(paste0(AUCages[1], " - ", AUCages[2]),
                         AUC_delta$AUC),
              do.call(cbind, AUC_delta$AUCCovs)
            )
)

rowname <- paste0("AUC (", traj, ")")
levelNames <- as.character(levels(as.factor(pull(modelDataEdit, condition))))
rownames(df) <- c("Age Range", paste0(rowname, " [", condition, ", level = ", levelNames, " ] (95% CIs)") )
df
```


```{r}
# -------
# Difference in AUC
# "Select two levels from your factor to calculate the difference in AUC between:"
levelsAUC <- c(0,3)

levelNames <- paste0(condition,levelsAUC) %>%
          str_remove_all("I|\\(|\\^|\\)|\\:")

coef <- summary(fit)$coefficients

ageOrig <- modelDataEdit %>%
          pull(age_original)
ageOrig <- ageOrig[!is.na(ageOrig)]
age1 <- AUCages[1] - mean(ageOrig)
age2 <- AUCages[2] - mean(ageOrig)

rowNames <- rownames(coef) %>%
          str_remove_all("I|\\(|\\^|\\)|\\:")

levelNames1 <- rowNames[str_detect(rowNames, paste0(levelNames, collapse = "|"))]

res <- deltaMethod(fit, c( paste0(" (( ((",age2,")*(",rowNames[1],")) + ((",rowNames[2],")*(",age2,")^2/2)   + ( (",levelNames1[1],")*(",age2,") ) + ((",levelNames1[2],")*(",age2,")^2/2) ) - ( ((",age1,")*(",rowNames[1],")) + ((",rowNames[2],")*(",age1,")^2/2)   + ( (",levelNames1[1],")*(",age1,") ) + ((",levelNames1[2],")*(",age1,")^2/2) )) - (( ((",age2,")*(",rowNames[1],")) + ((",rowNames[2],")*(",age2,")^2/2) ) - ( ((",age1,")*(",rowNames[1],")) + ((",rowNames[2],")*(",age1,")^2/2)     )) ") ) ,parameterNames = rowNames)
    
dif <- paste0( round(res$Estimate, 2), " 95% CI: (", round(res$`2.5 %`,2), " - ", round(res$`97.5 %`,2), ")")
statement <- paste0("The difference between the two factor levels (for the age ranges ", AUCages[1], " - ", AUCages[2] ,") is ", dif ,".")
statement

```



## {-}

# {-}

