---
title: "Simulation data - Survival n=100; p=20"
author: "Eunice"
output: html_document
---




# Covariates and time and censoring simulation

First, a matrix with three covariates is simulated following a normal distribution, with dimension n=100.



```{r}
set.seed(1980)
x <- mvtnorm::rmvnorm(n = 100 ,mean = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) , sigma = matrix(diag(20),nrow=20, ncol =20))
data.survival<-as.data.frame(x)
```

In second a simulation of the times and censoring is done using the following function, with the arguments:

- data        : data frame with random X
- n_outliers  : number of outlier observations
- beta_pure   : beta of pure population
- beta_noise  : beta of noisy population
- cens_amount : amount of censoring , typical = 0.2

```{r}
get_simdata_2coxes_bybook <- function(data,n_outliers,beta_pure,beta_noise,cens_amount,shape,scale)
  {
  times <- rep(x=0, nrow(data) )
  
  ## Generate the regular observations
  #
  for(i in 1:(nrow(data)-n_outliers)){
    times[i] <- ( - log( runif(n = 1 , min=0 , max= 1) )/( scale * exp( sum( beta_pure*data[i,] ) )  )) ^(1/shape) 
                     
  }
  ## Generate the outlier observations
  #
  for(i in (nrow(data)-n_outliers+1):(nrow(data)) ) {
    
    times[i] <- (  - log( runif(n = 1 , min=0 , max= 1) )/( scale * exp( sum( beta_noise*data[i,] ) )  ) )^(1/shape) 
    
  }
  ##### censoring
  #
  cens        <- rep(x=1, nrow(data))  # censoring vector , 0-alive 1-dead
  CENS_AMOUNT <- cens_amount
  decider_censoring <- runif(n=nrow(data) , min= 0, max= 1 )
  
  ## Random censoring
  #
  cens[decider_censoring<CENS_AMOUNT]  <-   0 
  
  
  return( cbind(times,cens) )
}
```


Bsed on this function, several scenarios for the parameters are presented.
Notice that in all the scenarios the amount of outliers is 10% of the total of observations. 

## 1st scenario 

Notice that in all the scenarios the amount of outliers is 10% of the total of observations

- beta_pure   : 1
- beta_noise  : -1
- cens_amount :0.2
- scale       : 1
- shape       : 1

```{r}
simulation_1<-get_simdata_2coxes_bybook(data.survival,10,1,-1,0.2,1,1)
data.time.status_1<-as.data.frame(simulation_1)
```


Libraries needed:
```{r}
library(survival)
library(coxrobust)
library(OutlierDC)
library(qvalue)
```


Package for outlier detection based on the concordance index.
```{r,include=FALSE}
# Package is retrieved from our group server
devtools::install_url('http://sels.tecnico.ulisboa.pt/gitlab/SOUND/bcsod/repository/archive.zip')
```
```{r}
library("BCSOD")
```



Functions used:

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


### Survival Analysis

```{r}
fit.cox_1<-coxph(Surv(data.time.status_1$times, data.time.status_1$cens)~.,data = data.survival)
fit.cox_1
```

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


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


The fitted model using exponential weights with 5% level of trimming is given bellow. Notice that the results show that all covariates are significant for a 5% level of significant.

```{r}
fit.coxrobust.H_1=rcoxph(data.time.status_1$times,data.time.status_1$cens,data.survival,wt.type="exponential",quant=.95)
fit.coxrobust.H_1$coefficients

```


```{r, echo=FALSE}
n<-dim(data.survival)[1]## number of observations of the data

par(mfrow=c(1,1))
K_1=quantile(fit.coxrobust.H_1$Lambda.rob*exp(fit.coxrobust.H_1$xbeta.rob),probs=c(.95))
nb=seq(1,n,1)
logwt0=-log(fit.coxrobust.H_1$wt[data.time.status_1$cens==0])*K_1
logwt1=-log(fit.coxrobust.H_1$wt[data.time.status_1$cens==1])*K_1
n0=nb[data.time.status_1$cens==0]
n1=nb[data.time.status_1$cens==1]
plot(-log(fit.coxrobust.H_1$wt)*K_1,ylab="-K log(weight)",xlab="Case number",xlim=c(0,n),ylim=c(0,20),type="n")
points(n0,logwt0,pch=21)
points(n1,logwt1,pch=19)
abline(h=K_1)
legend(10,20,c("Censored","Dead"),pch=c(21,19))
#K # value of K to put in the next text
text(0,2.3,"K=2.1",pos=4)
```


```{r}
sort(-log(fit.coxrobust.H_1$wt)*K_1,decreasing = TRUE)
```



### Outlier detection methods for survival analisys
In order to identify which of the observations might be outliers some methods for outlier detection for survival analysis are going to be performed.

#### Residuals
The first method used for outlier detection is the deviance residuals.


```{r}
res.dev_1<- resid(fit.cox_1,type="deviance")
```


```{r, echo=FALSE}
par(mfrow=c(1,1))
plot(res.dev_1,main="Deviance Residuals",xlab="Case number")
abline(h=0)
```


```{r}
head(sort(abs(res.dev_1),decreasing = TRUE))
```


#### Concordance c-index

```{r, include=FALSE}
set.seed(1980)
n<-dim(data.survival)[1]
outliers_dbht_1<- bcsod(surv.object = Surv(data.time.status_1$times,data.time.status_1$cens),covariate.data =data.survival ,sod.method = "dbht", B=1000 , B.N = n)
```



Here are the top observations with lowest p-values.
```{r}
outliers_dbht_1[1:10,]
```


```{r}
dbht_1<-as.data.frame(outliers_dbht_1)# convert to a dataframe (more easy to handle)
names(dbht_1)[1]="id"
dbht.order_1<-dbht_1[order(dbht_1$id),]
```

#### Censored quantile regression

The score algorithm arises from the outlier detection algorithm for the Cox's regression model for censored data proposed by Nardi (1999) but now based on quantile regression.

```{r}
odcS_1<-odc(formula=Surv((data.time.status_1$times),data.time.status_1$cens)~ data.survival$V1+data.survival$V2+data.survival$V3, data = data.survival,method="score",rq.model="Wang",h=0.05)
odcS_1
```




The top outliers are:
```{r}
head(sort(odcS_1@score,decreasing = TRUE))
```



Put in form of a vector the results obtained before
```{r}
odS.vec_1<-as.vector(odcS_1@score)
```

Define a matrix with the id and the values for the residual algorithm.

```{r}
n<-dim(data.survival)[1]#dimension of the data
id<-as.vector(seq(1:n))# vector with the id
odS.matrix_1<-as.data.frame(cbind(id,odS.vec_1))
```



### Rank Product test

First a rank matrix is going to be obtained based on the ranks of each outliers detection method.

```{r}
res.dev.vec_1<-as.vector(abs(res.dev_1))
rank_deviance_1<-rank(-res.dev.vec_1,ties.method = "first") 
rank_dbht_1<-rank(dbht.order_1$pvalue,ties.method = "first")
rank_score_1<-rank(-(odS.matrix_1$odS.vec_1),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_1<-cbind(rank_deviance_1,rank_dbht_1,rank_score_1)
```

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

