Classification Models.Rmd

---
title: "Investigating whether T1D models can identify T2-misdiagnoses."
author: "Ethan de Villiers"
date: "`r Sys.Date()`"
output: html_document
---
# Introduction


## Init

``` {r Packages}
library(tidyverse) # streamlining R coding since 2016
library(rio) # Import/Output master package
library(psych) # Essential for descriptive stats
library(pROC) # ROC curves and such.
library (rsq) # Used for R-squared values
library(ggplot2) # Advanced Graphs
library(tableone) # Nice descriptives for tables
library(rms) # Val.prob mainly
```

``` {r ImportingData}
filepath = 'startight Aug 22 comma separated.csv' # This must be the relative or absolute filepath pointing to the Startright August 2022 dataset
my_data <- import(filepath)

#> View(my_data)
```

## Prepping Data

Notes on Data
- age_atdiag is the exact patient age in double form
- AgeatDiagnosis is the patient age in int form (used)
- In equations where AGE is used, AgeatDiag was used instead - need to fix!
- no current clean GAD, ZNT8 code
- Islet autoantibodies were considered positive if: 
  - GADA ≥11 units/mL, IA2A ≥7.5 units/mL and ZNT8A ≥65 units/mL in those aged up to 30 years and ≥10 units/mL in those aged ≥30 years (17, 18)
- GRS Scores:
  - Correct GRS version = variable: "GRS"
  - No current, calibrated/standardised GRS scores against a type-1 enriched cohort

``` {r Cleaning and Prepping Data}
# Creating Anibody columns ####
#New column GADA_Status will default to NA, and if negative = 0, positive = 1
my_data = my_data %>%
  mutate(GAD_status = ifelse(GAD == 'negative', 0, NA)) %>%
  mutate(GAD_status = ifelse(grepl("[0-9]",GAD) & as.numeric(GAD) >= 11, 1, GAD_status))

#New column IA2_status will default to NA, and if negative = 0, positive = 1
my_data = my_data %>%
  mutate(IA2_status = ifelse(IA2 == 'negative', 0, NA)) %>%
  mutate(IA2_status = ifelse(grepl("[0-9]",IA2) & as.numeric(IA2) >= 7.5, 1, IA2_status))

#Filling in BMI for missing values
my_data = my_data %>%
  mutate(BMI = ifelse(is.na(BMI), Weight / (Height^2), BMI))

#Creating column for how many antibodies a patient has, sum()
my_data$Number_Antibodies = rowSums(my_data[,c("GAD_status", "IA2_status")], na.rm=TRUE)

#Creating binary column for if patients received any antibody testing
my_data = my_data %>%
  mutate(antibodies_tested = ifelse(!is.na(GAD_status) | !is.na(IA2_status), 1, 0 ) )

#Creating Binary outcome if patients were at least 1 antibody positive
my_data = my_data %>%
  mutate(Antibody_Positive = ifelse(Number_Antibodies >= 1, 1, 0))

#Creating Binary outcome if patients were MULTI antibody positive
my_data = my_data %>%
  mutate(Multi_Antibody_Positive = ifelse(Number_Antibodies >= 2, 1, 0))

#> Calibrating GRS-Scores
#Creating Type-1 enriched cohort
type_1_control_cohort = my_data %>%
  filter(Insulin == "Yes" & Number_Antibodies >= 2 & Type_of_diabetes == "Type 1") %>%
  filter(SCORE != '' & !is.na(SCORE)) # n = 456

#Calibrating GRS-Scores (pnorm())
my_data$Calibrated_GRS_Centiles = pnorm(my_data$SCORE, mean = mean(type_1_control_cohort$SCORE), sd = sd(type_1_control_cohort$SCORE))

#Calculating Time To Insulin - Date_of_diagnosis until Date_continuous insulin
# - Slightly problematic as only 134 patients have data for Date_continuous_insulin
my_data = my_data %>%
  mutate(Time_to_insulin = difftime(
    as.Date(Date_continuous_insulin, "%d-%b-%y"),
    as.Date(DateofDiagnosis, "%d-%b-%y"),
    units = "weeks")
  ) %>%
  mutate(Time_to_insulin = as.numeric(Time_to_insulin)) %>%
  mutate(Time_to_insulin = ifelse(Time_to_insulin < 0, 0, Time_to_insulin))

# Calculating Duration of Diabetes - Date of diagnosis until final clinic date = YEARS
my_data = my_data %>%
  mutate(Duration_diabetes = ifelse(
    V3Date_Contacted != "", 
    as.numeric(
      difftime(
        as.Date(V3Date_Contacted, "%m/%d/%Y"),
        as.Date(DateofDiagnosis, "%d-%b-%y"),
        units = "weeks"
      )/52
    ),
    as.numeric(
      difftime(
        as.Date(V2Date_Contacted, "%d-%b-%y"),
        as.Date(DateofDiagnosis, "%d-%b-%y"),
        units = "weeks"
      )/52
    ))
  )

hist(Type2_cohort$Duration_diabetes, na.rm = TRUE)
hist(Type2_cohort$Duration_followup, na.rm = TRUE)

# Calculating Median follow-up time from study
my_data = my_data %>%
  mutate(Duration_followup = ifelse(
    V3Date_Contacted != "", 
    as.numeric(
      difftime(
        as.Date(V3Date_Contacted, "%m/%d/%Y"),
        as.Date(Date_Visit, "%d-%b-%y"),
        units = "weeks"
      )/52
    ),
    as.numeric(
      difftime(
        as.Date(V2Date_Contacted, "%d-%b-%y"),
        as.Date(Date_Visit, "%d-%b-%y"),
        units = "weeks"
      )/52
    ))
  )

# Calculate Median time from diagnosis to recruitment
my_data = my_data %>%
  mutate(Duration_until_recruitment =
           as.numeric(
             difftime(
                as.Date(Date_Visit, "%d-%b-%y"),
                as.Date(DateofDiagnosis, "%d-%b-%y"),
                units = "weeks"
           )/52
           )
         )
```

