\documentclass[11pt]{article}

\usepackage{amsmath,amsfonts,graphicx} %Typical math resource packages
\usepackage{framed,color}
\definecolor{shadecolor}{rgb}{0.93, 0.93, 0.93}

%\usepackage{setspace}
\newtheorem{theorem}{Theorem}[section]

\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{corollary}[theorem]{Corollary}
%\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{lemma}{Lemma}[section]

\setlength{\hoffset} {-0pt}
\setlength{\textwidth} {6.0 in}
\setlength{\textheight} {8.6in}
\setlength{\oddsidemargin}{0.2 in}
\setlength{\evensidemargin}{0.2 in}
\setlength{\topmargin} {-0.3in}

\def\widenspacing{\multiply\baselineskip by 130
    \divide\baselineskip by 100}

\begin{document}

\begin{center}
{\Large \bf
Some Topics in Statistical Computation Using R \\
}

\medskip

\medskip

{\sc
Sta 695 002 Spring 2005\\
%\footnote{These are based on a class Sta 695. The actual materials covered may be different.} \\
}
\vspace{0.2in}
{\large
Department of Statistics \\
University of Kentucky \\
Lexington, KY 40536 \\
}

\end{center}


\medskip
\begin{center}
{\bf \small SUMMARY}
\end{center}
\begin{list}{}{\leftmargin 0.5in \rightmargin 0.5in} \item
{\small
These notes are based on a class of Sta 695, given by Dr. Mai Zhou in Spring 2005. 
The actual materials covered in class may be different.
}
\end{list}

{\sl AMS 2000 Subject Classification:}  Primary 60E15; secondary 60G30.
%\vspace{1ex} \newline

{\sl Key Words and Phrases: Truncated data, NPMLE, square error loss.}

\newpage

\widenspacing
%\vspace{-1ex}

\section{Introduction}

We aim to cover some topics that are {\em interesting} and {\em useful} and related to R.

We avoid topics that are too theoretical.
We try to make the level of the class
in between Sta 605 and Sta 705.
Finally, two semesters of Statistical Inference (based on Casella and Berger, for example)
is assumed.





\section{Errors in the R Numerical Computation}





\section{Examples}
\indent

We have developed computational codes for the
calculation of the Bayes estimator. These codes are
packaged as an R software package.
Available for downloading at
\begin{verbatim}
http://www.ms.uky.edu/~mai/NPBayesCensor.html
\end{verbatim}
First you need to install R on your computer, then install the package NPBayes.
The examples below are performed with these codes.

\textbf{Example 1} The following example contains one uncensored
observation and two right censored observations. Suppose the three
observations with the format $(X_i, Y_i, \delta_i)$ are $(9, 0.2,
1)$, $(13, 4, 0)$, $(15, 10, 0)$. Then in Table 1, we order all time points
(lifetimes and truncation times).

\begin{table}[h!tb] \centering
\caption{\label{1} {Ordered Data Points for Example 1}}
\begin{tabular}{lcccccc} \hline
\textbf{Time Points} & $y_1$ & $y_2$ & $x_1$ & $y_3$ & $x_2$ & $x_3$\\
\hline
\textbf{Time} & 0.2 & 4 & 9 & 10 & 13 & 15\\
 \hline
 \textbf{$p_i$} & -1 & -1 & 1 & -1 & 1 & 1 \\
 \hline\\
\end{tabular}
\end{table}

For $u=14$, after simplification based on equation
$(\ref{BayesES})$, we obtain
\begin{eqnarray}
\hat{S}_{D}(14)=\frac{\frac{(\alpha[14,\infty)+1)(\alpha[13,\infty)+2)}{(\alpha[10,\infty)+2)(\alpha[4,\infty)+2)(\alpha[0.2,\infty)+1)}}
{\frac{(\alpha[13,\infty)+1)\alpha[0,\infty)}{(\alpha[10,\infty)+1)(\alpha[4,\infty)+1)(\alpha[0.2,\infty)+0)}}~.
\nonumber
\end{eqnarray}

Let
\begin{eqnarray}
\alpha[\omega,\infty)=B\exp{(-\theta{\omega})}.
\end{eqnarray}
Table \ref{3data} illustrates the values of $\hat{S}_{D}(14)$
for $\theta=0.12$ and different values of $B$.

\begin{table}[h!tb] \centering
\caption{\label{3data} {Nonparametric Bayes Estimator of
$\hat{S}_{D}(14)$ for Example 1} }
\begin{tabular}{ccc} \hline
{B} & {$\theta$} & {$\hat{S}_{D}(14)$}\\
\hline
8 & 0.12 & 0.2509\\
1 & 0.12 & 0.3741\\
0.1 & 0.12 & 0.4687\\
0.001 & 0.12 & 0.4879\\\hline
\end{tabular}
\end{table}

Moreover, the product limit estimator at $u=14$ is
\begin{eqnarray}
\hat{S}_{PL}(14)=\prod_{x_i\leq{14}}\left(1-\frac{D_i}{R_i}\right)=(1-\frac{1}{2})=0.5
\nonumber
\end{eqnarray}
We can observe that when $\alpha(.)$ approaches $0$, the Bayes
estimator is very close to the product limit estimator.

Figure \ref{Qqplot400} presents the nonparametric Bayes estimator of
survival function for $B=8$ and $\theta=0.12$.


\begin{figure}[ht]
\begin{center}
\rotatebox{360}{ \resizebox{9cm}{9cm}
 {\includegraphics{NPBayes1.pdf}}}
\caption{Nonparametric Bayes Estimator of Survival Function for
Example 1} \label{Qqplot400}
\end{center}
\end{figure}

