---
title: "Outlier Ensemble Analysis of Triple-Negative Breast Cancer (TNBC) Data"
author: "Marta B. Lopes, André Veríssimo, Eunice Carrasquinha, Sandra Casimiro, Niko Beerenwinkel, Susana Vinga"
date: "16-01-2018"
output:
  html_document:
    number_sections: yes
    toc: yes
  output: null
  pdf_document:
    number_sections: yes"N"
    toc: yes
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning=FALSE)
```


The data analysis described below requires R version 3.3.2 and Bioconductor version 3.4.

# The BRCA dataset

The Cancer Genome Atlas Breast Invasive Carcinoma (TCGA-BRCA) data collection is part of a larger effort to build a research community focused on connecting cancer phenotypes to genotypes by providing clinical images matched to subjects from The Cancer Genome Atlas ([TCGA](https://cancergenome.nih.gov/)).

The BRCA data is publicly available (https://cancergenome.nih.gov/) and is decomposed into two datasets:
  
  1. the gene expression data, composed of 57251 variables for a total of 1222 samples with 1097 individuals. From those samples, 1102 with primary solid tumor, 7 metastatic and 113 with normal tissue;

2. the clinical data is composed of 114 variables obtained from the same individuals (much more is available in the gdc data variable, such as follow_up, drug, radiation, ...)


The BRCA datasets will be imported using the package 'brca.data'.

```{r}
# Package developed by our group
library(devtools)
devtools::install_url('https://github.com/averissimo/brca.data/releases/download/1.0/brca.data_1.0.tar.gz')
library(brca.data)
```

## TNBC dataset construction

Triple negative breast cancer (TNBC) refers to any BRCA that does not express the genes for the estrogen receptor (ESR1; ENSG00000091831), progesterone receptor (PGR; ENSG00000082175) and which does not overexpress the human epidermal growth factor receptor 2 (HER2; ENSG00000141736).

From the clinical data, the columns related to gene expression status of ESR1, PGR and HER2 must be selected. There are three possible variables for HER2:
  
  (a) ‘her2_immunohistochemistry_level_result’, which refers to the ImmunoHistoChemistry (IHC) HER2 result into levels ‘0’ (negative), ‘1’ (negative), ‘2+’ (indeterminate) and ‘3+’(positive);

(b) ‘lab_proc_her2_neu_immunohistochemistry_receptor_status’, corresponding to the IHC result into status indeterminate, equivocal, negative and positive;

(c) ‘lab_procedure_her2_neu_in_situ_hybrid_outcome_type’, the HER2 individuals’ classiﬁcation by the ﬂuorescence in situ hybridization (FISH) method, considered more accurate and often performed when IHC results do not show a clear classiﬁcation of cells into HER2-positive or negative.

```{r}
data("fpkm.per.tissue","fpkm.per.tissue.barcode", 'clinical', package = 'brca.data')

tnbc.vars <- c('breast_carcinoma_estrogen_receptor_status',
               'breast_carcinoma_progesterone_receptor_status',
               'her2_immunohistochemistry_level_result',
               'lab_proc_her2_neu_immunohistochemistry_receptor_status',
               'lab_procedure_her2_neu_in_situ_hybrid_outcome_type'
)
clinical.tnbc <- clinical$primary.solid.tumor[, tnbc.vars]

clinical.tnbc.aux <- clinical.tnbc
names(clinical.tnbc.aux) <- c("ESR1","PGR","HER2_level","HER2_status","HER2_FISH")

summary(clinical.tnbc.aux)
which(is.na(clinical.tnbc.aux[,1]))
clinical.tnbc.aux <- clinical.tnbc.aux[-(which(is.na(clinical.tnbc.aux[,1]))),]
```

Finding the concordance between the HER2 variables measured.

Concordance between HER2 (IHC) level and HER2 (IHC) status:
  
```{r}
table(as.factor(clinical.tnbc.aux[,3]),clinical.tnbc.aux[,4])
```

Creating a vector of HER2 (IHC) status based on HER2 (ICH) level:
  
```{r}
clinical.tnbc.aux$my_HER2_level <- as.character(clinical.tnbc.aux$HER2_level)

clinical.tnbc.aux[which(clinical.tnbc.aux$HER2_level == '0'),6] <- 'Negative'
clinical.tnbc.aux[which(clinical.tnbc.aux$HER2_level == '1+'),6] <- 'Negative'
clinical.tnbc.aux[which(clinical.tnbc.aux$HER2_level == '2+'),6] <- 'Indeterminate'
clinical.tnbc.aux[which(clinical.tnbc.aux$HER2_level == '3+'),6] <- 'Positive'
```

Individuals for which HER2 (IHC) status different from myHER2 (IHC) status (based on HER2 IHC level)

```{r}
HER2status.dif.HER2level <- clinical.tnbc.aux[c(rownames(clinical.tnbc.aux[which(clinical.tnbc.aux[,6] == 'Negative' & clinical.tnbc.aux[,4] == 'Positive'),]),rownames(clinical.tnbc.aux[which(clinical.tnbc.aux[,6] == 'Positive' & clinical.tnbc.aux[,4] == 'Negative'),])),]

# getting the individuals short name, as in the clinical data
clinical.tnbc.barcode.pattern <- '(TCGA-[A-Z0-9a-z]{2}-[a-zA-Z0-9]{4})-([0-9]{2}).+'
clinical.tnbc.shortname <- gsub(clinical.tnbc.barcode.pattern,'\\1',rownames(clinical.tnbc))

# individuals with discordant HER2 (IHC) status and HER2 (IHC) level
HER2status.dif.HER2level.names <- clinical.tnbc.shortname[which(clinical.tnbc.shortname %in% rownames(HER2status.dif.HER2level))]
length(HER2status.dif.HER2level.names)
```

Concordance between HER2 (IHC) status and HER2 (FISH):
  
```{r}
table(clinical.tnbc.aux[,5],clinical.tnbc.aux[,4])
```

Individuals for which HER2 (IHC) status different from HER2 (FISH):
  
```{r}
HER2status.dif.FISH <- clinical.tnbc.aux[c(rownames(clinical.tnbc.aux[which(clinical.tnbc.aux[,4] == 'Negative' & clinical.tnbc.aux[,5] == 'Positive'),]),rownames(clinical.tnbc.aux[which(clinical.tnbc.aux[,4] == 'Positive' & clinical.tnbc.aux[,5] == 'Negative'),])),]

HER2status.dif.HER2FISH.names <- clinical.tnbc.shortname[which(clinical.tnbc.shortname %in% rownames( HER2status.dif.FISH))]
length(HER2status.dif.HER2FISH.names)

# suspect individuals with discordante HER2 labels from different lab testing
suspect <- c(HER2status.dif.HER2level.names,HER2status.dif.HER2FISH.names)
length(suspect)
```

Replacing HER2 (IHC) status by HER2 (FISH) (a more accurate test) whenever available:
  
```{r}
clinical.tnbc.aux$HER2_status_plus_FISH <- clinical.tnbc.aux[,4]
clinical.tnbc.aux[clinical.tnbc.aux[,5] !='',7] <- clinical.tnbc.aux[clinical.tnbc.aux[,5] !='',5]
```

Concordance between HER2 (IHC) status and HER2 (IHC status + FISH):
  
```{r}
table(clinical.tnbc.aux[,4],clinical.tnbc.aux[,7])
```

Building the final clinical TNBC dataset, composed of all individuals and gene expression status of ESR1, PGR and HER2:
  
```{r}
clinical.tnbc <- clinical.tnbc.aux[,c(1,2,7)]

# replace all '' or NA with 'Indeterminate'
clinical.tnbc[is.na(clinical.tnbc) | clinical.tnbc == '' | clinical.tnbc == 'Equivocal'] <- 'Indeterminate'

# individuals with at least one 'Positive' and keep all are marked as not TNBC
tnbc.status.pos <- clinical.tnbc == 'Positive'
#tnbc.status.pos <- clinical.tnbc[,c(1,2,3)] == 'Positive'
not.tnbc.id     <- rownames(tnbc.status.pos[which(sapply(seq(nrow(tnbc.status.pos)), function(ix) { any(tnbc.status.pos[ix,]) })),])

# individuals with at least one 'Indeterminate' and excluded from dataset
tnbc.status.ind <- clinical.tnbc[!(rownames(clinical.tnbc) %in% not.tnbc.id),] == 'Indeterminate'
not.tnbc.id     <- rownames(tnbc.status.ind[which(sapply(seq(nrow(tnbc.status.ind)), function(ix) { any(tnbc.status.ind[ix,]) })),])
clinical.tnbc <- clinical.tnbc[!(rownames(clinical.tnbc) %in% not.tnbc.id), ]

tnbc.status.neg <- clinical.tnbc == 'Negative'
tnbc.ix <- sapply(seq(nrow(clinical.tnbc)), function(ix) { all(tnbc.status.neg[ix,]) })