Models Section:
This section generates and applies models that were previously developed. Although there are 7 models applied to this study, only 3 were included for the final abstract submission and presentation for Diabetes UK. To View the code, investigation and models that were NOT included, please view the "Master" branch on Github.

The models were taken from a previous study (found in Supplemental Material):
Lynam A, McDonald T, Hill A, Dennis J, Oram R, Pearson E, et al. Development and validation of multivariable clinical diagnostic models to identify type 1 diabetes requiring rapid insulin therapy in adults aged 18–50 years [Internet]. BMJ Open. British Medical Journal Publishing Group; 2019 [cited 2023Apr12]. Available from: https://bmjopen.bmj.com/content/9/9/e031586 

The models investigated that were INCLUDED in the presentation and abtract are as follows:
- Clinical Features Model:
  - = 37.94 + (-5.09 * log(age at diagnosis)) + (-6.34 * log(BMI))

- Clinical Features + Anti Model:
  - = 33.49649577 + (-4.665598345 * Log(Age)) + (-5.81137397 * Log(BMI)) + (3.082366 * AntiStatus1‡) + (3.494462 * AntiStatus2‡) + (4.350717 * AntiStatus3‡)
  - ‡where‡: AntiStatus1 = GAD +ve only, AntiStatus2 = IA-2 +ve only, AntiStatus3 = Both GADA and IA-2 +ve

- Clinical Features + Anti + GRS Model:
  - = 21.57649882 + (-4.086215772 * log(Age)) + (-5.096252172 * log(BMI)) + (2.702010666 * AntiStatus1‡) + (3.063255174 * AntiStatus2‡) + (3.813850704 * AntiStatus3‡) + (30.11052 * GRS)
  - ‡where‡: AntiStatus1 = GAD +ve only, AntiStatus2 = IA-2 +ve only, AntiStatus3 = Both GADA and IA-2 +ve

