12_prediction_modeling.R

rm(list = ls())
library(tibble)
library(tidyr)
library(readxl)
options(stringsAsFactors = FALSE)
options(width = 160)
set.seed(1)

###############################################################
##### Filter Non-Correlated Disease-specific FDR Proteins #####
###############################################################

Sumstats <- read_excel("WashU_pilot_data_explore/PredictionModel_03252025/Supplementary_tables.xlsx", sheet = 1, skip = 2)

new_colnames <- c(
  "Effect_size_AD", "SE_AD", "P_AD", "FDR_AD",
  "Effect_size_PD", "SE_PD", "P_PD", "FDR_PD",
  "Effect_size_FTD", "SE_FTD", "P_FTD", "FDR_FTD"
)

colnames(Sumstats)[6:(6 + length(new_colnames) - 1)] <- new_colnames

base_cols <- c("Aptamer", "UniProt", "Entrez ID", "Symbol", "Target")

# AD
AD_Sumstats <- Sumstats[, c(base_cols, "Effect_size_AD", "SE_AD", "P_AD", "FDR_AD")]
dim(AD_Sumstats) # 7289    9
AD_FDR <- na.omit(unique(AD_Sumstats[AD_Sumstats$FDR_AD < 0.05,]$Aptamer)) 

# PD
PD_Sumstats <- Sumstats[, c(base_cols, "Effect_size_PD", "SE_PD", "P_PD", "FDR_PD")]
dim(PD_Sumstats) # 7289    9
PD_FDR <- na.omit(unique(PD_Sumstats[PD_Sumstats$FDR_PD < 0.05,]$Aptamer))

# FTD
FTD_Sumstats <- Sumstats[, c(base_cols, "Effect_size_FTD", "SE_FTD", "P_FTD", "FDR_FTD")]
dim(FTD_Sumstats) # 7289    9
FTD_FDR <- na.omit(unique(FTD_Sumstats[FTD_Sumstats$FDR_FTD < 0.05,]$Aptamer))

length(AD_FDR) # 5187
length(FTD_FDR) # 2380
length(PD_FDR) # 3748
novel_proteins <- unique(na.omit(c(AD_FDR, FTD_FDR, PD_FDR))) 
length(novel_proteins) #6533

AD_expr <- read_csv("WashU_pilot_data_explore/PredictionModel_03252025/AD_expr.csv") # 7712 7299
FTD_expr <- read_csv("WashU_pilot_data_explore/PredictionModel_03252025/FTD_expr.csv") # 6428 7299
PD_expr <- read_csv("WashU_pilot_data_explore/PredictionModel_03252025/PD_expr.csv") # 6301 7299

rename_seq_columns <- function(df) {
  new_names <- names(df)
  new_names[11:length(new_names)] <- gsub("^seq\\.", "X", new_names[11:length(new_names)])
  names(df) <- new_names
  return(df)
}

AD_expr <- rename_seq_columns(AD_expr)
FTD_expr <- rename_seq_columns(FTD_expr)
PD_expr  <- rename_seq_columns(PD_expr)

table(novel_proteins %in% names(AD_expr)) # TRUE 6553 
table(novel_proteins %in% names(PD_expr)) # TRUE 6553 
table(novel_proteins %in% names(FTD_expr)) # TRUE 6553

analytes_corr_removal <- function(expression_data, novel_proteins) {
  expr_matrix <- expression_data[, novel_proteins]
  expr_matrix_nonNA <- as.data.frame(as.matrix(sapply(expr_matrix, as.numeric)))
  
  count.na <- apply(is.na(expr_matrix_nonNA), 2, sum)
  set.seed(1)
  for (i in which(count.na != 0)) {
    index <- is.na(expr_matrix_nonNA[, i])
    expr_matrix_nonNA[index, i] <- sample(expr_matrix_nonNA[!index, i], sum(index), replace = TRUE)
  }
  
  d2_cor <- expr_matrix_nonNA %>%
    as.matrix() %>%
    cor() %>%
    as.data.frame() %>%
    rownames_to_column(var = 'var1') %>%
    gather(var2, value, -var1)
  
  d2_cor <- na.omit(d2_cor)
  d2_cor <- d2_cor[d2_cor$var1 != d2_cor$var2, ]
  d2_cor <- d2_cor[order(d2_cor$value), ]
  
  print("Summary of correlation between analytes:")
  print(summary(d2_cor$value))
  
  correlated_Analytes <- d2_cor[abs(d2_cor$value) > 0.80, ]
  
  correlated_Analytes$temp <- apply(correlated_Analytes, 1, function(x) paste(sort(x), collapse = ""))
  correlated_Analytes <- correlated_Analytes[!duplicated(correlated_Analytes$temp), 1:3]
  
  analytes_noncor <- unique(novel_proteins[!(novel_proteins %in% correlated_Analytes$var1)])
  return(analytes_noncor)
}


# AD
dim(AD_expr) # 7712 7299
AD_FDR_nonCorr <- analytes_corr_removal(AD_expr, AD_FDR)
# [1] "Summary of correlation between analytes:"
# Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
# -0.752036 -0.164734  0.003602  0.014402  0.157071  0.996398 
length(AD_FDR_nonCorr) # 4021 (1166 analyte removed)

# FTD
dim(FTD_expr) # 6428 7299
FTD_FDR_nonCorr <- analytes_corr_removal(FTD_expr, FTD_FDR)
# [1] "Summary of correlation between analytes:"
# Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
# -0.740530 -0.187837 -0.001819  0.020311  0.172879  0.993367 
length(FTD_FDR_nonCorr) # 1791 (589 removed)

# PD
dim(PD_expr) # 6301 7299
PD_FDR_nonCorr <- analytes_corr_removal(PD_expr, PD_FDR)
# [1] "Summary of correlation between analytes:"
# Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
# -0.738255 -0.184609  0.001126  0.025016  0.173791  0.994858 
length(PD_FDR_nonCorr) # 2777 (971 removed)

save(AD_FDR_nonCorr, FTD_FDR_nonCorr, PD_FDR_nonCorr, file="WashU_pilot_data_explore/PredictionModel_03252025/Plasma_disease_specific_All_FDR_proteins_AD_PD_FTD_Overlapping_NonCorr.RData")


rm(list = ls())
library(dplyr) # dplyr_1.1.4
library(pROC) # pROC_1.18.5
library(ROCR) # ROCR_1.0-11
library(caret) # caret_6.0-94
library(glmnet) # glmnet_4.1-8
library(ggplot2) # ggplot2_3.5.1
library(readr)
library(readxl)

options(stringsAsFactors = FALSE)
options(width = 160)
set.seed(1)


Sumstats <- read_excel("WashU_pilot_data_explore/PredictionModel_03252025/Supplementary_tables.xlsx", sheet = 1, skip = 2)

new_colnames <- c(
  "Effect_size_AD", "SE_AD", "P_AD", "FDR_AD",
  "Effect_size_PD", "SE_PD", "P_PD", "FDR_PD",
  "Effect_size_FTD", "SE_FTD", "P_FTD", "FDR_FTD"
)

colnames(Sumstats)[6:(6 + length(new_colnames) - 1)] <- new_colnames

base_cols <- c("Aptamer", "UniProt", "Entrez ID", "Symbol", "Target")

# AD
AD_Sumstats <- Sumstats[, c(base_cols, "Effect_size_AD", "SE_AD", "P_AD", "FDR_AD")]
dim(AD_Sumstats) # 7289    9
AD_FDR <- na.omit(unique(AD_Sumstats[AD_Sumstats$FDR_AD < 0.05,]$Aptamer)) 

# PD
PD_Sumstats <- Sumstats[, c(base_cols, "Effect_size_PD", "SE_PD", "P_PD", "FDR_PD")]
dim(PD_Sumstats) # 7289    9
PD_FDR <- na.omit(unique(PD_Sumstats[PD_Sumstats$FDR_PD < 0.05,]$Aptamer))

# FTD
FTD_Sumstats <- Sumstats[, c(base_cols, "Effect_size_FTD", "SE_FTD", "P_FTD", "FDR_FTD")]
dim(FTD_Sumstats) # 7289    9
FTD_FDR <- na.omit(unique(FTD_Sumstats[FTD_Sumstats$FDR_FTD < 0.05,]$Aptamer))

length(AD_FDR) # 5187
length(FTD_FDR) # 2380
length(PD_FDR) # 3748
novel_proteins <- unique(na.omit(c(AD_FDR, FTD_FDR, PD_FDR))) 
length(novel_proteins) #6533

AD_expr <- read_csv("WashU_pilot_data_explore/PredictionModel_03252025/AD_expr.csv") # 7712 7299
FTD_expr <- read_csv("WashU_pilot_data_explore/PredictionModel_03252025/FTD_expr.csv") # 6428 7299
PD_expr <- read_csv("WashU_pilot_data_explore/PredictionModel_03252025/PD_expr.csv") # 6301 7299

rename_seq_columns <- function(df) {
  new_names <- names(df)
  new_names[11:length(new_names)] <- gsub("^seq\\.", "X", new_names[11:length(new_names)])
  names(df) <- new_names
  return(df)
}

AD_expr <- rename_seq_columns(AD_expr)
FTD_expr <- rename_seq_columns(FTD_expr)
PD_expr  <- rename_seq_columns(PD_expr)

convert_case_control <- function(df) {
  df$level <- ifelse(tolower(df$`Updated_Processed New Status`) == "co", 0, 1)
  colnames(df)[colnames(df) == "Updated_Processed New Status"] <- "status_for_analysis"
  return(df)
}

AD_expr <- convert_case_control(AD_expr)
PD_expr <- convert_case_control(PD_expr)
FTD_expr <- convert_case_control(FTD_expr)

# Read in the FDR proteins for each disease and are non-correlated
load("WashU_pilot_data_explore/PredictionModel_03252025/Plasma_disease_specific_All_FDR_proteins_AD_PD_FTD_Overlapping_NonCorr.RData")
length(AD_FDR_nonCorr) # 4021
length(FTD_FDR_nonCorr) # 1791
length(PD_FDR_nonCorr) # 2777

novel_proteins <- na.omit(unique(c(AD_FDR_nonCorr, FTD_FDR_nonCorr, PD_FDR_nonCorr)))
length(novel_proteins) # 5134

AD_expr <- AD_expr[,c("sample_id", "person_id", "contributor_code", "sequential_visit_number", "status_for_analysis","age_at_visit","sex", "level",novel_proteins)]
dim(AD_expr) # 7712 5142
row.names(AD_expr) <- AD_expr$sample_id

PD_expr <- PD_expr[,c("sample_id", "person_id", "contributor_code", "sequential_visit_number", "status_for_analysis","age_at_visit","sex", "level",novel_proteins)]
dim(PD_expr) # 6301 5142
row.names(PD_expr) <- PD_expr$sample_id

FTD_expr <- FTD_expr[,c("sample_id", "person_id", "contributor_code", "sequential_visit_number", "status_for_analysis","age_at_visit","sex", "level",novel_proteins)]
dim(FTD_expr) # 6428 5142
row.names(FTD_expr) <- FTD_expr$sample_id

# Impute the missing values
data_imputation <- function(expression_data) {
  expr_matrix <- expression_data[,9:ncol(expression_data)]
  expr_matrix_nonNA <- as.data.frame(as.matrix(sapply(expr_matrix, as.numeric)))
  count.na = apply(is.na(expr_matrix), 2, sum)
  set.seed(1)
  for (i in which(count.na != 0)) { # bootstrapping with replacement
    index = is.na(expr_matrix_nonNA[,i])
    expr_matrix_nonNA[index, i] = sample(expr_matrix_nonNA[!index, i], sum(index), replace = TRUE)
  }
  expression_data <- as.data.frame(cbind(expression_data[,1:8], expr_matrix_nonNA))
  print("any NA in the imputed data:")
  print(sum(is.na(expression_data[,9:ncol(expression_data)])))
  return(expression_data)
}

AD_expr_imputed <- data_imputation(AD_expr)
dim(AD_expr_imputed) # 7712 5142
any(is.na(AD_expr_imputed)) #False

PD_expr_imputed <- data_imputation(PD_expr)
dim(PD_expr_imputed) # 6301 5142
any(is.na(PD_expr_imputed)) 

FTD_expr_imputed <- data_imputation(FTD_expr)
FTD_expr_imputed <- FTD_expr_imputed[FTD_expr_imputed$age_at_visit != -1, ] # has one age == -1. removed
dim(FTD_expr_imputed) # 6427 5142
any(is.na(FTD_expr_imputed))

# AD_FDR , PD_FDR , FTD_FDR, AD_expr_imputed , PD_expr_imputed , FTD_expr_imputed

save(
  AD_FDR,
  PD_FDR,
  FTD_FDR,
  AD_expr_imputed,
  PD_expr_imputed,
  FTD_expr_imputed,
  file = "WashU_pilot_data_explore/PredictionModel_03252025/CSF_AD_PD_FTD_FDR_Sumstats_Expr_Imputed.RData"
)




######################################
######## Prediction Model ############
######################################

rm(list = ls())
library(dplyr) # dplyr_1.1.4
library(pROC) # pROC_1.18.5
library(ROCR) # ROCR_1.0-11
library(caret) # caret_6.0-94
library(glmnet) # glmnet_4.1-8
library(ggplot2) # ggplot2_3.5.1
library(readr)
library(readxl)

options(stringsAsFactors = FALSE)
options(width = 160)
set.seed(1)


Sumstats <- read_excel("WashU_pilot_data_explore/PredictionModel_03252025/Supplementary_tables.xlsx", sheet = 1, skip = 2)

new_colnames <- c(
  "Effect_size_AD", "SE_AD", "P_AD", "FDR_AD",
  "Effect_size_PD", "SE_PD", "P_PD", "FDR_PD",
  "Effect_size_FTD", "SE_FTD", "P_FTD", "FDR_FTD"
)