# set cases with TNBC
clinical.tnbc$tnbc <- 'NO_TNBC'
clinical.tnbc[tnbc.ix, 'tnbc'] <- 'TNBC'
```

Building the Y (binary data, i.e., TNBC vs. nonTNBC) response vector:
  
```{r}
# building the ydata (binary)
ydata <- data.frame(tnbc = clinical.tnbc$tnbc, row.names = rownames(clinical.tnbc))
ydata$tnbc <- factor(ydata$tnbc)

```

Protein coding genes catalog for building the X gene expression matrix

I. Install packages

```{r}
source("http://bioconductor.org/biocLite.R")  
biocLite(c('ensembldb', 'EnsDb.Hsapiens.v86'))

# Load libraries

suppressPackageStartupMessages({
  library(ensembldb)
  library(EnsDb.Hsapiens.v86)
  library(futile.logger)
  .Last.value <- flog.layout(layout.format('~m'))
})
```

II. Load ensembl data

Ensembl is a genome browser for vertebrate genomes that supports research in comparative genomics, evolution, sequence variation and transcriptional regulation. Ensembl annotate genes, computes multiple alignments, predicts regulatory function and collects disease data. Ensembl tools include BLAST, BLAT, BioMart and the Variant Effect Predictor (VEP) for all supported species.

- [Home page](http://www.ensembl.org/index.html)

```{r}
edb <- EnsDb.Hsapiens.v86
# hasProteinData(edb) # from devel only, as of 3.4
ensembl.protein.coding <- genes(edb, 
                                filter = list(GenebiotypeFilter('protein_coding')), 
                                columns = c('gene_id', 'gene_name'))
{
  flog.info('         Granges: %d', nrow(ensembl.protein.coding@elementMetadata))
  flog.info('Metadata columns: %d', ncol(ensembl.protein.coding@elementMetadata))
}
```

III. Load ccds data

The Consensus CDS (CCDS) project is a collaborative effort to identify a core set of human and mouse protein coding regions that are consistently annotated and of high quality. The long term goal is to support convergence towards a standard set of gene annotations. 

- [Home page](https://www.ncbi.nlm.nih.gov/projects/CCDS/CcdsBrowse.cgi) 
- [Statistics or release 20](https://www.ncbi.nlm.nih.gov/projects/CCDS/CcdsBrowse.cgi?REQUEST=SHOW_STATISTICS)

```{r}
ccds <- read.table(url('ftp://ftp.ncbi.nih.gov/pub/CCDS/current_human/CCDS.current.txt'),
                   sep = '\t', 
                   header = T, 
                   comment.char = "|", # necessary as header line has a # character
                   stringsAsFactors = FALSE)
flog.info('Size of ccds: %d x %d', nrow(ccds), ncol(ccds))
```

IV. Build unique gene set

```{r}
ensembl.genes         <- sort(unique(ensembl.protein.coding@elementMetadata$gene_name))
ensembl.genes.ensg.id <- sort(ensembl.protein.coding@elementMetadata$gene_id)

ccds.genes       <- sort(unique(ccds$gene))
ccds.extra.genes <- sort(ccds.genes[(!ccds.genes %in% ensembl.genes)])
ccds.extra.genes.ensg.id <- genes(edb, filter = list(GenenameFilter(ccds.extra.genes)), 
                                  columns = c('gene_id', 'gene_name'))

ensg.id <- sort(unique(c(ensembl.protein.coding@elementMetadata$gene_id, ensembl.genes.ensg.id)))

all.genes <- sort(unique(c(ensembl.genes, ccds.extra.genes)))
```

Building the X gene expression data matrix:
  
```{r}
xdata <- t(fpkm.per.tissue$primary.solid.tumor[which(rownames(fpkm.per.tissue$primary.solid.tumor) %in% ensg.id),-which(duplicated(fpkm.per.tissue.barcode$primary.solid.tumor))])
common.rows <- getParticipantCode(rownames(xdata)) %in% rownames(ydata)
xdata <- xdata[common.rows,]

# removing variables with sd=0
xdata.sd <- sapply(seq(ncol(xdata)), function(ix) {sd(xdata[,ix])})
xdata <- xdata[,xdata.sd != 0]

# Gene expression data for ESR1, PGR and HER2 from the suspect (potential outlier) cases identified above

# Clinical and gene expression data of individuals with discordant HER2 level and HER2 status
xdata.shortname <- getParticipantCode(rownames(xdata))

xdata.plus.clinical.HER2status.dif.HER2level <- cbind(xdata[which(xdata.shortname %in% HER2status.dif.HER2level.names), c("ENSG00000091831","ENSG00000082175","ENSG00000141736")],clinical.tnbc.aux[which(rownames(clinical.tnbc.aux) %in% HER2status.dif.HER2level.names),],ydata$tnbc[which(rownames(ydata) %in% HER2status.dif.HER2level.names)])

# clinical and gene expression data of individuals with discordant HER2 status and HER2 FISH
xdata.shortname <- getParticipantCode(rownames(xdata))
xdata.plus.clinical.HER2status.dif.HER2FISH <- cbind(xdata[which(xdata.shortname %in% HER2status.dif.HER2FISH.names), c("ENSG00000091831","ENSG00000082175","ENSG00000141736")],clinical.tnbc.aux[which(rownames(clinical.tnbc.aux) %in% HER2status.dif.HER2FISH.names),],ydata$tnbc[which(rownames(ydata) %in% HER2status.dif.HER2FISH.names)])
```

Data transformation to be used in further analysis:
  
```{r}
data.Y <- as.matrix(ydata)
data.Y[data.Y =="NO_TNBC"] <- 0
data.Y[data.Y =="TNBC"] <- 1
data.Y <- as.numeric(data.Y)

# log-transform xdata
data <- log2(xdata+1)

# normalizing xdata
data <- (data-matrix(apply(data,2,mean),nrow(data),ncol(data),byrow=TRUE))/matrix(apply(data,2,sd),nrow(data),ncol(data),byrow=TRUE)
```

# Exploratory data analysis

Summary of TNBC biomarkers:
  
```{r}
summary(xdata[which(data.Y=="1"),c("ENSG00000091831","ENSG00000082175","ENSG00000141736")])
summary(xdata[which(data.Y=="0"),c("ENSG00000091831","ENSG00000082175","ENSG00000141736")])

# Building the training and test sets; keep in the training set observations ESR1- and PGR- that would have the opposite classification if HER2 (FISH) had not replaced the HER2 (IHC) status; and observations ESR1- and PGR- for which FISH was not measured and have no concordance between HER2 (IHC) status and HER2 (IHC) level (having a potential impact on the final TNBC/nonTNBC classification)

# taking individuals ID from the 'HER2status.dif.FISH' and 'HER2status.dif.HER2level' created above
HER2status.dif.FISH
HER2status.dif.HER2level

set.seed(2010)

test.id <- sample((1:nrow(data))[-which(getParticipantCode(rownames(data)) %in% c("TCGA-A2-A0EQ","TCGA-BH-A18T","TCGA-AO-A0JL","TCGA-BH-A18Q","TCGA-A2-A04U","TCGA-LL-A73Y","TCGA-AN-A0FL","TCGA-AN-A0FX"))],round(nrow(data)/4),replace=FALSE)

data.train <- data[-test.id,]
data.Y.train <- data.Y[-test.id]

data.test <- data[test.id,]
data.Y.test <- data.Y[test.id]
```

Fitting a logistic model based on the three biomarker variables (ESR1, PGR and HER2):
  
```{r}
glm.data.train.3var <- data.frame(data.Y.train,data.train[,c("ENSG00000091831","ENSG00000082175","ENSG00000141736")])
colnames(glm.data.train.3var)[1] <- 'data.y'

# without interaction
glm.fit.3var <- glm(data.y~.,family=binomial,data=glm.data.train.3var,control = glm.control(maxit = 1000))
summary(glm.fit.3var)

# with interaction
glm.fit.3var <- glm(data.y~ENSG00000091831*ENSG00000082175*ENSG00000141736,family=binomial,data=glm.data.train.3var,control = glm.control(maxit = 1000))
summary(glm.fit.3var)

# predictions for the train and test sets
pred.train.glm.3var <- glm.fit.3var$fitted.values
pred.test.glm.3var <- predict(glm.fit.3var,as.data.frame(data.test[,c("ENSG00000091831","ENSG00000082175","ENSG00000141736")]),type="response")
```

Misclassifications in the training and test sets:
  
```{r}
# misclassified
glm.3var.miscl.train <- rownames(data.train[abs(data.Y.train - round(pred.train.glm.3var)) == 1,c("ENSG00000091831","ENSG00000082175","ENSG00000141736")])
length(glm.3var.miscl.train)

glm.3var.miscl.test <- rownames(data.test[abs(data.Y.test - round(pred.test.glm.3var)) == 1,c("ENSG00000091831","ENSG00000082175","ENSG00000141736")])
length(glm.3var.miscl.test)
```

45 misclassifications in the training set and 12 in the test set.

Dubious training cases which are misclassified:
  
```{r}
# misclassified and suspect in the train set
glm.3var.train.miscl.suspect <- glm.3var.miscl.train[which(getParticipantCode(glm.3var.miscl.train) %in% suspect)]
glm.3var.train.miscl.suspect