The models investigated that were NOT INCLUDED in the presentation and abtract are as follows:
- Clinical Features + GAD Model
  - = 34.8057844720 + (-4.801441792 * log (Age at diagnosis)) + (-5.980577792 * log(BMI)) + (2.937107976 * GADA)

- Clinical Features + IA2 Model
  - = 37.26905033 + (3.194096 * IA2† ) + (-5.047657308 * Log(Age)) + (-6.287258808 * Log(BMI))

- Clinical Features + GRS Model
  - = 24.46138054 + (-4.443506884 * Log(Age)) + (-5.534741384 * Log(BMI)) + (33.93968* T1D GRS)

- Clinical Features + Lipids Model
  - = 9.0034272  - (0.1915482 * BMI) – (0.1686227 * age at diagnosis) + (0.3026012 if female) – (0.2269216 * cholesterol) + (1.540850 * HDL) – (0.2784059 * triglycerides)

``` {r Generating Model LogOdds}
my_data = my_data %>% 
  mutate(
    logOR_clinF = 37.94 +
      (-5.09*log(AgeatDiagnosis)) +
      (-6.34 * log(BMI)) 
  ) %>% # First Model: Clinical Features alone (Age + BMI)
  mutate( 
    logOR_clinF_Antibodies = 33.49649577 +
            (-4.665598345 * log(AgeatDiagnosis)) +
            (-5.81137397 * log(BMI)) +
            (3.082366 * GAD_status) +
            (3.494462 * IA2_status) +
            (4.350717 * ifelse(Number_Antibodies >= 2, 1, 0))
  ) %>% # Second Model: ClinF and Antibodies (Age, BMI, GAD, IA2)
  mutate(
    logOR_clinF_antiGRS = 21.57649882 +
            (-4.086215772 * log(AgeatDiagnosis)) +
            (-5.096252172 * log(BMI)) +
            (2.702010666 * GAD_status) +
            (3.063255174 * IA2_status) +
            (3.813850704 * ifelse(Number_Antibodies >= 2, 1, 0)) +
            (30.11052 * GRS)
  ) # Third Model: ClinF + Anti + GRS (Age, BMI, GAD, IA2, GRS)
```

``` {r Calculating Model Probabilities from logOdds Ratio}
#> To calculate model probability from the previous LogOdds, the following formula must be applied:
#>  - exp(LogOR) / (1 + exp(LogOR)
my_data = my_data %>%
  mutate(model_clinF = exp(logOR_clinF) / (1 + exp(logOR_clinF))) %>%
  mutate(model_clinF_Anti = exp(logOR_clinF_Antibodies) / (1 + exp(logOR_clinF_Antibodies))) %>%
  mutate(model_clinF_AntiGRS = exp(logOR_clinF_antiGRS) / (1 + exp(logOR_clinF_antiGRS)))
```

Non-Insulin Cohort Building

This following section is where the initially non-insulin treated cohort is created.
 - Patients who are Initially treated != (not as) Insulin
 - Removing patients who have Date continuous Insulin within 2 weeks of their Diagnosis
 