rankproduct<-RankProduct(rankMat_1)
```


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.
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.survival)[1]
k_1<-dim(rankMat_1)[2]
```

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

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

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

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


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. 

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

```{r}
qobj_1<-qvalue(pvalues_1)
```

```{r}
qvalues_1<-qobj_1$qvalues
which(qvalues_1<0.01)
which(qvalues_1<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
outliers_1<-as.data.frame(cbind(id,rank_deviance_1,rank_dbht_1,rank_score_1,pvalues_1,qvalues_1))
sort.outliers_1 <- outliers_1[order(qvalues_1) , ]
sort.outliers_1[1:10,]
```

## 2nd scenario 

Notice that in all the scenarios the amount of outliers is 10% of the total of observations

- beta_pure   : 1
- beta_noise  : -1
- cens_amount :0.2
- scale       : 1.5
- shape       : 0.5

```{r}
simulation_2<-get_simdata_2coxes_bybook(data.survival,10,1,-1,0.2,1.5,0.5)
data.time.status_2<-as.data.frame(simulation_2)
```


### Survival Analysis

```{r}
fit.cox_2<-coxph(Surv(data.time.status_2$times, data.time.status_2$cens)~.,data = data.survival)
fit.cox_2
```

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


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


The fitted model using exponential weights with 5% level of trimming is given bellow. Notice that the results show that all covariates are significant for a 5% level of significant.

```{r}
fit.coxrobust.H_2=rcoxph(data.time.status_2$times,data.time.status_2$cens,data.survival,wt.type="exponential",quant=.95)
fit.coxrobust.H_2$coefficients

```




```{r}
K_2=quantile(fit.coxrobust.H_2$Lambda.rob*exp(fit.coxrobust.H_2$xbeta.rob),probs=c(.95))
sort(-log(fit.coxrobust.H_2$wt)*K_2,decreasing = TRUE)
```



### Outlier detection methods for survival analisys
In order to identify which of the observations might be outliers some methods for outlier detection for survival analysis are going to be performed.

#### Residuals
The first method used for outlier detection is the deviance residuals.


```{r}
res.dev_2<- resid(fit.cox_2,type="deviance")
```


```{r, echo=FALSE}
par(mfrow=c(1,1))
plot(res.dev_2,main="Deviance Residuals",xlab="Case number")
abline(h=0)
```


```{r}
head(sort(abs(res.dev_2),decreasing = TRUE))
```


#### Concordance c-index

```{r, include=FALSE}
set.seed(1980)
n<-dim(data.survival)[1]
outliers_dbht_2<- bcsod(surv.object = Surv(data.time.status_2$times,data.time.status_2$cens),covariate.data =data.survival ,sod.method = "dbht", B=1000 , B.N = n)
```



Here are the top observations with lowest p-values.
```{r}
outliers_dbht_2[1:10,]
```


```{r}
dbht_2<-as.data.frame(outliers_dbht_2)# convert to a dataframe (more easy to handle)
names(dbht_2)[1]="id"
dbht.order_2<-dbht_2[order(dbht_2$id),]
```

#### Censored quantile regression

The score algorithm arises from the outlier detection algorithm for the Cox's regression model for censored data proposed by Nardi (1999) but now based on quantile regression.

```{r}
odcS_2<-odc(formula=Surv((data.time.status_2$times),data.time.status_2$cens)~ data.survival$V1+data.survival$V2+data.survival$V3, data = data.survival,method="score",rq.model="Wang",h=0.05)
odcS_2
```




The top outliers are:
```{r}
head(sort(odcS_2@score,decreasing = TRUE))
```



Put in form of a vector the results obtained before
```{r}
odS.vec_2<-as.vector(odcS_2@score)
```

Define a matrix with the id and the values for the residual algorithm.

```{r}
n<-dim(data.survival)[1]#dimension of the data
id<-as.vector(seq(1:n))# vector with the id
odS.matrix_2<-as.data.frame(cbind(id,odS.vec_2))
```



### Rank Product test

First a rank matrix is going to be obtained based on the ranks of each outliers detection method.

```{r}
res.dev.vec_2<-as.vector(abs(res.dev_2))
rank_deviance_2<-rank(-res.dev.vec_2,ties.method = "first") 
rank_dbht_2<-rank(dbht.order_2$pvalue,ties.method = "first")
rank_score_2<-rank(-(odS.matrix_2$odS.vec_2),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_2<-cbind(rank_deviance_2,rank_dbht_2,rank_score_2)
```

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

rankproduct<-RankProduct(rankMat_2)
```


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.
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.survival)[1]
k_2<-dim(rankMat_2)[2]
```

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

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

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

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


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. 

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

```{r}
qobj_2<-qvalue(pvalues_2)
```

```{r}
qvalues_2<-qobj_2$qvalues
which(qvalues_2<0.01)
which(qvalues_2<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
outliers_2<-as.data.frame(cbind(id,rank_deviance_2,rank_dbht_2,rank_score_2,pvalues_2,qvalues_2))
sort.outliers_2 <- outliers_2[order(qvalues_2) , ]
sort.outliers_2[1:10,]
```


## 3rd scenario 

Notice that in all the scenarios the amount of outliers is 10% of the total of observations

- beta_pure   : 1
- beta_noise  : -1
- cens_amount : 0.2
- scale       : 0.5
- shape       : 1.5

```{r}
simulation_3<-get_simdata_2coxes_bybook(data.survival,10,1,-1,0.2,0.5,1.5)
data.time.status_3<-as.data.frame(simulation_3)
```


### Survival Analysis

```{r}
fit.cox_3<-coxph(Surv(data.time.status_3$times, data.time.status_3$cens)~.,data = data.survival)
fit.cox_3
```

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


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


The fitted model using exponential weights with 5% level of trimming is given bellow. Notice that the results show that all covariates are significant for a 5% level of significant.

```{r}
fit.coxrobust.H_3=rcoxph(data.time.status_3$times,data.time.status_3$cens,data.survival,wt.type="exponential",quant=.95)
fit.coxrobust.H_3$coefficients

```




```{r}
K_3=quantile(fit.coxrobust.H_3$Lambda.rob*exp(fit.coxrobust.H_3$xbeta.rob),probs=c(.95))
sort(-log(fit.coxrobust.H_3$wt)*K_3,decreasing = TRUE)
```



### Outlier detection methods for survival analisys
In order to identify which of the observations might be outliers some methods for outlier detection for survival analysis are going to be performed.

#### Residuals
The first method used for outlier detection is the deviance residuals.


```{r}
res.dev_3<- resid(fit.cox_3,type="deviance")
```


```{r, echo=FALSE}
par(mfrow=c(1,1))
plot(res.dev_3,main="Deviance Residuals",xlab="Case number")
abline(h=0)
```


```{r}
head(sort(abs(res.dev_3),decreasing = TRUE))
```


#### Concordance c-index