# misclassified and suspect in the test set
glm.3var.test.miscl.suspect <- glm.3var.miscl.test[which(getParticipantCode(glm.3var.miscl.test) %in% suspect)]
glm.3var.test.miscl.suspect
```

Several suspect cases with discordant HER2 (IHC) status and HER2 (FISH), as well as discordant HER2 (IHC) level and status are indeed influent and misclassified, based on the three biomarkers.


Looking for possible confounding variables in the clinical data:
  
```{r}
# categorical variables identified as potential interesting variables
clinical.interest <- c('age_at_initial_pathologic_diagnosis','race_list','history_of_neoadjuvant_treatment','gender','person_neoplasm_cancer_status','stage_event_pathologic_stage','stage_event_tnm_categories','menopause_status')
TNBC.glm.clinical <- clinical$primary.solid.tumor[, clinical.interest]

TNBC.glm.clinical <- data.frame(data.Y,TNBC.glm.clinical[rownames(clinical.tnbc), ])

# Univariate logistic models
TNBC.glm.fit <- glm(data.Y~age_at_initial_pathologic_diagnosis,family=binomial,data=TNBC.glm.clinical,control = glm.control(maxit = 1000))
summary(TNBC.glm.fit)# significant

TNBC.glm.fit <- glm(data.Y~gender,family=binomial,data=TNBC.glm.clinical,control = glm.control(maxit = 1000))
summary(TNBC.glm.fit)

TNBC.glm.fit <- glm(data.Y~race_list,family=binomial,data=TNBC.glm.clinical,control = glm.control(maxit = 1000))
summary(TNBC.glm.fit) # significant

TNBC.glm.fit <- glm(data.Y~history_of_neoadjuvant_treatment,family=binomial,data=TNBC.glm.clinical,control = glm.control(maxit = 1000))
summary(TNBC.glm.fit)

TNBC.glm.fit <- glm(data.Y~person_neoplasm_cancer_status,family=binomial,data=TNBC.glm.clinical,control = glm.control(maxit = 1000))

TNBC.glm.fit <- glm(data.Y~stage_event_pathologic_stage ,family=binomial,data=TNBC.glm.clinical,control = glm.control(maxit = 1000))

TNBC.glm.fit <- glm(data.Y~stage_event_tnm_categories ,family=binomial,data=TNBC.glm.clinical,control = glm.control(maxit = 1000))

TNBC.glm.fit <- glm(data.Y~menopause_status ,family=binomial,data=TNBC.glm.clinical,control = glm.control(maxit = 1000))
```


```{r}
# building a multivariate glm model accouting for all potential interesting categorical variables
TNBC.glm.fit <- glm(data.Y~age_at_initial_pathologic_diagnosis+gender+race_list+history_of_neoadjuvant_treatment+menopause_status+person_neoplasm_cancer_status+stage_event_pathologic_stage,family=binomial,data=TNBC.glm.clinical,control = glm.control(maxit = 1000))
#summary(TNBC.glm.fit)

# Chi-squared test between response and race (significant)
chisq.test(table(TNBC.glm.clinical$data.Y,TNBC.glm.clinical$race_list))

# Fisher test (more accurate - significant)
fisher.test(table(TNBC.glm.clinical$data.Y,TNBC.glm.clinical$race_list))


# building the final matrix (gene expression plus age and race)
data.plus.race.age <- data
names(TNBC.glm.clinical)
age.vect <- TNBC.glm.clinical[,2] 
race.vect <- as.character(TNBC.glm.clinical[,3])
race.vect[race.vect==""] <- ""
race.vect[race.vect=="AMERICAN INDIAN OR ALASKA NATIVE"] <- "0"
race.vect[race.vect=="WHITE"] <- "0"
race.vect[race.vect=="ASIAN"] <- "0"
race.vect[race.vect=="BLACK OR AFRICAN AMERICAN"] <- "1"
data.plus.race.age <- cbind(data.plus.race.age,as.factor(race.vect),age.vect,as.factor(data.Y))
dim(data.plus.race.age)

# changing binary code for race ("1" is "african american or black")
data.plus.race.age <- data.plus.race.age[data.plus.race.age[,19689]!="1",]# eliminates individuals without race
data.plus.race.age[data.plus.race.age[,19689]=="2",19689] <- 0
data.plus.race.age[data.plus.race.age[,19689]=="3",19689] <- 1

# changing binary code for Y (response)
data.plus.race.age[data.plus.race.age[,19691]=="1",19691] <- 0
data.plus.race.age[data.plus.race.age[,19691]=="2",19691] <- 1

summary(data.plus.race.age[,19690])
data.plus.race.age <- data.plus.race.age[-which(is.na(data.plus.race.age[,19690])),]
# eliminates individuals without age record

dim(data.plus.race.age)

data.plus.race.age.rank <- data.plus.race.age

# normalizing age
data.plus.race.age.rank[,19690] <- (data.plus.race.age.rank[,19690] - mean(data.plus.race.age.rank[,19690]))/ sd(data.plus.race.age.rank[,19690]) # normalized

data.Y.rank <- data.plus.race.age.rank[,19691]
length(data.Y.rank)

data.plus.race.age.rank <- data.plus.race.age.rank[,-19691]
dim(data.plus.race.age.rank)
```

# Ensemble outlier detection

## TNBC original data

```{r}
set.seed(2010)

library(glmnet)
my.alpha <- seq(0.1,0.9,0.1) 

## Optimizing alpha and lambda for logistic regression with Elastic net regularization

nvar.selected.EN <- matrix(0,1,length(my.alpha))
pred.EN <- matrix(0,dim(data.plus.race.age.rank)[1],length(my.alpha))
MSE.EN <- matrix(0,1,length(my.alpha))

# assigning samples to folds, to be used in cross-validation when tunning alpha
set.seed(2010)
foldid.tnbc <- sample(1:10,size=length(data.Y.rank),replace=TRUE)

for (j in 1:length(my.alpha)){
  
  # Logistic model fitting with 10-fold cross-validation for glmnet:
  fit.EN.cv <- cv.glmnet(as.matrix(data.plus.race.age.rank),as.factor(data.Y.rank),family="binomial",foldid=foldid.tnbc,alpha=my.alpha[j])
  
  var.selected.EN <- which(fit.EN.cv$glmnet.fit$beta[,which(fit.EN.cv$cvm == min(fit.EN.cv$cvm))] != 0)
  nvar.selected.EN[j] <- length(var.selected.EN)
  
  # Predictions obtained by model i
  pred.EN[,j] <- predict(fit.EN.cv,as.matrix(data.plus.race.age.rank),s="lambda.min",type="response")
  
  # Mean squared error of prediction (MSE)
  MSE.EN[j] <- mean((data.Y.rank-pred.EN)^2)
}

MSE.EN
nvar.selected.EN
my.alpha[which(MSE.EN == min(MSE.EN))]


## Model ensemble

n.models <- 3
nvar.selected <- matrix(0,1,n.models)
MSE <- matrix(0,1,n.models)
Miscl <- matrix(0,dim(data.plus.race.age.rank)[1],n.models)
CookD <- matrix(0,dim(data.plus.race.age.rank)[1],n.models)

# logistic regression with Elastic net regularization (LOGIT-EN)

i <- 1

my.alpha.opt <- my.alpha[which(MSE.EN == min(MSE.EN))]

fit.EN.cv <- cv.glmnet(as.matrix(data.plus.race.age.rank),as.factor(data.Y.rank),family="binomial",foldid=foldid.tnbc,alpha=my.alpha[which(MSE.EN == min(MSE.EN))])

var.selected.EN <- which(fit.EN.cv$glmnet.fit$beta[,which(fit.EN.cv$cvm == min(fit.EN.cv$cvm))] != 0)
nvar.selected[i] <- length(var.selected.EN)

assign(paste("var.selected", i, sep = ""), var.selected.EN)

# Predictions obtained by model i:
pred.EN <- predict(fit.EN.cv,as.matrix(data.plus.race.age.rank),s="lambda.min",type="response")

# Mean squared error of prediction (MSE):
MSE[i] <- mean((data.Y.rank-pred.EN)^2)

# Misclassified individuals based on model i:
Miscl[,i] <- abs(data.Y.rank - round(pred.EN))

table(data.Y.rank,round(pred.EN))

# The Cook’s distance for each individual based on model i:
V <- diag(as.vector(sqrt(pred.EN*(1-pred.EN))))
H <- V %*% data.plus.race.age.rank[,var.selected.EN] %*% (solve(t(data.plus.race.age.rank[,var.selected.EN]) %*% V %*% data.plus.race.age.rank[,var.selected.EN])) %*% t(data.plus.race.age.rank[,var.selected.EN]) %*% V # the hat matrix

CookD[,i] <- (data.Y.rank - pred.EN)^2 * diag(H) / ((pred.EN*(-pred.EN+1))*((-diag(H)+1))^2)


# SPLSDA

i <- 2

library(spls)