``` {r Non-Insulin Cohort Building}
#> Non-Insulin cohort referred to as "Type2_cohort" for remainder of script

#Type_2 = any patient initially NOT treated with Insulin
Type2_cohort = my_data %>%
  filter(Initial_diabetes_Insulin != "Insulin")

# Removing those who have a time_to_insulin of less than 2 weeks
Type2_cohort = Type2_cohort %>%
  filter(is.na(Time_to_insulin) | Time_to_insulin > 2)

# Filter patients with Visit 1 positive insulin, where visit 1 is >= 2 weeks from diagnosis
Type2_cohort = Type2_cohort %>%
  filter(as.numeric(
      difftime(
        as.Date(Date_Visit, "%d-%b-%y"),
        as.Date(DateofDiagnosis, "%d-%b-%y"),
        units = "weeks"
      )) > 2 | Insulin != "Yes"
  ) 

#> CREATING PROGRESSION OUTCOME:
#>  - Binary outcome if patients progressed to Insulin or Not
#>  - If Insulin == "Yes", patients have progressed to Insulin on Visit 1
#>  - If V3Insulin == "Yes", patients have progressed by the end of the study
#>  - If there is no available data for V3, Use V2Insulin
Type2_cohort = Type2_cohort %>%
  mutate(progressed = ifelse( # If on Insulin by study recruitment, progressed
    Insulin == "Yes",
    1,
    ifelse(
      V3Insulin == "Yes", # If on Insulin by Final Visit (V3) recruitment, progressed
      1,
      ifelse(
        V2Insulin == "Yes", # If no data on V3, if Insulin by V2, progressed
        1,
        0
      )
    )
  ))

# Looking at Cohort
# View(Type2_cohort)
```

``` {r Basic Descriptives for Cohort Creation Flowchart}
#> UNCOMMENT FOLLOWING LINE: See total number of progressed patients = 141
# nrow(Type2_cohort %>% filter(progressed == 1))
# nrow(Type2_cohort %>% filter(progressed == 0))

#> UNCOMMENT FOLLOWING LINES: See total number of progressed patients who had model outcomes = 
# nrow(Type2_cohort %>% filter(progressed == 1 & !is.na(model_clinF)))
# nrow(Type2_cohort %>% filter(progressed == 1 & !is.na(model_clinF_Anti)))
# nrow(Type2_cohort %>% filter(progressed == 1 & !is.na(model_clinF_AntiGRS)))
# nrow(Type2_cohort %>% filter(progressed == 1 & !is.na(model_clinF) & !is.na(model_clinF_Anti) & !is.na(model_clinF_AntiGRS)))

#> UNCOMMENT FOLLOWING LINES: See total number of progressed patients who had model outcomes = 
# nrow(Type2_cohort %>% filter(progressed == 0 & !is.na(model_clinF)))
# nrow(Type2_cohort %>% filter(progressed == 0 & !is.na(model_clinF_Anti)))
# nrow(Type2_cohort %>% filter(progressed == 0 & !is.na(model_clinF_AntiGRS)))
# nrow(Type2_cohort %>% filter(progressed == 0 & !is.na(model_clinF) & !is.na(model_clinF_Anti) & !is.na(model_clinF_AntiGRS)))

#> UNCOMMENT FOLLOWING LINES: Number of participants with model data
# nrow(Type2_cohort %>% filter(!is.na(model_clinF)))
# nrow(Type2_cohort %>% filter(!is.na(model_clinF_Anti)))
# nrow(Type2_cohort %>% filter(!is.na(model_clinF_AntiGRS)))
```

## Data Analysis

###Descriptive Statistics
``` {r Type-2 Cohort Baseline Table}
table = CreateTableOne(
    data=Type2_cohort,
    vars = c("BMI", "AgeatDiagnosis", "GRS", "HbA1c_at_diagnosis", "Duration_diabetes",
             "Duration_followup", "Duration_until_recruitment"),
    strata = "progressed"
  )

summary(table)

write.csv(table, file = "testTable.csv")
#> Take .CSV file and upload text to: 
#> https://www.becsv.com/csv-table.php
#> Generate and HTML table, and download to an .htm document
#> Open doc and copy/paste HTML table into word --> powerpoint :)
#> OR open .csv in Excel and copy paste across

#> Calculating How many were single/double antibody positive/negative
#> To use this section:
#>  - change progressed == either 1 or 0 to select different cohorts
#>  - change number_antibodies >= 2, == 1, or >=1 for different descriptives
nrow(
  Type2_cohort %>%
    filter(progressed == 0) %>%
    filter(Number_Antibodies > 1)
)

#> P-values for table
#> To use: Switch both Type2_cohort$[variable] for variable of choice (recommended vars outlined in above table)
ttest_normal = Type2_cohort$GRS[Type2_cohort$progressed==0]
ttest_insulin = Type2_cohort$GRS[Type2_cohort$progressed==1]
t.test(ttest_normal, y = ttest_insulin)
```

