%\VignetteIndexEntry{genefu An Introduction (HowTo)}
%\VignetteDepends{survcomp, mclust, rmeta, Biobase, xtable}
%\VignetteSuggests{GeneMeta, breastCancerVDX, breastCancerMAINZ, breastCancerTRANSBIG, breastCancerUPP, breastCancerUNT, breastCancerNKI}
%\VignetteImports{amap}
%\VignetteKeywords{Breast Cancer, Survival Analysis}
%\VignettePackage{genefu}

\documentclass[a4paper,11pt]{article}

\usepackage{amsmath}
\usepackage{times}
\usepackage{hyperref}
\usepackage[numbers]{natbib}
\usepackage[american]{babel}
\usepackage{authblk}
\renewcommand\Affilfont{\itshape\small}
\usepackage{Sweave}
\renewcommand{\topfraction}{0.85}
\renewcommand{\textfraction}{0.1}
\usepackage{graphicx}
\usepackage{tikz}


\textwidth=6.2in
\textheight=8.5in
%\parskip=.3cm
\oddsidemargin=.1in
\evensidemargin=.1in
\headheight=-.3in

%------------------------------------------------------------
% newcommand
%------------------------------------------------------------
\newcommand{\scscst}{\scriptscriptstyle}
\newcommand{\scst}{\scriptstyle}
\newcommand{\Robject}[1]{{\texttt{#1}}}
\newcommand{\Rfunction}[1]{{\texttt{#1}}}
\newcommand{\Rclass}[1]{\textit{#1}}
\newcommand{\Rpackage}[1]{\textit{#1}}
\newcommand{\Rexpression}[1]{\texttt{#1}}
\newcommand{\Rmethod}[1]{{\texttt{#1}}}
\newcommand{\Rfunarg}[1]{{\texttt{#1}}}


\begin{document}

%------------------------------------------------------------
\title{\Rpackage{genefu}: a package for breast cancer gene expression analysis}
%------------------------------------------------------------
\author[1]{Benjamin Haibe-Kains}
\author[2]{Markus Schr\"{o}der}
\author[3]{Christos Sotiriou}
\author[4]{Gianluca Bontempi}
\author[5,6]{John Quackenbush}


\affil[1]{Bioinformatics and Computational Genomics Laboratory, Princess Margaret Cancer Center, University Health Network, Toronto, Ontario, Canada}
\affil[2]{UCD School of Biomolecular and Biomedical Science, Conway Institute, University College Dublin, Belfield, Dublin, Ireland}
\affil[3]{Breast Cancer Translational Research Laboratory, Institut Jules Bordet, Universit\'{e} Libre de Bruxelles}
\affil[4]{Machine Learning Group, Universit\'{e} Libre de Bruxelles}
\affil[5]{Computational Biology and Functional Genomics Laboratory, Dana-Farber Cancer Institute, Harvard School of Public Health}
\affil[6]{Center for Cancer Computational Biology, Dana-Farber Cancer Institute}

\SweaveOpts{highlight=TRUE, tidy=TRUE, keep.space=TRUE, keep.blank.space=FALSE, keep.comment=TRUE}

%<<setup,echo=FALSE>>=
%library(pgfSweave)
%setCacheDir("cache")
%options(keep.source=TRUE)
%@

\maketitle
\tableofcontents

%------------------------------------------------------------
\section{Introduction}
%------------------------------------------------------------ 

The \Rpackage{genefu} package is providing relevant functions for gene expression analysis, especially in breast cancer.





%------------------------------------------------------------
\section{Case Study 1: comparing risk prediction models}
%------------------------------------------------------------ 
For computing the risk scores, estimates of the performance of the risk scores, combining the estimates and comparing the estimates we have to load the \Rpackage{genefu} and \Rpackage{survcomp} packages into the workspace. We also load the \Rpackage{xtable} package to display results insight this document.
<<loadPackages,results=hide>>=
library(genefu)
library(xtable)
library(rmeta)
@

The five data sets that we use in the case study are publicly available as experimental data packages on Bioconductor.org. In particular we used:

\begin{description}
\item[breastCancerMAINZ:]  \url{http://www.bioconductor.org/packages/release/data/experiment/html/breastCancerMAINZ.html}\\
\item[breastCancerUPP:]  \url{http://www.bioconductor.org/packages/release/data/experiment/html/breastCancerUPP.html}\\
\item[breastCancerUNT:]  \url{http://www.bioconductor.org/packages/release/data/experiment/html/breastCancerUNT.html}\\
\item[breastCancerNKI:]  \url{http://www.bioconductor.org/packages/release/data/experiment/html/breastCancerNKI.html}\\
\item[breastCancerTRANSBIG:]  \url{http://www.bioconductor.org/packages/release/data/experiment/html/breastCancerTRANSBIG.html}\\
\end{description}

We don't use the \Rpackage{breastCancerVDX} experimental package in the case study since it has been used as training data set for GENIUS \citep{HaibeKains2010}:

\begin{description}
\item[breastCancerVDX:]  \url{http://www.bioconductor.org/packages/release/data/experiment/html/breastCancerVDX.html}\\
\end{description}

These experimental data packages can be installed from Bioconductor version 2.8 or higher in R version 2.13.0 or higher. For the experimental data packages the commands for installing the data sets are:

<<installAndLoadMAINZ,results=hide,eval=FALSE>>=
source("http://www.bioconductor.org/biocLite.R")
biocLite("breastCancerMAINZ")
biocLite("breastCancerTRANSBIG")
biocLite("breastCancerUPP")
biocLite("breastCancerUNT")
biocLite("breastCancerNKI")
@

And for loading the data sets into the current workspace the commands are:
<<installAndLoadAllPackages,results=hide>>=
library(breastCancerMAINZ)
library(breastCancerTRANSBIG)
library(breastCancerUPP)
library(breastCancerUNT)
library(breastCancerNKI)
@

Table \ref{tab1} shows an overview of the data sets and the patients. From those 1123 breast cancer patients we selected only the patients that are node negative and didn't receive any treatment (except local radiotherapy), which results in 722 patients.

\begin{table}
\begin{center}
\caption{Detailed overview for the data sets used in the case study}
\label{tab1}
\begin{tabular}{ l | l*{4}{c}r }
Dataset & Patients [\#] & ER+ [\#] & HER2+ [\#] & Age [years] & Grade [1/2/3] & Platform \\ \hline
MAINZ & 200 & 155 & 23 & 25-90 & 29/136/35 & HGU133A \\
TRANSBIG & 198 & 123 & 35 & 24-60 & 30/83/83 & HGU133A \\
UPP & 251 & 175 & 46 & 28-93 & 67/128/54 & HGU133AB \\
UNT & 137 & 94 & 21 & 24-73 & 32/51/29 & HGU133AB \\
NKI & 337 & 212 & 53 & 26-62 & 79/109/149 & Agilent \\ \hline
Overall & 1123 & 759 & 178 & 25-73 & 237/507/350 & Affy/Agilent \\
\end{tabular}
\end{center}
\end{table}

%------------------------------------------------------------
%------------------------------------------------------------
%------------------------------------------------------------
%------------------------------------------------------------
%%% CODE START
%------------------------------------------------------------
%------------------------------------------------------------
%------------------------------------------------------------

Since there are duplicated patients in the five data sets, we have to identify the duplicated patients and we subsequently store them in a vector.
<<findDuplicatedPatients,results=hide, echo=FALSE>>=
library(Biobase)
data(breastCancerData)
cinfo <- colnames(pData(mainz7g))
data.all <- c("transbig7g"=transbig7g, "unt7g"=unt7g, "upp7g"=upp7g, "mainz7g"=mainz7g, "nki7g"=nki7g)

idtoremove.all <- NULL
duplres <- NULL

## UNT vs UPP vs TRANSBIG
demo.all <- rbind(pData(transbig7g), pData(unt7g), pData(upp7g))
dn2 <- c("TRANSBIG", "UNT", "UPP")

## Karolinska
## VDXKIU, KIU, UPPU
ds2 <- c("VDXKIU", "KIU", "UPPU")
demot <- demo.all[complete.cases(demo.all[ , c("series")]) & is.element(demo.all[ , "series"], ds2), ]

duplid <- sort(unique(demot[duplicated(demot[ , "id"]), "id"]))
duplrest <- NULL
for(i in 1:length(duplid)) {
  tt <- NULL
  for(k in 1:length(dn2)) {
    myx <- sort(row.names(demot)[complete.cases(demot[ , c("id", "dataset")]) & demot[ , "id"] == duplid[i] & demot[ , "dataset"] == dn2[k]])
    if(length(myx) > 0) { tt <- c(tt, myx) }
  }
  duplrest <- c(duplrest, list(tt))
}
names(duplrest) <- duplid
duplres <- c(duplres, duplrest)

## Oxford
## VVDXOXFU, OXFU
ds2 <- c("VDXOXFU", "OXFU")
demot <- demo.all[complete.cases(demo.all[ , c("series")]) & is.element(demo.all[ , "series"], ds2), ]

duplid <- sort(unique(demot[duplicated(demot[ , "id"]), "id"]))
duplrest <- NULL
for(i in 1:length(duplid)) {
  tt <- NULL
  for(k in 1:length(dn2)) {
    myx <- sort(row.names(demot)[complete.cases(demot[ , c("id", "dataset")]) & demot[ , "id"] == duplid[i] & demot[ , "dataset"] == dn2[k]])
    if(length(myx) > 0) { tt <- c(tt, myx) }
  }
  duplrest <- c(duplrest, list(tt))
}
names(duplrest) <- duplid
duplres <- c(duplres, duplrest)

## duplicated patients
duPL <- sort(unlist(lapply(duplres, function(x) { return(x[-1]) } )))
@
We then compute the risk scores for NPI, AURKA, GGI and GENIUS with the \Rfunction{npi()}, \Rfunction{sig.score()}, \Rfunction{ggi()} and \Rfunction{genius()} functions within \Rpackage{genefu}, respectively.
<<computeRiskScore>>=
dn <- c("transbig", "unt", "upp", "mainz", "nki")
dn.platform <- c("affy", "affy", "affy", "affy", "agilent")

res <- ddemo.all <- ddemo.coln <- NULL
for(i in 1:length(dn)) {

  ## load dataset
  dd <- get(data(list=dn[i]))

  ddata <- t(exprs(dd))
  ddemo <- phenoData(dd)@data
  dannot <- featureData(dd)@data
  ddemo.all <- c(ddemo.all, list(ddemo))
  if(is.null(ddemo.coln)) { ddemo.coln <- colnames(ddemo) } else { ddemo.coln <- intersect(ddemo.coln, colnames(ddemo)) }
  rest <- NULL

  ## NPI
  ss <- ddemo[ , "size"]
  gg <- ddemo[ , "grade"]
  nn <- rep(NA, nrow(ddemo))
  nn[complete.cases(ddemo[ , "node"]) & ddemo[ , "node"] == 0] <- 1 
  nn[complete.cases(ddemo[ , "node"]) & ddemo[ , "node"] == 1] <- 3
  names(ss) <- names(gg) <- names(nn) <- rownames(ddemo)
  rest <- cbind(rest, "NPI"=npi(size=ss, grade=gg, node=nn, na.rm=TRUE)$score)

  ## AURKA
  ## if affy platform consider the probe published in Desmedt et al., CCR, 2008
  if(dn.platform[i] == "affy") { domap <- FALSE } else { domap <- TRUE }
  modt <- scmgene.robust$mod$AURKA
  ## if agilent platform consider the probe published in Desmedt et al., CCR, 2008
  if(dn.platform[i] == "agilent") {
    domap <- FALSE
    modt[ , "probe"] <- "NM_003600"
  }
  rest <- cbind(rest, "AURKA"=sig.score(x=modt, data=ddata, annot=dannot, do.mapping=domap)$score)

  ## GGI
  if(dn.platform[i] == "affy") { domap <- FALSE } else { domap <- TRUE }
  rest <- cbind(rest, "GGI"=ggi(data=ddata, annot=dannot, do.mapping=domap)$score)
  ## GENIUS
  if(dn.platform[i] == "affy") { domap <- FALSE } else { domap <- TRUE }
  rest <- cbind(rest, "GENIUS"=genius(data=ddata, annot=dannot, do.mapping=domap)$score)
  res <- rbind(res, rest)
}
names(ddemo.all) <- dn
@

For further analysis and handling of the data we store all information in one object. We also remove the duplicated patients from the analysis and take only those patients into account, that have complete information for nodal, survival and treatment status.
<<simplifyAndRemoveDuplicatePatients>>=
ddemot <- NULL
for(i in 1:length(ddemo.all)) {
  ddemot <- rbind(ddemot, ddemo.all[[i]][ , ddemo.coln, drop=FALSE])
}
res[complete.cases(ddemot[ ,"dataset"]) & ddemot[ ,"dataset"] == "VDX", "GENIUS"] <- NA

## select only untreated node-negative patients with all risk predictions
myx <- complete.cases(res, ddemot[ , c("node", "treatment")]) & ddemot[ , "treatment"] == 0 & ddemot[ , "node"] == 0 & !is.element(rownames(ddemot), duPL)

res <- res[myx, , drop=FALSE]
ddemot <- ddemot[myx, , drop=FALSE]
@

To compare the risk score performances, we compute the concordance index\footnote{The same analysis could be performed with D index and hazard ratio by using the functions \Rfunction{D.index} and \Rfunction{hazard.ratio} from the \Rpackage{survcomp} respectively}:
<<cindexComputation>>=
cc.res <- complete.cases(res)
datasetList <- c("MAINZ","TRANSBIG","UPP","UNT","NKI")
riskPList <- c("NPI", "AURKA", "GGI", "GENIUS")
setT <- setE <- NULL
resMatrix <- as.list(NULL)

for(i in datasetList){
  dataset.only <- ddemot[,"dataset"] == i
  patientsAll <- cc.res & dataset.only
	
  ## set type of available survival data
  if(i == "UPP") {
    setT <- "t.rfs"
    setE <- "e.rfs"
  } else {
    setT <- "t.dmfs"
    setE <- "e.dmfs"
  }

  ## cindex computation
  cindexN <- t(apply(X=t(res[patientsAll,"NPI"]), MARGIN=1, function(x, y, z) {
    tt <- concordance.index(x=x, surv.time=y, surv.event=z, method="noether", na.rm=TRUE);
    return(c("cindex"=tt$c.index, "cindex.se"=tt$se, "lower"=tt$lower, "upper"=tt$upper)); },
    y=ddemot[patientsAll,setT], z=ddemot[patientsAll, setE]))
	
  cindexA <- t(apply(X=t(res[patientsAll,"AURKA"]), MARGIN=1, function(x, y, z) {
    tt <- concordance.index(x=x, surv.time=y, surv.event=z, method="noether", na.rm=TRUE);
    return(c("cindex"=tt$c.index, "cindex.se"=tt$se, "lower"=tt$lower, "upper"=tt$upper)); },
    y=ddemot[patientsAll,setT], z=ddemot[patientsAll, setE]))
	
  cindexM <- t(apply(X=t(res[patientsAll,"GGI"]), MARGIN=1, function(x, y, z) {
    tt <- concordance.index(x=x, surv.time=y, surv.event=z, method="noether", na.rm=TRUE);
    return(c("cindex"=tt$c.index, "cindex.se"=tt$se, "lower"=tt$lower, "upper"=tt$upper)); },
    y=ddemot[patientsAll, setT], z=ddemot[patientsAll, setE]))
	
  cindexG <- t(apply(X=t(res[patientsAll,"GENIUS"]), MARGIN=1, function(x, y, z) {
    tt <- concordance.index(x=x, surv.time=y, surv.event=z, method="noether", na.rm=TRUE);
    return(c("cindex"=tt$c.index, "cindex.se"=tt$se, "lower"=tt$lower, "upper"=tt$upper)); },
    y=ddemot[patientsAll, setT], z=ddemot[patientsAll, setE]))

	
  resMatrix[["NPI"]] <- rbind(resMatrix[["NPI"]], cindexN)
  resMatrix[["AURKA"]] <- rbind(resMatrix[["AURKA"]], cindexA)
  resMatrix[["GGI"]] <- rbind(resMatrix[["GGI"]], cindexM)
  resMatrix[["GENIUS"]] <- rbind(resMatrix[["GENIUS"]], cindexG)
}
@


Using a random-effects model we combine the dataset-specific performance estimated into overall estimates for each risk prediction model:
<<combineEstimations>>=
for(i in names(resMatrix)){
  ceData <- combine.est(x=resMatrix[[i]][,"cindex"], x.se=resMatrix[[i]][,"cindex.se"], hetero=TRUE)
  cLower <- ceData$estimate + qnorm(0.025, lower.tail=TRUE) * ceData$se
  cUpper <- ceData$estimate + qnorm(0.025, lower.tail=FALSE) * ceData$se
	
  cindexO <- cbind("cindex"=ceData$estimate, "cindex.se"=ceData$se, "lower"=cLower, "upper"=cUpper)
  resMatrix[[i]] <- rbind(resMatrix[[i]], cindexO)
  rownames(resMatrix[[i]]) <- c(datasetList, "Overall")
}
@


In order to compare the different risk prediction models we  compute one-sided p-values of the meta-estimates:
<<computePValues>>=
pv <- sapply(resMatrix, function(x) { return(x["Overall", c("cindex","cindex.se")]) })
pv <- apply(pv, 2, function(x) { return(pnorm((x[1] - 0.5) / x[2], lower.tail=x[1] < 0.5)) })
printPV <- matrix(pv,ncol=4)
rownames(printPV) <- "P-value"
colnames(printPV) <- names(pv)
@

<<printPvalue,results=tex>>==
xtable(printPV, digits=c(0, rep(-1,ncol(printPV))))
@

The following five figures represent the risk score performances measured by the concordance index for NPI, AURKA, GGI and GENIUS. The last figure represents the overall estimates.
<<forestplotNPI,fig=TRUE>>=
par(mfrow=c(2,2))
datasetListF <- c("MAINZ","TRANSBIG","UPP","UNT","NKI", NA, "Overall")
myspace <- "   "

## NPI Forestplot
tt <- rbind(resMatrix[["NPI"]][1:5,],
          "NA"=NA,
          "Overall"=resMatrix[["NPI"]][6,])

tt <- as.data.frame(tt)
labeltext <- cbind(c("Dataset", datasetListF))

r.mean <- c(NA,tt$cindex)
r.lower <- c(NA,tt$lower)
r.upper <- c(NA,tt$upper)

metaplot.surv(mn=r.mean, lower=r.lower, upper=r.upper, labels=labeltext, xlim=c(0.4,0.9), boxsize=0.5, zero=0.5, col=meta.colors(box="royalblue",line="darkblue",zero="firebrick"), main="NPI Concordance Index")
#@
#
#<<forestplotAURKA,fig=TRUE,echo=FALSE>>=
## AURKA Forestplot
tt <- rbind(resMatrix[["AURKA"]][1:5,],
          "NA"=NA,
          "Overall"=resMatrix[["AURKA"]][6,])

tt <- as.data.frame(tt)
labeltext <- cbind(c("Dataset", datasetListF))

r.mean <- c(NA,tt$cindex)
r.lower <- c(NA,tt$lower)
r.upper <- c(NA,tt$upper)

metaplot.surv(mn=r.mean, lower=r.lower, upper=r.upper, labels=labeltext, xlim=c(0.4,0.9), boxsize=0.5, zero=0.5, col=meta.colors(box="royalblue",line="darkblue",zero="firebrick"), main="AURKA Concordance Index")
#@
#
#<<forestplotGGI,fig=TRUE,echo=FALSE>>=
## GGI Forestplot
tt <- rbind(resMatrix[["GGI"]][1:5,],
          "NA"=NA,
          "Overall"=resMatrix[["GGI"]][6,])

tt <- as.data.frame(tt)
labeltext <- cbind(c("Dataset", datasetListF))

r.mean <- c(NA,tt$cindex)
r.lower <- c(NA,tt$lower)
r.upper <- c(NA,tt$upper)

metaplot.surv(mn=r.mean, lower=r.lower, upper=r.upper, labels=labeltext, xlim=c(0.4,0.9), boxsize=0.5, zero=0.5, col=meta.colors(box="royalblue",line="darkblue",zero="firebrick"), main="GGI Concordance Index")
#@
#
#<<forestplotGENIUS,fig=TRUE,echo=FALSE>>=
## GENIUS Forestplot
tt <- rbind(resMatrix[["GENIUS"]][1:5,],
          "NA"=NA,
          "Overall"=resMatrix[["GENIUS"]][6,])

tt <- as.data.frame(tt)
labeltext <- cbind(c("Dataset", datasetListF))

r.mean <- c(NA,tt$cindex)
r.lower <- c(NA,tt$lower)
r.upper <- c(NA,tt$upper)

metaplot.surv(mn=r.mean, lower=r.lower, upper=r.upper, labels=labeltext, xlim=c(0.4,0.9), boxsize=0.5, zero=0.5, col=meta.colors(box="royalblue",line="darkblue",zero="firebrick"), main="GENIUS Concordance Index")
@

<<forestplotOverall,fig=TRUE,height=5,echo=FALSE>>=
## Overall Forestplot
mybigspace <- "       "
tt <- rbind("OverallN"=resMatrix[["NPI"]][6,],
          "OverallA"=resMatrix[["AURKA"]][6,],
          "OverallM"=resMatrix[["GGI"]][6,],
          "OverallG"=resMatrix[["GENIUS"]][6,])

tt <- as.data.frame(tt)
labeltext <- cbind(c("Risk Prediction","NPI","AURKA","GGI","GENIUS"))

r.mean <- c(NA,tt$cindex)
r.lower <- c(NA,tt$lower)
r.upper <- c(NA,tt$upper)

metaplot.surv(mn=r.mean, lower=r.lower, upper=r.upper, labels=labeltext, xlim=c(0.45,0.75), boxsize=0.5, zero=0.5, col=meta.colors(box="royalblue",line="darkblue",zero="firebrick"), main="Overall Concordance Index")
@

In order to assess the difference between the risk scores, we compute the concordance indices with their p-values and compare the estimates with the \Rfunction{cindex.comp.meta} with a paired student t test.
<<computeCindexWithPvalue,echo=FALSE>>=
cc.res <- complete.cases(res)
datasetList <- c("MAINZ","TRANSBIG","UPP","UNT","NKI")
riskPList <- c("NPI", "AURKA", "GGI", "GENIUS")
setT <- setE <- NULL
resMatrixFull <- as.list(NULL)

for(i in datasetList){
  dataset.only <- ddemot[,"dataset"] == i
  patientsAll <- cc.res & dataset.only
	
  ## set type of available survival data
  if(i == "UPP") {
    setT <- "t.rfs"
    setE <- "e.rfs"
  } else {
    setT <- "t.dmfs"
    setE <- "e.dmfs"
  }
	
  ## cindex and p-value computation
  cindexN <- t(apply(X=t(res[patientsAll,"NPI"]), MARGIN=1, function(x, y, z) {
    tt <- concordance.index(x=x, surv.time=y, surv.event=z, method="noether", na.rm=TRUE); return(tt); },
    y=ddemot[patientsAll,setT], z=ddemot[patientsAll, setE]))
	
  cindexA <- t(apply(X=t(res[patientsAll,"AURKA"]), MARGIN=1, function(x, y, z) {
    tt <- concordance.index(x=x, surv.time=y, surv.event=z, method="noether", na.rm=TRUE);
    return(tt); },
    y=ddemot[patientsAll,setT], z=ddemot[patientsAll, setE]))
	
  cindexM <- t(apply(X=t(res[patientsAll,"GGI"]), MARGIN=1, function(x, y, z) {
    tt <- concordance.index(x=x, surv.time=y, surv.event=z, method="noether", na.rm=TRUE);
    return(tt); },
    y=ddemot[patientsAll, setT], z=ddemot[patientsAll, setE]))
	
  cindexG <- t(apply(X=t(res[patientsAll,"GENIUS"]), MARGIN=1, function(x, y, z) {
    tt <- concordance.index(x=x, surv.time=y, surv.event=z, method="noether", na.rm=TRUE);
    return(tt); },
    y=ddemot[patientsAll, setT], z=ddemot[patientsAll, setE]))
	
  resMatrixFull[["NPI"]] <- rbind(resMatrixFull[["NPI"]], cindexN)
  resMatrixFull[["AURKA"]] <- rbind(resMatrixFull[["AURKA"]], cindexA)
  resMatrixFull[["GGI"]] <- rbind(resMatrixFull[["GGI"]], cindexM)
  resMatrixFull[["GENIUS"]] <- rbind(resMatrixFull[["GENIUS"]], cindexG)
}

for(i in names(resMatrixFull)){
  rownames(resMatrixFull[[i]]) <- datasetList
}

ccmData <- tt <- rr <- NULL
for(i in 1:length(resMatrixFull)){
  tt <- NULL
  for(j in 1:length(resMatrixFull)){
    if(i != j) { rr <- cindex.comp.meta(list.cindex1=resMatrixFull[[i]], list.cindex2=resMatrixFull[[j]], hetero=TRUE)$p.value } else { rr <- 1 }
    tt <- cbind(tt, rr)
  }
  ccmData <- rbind(ccmData, tt)
}
ccmData <- as.data.frame(ccmData)
colnames(ccmData) <- riskPList
rownames(ccmData) <- riskPList
@

Table 2 displays the for multiple testing uncorrected p-values for the comparison of the different methods:
<<printCCM,results=tex>>=
xtable(ccmData, digits=c(0, rep(-1,ncol(ccmData))))
@

We correct the p-value with Holms method:
<<computeCCMPval>>=
ccmDataPval <- matrix(p.adjust(data.matrix(ccmData), method="holm"),ncol=4,dimnames=list(rownames(ccmData),colnames(ccmData)))
@

Table 3 displays the corrected p-values:
<<printCCMPval,results=tex>>=
xtable(ccmDataPval, digits=c(0, rep(-1,ncol(ccmDataPval))))
@


%------------------------------------------------------------
%------------------------------------------------------------
%------------------------------------------------------------
%------------------------------------------------------------
%%% CODE Stop
%------------------------------------------------------------
%------------------------------------------------------------
%------------------------------------------------------------


\newpage
%------------------------------------------------------------
\section{Session Info}
%------------------------------------------------------------ 
<<sessionInfo,echo=FALSE,results=tex>>==
toLatex(sessionInfo())
@

\end{document}