# Optimizing K and eta by cross-validation
#fit.splsda.LOGIT.cv <- cv.splsda(as.matrix(data.plus.race.age.rank),as.factor(data.Y.rank), fold=10, K = c(1:5), eta = c(0.9,0.8,0.7), kappa=0.5, classifier="logistic", scale.x=FALSE,n.core=10)

#fit.splsda.LOGIT.cv$K.opt
#fit.splsda.LOGIT.cv$eta.opt

# fixing the optimum K and eta obtained by cross-validation (for reproducibility)
fit.splsda.LOGIT <- splsda(as.matrix(data.plus.race.age.rank),as.factor(data.Y.rank), K = 4, eta=0.8, classifier="logistic", scale.x=FALSE)

var.selected.splsda.LOGIT <- fit.splsda.LOGIT$A
nvar.selected[i] <- length(var.selected.splsda.LOGIT)

assign(paste("var.selected", i, sep = ""), var.selected.splsda.LOGIT)

# Predictions obtained by model i:
pred.splsda.LOGIT <- predict(fit.splsda.LOGIT, as.matrix(data.plus.race.age.rank),type = "fit","coefficient", fit.type = "response")

table(data.Y.rank,round(pred.splsda.LOGIT))

# Mean squared error of prediction (MSE):
MSE[i] <- mean((data.Y.rank-pred.splsda.LOGIT)^2)

# Misclassified individuals based on model i:
Miscl[,i] <- abs(data.Y.rank - round(pred.splsda.LOGIT))

# The Cook’s distance for each individual based on model i:
V <- diag(as.vector(sqrt(pred.splsda.LOGIT*(1-pred.splsda.LOGIT))))
H <- V %*% (fit.splsda.LOGIT$T) %*% (solve(t(fit.splsda.LOGIT$T) %*% V %*% fit.splsda.LOGIT$T)) %*% t(fit.splsda.LOGIT$T) %*% V # the hat matrix

CookD[,i] <- (data.Y.rank - pred.splsda.LOGIT)^2 * diag(H) / ((pred.splsda.LOGIT*(-pred.splsda.LOGIT+1))*((-diag(H)+1))^2)


# SGPLS

i <- 3

library(spls)
library(parallel)

# Optimizing K and eta by cross-validation
#fit.sgpls.LOGIT.cv <- cv.sgpls(as.matrix(data.plus.race.age.rank),data.Y.rank, fold=10, K = c(1:5), eta = c(0.9,0.8,0.7), scale.x=FALSE,n.core=10)

# fixing the optimum K and eta obtained by cross-validation (for reproducibility)
fit.sgpls.LOGIT <- sgpls(as.matrix(data.plus.race.age.rank),data.Y.rank, K = 4, eta=0.7, scale.x=FALSE)

var.selected.sgpls.LOGIT <- fit.sgpls.LOGIT$A
nvar.selected[i] <- length(var.selected.sgpls.LOGIT)

assign(paste("var.selected", i, sep = ""), var.selected.sgpls.LOGIT)

# Predictions obtained by model i:
pred.sgpls.LOGIT <- predict(fit.sgpls.LOGIT, as.matrix(data.plus.race.age.rank),type="fit",fit.type="response")

table(data.Y.rank,round(pred.sgpls.LOGIT))

# Mean squared error of prediction (MSE):
MSE[i] <- mean((data.Y.rank-pred.sgpls.LOGIT)^2)

# Misclassified individuals based on model i:
Miscl[,i] <- abs(data.Y.rank - round(pred.sgpls.LOGIT))

# The Cook’s distance for each individual based on model i:
V <- diag(as.vector(sqrt(pred.sgpls.LOGIT*(1-pred.sgpls.LOGIT))))
H <- V %*% (data.plus.race.age.rank[,fit.sgpls.LOGIT$A]%*%fit.sgpls.LOGIT$W) %*% (solve(t(data.plus.race.age.rank[,fit.sgpls.LOGIT$A]%*%fit.sgpls.LOGIT$W) %*% V %*% data.plus.race.age.rank[,fit.sgpls.LOGIT$A]%*%fit.sgpls.LOGIT$W)) %*% t(data.plus.race.age.rank[,fit.sgpls.LOGIT$A]%*%fit.sgpls.LOGIT$W) %*% V 

CookD[,i] <- (data.Y.rank - pred.sgpls.LOGIT)^2 * diag(H) / ((pred.sgpls.LOGIT*(-pred.sgpls.LOGIT+1))*((-diag(H)+1))^2)

# Mean squared error for the three models:
MSE

# Misclassifications for the three models:
apply(Miscl,2,sum)

# Number of variables selected by the three models:
nvar.selected

# Variables selected by the three models:
var.selected=list(var.selected1,var.selected2,var.selected3)


# The rank product:
rank.matrix <- apply(CookD,2, function(CookD) rank(-(CookD),ties.method = "first"))
rank.product <- apply(rank.matrix, 1, prod)

rho <-rank.product
n <- dim(data.plus.race.age.rank)[1]
k <- dim(rank.matrix)[2]

source("~/Heskes_pvalues.R")

# The p-values:
pvalues <- as.vector(rankprodbounds(rho,n,k,Delta ='geometric'))

# The q-values:
library(qvalue)
qobj <- qvalue(pvalues)
qvalues <- qobj$qvalues

# Misclassifications (%) by the three models
Miscl.percent <- apply(Miscl,1, function(Miscl) round((sum(Miscl)/k)*100))

# suspect individuals within the misclassified by model
Miscl.suspect.LOGIT.EN <- rownames(data.plus.race.age.rank[which(Miscl[,1]==1),])[which(getParticipantCode(rownames(data.plus.race.age.rank[which(Miscl[,1]==1),])) %in% suspect)]
Miscl.suspect.LOGIT.EN

Miscl.suspect.SPLSDA <- rownames(data.plus.race.age.rank[which(Miscl[,2]==1),])[which(getParticipantCode(rownames(data.plus.race.age.rank[which(Miscl[,2]==1),])) %in% suspect)]
Miscl.suspect.SPLSDA

Miscl.suspect.SGPLS <- rownames(data.plus.race.age.rank[which(Miscl[,3]==1),])[which(getParticipantCode(rownames(data.plus.race.age.rank[which(Miscl[,3]==1),])) %in% suspect)]
Miscl.suspect.SGPLS

# Building the rank matrix
id <- as.vector(seq(1:dim(data.plus.race.age.rank)[1]))

outliers.rank <- as.data.frame(cbind(id, rank.matrix,rank.product,pvalues,qvalues,Miscl.percent))

# suspect in the top 25 logistic with EN model rank
rank.LOGIT.EN <- outliers.rank[order(outliers.rank[,2]),]
rownames(data.plus.race.age.rank[rank.LOGIT.EN[1:25,1],])[which(getParticipantCode(rownames(data.plus.race.age.rank[rank.LOGIT.EN[1:25,1],])) %in% suspect)]

# suspect in the top 25 SPLSDA model rank
rank.SPLSDA <- outliers.rank[order(outliers.rank[,3]),]
rownames(data.plus.race.age.rank[rank.SPLSDA[1:25,1],])[which(getParticipantCode(rownames(data.plus.race.age.rank[rank.SPLSDA[1:25,1],])) %in% suspect)]

# suspect in the top 25 SGPLS model rank
rank.SGPLS <- outliers.rank[order(outliers.rank[,4]),]
rownames(data.plus.race.age.rank[rank.SGPLS[1:25,1],])[which(getParticipantCode(rownames(data.plus.race.age.rank[rank.SGPLS[1:25,1],])) %in% suspect)]

outliers.rank <- outliers.rank[order(qvalues),]

# influent patients
TNBC.influent <- outliers.rank[which(outliers.rank[,7] < 0.05),1]
length(TNBC.influent)

# brca type of the influent patients
data.Y.rank[TNBC.influent]
sum(data.Y.rank[TNBC.influent]) # 10 TNBC and 14 non-TNBC patients

# influent patients which are suspect
TNBC.influent.suspect <- rownames(data.plus.race.age.rank[TNBC.influent,])[which(getParticipantCode(rownames(data.plus.race.age.rank[TNBC.influent,])) %in% suspect)]

# brca type of the influent, suspect patients
data.Y.rank[TNBC.influent][which(getParticipantCode(rownames(data.plus.race.age.rank[TNBC.influent,])) %in% suspect)]

# patients which are always misclassified
TNBC.alw.miscl <- which(apply(Miscl,1,sum)==3)
length(TNBC.alw.miscl)

# which are suspect
TNBC.miscl.suspect <- rownames(data.plus.race.age.rank[TNBC.alw.miscl,])[which(getParticipantCode(rownames(data.plus.race.age.rank[TNBC.alw.miscl,])) %in% suspect)]
length(TNBC.miscl.suspect)

# TNBC-associated variables (ER, PG and HER2) selected by the three models:
which(colnames(data.plus.race.age.rank)=="ENSG00000091831") %in% var.selected1
which(colnames(data.plus.race.age.rank)=="ENSG00000091831") %in% var.selected2
which(colnames(data.plus.race.age.rank)=="ENSG00000091831") %in% var.selected3
# ER selected by the three models