### Models Section

``` {r Creating ROC_DATA dataframe for ROC Cohort}
# Creating a new dataframe to work with cohort that has ALL Model outcomes
roc_data = Type2_cohort %>%
  filter(!is.na(progressed) &!is.na(model_clinF) & !is.na(model_clinF_Anti) &
           !is.na(model_clinF_AntiGRS) )
```

``` {r 1: Clinical Features Model}
#> This section dedicated to the Clin F slide
#> Incorporates: ROC and ggBoxPlot
#> To view/generate Calibration plot, please proceed to Chunk 12, titled: "Model Calibration Plots"

# ROC AUC
glm.fit = glm(roc_data$progressed ~ roc_data$model_clinF, family = binomial())
plot.new()
lines(roc_data$model_clinF, glm.fit$fitted.values)

(roc = roc(roc_data$progressed, glm.fit$fitted.values, auc=TRUE, plot=TRUE))

#> To View Threshold and specificity/sensitivity UNCOMMENT the following line:
# coords(roc, x='best')

#Distribution 2.0: Boxplot
ggplot(Type2_cohort %>%
         mutate(progressed = ifelse(progressed == 0,
                                    "Did not Progress",
                                    "Progressed")
                )
        ,aes(x=as.factor(progressed),y=model_clinF)
       ) + 
  geom_boxplot() +
  ylab("Clinical Features Model Probability") +
  xlab("Progression to Insulin <= 3 years")

```

``` {r 2: Clin F + GAD and IA2 Model}
#> This section dedicated to the Clin F + Antibodies slide
#> Incorporates: ROC and ggBoxPlot
#> To view/generate Calibration plot, please proceed to Chunk 12, titled: "Model Calibration Plots"

# ROC AUC
glm.fit = glm(roc_data$progressed ~ roc_data$model_clinF_Anti, family = binomial())
plot.new()
lines(roc_data$model_clinF_Anti, glm.fit$fitted.values)

(roc = roc(roc_data$progressed, glm.fit$fitted.values, auc=TRUE, plot=TRUE))

#> To View Threshold and specificity/sensitivity UNCOMMENT the following line:
# coords(roc, x='best')

#Distribution 2.0: Boxplot
ggplot(Type2_cohort %>%
         mutate(progressed = ifelse(progressed == 0,
                                    "Did not Progress",
                                    "Progressed")
                )
        , aes(as.factor(progressed), model_clinF_Anti)) +
  geom_boxplot() +
  xlab("Progression to Insulin <= 3 years") +
  ylab("Clinical Features + Antibodies Model Probability")
```

``` {r 3: Clin F + Anti + GRS Model}
#> This section dedicated to the Clin F + Antibodies + GRS slide
#> Incorporates: ROC and ggBoxPlot
#> To view/generate Calibration plot, please proceed to Chunk 12, titled: "Model Calibration Plots"

#ROC AUC
glm.fit = glm(roc_data$progressed ~ roc_data$model_clinF_AntiGRS, family = binomial())
plot.new()
lines(roc_data$model_clinF_AntiGRS, glm.fit$fitted.values)

(roc = roc(roc_data$progressed, glm.fit$fitted.values, auc=TRUE, plot=TRUE))

#> To View Threshold and specificity/sensitivity UNCOMMENT the following line:
# coords(roc, x='best')

#Distribution 2.0: Boxplot
ggplot(Type2_cohort %>%
         mutate(progressed = ifelse(progressed == 0,
                                    "Did not Progress",
                                    "Progressed")
                )
        , aes(as.factor(progressed), model_clinF_AntiGRS)) +
  geom_boxplot() +
  xlab("Progression to Insulin <= 3 years") +
  ylab("Clinical Features + Anti + GRS Model Probability")
```