```{r, include=FALSE}
set.seed(1980)
n<-dim(data.survival)[1]
outliers_dbht_3<- bcsod(surv.object = Surv(data.time.status_3$times,data.time.status_3$cens),covariate.data =data.survival ,sod.method = "dbht", B=1000 , B.N = n)
```



Here are the top observations with lowest p-values.
```{r}
outliers_dbht_3[1:10,]
```


```{r}
dbht_3<-as.data.frame(outliers_dbht_3)# convert to a dataframe (more easy to handle)
names(dbht_3)[1]="id"
dbht.order_3<-dbht_3[order(dbht_3$id),]
```

#### Censored quantile regression

The score algorithm arises from the outlier detection algorithm for the Cox's regression model for censored data proposed by Nardi (1999) but now based on quantile regression.

```{r}
odcS_3<-odc(formula=Surv((data.time.status_3$times),data.time.status_3$cens)~ data.survival$V1+data.survival$V2+data.survival$V3, data = data.survival,method="score",rq.model="Wang",h=0.05)
odcS_3
```




The top outliers are:
```{r}
head(sort(odcS_3@score,decreasing = TRUE))
```



Put in form of a vector the results obtained before
```{r}
odS.vec_3<-as.vector(odcS_3@score)
```

Define a matrix with the id and the values for the residual algorithm.

```{r}
n<-dim(data.survival)[1]#dimension of the data
id<-as.vector(seq(1:n))# vector with the id
odS.matrix_3<-as.data.frame(cbind(id,odS.vec_3))
```



### Rank Product test

First a rank matrix is going to be obtained based on the ranks of each outliers detection method.

```{r}
res.dev.vec_3<-as.vector(abs(res.dev_3))
rank_deviance_3<-rank(-res.dev.vec_3,ties.method = "first") 
rank_dbht_3<-rank(dbht.order_3$pvalue,ties.method = "first")
rank_score_3<-rank(-(odS.matrix_3$odS.vec_3),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_3<-cbind(rank_deviance_3,rank_dbht_3,rank_score_3)
```

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

rankproduct<-RankProduct(rankMat_3)
```


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.
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.survival)[1]
k_3<-dim(rankMat_3)[2]
```

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

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

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

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


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. 

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

```{r}
qobj_3<-qvalue(pvalues_3)
```

```{r}
qvalues_3<-qobj_3$qvalues
which(qvalues_3<0.01)
which(qvalues_3<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
outliers_3<-as.data.frame(cbind(id,rank_deviance_3,rank_dbht_3,rank_score_3,pvalues_3,qvalues_3))
sort.outliers_3 <- outliers_3[order(qvalues_3) , ]
sort.outliers_3[1:10,]
```





## 4th scenario 

Notice that in all the scenarios the amount of outliers is 10% of the total of observations

- beta_pure   : 1
- beta_noise  : -0.2
- cens_amount : 0.2
- scale       : 1
- shape       : 1

```{r}
simulation_4<-get_simdata_2coxes_bybook(data.survival,10,1,-0.2,0.2,1,1)
data.time.status_4<-as.data.frame(simulation_4)
```


### Survival Analysis

```{r}
fit.cox_4<-coxph(Surv(data.time.status_4$times, data.time.status_4$cens)~.,data = data.survival)
fit.cox_4
```

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


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


The fitted model using exponential weights with 5% level of trimming is given bellow. Notice that the results show that all covariates are significant for a 5% level of significant.

```{r}
fit.coxrobust.H_4=rcoxph(data.time.status_4$times,data.time.status_4$cens,data.survival,wt.type="exponential",quant=.95)
fit.coxrobust.H_4$coefficients

```




```{r}
K_4=quantile(fit.coxrobust.H_4$Lambda.rob*exp(fit.coxrobust.H_4$xbeta.rob),probs=c(.95))
sort(-log(fit.coxrobust.H_4$wt)*K_4,decreasing = TRUE)
```



### Outlier detection methods for survival analisys
In order to identify which of the observations might be outliers some methods for outlier detection for survival analysis are going to be performed.

#### Residuals
The first method used for outlier detection is the deviance residuals.


```{r}
res.dev_4<- resid(fit.cox_4,type="deviance")
```


```{r, echo=FALSE}
par(mfrow=c(1,1))
plot(res.dev_4,main="Deviance Residuals",xlab="Case number")
abline(h=0)
```


```{r}
head(sort(abs(res.dev_4),decreasing = TRUE))
```


#### Concordance c-index

```{r, include=FALSE}
set.seed(1980)
n<-dim(data.survival)[1]
outliers_dbht_4<- bcsod(surv.object = Surv(data.time.status_4$times,data.time.status_4$cens),covariate.data =data.survival ,sod.method = "dbht", B=1000 , B.N = n)
```



Here are the top observations with lowest p-values.
```{r}
outliers_dbht_4[1:10,]
```


```{r}
dbht_4<-as.data.frame(outliers_dbht_4)# convert to a dataframe (more easy to handle)
names(dbht_4)[1]="id"
dbht.order_4<-dbht_4[order(dbht_4$id),]
```

#### Censored quantile regression

The score algorithm arises from the outlier detection algorithm for the Cox's regression model for censored data proposed by Nardi (1999) but now based on quantile regression.

```{r}
odcS_4<-odc(formula=Surv((data.time.status_4$times),data.time.status_4$cens)~ data.survival$V1+data.survival$V2+data.survival$V3, data = data.survival,method="score",rq.model="Wang",h=0.05)
odcS_4
```




The top outliers are:
```{r}
head(sort(odcS_4@score,decreasing = TRUE))
```



Put in form of a vector the results obtained before
```{r}
odS.vec_4<-as.vector(odcS_4@score)
```

Define a matrix with the id and the values for the residual algorithm.

```{r}
n<-dim(data.survival)[1]#dimension of the data
id<-as.vector(seq(1:n))# vector with the id
odS.matrix_4<-as.data.frame(cbind(id,odS.vec_4))
```



### Rank Product test

First a rank matrix is going to be obtained based on the ranks of each outliers detection method.

```{r}
res.dev.vec_4<-as.vector(abs(res.dev_4))
rank_deviance_4<-rank(-res.dev.vec_4,ties.method = "first") 
rank_dbht_4<-rank(dbht.order_4$pvalue,ties.method = "first")
rank_score_4<-rank(-(odS.matrix_4$odS.vec_4),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_4<-cbind(rank_deviance_4,rank_dbht_4,rank_score_4)
```

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

rankproduct<-RankProduct(rankMat_4)
```


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.
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.survival)[1]
k_4<-dim(rankMat_4)[2]
```

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

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

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

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


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. 

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

```{r}
qobj_4<-qvalue(pvalues_4)
```

```{r}
qvalues_4<-qobj_4$qvalues
which(qvalues_4<0.01)
which(qvalues_4<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
outliers_4<-as.data.frame(cbind(id,rank_deviance_4,rank_dbht_4,rank_score_4,pvalues_4,qvalues_4))
sort.outliers_4 <- outliers_4[order(qvalues_4) , ]
sort.outliers_4[1:10,]
```



## 5th scenario 

Notice that in all the scenarios the amount of outliers is 10% of the total of observations