\textbf{Example 2} Presented in Table \ref{3} is an example with two
uncensored observations and four right censored observations.
Suppose the six observations with the format of $(X_i, Y_i,
\delta_i)$ are $(0.6, 0.1, 1)$, $(1.5, 0.3, 1)$, $(2.9, 0.5, 0)$,
$(3.1, 0.9, 0)$, $(3.7, 3.2, 0)$, and $(4.3, 4.2, 0)$. Then we
order all time points (lifetimes and truncation times).

\begin{table}[htb] \centering
\caption{\label{3} {Ordered Data Points for Example 2}}
\begin{tabular}{lcccccccccccc} \hline
\textbf{Time Points} & $y_1$ & $y_2$ & $y_3$ & $x_1$ & $y_4$ & $x_2$ &
                       $x_3$ & $y_5$ & $x_4$ & $x_5$ & $y_6$ & $x_6$\\
\hline
\textbf{Time} & 0.1 & 0.3 & 0.5 & 0.6 & 0.9 & 1.5 & 2.9 & 3.1 & 3.2 & 3.7 & 4.2 & 4.3\\
 \hline
 \textbf{$p_i$} & -1 & -1 & -1 & 1 & -1 & 1 & 1 & -1 & 1 & 1 & -1 & 1 \\
 \hline\\
 %\caption{Data for Example 1}
\end{tabular}
\end{table}

For $u=3.9$, after simplification based on equation
$(\ref{BayesES})$, we have
\begin{eqnarray}
\hat{S}_{D}(3.9)=\frac{\frac{(\alpha[3.9, \infty)+1)(\alpha[3.7,
\infty)+2)(\alpha[3.2, \infty)+3)(\alpha[2.9, \infty)+3)}
{(\alpha[3.1, \infty)+2)(\alpha[0.9, \infty)+3)(\alpha[0.5,
\infty)+3)(\alpha[0.3, \infty)+2)(\alpha[0.1, \infty)+1)}}
{\frac{(\alpha[3.7, \infty)+1)(\alpha[3.2, \infty)+2)(\alpha[2.9,
\infty)+2)(\alpha[0, \infty))} {(\alpha[3.1,
\infty)+1)(\alpha[0.9, \infty)+2)(\alpha[0.5,
\infty)+2)(\alpha[0.3, \infty)+1)(\alpha[0.1, \infty)+0)}}.
\nonumber
\end{eqnarray}

Again, let
\begin{eqnarray}
\alpha[\omega,\infty)=B\exp{(-\theta{\omega})}.
\end{eqnarray}

The $\hat{S}_{D}(3.9)$ for $\theta=0.12$ and different values of
$B$ are presented in Table \ref{table4}.

\begin{table}[htb] \centering
\caption{\label{table4} {Nonparametric Bayes Estimator of
$\hat{S}_{D}(3.9)$ for Example 2 }}
\begin{tabular}{ccc} \hline
{B} & {$\theta$} & {$\hat{S}_{D}(3.9)$}\\
\hline
8 & 0.12 & 0.5585\\
1 & 0.12 & 0.4612\\
0.1 & 0.12 & 0.4277\\\hline
\end{tabular}
\end{table}

Moreover, the product limit estimator is
\begin{eqnarray}
\hat{S}_{PL}(3.9)=\prod_{x_i\leq{3.9}}\left(1-\frac{D_i}{R_i}\right)=(1-\frac{1}{3})(1-\frac{1}{3})=0.44.
\nonumber
\end{eqnarray}
We can see that when $\alpha(.)$ approaches $0$, the Bayes
estimator is very close to product limit estimator.

The nonparametric Bayes estimator of survival function is plotted
in Figure \ref{Qqplot401} for $B=8$ and $\theta=0.12$.


{\bf Example 3} Finally we shall use a data set from the book of Klein and Moeschberger (1997).
We take all the male residents from the Channing housing data. There are 97 such cases.
There is one case that the subject
enters and dies at the same month. Since our program does not allow truncation and survival times
to be equal, we just add 0.1 to all survival times to break tie. Below are the codes and result:
$\hat S(1050) = 0.2802438$.

\begin{snugshade}
\begin{verbatim}
library(NPBayes)
library(KMsurv)
data(channing)

cha1 <- channing[,(channing$gender==1)]
NPBayesT(B=10, theta=0.0002, u=1050, y=cha1[,3],
                   x=cha1[,4]+0.1, status=cha1[,2])
## [1] 0.2802438
\end{verbatim}
\end{snugshade}

The measure that is the parameter of the Dirichlet process prior is taken to be
the same exponential measure as in previous examples. But this can easily be changed by
rewrite the {\tt MeasureB( )} function inside the {\tt NPBayes} package.

{\bf Remark 3}: The same package {\tt NPBayes}
also contains function that can handle interval censored data.
However, the computation time increases exponentially with
the number of truly interval censored data points.
On the other hand, the speed of the computation is quite fast when dealing with
right censored and left truncated data.

\addcontentsline{toc}{chapter}{\protect\numberline{ }{References}}
\begin{thebibliography}{plain}

\bibitem{KleinMoeschberger}
Klein, J. and Moeschberger, M. (1997),
{\em Survival Analysis}
Springer, New York

\bibitem{R}
R Development Core Team (2011). 
{\em R: A language and environment for statistical computing}
R Foundation for Statistical Computing,
Vienna, Austria. 
ISBN 3-900051-07-0, URL http://www.R-project.org/.

\bibitem{Zhou2}
Zhou, M. (2004). Nonparametric Bayes estimator of survival
functions for doubly/interval censored data. \textit{Statistica
Sinica} \textbf{14}, 533-546.
\end{thebibliography}

\end{document}