``` {r Model Calibration Plots}
#> To Use following section:
#>  - Reassign model_probabilities variable to your chosen Model Probs
#>  - Change the labs( title = "") on line 430 to your title
model_probabilities = Type2_cohort$model_clinF_AntiGRS

quantiles <- quantile(model_probabilities, probs = seq(0, 1, by = 0.1), na.rm = TRUE)
quantiles[1] <- 0.99 * quantiles[1]
quantiles[length(quantiles)] <- 1.1 * quantiles[length(quantiles)]
processed_quantiles <- cut(
  model_probabilities,
  breaks = quantiles,
  include_lowest = TRUE
)
processed_quantiles <- data.frame(y = Type2_cohort$progressed, pred = model_probabilities, dec = processed_quantiles) %>%
  group_by(dec) %>%
  mutate(prob_obs = sum(y) / n(),
    obs = sum(y),
    n_group = n(),
    mnpred = mean(pred),
    lower = lapply(sum(y), prop.test, n = n()),
    upper = sapply(lower, function(x) x$conf.int[2]),
    lower = sapply(lower, function(x) x$conf.int[1]))
## plot
ggplot(processed_quantiles, aes(x = mnpred, y = prob_obs)) +
  geom_point() +
  xlab("Mean predicted probability in each decile") +
  ylab("Observed probability in each decile") +
  labs(title = "Clin F + Antibodies + GRS Model Cal. Plot") +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  ylim(c(0, 1)) + xlim(c(0, 1)) +
  geom_errorbar(aes(ymin = lower, ymax = upper))

```


``` {r Antibodies predicting Progression to Insulin}
#> This section dedicated to the Antibodies comparative slide
# ROC
glm.fit = glm(roc_data$progressed ~ roc_data$Number_Antibodies, family = binomial())
plot.new()
lines(roc_data$Number_Antibodies, glm.fit$fitted.values)

(roc = roc(roc_data$progressed, glm.fit$fitted.values, auc=TRUE, plot=TRUE))
# coords(roc, x='best')
```


## Clinical Application Section

!!!!!       TEMPORARILY REMOVED FROM PRESENTATION PENDING FURTHER DECISION        !!!!!
!!!!!       TEMPORARILY REMOVED FROM PRESENTATION PENDING FURTHER DECISION        !!!!!
!!!!!       TEMPORARILY REMOVED FROM PRESENTATION PENDING FURTHER DECISION        !!!!!

The following section is focused around illistrating how these models can be used in a clinical context, to aid diagnosis.
``` {r Models as Clinical Diagnostic Tools} 
#> This section is dedicated to the The clinical tools slide:
#>  - Identify 3 patients that are phenotypically similar, but have drastic different pathologies

#> True, Type-1 patient (Antibodies and Progression): SR0280
#> Type-2 patient: SR2462
#> False Type-1 patient (antibody, no progression): SR1707
Type2_cohort %>%
  filter(StartRightID == "SR0280" | StartRightID == "SR2462" | StartRightID == "SR1707") %>%
  select(StartRightID, BMI, AgeatDiagnosis, model_clinF, Number_Antibodies, model_clinF_Anti, progressed)
```

