MAPAcor_demo.Rmd

---
title: "MAPA Functions"
author: "Ishaan Gadiyar"
date: "`r Sys.Date()`"
output:
  html_document: 
    code_folding: hide
    toc: true
    toc_float: true
    number_sections: true
    theme: spacelab
    highlight: haddock
---

<!-- MJ: General feedback: 
Add text context about the results and the format of the results. 
Would show the sample split level first before the summary.
Check analysis conclusions - looks like random forest also performs better with struct?
Make sure the interpretation of ratio is clear - it’s on the log scale of the effect size (f^2) -->

# Introduction

This tutorial illustrates how to download and use the *MAPA* package through analyses of the *Reproducible Brain Charts (RBC)* dataset (https://reprobrainchart.github.io/) to train and evaluate random forest and ridge regression models.
MAPA stands for maximum achievable predictive accuracy and represents the best predictive accuracy that can be achieved by for a given model class, outcome, and input features.
Pearson's correlation has slower than root-n convergence when estimating achievable prediction accuracy for machine learning prediction models.
This causes bias and failing confidence intervals.
The MAPA estimator is based on semiparametric theory and we showed in our preprint that the the estimator is less biased than the classical estimate of Pearson's correlation ([Jones et al., 2025](https://www.biorxiv.org/content/10.1101/2025.11.26.690778v1)).
*MAPA* can be used across a variety of black-box models that predict continuous outcomes if the convergence criteria are appropriately satisfied.
We ensure the convergence criteria by hyperparameter tuning with nested cross-validation.
Similar analyses and details of the methods are available in our preprint ([Jones et al., 2025](https://www.biorxiv.org/content/10.1101/2025.11.26.690778v1)).

# Loading Packages/Libraries
```{r, message = FALSE}
library(tidyverse)
library(knitr) 
library(glmnet)
library(locfit)
library(ggplot2)
library(ranger)
library(tuneRanger)
library(mlr)    

# Installing and loading MAPA package
if (!"MAPA" %in% rownames(installed.packages())) {
  devtools::install_github("statimagcoll/MAPA")
}
library(MAPA)

# Setting the seed
set.seed(123)
```


# Installing Data 

We analyze the Philadelphia Neurodevelopmental Cohort (PNC) from the *Reproducible Brain Charts (RBC)* dataset, accessible via the *datalad* package in R (https://reprobrainchart.github.io/docs/get_data).

In *pnc_data*, the pertinent data that we use are: 

- `age`: a continuous measure in years that represents how old a subject’s brain is. We predict this measure in both model examples. 
- `17Networks`: columns that contain the string '17Networks' are brain volumes within the Yeo 17 Network brain parcellation.
- _to_: columns that contain '_to_' represent functional connectivity measures, or variables that give a measure of correlation of brain activity between specific regions of the brain.  
- `sex_bin`: the sex of the subject, where 1 codes for male. 
- `euler`: the Euler number, which is used to assess the quality of structural MRI scans. 
- `meanFD`: the average mean frame displacement across the resting state fMRI timeseries. 

The imaging data are not site corrected (unharmonized) because all data for the PNC was collected at a single site.
 
```{r}
# Loading in the structural data - denoted by column names with '17Networks'
rbc <- read_csv("/media/disk2/RBC_version0.1/regional_GMV_all_sites_unharmonized.csv")

# Loading in the functional data - denoted by column names with '_to_'
func <- read_csv("/media/disk2/RBC_version0.1/FC_network17_unharmonized.csv")

# Setting each dataset to only include patients that overlap 
rbc = rbc[rbc$participant_id %in% func$participant_id,]
func = func[func$participant_id %in% rbc$participant_id,]

# Checking that only overlapping participants exist 
sum(rbc$participant_id == func$participant_id) == nrow(func)

# Merging datasets
cols_to_add <- setdiff(names(func), names(rbc))
total_data <- cbind(rbc, func[cols_to_add])

# Subsetting study site to PNC
pnc_data = total_data %>% filter(study_site == "PNC1")

# Turning the sex column into a binary stored in 'sex_bin'
pnc_data <- pnc_data %>%
  mutate(sex_bin = ifelse(pnc_data$sex == "Male", 1, 0))
```

# Ridge Regression Model Examples

This section shows how the *MAPA* package can be used to evaluate ridge regression models with single five-fold cross-fitting and with 20 repeated iterations ("sample splits") of five-fold cross-fitting.
The second approach takes longer, but ensures the results are not dependent on a single permutation of the fold grouping.

## Ridge Regression with Single Cross-Fitting

### Getting Predicted Values

Here, we investigate how the inclusion of either structural or functional MRI data in a ridge regression model impacts model accuracy when predicting age.
Five-fold cross-fitting is used to obtain training predictions for the left-out fold.
The tuning parameter for ridge regression model is trained within the four folds nested 10-fold cross-validation and the minimum lambda value is selected for prediction.
The model out-of-sample predictions are stored as a separate column in the original data frame.
We train two models, one using the structural data and one using the functional connectivity data.

```{r}
# Cleaning data - drop NA before cross-fold validation
data_clean = pnc_data %>% drop_na(age) 

# Instantiating number of folds and assigning fold number randomly 
folds = 5
data_clean$fold <- sample(rep(1:folds, length.out = nrow(data_clean)))

# Sanity check for fold distribution 
table(data_clean$fold)

# Creating storage variables for lambda (a finetuning variable)
lambda_struct = c()
lambda_func = c()

# Creating storage columns for model predictions 
data_clean$preds_struct = NA
data_clean$preds_func = NA

# Creating a for loop for 5-fold cross validation
for(i in 1:folds){
  # Assigning data
  test = data_clean %>% filter(fold == i)
  train = data_clean %>% filter(fold != i)
  
  # Creating model data variables for training  
  x_var_struct = as.matrix(train[, grepl('17Networks', colnames(train)) | colnames(train) %in% c("sex_bin", "euler", "meanFD")])
  y_var = as.matrix(train[, colnames(train) %in% c("age")])
  x_var_func = as.matrix(train[, grepl('_to_', colnames(train)) | colnames(train) %in% c("sex_bin", "euler", "meanFD")])

  # Creating model variables for prediction
  newx_struct = as.matrix(test[, grepl('17Networks', colnames(test)) | colnames(test) %in% c("sex_bin", "euler", "meanFD")])
  newx_func = as.matrix(test[, grepl('_to_', colnames(test)) | colnames(test) %in% c("sex_bin", "euler", "meanFD")])
  newy = as.matrix(test[, colnames(test) %in% c("age")])
  
  # Fitting a ridge regression model
  lambda_seq <- c(
    seq(10, 1, -1),        # From 10 to 1 (step of 1)
    seq(0.9, 0.1, -0.1),   # From 0.9 to 0.1 (step of 0.1)
    seq(0.09, 0.01, -0.01),# From 0.09 to 0.01 (step of 0.01)
    seq(0.009, 0.001, -0.001), # From 0.009 to 0.001 (step of 0.001)
    seq(0.0009, 0.0001, -0.0001), # From 0.0009 to 0.0001 (step of 0.0001)
    seq(0.00009, 0.00001, -0.00001), # From 0.00009 to 0.00001
    seq(0.000009, 0.000001, -0.000001), # From 0.000009 to 0.000001
    seq(0.0000009, 0.0000001, -0.0000001), # From 0.0000009 to 0.0000001
    0
  )
  
  # Training respective models
  ridgeModel_struct <- cv.glmnet(x_var_struct, y_var, alpha = 0, lambda = lambda_seq)
  ridgeModel_func <- cv.glmnet(x_var_func, y_var, alpha = 0, lambda = lambda_seq)
  
  # Getting predictions
  test$preds_struct = predict(ridgeModel_struct, newx = newx_struct, s = "lambda.min")[,]
  test$preds_func = predict(ridgeModel_func, newx = newx_func, s = "lambda.min")[,]
  
  # Storing prediction values
  ## Isolating to a singular data frame with only the prediction values and IDs
  test_collapsed <- test %>%
  group_by(participant_id) %>%
  summarise(
    preds_struct = first(na.omit(preds_struct)),
    preds_func = first(na.omit(preds_func)),
    .groups = "drop"
  )
  
  ## Matching the predictions by ID to the original data frame
  data_clean <- data_clean %>%
  left_join(test_collapsed, by = "participant_id", suffix = c("", "_new"))%>%
  mutate(
    preds_struct = ifelse(is.na(preds_struct_new), preds_struct, preds_struct_new),
    preds_func = ifelse(is.na(preds_func_new), preds_func, preds_func_new)
  ) %>%
  select(-ends_with("_new"))
  
  
  # Storing minimum lambda used for each ridge regression model
  lambda_struct = c(lambda_struct, ridgeModel_struct$lambda.min)
  lambda_func = c(lambda_func, ridgeModel_func$lambda.min)
}
```





### Applying `mapa` to Data Predictions


The `mapa` function is used to obtain the Pearson correlation MAPA.
The function takes in the true outcome values, the predictions from the first model (here, the structural predicted values for age), the fold IDs, and (optionally) the predictions from the second model (here, the functional predicted values).

When two sets of predicted values are given, the `mapa` function computes confidence intervals for both estimators as well as a confidence interval for the difference.
The confidence interval is constructed on the logit squared scale and back-transformed to Pearson's correlation interval for predicted values, [0,1].

<!-- SNV: Ishaan, describe the output structure and interpretation -->
The output returns XX.

```{r}
# Gathering the estimates with one sample split
results = mapa(as.matrix(data_clean$age), as.matrix(data_clean$preds_struct), as.matrix(data_clean$fold), as.matrix(data_clean$preds_func))

# Here, yhat1 represents the prediction accuracy of the structural model and 
# yhat2 represents the prediction accuracy of the functional model 
results$est <- results$est %>% mutate(metric = recode(metric, yhat1 = "Structural Data", yhat2 = "Functional Data"))
results$est_summary <- results$est_summary %>% mutate(metric = recode(metric, yhat1 = "Structural Data", yhat2 = "Functional Data"))

# Previewing summarized results
results$est_summary
head(results$est)
```

### Graphing Outputs

The results from the `mapa` function can be used to plot estimates and 95% confidence intervals for the models and their difference.
Here, we see that the ridge regression model trained with structural MRI data has better age prediction than the model trained with functional MRI data when predicting the continuous variable age.

```{r}
# Subsetting results to only estimates
toPlot = results$est_summary %>% filter((metric == "Structural Data" | metric == "Functional Data") & scale == "r")
toPlot_Diff = results$est_summary %>% filter(metric == "Difference")

# Plotting one-step results 
ggplot(toPlot, aes(x = metric, y = est)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = LB, ymax = UB),
                width = 0.2,
                position = position_dodge(width = 0.5)) +
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust=1)) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "black")  + 
  labs(x = "Ridge Regression Model Data", y = "One-Step Estimator Value") 

# Plotting one-step difference results 
ggplot(toPlot_Diff, aes(x = metric, y = est)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = LB, ymax = UB),
                width = 0.2,
                position = position_dodge(width = 0.5)) +
  theme_bw() + 
  theme(axis.ticks.y = element_blank(), axis.text.y = element_blank()) + 
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "black")  + 
  ylim(-1,1) + 
  labs(x = "One-Step Difference", y = "Difference Magnitude")
```



## Ridge Regression using Sample Splits

### Getting Predicted Values

Here, like section 4.1, we test how the inclusion of either structural or functional MRI data in a ridge regression model impacts model accuracy when predicting age. However, this time we use 20 repeated sample splits to ensure stability across random configurations of the five-fold cross-fitting.
The procedure is similar, except now 20 sets of out-of-sample predictions are generated one for each configuration of the five-fold cross-fitting.

The resulting value is the estimate and confidence interval obtained from the median estimate and standard error across all the sample splits.
```{r}
# Cleaning data
data_clean = pnc_data %>% drop_na(age)

# Assigning sample split number and fold number 
ss_num = 20
folds = 5

# Creating storage variables for hypertuning
lambda_struct = c()
lambda_func = c()

# Assigning storage matrices within data frame columns 
data_clean$struct_ss  <- matrix(nrow = nrow(data_clean), ncol = 0)
data_clean$func_ss  <- matrix(nrow = nrow(data_clean), ncol = 0)
data_clean$fold_ss  <- matrix(nrow = nrow(data_clean), ncol = 0)


for(j in 1:ss_num){
  # Creating a temporary data frame
  data_temp <- data_clean
  
  # Creating storage columns
  data_temp$preds_struct = NA
  data_temp$preds_func = NA
  
  # Randomly assigning fold numbers within the temporary data frame 
  data_temp$fold <- sample(rep(1:folds, length.out = nrow(data_temp)))
  
  # Creating a for loop for 5-fold cross validation
  for(i in 1:folds){
    # Assigning data
    test = data_temp %>% filter(fold == i)
    train = data_temp %>% filter(fold != i)
    
    # Creating model data variables for training  
    x_var_struct = as.matrix(train[, grepl('17Networks', colnames(train)) | colnames(train) %in% c("sex_bin", "euler", "meanFD")])
    y_var = as.matrix(train[, colnames(train) %in% c("age")])
    x_var_func = as.matrix(train[, grepl('_to_', colnames(train)) | colnames(train) %in% c("sex_bin", "euler", "meanFD")])
  
    # Creating model variables for prediction
    newx_struct = as.matrix(test[, grepl('17Networks', colnames(test)) | colnames(test) %in% c("sex_bin", "euler", "meanFD")])
    newx_func = as.matrix(test[, grepl('_to_', colnames(test)) | colnames(test) %in% c("sex_bin", "euler", "meanFD")])
    newy = as.matrix(test[, colnames(test) %in% c("age")])
    
    # Fitting a ridge regression model
    lambda_seq <- c(
      seq(10, 1, -1),        # From 10 to 1 (step of 1)
      seq(0.9, 0.1, -0.1),   # From 0.9 to 0.1 (step of 0.1)
      seq(0.09, 0.01, -0.01),# From 0.09 to 0.01 (step of 0.01)
      seq(0.009, 0.001, -0.001), # From 0.009 to 0.001 (step of 0.001)
      seq(0.0009, 0.0001, -0.0001), # From 0.0009 to 0.0001 (step of 0.0001)
      seq(0.00009, 0.00001, -0.00001), # From 0.00009 to 0.00001
      seq(0.000009, 0.000001, -0.000001), # From 0.000009 to 0.000001
      seq(0.0000009, 0.0000001, -0.0000001), # From 0.0000009 to 0.0000001
      0
    )
    
    # Training respective models
    ridgeModel_struct <- cv.glmnet(x_var_struct, y_var, alpha = 0, lambda = lambda_seq)
    ridgeModel_func <- cv.glmnet(x_var_func, y_var, alpha = 0, lambda = lambda_seq)
    
    # Getting predictions
    test$preds_struct = predict(ridgeModel_struct, newx = newx_struct, s = "lambda.min")[,]
    test$preds_func = predict(ridgeModel_func, newx = newx_func, s = "lambda.min")[,]
    
    # Storing prediction values
    ## Isolating to a singular data frame with only the prediction values and IDs
    test_collapsed <- test %>%
      group_by(participant_id) %>%
      summarise(
        preds_struct = first(na.omit(preds_struct)),
        preds_func = first(na.omit(preds_func)),
        .groups = "drop"
      )
    
    ## Matching the predictions by ID to the original data frame
    data_temp <- data_temp %>%
      left_join(test_collapsed, by = "participant_id", suffix = c("", "_new"))%>%
      mutate(
        preds_struct = ifelse(is.na(preds_struct_new), preds_struct, preds_struct_new),
        preds_func = ifelse(is.na(preds_func_new), preds_func, preds_func_new)
      ) %>%
      select(-ends_with("_new"))
    
    
    # Storing minimum lambda used for each ridge regression model
    lambda_struct = c(lambda_struct, ridgeModel_struct$lambda.min)
    lambda_func = c(lambda_func, ridgeModel_func$lambda.min)
  }
  
  # Storing results from each model 
  data_clean$struct_ss <- cbind(data_clean$struct_ss, as.matrix(data_temp$preds_struct))
  data_clean$func_ss <- cbind(data_clean$func_ss, as.matrix(data_temp$preds_func))
  data_clean$fold_ss <- cbind(data_clean$fold_ss, as.matrix(data_temp$fold))
}
```




### Applying `mapa` to Data Predictions


After fitting the results from the 20 sample splits, we have 20 sets of predictions for each model and fold assignment.
The `mapa` function takes matrix valued arguments for the structural predictions, fold ID, and functional predictions to compute the MAPA estimates and confidence intervals across the sample splits.
Each column of the prediction matrices corresponds to one configuration of five-fold cross-fitting.

```{r}
# Gathering the estimates with four sample splits
results = mapa(data_clean$age, data_clean$struct_ss, data_clean$fold_ss, data_clean$func_ss)

# Here, yhat1 represents the prediction accuracy of the structural model and 
# yhat2 represents the prediction accuracy of the functional model 
results$est <- results$est %>% mutate(metric = recode(metric, yhat1 = "Structural Data", yhat2 = "Functional Data"))
results$est_summary <- results$est_summary %>% mutate(metric = recode(metric, yhat1 = "Structural Data", yhat2 = "Functional Data"))

# Previewing summarized results
results$est_summary
head(results$est)
```

### Graphing Outputs

The results here align with the previous section; the model trained on structural data performs better than the model trained on functional data when predicting participant age.

```{r}
# Filtering to only the estimates for the structural and functional estimates
toPlot <- results$est_summary %>% filter((metric == "Structural Data" | metric == "Functional Data") & scale == "r")
toPlot_Diff <- results$est_summary %>% filter(metric == "Difference")

# Plotting results
ggplot(toPlot, aes(x = metric, y = est)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = LB, ymax = UB),
                width = 0.2,
                position = position_dodge(width = 0.5)) +
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust=1)) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "black")  + 
  labs(x = "Ridge Regression Model Data", y = "One-Step Estimator Value")


# Plotting one-step difference results 
ggplot(toPlot_Diff, aes(x = metric, y = est)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = LB, ymax = UB),
                width = 0.2,
                position = position_dodge(width = 0.5)) +
  theme_bw() + 
  theme(axis.ticks.y = element_blank(), axis.text.y = element_blank()) + 
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "black")  + 
  ylim(-1,1) + 
  labs(x = "One-Step Difference", y = "Difference Magnitude")
```

 

# Random Forest Model Examples

This section parallels the previous analyses using random forest instead of ridge regression. 

## Methods

Here, outside of the random forest model package itself, we create the *helper_ranger* function to help tune and check the validity of our model tuning.

```{r}
# tuning function
helper_ranger = function(trainDat, outcome, testDat, task){
  # If the tuned model already gives non-constant preds, keep it.
  if (sd(predict(task$model, newdata = data.frame(testDat))$data$response) > 1e-5) {
    return(task$model)
  }

  # Otherwise walk the tried hyperparams (best MSE first) and re-train via mlr
  tune_table <- task$results[order(task$results$mse), ]

  for (i in seq_len(nrow(tune_table))) {
    params <- tune_table[i, ]

    # build an mlr learner with these params
    lrn <- mlr::makeLearner(
      "regr.ranger",
      num.trees       = 250,
      mtry            = params$mtry,
      min.node.size   = params$min.node.size,
      sample.fraction = params$sample.fraction,
      num.threads     = 1,
      respect.unordered.factors = "order"  # or as needed
    )

    # train an mlr model so predict() keeps the same API
    tr_task <- mlr::makeRegrTask(id = "rf_task_retry", data = trainDat, target = outcome)
    wm <- mlr::train(lrn, tr_task)

    # check for non-constant predictions
    if (sd(predict(wm, newdata = data.frame(testDat))$data$response) > 1e-5) {
      return(wm)
    }
  }

  # if everything fails, just return the original
  return(task$model)
}
```

## Random Forest with Single Cross-Fitting

### Getting Predicted Values
Here, we test how the inclusion of either structural or functional data in a random forest model impacts model accuracy when predicting age. Five-fold cross-fitting is used to obtain training predictions for the left-out fold.

```{r}
# Cleaning data - drop NA before CV
data_clean = pnc_data %>% drop_na(age) 

# Assigning folds randomly
folds = 5
data_clean$fold <- sample(rep(1:folds, length.out = nrow(data_clean)))

# Selecting an outcome variable 
outcome = "age"

# Sanity check for fold distribution 
table(data_clean$fold)

# Creating storage columns
data_clean$preds_struct = NA
data_clean$preds_func = NA

# Creating a for loop for 5-fold cross validation
for(i in 1:folds){
  # Assigning data
  test = data_clean %>% filter(fold == i)
  train = data_clean %>% filter(fold != i)
  
  # Creating model data variables for training  
  x_var_struct = data.frame(train[, grepl('17Networks', colnames(train)) | colnames(train) %in% c(outcome, "sex_bin", "euler", "meanFD")])
  y_var = data.frame(train[, colnames(train) %in% c("age")])
  x_var_func = data.frame(train[, grepl('_to_', colnames(train)) | colnames(train) %in% c(outcome, "sex_bin", "euler", "meanFD")])

  # Creating model variables for prediction
  newx_struct = data.frame(test[, grepl('17Networks', colnames(test)) | colnames(test) %in% c("sex_bin", "euler", "meanFD")])
  newx_func = data.frame(test[, grepl('_to_', colnames(test)) | colnames(test) %in% c("sex_bin", "euler", "meanFD")])
  newy = data.frame(test[, colnames(test) %in% c("age")])
  
  # Creating random forest tasks
  task_struct = makeRegrTask(id = "rf_task", data = x_var_struct, target = outcome)
  task_func = makeRegrTask(id = "rf_task", data = x_var_func, target = outcome)
  
  # Tuning each model
  rfTune_struct = tuneRanger(task_struct, num.trees = 250, iters = 25, iters.warmup = 10, show.info = FALSE)
  rfTune_func = tuneRanger(task_func, num.trees = 250, iters = 25, iters.warmup = 10, show.info = FALSE)
  
  # Checking the tuning of each model
  rfTune_struct$model = helper_ranger(x_var_struct, outcome, newx_struct, rfTune_struct) # Calling for checking of tuning 
  rfTune_func$model = helper_ranger(x_var_func, outcome, newx_func, rfTune_func) # Calling for checking of tuning 
  
  # Getting predictions
  test$preds_struct = predict(rfTune_struct$model, newdata = newx_struct)$data$response
  test$preds_func = predict(rfTune_func$model, newdata = newx_func)$data$response
  
  # Storing prediction values
  ## Isolating to a singular data frame with only the prediction values and IDs
  test_collapsed <- test %>%
  group_by(participant_id) %>%
  summarise(
    preds_struct = first(na.omit(preds_struct)),
    preds_func= first(na.omit(preds_func)),
    .groups = "drop"
  )
  
  ## Matching the predictions by ID to the original data frame
  data_clean <- data_clean %>%
  left_join(test_collapsed, by = "participant_id", suffix = c("", "_new"))%>%
  mutate(
    preds_struct = ifelse(is.na(preds_struct_new), preds_struct, preds_struct_new),
    preds_func = ifelse(is.na(preds_func_new), preds_func, preds_func_new)
  ) %>%
  select(-ends_with("_new"))
}
```

### Applying `mapa` to Data Predictions

The approach is the same regardless of the model used, assuming it has been tuned appropriately to satisfy the convergence criteria for the MAPA estimator.

```{r}
# Gathering the estimates with one sample split
results = mapa(as.matrix(data_clean$age), as.matrix(data_clean$preds_struct), as.matrix(data_clean$fold), as.matrix(data_clean$preds_func))

# Here, yhat1 represents the prediction accuracy of the structural model and 
# yhat2 represents the prediction accuracy of the functional model 
results$est <- results$est %>% mutate(metric = recode(metric, yhat1 = "Structural Data", yhat2 = "Functional Data"))
results$est_summary <- results$est_summary %>% mutate(metric = recode(metric, yhat1 = "Structural Data", yhat2 = "Functional Data"))

# Previewing summarized results
results$est_summary
head(results$est)
```

### Graphing Outputs

As with the ridge regression model, we see that the model trained on structural data performs significantly better than the model trained on functional data when predicting subject age. 

```{r}
# Subsetting results to only estimates
toPlot = results$est_summary %>% filter((metric == "Structural Data" | metric == "Functional Data") & scale == "r")
toPlot_Diff =  results$est_summary %>% filter(metric == "Difference")

# Plotting one-step results 
ggplot(toPlot, aes(x = metric, y = est)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = LB, ymax = UB),
                width = 0.2,
                position = position_dodge(width = 0.5)) +
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust=1)) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "black")  + 
  labs(x = "Random Forest Model Data", y = "One-Step Estimator Value")


# Plotting one-step difference results 
ggplot(toPlot_Diff, aes(x = metric, y = est)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = LB, ymax = UB),
                width = 0.2,
                position = position_dodge(width = 0.5)) +
  theme_bw() + 
  theme(axis.ticks.y = element_blank(), axis.text.y = element_blank()) + 
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "black")  + 
  ylim(-1,1) + 
  labs(x = "One-Step Difference", y = "Difference Magnitude")
```


## Random Forest using Sample Splits

### Getting Predicted Values
```{r}
# Cleaning data
data_clean = pnc_data %>% drop_na(age)

# Assigning sample split number and fold number 
ss_num = 20
folds = 5

# Creating storage variables for hypertuning
lambda_struct = c()
lambda_func = c()

# Assigning storage matrices within data frame columns 
data_clean$struct_ss  <- matrix(nrow = nrow(data_clean), ncol = 0)
data_clean$func_ss  <- matrix(nrow = nrow(data_clean), ncol = 0)
data_clean$fold_ss  <- matrix(nrow = nrow(data_clean), ncol = 0)

# Running the sample splits
for(j in 1:ss_num){
  # Creating a temporary data frame
  data_temp <- data_clean
  
  # Creating storage columns
  data_temp$preds_struct = NA
  data_temp$preds_func = NA
  
  # Setting fold numbers for temporary data frame
  data_temp$fold <- sample(rep(1:folds, length.out = nrow(data_temp)))
  
  
  # Creating a for loop for 5-fold cross validation
  for(i in 1:folds){
    # Assigning data
    test = data_temp %>% filter(fold == i)
    train = data_temp %>% filter(fold != i)
    
    # Creating model data variables for training  
    x_var_struct = data.frame(train[, grepl('17Networks', colnames(train)) | colnames(train) %in% c(outcome, "sex_bin", "euler", "meanFD")])
    y_var = data.frame(train[, colnames(train) %in% c("age")])
    x_var_func = data.frame(train[, grepl('_to_', colnames(train)) | colnames(train) %in% c(outcome, "sex_bin", "euler", "meanFD")])
  
    # Creating model variables for prediction
    newx_struct = data.frame(test[, grepl('17Networks', colnames(test)) | colnames(test) %in% c("sex_bin", "euler", "meanFD")])
    newx_func = data.frame(test[, grepl('_to_', colnames(test)) | colnames(test) %in% c("sex_bin", "euler", "meanFD")])
    newy = data.frame(test[, colnames(test) %in% c("age")])
    
    # Creating random forest tasks
    task_struct = makeRegrTask(id = "rf_task", data = x_var_struct, target = outcome)
    task_func = makeRegrTask(id = "rf_task", data = x_var_func, target = outcome)
    
    # Tuning each model
    rfTune_struct = tuneRanger(task_struct, num.trees = 250, iters = 25, iters.warmup = 10, show.info = FALSE)
    rfTune_func = tuneRanger(task_func, num.trees = 250, iters = 25, iters.warmup = 10, show.info = FALSE)
    
    # Checking the tuning of each model
    rfTune_struct$model = helper_ranger(x_var_struct, outcome, newx_struct, rfTune_struct) # Calling for checking of tuning 
    rfTune_func$model = helper_ranger(x_var_func, outcome, newx_func, rfTune_func) # Calling for checking of tuning 
    
    # Getting predictions
    test$preds_struct = predict(rfTune_struct$model, newdata = newx_struct)$data$response
    test$preds_func = predict(rfTune_func$model, newdata = newx_func)$data$response
    
    # Storing prediction values
    ## Isolating to a singular data frame with only the prediction values and IDs
    test_collapsed <- test %>%
      group_by(participant_id) %>%
      summarise(
        preds_struct = first(na.omit(preds_struct)),
        preds_func = first(na.omit(preds_func)),
        .groups = "drop"
      )
    
    ## Matching the predictions by ID to the original data frame
    data_temp <- data_temp %>%
      left_join(test_collapsed, by = "participant_id", suffix = c("", "_new"))%>%
      mutate(
        preds_struct = ifelse(is.na(preds_struct_new), preds_struct, preds_struct_new),
        preds_func = ifelse(is.na(preds_func_new), preds_func, preds_func_new)
      ) %>%
      select(-ends_with("_new"))
  }
  
  # Storing results from each model 
  data_clean$struct_ss <- cbind(data_clean$struct_ss, as.matrix(data_temp$preds_struct))
  data_clean$func_ss <- cbind(data_clean$func_ss, as.matrix(data_temp$preds_func))
  data_clean$fold_ss <- cbind(data_clean$fold_ss, as.matrix(data_temp$fold))
}
```

### Applying `mapa` to Data Predictions
```{r}
# Gathering the estimates with four sample splits
results = mapa(data_clean$age, data_clean$struct_ss, data_clean$fold_ss, data_clean$func_ss)

# Here, yhat1 represents the prediction accuracy of the structural model and 
# yhat2 represents the prediction accuracy of the functional model 
results$est <- results$est %>% mutate(metric = recode(metric, yhat1 = "Structural Data", yhat2 = "Functional Data"))
results$est_summary <- results$est_summary %>% mutate(metric = recode(metric, yhat1 = "Structural Data", yhat2 = "Functional Data"))

# Previewing summarized results
results$est_summary
head(results$est)
```

### Graphing Outputs


Here, with the inclusion of sample splits to test other scenarios where fold assignment randomly differs, we see that the model trained on structural data does not perform better than the model trained on functional data when predicting subject age (the confidence interval for the difference includes zero). Further, the plausible range of difference in correlation is between 0.07 and 0.14 between the structural and functional models. 


```{r}
# Filtering to only the estimates for the eTIV and nWBV estimates
toPlot <- results$est_summary %>% filter((metric == "Structural Data" | metric == "Functional Data") & scale == "r")
toPlot_Diff <- results$est_summary %>% filter(metric == "Difference")

# Plotting one-step results 
ggplot(toPlot, aes(x = metric, y = est)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = LB, ymax = UB),
                width = 0.2,
                position = position_dodge(width = 0.5)) +
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust=1)) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "black")  + 
  labs(x = "Random Forest Model Data", y = "One-Step Estimator Value")



# Plotting one-step difference results 
ggplot(toPlot_Diff, aes(x = metric, y = est)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = LB, ymax = UB),
                width = 0.2,
                position = position_dodge(width = 0.5)) +
  theme_bw() + 
  theme(axis.ticks.y = element_blank(), axis.text.y = element_blank()) + 
  coord_flip() + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "black")  + 
  ylim(-1,1) + 
  labs(x = "One-Step Difference", y = "Difference Magnitude")
```