which(colnames(data.plus.race.age.rank)=="ENSG00000082175") %in% var.selected1
which(colnames(data.plus.race.age.rank)=="ENSG00000082175") %in% var.selected2
which(colnames(data.plus.race.age.rank)=="ENSG00000082175") %in% var.selected3
# PR only selected by SPLSDA

which(colnames(data.plus.race.age.rank)=="ENSG00000141736") %in% var.selected1
which(colnames(data.plus.race.age.rank)=="ENSG00000141736") %in% var.selected2
which(colnames(data.plus.race.age.rank)=="ENSG00000141736") %in% var.selected3
# HER2 selected by SPLSDA and SGPLS 

# Table 8 (Lopes et al.)
cbind(xdata[rownames(data.plus.race.age.rank[TNBC.influent,]),c("ENSG00000091831","ENSG00000082175","ENSG00000141736")],clinical.tnbc.aux[getParticipantCode(rownames(data.plus.race.age.rank[TNBC.influent,])),c(1:5)],clinical.tnbc[getParticipantCode(rownames(data.plus.race.age.rank[TNBC.influent,])),4],outliers.rank[1:length(TNBC.influent),2:8])


# plot SPLSDA-scores (Lopes et al.)
TNBC.matrix.splsda12.plot <- data.frame(LV1 = fit.splsda.LOGIT$T[,1],LV2 = fit.splsda.LOGIT$T[,2], Class = as.character(data.Y.rank))

TNBC.influent.matrix.splsda12.plot <- data.frame(LV1 =  fit.splsda.LOGIT$T[TNBC.influent,1],LV2 = fit.splsda.LOGIT$T[TNBC.influent,2], Class = as.character(data.Y.rank[TNBC.influent]))

TNBC.suspect.matrix.splsda12.plot <- data.frame(LV1 = fit.splsda.LOGIT$T[which(rownames(data.plus.race.age.rank)==TNBC.influent.suspect),1],LV2 = fit.splsda.LOGIT$T[which(rownames(data.plus.race.age.rank)==TNBC.influent.suspect),2], Class = as.character(TNBC.matrix.splsda12.plot[which(rownames(data.plus.race.age.rank)==TNBC.influent.suspect),3]))

library(ggplot2)
#postscript("splsda.ps")
#tiff("splsda.tiff")
ggplot(TNBC.matrix.splsda12.plot, aes(LV1, LV2)) + geom_point(aes(color = Class,shape = Class),size=3) + scale_shape_manual(values=c(21,24)) + scale_colour_manual(values=c("black", "black","black","black")) + theme(legend.text=element_text(size=2.5)) + theme_minimal() + theme(axis.text=element_text(size=14)) + theme(axis.title.y = element_text(size = 20, angle = 90)) + theme(axis.title.x = element_text(size = 20, angle = 0)) + 
  geom_point(data=TNBC.influent.matrix.splsda12.plot,aes(x=LV1, y=LV2,shape = Class,col="black",bg="red"),size=3) + 
  geom_point(data=TNBC.suspect.matrix.splsda12.plot,aes(x=LV1, y=LV2,shape = Class,col="black",bg="blue"),size=3) + guides(colour=FALSE)
#dev.off()

# PCA on TNBC data
ytnbc.col <- data.Y.rank
ytnbc.col[ytnbc.col==1] <- "black"
ytnbc.col[ytnbc.col==0] <- "black"
ytnbc.col[TNBC.influent] <- "blue"
ytnbc.col[which(rownames(data.plus.race.age.rank)==TNBC.influent.suspect)] <- "red"
 
ytnbc.shape <- data.Y.rank
ytnbc.shape[ytnbc.shape==1] <- 2
ytnbc.shape[ytnbc.shape==0] <- 1

ytnbc.shape <- as.factor(ytnbc.shape)
ytnbc.col <- as.factor(ytnbc.col)

tnbc.pca.matrix <- as.data.frame(cbind(data.plus.race.age.rank,ytnbc.shape,ytnbc.col))

tnbc.pca <- prcomp(tnbc.pca.matrix[,1:19690], center=FALSE,scale=FALSE)

library(ggfortify)
#postscript("pca.ps")
#tiff("pca.tiff")
autoplot(tnbc.pca, data = tnbc.pca.matrix, colour = ytnbc.col, shape= ytnbc.shape, size=3)  + theme_minimal()
#dev.off()

```

Evaluating the ensemble performance for synthetic datasets built based on mean and covariance matrices from the TNBC and non-TNBC classes

```{r}
n.synthetic <- 1000

# taking the covariance from the original TNBC data, to make simulation more realistic
xdata.synthetic.cov <- cov(data.plus.race.age.rank)
xdata.synthetic.cov.event <- cov(data.plus.race.age.rank[which(data.Y.rank==1),])
xdata.synthetic.cov.noevent <- cov(data.plus.race.age.rank[which(data.Y.rank==0),])

# for p = 5000 randomly chosen variables
p.synthetic <- 5000
set.seed(7)
cov.vars <- sample(c(1:dim(xdata.synthetic.cov)[2]),p.synthetic)
xdata.synthetic.cov.event5000 <- xdata.synthetic.cov.event[cov.vars,cov.vars]
xdata.synthetic.cov.noevent5000 <- xdata.synthetic.cov.noevent[cov.vars,cov.vars]

# for p = full TNBC dimensionality
p.synthetic <- dim(data.plus.race.age.rank[,1:19688])[2]
xdata.synthetic.cov.eventfull <- xdata.synthetic.cov.event[1:p.synthetic,1:p.synthetic]
xdata.synthetic.cov.noeventfull <- xdata.synthetic.cov.noevent[1:p.synthetic,1:p.synthetic]

# taking the covariance from the original TNBC data, to make simulation more realistic
# for p = 5000 randomly chosen variables
mean.synthetic.event5000 <- apply(data.plus.race.age.rank[which(data.Y.rank==1),],2,mean)[cov.vars]
mean.synthetic.noevent5000 <- apply(data.plus.race.age.rank[which(data.Y.rank==0),],2,mean)[cov.vars]

# for p = full TNBC dimensionality
mean.synthetic.eventfull <- apply(data.plus.race.age.rank[which(data.Y.rank==1),],2,mean)[1:p.synthetic]
mean.synthetic.noeventfull <- apply(data.plus.race.age.rank[which(data.Y.rank==0),],2,mean)[1:p.synthetic]

# generate  a random X
library(MASS)
set.seed(8)
xdata.synthetic.event5000 <- mvrnorm(n = 500, mean.synthetic.event5000, xdata.synthetic.cov.event5000, tol = 1e-6, empirical = FALSE)
set.seed(9)
xdata.synthetic.noevent5000 <- mvrnorm(n = 500, mean.synthetic.noevent5000, xdata.synthetic.cov.noevent5000, tol = 1e-6, empirical = FALSE)

xdata.synthetic <- rbind(xdata.synthetic.event5000,xdata.synthetic.noevent5000)
ydata.synthetic <- c(rep(1,500), rep(0,500))
data.Y.synthetic.out <- ydata.synthetic
set.seed(9)
data.synthetic.outliers <- sample(c(1:1000),25, replace=FALSE) # simulating 25 outliers
data.synthetic.outliers 
data.Y.synthetic.out[data.synthetic.outliers] <- abs(1-data.Y.synthetic.out[data.synthetic.outliers])

data.synthetic <- xdata.synthetic
# normalized X data
data.synthetic <- (data.synthetic-matrix(apply(data.synthetic,2,mean),nrow(data.synthetic),ncol(data.synthetic),byrow=TRUE))/matrix(apply(data.synthetic,2,sd),nrow(data.synthetic),ncol(data.synthetic),byrow=TRUE)

# Plots

#postscript("Fig_synthetic_outliers_PCA_n1000_p5000_muTNBC_covTNBC_out25.ps")
#tiff("Fig_synthetic_outliers_PCA_n1000_p5000_muTNBC_covTNBC_out25.tiff")
ysynthetic.col <- data.Y.synthetic.out
ysynthetic.col[ysynthetic.col==1] <- "black"
ysynthetic.col[ysynthetic.col==0] <- "gray"
ysynthetic.col[data.synthetic.outliers] <- "red"
 
synthetic.pca.matrix <- as.data.frame(data.synthetic)
synthetic.pca.matrix.groups <- factor(ysynthetic.col)
synthetic.pca.matrix <- as.matrix(synthetic.pca.matrix)

synthetic.pca <- prcomp(synthetic.pca.matrix, center=FALSE,scale=FALSE)
 
autoplot(synthetic.pca, data = synthetic.pca.matrix, colour = synthetic.pca.matrix.groups)
#dev.off()


# Ensemble outlier detection

# e.g., for p=5000 and % of outliers = 2.5

set.seed(2010)

my.alpha.synthetic <- seq(0.1,0.9,0.1)