colnames(Sumstats)[6:(6 + length(new_colnames) - 1)] <- new_colnames

base_cols <- c("Aptamer", "UniProt", "Entrez ID", "Symbol", "Target")

# AD
AD_Sumstats <- Sumstats[, c(base_cols, "Effect_size_AD", "SE_AD", "P_AD", "FDR_AD")]
dim(AD_Sumstats) # 7289    9
AD_FDR <- na.omit(unique(AD_Sumstats[AD_Sumstats$FDR_AD < 0.05,]$Aptamer)) 

# PD
PD_Sumstats <- Sumstats[, c(base_cols, "Effect_size_PD", "SE_PD", "P_PD", "FDR_PD")]
dim(PD_Sumstats) # 7289    9
PD_FDR <- na.omit(unique(PD_Sumstats[PD_Sumstats$FDR_PD < 0.05,]$Aptamer))

# FTD
FTD_Sumstats <- Sumstats[, c(base_cols, "Effect_size_FTD", "SE_FTD", "P_FTD", "FDR_FTD")]
dim(FTD_Sumstats) # 7289    9
FTD_FDR <- na.omit(unique(FTD_Sumstats[FTD_Sumstats$FDR_FTD < 0.05,]$Aptamer))

length(AD_FDR) # 5187
length(FTD_FDR) # 2380
length(PD_FDR) # 3748
novel_proteins <- unique(na.omit(c(AD_FDR, FTD_FDR, PD_FDR))) 
length(novel_proteins) #6533

AD_expr <- read_csv("WashU_pilot_data_explore/PredictionModel_03252025/AD_expr.csv") # 7712 7299
FTD_expr <- read_csv("WashU_pilot_data_explore/PredictionModel_03252025/FTD_expr.csv") # 6428 7299
PD_expr <- read_csv("WashU_pilot_data_explore/PredictionModel_03252025/PD_expr.csv") # 6301 7299

rename_seq_columns <- function(df) {
  new_names <- names(df)
  new_names[11:length(new_names)] <- gsub("^seq\\.", "X", new_names[11:length(new_names)])
  names(df) <- new_names
  return(df)
}

AD_expr <- rename_seq_columns(AD_expr)
FTD_expr <- rename_seq_columns(FTD_expr)
PD_expr  <- rename_seq_columns(PD_expr)

convert_case_control <- function(df) {
  df$level <- ifelse(tolower(df$`Updated_Processed New Status`) == "co", 0, 1)
  colnames(df)[colnames(df) == "Updated_Processed New Status"] <- "status_for_analysis"
  return(df)
}

AD_expr <- convert_case_control(AD_expr)
PD_expr <- convert_case_control(PD_expr)
FTD_expr <- convert_case_control(FTD_expr)

# Read in the FDR proteins for each disease and are non-correlated
load("WashU_pilot_data_explore/PredictionModel_03252025/Plasma_disease_specific_All_FDR_proteins_AD_PD_FTD_Overlapping_NonCorr.RData")
length(AD_FDR_nonCorr) # 4021
length(FTD_FDR_nonCorr) # 1791
length(PD_FDR_nonCorr) # 2777

novel_proteins <- na.omit(unique(c(AD_FDR_nonCorr, FTD_FDR_nonCorr, PD_FDR_nonCorr)))
length(novel_proteins) # 5134

AD_expr <- AD_expr[,c("sample_id", "person_id", "contributor_code", "sequential_visit_number", "status_for_analysis","age_at_visit","sex", "level",novel_proteins)]
dim(AD_expr) # 7712 5142
row.names(AD_expr) <- AD_expr$sample_id

PD_expr <- PD_expr[,c("sample_id", "person_id", "contributor_code", "sequential_visit_number", "status_for_analysis","age_at_visit","sex", "level",novel_proteins)]
dim(PD_expr) # 6301 5142
row.names(PD_expr) <- PD_expr$sample_id

FTD_expr <- FTD_expr[,c("sample_id", "person_id", "contributor_code", "sequential_visit_number", "status_for_analysis","age_at_visit","sex", "level",novel_proteins)]
dim(FTD_expr) # 6428 5142
row.names(FTD_expr) <- FTD_expr$sample_id

# Impute the missing values
data_imputation <- function(expression_data) {
  expr_matrix <- expression_data[,9:ncol(expression_data)]
  expr_matrix_nonNA <- as.data.frame(as.matrix(sapply(expr_matrix, as.numeric)))
  count.na = apply(is.na(expr_matrix), 2, sum)
  set.seed(1)
  for (i in which(count.na != 0)) { # bootstrapping with replacement
    index = is.na(expr_matrix_nonNA[,i])
    expr_matrix_nonNA[index, i] = sample(expr_matrix_nonNA[!index, i], sum(index), replace = TRUE)
  }
  expression_data <- as.data.frame(cbind(expression_data[,1:8], expr_matrix_nonNA))
  print("any NA in the imputed data:")
  print(sum(is.na(expression_data[,9:ncol(expression_data)])))
  return(expression_data)
}

AD_expr_imputed <- data_imputation(AD_expr)
dim(AD_expr_imputed) # 7712 5142
any(is.na(AD_expr_imputed)) #False

PD_expr_imputed <- data_imputation(PD_expr)
dim(PD_expr_imputed) # 6301 5142
any(is.na(PD_expr_imputed)) 

FTD_expr_imputed <- data_imputation(FTD_expr)
FTD_expr_imputed <- FTD_expr_imputed[FTD_expr_imputed$age_at_visit != -1, ] # has one age == -1. removed
dim(FTD_expr_imputed) # 6427 5142
any(is.na(FTD_expr_imputed))

# AD_FDR , PD_FDR , FTD_FDR, AD_expr_imputed , PD_expr_imputed , FTD_expr_imputed

save(
  AD_FDR,
  PD_FDR,
  FTD_FDR,
  AD_expr_imputed,
  PD_expr_imputed,
  FTD_expr_imputed,
  file = "WashU_pilot_data_explore/PredictionModel_03252025/CSF_AD_PD_FTD_FDR_Sumstats_Expr_Imputed.RData"
)


#################
##### AD ########
#################

#####################################################
#### AD Lasso Model Generation (Iterative; n=50) ####
#####################################################

# subset for AD novel proteins
table(AD_expr_imputed$level)
#    0    1
#  5776 1936 
table(AD_expr_imputed$sex)
#   1    2 
# 3297 4415 

selected_cols <- c("sample_id", "level", "age_at_visit", "sex", AD_FDR_nonCorr)
length(selected_cols) # 4025 = 4 + 4021
AD_expr_subset <- AD_expr_imputed[, selected_cols]
dim(AD_expr_subset) # 7712 4025
row.names(AD_expr_subset) <- AD_expr_subset$sample_id

# comparison <- "ADvsCO"
AD_expr_subset[1:2,1:5]
#                           sample_id                               level     age_at_visit  sex   X8997.4
# 000323e3-4113-4404-9088-833c1fb08bd3 000323e3-4113-4404-9088-833c1fb08bd3     1           78   1 2.7534588
# 000bd2bf-2566-44d5-af4d-ef86e423fe82 000bd2bf-2566-44d5-af4d-ef86e423fe82     1           90   2 0.9741967

iterative_lasso_modeling <- function(expr_data, comparison) {
  #Define a "not in" opperator
  `%!in%` <- Negate(`%in%`)
  
  # Expression matrix
  quantMatrix <- expr_data[,5:ncol(expr_data)]
  dim(quantMatrix) 
  quantMatrix <- as.matrix(quantMatrix)
  quantMatrix[1:4,1:4]
  
  # Phenotype
  pheno <- expr_data[,1:4]
  colnames(pheno)[1] <- "sample_ID"
  dim(pheno) 
  head(pheno)
  
  #List the proteins
  proteins <- colnames(quantMatrix)
  length(proteins) 
  prot <- as.data.frame(proteins)
  names(prot) <- "Protein"
  protWeight_df <- as.data.frame(prot)
  colnames(protWeight_df) <- c("Proteins")
  dim(protWeight_df) 
  aucdf <- data.frame()
  
  for (i in 1:50)
  {
    set.seed(i)
    #Set up testing and training data
    trainSamples <- c(sample(pheno$sample_ID[pheno$level=='1'], round(length(pheno$level[pheno$level=='1'])*0.7)),
                      sample(pheno$sample_ID[pheno$level=='0'], round(length(pheno$level[pheno$level=='0'])*0.7)))
    trainX <- quantMatrix[rownames(quantMatrix)%in%trainSamples,]
    trainX <- as.matrix(trainX)
    trainY <- as.numeric(pheno$level[pheno$sample_ID %in% trainSamples])
    #Build a lasso model
    set.seed(567)
    cvOut <- cv.glmnet(trainX, trainY, family="binomial",alpha=1)
    lamMin <- cvOut$lambda.min
    #Rebuild the model
    lassoBest <- glmnet(trainX, trainY, family="binomial",alpha=1, lambda = lamMin)
    coefProt <- coef(lassoBest) %>% as.matrix() %>% as_tibble(rownames = "Proteins")
    sigProt <- coefProt[coefProt$s0 != 0,]
    protWeight_df <- merge(protWeight_df, sigProt, by="Proteins", all.x = T)
    testX <- quantMatrix[rownames(quantMatrix)%!in%trainSamples,]
    testY <- as.numeric(pheno$level[pheno$sample_ID %!in% trainSamples])
    preds <- as.data.frame(predict(lassoBest, newx = testX, type = "response"))
    t <- ci.auc(testY, as.numeric(preds$s0), conf.level = 0.9)
    t <- data.frame(low = round(t[1],3), AUC = round(t[2],3), high = round(t[3],3), Proteins = (nrow(sigProt)-1))
    aucdf <- rbind(aucdf,t)
    print(i)
  }
  # Save predictions and performance obtained from lasso model training
  dim(protWeight_df) 
  dim(aucdf) 
  # hist(aucdf$Proteins)
  print("summary of proteins # selected in each of the 50 iterations:")
  print(summary(aucdf$Proteins))
  write.csv(protWeight_df, paste0("ProteinWeights_50runs_", comparison, "_AgeSexPlatePCs.csv", sep=""), row.names = F)
  write.csv(aucdf, paste0("ModelPerformance_50runs_", comparison, "_AgeSexPlatePCs.csv", sep=""), row.names = F)
  
  col_names <- c("Proteins", seq(from=1,to=50))
  colnames(protWeight_df) <- col_names
  protWeight_df$count <- rowSums(!is.na(select(protWeight_df, 2:51)))
  dim(protWeight_df) 
  protWeight_df <- protWeight_df[order(protWeight_df$count, decreasing=T),]
  protWeight_df <<- protWeight_df[order(protWeight_df$count, decreasing=T),]
  print("summary of most important protein selection:")
  print(summary(protWeight_df$count))
  write.csv(protWeight_df, paste0("ProteinWeights_Sorted_50runs_", comparison, "_AgeSexPlatePCs.csv", sep=""), row.names = F)
  
  # Carlos asked me to try the model with top 10 analytes as they appeared in 45 of the 50 iterations.
  iterativeModel_10analytes <<- protWeight_df[1:10,]$Proteins 
  print("Top 10 proteins selected for final Lasso model:")
  print(iterativeModel_10analytes)
  save(iterativeModel_10analytes, file=paste0("iterativeModel_10analytes_", comparison,"_AgeSexPlatePCs.RData", sep=""))
}   

iterative_lasso_modeling(AD_expr_subset, "ADvsCO")
iterativeModel_10analytes <- c(
  "X15321.8",
  "X22782.80",
  "X10980.11",
  "X8997.4",
  "X15539.15",
  "X3461.58",
  "X3204.2",
  "X14136.234",
  "X19553.14",
  "X5383.14"
)


comparison <- "ADvsCO"

save(iterativeModel_10analytes, 
     file = paste0("iterativeModel_10analytes_", comparison, "_AgeSexPlatePCs.RData"))

###################################################################################
##### AD - Training and Testing (70 - 30%) in AD data using 10 analytes model #####
###################################################################################

library(dplyr) # dplyr_1.1.4
library(pROC) # pROC_1.18.5
library(ROCR) # ROCR_1.0-11
library(caret) # caret_6.0-94
library(glmnet) # glmnet_4.1-8
library(ggplot2) # ggplot2_3.5.1
options(stringsAsFactors = FALSE)
options(width = 160)
set.seed(1)

# comparison <- "ADvsCO"
AD_expr_subset[1:2,1:5]

#                           sample_id                                 level   age_at_visit sex   X8997.4
# 000323e3-4113-4404-9088-833c1fb08bd3 000323e3-4113-4404-9088-833c1fb08bd3     1           78   1 2.7534588
# 000bd2bf-2566-44d5-af4d-ef86e423fe82 000bd2bf-2566-44d5-af4d-ef86e423fe82     1           90   2 0.9741967