- beta_pure   : 1
- beta_noise  : -0.2
- cens_amount : 0.2
- scale       : 1.5
- shape       : 0.5

```{r}
simulation_5<-get_simdata_2coxes_bybook(data.survival,10,1,-0.2,0.2,1.5,0.5)
data.time.status_5<-as.data.frame(simulation_5)
```


### Survival Analysis

```{r}
fit.cox_5<-coxph(Surv(data.time.status_5$times, data.time.status_5$cens)~.,data = data.survival)
fit.cox_5
```

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


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


The fitted model using exponential weights with 5% level of trimming is given bellow. Notice that the results show that all covariates are significant for a 5% level of significant.

```{r}
fit.coxrobust.H_5=rcoxph(data.time.status_5$times,data.time.status_5$cens,data.survival,wt.type="exponential",quant=.95)
fit.coxrobust.H_5$coefficients

```




```{r}
K_5=quantile(fit.coxrobust.H_5$Lambda.rob*exp(fit.coxrobust.H_5$xbeta.rob),probs=c(.95))
sort(-log(fit.coxrobust.H_5$wt)*K_5,decreasing = TRUE)
```



### Outlier detection methods for survival analisys
In order to identify which of the observations might be outliers some methods for outlier detection for survival analysis are going to be performed.

#### Residuals
The first method used for outlier detection is the deviance residuals.


```{r}
res.dev_5<- resid(fit.cox_5,type="deviance")
```


```{r, echo=FALSE}
par(mfrow=c(1,1))
plot(res.dev_5,main="Deviance Residuals",xlab="Case number")
abline(h=0)
```


```{r}
head(sort(abs(res.dev_5),decreasing = TRUE))
```


#### Concordance c-index

```{r, include=FALSE}
set.seed(1980)
n<-dim(data.survival)[1]
outliers_dbht_5<- bcsod(surv.object = Surv(data.time.status_5$times,data.time.status_5$cens),covariate.data =data.survival ,sod.method = "dbht", B=1000 , B.N = n)
```



Here are the top observations with lowest p-values.
```{r}
outliers_dbht_5[1:10,]
```


```{r}
dbht_5<-as.data.frame(outliers_dbht_5)# convert to a dataframe (more easy to handle)
names(dbht_5)[1]="id"
dbht.order_5<-dbht_5[order(dbht_5$id),]
```

#### Censored quantile regression

The score algorithm arises from the outlier detection algorithm for the Cox's regression model for censored data proposed by Nardi (1999) but now based on quantile regression.

```{r}
odcS_5<-odc(formula=Surv((data.time.status_5$times),data.time.status_5$cens)~ data.survival$V1+data.survival$V2+data.survival$V3, data = data.survival,method="score",rq.model="Wang",h=0.05)
odcS_5
```




The top outliers are:
```{r}
head(sort(odcS_5@score,decreasing = TRUE))
```



Put in form of a vector the results obtained before
```{r}
odS.vec_5<-as.vector(odcS_5@score)
```

Define a matrix with the id and the values for the residual algorithm.

```{r}
n<-dim(data.survival)[1]#dimension of the data
id<-as.vector(seq(1:n))# vector with the id
odS.matrix_5<-as.data.frame(cbind(id,odS.vec_5))
```



### Rank Product test

First a rank matrix is going to be obtained based on the ranks of each outliers detection method.

```{r}
res.dev.vec_5<-as.vector(abs(res.dev_5))
rank_deviance_5<-rank(-res.dev.vec_5,ties.method = "first") 
rank_dbht_5<-rank(dbht.order_5$pvalue,ties.method = "first")
rank_score_5<-rank(-(odS.matrix_5$odS.vec_5),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_5<-cbind(rank_deviance_5,rank_dbht_5,rank_score_5)
```

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

rankproduct<-RankProduct(rankMat_5)
```


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.
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.survival)[1]
k_5<-dim(rankMat_5)[2]
```

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

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

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

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


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. 

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

```{r}
qobj_5<-qvalue(pvalues_5)
```

```{r}
qvalues_5<-qobj_5$qvalues
which(qvalues_5<0.01)
which(qvalues_5<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
outliers_5<-as.data.frame(cbind(id,rank_deviance_5,rank_dbht_5,rank_score_5,pvalues_5,qvalues_5))
sort.outliers_5 <- outliers_5[order(qvalues_5) , ]
sort.outliers_5[1:10,]
```


## 6th scenario 

Notice that in all the scenarios the amount of outliers is 10% of the total of observations

- beta_pure   : 1
- beta_noise  : -0.2
- cens_amount : 0.2
- scale       : 0.5
- shape       : 1.5

```{r}
simulation_6<-get_simdata_2coxes_bybook(data.survival,10,1,-0.2,0.2,0.5,1.5)
data.time.status_6<-as.data.frame(simulation_6)
```


### Survival Analysis

```{r}
fit.cox_6<-coxph(Surv(data.time.status_6$times, data.time.status_6$cens)~.,data = data.survival)
fit.cox_6
```

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


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


The fitted model using exponential weights with 5% level of trimming is given bellow. Notice that the results show that all covariates are significant for a 5% level of significant.

```{r}
fit.coxrobust.H_6=rcoxph(data.time.status_6$times,data.time.status_6$cens,data.survival,wt.type="exponential",quant=.95)
fit.coxrobust.H_6$coefficients

```




```{r}
K_6=quantile(fit.coxrobust.H_6$Lambda.rob*exp(fit.coxrobust.H_6$xbeta.rob),probs=c(.95))
sort(-log(fit.coxrobust.H_6$wt)*K_6,decreasing = TRUE)
```



### Outlier detection methods for survival analisys
In order to identify which of the observations might be outliers some methods for outlier detection for survival analysis are going to be performed.

#### Residuals
The first method used for outlier detection is the deviance residuals.


```{r}
res.dev_6<- resid(fit.cox_6,type="deviance")
```


```{r, echo=FALSE}
par(mfrow=c(1,1))
plot(res.dev_6,main="Deviance Residuals",xlab="Case number")
abline(h=0)
```


```{r}
head(sort(abs(res.dev_6),decreasing = TRUE))
```


#### Concordance c-index

```{r, include=FALSE}
set.seed(1980)
n<-dim(data.survival)[1]
outliers_dbht_6<- bcsod(surv.object = Surv(data.time.status_6$times,data.time.status_6$cens),covariate.data =data.survival ,sod.method = "dbht", B=1000 , B.N = n)
```



Here are the top observations with lowest p-values.
```{r}
outliers_dbht_6[1:10,]
```


```{r}
dbht_6<-as.data.frame(outliers_dbht_6)# convert to a dataframe (more easy to handle)
names(dbht_6)[1]="id"
dbht.order_6<-dbht_6[order(dbht_6$id),]
```

#### Censored quantile regression

The score algorithm arises from the outlier detection algorithm for the Cox's regression model for censored data proposed by Nardi (1999) but now based on quantile regression.

```{r}
odcS_6<-odc(formula=Surv((data.time.status_6$times),data.time.status_6$cens)~ data.survival$V1+data.survival$V2+data.survival$V3, data = data.survival,method="score",rq.model="Wang",h=0.05)
odcS_6
```




