R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

# Load required packages
library(randomForestSRC)
library(survival)
library(dplyr)
library(parallel)

# Assuming your data frame is named 'data' with:
# - time: survival time
# - status: event indicator (1=event, 0=censored)
# - other columns: predictor variables

# Example data preparation (replace with your actual data)
# data <- your_dataframe
# Make sure it has columns: time, status, and predictor variables

Random Grid Generation (200 Combinations) The code samples 200 random hyperparameter combinations based on the ranges specified in your Eye-Risk project document

set.seed(123)

# Generate 200 random hyperparameter combinations
n_combinations <- 200

hyper_grid <- data.frame(
  ntree    = sample(800:1500, n_combinations, replace = TRUE),
  mtry     = sample(1:10, n_combinations, replace = TRUE),
  nodesize = sample(1:10, n_combinations, replace = TRUE),
  nsplit   = sample(5:15, n_combinations, replace = TRUE)
) %>%
  distinct()

# Display default values mentioned in document (ntree=1000, mtry=6, nodesize=3, nsplit=10)
cat("Default values: ntree=1000, mtry=6, nodesize=3, nsplit=10\n")
cat("Number of unique combinations:", nrow(hyper_grid), "\n")

5-Fold Cross-Validation Stratified on Events This stratifies the folds to ensure balanced event distribution across folds

# Create stratified folds based on event status
set.seed(456)

# Stratify on number of events
create_stratified_folds <- function(data, k = 5) {
  events <- which(data$status == 1)
  censored <- which(data$status == 0)
  
  # Create folds for events and censored separately
  event_folds <- split(sample(events), 1:k)
  censored_folds <- split(sample(censored), 1:k)
  
  # Combine to create balanced folds
  folds <- lapply(1:k, function(i) {
    c(event_folds[[i]], censored_folds[[i]])
  })
  
  return(folds)
}

folds <- create_stratified_folds(data, k = 5)

# Verify stratification
cat("Events per fold:\n")
sapply(folds, function(fold) sum(data$status[fold]))

Cross-Validation Function with C-Index The performance metric is Harrell’s C-index, calculated as 1 - err.rate in randomForestSRC

# Function to evaluate one hyperparameter combination
cv_evaluate <- function(params, data, folds) {
  
  c_indices <- numeric(length(folds))
  
  for (i in seq_along(folds)) {
    
    # Split data
    test_idx <- folds[[i]]
    train_data <- data[-test_idx, ]
    test_data <- data[test_idx, ]
    
    # Fit Random Survival Forest
    rsf_model <- rfsrc(
      Surv(time, status) ~ .,
      data = train_data,
      ntree = params$ntree,
      mtry = params$mtry,
      nodesize = params$nodesize,
      nsplit = params$nsplit,
      splitrule = "logrank",
      importance = "none",
      seed = -1
    )
    
    # Predict on test set
    rsf_pred <- predict(rsf_model, newdata = test_data)
    
    # Calculate C-index (1 - error rate)
    c_indices[i] <- 1 - rsf_pred$err.rate[rsf_pred$ntree]
  }
  
  # Return mean C-index across folds
  return(mean(c_indices, na.rm = TRUE))
}

Grid Search with Parallel Processing Using parallel processing as mentioned in your document (parallelMap):

# Setup parallel processing
n_cores <- detectCores() - 1
cl <- makeCluster(n_cores)
clusterEvalQ(cl, {
  library(randomForestSRC)
  library(survival)
})
clusterExport(cl, c("data", "folds", "cv_evaluate"))

# Perform grid search
cat("Starting grid search with", nrow(hyper_grid), "combinations...\n")

cv_results <- parLapply(cl, 1:nrow(hyper_grid), function(i) {
  c_index <- cv_evaluate(hyper_grid[i, ], data, folds)
  return(data.frame(
    combo_id = i,
    ntree = hyper_grid$ntree[i],
    mtry = hyper_grid$mtry[i],
    nodesize = hyper_grid$nodesize[i],
    nsplit = hyper_grid$nsplit[i],
    mean_c_index = c_index
  ))
})

stopCluster(cl)

# Combine results
results_df <- do.call(rbind, cv_results) %>%
  arrange(desc(mean_c_index))

# Display top 10 combinations
head(results_df, 10)

Final Model with Best Hyperparameters Train the final model using the best hyperparameters on the full dataset:

# Extract best hyperparameters
best_params <- results_df[1, ]

cat("\nBest hyperparameters:\n")
print(best_params)

# Train final model on full dataset
final_model <- rfsrc(
  Surv(time, status) ~ .,
  data = data,
  ntree = best_params$ntree,
  mtry = best_params$mtry,
  nodesize = best_params$nodesize,
  nsplit = best_params$nsplit,
  splitrule = "logrank",
  importance = "permute",  # Get variable importance
  seed = 123
)

# Model summary
print(final_model)

# Variable importance
var_imp <- data.frame(
  variable = names(final_model$importance),
  importance = final_model$importance
) %>%
  arrange(desc(importance))

print(var_imp)

# Final C-index on training data
final_c_index <- 1 - final_model$err.rate[final_model$ntree]
cat("\nFinal C-index:", round(final_c_index, 4), "\n")

Prediction on New Data To make predictions on new observations

# Example: predict on new data
# new_data <- your_new_dataframe

# predictions <- predict(final_model, newdata = new_data)

# Extract predicted mortality (cumulative hazard)
# predicted_risk <- predictions$predicted

# Extract survival probabilities at specific times
# survival_prob <- predictions$survival