```{r Clinical Features Model as a First-Line Screening Tool}
#> This chunk is used to investigate the Clinical Features model as a screening tool

# ROC AUC for Clinical Features Model discerning Single-Antibody Positivity
glm.fit = glm(roc_data$Antibody_Positive ~ roc_data$model_clinF, family = binomial())
plot.new()
lines(roc_data$model_clinF, glm.fit$fitted.values)

(roc = roc(roc_data$Antibody_Positive, glm.fit$fitted.values, auc=TRUE, plot=TRUE))
# coords(roc, x='best')

# ROC AUC for Clinical Features Model discerning MULTI-Antibody Positivity
glm.fit = glm(roc_data$Multi_Antibody_Positive ~ roc_data$model_clinF, family = binomial())
plot.new()
lines(roc_data$model_clinF, glm.fit$fitted.values)

(roc = roc(roc_data$Multi_Antibody_Positive, glm.fit$fitted.values, auc=TRUE, plot=TRUE))
# coords(roc, x='best')

#> CALIBRATION PLOT for Clin F Model Calibrating against Antibody Positivity
#>  - Change "$Antibody_Positive" to "$Multi
quantiles <- quantile(Type2_cohort$model_clinF, probs = seq(0, 1, by = 0.1), na.rm = TRUE)
quantiles[1] <- 0.99 * quantiles[1]
quantiles[length(quantiles)] <- 1.1 * quantiles[length(quantiles)]
processed_quantiles <- cut(
  Type2_cohort$model_clinF,
  breaks = quantiles,
  include_lowest = TRUE
)
processed_quantiles <- data.frame(y = Type2_cohort$Multi_Antibody_Positive, pred = Type2_cohort$model_clinF, dec = processed_quantiles) %>%
  group_by(dec) %>%
  mutate(prob_obs = sum(y) / n(),
    obs = sum(y),
    n_group = n(),
    mnpred = mean(pred),
    lower = lapply(sum(y), prop.test, n = n()),
    upper = sapply(lower, function(x) x$conf.int[2]),
    lower = sapply(lower, function(x) x$conf.int[1]))
## plot
ggplot(processed_quantiles, aes(x = mnpred, y = prob_obs)) +
  geom_point() +
  xlab("Mean predicted probability in each decile") +
  ylab("Observed probability in each decile") +
  labs(title = "Clinical Features on Multi-Antibody Positivity Cal. Plot") +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  ylim(c(0, 1)) + xlim(c(0, 1)) +
  geom_errorbar(aes(ymin = lower, ymax = upper))
```

Clinical Features + Antibodies Model Adds to Antibodies Interpretation
 - Section use for final 2 results slides in presentation:
 - Does Antibodies Model perform well/is able to discriminate in those who HAVE BEEN antibody tested?
 - DO they offer clinical use after antibodies?

``` {r Antibodies Model in Antibody tested}
#> This chunk was used to create all 6 graphs on final 2 results slides: 2 ROC, 2 Calibration, 2 Box.
#> To use this section:
#>  - Set temp_roc_data filter to either be Antibody_Positive == 1 OR Antibody_Positive == 0 for +ves/-ves
#>  - Change model_probabilities to the specified model (either "model_clinF_Anti" OR "model_clinF_AntiGRS") 
#>  - Change model probabilities on line 569 to produce associated box plot

#Creating temporary data frame to avoid messing with other ROC data and curves
#This will be re-written later for second slide in this section
temp_roc_data = roc_data %>%
  filter(Antibody_Positive == 1)

#How many patients of this group progressed?
nrow(
  temp_roc_data %>%
  filter(progressed == 1)
)

# Reassign model_probabilities
model_probabilities = temp_roc_data$model_clinF_Anti

#> ROC AUC
glm.fit = glm(temp_roc_data$progressed ~ model_probabilities, family = binomial())
plot.new()
lines(model_probabilities, glm.fit$fitted.values)

(roc = roc(temp_roc_data$progressed, glm.fit$fitted.values, auc=TRUE, plot=TRUE))
# coords(roc, x='best')

#Distribution 2.0: Boxplot
ggplot(
  temp_roc_data %>%
    mutate(progressed = ifelse(progressed == 1, "Progressed", "Did not progress")),
  aes(progressed, model_clinF_AntiGRS)
) +
  geom_boxplot() +
  xlab("Progression to Insulin") +
  ylab("Clin F + Antibodies + GRS Model Probabilities") +
  ylim(c(0,1))

#> CALIBRATION PLOT
quantiles <- quantile(model_probabilities, probs = seq(0, 1, by = 0.1), na.rm = TRUE)
quantiles[1] <- 0.99 * quantiles[1]
quantiles[length(quantiles)] <- 1.1 * quantiles[length(quantiles)]
processed_quantiles <- cut(
  model_probabilities,
  breaks = quantiles,
  include_lowest = TRUE
)
processed_quantiles <- data.frame(y = temp_roc_data$progressed, pred = model_probabilities, dec = processed_quantiles) %>%
  group_by(dec) %>%
  mutate(prob_obs = sum(y) / n(),
    obs = sum(y),
    n_group = n(),
    mnpred = mean(pred),
    lower = lapply(sum(y), prop.test, n = n()),
    upper = sapply(lower, function(x) x$conf.int[2]),
    lower = sapply(lower, function(x) x$conf.int[1]))
## plot
ggplot(processed_quantiles, aes(x = mnpred, y = prob_obs)) +
  geom_point() +
  xlab("Mean predicted probability in each decile") +
  ylab("Observed probability in each decile") +
  labs(title = "Clin F + Anti + GRS Model Cal. Plot") +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  ylim(c(0, 1)) + xlim(c(0, 1)) +
  geom_errorbar(aes(ymin = lower, ymax = upper))

#Cleaning Environment
remove(temp_roc_data)

```