The top outliers are:
```{r}
head(sort(odcS_6@score,decreasing = TRUE))
```



Put in form of a vector the results obtained before
```{r}
odS.vec_6<-as.vector(odcS_6@score)
```

Define a matrix with the id and the values for the residual algorithm.

```{r}
n<-dim(data.survival)[1]#dimension of the data
id<-as.vector(seq(1:n))# vector with the id
odS.matrix_6<-as.data.frame(cbind(id,odS.vec_6))
```



### Rank Product test

First a rank matrix is going to be obtained based on the ranks of each outliers detection method.

```{r}
res.dev.vec_6<-as.vector(abs(res.dev_6))
rank_deviance_6<-rank(-res.dev.vec_6,ties.method = "first") 
rank_dbht_6<-rank(dbht.order_6$pvalue,ties.method = "first")
rank_score_6<-rank(-(odS.matrix_6$odS.vec_6),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_6<-cbind(rank_deviance_6,rank_dbht_6,rank_score_6)
```

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

rankproduct<-RankProduct(rankMat_6)
```


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.
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.survival)[1]
k_6<-dim(rankMat_6)[2]
```

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

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

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

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


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. 

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

```{r}
qobj_6<-qvalue(pvalues_6)
```

```{r}
qvalues_6<-qobj_6$qvalues
which(qvalues_6<0.01)
which(qvalues_6<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
outliers_6<-as.data.frame(cbind(id,rank_deviance_6,rank_dbht_6,rank_score_6,pvalues_6,qvalues_6))
sort.outliers_6 <- outliers_6[order(qvalues_6) , ]
sort.outliers_6[1:10,]
```




## 7th scenario 

Notice that in all the scenarios the amount of outliers is 10% of the total of observations

- beta_pure   : 1
- beta_noise  : -1
- cens_amount : 0.3
- scale       : 1
- shape       : 1

```{r}
simulation_7<-get_simdata_2coxes_bybook(data.survival,10,1,-1,0.3,1,1)
data.time.status_7<-as.data.frame(simulation_7)
```


### Survival Analysis

```{r}
fit.cox_7<-coxph(Surv(data.time.status_7$times, data.time.status_7$cens)~.,data = data.survival)
fit.cox_7
```

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


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


The fitted model using exponential weights with 5% level of trimming is given bellow. Notice that the results show that all covariates are significant for a 5% level of significant.

```{r}
fit.coxrobust.H_7=rcoxph(data.time.status_7$times,data.time.status_7$cens,data.survival,wt.type="exponential",quant=.95)
fit.coxrobust.H_7$coefficients

```




```{r}
K_7=quantile(fit.coxrobust.H_7$Lambda.rob*exp(fit.coxrobust.H_7$xbeta.rob),probs=c(.95))
sort(-log(fit.coxrobust.H_7$wt)*K_7,decreasing = TRUE)
```



### Outlier detection methods for survival analisys
In order to identify which of the observations might be outliers some methods for outlier detection for survival analysis are going to be performed.

#### Residuals
The first method used for outlier detection is the deviance residuals.


```{r}
res.dev_7<- resid(fit.cox_7,type="deviance")
```


```{r, echo=FALSE}
par(mfrow=c(1,1))
plot(res.dev_7,main="Deviance Residuals",xlab="Case number")
abline(h=0)
```


```{r}
head(sort(abs(res.dev_7),decreasing = TRUE))
```


#### Concordance c-index

```{r, include=FALSE}
set.seed(1980)
n<-dim(data.survival)[1]
outliers_dbht_7<- bcsod(surv.object = Surv(data.time.status_7$times,data.time.status_7$cens),covariate.data =data.survival ,sod.method = "dbht", B=1000 , B.N = n)
```



Here are the top observations with lowest p-values.
```{r}
outliers_dbht_7[1:10,]
```


```{r}
dbht_7<-as.data.frame(outliers_dbht_7)# convert to a dataframe (more easy to handle)
names(dbht_7)[1]="id"
dbht.order_7<-dbht_7[order(dbht_7$id),]
```

#### Censored quantile regression

The score algorithm arises from the outlier detection algorithm for the Cox's regression model for censored data proposed by Nardi (1999) but now based on quantile regression.

```{r}
odcS_7<-odc(formula=Surv((data.time.status_7$times),data.time.status_7$cens)~ data.survival$V1+data.survival$V2+data.survival$V3, data = data.survival,method="score",rq.model="Wang",h=0.05)
odcS_7
```




The top outliers are:
```{r}
head(sort(odcS_7@score,decreasing = TRUE))
```



Put in form of a vector the results obtained before
```{r}
odS.vec_7<-as.vector(odcS_7@score)
```

Define a matrix with the id and the values for the residual algorithm.

```{r}
n<-dim(data.survival)[1]#dimension of the data
id<-as.vector(seq(1:n))# vector with the id
odS.matrix_7<-as.data.frame(cbind(id,odS.vec_7))
```



### Rank Product test

First a rank matrix is going to be obtained based on the ranks of each outliers detection method.

```{r}
res.dev.vec_7<-as.vector(abs(res.dev_7))
rank_deviance_7<-rank(-res.dev.vec_7,ties.method = "first") 
rank_dbht_7<-rank(dbht.order_7$pvalue,ties.method = "first")
rank_score_7<-rank(-(odS.matrix_7$odS.vec_7),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_7<-cbind(rank_deviance_7,rank_dbht_7,rank_score_7)
```

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

rankproduct<-RankProduct(rankMat_7)
```


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.
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.survival)[1]
k_7<-dim(rankMat_7)[2]
```

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

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

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

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


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. 

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

```{r}
qobj_7<-qvalue(pvalues_7)
```

```{r}
qvalues_7<-qobj_7$qvalues
which(qvalues_7<0.01)
which(qvalues_7<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
outliers_7<-as.data.frame(cbind(id,rank_deviance_7,rank_dbht_7,rank_score_7,pvalues_7,qvalues_7))
sort.outliers_7 <- outliers_7[order(qvalues_7) , ]
sort.outliers_7[1:10,]
```












## 8th scenario 

Notice that in all the scenarios the amount of outliers is 10% of the total of observations

- beta_pure   : 1
- beta_noise  : -1
- cens_amount : 0.3
- scale       : 1.5
- shape       : 0.5

```{r}
simulation_8<-get_simdata_2coxes_bybook(data.survival,10,1,-1,0.3,1.5,0.5)
data.time.status_8<-as.data.frame(simulation_8)
```


### Survival Analysis

```{r}
fit.cox_8<-coxph(Surv(data.time.status_8$times, data.time.status_8$cens)~.,data = data.survival)
fit.cox_8
```

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


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


The fitted model using exponential weights with 5% level of trimming is given bellow. Notice that the results show that all covariates are significant for a 5% level of significant.

```{r}
fit.coxrobust.H_8=rcoxph(data.time.status_8$times,data.time.status_8$cens,data.survival,wt.type="exponential",quant=.95)
fit.coxrobust.H_8$coefficients

```




