---
title: "Outlier detection, RP"
author: "IDMEC Team 2017"
date: "Dataset -Ovarian"
output:
  html_document:
    number_sections: yes
    toc: yes
---


# Description of the dataset
The ovarian cancer dataset is based on gene expression data of oncological patients and is constituted by 517 observations over 12042 covariates. This data was obtained from The Cancer Genome Atlas (TCGA) (http://cancergenome.nih.gov/).

The dataset is publicly available (https://gdc-portal.nci.nih.gov/) and it was normalized and aggregated by the TCGA consortium allowing for the analysis to be reproducible with the original dataset. We only cleaned corrupted case records from the data, where correlated attributes had invalid information.

The clinical data was cleaned using "Days to last follow-up" and "Days to death" attributes to detect inconsistencies between them. Only the cases where the number of days matched were included in the analysis. The same process was performed for the attributes "Days to death" and "Vital status", where some cases had as status "deceased", but a missing "Days to death".

This dataset was analysed in three different ways. 

Lets read each dataset. Do not forget to change the filepath for the correct directory of your computer.

The covariate matrix has to be defined. The data is stored in a txt file.


# TCGA ovarian cancer 18 genes

```{r}
data.survival18<- read.table("./tcga18.txt",dec = ",", header = TRUE)

head(data.survival18)
```


And the matrix for the time and status is given by:

```{r}

data.time.status<-read.table("./timestatus.txt",dec = ",", header = TRUE)
head(data.time.status)
```



## Survival analysis
The Cox's proportional hazard was performed in all variables. 

### Cox Regression model

Let's start the analysis for the Cox regression model.
```{r}
library(survival)
```

Fit the Cox proportional hazards to the dataset. 

```{r}
fit.cox.all18<-coxph(Surv(data.time.status$time,data.time.status$status) ~ .,data=data.survival18)
fit.cox.all18
```


Next, see if the hypothesis of proportional hazards is not violated. By the results obtained, the proportional hazards hypothesis is not violated.
```{r}
pph18<-cox.zph(fit.cox.all18, transform="km", global=TRUE)
pph18
```



### Cox Robust

In order to see if the results obtained before are consistent, a robust version of the Cox regression model is going to be performed.
First call the following package.

```{r}
library(coxrobust)
```

To see the results obtained by fitting the robust Cox using exponential weights with a 5% level of trimming, do the following:
```{r}
fit.coxrobust18<-coxr(Surv(data.time.status$time,data.time.status$status) ~ .,data=data.survival18,trunc=0.95, f.weight="exponential",singular.ok=T,model=F)
fit.coxrobust18

```


## Cox Robust based on Heritier

Another proposal of the robust Cox is given by Heritier (2009), and is going to be performed next.

```{r}

filepath="../functions/"
source(paste(filepath,"Chapter7_functions.r",sep=""))

```

The fitted model using exponential weights with 5% level of trimming is given bellow. 
```{r}
fit.coxrobust.H18=rcoxph(data.time.status$time,data.time.status$status,data.survival18,wt.type="exponential",quant=.95)
fit.coxrobust.H18$coefficients

```

The only genes statistically significant were: CRYAB and SPARC, for Cox's and Cox's robust. 

The next figure shows that observations 113 and 219 are identified as influential observations in the sense that the weight given for each one is the lowest among the others.
```{r, echo=FALSE}
n<-dim(data.survival18)[1]## number of observations of the data
#n
par(mfrow=c(1,1))
K=quantile(fit.coxrobust.H18$Lambda.rob*exp(fit.coxrobust.H18$xbeta.rob),probs=c(.95))
nb=seq(1,n,1)
logwt0=-log(fit.coxrobust.H18$wt[data.time.status$status==0])*K
logwt1=-log(fit.coxrobust.H18$wt[data.time.status$status==1])*K
n0=nb[data.time.status$status==0]
n1=nb[data.time.status$status==1]
plot(-log(fit.coxrobust.H18$wt)*K,ylab="-K log(weight)",xlab="Case number",xlim=c(0,n),ylim=c(0,8),type="n")
points(n0,logwt0,pch=21)
points(n1,logwt1,pch=19)
abline(h=K)


legend(400,7,c("Censored","Dead"),pch=c(21,19))
#K # value of K to put in the next text
text(-3.75,2.5,"K=1.81",pos=4)

## In order to identify the observation that have the lowest weight (in this case the highest because the scale was changed)
#sort(-log(fit.coxrobust.H$wt)*K,decreasing = TRUE) ## para saber quem são as aobservações com valor mais elevado
text(114,4,"case 113",pos=4)
text(220,4,"case 219",pos=4)

```

In the next section it is going to be evaluated if the observations identified before are or not influential.

## Residual for survival analisys
In order to identify which of the observations might be outliers the martingale residual is going to be performed.


### Martingale
The results regarding the martingale residuals are shown bellow. Observation 219 presented the lowest values for the martingale residual.

```{r}
res.mart.18<- resid(fit.cox.all18,type="martingale")
```


```{r, echo=FALSE}
par(mfrow=c(1,1))
plot(res.mart.18,main="Martingale Residuals",xlab="Case number")
abline(h=0)
#identify(res.mart.18)
```
```{r}
head(sort(res.mart.18))
```

NOTE: I can not find a way to identify in R markdown a way of identify the observations.





# TCGA ovarian cancer 22 genes

```{r}
data.survival22<- read.table("./tcga22.txt",dec = ",", header = TRUE)

head(data.survival22)
```


And the matrix for the time and status is given by:

```{r}
data.time.status<-read.table("./timestatus.txt",dec = ",", header = TRUE)
head(data.time.status)
```



## Survival analysis
The Cox's proportional hazard was performed in all variables. 

### Cox Regression model

Let's start the analysis for the Cox regression model.
```{r}
library(survival)
```

Fit the Cox proportional hazards to the dataset. 

```{r}
fit.cox.all22<-coxph(Surv(data.time.status$time,data.time.status$status) ~ .,data=data.survival22)
fit.cox.all22
```


Next, see if the hypothesis of proportional hazards is not violated. By the results obtained, the proportional hazards hypothesis is not violated.
```{r}
pph22<-cox.zph(fit.cox.all22, transform="km", global=TRUE)
pph22
```



### Cox Robust

In order to see if the results obtained before are consistent, a robust version of the Cox regression model is going to be performed.
First call the following package.

```{r}
library(coxrobust)
```

To see the results obtained by fitting the robust Cox using exponential weights with a 5% level of trimming, do the following:
```{r}
fit.coxrobust22<-coxr(Surv(data.time.status$time,data.time.status$status) ~ .,data=data.survival22,trunc=0.95, f.weight="exponential",singular.ok=T,model=F)
fit.coxrobust22

```


### Cox Robust based on Heritier

Another proposal of the robust Cox is given by Heritier (2009), and is going to be performed next.

```{r}

filepath="../functions/"
source(paste(filepath,"Chapter7_functions.r",sep=""))

```

The fitted model using exponential weights with 5% level of trimming is given bellow. 
```{r}
fit.coxrobust.H22=rcoxph(data.time.status$time,data.time.status$status,data.survival22,wt.type="exponential",quant=.95)
fit.coxrobust.H22$coefficients

```

The only genes statistically significant were: BRCA2 and PALB2 for the Cox's and PALB2 for the robust version of Heritier.

The next figure shows that observations 114 and 211 are identified as influential observations in the sense that the weight given for each one is the lowest among the others.
```{r, echo=FALSE}
n<-dim(data.survival22)[1]## number of observations of the data
n
par(mfrow=c(1,1))
K=quantile(fit.coxrobust.H22$Lambda.rob*exp(fit.coxrobust.H22$xbeta.rob),probs=c(.95))
nb=seq(1,n,1)
logwt0=-log(fit.coxrobust.H22$wt[data.time.status$status==0])*K
logwt1=-log(fit.coxrobust.H22$wt[data.time.status$status==1])*K
n0=nb[data.time.status$status==0]
n1=nb[data.time.status$status==1]
plot(-log(fit.coxrobust.H22$wt)*K,ylab="-K log(weight)",xlab="Case number",xlim=c(0,n),ylim=c(0,6),type="n")
points(n0,logwt0,pch=21)
points(n1,logwt1,pch=19)
abline(h=K)


legend(430,5,c("Censored","Dead"),pch=c(21,19))
K # value of K to put in the next text
text(-3.75,2,"K=1.71",pos=4)

## In order to identify the observation that have the lowest weight (in this case the highest because the scale was changed)
sort(-log(fit.coxrobust.H22$wt)*K,decreasing = TRUE) ## para saber quem são as aobservações com valor mais elevado
text(115,5,"case 114",pos=4)
text(212,3.5,"case 211",pos=4)


```

In the next section it is going to be evaluated if the observations identified before are or not influential.

## Outlier detection methods for survival analisys
In order to identify which of the observations might be outliers the martingale residual is going to be performed.


#### Martingale
The results regarding the martingale residuals are shown bellow. Observations 114 and 211 presented the lowest values for the martingale residual.

```{r}
res.mart.22<- resid(fit.cox.all22,type="martingale")
```


```{r, echo=FALSE}
par(mfrow=c(1,1))
plot(res.mart.22,main="Martingale Residuals",xlab="Case number")
abline(h=0)
#identify(res.mart.22)
```
```{r}
head(sort(res.mart.22))
```



# TCGA ovarian cancer 63 genes

```{r}
data.survival63<- read.table("./tcga63.txt",dec = ",", header = TRUE)

head(data.survival63)
```


And the matrix for the time and status is given by:

```{r}
data.time.status<-read.table("./timestatus.txt",dec = ",", header = TRUE)
head(data.time.status)
```



## Survival analysis
The Cox's proportional hazard was performed in all variables. 

### Cox Regression model

Let's start the analysis for the Cox regression model.
```{r}
library(survival)
```

Fit the Cox proportional hazards to the dataset. 

```{r}
fit.cox.all63<-coxph(Surv(data.time.status$time,data.time.status$status) ~ .,data=data.survival63)
fit.cox.all63
```


Next, see if the hypothesis of proportional hazards is not violated. By the results obtained, the proportional hazards hypothesis is not violated.
```{r}
pph63<-cox.zph(fit.cox.all63, transform="km", global=TRUE)
pph63
```

Genes not significant: TGM2, MGAT4C, TREML2, FZD10, ELA3A, SLC37A4, HSPB7, FOXE3, CLTCL1, GAPDHS, UBE2J1.


### Cox Robust

In order to see if the results obtained before are consistent, a robust version of the Cox regression model is going to be performed.
First call the following package.

```{r}
library(coxrobust)
```

To see the results obtained by fitting the robust Cox using exponential weights with a 5% level of trimming, do the following:
```{r}
fit.coxrobust63<-coxr(Surv(data.time.status$time,data.time.status$status) ~ .,data=data.survival63,trunc=0.95, f.weight="exponential",singular.ok=T,model=F)
fit.coxrobust63

```


### Cox Robust based on Heritier

Another proposal of the robust Cox is given by Heritier (2009), and is going to be performed next.

```{r}

filepath="../functions/"
source(paste(filepath,"Chapter7_functions.r",sep=""))

```

The fitted model using exponential weights with 5% level of trimming is given bellow. 
```{r}
fit.coxrobust.H63=rcoxph(data.time.status$time,data.time.status$status,data.survival63,wt.type="exponential",quant=.95)
fit.coxrobust.H63$coefficients

```

Genes not significant: TGM2, ALG8, MGAT4, TREML2, FZD10, SLC37A4, HSPB7, FOXE3, CLTCL1, PRR16, GAPDHS, SDF2L1, and UBE2J1.


The next figure shows that observations 39 and 350 are identified as influential observations in the sense that the weight given for each one is the lowest among the others.
```{r, echo=FALSE}
n<-dim(data.survival63)[1]## number of observations of the data
n

par(mfrow=c(1,1))
K=quantile(fit.coxrobust.H63$Lambda.rob*exp(fit.coxrobust.H63$xbeta.rob),probs=c(.95))
nb=seq(1,n,1)
logwt0=-log(fit.coxrobust.H63$wt[data.time.status$status==0])*K
logwt1=-log(fit.coxrobust.H63$wt[data.time.status$status==1])*K
n0=nb[data.time.status$status==0]
n1=nb[data.time.status$status==1]
plot(-log(fit.coxrobust.H63$wt)*K,ylab="-K log(weight)",xlab="Case number",xlim=c(0,n),ylim=c(0,15),type="n")
points(n0,logwt0,pch=21)
points(n1,logwt1,pch=19)
abline(h=K)


legend(450,12,c("Censored","Dead"),pch=c(21,19))
K # value of K to put in the next text
text(-3.75,2.5,"K=2.03",pos=4)

## In order to identify the observation that have the lowest weight (in this case the highest because the scale was changed)
sort(-log(fit.coxrobust.H63$wt)*K,decreasing = TRUE) ## para saber quem são as aobservações com valor mais elevado
text(40,13,"case 39",pos=4)
text(351,5,"case 350",pos=4)


```

In the next section it is going to be evaluated if the observations identified before are or not influential.

## Outlier detection methods for survival analisys
In order to identify which of the observations might be outliers the martingale residual is going to be performed.


#### Martingale
The results regarding the martingale residuals are shown bellow. Observations 39 and 350 presented the lowest values for the martingale residual.

```{r}
res.mart.63<- resid(fit.cox.all63,type="martingale")
```


```{r, echo=FALSE}
par(mfrow=c(1,1))
plot(res.mart.63,main="Martingale Residuals",xlab="Case number")
abline(h=0)
#identify(res.mart.63)
```
```{r}
head(sort(res.mart.63))
```






# Rank Product test

In order to obtain a global result concerning each one of the techniques performed in the previous section, a rank product test is going to be used.
First a rank matrix is going to be obtained based on the ranks of each outliers detection method.

```{r}
res.mart.vec18<-as.vector(res.mart.18)
rank_martingale18<-rank(res.mart.vec18,ties.method = "first")
```


```{r}
res.mart.vec22<-as.vector(res.mart.22)
rank_martingale22<-rank(res.mart.vec22,ties.method = "first")
```


```{r}
res.mart.vec63<-as.vector(res.mart.63)
rank_martingale63<-rank(res.mart.vec63,ties.method = "first")
```


Now let us obtain based on the ranks above the rank product matrix.
First let's consider a rank matrix, with all the ranks obtained by each method.

```{r}
rankMat<-cbind(rank_martingale18,rank_martingale22,rank_martingale63)
```

In second obtain the rank product, which is obtained based on the input of the rank matrix.
```{r}
RankProduct = function(rankMat)
{
  return(apply(rankMat, 1, prod))
}

rankproduct<-RankProduct(rankMat)
```

## p-values

In order to obtain the p-values the algorithm proposed by Heskes (2014) was used. Notice that the p-values are based on the geometric mean of upper and lower bounds, defined recursively.

```{r}
filepath="../functions/"
source(paste(filepath,"Heskes_pvalues.R",sep=""))

```

The input to obtain the p-values is the following:

* rho: the rank product matrix;
* the number of observations, n;
* the number of methods used, k.

```{r}
rho=rankproduct
n<-dim(data.survival18)[1]
k<-dim(rankMat)[2]
```

The p-values are obtained by the following, where Delta option is the geometric mean.

```{r}
pvalues<-as.vector(rankprodbounds(rho,n,k,Delta ='geometric'))
```

In order to obtain the observations with the lowest p-values we can do the lines:

```{r}
pvalues.matrix<-cbind(pvalues)
which(pvalues.matrix<0.01)## pvalues < to 1%
which(pvalues.matrix<0.05) ## pvalues < to 5%
```



## q-values
Another key issue when performing Rank product tests is related with the  multiple testing problem. In fact, since many observations are tested, type-I errors (false positives) will increase.

The FDR, which is the expected proportion of false positives among all tests that are significant, sorts in an ascendant order the p-values and divides them by their percentile rank.  The measure used to determine the FDR is the q-value.  For the p-value: 0.05 implies that 5% of all tests will result in false positives, instead, for the q-value: 0.05 implies that 5% of significant tests will result in false positives. In this sense the q-value is able to control the number of false discoveries in those tests. For this reason it has the ability of finding truly significant results. 

The q-value package must be installed.

```{r,include=FALSE}
source("https://bioconductor.org/biocLite.R")
biocLite("qvalue")

```

```{r}
library(qvalue)
```


Here the only input that is need is the p-values obtained before for each observation.

```{r}
qobj<-qvalue(pvalues)
```

```{r}
qvalues<-qobj$qvalues
which(qvalues<0.01)
which(qvalues<0.05)
```



In order to resume all the information, let us get an id vector and a data frame that combines the ranks the pvalues and the qvalues.

```{r}
id<-as.vector(seq(1:n))# vector with the id
tcga.outliers<-as.data.frame(cbind(id,rank_martingale18,rank_martingale22,rank_martingale63,pvalues,qvalues))
```

Next we can obtain the top 15 outlier observations, with the qvalues sorted in increasing order. 

```{r}
sort.tcga.outliers <- tcga.outliers[order(qvalues) , ]
sort.tcga.outliers[1:15,]
```

Notice the observations the with the lowest qvalues lead us to believe that they are influential observations and should taken for further study.


```{r}
tcga_outliers_time<-cbind(id,data.time.status$time,data.time.status$status,rank_martingale18,rank_martingale22,rank_martingale63,pvalues,qvalues)
sort.tcga.outliers_time <- tcga_outliers_time[order(qvalues) , ]
sort.tcga.outliers_time[1:15,]
```