``` {r Antibodies Model in Older Agegroup}
# This section dedicated to investigating the antibodies model in those over 50yo
#Re-assigning temp-roc data
temp_roc_data = roc_data %>%
  filter(AgeatDiagnosis > 50)

#How many patients of this group progressed?
nrow(
  temp_roc_data %>%
  filter(progressed == 1)
)

# ROC AUC
glm.fit = glm(temp_roc_data$progressed ~ temp_roc_data$model_clinF_Anti, family = binomial())
plot.new()
lines(temp_roc_data$model_clinF_Anti, glm.fit$fitted.values)

(roc = roc(temp_roc_data$progressed, glm.fit$fitted.values, auc=TRUE, plot=TRUE))

#> To View Threshold and specificity/sensitivity UNCOMMENT the following line:
# coords(roc, x='best')

#Distribution 2.0: Boxplot
ggplot(temp_roc_data %>%
         mutate(progressed = ifelse(progressed == 0,
                                    "Did not Progress",
                                    "Progressed")
                )
        , aes(as.factor(progressed), model_clinF_Anti)) +
  geom_boxplot() +
  xlab("Progression to Insulin <= 3 years") +
  ylab("Clinical Features + Antibodies Model Probability")

#Calibration
quantiles <- quantile(temp_roc_data$model_clinF_Anti, probs = seq(0, 1, by = 0.1), na.rm = TRUE)
quantiles[1] <- 0.99 * quantiles[1]
quantiles[length(quantiles)] <- 1.1 * quantiles[length(quantiles)]
processed_quantiles <- cut(
  temp_roc_data$model_clinF_Anti,
  breaks = quantiles,
  include_lowest = TRUE
)
processed_quantiles <- data.frame(y = temp_roc_data$progressed, pred = temp_roc_data$model_clinF_Anti, dec = processed_quantiles) %>%
  group_by(dec) %>%
  mutate(prob_obs = sum(y) / n(),
    obs = sum(y),
    n_group = n(),
    mnpred = mean(pred),
    lower = lapply(sum(y), prop.test, n = n()),
    upper = sapply(lower, function(x) x$conf.int[2]),
    lower = sapply(lower, function(x) x$conf.int[1]))
## plot
ggplot(processed_quantiles, aes(x = mnpred, y = prob_obs)) +
  geom_point() +
  xlab("Mean predicted probability in each decile") +
  ylab("Observed probability in each decile") +
  labs(title = "Clin F + Antibodies Model Cal. Plot") +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  ylim(c(0, 1)) + xlim(c(0, 1)) +
  geom_errorbar(aes(ymin = lower, ymax = upper))


```