iterative_model_testing <- function(expr_data, comparison) {
  
  selected_cols <<- c("level", "age_at_visit", "sex", iterativeModel_10analytes)
  AD_expr_final <- expr_data[, selected_cols]
  
  print("Overall sample case-control distribution:")
  print(table(AD_expr_final$level))
  
  set.seed(1)
  ones <- AD_expr_final %>% filter(level == 1)
  zeros <- AD_expr_final %>% filter(level == 0)
  ones <- ones[sample(nrow(ones)), ]
  zeros <- zeros[sample(nrow(zeros)), ]
  
  train_ones <- ones[1:round(nrow(ones) * 0.70, 0), ]
  test_ones <- ones[(round(nrow(ones) * 0.70, 0) + 1):nrow(ones), ]
  train_zeros <- zeros[1:round(nrow(zeros) * 0.70, 0), ]
  test_zeros <- zeros[(round(nrow(zeros) * 0.70, 0) + 1):nrow(zeros), ]
  
  train_discovery <- bind_rows(train_ones, train_zeros)
  test_discovery <- bind_rows(test_ones, test_zeros)
  train_discovery <<- train_discovery[sample(nrow(train_discovery)), ]
  test_discovery <<- test_discovery[sample(nrow(test_discovery)), ]
  
  print("Any train sample present in test data:")
  print(table(row.names(train_discovery) %in% row.names(test_discovery)))
  
  print("Dim of train data:")
  print(dim(train_discovery))
  print(table(train_discovery$level))
  print("Dim of test data:")
  print(dim(test_discovery))
  print(table(test_discovery$level))
  
  bl_disc <<- glm(as.formula(level ~ age_at_visit + sex), data = train_discovery, family = 'binomial')
  out_prob_bl <- predict(bl_disc, newdata = test_discovery, type = 'response')
  auc_test_bl <<- roc(as.formula(test_discovery$level ~ out_prob_bl), plot = FALSE, print.auc = FALSE)
  npv_ppv_best <- coords(auc_test_bl, x = "best", best.method = "youden", input = "threshold", ret = c("threshold", "acc", "npv", "ppv", "sensitivity", "specificity"))
  print("Performance for baseline model:")
  print(npv_ppv_best)
  print(round(auc_test_bl$auc, 2))
  
  all_disc <<- glm(level ~ ., data = train_discovery, family = 'binomial')
  out_prob <- predict(all_disc, newdata = test_discovery, type = 'response')
  auc_test_10analytes <<- roc(as.formula(test_discovery$level ~ out_prob), plot = FALSE, print.auc = FALSE)
  npv_ppv_best <- coords(auc_test_10analytes, x = "best", best.method = "youden", input = "threshold", ret = c("threshold", "acc", "npv", "ppv", "sensitivity", "specificity"))
  print("Performance for 10-analyte model:")
  print(npv_ppv_best)
  print(round(auc_test_10analytes$auc, 2))
  print("Fixed threshold:")
  fixed_threshold <<- npv_ppv_best$threshold
  print(fixed_threshold)
  
  save(all_disc, bl_disc, train_discovery, test_discovery, iterativeModel_10analytes, out_prob, fixed_threshold, file=paste0(comparison, "_10Analytes_Model_AgeSexPlatePCs.RData", sep=""))
} 

iterative_model_testing(AD_expr_subset, "ADvsCO")
# [1] "Overall sample case-control distribution:"
# 
# 0    1 
# 5776 1936 
# [1] "Any train sample present in test data:"
# 
# FALSE 
# 5398 
# [1] "Dim of train data:"
# [1] 5398   13
# 
# 0    1 
# 4043 1355 
# [1] "Dim of test data:"
# [1] 2314   13
# 
# 0    1 
# 1733  581 
# Setting levels: control = 0, case = 1
# Setting direction: controls < cases
# [1] "Performance for baseline model:"
#           threshold  accuracy       npv      ppv sensitivity specificity
# threshold 0.2100544 0.5090752 0.8754717 0.317314   0.8296041   0.4016157
# [1] 0.63
# Setting levels: control = 0, case = 1
# Setting direction: controls < cases
# [1] "Performance for 10-analyte model:"
#           threshold  accuracy       npv       ppv sensitivity specificity
# threshold 0.2669046 0.7532411 0.8930988 0.5059809   0.7280551   0.7616849
# [1] 0.82
# [1] "Fixed threshold:"
# [1] 0.2669046

# Entire figure at one place
png("AUC_Lasso_10_IterativeModel_ADvsCO_AgeSexPlatePCs.png", units="mm", width=120, height=120, res=1000)
plot_curve <- plot.roc(auc_test_bl, xlim=c(1,0), ylim=c(0, 1), col=alpha("blue", 0.5), lty=2, 
                       main = "ADvsCO")
text(0.6, 0.20, paste("Baseline: ", round(auc_test_bl$auc, 2), sep=""), col=alpha("blue", 0.5), cex=0.75, pos=4)
plot_curve <- plot.roc(auc_test_10analytes, xlim=c(1,0), ylim=c(0, 1), lty=1, lwd=6, col=alpha("blue", 0.5), add=TRUE)
text(0.6, 0.25, paste("Test Data (30%): ",  round(auc_test_10analytes$auc, 2), sep=""), col=("blue"), cex=0.85, pos=4)
dev.off()
save(auc_test_bl, auc_test_10analytes, file="AUC_Lasso_10_IterativeModel_ADvsCO_AgeSexPlatePCs.RData")

# get AD probablity for all individuals
out_prob_AD <- predict(all_disc, newdata = rbind(test_discovery, train_discovery), type='response')
out_prob_AD <- as.data.frame(out_prob_AD)
out_prob_AD$sample_id <- row.names(out_prob_AD)
dim(out_prob_AD) # 7712    2
out_prob_AD <- inner_join(out_prob_AD, AD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_AD, sample_id, status_for_analysis)
dim(out_prob_AD) # 7712    3
table(out_prob_AD$out_prob_AD > fixed_threshold, out_prob_AD$status_for_analysis)
#         ad   CO
# FALSE  538 4434
# TRUE  1398 1342


###########################################
##### AD - Testing in other disorders #####
###########################################
## PD ##
dim(PD_expr_imputed) # 6301 5142
PD_expr_imputed[1:4,1:13]
table(PD_expr_imputed$level)
# 0    1 
# 5776  525 
table(PD_expr_imputed$sex)
#   1    2 
# 2696 3605 
PD_expr_subset <- PD_expr_imputed[, selected_cols]
row.names(PD_expr_subset) <- PD_expr_imputed$sample_id
dim(PD_expr_subset) # 6301   13
# Testing for all analytes in PD (N=10)
out_prob_PD <- predict(all_disc, newdata = PD_expr_subset, type='response')
auc_test_PD_10analytes <- roc(as.formula(PD_expr_subset$level ~ out_prob_PD), plot = FALSE, print.auc = FALSE)
npv_ppv_best <- coords(auc_test_PD_10analytes, x = "all", best.method="youden", input = "threshold", ret = c("threshold","acc","npv","ppv","sensitivity","specificity"))
npv_ppv_best[which.min(abs(fixed_threshold-npv_ppv_best$threshold)),]
#       threshold  accuracy       npv        ppv sensitivity specificity
# 1420 0.2668937 0.2796382 0.9358703 0.08889799   0.8266667   0.2299169
round(auc_test_PD_10analytes$auc, 2) # 0.53
out_prob_PD <- as.data.frame(out_prob_PD)
sum(is.na(out_prob_PD)) # 0
out_prob_PD$sample_id <- row.names(out_prob_PD)
out_prob_PD <- inner_join(out_prob_PD, PD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_PD, sample_id, status_for_analysis)
table(out_prob_PD$out_prob_PD > fixed_threshold, out_prob_PD$status_for_analysis)
#         CO   pd
# FALSE 4449  434
# TRUE  1327   91

## FTD ##
dim(FTD_expr_imputed) #  6427 5142
FTD_expr_imputed[1:4,1:13]
table(FTD_expr_imputed$level)
#   0    1 
# 6264  163 
table(FTD_expr_imputed$sex)
#   1    2 
# 2734 3693 
sum(is.na(FTD_expr_imputed)) # 0
FTD_expr_subset <- FTD_expr_imputed[, selected_cols]
row.names(FTD_expr_subset) <- FTD_expr_imputed$sample_id
dim(FTD_expr_subset) # 6427   13
# Testing for all analytes in FTD (N=10)
out_prob_FTD <- predict(all_disc, newdata = FTD_expr_subset, type='response')
auc_test_FTD_10analytes <- roc(as.formula(FTD_expr_subset$level ~ out_prob_FTD), plot = FALSE, print.auc = FALSE)
npv_ppv_best <- coords(auc_test_FTD_10analytes, x = "all", best.method="youden", input = "threshold", ret = c("threshold","acc","npv","ppv","sensitivity","specificity"))
npv_ppv_best[which.min(abs(fixed_threshold-npv_ppv_best$threshold)),]
#     threshold  accuracy       npv        ppv sensitivity specificity
# 5005 0.2668937 0.7662984 0.9758193 0.02951511   0.2576687   0.7795338
round(auc_test_FTD_10analytes$auc, 2) # 0.57
out_prob_FTD <- as.data.frame(out_prob_FTD)
sum(is.na(out_prob_FTD)) # 0
out_prob_FTD$sample_id <- row.names(out_prob_FTD)
out_prob_FTD <- inner_join(out_prob_FTD, FTD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_FTD, sample_id, status_for_analysis)
table(out_prob_FTD$out_prob_FTD > fixed_threshold, out_prob_FTD$status_for_analysis)
#         CO  ftd
# FALSE 4884  121
# TRUE  1380   42

save(AD_expr_imputed, PD_expr_imputed, FTD_expr_imputed, file="CSF_Complete_DiseaseSpecific_Imputed_Expression_DFs_For_CSFPredModel_AgeSexPlatePCs.RData")

colnames(out_prob_PD)[1] <- "out_prob_AD"
colnames(out_prob_FTD)[1] <- "out_prob_AD"
AD_Probablities <- rbind(out_prob_AD, out_prob_PD, out_prob_FTD)
dim(AD_Probablities) # 20440     3
names(AD_Probablities)[3] <- "Final_Status"
length(unique(AD_Probablities$sample_id)) # 8888

AD_Probablities <- AD_Probablities[order(AD_Probablities$sample_id, AD_Probablities$out_prob_AD),]
AD_Probablities <- AD_Probablities[!duplicated(AD_Probablities$sample_id),]
AD_Probablities <- na.omit(AD_Probablities)
dim(AD_Probablities) # 8888    3
table(AD_Probablities$Final_Status)
# ad   CO  ftd   pd 
# 1936 6264  163  525 
save(AD_Probablities, file="CSF_AD_Probability_N8888_AgeSexPlatePCs.RData")

png("Probability_Violinplot_for_CSF_AD_Signature_AgeSexPlatePCs.png", units="mm", width=100, height=100, res=1000)
ggplot(AD_Probablities, aes(x=Final_Status, y=out_prob_AD, color=Final_Status)) + geom_violin(trim = TRUE) + 
  theme(legend.position = "none", panel.background = element_blank(), axis.line = element_line(colour = "black")) +
  labs(x ="Clinical Diagnosis", y = "Predicted AD Probability") + ggtitle("Probability Distribution for AD signature")
dev.off()

png("AUC_Lasso_10_IterativeModel_ADvsCO_Testing_Other_NDs_AgeSexPlatePCs.png", units="mm", width=120, height=120, res=1000)
plot_curve <- plot.roc(auc_test_PD_10analytes, xlim=c(1,0), ylim=c(0, 1), lty=1, lwd=6, col=alpha("red", 0.5), 
                       main = c("AD Model Testing across other dementia"))
text(0.6, 0.15, paste("PD: ",  round(auc_test_PD_10analytes$auc, 2), sep=""), col=("red"), cex=0.85, pos=4)
plot_curve <- plot.roc(auc_test_FTD_10analytes, xlim=c(1,0), ylim=c(0, 1), lty=1, lwd=6, col=alpha("green4", 0.5), add=TRUE)
text(0.6, 0.10, paste("FTD: ",  round(auc_test_FTD_10analytes$auc, 2), sep=""), col=("green4"), cex=0.80, pos=4)
dev.off()
save(auc_test_PD_10analytes, auc_test_FTD_10analytes, file="AUC_Lasso_10_IterativeModel_ADvsCO_Testing_Other_NDs_AgeSexPlatePCs.RData")


##########################################################################################
##### AD Prediction Model (Iterative Class Balanced Training {70%} and Testing{30%}) #####
##########################################################################################
rm(list = ls())
library(dplyr) # dplyr_1.1.4
library(pROC) # pROC_1.18.5
library(ROCR) # ROCR_1.0-11
library(caret) # caret_6.0-94
library(glmnet) # glmnet_4.1-8
options(stringsAsFactors = FALSE)
options(width = 160)
set.seed(1)

load("~/files/WashU_pilot_data_explore/PredictionModel_03252025/AD/CSF_Complete_DiseaseSpecific_Imputed_Expression_DFs_For_CSFPredModel_AgeSexPlatePCs.RData")

selected_cols <- c("level", "age_at_visit", "sex", iterativeModel_10analytes)
length(selected_cols) # 13

row.names(AD_expr_imputed) <- AD_expr_imputed$sample_id
AD_expr_subset <- AD_expr_imputed[, selected_cols]
dim(AD_expr_subset) # 7712   13

row.names(FTD_expr_imputed) <- FTD_expr_imputed$sample_id
FTD_expr_subset <- FTD_expr_imputed[, selected_cols]
dim(FTD_expr_subset) # 6427   13

row.names(PD_expr_imputed) <- PD_expr_imputed$sample_id
PD_expr_subset <- PD_expr_imputed[, selected_cols]
dim(PD_expr_subset) # 6301   13


# for saving performance info from the model
perf_df_ad <- c()
perf_df_ad_base <- c()
# for saving predictions from the model
iter_pred_ad_base <- list()
iter_pred_ad <- list()
iter_y_ad <- list()

# other dementias
perf_df_DLB <- c()
iter_pred_DLB <- list()
iter_y_DLB <- list()
perf_df_FTD <- c()
iter_pred_FTD <- list()
iter_y_FTD <- list()
perf_df_PD <- c()
iter_pred_PD <- list()
iter_y_PD <- list()

# round(table(AD_expr_subset$level)[1]*0.70, 0) # 4043
# round(table(AD_expr_subset$level)[1]*0.30, 0) # 1733
table(AD_expr_subset$level)
# 0    1 
# 5776 1936 
min(table(AD_expr_subset$level)) # 1936