```{r}
K_8=quantile(fit.coxrobust.H_8$Lambda.rob*exp(fit.coxrobust.H_8$xbeta.rob),probs=c(.95))
sort(-log(fit.coxrobust.H_8$wt)*K_8,decreasing = TRUE)
```



### Outlier detection methods for survival analisys
In order to identify which of the observations might be outliers some methods for outlier detection for survival analysis are going to be performed.

#### Residuals
The first method used for outlier detection is the deviance residuals.


```{r}
res.dev_8<- resid(fit.cox_8,type="deviance")
```


```{r, echo=FALSE}
par(mfrow=c(1,1))
plot(res.dev_8,main="Deviance Residuals",xlab="Case number")
abline(h=0)
```


```{r}
head(sort(abs(res.dev_8),decreasing = TRUE))
```


#### Concordance c-index

```{r, include=FALSE}
set.seed(1980)
n<-dim(data.survival)[1]
outliers_dbht_8<- bcsod(surv.object = Surv(data.time.status_8$times,data.time.status_8$cens),covariate.data =data.survival ,sod.method = "dbht", B=1000 , B.N = n)
```



Here are the top observations with lowest p-values.
```{r}
outliers_dbht_8[1:10,]
```


```{r}
dbht_8<-as.data.frame(outliers_dbht_8)# convert to a dataframe (more easy to handle)
names(dbht_8)[1]="id"
dbht.order_8<-dbht_8[order(dbht_8$id),]
```

#### Censored quantile regression

The score algorithm arises from the outlier detection algorithm for the Cox's regression model for censored data proposed by Nardi (1999) but now based on quantile regression.

```{r}
odcS_8<-odc(formula=Surv((data.time.status_8$times),data.time.status_8$cens)~ data.survival$V1+data.survival$V2+data.survival$V3, data = data.survival,method="score",rq.model="Wang",h=0.05)
odcS_8
```




The top outliers are:
```{r}
head(sort(odcS_8@score,decreasing = TRUE))
```



Put in form of a vector the results obtained before
```{r}
odS.vec_8<-as.vector(odcS_8@score)
```

Define a matrix with the id and the values for the residual algorithm.

```{r}
n<-dim(data.survival)[1]#dimension of the data
id<-as.vector(seq(1:n))# vector with the id
odS.matrix_8<-as.data.frame(cbind(id,odS.vec_8))
```



### Rank Product test

First a rank matrix is going to be obtained based on the ranks of each outliers detection method.

```{r}
res.dev.vec_8<-as.vector(abs(res.dev_8))
rank_deviance_8<-rank(-res.dev.vec_8,ties.method = "first") 
rank_dbht_8<-rank(dbht.order_8$pvalue,ties.method = "first")
rank_score_8<-rank(-(odS.matrix_8$odS.vec_8),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_8<-cbind(rank_deviance_8,rank_dbht_8,rank_score_8)
```

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

rankproduct<-RankProduct(rankMat_8)
```


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.
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.survival)[1]
k_8<-dim(rankMat_8)[2]
```

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

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

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

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


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. 

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

```{r}
qobj_8<-qvalue(pvalues_8)
```

```{r}
qvalues_8<-qobj_8$qvalues
which(qvalues_8<0.01)
which(qvalues_8<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
outliers_8<-as.data.frame(cbind(id,rank_deviance_8,rank_dbht_8,rank_score_8,pvalues_8,qvalues_8))
sort.outliers_8 <- outliers_8[order(qvalues_8) , ]
sort.outliers_8[1:10,]
```





## 9th scenario 

Notice that in all the scenarios the amount of outliers is 10% of the total of observations

- beta_pure   : 1
- beta_noise  : -1
- cens_amount : 0.3
- scale       : 0.5
- shape       : 1.5

```{r}
simulation_9<-get_simdata_2coxes_bybook(data.survival,10,1,-1,0.3,0.5,1.5)
data.time.status_9<-as.data.frame(simulation_9)
```


### Survival Analysis

```{r}
fit.cox_9<-coxph(Surv(data.time.status_9$times, data.time.status_9$cens)~.,data = data.survival)
fit.cox_9
```

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


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


The fitted model using exponential weights with 5% level of trimming is given bellow. Notice that the results show that all covariates are significant for a 5% level of significant.

```{r}
fit.coxrobust.H_9=rcoxph(data.time.status_9$times,data.time.status_9$cens,data.survival,wt.type="exponential",quant=.95)
fit.coxrobust.H_9$coefficients

```




```{r}
K_9=quantile(fit.coxrobust.H_9$Lambda.rob*exp(fit.coxrobust.H_9$xbeta.rob),probs=c(.95))
sort(-log(fit.coxrobust.H_9$wt)*K_9,decreasing = TRUE)
```



### Outlier detection methods for survival analisys
In order to identify which of the observations might be outliers some methods for outlier detection for survival analysis are going to be performed.

#### Residuals
The first method used for outlier detection is the deviance residuals.


```{r}
res.dev_9<- resid(fit.cox_9,type="deviance")
```


```{r, echo=FALSE}
par(mfrow=c(1,1))
plot(res.dev_9,main="Deviance Residuals",xlab="Case number")
abline(h=0)
```


```{r}
head(sort(abs(res.dev_9),decreasing = TRUE))
```


#### Concordance c-index

```{r, include=FALSE}
set.seed(1980)
n<-dim(data.survival)[1]
outliers_dbht_9<- bcsod(surv.object = Surv(data.time.status_9$times,data.time.status_9$cens),covariate.data =data.survival ,sod.method = "dbht", B=1000 , B.N = n)
```



Here are the top observations with lowest p-values.
```{r}
outliers_dbht_9[1:10,]
```


```{r}
dbht_9<-as.data.frame(outliers_dbht_9)# convert to a dataframe (more easy to handle)
names(dbht_9)[1]="id"
dbht.order_9<-dbht_9[order(dbht_9$id),]
```

#### Censored quantile regression

The score algorithm arises from the outlier detection algorithm for the Cox's regression model for censored data proposed by Nardi (1999) but now based on quantile regression.

```{r}
odcS_9<-odc(formula=Surv((data.time.status_9$times),data.time.status_9$cens)~ data.survival$V1+data.survival$V2+data.survival$V3, data = data.survival,method="score",rq.model="Wang",h=0.05)
odcS_9
```




The top outliers are:
```{r}
head(sort(odcS_9@score,decreasing = TRUE))
```



Put in form of a vector the results obtained before
```{r}
odS.vec_9<-as.vector(odcS_9@score)
```

Define a matrix with the id and the values for the residual algorithm.

```{r}
n<-dim(data.survival)[1]#dimension of the data
id<-as.vector(seq(1:n))# vector with the id
odS.matrix_9<-as.data.frame(cbind(id,odS.vec_9))
```



### Rank Product test

First a rank matrix is going to be obtained based on the ranks of each outliers detection method.

```{r}
res.dev.vec_9<-as.vector(abs(res.dev_9))
rank_deviance_9<-rank(-res.dev.vec_9,ties.method = "first") 
rank_dbht_9<-rank(dbht.order_9$pvalue,ties.method = "first")
rank_score_9<-rank(-(odS.matrix_9$odS.vec_9),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_9<-cbind(rank_deviance_9,rank_dbht_9,rank_score_9)
```

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

