
\documentclass{article}

\usepackage{graphics}

%\VignetteIndexEntry{pensim package user guide}

\begin{document}

<<foo,echo=FALSE>>=
options(keep.source = TRUE)
options(width = 60)
options(eps = FALSE)
foo <- packageDescription("pensim")
@


\title{pensim Package Example (Version \Sexpr{foo$Version})}
\author{Levi Waldron}
\maketitle
\tableofcontents 

\section{Introduction}

This package acts as a wrapper to the \emph{penalized} R package to
add the following functionality to that package:

\begin{itemize}
\item repeated tuning on different data folds, with parallelization to
  take advantage of multiple processors
\item two-dimensional tuning of the Elastic Net
\item calculation of cross-validated risk scores through a nested
  cross-validation procedure, with parallelization to take advantage
  of multiple processors
\end{itemize}

It also provides a function for simulation of collinear
high-dimensional data with survival or binary response.

This package was developed in support of the study by Waldron et
al.\cite{waldron_optimized_2011}.  This paper contains greater detail
on proper application of the methods provided here.  Please cite this
paper when using the pensim package in your research, as well as the
penalized package\cite{goeman_l1_2010}.

\section{Example data}

pensim provides example data from a microarray experiment
investigating survival of cancer patients with lung adenocarcinomas
(Beer et al., 2002)\cite{beer_gene-expression_2002}.  Load this data
and do an initial pre-filter of genes with low IQR:

@ 
<<dataprep>>=
library(pensim)
data(beer.exprs)
data(beer.survival)
##select just 100 genes to speed computation, just for the sake of example:
set.seed(1)
beer.exprs.sample <- beer.exprs[sample(1:nrow(beer.exprs), 100), ]
#
gene.quant <- apply(beer.exprs.sample, 1, quantile, probs=0.75)
dat.filt <- beer.exprs.sample[gene.quant>log2(100), ]
gene.iqr <- apply(dat.filt, 1, IQR)
dat.filt <- as.matrix(dat.filt[gene.iqr>0.5, ])
dat.filt <- t(dat.filt)
dat.filt <- data.frame(dat.filt)
#
library(survival)
surv.obj <- Surv(beer.survival$os, beer.survival$status)
@ %def 

Note that the expression data are in ``wide'' format, with one column
per predictor (gene).  It is recommended to put covariate data in a
dataframe object, rather than a matrix.

\section{Nested cross-validation}

Unbiased estimation of prediction accuracy involves two levels of
cross-validation: an outer level for estimating prediction accuracy,
and an inner level for model tuning.  This process is simplified by
the opt.nested.crossval function.


It is recommended first to establish the arguments for the penalized
regression by testing on the \emph{penalized} package functions optL1
(for LASSO), optL2 (for Ridge) or cvl (for Elastic Net).  Here we use
LASSO.  Setting maxlambda1=5 is not a generally recommended procedure,
but is useful in this toy example to avoid converging on the null
model.

@ 
<<lassotest>>=
library(penalized)
testfit <- optL1(response=surv.obj,
                 penalized=dat.filt,
                 fold=5,
                 maxlambda1=5,
                 positive=FALSE,
                 standardize=TRUE,
                 trace=FALSE)
@ %def 


Now pass these arguments to opt.nested.crossval() for cross-validated
calculation and assessment of risk scores, with the additional
arguments:

\begin{itemize}
\item outerfold and nprocessors: number of folds for the outer level of cross-validation, and the number of processors to use for the outer level of cross-validation  (see ?opt.nested.crossval)
\item optFUN and scaling: which optimization function to use (opt1D for LASSO or Ridge, opt2D for Elastic Net) - see ?opt.splitval.
\item   setpen and nsim: setpen defines whether to do LASSO (``L1'') or Ridge (``L2''), nsim defines the number of times to repeat tuning (see ?opt1D.  opt2D has different required arguments.)
\end{itemize}