for (i in 1:100) {
  set.seed(i)
  print(i)
  # Separate data by levels # In each iteration 70% data is training (equal number of AD and CO) and 30% data is testing (equal number of AD and CO)
  level_0 <- AD_expr_subset[AD_expr_subset$level == 0, ]
  level_1 <- AD_expr_subset[AD_expr_subset$level == 1, ]
  # Sample rows for level = 1
  train_level_1 <- level_1[sample(1:nrow(level_1), size = round(0.7 * min(table(AD_expr_subset$level)))), ]
  test_level_1 <- level_1[!rownames(level_1) %in% rownames(train_level_1), ] # maker sure test has no train samples
  test_level_1 <- test_level_1[sample(1:nrow(test_level_1), size = round(0.3 * min(table(AD_expr_subset$level)))), ] # downsample test dataset
  # Sample rows for level = 0
  train_level_0 <- level_0[sample(1:nrow(level_0), size = nrow(train_level_1)), ]
  test_level_0 <- level_0[!rownames(level_0) %in% rownames(train_level_0), ] # maker sure test has no train samples
  test_level_0 <- test_level_0[sample(1:nrow(test_level_0), size = nrow(test_level_1)), ] # downsample test dataset
  # Combine training and testing sets
  Train_set <- rbind(train_level_1, train_level_0)
  Test_set <- rbind(test_level_1, test_level_0)
  
  # Baseline model
  bl_disc <- glm(as.formula(level ~ age_at_visit + sex), data = Train_set, family = 'binomial')
  base_pred_val_ad <- predict(bl_disc, newdata = Test_set, type='response')
  ROC_base_val_ad <- roc(as.formula(Test_set$level ~ base_pred_val_ad), plot = FALSE, print.auc = FALSE)
  
  # 10-protein prediction model
  all_disc <<- glm(level ~ ., data = Train_set, family = 'binomial')
  pred_val_glm_ad <- predict(all_disc, newdata = Test_set, type='response')
  ROC_val_ad <- roc(as.formula(Test_set$level ~ pred_val_glm_ad), plot = FALSE, print.auc = FALSE)
  
  iter_pred_ad_base[[i]] <- as.numeric(base_pred_val_ad)
  iter_pred_ad[[i]] <- as.numeric(pred_val_glm_ad)
  iter_y_ad[[i]] <- Test_set$level
  
  # each performance result for iterations, if there is any ties, get the first one only
  perf_df_ad <- rbind(perf_df_ad, 
                      coords(ROC_val_ad, x = "best", best.method="youden", input = "threshold", 
                             ret = c("threshold","auc", "acc","npv","ppv","sensitivity","specificity")) %>% 
                        mutate(iteration = i, auc = ROC_val_ad$auc, .before = 1) %>% 
                        summarise_all(first))
  perf_df_ad_base <- rbind(perf_df_ad_base, 
                           coords(ROC_base_val_ad, x = "best", best.method="youden", input = "threshold", 
                                  ret = c("threshold","auc", "acc","npv","ppv","sensitivity","specificity")) %>% 
                             mutate(iteration = i, auc = ROC_base_val_ad$auc, .before = 1) %>% 
                             summarise_all(first))
  # FTD
  pred_val_glm_FTD <- predict(all_disc, newdata = FTD_expr_subset, type='response')
  ROC_val_FTD <- roc(as.formula(FTD_expr_subset$level ~ pred_val_glm_FTD), plot = FALSE, print.auc = FALSE)
  iter_pred_FTD[[i]] <- as.numeric(pred_val_glm_FTD)
  iter_y_FTD[[i]] <- FTD_expr_subset$level
  perf_df_FTD <- rbind(perf_df_FTD, 
                       coords(ROC_val_FTD, x = "best", best.method="youden", input = "threshold", 
                              ret = c("threshold","auc", "acc","npv","ppv","sensitivity","specificity")) %>% 
                         mutate(iteration = i, auc = ROC_val_FTD$auc, .before = 1) %>%  summarise_all(first))
  # PD
  pred_val_glm_PD <- predict(all_disc, newdata = PD_expr_subset, type='response')
  ROC_val_PD <- roc(as.formula(PD_expr_subset$level ~ pred_val_glm_PD), plot = FALSE, print.auc = FALSE)
  iter_pred_PD[[i]] <- as.numeric(pred_val_glm_PD)
  iter_y_PD[[i]] <- PD_expr_subset$level
  perf_df_PD <- rbind(perf_df_PD, 
                      coords(ROC_val_PD, x = "best", best.method="youden", input = "threshold", 
                             ret = c("threshold","auc", "acc","npv","ppv","sensitivity","specificity")) %>% 
                        mutate(iteration = i, auc = ROC_val_PD$auc, .before = 1) %>% summarise_all(first))
}
roc.test(ROC_base_val_ad, ROC_val_ad, method="delong") # Z = -10.214, p-value < 2.2e-16, 95 percent confidence interval: -0.2107188 -0.1435639
roc.test(ROC_val_ad, ROC_val_FTD, method="delong") # p-value < 2.2e-16
roc.test(ROC_val_ad, ROC_val_PD, method="delong") # p-value < 2.2e-16

pred_ad <- prediction(iter_pred_ad, iter_y_ad)
pred_ad_base <- prediction(iter_pred_ad_base, iter_y_ad)
perf_ad <- performance(pred_ad, "tpr", "fpr")
perf_ad_base <- performance(pred_ad_base, "tpr", "fpr")
mean(perf_df_ad$auc) # 0.8110871
sd(perf_df_ad$auc) # 0.009476067
mean(perf_df_ad_base$auc) # 0.6295539
sd(perf_df_ad_base$auc) # 0.01653251
summary(perf_df_ad$auc)
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 0.7902  0.8048  0.8100  0.8111  0.8171  0.8361 

dim(perf_df_ad) # 100   8
dim(perf_df_ad_base) # 100   8
# get AD probablity for all individuals
out_prob_AD <- predict(all_disc, newdata = AD_expr_subset, type='response')
out_prob_AD <- as.data.frame(out_prob_AD)
out_prob_AD$sample_id <- row.names(out_prob_AD)
dim(out_prob_AD) # 7712    2

out_prob_AD_mean <- inner_join(out_prob_AD, AD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_AD, sample_id, status_for_analysis)
dim(out_prob_AD_mean) # 7712    3

plot(perf_ad, avg = "vertical", 
     spread.estimate="stderror", lwd=6, lty=1, plotCI.col = alpha("blue", 0.5), col = alpha("blue", 0.5), main = c("ADvsCO"))
text(0.4, 0.25, paste("Test Data (30%): ",  round(mean(perf_df_ad$auc), 2), sep=""), col=("blue"), cex=0.85, pos=4)
plot(perf_ad_base, avg = "vertical", 
     spread.estimate="stderror", lty=2, add = T, plotCI.col = alpha("blue", 0.5), col = alpha("blue", 0.5))
text(0.4, 0.20, paste("Baseline: ", round(mean(perf_df_ad_base$auc), 2), sep=""), col=alpha("blue", 0.5), cex=0.75, pos=4)
save(perf_ad, perf_df_ad, perf_ad_base, perf_df_ad_base, 
     file="AUC_Lasso_10_IterativeModel_ADvsCO_AgeSexProjectPCs_Iterative_ClassBalanced.RData")

# FTD
pred_FTD <- prediction(iter_pred_FTD, iter_y_FTD)
perf_FTD <- performance(pred_FTD, "tpr", "fpr")
mean(perf_df_FTD$auc) # 0.5418467
sd(perf_df_FTD$auc) # 0.01139694
dim(perf_df_FTD) # 100   8
out_prob_FTD_mean <- sapply(1:length(iter_pred_FTD[[1]]), function(i) {
  mean(sapply(iter_pred_FTD, function(x) x[i]))
})
out_prob_FTD_mean <- as.data.frame(cbind(row.names(FTD_expr_subset), out_prob_FTD_mean))
colnames(out_prob_FTD_mean) <- c("sample_id", "out_prob_FTD")
out_prob_FTD_mean <- inner_join(out_prob_FTD_mean, FTD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_FTD, sample_id, status_for_analysis)

# PD
pred_PD <- prediction(iter_pred_PD, iter_y_PD)
perf_PD <- performance(pred_PD, "tpr", "fpr")
mean(perf_df_PD$auc) # 0.5396055
sd(perf_df_PD$auc) # 0.005534958
dim(perf_df_PD) # 100   8
out_prob_PD_mean <- sapply(1:length(iter_pred_PD[[1]]), function(i) {
  mean(sapply(iter_pred_PD, function(x) x[i]))
})
out_prob_PD_mean <- as.data.frame(cbind(row.names(PD_expr_subset), out_prob_PD_mean))
colnames(out_prob_PD_mean) <- c("sample_id", "out_prob_PD")
out_prob_PD_mean <- inner_join(out_prob_PD_mean, PD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_PD, sample_id, status_for_analysis)

png("AUC_ADvsCO_Testing_FTD_PD_Iterative.png", units = "mm", width = 120, height = 120, res = 1000)

plot(perf_FTD, avg = "vertical", spread.estimate = "stderror", lwd = 6,
     plotCI.col = alpha("green4", 0.5), col = alpha("green4", 0.5),
     main = "AD model testing across other dementia")
text(0.4, 0.10, paste("FTD: ", round(mean(perf_df_FTD$auc), 2)), col = "green4", cex = 0.80, pos = 4)
plot(perf_PD, avg = "vertical", spread.estimate = "stderror", lwd = 6, add = TRUE,
     plotCI.col = alpha("red", 0.5), col = alpha("red", 0.5))
text(0.4, 0.15, paste("PD: ", round(mean(perf_df_PD$auc), 2)), col = "red", cex = 0.80, pos = 4)
dev.off()
save(perf_FTD, perf_df_FTD, perf_PD, perf_df_PD, 
     file="AUC_Lasso_10_IterativeModel_ADvsCO_Testing_Other_NDs_AgeSexProjectPCs_Iterative_ClassBalanced.RData")

colnames(out_prob_FTD_mean)[1] <- "out_prob_AD"
colnames(out_prob_PD_mean)[1] <- "out_prob_AD"
AD_Probablities <- rbind(out_prob_AD_mean, out_prob_FTD_mean, out_prob_PD_mean)
dim(AD_Probablities) # 20440     3
names(AD_Probablities)[3] <- "Final_Status"
length(unique(AD_Probablities$sample_id)) # 8888
AD_Probablities <- AD_Probablities[order(AD_Probablities$sample_id, AD_Probablities$out_prob_AD),]
AD_Probablities <- AD_Probablities[!duplicated(AD_Probablities$sample_id),]
AD_Probablities <- na.omit(AD_Probablities)
AD_Probablities$out_prob_AD <- as.numeric(as.character(AD_Probablities$out_prob_AD))
dim(AD_Probablities) # 8888    3
table(AD_Probablities$Final_Status)
# ad   CO  ftd   pd 
# 1936 6264  163  525 

save(AD_Probablities, file="Plasma_AD_Probability_N8888_AgeSexPlatePCs_Iterative_ClassBalanced.RData")
png("Probability_Violinplot_for_Plasma_AD_Signature_AgeSexProjectPCs_Iterative_ClassBalanced.png", units="mm", width=100, height=100, res=1000)
ggplot(AD_Probablities, aes(x=Final_Status, y=out_prob_AD, color=Final_Status)) + geom_violin(trim = TRUE) + 
  theme(legend.position = "none", panel.background = element_blank(), axis.line = element_line(colour = "black")) +
  labs(x ="Clinical Diagnosis", y = "Predicted AD Probability") + ggtitle("Probability Distribution for AD signature")
dev.off()

# Calculate mean and standard deviation for each column
calculate_mean_sd <- function(df) {
  df <- df[,-1]
  # Ensure only numeric columns are included
  numeric_cols <- df[sapply(df, is.numeric)]
  # Calculate mean and SD for each numeric column
  formatted_stats <- sapply(numeric_cols, function(x) {
    mean_val <- mean(x, na.rm = TRUE)
    sd_val <- sd(x, na.rm = TRUE)
    paste0(round(mean_val, 2), "(", round(sd_val, 2), ")")
  })
  # Convert to a single-row dataframe
  result <- as.data.frame(t(formatted_stats))
  colnames(result) <- colnames(numeric_cols)
  return(result)
}
calculate_mean_sd(perf_df_ad)
#           auc  threshold   accuracy        npv        ppv sensitivity specificity
# 1 0.81(0.01) 0.52(0.04) 0.74(0.01) 0.74(0.02) 0.76(0.02)  0.72(0.04)  0.76(0.04)
calculate_mean_sd(perf_df_ad_base)
#           auc  threshold   accuracy        npv        ppv sensitivity specificity
# 1 0.63(0.02) 0.47(0.03) 0.61(0.01) 0.66(0.03) 0.59(0.01)  0.75(0.07)  0.47(0.07)
calculate_mean_sd(perf_df_FTD)
#           auc threshold   accuracy     npv     ppv sensitivity specificity
# 1 0.54(0.01) 0.2(0.03) 0.34(0.04) 0.98(0) 0.03(0)  0.79(0.04)  0.32(0.04)
calculate_mean_sd(perf_df_PD)
#           auc  threshold   accuracy     npv     ppv sensitivity specificity
# 1 0.54(0.01) 0.39(0.05) 0.41(0.06) 0.93(0) 0.09(0)  0.69(0.07)  0.39(0.07)


#################
##### PD ########
#################

#####################################################
#### PD Lasso Model Generation (Iterative; n=50) ####
#####################################################
rm(list = ls())
library(dplyr) # dplyr_1.1.4
library(pROC) # pROC_1.18.5
library(ROCR) # ROCR_1.0-11
library(caret) # caret_6.0-94
library(glmnet) # glmnet_4.1-8
options(stringsAsFactors = FALSE)
options(width = 160)
set.seed(1)

