ROCAD analysis

make_package.R

# Time-stamp: <2010-04-05 17:11:20 chl>
#
# Automaticaly collate R file and try to build a package with
# documentation generated by Roxygen.
#
require(roxygen)
rm(list=ls())
package.skeleton('rocad',code_files=c('rocad.R','rocad-package.R'),force=T)
roxygenize('rocad',roxygen.dir='rocad',copy.package=F,unlink.target=F)
try(system("R CMD build rocad"))
try(system("R CMD Rd2dvi --pdf rocad"))

rocad.nw

% Time-stamp: <2010-04-05 17:13:26 chl>

%
% Small illustrations around ROCAD analysis, using rocad.R utilities
% functions.
%

% To produce a pdf, use something like:
%   $ noweave -index -delay rocad.nw > rocad.tex
%   $ pdflatex rocad.tex (x 2)
%
% Html output may be produced using
%   $ noweave -filter l2h -index -html rocad.nw | htmltoc > rocad.html
%
% To extract relevant code, use e.g.
%   $ notangle -L -Rrocad.R rocad.nw > rocad.R
\documentclass[10pt]{article}
\usepackage{noweb}
\noweboptions{smallcode,longchunks}
\usepackage{subfigure}
\usepackage{graphicx}
\usepackage[usenames,dvipsnames]{color}
\usepackage{wrapfig}
\usepackage{hyperref}
\hypersetup{colorlinks=true,citecolor=red,urlcolor=NavyBlue}
\makeatletter
\renewcommand{\@cite}[1]{#1} % remove bracket around citations
\makeatother
\newcommand{\marcite}[1]{\marginpar{\raggedright\small\cite{#1}}}

\setlength{\parindent}{0pt}
\setlength{\parskip}{1ex plus 0.5ex minus 0.2ex}

\title{ROCAD analysis}
\author{Christophe Lalanne\thanks{Email: \href{mailto:ch.lalanne@gmail.com}{ch.lalanne@gmail.com}}}

\begin{document}
\maketitle
\pagestyle{noweb}
\tableofcontents

\section{What are ROC area discimination curve}
\label{sec:intro}

\section{How it may be implemented}
\label{sec:code}
@
We first need to compute the area under the curve, which is readily implemented as shown below:
<<rocad.R>>=
auc <- function(se, sp) {
	trap.rule <- function(x,y) { sum(diff(x)*(y[-1]+y[-length(y)]))/2 }
	return(trap.rule(1-sp,se))	
}
@ %def auc

@
We then define a funciton that iteratively compute Sensibility and Specificity given item score and a binary gold standard.
<<rocad.R>>=
tr.iterate <- function(item, score) {
	tpr <- function(tab) { tab["1","1"]/apply(tab,2,sum)["1"] }
	tnr <- function(tab) { tab["0","0"]/apply(tab,2,sum)["0"] }
	if (length(unique(item))>1 & !is.factor(item)) item <- as.factor(item)
	se <- sp <- numeric(nlevels(item)-1)
	for (k in 1:(nlevels(item)-1)) {
		tmp <- factor(ifelse(item==k,0,1), levels=c("0","1"))
		tab <- table(tmp,score)
		se[k] <- tpr(tab)
		sp[k] <- tnr(tab)
	}
	auc <- auc(se,sp)
	return(list(se=se,sp=sp,auc=auc))
}
@ %def tr.iterate

@
It can be tested using a simple simulation:
<<rocad-testing.R>>=
res <- tr.iterate(test[,1],tot.score.bin)
plot(1-res$sp,res$se,type="b")
@

%\bibliographystyle{alpha}
%\begin{small}
%\bibliography{new}
%\end{small}


\clearpage
\subsection*{Defined chunks}
\par\noindent
\nowebchunks
@

\subsection*{Index}
\par\noindent
\nowebindex
@
\end{document}

rocad.R

# Time-stamp: <2010-04-05 11:08:15 chl>

#          +--------+-------+
#          |  + (1) | - (0) |
#  +-------+--------+-------+
#  | + (1) |   TP   |  FP   |
#  +-------+--------+-------+
#  | - (0) |   FN   |  TN   |
#  +-------+--------+-------+
#

#' A simulated data set.
#'
#' Simulate a set of categorical responses.
#' 
#' @param None
#' @return A matrix with items responses
#' @seealso \code{sim.item} in the \code{psych} package.
#' @author chl
test.data <- function() {
	require(psych)
	test <- sim.item(nvar=20, nsub=300, categorical=TRUE, low=0, high=5)
	return(test)
}

#' Calculate Area Under the Curve (AUC).
#'
#' Given a value for the specificity and sensibility of a test, 
#' calculate the AUC using the trapezoidal method.
#'
#' @param se sensibility
#' @param sp specificity
#' @return The value of AUC
#' @note Should add the 95% CIs; should also look at colAUC {caTools}
#'       which provide column-wise computation
#' @references Centor RM, Schwartz JS. An evaluation of methods for
#'             estimating the area under the receiver operating
#'             characteristic (ROC) curve. 
#'             Med Decis Making (1985) 5:149-56.
#' @author chl
auc <- function(se, sp) {
	trap.rule <- function(x,y) { sum(diff(x)*(y[-1]+y[-length(y)]))/2 }
	return(trap.rule(1-sp,se))	
}

#' Calculate Se and Sp values for an ordinal variable.
#'
#' Assuming an item with \eqn{k} categories, and a binary variable (gold 
#' standard), compute Sensibility and Specificity values.
#'
#' @param item a vector of responses
#' @param score a binary vector of classification
#' @return A list with se and sp computed for the \eqn{k-1} categories
#' @seealso other func in package
#' @author chl
tr.iterate <- function(item, score) {
	tpr <- function(tab) { tab["1","1"]/apply(tab,2,sum)["1"] }
	tnr <- function(tab) { tab["0","0"]/apply(tab,2,sum)["0"] }
	if (length(unique(item))>1 & !is.factor(item)) item <- as.factor(item)
	se <- sp <- numeric(nlevels(item)-1)
	for (k in 1:(nlevels(item)-1)) {
		tmp <- factor(ifelse(item==k,0,1), levels=c("0","1"))
		tab <- table(tmp,score)
		se[k] <- tpr(tab)
		sp[k] <- tnr(tab)
	}
	auc <- auc(se,sp)
	return(list(se=se,sp=sp,auc=auc))
}

Archiving (very) old stuff.