rankproduct<-RankProduct(rankMat_9)
```


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.
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.survival)[1]
k_9<-dim(rankMat_9)[2]
```

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

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

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

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


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. 

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

```{r}
qobj_9<-qvalue(pvalues_9)
```

```{r}
qvalues_9<-qobj_9$qvalues
which(qvalues_9<0.01)
which(qvalues_9<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
outliers_9<-as.data.frame(cbind(id,rank_deviance_9,rank_dbht_9,rank_score_9,pvalues_9,qvalues_9))
sort.outliers_9 <- outliers_9[order(qvalues_9) , ]
sort.outliers_9[1:10,]
```









## 10th scenario 

Notice that in all the scenarios the amount of outliers is 10% of the total of observations

- beta_pure   : 1
- beta_noise  : -0.2
- cens_amount : 0.3
- scale       : 1
- shape       : 1

```{r}
simulation_10<-get_simdata_2coxes_bybook(data.survival,10,1,-0.2,0.3,1,1)
data.time.status_10<-as.data.frame(simulation_10)
```


### Survival Analysis

```{r}
fit.cox_10<-coxph(Surv(data.time.status_10$times, data.time.status_10$cens)~.,data = data.survival)
fit.cox_10
```

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


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


The fitted model using exponential weights with 5% level of trimming is given bellow. Notice that the results show that all covariates are significant for a 5% level of significant.

```{r}
fit.coxrobust.H_10=rcoxph(data.time.status_10$times,data.time.status_10$cens,data.survival,wt.type="exponential",quant=.95)
fit.coxrobust.H_10$coefficients

```




```{r}
K_10=quantile(fit.coxrobust.H_10$Lambda.rob*exp(fit.coxrobust.H_10$xbeta.rob),probs=c(.95))
sort(-log(fit.coxrobust.H_10$wt)*K_10,decreasing = TRUE)
```



### Outlier detection methods for survival analisys
In order to identify which of the observations might be outliers some methods for outlier detection for survival analysis are going to be performed.

#### Residuals
The first method used for outlier detection is the deviance residuals.


```{r}
res.dev_10<- resid(fit.cox_10,type="deviance")
```


```{r, echo=FALSE}
par(mfrow=c(1,1))
plot(res.dev_10,main="Deviance Residuals",xlab="Case number")
abline(h=0)
```


```{r}
head(sort(abs(res.dev_10),decreasing = TRUE))
```


#### Concordance c-index

```{r, include=FALSE}
set.seed(1980)
n<-dim(data.survival)[1]
outliers_dbht_10<- bcsod(surv.object = Surv(data.time.status_10$times,data.time.status_10$cens),covariate.data =data.survival ,sod.method = "dbht", B=1000 , B.N = n)
```



Here are the top observations with lowest p-values.
```{r}
outliers_dbht_10[1:10,]
```


```{r}
dbht_10<-as.data.frame(outliers_dbht_10)# convert to a dataframe (more easy to handle)
names(dbht_10)[1]="id"
dbht.order_10<-dbht_10[order(dbht_10$id),]
```

#### Censored quantile regression

The score algorithm arises from the outlier detection algorithm for the Cox's regression model for censored data proposed by Nardi (1999) but now based on quantile regression.

```{r}
odcS_10<-odc(formula=Surv((data.time.status_10$times),data.time.status_10$cens)~ data.survival$V1+data.survival$V2+data.survival$V3, data = data.survival,method="score",rq.model="Wang",h=0.05)
odcS_10
```




The top outliers are:
```{r}
head(sort(odcS_10@score,decreasing = TRUE))
```



Put in form of a vector the results obtained before
```{r}
odS.vec_10<-as.vector(odcS_10@score)
```

Define a matrix with the id and the values for the residual algorithm.

```{r}
n<-dim(data.survival)[1]#dimension of the data
id<-as.vector(seq(1:n))# vector with the id
odS.matrix_10<-as.data.frame(cbind(id,odS.vec_10))
```



### Rank Product test

First a rank matrix is going to be obtained based on the ranks of each outliers detection method.

```{r}
res.dev.vec_10<-as.vector(abs(res.dev_10))
rank_deviance_10<-rank(-res.dev.vec_10,ties.method = "first") 
rank_dbht_10<-rank(dbht.order_10$pvalue,ties.method = "first")
rank_score_10<-rank(-(odS.matrix_10$odS.vec_10),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_10<-cbind(rank_deviance_10,rank_dbht_10,rank_score_10)
```

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

rankproduct<-RankProduct(rankMat_10)
```


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.
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.survival)[1]
k_10<-dim(rankMat_10)[2]
```

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

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

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

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


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. 

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

```{r}
qobj_10<-qvalue(pvalues_10)
```

```{r}
qvalues_10<-qobj_10$qvalues
which(qvalues_10<0.01)
which(qvalues_10<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
outliers_10<-as.data.frame(cbind(id,rank_deviance_10,rank_dbht_10,rank_score_10,pvalues_10,qvalues_10))
sort.outliers_10 <- outliers_10[order(qvalues_10) , ]
sort.outliers_10[1:10,]
```






## 11th scenario 

Notice that in all the scenarios the amount of outliers is 10% of the total of observations

- beta_pure   : 1
- beta_noise  : -0.2
- cens_amount : 0.3
- scale       : 1.5
- shape       : 0.5

```{r}
simulation_11<-get_simdata_2coxes_bybook(data.survival,10,1,-0.2,0.3,1.5,0.5)
data.time.status_11<-as.data.frame(simulation_11)
```


### Survival Analysis

```{r}
fit.cox_11<-coxph(Surv(data.time.status_11$times, data.time.status_11$cens)~.,data = data.survival)
fit.cox_11
```

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


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


The fitted model using exponential weights with 5% level of trimming is given bellow. Notice that the results show that all covariates are significant for a 5% level of significant.

```{r}
fit.coxrobust.H_11=rcoxph(data.time.status_11$times,data.time.status_11$cens,data.survival,wt.type="exponential",quant=.95)
fit.coxrobust.H_11$coefficients

```




```{r}
K_11=quantile(fit.coxrobust.H_11$Lambda.rob*exp(fit.coxrobust.H_11$xbeta.rob),probs=c(.95))
sort(-log(fit.coxrobust.H_11$wt)*K_11,decreasing = TRUE)
```



### Outlier detection methods for survival analisys
In order to identify which of the observations might be outliers some methods for outlier detection for survival analysis are going to be performed.

#### Residuals
The first method used for outlier detection is the deviance residuals.


```{r}
res.dev_11<- resid(fit.cox_11,type="deviance")
```


```{r, echo=FALSE}
par(mfrow=c(1,1))
plot(res.dev_11,main="Deviance Residuals",xlab="Case number")
abline(h=0)
```


```{r}
head(sort(abs(res.dev_11),decreasing = TRUE))
```


#### Concordance c-index