# subset for PD novel proteins
table(PD_expr_imputed$level)
#    0    1
#  5776 525 
table(PD_expr_imputed$sex)
#   1    2 
# 2696 3605

# subset for PD novel proteins
selected_cols <- c("sample_id", "level", "age_at_visit", "sex", PD_FDR_nonCorr)
length(selected_cols) # 2781
PD_expr_subset <- PD_expr_imputed[, selected_cols]
dim(PD_expr_subset) # 6301 2781
row.names(PD_expr_subset) <- PD_expr_subset$sample_id

# comparison <- "PDvsCO"
PD_expr_subset[1:2,1:5]
#                                                     sample_id              level age_at_visit sex   X8997.4
# 0014d22a-aa65-4147-90f0-f6509c875e10 0014d22a-aa65-4147-90f0-f6509c875e10     1           77   2 0.9984053
# 00296564-0852-46f1-9de9-1fcb7b461e44 00296564-0852-46f1-9de9-1fcb7b461e44     0           89   2 1.7910329

iterative_lasso_modeling(PD_expr_subset, "PDvsCO")   

# [1] "summary of proteins # selected in each of the 50 iterations:"
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 228.0   267.2   292.5   292.4   312.2   365.0 
# [1] "summary of most important protein selection:"
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 0.000   0.000   0.000   5.264   4.000  50.000 

# [1] "Top 10 proteins selected for final Lasso model:"
# [1] "X10064.12"  "X10967.12"  "X12352.70"  "X12549.33"  "X12881.17"  "X15522.2"   "X16859.100" "X17796.15"  "X18310.26" 
# [10] "X23580.2" 

###################################################################################
##### PD - Training and Testing (70 - 30%) in AD data using 10 analytes model #####
###################################################################################
iterative_model_testing(PD_expr_subset, "PDvsCO")
# 
# [1] "Overall sample case-control distribution:"
# 
# 0    1 
# 5776  525 
# [1] "Any train sample present in test data:"
# 
# FALSE 
# 4411 
# [1] "Dim of train data:"
# [1] 4411   13
# 
# 0    1 
# 4043  368 
# [1] "Dim of test data:"
# [1] 1890   13
# 
# 0    1 
# 1733  157 
# Setting levels: control = 0, case = 1
# Setting direction: controls < cases
# [1] "Performance for baseline model:"
# threshold  accuracy       npv        ppv sensitivity specificity
# threshold 0.08561892 0.5820106 0.9290264 0.09987358   0.5031847   0.5891518
# [1] 0.56
# Setting levels: control = 0, case = 1
# Setting direction: controls < cases
# [1] "Performance for 10-analyte model:"
# threshold  accuracy       npv       ppv sensitivity specificity
# threshold 0.0486701 0.6433862 0.9783198 0.1698595   0.8471338   0.6249279
# [1] 0.82
# [1] "Fixed threshold:"
# [1] 0.0486701

# Entire figure at one place
png("AUC_Lasso_10_IterativeModel_PDvsCO_AgeSexPlatePCs.png", units="mm", width=120, height=120, res=1000)
plot_curve <- plot.roc(auc_test_bl, xlim=c(1,0), ylim=c(0, 1), col=alpha("red", 0.5), lty=2, 
                       main = "PDvsCO")
text(0.6, 0.20, paste("Baseline: ", round(auc_test_bl$auc, 2), sep=""), col=alpha("red", 0.5), cex=0.75, pos=4)
plot_curve <- plot.roc(auc_test_10analytes, xlim=c(1,0), ylim=c(0, 1), lty=1, lwd=6, col=alpha("red", 0.5), add=TRUE)
text(0.6, 0.25, paste("Test Data (30%): ",  round(auc_test_10analytes$auc, 2), sep=""), col=("red"), cex=0.85, pos=4)
dev.off()
save(auc_test_bl, auc_test_10analytes, file="AUC_Lasso_10_IterativeModel_PDvsCO_AgeSexPlatePCs.RData")

# get PD probablity for all individuals
out_prob_PD <- predict(all_disc, newdata = rbind(test_discovery, train_discovery), type='response')
out_prob_PD <- as.data.frame(out_prob_PD)
out_prob_PD$sample_id <- row.names(out_prob_PD)
dim(out_prob_PD) # 6301    2
out_prob_PD <- inner_join(out_prob_PD, PD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_PD, sample_id, status_for_analysis)
dim(out_prob_PD) # 6301    3
table(out_prob_PD$out_prob_PD > fixed_threshold, out_prob_PD$status_for_analysis)
#         CO   pd
# FALSE 3620   79
# TRUE  2156  446


###########################################
##### PD - Testing in other disorders #####
###########################################
## AD
AD_expr_subset <- AD_expr_imputed[, selected_cols]
row.names(AD_expr_subset) <- AD_expr_imputed$sample_id
dim(AD_expr_subset) # 7712   13
# Testing for all analytes in AD (N=10)
out_prob_AD <- predict(all_disc, newdata = AD_expr_subset, type='response')
auc_test_AD_10analytes <- roc(as.formula(AD_expr_subset$level ~ out_prob_AD), plot = FALSE, print.auc = FALSE)
npv_ppv_best <- coords(auc_test_AD_10analytes, x = "all", best.method="youden", input = "threshold", ret = c("threshold","acc","npv","ppv","sensitivity","specificity"))
npv_ppv_best[which.min(abs(fixed_threshold-npv_ppv_best$threshold)),]
#        threshold  accuracy       npv       ppv sensitivity specificity
# 4616 0.04865769 0.5927127 0.7854821 0.3054569   0.4886364    0.627597
round(auc_test_AD_10analytes$auc, 2) # 0.57
out_prob_AD <- as.data.frame(out_prob_AD)
sum(is.na(out_prob_AD)) # 0
out_prob_AD$sample_id <- row.names(out_prob_AD)
out_prob_AD <- inner_join(out_prob_AD, AD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_AD, sample_id, status_for_analysis)
table(out_prob_AD$out_prob_AD > fixed_threshold, out_prob_AD$status_for_analysis)
#         ad   CO
# FALSE  991 3625
# TRUE   945 2151

## FTD ##
dim(FTD_expr_imputed) #  6427 5142
FTD_expr_imputed[1:4,1:13]
table(FTD_expr_imputed$level)
#   0    1 
# 6264  163 
table(FTD_expr_imputed$sex)
#   1    2 
# 2734 3693 
sum(is.na(FTD_expr_imputed)) # 0
FTD_expr_subset <- FTD_expr_imputed[, selected_cols]
row.names(FTD_expr_subset) <- FTD_expr_imputed$sample_id
dim(FTD_expr_subset) # 6427   13
# Testing for all analytes in FTD (N=10)
out_prob_FTD <- predict(all_disc, newdata = FTD_expr_subset, type='response')
auc_test_FTD_10analytes <- roc(as.formula(FTD_expr_subset$level ~ out_prob_FTD), plot = FALSE, print.auc = FALSE)
npv_ppv_best <- coords(auc_test_FTD_10analytes, x = "all", best.method="youden", input = "threshold", ret = c("threshold","acc","npv","ppv","sensitivity","specificity"))
npv_ppv_best[which.min(abs(fixed_threshold-npv_ppv_best$threshold)),]
#      threshold  accuracy       npv        ppv sensitivity specificity
# 4012 0.0486824 0.6279757 0.9827973 0.03890728   0.5766871   0.6293103
round(auc_test_FTD_10analytes$auc, 2) # 0.65
out_prob_FTD <- as.data.frame(out_prob_FTD)
sum(is.na(out_prob_FTD)) # 0
out_prob_FTD$sample_id <- row.names(out_prob_FTD)
out_prob_FTD <- inner_join(out_prob_FTD, FTD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_FTD, sample_id, status_for_analysis)
table(out_prob_FTD$out_prob_FTD > fixed_threshold, out_prob_FTD$status_for_analysis)
#         CO  ftd
# FALSE 3942   69
# TRUE  2322   94

## AD
AD_expr_subset <- AD_expr_imputed[, selected_cols]
row.names(AD_expr_subset) <- AD_expr_imputed$sample_id
dim(AD_expr_subset) # 7712   13
# Testing for all analytes in AD (N=10)
out_prob_AD <- predict(all_disc, newdata = AD_expr_subset, type='response')
auc_test_AD_10analytes <- roc(as.formula(AD_expr_subset$level ~ out_prob_AD), plot = FALSE, print.auc = FALSE)
npv_ppv_best <- coords(auc_test_AD_10analytes, x = "all", best.method="youden", input = "threshold", ret = c("threshold","acc","npv","ppv","sensitivity","specificity"))
npv_ppv_best[which.min(abs(fixed_threshold-npv_ppv_best$threshold)),]
#        threshold  accuracy       npv       ppv sensitivity specificity
# 4616 0.04865769 0.5927127 0.7854821 0.3054569   0.4886364    0.627597
round(auc_test_AD_10analytes$auc, 2) # 0.57
out_prob_AD <- as.data.frame(out_prob_AD)
sum(is.na(out_prob_AD)) # 0
out_prob_AD$sample_id <- row.names(out_prob_AD)
out_prob_AD <- inner_join(out_prob_AD, AD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_AD, sample_id, status_for_analysis)
table(out_prob_AD$out_prob_AD > fixed_threshold, out_prob_AD$status_for_analysis)
#         ad   CO
# FALSE  991 3625
# TRUE   945 2151

colnames(out_prob_AD)[1] <- "out_prob_PD"
colnames(out_prob_FTD)[1] <- "out_prob_PD"
PD_Probablities <- rbind(out_prob_PD, out_prob_AD, out_prob_FTD)
dim(PD_Probablities) # 20440    3
names(PD_Probablities)[3] <- "Final_Status"
length(unique(PD_Probablities$sample_id)) # 8888

PD_Probablities <- PD_Probablities[order(PD_Probablities$sample_id, PD_Probablities$out_prob_PD),]
PD_Probablities <- PD_Probablities[!duplicated(PD_Probablities$sample_id),]
PD_Probablities <- na.omit(PD_Probablities)
dim(PD_Probablities) # 8888    3
table(PD_Probablities$Final_Status)
#   ad   CO  ftd   pd 
# 1936 6264  163  525 
save(PD_Probablities, file="Plasma_PD_Probability_N8888_AgeSexPlatePCs.RData")

png("Probability_Violinplot_for_Plasma_PD_Signature_AgeSexPlatePCs.png", units="mm", width=100, height=100, res=1000)
ggplot(PD_Probablities, aes(x=Final_Status, y=out_prob_PD, color=Final_Status)) + geom_violin(trim = TRUE) + 
  theme(legend.position = "none", panel.background = element_blank(), axis.line = element_line(colour = "black")) +
  labs(x ="Clinical Diagnosis", y = "Predicted PD Probability") + ggtitle("Probability Distribution for PD signature")
dev.off()

png("AUC_Lasso_10_IterativeModel_PDvsCO_Testing_Other_NDs_AgeSexPlatePCs.png", units="mm", width=120, height=120, res=1000)
plot_curve <- plot.roc(auc_test_AD_10analytes, xlim=c(1,0), ylim=c(0, 1), lty=1, lwd=6, col=alpha("blue", 0.5), 
                       main = c("PD Model Testing across other dementia"))
text(0.6, 0.15, paste("AD: ",  round(auc_test_AD_10analytes$auc, 2), sep=""), col=("blue"), cex=0.85, pos=4)
plot_curve <- plot.roc(auc_test_FTD_10analytes, xlim=c(1,0), ylim=c(0, 1), lty=1, lwd=6, col=alpha("green4", 0.5), add=TRUE)
text(0.6, 0.10, paste("FTD: ",  round(auc_test_FTD_10analytes$auc, 2), sep=""), col=("green4"), cex=0.80, pos=4)
dev.off()
save(auc_test_AD_10analytes, auc_test_FTD_10analytes, file="AUC_Lasso_10_IterativeModel_PDvsCO_Testing_Other_NDs_AgeSexPlatePCs.RData")

##########################################################################################
##### PD Prediction Model (Iterative Class Balanced Training {70%} and Testing{30%}) #####
##########################################################################################
rm(list = ls())
library(dplyr) # dplyr_1.1.4
library(pROC) # pROC_1.18.5
library(ROCR) # ROCR_1.0-11
library(caret) # caret_6.0-94
library(glmnet) # glmnet_4.1-8
options(stringsAsFactors = FALSE)
options(width = 160)
set.seed(1)

load("~/files/WashU_pilot_data_explore/PredictionModel_03252025/AD/CSF_Complete_DiseaseSpecific_Imputed_Expression_DFs_For_CSFPredModel_AgeSexPlatePCs.RData")
protWeight_df <- read.csv("~/files/WashU_pilot_data_explore/PredictionModel_03252025/PD/ProteinWeights_Sorted_50runs_PDvsCO_AgeSexPlatePCs.csv")
protWeight_df <- protWeight_df[order(protWeight_df$count, decreasing=T),]
protWeight_df$count[1:15]
# [1] 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50

iterativeModel_10analytes <- protWeight_df$Proteins[1:10]
selected_cols <- c("level", "age_at_visit", "sex", iterativeModel_10analytes)
length(selected_cols) # 13

row.names(AD_expr_imputed) <- AD_expr_imputed$sample_id
AD_expr_subset <- AD_expr_imputed[, selected_cols]
dim(AD_expr_subset) # 7712   13

row.names(FTD_expr_imputed) <- FTD_expr_imputed$sample_id
FTD_expr_subset <- FTD_expr_imputed[, selected_cols]
dim(FTD_expr_subset) # 6427   13

row.names(PD_expr_imputed) <- PD_expr_imputed$sample_id
PD_expr_subset <- PD_expr_imputed[, selected_cols]
dim(PD_expr_subset) # 6301   13

