rpotter6298
96 lines
3.8 KiB
R
Executable File
96 lines
3.8 KiB
R
Executable File
loocv_validation <- function(df, logistic_model) {
|
|
n_samples <- nrow(df)
|
|
predicted_probabilities <- numeric(n_samples)
|
|
|
|
for (i in 1:n_samples) {
|
|
test_set <- df[i,]
|
|
test_prob <- predict(logistic_model, newdata = test_set[-1], type = "response")
|
|
predicted_probabilities[i] <- test_prob
|
|
}
|
|
|
|
return(predicted_probabilities)
|
|
|
|
}
|
|
|
|
roc_curve <- function(df, report, P.Val = 0.05, stepwise = FALSE, validation = "none") {
|
|
lrep_sigs <- report[report$adj.P.Val < P.Val,]
|
|
analytes <- row.names(lrep_sigs)
|
|
df_selected <- df[, c("group", analytes)]
|
|
df_selected$group = as.numeric(df_selected$group)
|
|
logistic_model <- glm(group ~ ., data = df_selected, family = "binomial")
|
|
|
|
if (stepwise) {
|
|
logistic_model <- step(logistic_model, direction = "backward")
|
|
}
|
|
|
|
model_summary <- summary(logistic_model)
|
|
accuracy <- NULL
|
|
if (validation == "LOOCV") {
|
|
predicted_probabilities <- loocv_validation(df, logistic_model)
|
|
true_labels <- as.numeric(df$group)
|
|
threshold <- 0.5
|
|
predicted_labels <- ifelse(predicted_probabilities >= threshold, 1, 0)
|
|
correct_predictions <- predicted_labels == true_labels
|
|
accuracy <- mean(correct_predictions)
|
|
roc_obj <- roc(true_labels, predicted_probabilities)
|
|
} else {
|
|
predicted_probabilities <- predict(logistic_model, type = "response")
|
|
roc_obj <- roc(df_selected$group, predicted_probabilities)
|
|
}
|
|
|
|
auc <- auc(roc_obj)
|
|
plot(roc_obj, main = paste("ROC Curve (AUC =", round(auc, 3), ")"),asp=1)
|
|
roc_plot <- recordPlot()
|
|
|
|
logistic_regression <- list(model_summary = model_summary, logistic_model = logistic_model)
|
|
validation_results <- list(validation_method = validation, predicted_probabilities = predicted_probabilities, accuracy = accuracy)
|
|
roc_results <- list(roc_obj = roc_obj, auc = auc, plot = roc_plot)
|
|
|
|
if (validation == "none") {
|
|
results <- list(logistic_regression = logistic_regression, roc = roc_results)
|
|
} else {
|
|
results <- list(logistic_regression = logistic_regression, validation = validation_results, roc = roc_results)
|
|
}
|
|
|
|
return(results)
|
|
}
|
|
|
|
fancy_roc <- function(df, report, P.Val = 0.05, validation = "none") {
|
|
library(pROC)
|
|
lrep_sigs <- report[report$adj.P.Val < P.Val,]
|
|
analytes <- row.names(lrep_sigs)
|
|
df_selected <- df[, c("group", analytes)]
|
|
df_selected$group = as.numeric(df_selected$group)
|
|
logistic_model <- glm(group ~ ., data = df_selected, family = "binomial")
|
|
|
|
logistic_model_stepwise <- NULL
|
|
auc_stepwise <- NULL
|
|
predicted_probabilities_stepwise <- NULL
|
|
if (length(analytes) > 1) {
|
|
logistic_model_stepwise <- step(logistic_model, direction = "backward")
|
|
predicted_probabilities_stepwise <- predict(logistic_model_stepwise, type = "response")
|
|
roc_obj_stepwise <- roc(df_selected$group, predicted_probabilities_stepwise)
|
|
auc_stepwise <- auc(roc_obj_stepwise)
|
|
}
|
|
|
|
predicted_probabilities <- predict(logistic_model, type = "response")
|
|
roc_obj <- roc(df_selected$group, predicted_probabilities)
|
|
auc <- auc(roc_obj)
|
|
|
|
par(cex.axis = 1.5)
|
|
plot(roc_obj, col="red", main = "", xlim=c(1,0), ylim=c(0,1))
|
|
|
|
if (!is.null(logistic_model_stepwise)) {
|
|
predicted_probabilities_stepwise <- predict(logistic_model_stepwise, type = "response")
|
|
roc_obj_stepwise <- roc(df_selected$group, predicted_probabilities_stepwise)
|
|
auc_stepwise <- auc(roc_obj_stepwise)
|
|
lines(roc_obj_stepwise, col="blue")
|
|
legend("bottomright", legend=c("Logistic Regression", "Logistic Regression (Backwards Step)"), col=c("red", "blue"), lty=1, cex=0.8)
|
|
}
|
|
else {
|
|
legend("bottomright", legend="Logistic Regression", col="red", lty=1, cex=0.8)
|
|
}
|
|
text(x=0.4, y=0.4, labels=paste("AUC: ", paste0(round(auc, 3))), pos=1, cex=1, col="red")
|
|
text(x=0.4, y=0.35, labels=paste("AUC: ", paste0(round(auc_stepwise, 3))), pos=1, cex=1, col="blue")}
|
|
|