```{r, include=FALSE}
set.seed(1980)
n<-dim(data.survival)[1]
outliers_dbht_11<- bcsod(surv.object = Surv(data.time.status_11$times,data.time.status_11$cens),covariate.data =data.survival ,sod.method = "dbht", B=1000 , B.N = n)
```



Here are the top observations with lowest p-values.
```{r}
outliers_dbht_11[1:10,]
```


```{r}
dbht_11<-as.data.frame(outliers_dbht_11)# convert to a dataframe (more easy to handle)
names(dbht_11)[1]="id"
dbht.order_11<-dbht_11[order(dbht_11$id),]
```

#### Censored quantile regression

The score algorithm arises from the outlier detection algorithm for the Cox's regression model for censored data proposed by Nardi (1999) but now based on quantile regression.

```{r}
odcS_11<-odc(formula=Surv((data.time.status_11$times),data.time.status_11$cens)~ data.survival$V1+data.survival$V2+data.survival$V3, data = data.survival,method="score",rq.model="Wang",h=0.05)
odcS_11
```




The top outliers are:
```{r}
head(sort(odcS_11@score,decreasing = TRUE))
```



Put in form of a vector the results obtained before
```{r}
odS.vec_11<-as.vector(odcS_11@score)
```

Define a matrix with the id and the values for the residual algorithm.

```{r}
n<-dim(data.survival)[1]#dimension of the data
id<-as.vector(seq(1:n))# vector with the id
odS.matrix_11<-as.data.frame(cbind(id,odS.vec_11))
```



### Rank Product test

First a rank matrix is going to be obtained based on the ranks of each outliers detection method.

```{r}
res.dev.vec_11<-as.vector(abs(res.dev_11))
rank_deviance_11<-rank(-res.dev.vec_11,ties.method = "first") 
rank_dbht_11<-rank(dbht.order_11$pvalue,ties.method = "first")
rank_score_11<-rank(-(odS.matrix_11$odS.vec_11),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_11<-cbind(rank_deviance_11,rank_dbht_11,rank_score_11)
```

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

rankproduct<-RankProduct(rankMat_11)
```


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.
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.survival)[1]
k_11<-dim(rankMat_11)[2]
```

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

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

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

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


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. 

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

```{r}
qobj_11<-qvalue(pvalues_11)
```

```{r}
qvalues_11<-qobj_11$qvalues
which(qvalues_11<0.01)
which(qvalues_11<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
outliers_11<-as.data.frame(cbind(id,rank_deviance_11,rank_dbht_11,rank_score_11,pvalues_11,qvalues_11))
sort.outliers_11 <- outliers_11[order(qvalues_11) , ]
sort.outliers_11[1:10,]
```







## 12th scenario 

Notice that in all the scenarios the amount of outliers is 10% of the total of observations

- beta_pure   : 1
- beta_noise  : -0.2
- cens_amount : 0.3
- scale       : 0.5
- shape       : 1.5

```{r}
simulation_12<-get_simdata_2coxes_bybook(data.survival,10,1,-0.2,0.3,0.5,1.5)
data.time.status_12<-as.data.frame(simulation_12)
```


### Survival Analysis

```{r}
fit.cox_12<-coxph(Surv(data.time.status_12$times, data.time.status_12$cens)~.,data = data.survival)
fit.cox_12
```

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


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


The fitted model using exponential weights with 5% level of trimming is given bellow. Notice that the results show that all covariates are significant for a 5% level of significant.

```{r}
fit.coxrobust.H_12=rcoxph(data.time.status_12$times,data.time.status_12$cens,data.survival,wt.type="exponential",quant=.95)
fit.coxrobust.H_12$coefficients

```




```{r}
K_12=quantile(fit.coxrobust.H_12$Lambda.rob*exp(fit.coxrobust.H_12$xbeta.rob),probs=c(.95))
sort(-log(fit.coxrobust.H_12$wt)*K_12,decreasing = TRUE)
```



### Outlier detection methods for survival analisys
In order to identify which of the observations might be outliers some methods for outlier detection for survival analysis are going to be performed.

#### Residuals
The first method used for outlier detection is the deviance residuals.


```{r}
res.dev_12<- resid(fit.cox_12,type="deviance")
```


```{r, echo=FALSE}
par(mfrow=c(1,1))
plot(res.dev_12,main="Deviance Residuals",xlab="Case number")
abline(h=0)
```


```{r}
head(sort(abs(res.dev_12),decreasing = TRUE))
```


#### Concordance c-index

```{r, include=FALSE}
set.seed(1980)
n<-dim(data.survival)[1]
outliers_dbht_12<- bcsod(surv.object = Surv(data.time.status_12$times,data.time.status_12$cens),covariate.data =data.survival ,sod.method = "dbht", B=1000 , B.N = n)
```



Here are the top observations with lowest p-values.
```{r}
outliers_dbht_12[1:10,]
```


```{r}
dbht_12<-as.data.frame(outliers_dbht_12)# convert to a dataframe (more easy to handle)
names(dbht_12)[1]="id"
dbht.order_12<-dbht_12[order(dbht_12$id),]
```

#### Censored quantile regression

The score algorithm arises from the outlier detection algorithm for the Cox's regression model for censored data proposed by Nardi (1999) but now based on quantile regression.

```{r}
odcS_12<-odc(formula=Surv((data.time.status_12$times),data.time.status_12$cens)~ data.survival$V1+data.survival$V2+data.survival$V3, data = data.survival,method="score",rq.model="Wang",h=0.05)
odcS_12
```




The top outliers are:
```{r}
head(sort(odcS_12@score,decreasing = TRUE))
```



Put in form of a vector the results obtained before
```{r}
odS.vec_12<-as.vector(odcS_12@score)
```

Define a matrix with the id and the values for the residual algorithm.

```{r}
n<-dim(data.survival)[1]#dimension of the data
id<-as.vector(seq(1:n))# vector with the id
odS.matrix_12<-as.data.frame(cbind(id,odS.vec_12))
```



### Rank Product test

First a rank matrix is going to be obtained based on the ranks of each outliers detection method.

```{r}
res.dev.vec_12<-as.vector(abs(res.dev_12))
rank_deviance_12<-rank(-res.dev.vec_12,ties.method = "first") 
rank_dbht_12<-rank(dbht.order_12$pvalue,ties.method = "first")
rank_score_12<-rank(-(odS.matrix_12$odS.vec_12),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_12<-cbind(rank_deviance_12,rank_dbht_12,rank_score_12)
```

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

rankproduct<-RankProduct(rankMat_12)
```


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.
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.survival)[1]
k_12<-dim(rankMat_12)[2]
```

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

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

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

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


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. 

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

```{r}
qobj_12<-qvalue(pvalues_12)
```

```{r}
qvalues_12<-qobj_12$qvalues
which(qvalues_12<0.01)
which(qvalues_12<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
outliers_12<-as.data.frame(cbind(id,rank_deviance_12,rank_dbht_12,rank_score_12,pvalues_12,qvalues_12))
sort.outliers_12 <- outliers_12[order(qvalues_12) , ]
sort.outliers_12[1:10,]
```