# for saving performance info from the model
perf_df_pd <- c()
perf_df_pd_base <- c()
# for saving predictions from the model
iter_pred_pd_base <- list()
iter_pred_pd <- list()
iter_y_pd <- list()

# other dementias
perf_df_DLB <- c()
iter_pred_DLB <- list()
iter_y_DLB <- list()
perf_df_FTD <- c()
iter_pred_FTD <- list()
iter_y_FTD <- list()
perf_df_AD <- c()
iter_pred_AD <- list()
iter_y_AD <- list()

table(PD_expr_subset$level)
#   0    1 
# 5776  525 
min(table(PD_expr_subset$level)) # 525

for (i in 1:100) {
  set.seed(i)
  print(i)  
  # Separate data by levels # In each iteration 70% data is training (equal number of AD and CO) and 30% data is testing (equal number of AD and CO)
  level_0 <- PD_expr_subset[PD_expr_subset$level == 0, ]
  level_1 <- PD_expr_subset[PD_expr_subset$level == 1, ]
  
  N <- min(nrow(level_0), nrow(level_1))
  train_n <- round(0.7 * N)
  test_n <- round(0.3 * N)
  
  train_level_1 <- level_1[sample(1:nrow(level_1), size = min(train_n, nrow(level_1))), ]
  remain_level_1 <- level_1[!rownames(level_1) %in% rownames(train_level_1), ]
  test_level_1 <- remain_level_1[sample(1:nrow(remain_level_1), size = min(test_n, nrow(remain_level_1))), ]
  
  train_level_0 <- level_0[sample(1:nrow(level_0), size = nrow(train_level_1)), ]
  remain_level_0 <- level_0[!rownames(level_0) %in% rownames(train_level_0), ]
  test_level_0 <- remain_level_0[sample(1:nrow(remain_level_0), size = nrow(test_level_1)), ]
  
  # Combine training and testing sets
  Train_set <- rbind(train_level_1, train_level_0)
  Test_set <- rbind(test_level_1, test_level_0)
  print(nrow(Train_set))
  print(nrow(Test_set))
  
  # Baseline model
  bl_disc <- glm(as.formula(level ~ age_at_visit + sex), data = Train_set, family = 'binomial')
  base_pred_val_pd <- predict(bl_disc, newdata = Test_set, type='response')
  ROC_base_val_pd <- roc(as.formula(Test_set$level ~ base_pred_val_pd), plot = FALSE, print.auc = FALSE)
  
  # 10-protein prediction model
  all_disc <<- glm(level ~ ., data = Train_set, family = 'binomial')
  pred_val_glm_pd <- predict(all_disc, newdata = Test_set, type='response')
  ROC_val_pd <- roc(as.formula(Test_set$level ~ pred_val_glm_pd), plot = FALSE, print.auc = FALSE)
  
  iter_pred_pd_base[[i]] <- as.numeric(base_pred_val_pd)
  iter_pred_pd[[i]] <- as.numeric(pred_val_glm_pd)
  iter_y_pd[[i]] <- Test_set$level
  
  # each performance result for iterations, if there is any ties, get the first one only
  perf_df_pd <- rbind(perf_df_pd, 
                      coords(ROC_val_pd, x = "best", best.method="youden", input = "threshold", 
                             ret = c("threshold","auc", "acc","npv","ppv","sensitivity","specificity")) %>% 
                        mutate(iteration = i, auc = ROC_val_pd$auc, .before = 1) %>% 
                        summarise_all(first))
  perf_df_pd_base <- rbind(perf_df_pd_base, 
                           coords(ROC_base_val_pd, x = "best", best.method="youden", input = "threshold", 
                                  ret = c("threshold","auc", "acc","npv","ppv","sensitivity","specificity")) %>% 
                             mutate(iteration = i, auc = ROC_base_val_pd$auc, .before = 1) %>% 
                             summarise_all(first))
  
  # Predictions in other dementia
  # FTD
  pred_val_glm_FTD <- predict(all_disc, newdata = FTD_expr_subset, type='response')
  ROC_val_FTD <- roc(as.formula(FTD_expr_subset$level ~ pred_val_glm_FTD), plot = FALSE, print.auc = FALSE)
  iter_pred_FTD[[i]] <- as.numeric(pred_val_glm_FTD)
  iter_y_FTD[[i]] <- FTD_expr_subset$level
  perf_df_FTD <- rbind(perf_df_FTD, 
                       coords(ROC_val_FTD, x = "best", best.method="youden", input = "threshold", 
                              ret = c("threshold","auc", "acc","npv","ppv","sensitivity","specificity")) %>% 
                         mutate(iteration = i, auc = ROC_val_FTD$auc, .before = 1) %>%  summarise_all(first))
  # AD
  pred_val_glm_AD <- predict(all_disc, newdata = AD_expr_subset, type='response')
  ROC_val_AD <- roc(as.formula(AD_expr_subset$level ~ pred_val_glm_AD), plot = FALSE, print.auc = FALSE)
  iter_pred_AD[[i]] <- as.numeric(pred_val_glm_AD)
  iter_y_AD[[i]] <- AD_expr_subset$level
  perf_df_AD <- rbind(perf_df_AD, 
                      coords(ROC_val_AD, x = "best", best.method="youden", input = "threshold", 
                             ret = c("threshold","auc", "acc","npv","ppv","sensitivity","specificity")) %>% 
                        mutate(iteration = i, auc = ROC_val_AD$auc, .before = 1) %>% summarise_all(first))
}
# Train_set = 736 ( 368 PD, 368 CO)
# Test_set = 314 ( 157 PD, 157 CO)

pred_pd <- prediction(iter_pred_pd, iter_y_pd)
pred_pd_base <- prediction(iter_pred_pd_base, iter_y_pd)
perf_pd <- performance(pred_pd, "tpr", "fpr")
perf_pd_base <- performance(pred_pd_base, "tpr", "fpr")

mean(perf_df_pd$auc) # 0.8344704
sd(perf_df_pd$auc) # 0.02161382
mean(perf_df_pd_base$auc) # 0.5590008
sd(perf_df_pd_base$auc) # 0.02993082
dim(perf_df_pd) # 100   8
dim(perf_df_pd_base) # 100   8

out_prob_PD <- predict(all_disc, newdata = PD_expr_subset, type='response')
out_prob_PD <- as.data.frame(out_prob_PD)
out_prob_PD$sample_id <- row.names(out_prob_PD)
dim(out_prob_PD) # 6301    2
out_prob_PD_mean <- inner_join(out_prob_PD, PD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_PD, sample_id, status_for_analysis)
dim(out_prob_PD_mean) # 6301    3

png("AUC_Lasso_10_IterativeModel_PDvsCO_AgeSexProjectPCs_Iterative_ClassBalanced.png", units="mm", width=120, height=120, res=1000)
plot(perf_pd, avg = "vertical", 
     spread.estimate="stderror", lwd=6, lty=1, plotCI.col = alpha("red", 0.5), col = alpha("red", 0.5), main = c("PDvsCO"))
text(0.4, 0.25, paste("Test Data (30%): ",  round(mean(perf_df_pd$auc), 2), sep=""), col=("red"), cex=0.85, pos=4)
plot(perf_pd_base, avg = "vertical", 
     spread.estimate="stderror", lty=2, add = T, plotCI.col = alpha("red", 0.5), col = alpha("red", 0.5))
text(0.4, 0.20, paste("Baseline: ", round(mean(perf_df_pd_base$auc), 2), sep=""), col=alpha("red", 0.5), cex=0.75, pos=4)
dev.off()
save(perf_pd, perf_df_pd, perf_pd_base, perf_df_pd_base, 
     file="AUC_Lasso_10_IterativeModel_PDvsCO_AgeSexProjectPCs_Iterative_ClassBalanced.RData")

## prediction/performance in other dementia ##
# FTD
pred_FTD <- prediction(iter_pred_FTD, iter_y_FTD)
perf_FTD <- performance(pred_FTD, "tpr", "fpr")
mean(perf_df_FTD$auc) # 0.6193299
sd(perf_df_FTD$auc) # 0.01712578
dim(perf_df_FTD) # 100   8
out_prob_FTD_mean <- sapply(1:length(iter_pred_FTD[[1]]), function(i) {
  mean(sapply(iter_pred_FTD, function(x) x[i]))
})
out_prob_FTD_mean <- as.data.frame(cbind(row.names(FTD_expr_subset), out_prob_FTD_mean))
colnames(out_prob_FTD_mean) <- c("sample_id", "out_prob_FTD")
out_prob_FTD_mean <- inner_join(out_prob_FTD_mean, FTD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_FTD, sample_id, status_for_analysis)

# AD
pred_AD <- prediction(iter_pred_AD, iter_y_AD)
perf_AD <- performance(pred_AD, "tpr", "fpr")
mean(perf_df_AD$auc) # 0.5891634
sd(perf_df_AD$auc) # 0.01050784
dim(perf_df_AD) # 100   8
out_prob_AD_mean <- sapply(1:length(iter_pred_AD[[1]]), function(i) {
  mean(sapply(iter_pred_AD, function(x) x[i]))
})
out_prob_AD_mean <- as.data.frame(cbind(row.names(AD_expr_subset), out_prob_AD_mean))
colnames(out_prob_AD_mean) <- c("sample_id", "out_prob_AD")
out_prob_AD_mean <- inner_join(out_prob_AD_mean, AD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_AD, sample_id, status_for_analysis)

png("AUC_Lasso_10_IterativeModel_PDvsCO_Testing_Other_NDs_AgeSexProjectPCs_Iterative_ClassBalanced.png", units="mm", width=120, height=120, res=1000)
# xlab='Average Specificity', ylab='Average Sensitivity', xlim=c(1,0), 
plot(perf_FTD, avg = "vertical", 
     spread.estimate = "stderror", lwd = 6, 
     plotCI.col = alpha("green4", 0.5), col = alpha("green4", 0.5),
     main = "PD model testing across other dementia")
text(0.4, 0.10, paste("FTD: ", round(mean(perf_df_FTD$auc), 2)), 
     col = "green4", cex = 0.80, pos = 4)

plot(perf_AD, avg = "vertical", 
     spread.estimate = "stderror", lwd = 6, add = TRUE,
     plotCI.col = alpha("blue", 0.5), col = alpha("blue", 0.5))
text(0.4, 0.15, paste("AD: ", round(mean(perf_df_AD$auc), 2)), 
     col = "blue", cex = 0.80, pos = 4)
dev.off()
save(perf_FTD, perf_df_FTD, perf_AD, perf_df_AD, 
     file="AUC_Lasso_10_IterativeModel_PDvsCO_Testing_Other_NDs_AgeSexProjectPCs_Iterative_ClassBalanced.RData")


colnames(out_prob_FTD_mean)[1] <- "out_prob_PD"
colnames(out_prob_AD_mean)[1] <- "out_prob_PD"
PD_Probablities <- rbind(out_prob_PD_mean, out_prob_FTD_mean, out_prob_AD_mean)
dim(PD_Probablities) # 20440    3
names(PD_Probablities)[3] <- "Final_Status"
length(unique(PD_Probablities$sample_id)) # 8888
PD_Probablities <- PD_Probablities[order(PD_Probablities$sample_id, PD_Probablities$out_prob_PD),]
PD_Probablities <- PD_Probablities[!duplicated(PD_Probablities$sample_id),]
PD_Probablities <- na.omit(PD_Probablities)
PD_Probablities$out_prob_PD <- as.numeric(as.character(PD_Probablities$out_prob_PD))
dim(PD_Probablities) # 8888    3
table(PD_Probablities$Final_Status)
# ad   CO  ftd   pd 
# 1936 6264  163  525 
save(PD_Probablities, file="Plasma_PD_Probability_N8888_AgeSexPlatePCs_Iterative_ClassBalanced.RData")
png("Probability_Violinplot_for_Plasma_PD_Signature_AgeSexProjectPCs_Iterative_ClassBalanced.png", units="mm", width=100, height=100, res=1000)
ggplot(PD_Probablities, aes(x=Final_Status, y=out_prob_PD, color=Final_Status)) + geom_violin(trim = TRUE) + 
  theme(legend.position = "none", panel.background = element_blank(), axis.line = element_line(colour = "black")) +
  labs(x ="Clinical Diagnosis", y = "Predicted PD Probability") + ggtitle("Probability Distribution for PD signature")
dev.off()


# Calculate mean and standard deviation for each column
calculate_mean_sd <- function(df) {
  df <- df[,-1]
  # Ensure only numeric columns are included
  numeric_cols <- df[sapply(df, is.numeric)]
  # Calculate mean and SD for each numeric column
  formatted_stats <- sapply(numeric_cols, function(x) {
    mean_val <- mean(x, na.rm = TRUE)
    sd_val <- sd(x, na.rm = TRUE)
    paste0(round(mean_val, 2), "(", round(sd_val, 2), ")")
  })
  # Convert to a single-row dataframe
  result <- as.data.frame(t(formatted_stats))
  colnames(result) <- colnames(numeric_cols)
  return(result)
}
calculate_mean_sd(perf_df_pd)
#          auc  threshold   accuracy        npv        ppv sensitivity specificity
# 1 0.83(0.02) 0.48(0.08) 0.77(0.02) 0.77(0.04) 0.78(0.05)  0.76(0.07)  0.78(0.08)
calculate_mean_sd(perf_df_pd_base)
#          auc  threshold   accuracy       npv        ppv sensitivity specificity
#1 0.56(0.03) 0.49(0.04) 0.56(0.02) 0.6(0.07) 0.57(0.04)   0.64(0.2)  0.49(0.21)
calculate_mean_sd(perf_df_FTD)
#          auc  threshold   accuracy     npv     ppv sensitivity specificity
# 1 0.62(0.02) 0.37(0.09) 0.63(0.11) 0.98(0) 0.04(0)  0.56(0.12)  0.63(0.12)
calculate_mean_sd(perf_df_AD)
#          auc  threshold   accuracy        npv        ppv sensitivity specificity
# 1 0.59(0.01) 0.41(0.04) 0.62(0.02) 0.79(0.01) 0.32(0.01)  0.45(0.05)  0.68(0.05)