@ 
<<opt.nested.crossval>>=
set.seed(1)
preds <- opt.nested.crossval(outerfold=5, nprocessors=1, #opt.nested.crossval arguments
             optFUN="opt1D", scaling=FALSE,                     #opt.splitval arguments
             setpen="L1", nsim=1,                                      #opt1D arguments
             response=surv.obj,                    #rest are penalized::optl1 arguments
             penalized=dat.filt,
             fold=5,
             positive=FALSE,
             standardize=TRUE,
             trace=FALSE)
@ %def 

Ideally nsim would be 50, and outerfold and fold would be 10, but the
values below speed computation 200x compared to these recommended
values.  Note that here we are using the standardize=TRUE argument of
optL1 rather than the scaling=TRUE argument of opt.splitval.  These
two approaches to scaling are roughly equivalent, but the scaling
approaches are not the same (scaling=TRUE does z-score,
standardize=TRUE scales to unit central L2 norm), and results will not
be identical.  Also, using standardize=TRUE scales variables but
provides coeffients for the original scale, whereas using scaling=TRUE
scales variables in the training set then applies the same scales to
the test set.


\subsection{Summarization and plotting}

Cox fit on the continuous risk predictions:

@ 
<<coxfit>>=
coxfit.continuous <- coxph(surv.obj~preds)
summary(coxfit.continuous)
@ %def 


Dichotomize the cross-validated risk predictions at the median, for visualization:

@ 
<<dichot>>=
preds.dichot <- preds > median(preds)
@ %def 

Plot the ROC curve:

@ 
<<prelim,echo=FALSE>>=
options(device="pdf")
@ %def 

