Train PLS Smoking

knitr::opts_chunk$set(echo = TRUE)

library(caret)
library(R.matlab)
library(tidyverse)
library(pROC)
library(broom)
library(kableExtra)
library(rlang)
library(testthat)
library(doParallel)

source("train_pls_formulae.R")

n_cores <- detectCores()
message(sprintf("Using %d cores", n_cores))

## Using 20 cores

Note: the output HTML & the training RDS data are saved in PLS.

Parameters

name	value
Cohort name	UKBB
Risk factor	Smoking
Maximum number of PLS components	30
Pct PCA variance explained (%)	99.9
Data file	MESA_UKBB_FinalData.csv

Get the data

1. Read the variable and select only cases within the cohort UKBB
2. Select ID, Age, Gender & the risk factor Smoking.
3. Recode the risk factor values where TRUE -> RISK and FALSE -> NO_RISK

dt.ori <- read_csv(datafile, show_col_types = FALSE) %>% 
  filter(Cohort == params$cohort) %>%
  select(c(ID, Age, Sex, NoRisk, sym(params$rf_name))) %>%
  mutate(
    !!sym(params$rf_name) := recode_factor(as.factor(!!sym(params$rf_name)), "FALSE" = "NO_RISK", "TRUE" = "RISK"),
    Sex = as.factor(Sex)
  )

message(sprintf("Original %s: %d rows", params$cohort, nrow(dt.ori)))

## Original UKBB: 2960 rows

4. Define the NO RISK sub-cohort, which is all cases without any risk factor at all

dt.norf <- filter(dt.ori, NoRisk)

message(sprintf("No risk sub-cohort %s: %d rows", params$cohort, nrow(dt.norf)))

## No risk sub-cohort UKBB: 645 rows

5. Filter cases having Smoking

dt.rf <- filter(dt.ori,  !!sym(params$rf_name) == "RISK")

message(sprintf("Risk sub-cohort %s: %d rows", params$cohort, nrow(dt.rf)))

## Risk sub-cohort UKBB: 1182 rows

Sanity check: the intersection between NO_RISK and RISK id’s should be zero

test_that("Unique IDs", {
  expect_equal(length(intersect(dt.rf$ID, dt.norf$ID)), 0)
})

## Test passed 🎊

6. Final data

dt <- bind_rows(dt.norf, dt.rf) %>% 
  select(-c(NoRisk))

message(sprintf("Training data: %d rows", nrow(dt)))

## Training data: 1827 rows

dt %>% select(c(Sex, sym(params$rf_name))) %>% table() %>% 
  kbl(caption = "Risk factor distributions") %>% kable_styling("hover", full_width = F, position="left")

Risk factor distributions

	NO_RISK	RISK
female	395	560
male	250	622

just for fixing plot title

cohort_title <- if_else(params$cohort == "UKBB", "UK Biobank", params$cohort)  

Prepare the PC scores

1. Determine how many PCA components to use

model <- list()

expl_vars <- readMat(paste("PCA", params$cohort, "PCA_explained.mat", sep = "/"))$explained
model$npcs <- length(expl_vars[cumsum(expl_vars) < params$pca_var])
message(sprintf("Number of PC components with %.2f%% variance explained: %d", params$pca_var, model$npcs))

## Number of PC components with 99.90% variance explained: 210

rm(expl_vars)

2. Read PCA components
3. Combine with Age, Sex & the risk factor using the same ID

# this is a combined X (Age + Sex + PCScores) & Y (Risk Factor)
model$dt.train <- inner_join(
  dt %>% 
    select(c(ID, Age, Sex, sym(params$rf_name))) %>% 
    mutate(
      Sex = as.numeric(Sex)
    ),
  cbind(
    read.csv(paste("PCA", sprintf("%s_ids.csv", params$cohort), sep="/")),
    readMat(paste("PCA", params$cohort, "PCA_score.mat", sep="/"))$score[, 1:model$npcs]
  ),
  by = c("ID" = "ID"))

# check the cohort
message(sprintf("Number of cases = %d, columns = %d", nrow(model$dt.train), ncol(model$dt.train)))

## Number of cases = 1827, columns = 214

4. Sanity checks.

This should only include that specific cohort

test_that("Check cohort", {
  expect_equal(substr(params$cohort, 1, 3), (model$dt.train %>% transmute(cohort = substr(ID, 1, 3)) %>% distinct())$cohort)
})

## Test passed 😀

No NA rows

test_that("No NA rows", {
  expect_equal(0, sum(is.na(model$dt.train)))
})

## Test passed 😸

Find optimal number of PLS components

We’re going to use 5-fold cross validation to determine the optimal parameter for PLS, which is the ncomp.

(m_kcv <- train_pls(as.formula(sprintf("%s ~ Age + Sex + .", params$rf_name)), dt = model$dt.train %>% select(-c(ID)), n_comps=params$max_ncomps))

## + Fold1: ncomp=30 
## - Fold1: ncomp=30 
## + Fold2: ncomp=30 
## - Fold2: ncomp=30 
## + Fold3: ncomp=30 
## - Fold3: ncomp=30 
## + Fold4: ncomp=30 
## - Fold4: ncomp=30 
## + Fold5: ncomp=30 
## - Fold5: ncomp=30 
## Aggregating results
## Selecting tuning parameters
## Fitting ncomp = 6 on full training set