##################
##### FTD ########
##################

######################################################
#### FTD Lasso Model Generation (Iterative; n=50) ####
######################################################
# subset for FTD novel proteins
table(FTD_expr_imputed$level)
#    0    1
#  6264 163 
table(FTD_expr_imputed$sex)
#   1    2 
# 2734 3693

# subset for FTD novel proteins
selected_cols <- c("sample_id", "level", "age_at_visit", "sex", FTD_FDR_nonCorr)
length(selected_cols) # 1795
FTD_expr_subset <- FTD_expr_imputed[, selected_cols]
dim(FTD_expr_subset) # 6427 1795
row.names(FTD_expr_subset) <- FTD_expr_subset$sample_id

# comparison <- "FTDvsCO"
FTD_expr_subset[1:2,1:5]
#                                                                 sample_id level age_at_visit sex  X8997.4
# 00296564-0852-46f1-9de9-1fcb7b461e44 00296564-0852-46f1-9de9-1fcb7b461e44     0           89   2 1.791033
# 003314ef-0ad6-4baf-8dbb-8a2e5eecb7df 003314ef-0ad6-4baf-8dbb-8a2e5eecb7df     0           49   2 1.346433

iterative_lasso_modeling(FTD_expr_subset, "FTDvsCO")   

# [1] "summary of proteins # selected in each of the 50 iterations:"
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 134.0   152.0   165.5   169.6   184.8   223.0 
# [1] "summary of most important protein selection:"
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 0.000   0.000   0.000   4.735   4.000  50.000 
# [1] "Top 10 proteins selected for final Lasso model:"
# [1] "X18175.65" "X18896.23" "X21909.10" "X22141.59" "X2962.50"  "X4930.21"  "X5109.24"  "X8997.4"   "X15631.18" "X19279.42"


####################################################################################
##### FTD - Training and Testing (70 - 30%) in AD data using 10 analytes model #####
####################################################################################
rm(list = ls())
library(dplyr) # dplyr_1.1.4
library(pROC) # pROC_1.18.5
library(ROCR) # ROCR_1.0-11
library(caret) # caret_6.0-94
library(glmnet) # glmnet_4.1-8
options(stringsAsFactors = FALSE)
options(width = 160)
set.seed(1)

load("~/files/WashU_pilot_data_explore/PredictionModel_03252025/AD/CSF_Complete_DiseaseSpecific_Imputed_Expression_DFs_For_CSFPredModel_AgeSexPlatePCs.RData")
iterative_model_testing(FTD_expr_subset, "FTDvsCO")
# [1] "Overall sample case-control distribution:"
# 
# 0    1 
# 6264  163 
# [1] "Any train sample present in test data:"
# 
# FALSE 
# 4499 
# [1] "Dim of train data:"
# [1] 4499   13
# 
# 0    1 
# 4385  114 
# [1] "Dim of test data:"
# [1] 1928   13
# 
# 0    1 
# 1879   49 
# Setting levels: control = 0, case = 1
# Setting direction: controls < cases
# [1] "Performance for baseline model:"
# threshold accuracy       npv        ppv sensitivity specificity
# threshold 0.0327656  0.80861 0.9864691 0.07446809   0.5714286   0.8147951
# [1] 0.73
# Setting levels: control = 0, case = 1
# Setting direction: controls < cases
# [1] "Performance for 10-analyte model:"
# threshold  accuracy       npv       ppv sensitivity specificity
# threshold 0.02070398 0.7852697 0.9945982 0.0917226   0.8367347   0.7839276
# [1] 0.88
# [1] "Fixed threshold:"
# [1] 0.02070398

# Entire figure at one place
png("AUC_Lasso_10_IterativeModel_FTDvsCO_AgeSexPlatePCs.png", units="mm", width=120, height=120, res=1000)
plot_curve <- plot.roc(auc_test_bl, xlim=c(1,0), ylim=c(0, 1), col=alpha("green4", 0.5), lty=2, 
                       main = "FTDvsCO")
text(0.6, 0.20, paste("Baseline: ", round(auc_test_bl$auc, 2), sep=""), col=alpha("green4", 0.5), cex=0.75, pos=4)
plot_curve <- plot.roc(auc_test_10analytes, xlim=c(1,0), ylim=c(0, 1), lty=1, lwd=6, col=alpha("green4", 0.5), add=TRUE)
text(0.6, 0.25, paste("Test Data (30%): ",  round(auc_test_10analytes$auc, 2), sep=""), col=("green4"), cex=0.85, pos=4)
dev.off()
save(auc_test_bl, auc_test_10analytes, file="AUC_Lasso_10_IterativeModel_FTDvsCO_AgeSexPlatePCs.RData")

# get FTD probablity for all individuals
out_prob_FTD <- predict(all_disc, newdata = rbind(test_discovery, train_discovery), type='response')
out_prob_FTD <- as.data.frame(out_prob_FTD)
out_prob_FTD$sample_id <- row.names(out_prob_FTD)
dim(out_prob_FTD) # 6427    2
out_prob_FTD <- inner_join(out_prob_FTD, FTD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_FTD, sample_id, status_for_analysis)
dim(out_prob_FTD) # 6427    3
table(out_prob_FTD$out_prob_FTD > fixed_threshold, out_prob_FTD$status_for_analysis)
#         CO  ftd
# FALSE 4863   22
# TRUE  1401  141


############################################
##### FTD - Testing in other disorders #####
############################################
## AD
AD_expr_subset <- AD_expr_imputed[, selected_cols]
row.names(AD_expr_subset) <- AD_expr_imputed$sample_id
dim(AD_expr_subset) # 7712   13
# Testing for all analytes in AD (N=10)
out_prob_AD <- predict(all_disc, newdata = AD_expr_subset, type='response')
auc_test_AD_10analytes <- roc(as.formula(AD_expr_subset$level ~ out_prob_AD), plot = FALSE, print.auc = FALSE)
npv_ppv_best <- coords(auc_test_AD_10analytes, x = "all", best.method="youden", input = "threshold", ret = c("threshold","acc","npv","ppv","sensitivity","specificity"))
npv_ppv_best[which.min(abs(fixed_threshold-npv_ppv_best$threshold)),]
#       threshold  accuracy       npv       ppv sensitivity specificity
# 5945 0.02069939 0.6807573 0.7787685 0.3512443   0.3207645   0.8014197
round(auc_test_AD_10analytes$auc, 2) # 0.59
out_prob_AD <- as.data.frame(out_prob_AD)
sum(is.na(out_prob_AD)) # 0
out_prob_AD$sample_id <- row.names(out_prob_AD)
out_prob_AD <- inner_join(out_prob_AD, AD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_AD, sample_id, status_for_analysis)
table(out_prob_AD$out_prob_AD > fixed_threshold, out_prob_AD$status_for_analysis)
#         ad   CO
# FALSE 1315 4630
# TRUE   621 1146


## PD
PD_expr_subset <- PD_expr_imputed[, selected_cols]
row.names(PD_expr_subset) <- PD_expr_imputed$sample_id
dim(PD_expr_subset) # 6301   13
# Testing for all analytes in PD (N=10)
out_prob_PD <- predict(all_disc, newdata = PD_expr_subset, type='response')
auc_test_PD_10analytes <- roc(as.formula(PD_expr_subset$level ~ out_prob_PD), plot = FALSE, print.auc = FALSE)
npv_ppv_best <- coords(auc_test_PD_10analytes, x = "all", best.method="youden", input = "threshold", ret = c("threshold","acc","npv","ppv","sensitivity","specificity"))
npv_ppv_best[which.min(abs(fixed_threshold-npv_ppv_best$threshold)),]
#       threshold  accuracy       npv       ppv sensitivity specificity
# 4951 0.02069939 0.7670211 0.9351515 0.1509993   0.3885714   0.8014197
round(auc_test_PD_10analytes$auc, 2) # 0.66
out_prob_PD <- as.data.frame(out_prob_PD)
sum(is.na(out_prob_PD)) # 0
out_prob_PD$sample_id <- row.names(out_prob_PD)
out_prob_PD <- inner_join(out_prob_PD, PD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_PD, sample_id, status_for_analysis)
table(out_prob_PD$out_prob_PD > fixed_threshold, out_prob_PD$status_for_analysis)
#         CO   pd
# FALSE 4630  321
# TRUE  1146  204

colnames(out_prob_AD)[1] <- "out_prob_FTD"
colnames(out_prob_PD)[1] <- "out_prob_FTD"
FTD_Probablities <- rbind(out_prob_FTD, out_prob_AD, out_prob_PD)
dim(FTD_Probablities) # 20440    3
names(FTD_Probablities)[3] <- "Final_Status"
length(unique(FTD_Probablities$sample_id)) # 8888
FTD_Probablities <- FTD_Probablities[order(FTD_Probablities$sample_id, FTD_Probablities$out_prob_FTD),]
FTD_Probablities <- FTD_Probablities[!duplicated(FTD_Probablities$sample_id),]
FTD_Probablities <- na.omit(FTD_Probablities)
dim(FTD_Probablities) # 8888    3
table(FTD_Probablities$Final_Status)
# ad   CO  ftd   pd 
# 1936 6264  163  525 
save(FTD_Probablities, file="CSF_FTD_Probability_N8888_AgeSexPlatePCs.RData")
png("Probability_Violinplot_for_CSF_FTD_Signature_AgeSexPlatePCs.png", units="mm", width=100, height=100, res=1000)
ggplot(FTD_Probablities, aes(x=Final_Status, y=out_prob_FTD, color=Final_Status)) + geom_violin(trim = TRUE) + 
  theme(legend.position = "none", panel.background = element_blank(), axis.line = element_line(colour = "black")) +
  labs(x ="Clinical Diagnosis", y = "Predicted FTD Probability") + ggtitle("Probability Distribution for FTD signature")
dev.off()

png("AUC_Lasso_10_IterativeModel_FTDvsCO_Testing_Other_NDs_AgeSexPlatePCs.png", units="mm", width=120, height=120, res=1000)
plot_curve <- plot.roc(auc_test_AD_10analytes, xlim=c(1,0), ylim=c(0, 1), lty=1, lwd=6, col=alpha("blue", 0.5), 
                       main = c("FTD Model Testing across other dementia"))
text(0.6, 0.15, paste("AD: ",  round(auc_test_AD_10analytes$auc, 2), sep=""), col=("blue"), cex=0.85, pos=4)
plot_curve <- plot.roc(auc_test_PD_10analytes, xlim=c(1,0), ylim=c(0, 1), lty=1, lwd=6, col=alpha("red", 0.5), add=TRUE)
text(0.6, 0.05, paste("PD: ",  round(auc_test_PD_10analytes$auc, 2), sep=""), col=("red"), cex=0.80, pos=4)
dev.off()
save(auc_test_AD_10analytes, auc_test_PD_10analytes, file="AUC_Lasso_10_IterativeModel_FTDvsCO_Testing_Other_NDs_AgeSexPlatePCs.RData")


###########################################################################################
##### FTD Prediction Model (Iterative Class Balanced Training {70%} and Testing{30%}) #####
###########################################################################################
rm(list = ls())
library(dplyr) # dplyr_1.1.4
library(pROC) # pROC_1.18.5
library(ROCR) # ROCR_1.0-11
library(caret) # caret_6.0-94
library(glmnet) # glmnet_4.1-8
options(stringsAsFactors = FALSE)
options(width = 160)
set.seed(1)

load("~/files/WashU_pilot_data_explore/PredictionModel_03252025/AD/CSF_Complete_DiseaseSpecific_Imputed_Expression_DFs_For_CSFPredModel_AgeSexPlatePCs.RData")

protWeight_df <- read.csv("~/files/WashU_pilot_data_explore/PredictionModel_03252025/FTD/ProteinWeights_Sorted_50runs_FTDvsCO_AgeSexPlatePCs.csv")
protWeight_df <- protWeight_df[order(protWeight_df$count, decreasing=T),]
protWeight_df$count[1:15]
# [1] 50 50 50 50 50 50 50 50 49 49 49 49 48 48 48

iterativeModel_10analytes <- protWeight_df$Proteins[1:10]
selected_cols <- c("level", "age_at_visit", "sex", iterativeModel_10analytes)
length(selected_cols) # 13

row.names(AD_expr_imputed) <- AD_expr_imputed$sample_id
AD_expr_subset <- AD_expr_imputed[, selected_cols]
dim(AD_expr_subset) # 7712   13

row.names(FTD_expr_imputed) <- FTD_expr_imputed$sample_id
FTD_expr_subset <- FTD_expr_imputed[, selected_cols]
dim(FTD_expr_subset) # 6427   13

row.names(PD_expr_imputed) <- PD_expr_imputed$sample_id
PD_expr_subset <- PD_expr_imputed[, selected_cols]
dim(PD_expr_subset) # 6301   13

# for saving performance info from the model
perf_df_ftd <- c()
perf_df_ftd_base <- c()
# for saving predictions from the model
iter_pred_ftd_base <- list()
iter_pred_ftd <- list()
iter_y_ftd <- list()

# other dementias
perf_df_AD <- c()
iter_pred_AD <- list()
iter_y_AD <- list()
perf_df_DLB <- c()
iter_pred_DLB <- list()
iter_y_DLB <- list()
perf_df_PD <- c()
iter_pred_PD <- list()
iter_y_PD <- list()

table(FTD_expr_subset$level)
#   0    1 
# 6264  163 
min(table(FTD_expr_subset$level)) # 163