# Optimizing alpha and lambda for logistic regression with Elastic net regularization
nvar.selected.EN.synthetic <- matrix(0,1,length(my.alpha.synthetic) )
MSE.EN.synthetic <- matrix(0,1,length(my.alpha.synthetic))
MSE.EN.synthetic <- matrix(0,1,length(my.alpha.synthetic))

# assigning samples to folds, to be used in cross-validation when tunning alpha
set.seed(2010)
foldid.synthetic <-sample(1:10,size=length(data.Y.synthetic.out),replace=TRUE)

library(glmnet)

for (j in 1:length(my.alpha.synthetic)){
  
  # Logistic model fitting with 10-fold cross-validation for glmnet:
  fit.EN.cv.synthetic <- cv.glmnet(as.matrix(data.synthetic),as.factor(data.Y.synthetic.out),family="binomial",foldid=foldid.synthetic,alpha=my.alpha.synthetic[j])
  var.selected.EN.synthetic <- which(fit.EN.cv.synthetic$glmnet.fit$beta[,which(fit.EN.cv.synthetic$cvm == min(fit.EN.cv.synthetic$cvm))] != 0)
  nvar.selected.EN.synthetic[j] <- length(var.selected.EN.synthetic)
  
  # Predictions obtained by model i:
  pred.EN.cv.synthetic <- predict(fit.EN.cv.synthetic,as.matrix(data.synthetic),s="lambda.min",type="response")
  
  # Mean squared error of prediction (MSE):
  MSE.EN.synthetic[j] <- mean((data.Y.synthetic.out-pred.EN.cv.synthetic)^2)
}

MSE.EN.synthetic
nvar.selected.EN.synthetic
my.alpha.synthetic[which(MSE.EN.synthetic == min(MSE.EN.synthetic))]


## Model ensemble

n.models.synthetic <- 3
nvar.selected.synthetic <- matrix(0,1,n.models.synthetic)
MSE.synthetic <- matrix(0,1,n.models.synthetic)
Miscl.synthetic <- matrix(0,dim(data.synthetic)[1],n.models.synthetic)
CookD.synthetic <- matrix(0,dim(data.synthetic)[1],n.models.synthetic)

# logistic regression with Elastic net regularization (LOGIT-EN)

i <- 1

my.alpha.opt.synthetic <- my.alpha.synthetic[which(MSE.EN.synthetic == min(MSE.EN.synthetic))]

fit.EN.cv.synthetic <- cv.glmnet(as.matrix(data.synthetic),as.factor(data.Y.synthetic.out),family="binomial",foldid=foldid.synthetic,alpha=my.alpha.opt.synthetic)    
var.selected.EN.synthetic <- which(fit.EN.cv.synthetic$glmnet.fit$beta[,which(fit.EN.cv.synthetic$cvm == min(fit.EN.cv.synthetic$cvm))] != 0)
nvar.selected.synthetic[i] <- length(var.selected.EN.synthetic)

assign(paste("var.selected.synthetic", i, sep = ""), which(fit.EN.cv.synthetic$glmnet.fit$beta[,which(fit.EN.cv.synthetic$cvm == min(fit.EN.cv.synthetic$cvm))] != 0))

# Predictions obtained by model i:
pred.EN.cv.synthetic <- predict(fit.EN.cv.synthetic,as.matrix(data.synthetic),s="lambda.min",type="response")

# Mean squared error of prediction (MSE)
MSE.synthetic[i] <- mean((data.Y.synthetic.out-pred.EN.cv.synthetic)^2)

# Misclassified individuals (‘1’) based on method i:
Miscl.synthetic[,i] <- abs(data.Y.synthetic.out - round(pred.EN.cv.synthetic))

table(data.Y.synthetic.out,round(pred.EN.cv.synthetic))

# The Cook’s distance for each individual based on method i:
V <- diag(as.vector(sqrt(pred.EN.cv.synthetic*(1-pred.EN.cv.synthetic))))
H <- V %*% data.synthetic[,var.selected.EN.synthetic] %*% (solve(t(data.synthetic[,var.selected.EN.synthetic]) %*% V %*% data.synthetic[,var.selected.EN.synthetic])) %*% t(data.synthetic[,var.selected.EN.synthetic]) %*% V 

CookD.synthetic[,i] <- (data.Y.synthetic.out - pred.EN.cv.synthetic)^2 * diag(H) / ((pred.EN.cv.synthetic*(-pred.EN.cv.synthetic+1))*((-diag(H)+1))^2)

which(Miscl.synthetic[,i]==1)
which(data.synthetic.outliers %in%  which(Miscl.synthetic[,i]==1))

# SPLS-DA

library(spls)

i <- 2

# Optimizing the parameters by cross-validation
#fit.splsda.synthetic.cv <- cv.splsda(as.matrix(data.synthetic),as.factor(data.Y.synthetic.out), fold=10, K = c(1:5), eta = c(0.9,0.8,0.7), kappa=0.5, classifier="logistic", scale.x=FALSE)

# using the parameters optimized above
fit.splsda.synthetic <- splsda(as.matrix(data.synthetic),as.factor(data.Y.synthetic.out), K = 5, eta=0.7, classifier="logistic", scale.x=FALSE)

var.selected.splsda.synthetic <- fit.splsda.synthetic$A
nvar.selected.synthetic[i] <- length(var.selected.splsda.synthetic)

assign(paste("var.selected.synthetic", i, sep = ""), var.selected.splsda.synthetic)

# Predictions obtained by model i:
pred.splsda.synthetic<- predict(fit.splsda.synthetic, as.matrix(data.synthetic),type = "fit","coefficient", fit.type = "response")

# Mean squared error of prediction (MSE)
MSE.synthetic[i] <- mean((data.Y.synthetic.out-pred.splsda.synthetic)^2)

# Misclassified individuals (‘1’) based on method i:
Miscl.synthetic[,i] <- abs(data.Y.synthetic.out - round(pred.splsda.synthetic))

table(data.Y.synthetic.out,round(pred.splsda.synthetic))

# The Cook’s distance for each individual based on method i:
V <- diag(as.vector(sqrt(pred.splsda.synthetic*(1-pred.splsda.synthetic))))
H <- V %*% (fit.splsda.synthetic$T) %*% (solve(t(fit.splsda.synthetic$T) %*% V %*% fit.splsda.synthetic$T)) %*% t(fit.splsda.synthetic$T) %*% V # the hat matrix

CookD.synthetic[,i] <- (data.Y.synthetic.out - pred.splsda.synthetic)^2 * diag(H) / ((pred.splsda.synthetic*(-pred.splsda.synthetic+1))*((-diag(H)+1))^2)

which(Miscl.synthetic[,i]==1)
which(data.synthetic.outliers %in%  which(Miscl.synthetic[,i]==1))

# SGPLS

i <- 3

# Optimizing the parameters by cross-validation
#fit.sgpls.synthetic.cv <- cv.sgpls(as.matrix(data.synthetic),data.Y.synthetic.out, fold=10, K = c(1:5), eta = c(0.9,0.8,0.7), scale.x=FALSE)

# using the parameters optimized above
fit.sgpls.synthetic <- sgpls(as.matrix(data.synthetic),data.Y.synthetic.out, K = 4, eta=0.8, scale.x=FALSE)

var.selected.sgpls.synthetic <- fit.sgpls.synthetic$A
nvar.selected.synthetic[i] <- length(var.selected.sgpls.synthetic)

assign(paste("var.selected.synthetic", i, sep = ""), var.selected.sgpls.synthetic)

# Predictions obtained for the train set:
pred.sgpls.synthetic <- predict(fit.sgpls.synthetic,as.matrix(data.synthetic),type="fit",fit.type="response")

# Mean squared error of prediction (MSE)
MSE.synthetic[i] <- mean((data.Y.synthetic.out-pred.sgpls.synthetic)^2)

# Misclassified individuals (‘1’) based on method i:
Miscl.synthetic[,i] <- abs(data.Y.synthetic.out - round(pred.sgpls.synthetic))

table(data.Y.synthetic.out,round(pred.sgpls.synthetic))

# The Cook’s distance for each individual based on method i:
V <- diag(as.vector(sqrt(pred.sgpls.synthetic*(1-pred.sgpls.synthetic))))
H <- V %*% (data.synthetic[,var.selected.sgpls.synthetic]%*%fit.sgpls.synthetic$W) %*% (solve(t(data.synthetic[,var.selected.sgpls.synthetic]%*%fit.sgpls.synthetic$W) %*% V %*% data.synthetic[,var.selected.sgpls.synthetic]%*%fit.sgpls.synthetic$W)) %*% t(data.synthetic[,var.selected.sgpls.synthetic]%*%fit.sgpls.synthetic$W) %*% V # the hat matrix

CookD.synthetic[,i] <- (data.Y.synthetic.out - pred.sgpls.synthetic)^2 * diag(H) / ((pred.sgpls.synthetic*(-pred.sgpls.synthetic+1))*((-diag(H)+1))^2)

which(Miscl.synthetic[,i]==1)
which(data.synthetic.outliers %in%  which(Miscl.synthetic[,i]==1))