## Partial Least Squares 
## 
## 1827 samples
##  212 predictor
##    2 classes: 'NO_RISK', 'RISK' 
## 
## Pre-processing: centered (212) 
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 1461, 1462, 1461, 1462, 1462 
## Resampling results across tuning parameters:
## 
##   ncomp  ROC        Sens        Spec     
##    1     0.6074807  0.02945736  0.9754380
##    2     0.6556874  0.23100775  0.9018022
##    3     0.6748175  0.26356589  0.8908317
##    4     0.6854824  0.30387597  0.8688407
##    5     0.6873838  0.31782946  0.8654509
##    6     0.6874211  0.32558140  0.8603840
##    7     0.6827327  0.33333333  0.8553243
##    8     0.6771454  0.33643411  0.8434778
##    9     0.6724043  0.34418605  0.8384074
##   10     0.6681054  0.35658915  0.8307874
##   11     0.6689635  0.36434109  0.8274262
##   12     0.6712258  0.36434109  0.8223343
##   13     0.6728020  0.38449612  0.8104663
##   14     0.6713925  0.38449612  0.8113280
##   15     0.6738777  0.39534884  0.8087857
##   16     0.6742345  0.39069767  0.8071051
##   17     0.6752221  0.39844961  0.8037188
##   18     0.6704259  0.39689922  0.8011979
##   19     0.6702714  0.39534884  0.7986519
##   20     0.6629029  0.39534884  0.7918723
##   21     0.6590430  0.39379845  0.7901809
##   22     0.6584372  0.40310078  0.7901809
##   23     0.6584101  0.40620155  0.7910355
##   24     0.6587308  0.40000000  0.7884932
##   25     0.6583679  0.39379845  0.7935636
##   26     0.6560695  0.39844961  0.7808839
##   27     0.6559729  0.40465116  0.7850998
##   28     0.6553352  0.40930233  0.7791819
##   29     0.6545037  0.40775194  0.7758028
##   30     0.6540642  0.41395349  0.7800365
## 
## ROC was used to select the optimal model using the largest value.
## The final value used for the model was ncomp = 6.

Get the number of component

(model$ncomp <- m_kcv$bestTune$ncomp)

## [1] 6

Plot for analysis

plot(m_kcv, main=sprintf("PLS components (%s - %s)", params$rf_name, params$cohort))

Train PLS

After we get the optimal parameter, we train the PLS using leave-one-out cross validation. Since this will take time, we need multiple cores to compute.

Create train control with LOOCV

tc_loo <- trainControl(method="LOOCV", savePredictions = "all", classProbs = TRUE, verboseIter=FALSE,
                       summaryFunction = twoClassSummary, allowParallel = TRUE)

# parallel cores
if( n_cores > 1 ) {
  cl <- makeCluster(n_cores - 1, outfile = "")
  registerDoParallel(cl)
  getDoParWorkers()
}

## [1] 19

# train (take a while)
model$pls <- train(
  form = sprintf("%s ~ Age + Sex + .", params$rf_name) %>% as.formula(),
  data = model$dt.train %>% select(-ID),
  method = "pls", 
  metric = "ROC",
  preProc = c("center"),
  tuneGrid = data.frame(ncomp=model$ncomp),
  probMethod = "softmax",
  trControl = tc_loo
)

# unregister cores
if( n_cores > 1 ) {
  stopCluster(cl)
  registerDoSEQ()
}

model$pls

## Partial Least Squares 
## 
## 1827 samples
##  212 predictor
##    2 classes: 'NO_RISK', 'RISK' 
## 
## Pre-processing: centered (212) 
## Resampling: Leave-One-Out Cross-Validation 
## Summary of sample sizes: 1826, 1826, 1826, 1826, 1826, 1826, ... 
## Resampling results:
## 
##   ROC        Sens       Spec     
##   0.6935073  0.3069767  0.8663283
## 
## Tuning parameter 'ncomp' was held constant at a value of 6

Or use getTrainPerf to get the results

getTrainPerf(model$pls)

##    TrainROC TrainSens TrainSpec method
## 1 0.6935073 0.3069767 0.8663283    pls

The prediction is in model$pls$pred data frame. Let’s plot the training ROC

Check the level order, because RISK should be the first, since that’s the prediction target.

levels(model$pls$pred$obs)

## [1] "NO_RISK" "RISK"

Calculate the ROC. Note we need to reverse the category level.

(m_roc <- roc(model$pls$pred$obs, model$pls$pred$RISK, levels = rev(levels(model$pls$pred$obs))))

## Setting direction: controls > cases

## 
## Call:
## roc.default(response = model$pls$pred$obs, predictor = model$pls$pred$RISK,     levels = rev(levels(model$pls$pred$obs)))
## 
## Data: model$pls$pred$RISK in 1182 controls (model$pls$pred$obs RISK) > 645 cases (model$pls$pred$obs NO_RISK).
## Area under the curve: 0.6935

plot(m_roc, legacy_axes=TRUE, main = sprintf("Training ROC: %s (%s)", params$rf_name, params$cohort))

Show the first 10 coefficients

coef(model$pls$finalModel)[1:10,,]

##           NO_RISK          RISK
## Age -1.248862e-03  1.248862e-03
## Sex  8.620371e-06 -8.620371e-06
## `1`  3.675642e-04 -3.675642e-04
## `2`  2.642924e-04 -2.642924e-04
## `3` -3.808875e-04  3.808875e-04
## `4` -7.600182e-04  7.600182e-04
## `5`  8.175104e-04 -8.175104e-04
## `6`  1.004172e-03 -1.004172e-03
## `7`  3.190593e-04 -3.190593e-04
## `8`  6.972424e-04 -6.972424e-04

Save the results

Add some other information too in the model

model$cohort <- params$cohort
model$risk_factor <- params$rf_name
model$pca_var <- params$pca_var

model_filename <- paste(params$output_folder, sprintf("PLS_%s_%s.rds", params$cohort, params$rf_name), sep="/")
write_rds(model, file=model_filename, compress = "gz")
message(sprintf("Saved model to %s", model_filename))

## Saved model to PLS/PLS_UKBB_Smoking.rds