for (i in 1:100) {
  set.seed(i)
  print(i)
  # Separate data by levels # In each iteration 70% data is training (equal number of AD and CO) and 30% data is testing (equal number of AD and CO)
  level_0 <- FTD_expr_subset[FTD_expr_subset$level == 0, ]
  level_1 <- FTD_expr_subset[FTD_expr_subset$level == 1, ]
  # Sample rows for level = 1
  train_level_1 <- level_1[sample(1:nrow(level_1), size = round(0.7 * min(table(FTD_expr_subset$level)))), ]
  test_level_1 <- level_1[!rownames(level_1) %in% rownames(train_level_1), ] # maker sure test has no train samples
  test_level_1 <- test_level_1[sample(1:nrow(test_level_1), size = round(0.3 * min(table(FTD_expr_subset$level)))), ] # downsample test dataset
  # Sample rows for level = 0
  train_level_0 <- level_0[sample(1:nrow(level_0), size = nrow(train_level_1)), ]
  test_level_0 <- level_0[!rownames(level_0) %in% rownames(train_level_0), ] # maker sure test has no train samples
  test_level_0 <- test_level_0[sample(1:nrow(test_level_0), size = nrow(test_level_1)), ] # downsample test dataset
  # Combine training and testing sets
  Train_set <- rbind(train_level_1, train_level_0)
  Test_set <- rbind(test_level_1, test_level_0)
  print(nrow(Train_set))
  print(nrow(Test_set))
  
  # Baseline model
  bl_disc <- glm(as.formula(level ~ age_at_visit + sex), data = Train_set, family = 'binomial')
  base_pred_val_ftd <- predict(bl_disc, newdata = Test_set, type='response')
  ROC_base_val_ftd <- roc(as.formula(Test_set$level ~ base_pred_val_ftd), plot = FALSE, print.auc = FALSE)
  
  # 10-protein prediction model
  all_disc <<- glm(level ~ ., data = Train_set, family = 'binomial')
  pred_val_glm_ftd <- predict(all_disc, newdata = Test_set, type='response')
  ROC_val_ftd <- roc(as.formula(Test_set$level ~ pred_val_glm_ftd), plot = FALSE, print.auc = FALSE)
  
  iter_pred_ftd_base[[i]] <- as.numeric(base_pred_val_ftd)
  iter_pred_ftd[[i]] <- as.numeric(pred_val_glm_ftd)
  iter_y_ftd[[i]] <- Test_set$level
  
  # each performance result for iterations, if there is any ties, get the first one only
  perf_df_ftd <- rbind(perf_df_ftd, 
                       coords(ROC_val_ftd, x = "best", best.method="youden", input = "threshold", 
                              ret = c("threshold","auc", "acc","npv","ppv","sensitivity","specificity")) %>% 
                         mutate(iteration = i, auc = ROC_val_ftd$auc, .before = 1) %>% 
                         summarise_all(first))
  perf_df_ftd_base <- rbind(perf_df_ftd_base, 
                            coords(ROC_base_val_ftd, x = "best", best.method="youden", input = "threshold", 
                                   ret = c("threshold","auc", "acc","npv","ppv","sensitivity","specificity")) %>% 
                              mutate(iteration = i, auc = ROC_base_val_ftd$auc, .before = 1) %>% 
                              summarise_all(first))
  
  # Predictions in other dementia
  # AD
  pred_val_glm_AD <- predict(all_disc, newdata = AD_expr_subset, type='response')
  ROC_val_AD <- roc(as.formula(AD_expr_subset$level ~ pred_val_glm_AD), plot = FALSE, print.auc = FALSE)
  iter_pred_AD[[i]] <- as.numeric(pred_val_glm_AD)
  iter_y_AD[[i]] <- AD_expr_subset$level
  perf_df_AD <- rbind(perf_df_AD, 
                      coords(ROC_val_AD, x = "best", best.method="youden", input = "threshold", 
                             ret = c("threshold","auc", "acc","npv","ppv","sensitivity","specificity")) %>% 
                        mutate(iteration = i, auc = ROC_val_AD$auc, .before = 1) %>%  summarise_all(first))
  # PD
  pred_val_glm_PD <- predict(all_disc, newdata = PD_expr_subset, type='response')
  ROC_val_PD <- roc(as.formula(PD_expr_subset$level ~ pred_val_glm_PD), plot = FALSE, print.auc = FALSE)
  iter_pred_PD[[i]] <- as.numeric(pred_val_glm_PD)
  iter_y_PD[[i]] <- PD_expr_subset$level
  perf_df_PD <- rbind(perf_df_PD, 
                      coords(ROC_val_PD, x = "best", best.method="youden", input = "threshold", 
                             ret = c("threshold","auc", "acc","npv","ppv","sensitivity","specificity")) %>% 
                        mutate(iteration = i, auc = ROC_val_PD$auc, .before = 1) %>% summarise_all(first))
}
# Train_set = 228 (114 FTD, 114 CO)
# Test_set = 98 (49 FTD, 49 CO)

pred_ftd <- prediction(iter_pred_ftd, iter_y_ftd)
pred_ftd_base <- prediction(iter_pred_ftd_base, iter_y_ftd)
perf_ftd <- performance(pred_ftd, "tpr", "fpr")
perf_ftd_base <- performance(pred_ftd_base, "tpr", "fpr")

mean(perf_df_ftd$auc) # 0.8843815
sd(perf_df_ftd$auc) # 0.03494659
mean(perf_df_ftd_base$auc) # 0.7075281
sd(perf_df_ftd_base$auc) # 0.05194398
dim(perf_df_ftd) # 100   8
dim(perf_df_ftd_base) # 100   8

# get FTD probablity for all individuals
out_prob_FTD <- predict(all_disc, newdata = FTD_expr_subset, type='response')
out_prob_FTD <- as.data.frame(out_prob_FTD)
out_prob_FTD$sample_id <- row.names(out_prob_FTD)
dim(out_prob_FTD) # 6427    2
out_prob_FTD_mean <- inner_join(out_prob_FTD, FTD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_FTD, sample_id, status_for_analysis)
dim(out_prob_FTD_mean) # 6427   3

png("AUC_Lasso_10_IterativeModel_FTDvsCO_AgeSexPlatePCs_Iterative_ClassBalanced.png", units="mm", width=120, height=120, res=1000)
plot(perf_ftd, avg = "vertical", 
     spread.estimate="stderror", lwd=6, lty=1, plotCI.col = alpha("green4", 0.5), col = alpha("green4", 0.5), main = c("FTDvsCO"))
text(0.4, 0.25, paste("Test Data (30%): ",  round(mean(perf_df_ftd$auc), 2), sep=""), col=("green4"), cex=0.85, pos=4)
plot(perf_ftd_base, avg = "vertical", 
     spread.estimate="stderror", lty=2, add = T, plotCI.col = alpha("green4", 0.5), col = alpha("green4", 0.5))
text(0.4, 0.20, paste("Baseline: ", round(mean(perf_df_ftd_base$auc), 2), sep=""), col=alpha("green4", 0.5), cex=0.75, pos=4)
dev.off()
save(perf_ftd, perf_df_ftd, perf_ftd_base, perf_df_ftd_base,
     file="AUC_Lasso_10_IterativeModel_FTDvsCO_AgeSexPlatePCs_Iterative_ClassBalanced.RData")

############## prediction/performance in other dementias ############## 
# AD
pred_AD <- prediction(iter_pred_AD, iter_y_AD)
perf_AD <- performance(pred_AD, "tpr", "fpr")
mean(perf_df_AD$auc) # 0.5578715
sd(perf_df_AD$auc) # 0.01972435
dim(perf_df_AD) # 100   8
out_prob_AD_mean <- sapply(1:length(iter_pred_AD[[1]]), function(i) {
  mean(sapply(iter_pred_AD, function(x) x[i]))
})
out_prob_AD_mean <- as.data.frame(cbind(row.names(AD_expr_subset), out_prob_AD_mean))
colnames(out_prob_AD_mean) <- c("sample_id", "out_prob_AD")
out_prob_AD_mean <- inner_join(out_prob_AD_mean, AD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_AD, sample_id, status_for_analysis)
dim(out_prob_AD_mean)

# PD
pred_PD <- prediction(iter_pred_PD, iter_y_PD)
perf_PD <- performance(pred_PD, "tpr", "fpr")
mean(perf_df_PD$auc) # 0.6532571
sd(perf_df_PD$auc) # 0.007372025
dim(perf_df_PD) # 100   8
out_prob_PD_mean <- sapply(1:length(iter_pred_PD[[1]]), function(i) {
  mean(sapply(iter_pred_PD, function(x) x[i]))
})
out_prob_PD_mean <- as.data.frame(cbind(row.names(PD_expr_subset), out_prob_PD_mean))
colnames(out_prob_PD_mean) <- c("sample_id", "out_prob_PD")
out_prob_PD_mean <- inner_join(out_prob_PD_mean, PD_expr_imputed[,1:5], by="sample_id") %>% select(out_prob_PD, sample_id, status_for_analysis)
dim(out_prob_PD_mean)


png("AUC_Lasso_10_IterativeModel_FTDvsCO_Testing_Other_NDs_AgeSexPlatePCs_Iterative_ClassBalanced.png", units="mm", width=120, height=120, res=1000)
# xlab='Average Specificity', ylab='Average Sensitivity', xlim=c(1,0), 
plot(perf_AD, avg = "vertical",  
     spread.estimate="stderror", lwd=6, plotCI.col = alpha("blue", 0.5), col = alpha("blue", 0.5), main = c("FTD model testing across other dementia"))
text(0.4, 0.15, paste("AD: ",  round(mean(perf_df_AD$auc), 2), sep=""), col=("blue"), cex=0.85, pos=4)
plot(perf_PD, avg = "vertical", spread.estimate="stderror", lwd=6, add = T, plotCI.col = alpha("red", 0.5), col = alpha("red", 0.5))
text(0.4, 0.05, paste("PD: ",  round(mean(perf_df_PD$auc), 2), sep=""), col=("red"), cex=0.80, pos=4)
dev.off()
save(perf_AD, perf_df_AD, perf_PD, perf_df_PD, 
     file="AUC_Lasso_10_IterativeModel_FTDvsCO_Testing_Other_NDs_AgeSexPlatePCs_Iterative_ClassBalanced.RData")

colnames(out_prob_AD_mean)[1] <- "out_prob_FTD"
colnames(out_prob_PD_mean)[1] <- "out_prob_FTD"
FTD_Probablities <- rbind(out_prob_FTD_mean, out_prob_AD_mean, out_prob_PD_mean)
dim(FTD_Probablities) # 20440     3
names(FTD_Probablities)[3] <- "Final_Status"
length(unique(FTD_Probablities$sample_id)) # 8888
FTD_Probablities <- FTD_Probablities[order(FTD_Probablities$sample_id, FTD_Probablities$out_prob_FTD),]
FTD_Probablities <- FTD_Probablities[!duplicated(FTD_Probablities$sample_id),]
FTD_Probablities <- na.omit(FTD_Probablities)
FTD_Probablities$out_prob_FTD <- as.numeric(as.character(FTD_Probablities$out_prob_FTD))
dim(FTD_Probablities) # 8888    3
table(FTD_Probablities$Final_Status)
# ad   CO  ftd   pd 
# 1936 6264  163  525 
save(FTD_Probablities, file="CSF_FTD_Probability_N8888_AgeSexPlatePCs_Iterative_ClassBalanced.RData")
png("Probability_Violinplot_for_CSF_FTD_Signature_AgeSexPlatePCs_Iterative_ClassBalanced.png", units="mm", width=100, height=100, res=1000)
ggplot(FTD_Probablities, aes(x=Final_Status, y=out_prob_FTD, color=Final_Status)) + geom_violin(trim = TRUE) + 
  theme(legend.position = "none", panel.background = element_blank(), axis.line = element_line(colour = "black")) +
  labs(x ="Clinical Diagnosis", y = "Predicted FTD Probability") + ggtitle("Probability Distribution for FTD signature")
dev.off()

# Calculate mean and standard deviation for each column
calculate_mean_sd <- function(df) {
  df <- df[,-1]
  # Ensure only numeric columns are included
  numeric_cols <- df[sapply(df, is.numeric)]
  # Calculate mean and SD for each numeric column
  formatted_stats <- sapply(numeric_cols, function(x) {
    mean_val <- mean(x, na.rm = TRUE)
    sd_val <- sd(x, na.rm = TRUE)
    paste0(round(mean_val, 2), "(", round(sd_val, 2), ")")
  })
  # Convert to a single-row dataframe
  result <- as.data.frame(t(formatted_stats))
  colnames(result) <- colnames(numeric_cols)
  return(result)
}
calculate_mean_sd(perf_df_ftd)
#         auc  threshold   accuracy        npv        ppv sensitivity specificity
# 1 0.88(0.03) 0.45(0.13) 0.83(0.04) 0.86(0.06) 0.81(0.06)  0.86(0.07)   0.8(0.08)
calculate_mean_sd(perf_df_ftd_base)
#         auc  threshold   accuracy        npv        ppv sensitivity specificity
# 1 0.71(0.05) 0.46(0.08) 0.69(0.04) 0.73(0.08) 0.68(0.06)  0.74(0.14)  0.63(0.15)
calculate_mean_sd(perf_df_PD)
#         auc  threshold   accuracy     npv        ppv sensitivity specificity
# 1 0.65(0.01) 0.14(0.06) 0.55(0.06) 0.95(0) 0.12(0.01)   0.7(0.07)  0.53(0.07)
calculate_mean_sd(perf_df_AD)
#           auc threshold   accuracy        npv        ppv sensitivity specificity
# 1 0.56(0.02) 0.29(0.1) 0.62(0.04) 0.78(0.01) 0.31(0.02)  0.41(0.07)  0.69(0.07)