\begin{figure}[htbp]
\begin{center}
@ 
<<label=ROCplot,fig=TRUE,echo=TRUE,include=TRUE,eval=TRUE>>=
if(require(survivalROC)){
  nobs <- length(preds)
  cutoff <- 12
  preds.roc <- survivalROC(Stime=beer.survival$os, status=beer.survival$status,
                         marker=preds, predict.time=cutoff, span = 0.01*nobs^(-0.20))
  plot(preds.roc$FP, preds.roc$TP, type="l", xlim=c(0, 1), ylim=c(0, 1),   
       xlab=paste( "FP", "\n", "AUC = ", round(preds.roc$AUC, 3)), 
       lty=2,
       ylab="TP", main="LASSO predictions\n ROC curve at 12 months")
  abline(0, 1)
}
@ %def 
\end{center}
\caption{ ROC plot of cross-validated continuous risk predictions at 12 months. Note that the predictions are better if you don't randomly select 250 genes to start with!  We only did this to ease the load on the CRAN checking servers. }
\label{fig:fig1}
\end{figure}

\section{Getting coefficients for model fit on all the data}

Finally, we can get coefficients for the model fit on all the data,
for future use.  Note that nsim should ideally be greater than 1, to
train the model using multiple foldings for cross-validation.  The
output of opt1D or opt2D will be a matrix with one row per simulation.
The default behavior in opt.nested.crossval() is to take the
simulation with highest cross-validated partial log likelihood (CVL),
which is the recommended way to select a model from the multiple
simulations.

@ 
<<full.model>>=
beer.coefs <- opt1D(setpen="L1",nsim=1,                        
             response=surv.obj,                         
             penalized=dat.filt,
             fold=5,
             maxlambda1=5,
             positive=FALSE,
             standardize=TRUE,
             trace=FALSE)
@ %def 

We can also include unpenalized covariates, if desired.  Note that
when keeping only one variable for a penalized or unpenalized
covariate, indexing a dataframe like [1] instead of doing [, 1]
preserves the variable name.  With [, 1] the variable name gets
converted to "".

@ 
<<unpenalized.eg>>=
beer.coefs.unpen <- opt1D(setpen="L1", nsim=1,                        
             response=surv.obj,                         
             penalized=dat.filt[-1],   # This is equivalent to dat.filt[,-1]
             unpenalized=dat.filt[1],  
             fold=5,
             maxlambda1=5,
             positive=FALSE,
             standardize=TRUE,
             trace=FALSE)
@ %def 

Note the non-zero first coefficient this time, due to it being unpenalized:

@ 
<<lookatcoefs>>=
beer.coefs[1, 1:5]        #example output with no unpenalized covariates
beer.coefs.unpen[1, 1:5]  #example output with first covariate unpenalized
@ %def 

\section{Simulation of collinear high-dimensional data with survival or binary outcome}

The \emph{pensim} also provides a convenient means to simulation
high-dimensional expression data with (potentially censored) survival
outcome or binary outcome which is dependent on specified covariates.

\subsection{Binary outcome example}

First, generate the data.  Here we simulate 20 variables.  The first
15 (group ``a'') are uncorrelated, and have no association with
outcome.  The final five (group ``b'') have covariance of 0.8 to each
other variable in that group.  The response variable is associated
with the first variable group ``b'' (firstonly=TRUE) with a
coefficient of 2.

Binary outcomes for $n_s= 50$ samples are simulated as a Bernoulli
distribution with probability for patient s:

\begin{equation}
p_{s} =\frac{1}{1 + exp(-\beta X_{s})}
\end{equation}

with $\beta_{s,16} = 0.5$ and all other $\beta_{s,i}$ equal to zero.

The code for this simulation is as follows:

@ 
<<genbinary>>=
set.seed(9)
x <- create.data(nvars=c(15, 5),
                 cors=c(0, 0.8),
                 associations=c(0, 2),
                 firstonly=c(TRUE, TRUE),
                 nsamples=50,
                 response="binary",
                 logisticintercept=0.5)
@ %def 

Take a look at the simulated data:

@ 
<<lookbinary>>=
summary(x)
x$summary
@ %def 

A simple logistic model fails at variable selection in this case:

@ 
<<fitmodel>>=
simplemodel <- glm(outcome ~ ., data=x$data, family=binomial)
summary(simplemodel)
@ %def 

But LASSO does a better job, selecting several of the collinear
variables in the ``b'' group of variables which are associated with
outcome:

@ 
<<binarylassodemo>>=
lassofit <- opt1D(nsim=3, nprocessors=1, setpen="L1", penalized=x$data[1:20],
                 response=x$data[, "outcome"], trace=FALSE, fold=10)
print(lassofit)
@ %def 

And visualize the data as a heatmap:


\begin{figure}[htbp]
\begin{center}
@ 
<<label=heatmap, fig=TRUE, echo=TRUE, include=TRUE, eval=TRUE>>=
dat <- t(as.matrix(x$data[, -match("outcome", colnames(x$data))]))
heatmap(dat, ColSideColors=ifelse(x$data$outcome==0, "black", "white"))
@ %def 
\end{center}
\caption{Heatmap of simulated data with binary response.}
\label{fig:fig2}
\end{figure}

\subsection{Survival outcome example}

We simulate these data in the same way, but with
response=''timetoevent''.  Here censoring is uniform random between
times 2 and 10, generating approximately 34\% censoring:

@ 
<<survoutcome>>=
set.seed(1)
x <- create.data(nvars=c(15, 5),
                 cors=c(0, 0.8),
                 associations=c(0, 0.5),
                 firstonly=c(TRUE, TRUE),
                 nsamples=50,
                 censoring=c(2, 10),
                 response="timetoevent")
@ %def 

How many events are censored?

@ 
<<howmanycensored>>=
sum(x$data$cens==0)/nrow(x$data)
@ %def 

Kaplan-Meier plot of this simulated cohort:

\begin{figure}[htbp]
\begin{center}
@ 
<<label=simulatedKM, fig=TRUE, echo=TRUE, include=TRUE, eval=TRUE>>=
library(survival)
surv.obj <- Surv(x$data$time, x$data$cens)
plot(survfit(surv.obj~1), ylab="Survival probability", xlab="time")
@ %def 
\end{center}
\caption{Kaplan-Meier plot of survival of simulated cohort.}
\label{fig:fig3}
\end{figure}


\clearpage

\bibliographystyle{plain}
\bibliography{pensim}

@ 
<<>>=
sessionInfo()
@ %def 

\end{document}