nvar.selected.synthetic
var.selected.synthetic=list(var.selected.synthetic1,var.selected.synthetic2,var.selected.synthetic3)

source("~/Heskes_pvalues.R")

# The rank product:
rank.matrix.synthetic <- apply(CookD.synthetic,2, function(CookD.synthetic) rank(-(CookD.synthetic),ties.method = "first"))
rank.product.synthetic <- apply(rank.matrix.synthetic, 1, prod)

rho.synthetic <-rank.product.synthetic
n.synthetic <- dim(data.synthetic)[1]
k.synthetic <- dim(rank.matrix.synthetic)[2]

# The p-values:
pvalues.synthetic <- as.vector(rankprodbounds(rho.synthetic,n.synthetic,k.synthetic,Delta ='geometric'))

library(qvalue)
# The q-values:
qobj.synthetic <- qvalue(pvalues.synthetic)
qvalues.synthetic <- qobj.synthetic$qvalues

# Misclassifications
Miscl.train.percent.synthetic <- apply(Miscl.synthetic,1, function(Miscl.synthetic) round((sum(Miscl.synthetic)/k.synthetic)*100))
Miscl.train.by.method.synthetic <- apply(Miscl.synthetic,2, sum)

# Building the rank matrix
id.synthetic <- as.vector(seq(1:n.synthetic))

outliers.rank.synthetic <-as.data.frame(cbind(id.synthetic, rank.matrix.synthetic,rank.product.synthetic,pvalues.synthetic,qvalues.synthetic,Miscl.train.percent.synthetic))
outliers.rank.synthetic <- outliers.rank.synthetic[order(qvalues.synthetic),]
outliers.rank.synthetic[1:25,]

synthetic.influent <- outliers.rank.synthetic[which(outliers.rank.synthetic[,7] < 0.05),1]

which(data.synthetic.outliers %in% synthetic.influent)
length(which(data.synthetic.outliers %in% synthetic.influent))

length(which(outliers.rank.synthetic[1:length(data.synthetic.outliers),1] %in% data.synthetic.outliers))
length(which(outliers.rank.synthetic[which(outliers.rank.synthetic[,2]<=length(data.synthetic.outliers)),1] %in% data.synthetic.outliers))
length(which(outliers.rank.synthetic[which(outliers.rank.synthetic[,3]<=length(data.synthetic.outliers)),1] %in% data.synthetic.outliers))
length(which(outliers.rank.synthetic[which(outliers.rank.synthetic[,4]<=length(data.synthetic.outliers)),1] %in% data.synthetic.outliers))

```

## Genes selected

Genes selected by LOGIT-EN, SPLS-DA and SGPLS

```{r}
# variables in common between EN and sPLSDA and sgPLS
TNBC.vars.in.common <- as.data.frame(table(c(var.selected1,var.selected2,var.selected3)))
TNBC.vars.in.common <- as.numeric(as.character(TNBC.vars.in.common[which(TNBC.vars.in.common[,2] == 3),1]))
length(TNBC.vars.in.common)
# 26 vars in common

TNBC.vars.in.common.code <- as.matrix(mcols(ensembl.protein.coding[colnames(data.plus.race.age.rank[,TNBC.vars.in.common]),2]))
TNBC.vars.in.common.code

# summary of gene expression in TNBC and non-TNBC for the 26 genes selected
TNBC.vars.in.common.expression.TNBC <- xdata[ydata[,1]=="TNBC",TNBC.vars.in.common]
colnames(TNBC.vars.in.common.expression.TNBC) <- as.vector(TNBC.vars.in.common.code)
summary(TNBC.vars.in.common.expression.TNBC)

TNBC.vars.in.common.expression.nonTNBC <- xdata[ydata[,1]=="NO_TNBC",TNBC.vars.in.common]
colnames(TNBC.vars.in.common.expression.nonTNBC) <- as.vector(TNBC.vars.in.common.code)
summary(TNBC.vars.in.common.expression.nonTNBC)


# ranking coefficients of the variables selected in common by the three methods

EN.coef <- abs(fit.EN.cv$glmnet.fit$beta[,which(fit.EN.cv$cvm == min(fit.EN.cv$cvm))])
EN.coef.matrix <- cbind(c(1:dim(data.plus.race.age.rank)[2]), EN.coef)
rank.EN.coef <- cbind(EN.coef.matrix[order(EN.coef.matrix[,2], decreasing = TRUE),],c(1:dim(data.plus.race.age.rank)[2]))
rank.EN.coef <- rank.EN.coef[which(rank.EN.coef[,1] %in% TNBC.vars.in.common),]

splsda.coef <- abs(coef(fit.splsda.LOGIT))[-1]
splsda.coef.matrix <- cbind(c(1: dim(data.plus.race.age.rank)[2]), splsda.coef)
rank.splsda.coef <- cbind(splsda.coef.matrix[order(splsda.coef.matrix[,2], decreasing = TRUE),],c(1:dim(data.plus.race.age.rank)[2]))
rank.splsda.coef <- rank.splsda.coef[which(rank.splsda.coef[,1] %in% TNBC.vars.in.common),]

sgpls.coef <- abs(coef(fit.sgpls.LOGIT))[-1]
sgpls.coef.matrix <- cbind(c(1: dim(data.plus.race.age.rank)[2]), sgpls.coef)
rank.sgpls.coef <- cbind(sgpls.coef.matrix[order(sgpls.coef.matrix[,2] , decreasing = TRUE),],c(1:dim(data.plus.race.age.rank)[2]))
rank.sgpls.coef <- rank.sgpls.coef[which(rank.sgpls.coef[,1] %in% TNBC.vars.in.common),]

variables.rank.matrix <- cbind(sort(rank.EN.coef[,1]),rank.EN.coef[order(rank.EN.coef[,1]),3], rank.splsda.coef[order(rank.splsda.coef[,1]),3], rank.sgpls.coef[order(rank.sgpls.coef[,1]),3])
colnames(variables.rank.matrix) <- c("variable","LOGIT-EN","SPLSDA","SGPLS")
#variables.rank.matrix <- cbind(as.matrix(mcols(ensembl.protein.coding[rownames(variables.rank.matrix),2])),variables.rank.matrix)

combined.variables.rank.matrix <- cbind(TNBC.vars.in.common.code,apply(variables.rank.matrix[,2:4],1,prod))
combined.variables.rank.matrix

combined.variables.rank.matrix[order(as.numeric(combined.variables.rank.matrix[,2])),1]
```

Genes selected by LOGIT-EN, SPLS-DA and SGPLS when changing the label of the influential observations misclassified in more than 2 methods 

```{r}
data.Y.rank.change <- data.Y.rank
data.Y.rank.change[outliers.rank[which(outliers.rank[1:24,8] >= 67),1]
] <- 1-data.Y.rank.change[outliers.rank[which(outliers.rank[1:24,8] >= 67),1]]

# LOGIT-EN

my.alpha <- seq(0.1,0.9,0.1) 

## Optimizing alpha and lambda for logistic regression with Elastic net regularization

nvar.selected.EN.outliers.change <- matrix(0,1,length(my.alpha))
pred.EN.outliers.change <- matrix(0,dim(data.plus.race.age.rank)[1],length(my.alpha))
MSE.EN.outliers.change <- matrix(0,1,length(my.alpha))

set.seed(2010)

for (j in 1:length(my.alpha)){
  
  # Logistic model fitting with 10-fold cross-validation for glmnet:
  fit.EN.cv.outliers.change <- cv.glmnet(as.matrix(data.plus.race.age.rank),as.factor(data.Y.rank.change),family="binomial",foldid=foldid.tnbc,alpha=my.alpha[j])
  
  var.selected.EN.outliers.change <- which(fit.EN.cv.outliers.change$glmnet.fit$beta[,which(fit.EN.cv.outliers.change$cvm == min(fit.EN.cv.outliers.change$cvm))] != 0)
  nvar.selected.EN.outliers.change[j] <- length(var.selected.EN.outliers.change)
  
  # Predictions obtained by model i
  pred.EN.outliers.change[,j] <- predict(fit.EN.cv.outliers.change,as.matrix(data.plus.race.age.rank),s="lambda.min",type="response")
  
  # Mean squared error of prediction (MSE)
  MSE.EN.outliers.change[j] <- mean((data.Y.rank.change-pred.EN.outliers.change)^2)
}

MSE.EN.outliers.change
nvar.selected.EN.outliers.change
my.alpha[which(MSE.EN.outliers.change == min(MSE.EN.outliers.change))]

my.alpha.opt <- my.alpha[which(MSE.EN.outliers.change == min(MSE.EN.outliers.change))]

fit.EN.cv.outliers.change <- cv.glmnet(as.matrix(data.plus.race.age.rank),as.factor(data.Y.rank.change),family="binomial",foldid=foldid.tnbc,alpha=my.alpha.opt)
var.selected.EN.outliers.change <- which(fit.EN.cv.outliers.change$glmnet.fit$beta[,which(fit.EN.cv.outliers.change$cvm == min(fit.EN.cv.outliers.change$cvm))] != 0)

# Predictions
pred.EN.outliers.change <- predict(fit.EN.cv.outliers.change,as.matrix(data.plus.race.age.rank),s="lambda.min",type="response")

# Mean squared error of prediction (MSE):
MSE.EN.outliers.change <- mean((data.Y.rank.change-pred.EN.outliers.change)^2)
MSE.EN.outliers.change

# SPLSDA

#fit.splsda.LOGIT.cv.outliers.change <- cv.splsda(as.matrix(data.plus.race.age.rank),as.factor(data.Y.rank.change), fold=10, K = c(1:5), eta = c(0.9,0.8,0.7), kappa=0.5, classifier="logistic", scale.x=FALSE,n.core=10)
fit.splsda.LOGIT.outliers.change <- splsda(as.matrix(data.plus.race.age.rank),as.factor(data.Y.rank.change), K = 5, eta=0.7, classifier="logistic", scale.x=FALSE) # parameters K and eta obtained by cross-validation above
var.selected.splsda.LOGIT.outliers.change <- fit.splsda.LOGIT.outliers.change$A

# Predictions
pred.splsda.LOGIT.outliers.change <- predict(fit.splsda.LOGIT.outliers.change, as.matrix(data.plus.race.age.rank),type = "fit","coefficient", fit.type = "response")

# Mean squared error of prediction (MSE):
MSE.splsda.LOGIT.outliers.change <- mean((data.Y.rank.change-pred.splsda.LOGIT.outliers.change)^2)
MSE.splsda.LOGIT.outliers.change

# SGPLS

#fit.sgpls.LOGIT.cv.outliers.out <- cv.sgpls(as.matrix(data.plus.race.age.rank),data.Y.rank.change, fold=10, K = c(1:5), eta = c(0.9,0.8,0.7), scale.x=FALSE,n.core=10)
fit.sgpls.LOGIT.outliers.change <- sgpls(as.matrix(data.plus.race.age.rank),data.Y.rank.change, K = 5, eta=0.8, scale.x=FALSE)
var.selected.sgpls.LOGIT.outliers.change <- fit.sgpls.LOGIT.outliers.change$A

# Predictions
pred.sgpls.LOGIT.outliers.change <- predict(fit.sgpls.LOGIT.outliers.change, as.matrix(data.plus.race.age.rank),type="fit",fit.type="response")

# Mean squared error of prediction (MSE):
MSE.sgpls.LOGIT.outliers.change <- mean((data.Y.rank.change-pred.sgpls.LOGIT.outliers.change)^2)
MSE.sgpls.LOGIT.outliers.change


# variables in common between EN and SPLSDA and SGPLS
outliers.change.vars.in.common <- as.data.frame(table(c(var.selected.EN.outliers.change,var.selected.splsda.LOGIT.outliers.change,var.selected.sgpls.LOGIT.outliers.change)))
outliers.change.vars.in.common <- as.numeric(as.character(outliers.change.vars.in.common[which(outliers.change.vars.in.common[,2] == 3),1]))
length(outliers.change.vars.in.common)

outliers.change.vars.in.common.code <- as.matrix(mcols(ensembl.protein.coding[colnames(data.plus.race.age.rank[,outliers.change.vars.in.common]),2]))
outliers.change.vars.in.common.code

# variables in common with the previous 26 relevant selected by the three methods
outliers.change.vars.in.common.TNBC <- outliers.change.vars.in.common[which(outliers.change.vars.in.common %in% TNBC.vars.in.common)]
outliers.change.vars.in.common.TNBC.code <- as.matrix(mcols(ensembl.protein.coding[colnames(data.plus.race.age.rank[,outliers.change.vars.in.common.TNBC]),2]))
outliers.change.vars.in.common.TNBC.code 

# summary of gene expression in TNBC and non-TNBC for the 10 genes selected

dim(xdata)

outliers.change.vars.in.common.name <- colnames(xdata[,outliers.change.vars.in.common])

outliers.change.vars.in.common.expression.TNBC <- xdata[data.Y=="1",outliers.change.vars.in.common.name]
colnames(outliers.change.vars.in.common.expression.TNBC) <- as.vector(outliers.change.vars.in.common.code)
summary(outliers.change.vars.in.common.expression.TNBC)

outliers.change.vars.in.common.expression.nonTNBC <- xdata[data.Y=="0",outliers.change.vars.in.common.name]
colnames(outliers.change.vars.in.common.expression.nonTNBC) <- as.vector(outliers.change.vars.in.common.code)
summary(outliers.change.vars.in.common.expression.nonTNBC)
```

## TNBC data with varying sets of patients

The different models used for ensemble analysis are sparse logistic model with alpha = 0.7 (EN) based on randonly selected 80 % of patients.

```{r}
source("~/ensembl.rand.patients.R")
library('qvalue')

set.seed(2010)

n.models.rand <- 100
outliers.rand.patients.EN <- ensembl.rand.patients(data.plus.race.age.rank,data.Y.rank,n.models.rand)

# mean number of variables selected in the 100 runs
mean(outliers.rand.patients.EN$nvariables.selected)
# mean MSE for the 100 runs
mean(outliers.rand.patients.EN$MSE)
# mean number of misclassification in the 100 runs
mean(outliers.rand.patients.EN$Misclassifications.by.method)

# influent patients
TNBC.rand.patient.influent <- outliers.rand.patients.EN$outliers.rank[which(outliers.rand.patients.EN$outliers.rank[,n.models.rand+4] < 0.05),1]
length(TNBC.rand.patient.influent)

# influent patients which are suspect
rownames(data.plus.race.age.rank[TNBC.rand.patient.influent,])[which(getParticipantCode(rownames(data.plus.race.age.rank[TNBC.rand.patient.influent,])) %in% suspect)]

TNBC.rand.patient.var.select <- as.data.frame(table(outliers.rand.patients.EN$var.selected))
TNBC.rand.patient.var.select$Var1 <- as.numeric(as.character(TNBC.rand.patient.var.select[,1]))
TNBC.rand.patient.var.select$var.name <- colnames(data.plus.race.age.rank[,TNBC.rand.patient.var.select[,1]])

# variabes selected in more that 75% of the runs
TNBC.rand.patient.var.select.more75 <- TNBC.rand.patient.var.select[which(TNBC.rand.patient.var.select[,2]>75),3]
length(TNBC.rand.patient.var.select.more75)

# variables selected in common by the ensembles for the original TNBC and for the 'random patient' strategy
TNBC.rand.patient.var.select.common.originalTNBC <- TNBC.rand.patient.var.select.more75[which(TNBC.rand.patient.var.select.more75 %in% colnames(data.plus.race.age.rank[,TNBC.vars.in.common]))]
length(TNBC.rand.patient.var.select.common.originalTNBC)

# variables code
ensembl.protein.coding[TNBC.rand.patient.var.select.common.originalTNBC,2][,1]
as.matrix(mcols(ensembl.protein.coding[TNBC.rand.patient.var.select.common.originalTNBC,2]))
```

## TNBC data with varying sets of variables

The different models used for ensemble analysis are sparse logistic model with alpha = 0.7 (EN), departing from a set of 1000 randomly selected variables.

```{r}
source("~/ensembl.rand.variables.R")
times <- 100
vars.rand <- 1000
outliers.rand.variables.EN <- ensembl.rand.variables(data.plus.race.age.rank,data.Y.rank, times, vars.rand)

# mean number of variables selected in the 100 runs
mean(outliers.rand.variables.EN$nvariables.selected)
# mean MSE for the 100 runs
mean(outliers.rand.variables.EN$MSE)
# mean number of misclassification in the 100 runs
mean(outliers.rand.variables.EN$Misclassifications.by.method)

# influent patients
TNBC.rand.variables.influent <- outliers.rand.variables.EN$outliers.rank[which(outliers.rand.variables.EN$outliers.rank[,times+4] < 0.05),1]
length(TNBC.rand.variables.influent)

# influent patients which are suspect
rownames(data.plus.race.age.rank[TNBC.rand.variables.influent,])[which(getParticipantCode(rownames(data.plus.race.age.rank[TNBC.rand.variables.influent,])) %in% suspect)]

# which influent patients in common with the ensemble for the original TNBC data and by the 'random patient' and 'random variable' strategies
TNBC.rand.outliers.in.common <- as.data.frame(table(c(TNBC.rand.variables.influent,TNBC.rand.patient.influent,TNBC.influent)))
TNBC.rand.outliers.in.common <- as.numeric(as.character(TNBC.rand.outliers.in.common[which(TNBC.rand.outliers.in.common[,2] == 3),1]))
length(TNBC.rand.outliers.in.common)

getParticipantCode(rownames(data.plus.race.age.rank[TNBC.rand.outliers.in.common,]))
# 18 patients in common (TNBC original + rand patients + rand variables) detected as outliers 

# influent patients in common for the three emsemble strategies which are suspect
rownames(data.plus.race.age.rank[TNBC.rand.outliers.in.common,])[which(getParticipantCode(rownames(data.plus.race.age.rank[TNBC.rand.outliers.in.common,])) %in% suspect)]
```
