% \documentclass[serif]{beamer} % Serif for Computer Modern math font.
\documentclass[serif, handout]{beamer} % Handout mode to ignore pause statements
\hypersetup{colorlinks,linkcolor=,urlcolor=red}

\usefonttheme{serif} % Looks like Computer Modern for non-math text -- nice!
\setbeamertemplate{navigation symbols}{} % Supress navigation symbols
\usetheme{AnnArbor}  % CambridgeUS Blue and yellow, Shows current section title
% \usetheme{Berlin} % Displays sections on top
\usepackage[english]{babel}
% \definecolor{links}{HTML}{2A1B81}
% \definecolor{links}{red}

\setbeamertemplate{footline}[frame number] 

\mode<presentation>
% \mode<handout>{\setbeamercolor{background canvas}{bg=black!5}}

\title{Model Diagnostics\footnote{See last slide for copyright information.}}
\subtitle{STA312 Spring 2019}
\date{} % To suppress date

\begin{document}

\begin{frame}
  \titlepage
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Background Reading}
%\framesubtitle{} 

\begin{itemize}
    \item Chapter 7 in \emph{Applied Survival Analysis Using R} by Dirk Moore
    \item[]
    \item \emph{Modeling Survival Data: Extending the Cox Model} (2000) by Terry Thereau and Patricia Grambsch
\end{itemize}

\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{Overview}
\tableofcontents
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{What could go wrong?} \pause
%\framesubtitle{} 

\begin{itemize}
    \item Proportional hazards assumption could be incorrect. \pause The log-normal model is an example. \pause
    \item Relationships might not be straight-line.  \pause For example, \pause
\begin{displaymath}
    h(t) = h_0(t) \, \exp\{ \beta_1 \cos(\beta_2 x) \}
\end{displaymath} \pause
    \item Some individual observations may have too much influence on the results. \pause
    \item Look at residuals. \pause
    \item \emph{Martingale} residuals?
\end{itemize}
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Stochastic processes}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Stochastic Processes} \pause
%\framesubtitle{} 
\begin{itemize}
    \item A \emph{stochastic process} is an infinite collection of random variables. \pause
    \item A \emph{counting process} $N(t)$ counts the number of events up to and including time $t$.  \pause
    \item Let $N_i(t)$ be the number of deaths for patient $i$, in the interval $(0,t]$ \pause
    \item This means more general counts are possible (and useful). \pause
        \begin{itemize}
            \item Number of heart attacks.  \pause
            \item Number of major auto repairs.  \pause
            \item Number of admissions to hospital.  \pause
            \item Number of lectures missed. \pause
            \item Number of times a sexually transmitted disease was diagnosed (for one person).  \pause
        \end{itemize}
    \item These all are in the category of \emph{recurrent risks}. \pause
    \item Being at risk is also a stochastic process that can turn on or off.
\end{itemize}
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Stochastic processes formulation for survival analysis} \pause
%\framesubtitle{} 
The pair $(T_i,\delta_i)$ \pause is replaced by \pause

\begin{itemize}
    \item $N_i(t)$: Number of observed events in $(0,t]$ for unit $i$. \pause
    \item $Y_i(t) = \left\{ \begin{array}{ll} % ll means left left
   1     & \mbox{if unit $i$ is at risk at time $t$}  \\
   0     & \mbox{otherwise} 
   \end{array}  \right.$. \pause %
 \\  This is called the \emph{risk process}. \pause
\end{itemize}
And the probability distribution is determined by the hazard function \pause
\begin{displaymath}
    h_i(t) = h_0(t) e^{\mathbf{x}_i(t)^\top \boldsymbol{\beta}}
\end{displaymath} \pause
Note this is a conditional model, in which $\mathbf{x}_i$ is a fixed function of $t$. 
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Martingales} \pause
\framesubtitle{} 
A \emph{discrete-time martingale} \pause is a sequence of random variables $X_1, X_2, \ldots$
that satisfies \pause
\begin{itemize}
    \item $E(|X_n|)<\infty$ \pause
    \item $E(X_{n+1}|X_1, \ldots, X_n) = X_n$ \pause
\end{itemize}
Examples: 
\begin{itemize}
    \item An unbiased random walk. \pause
    \item A gambler's current fortune if the game is fair.
\end{itemize}
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{frame}
\frametitle{Martingale sequence with respect to another sequence} 
\framesubtitle{Still discrete time}  \pause
The sequence $Y_1, Y_2, \ldots$ is a martingale with respect to $X_1, X_2, \ldots$ if \pause
\begin{itemize}
    \item $E(|Y_n|)<\infty$ \pause
    \item $E(Y_{n+1}|X_1, \ldots, X_n) = Y_n$ \pause
\end{itemize} \vspace{4mm}
Example: Likelihood ratio. \pause
Let $L_n = \displaystyle \prod_{i=1}^n\frac{g(X_i)}{f(X_i)}$. \pause If $X_1, X_2, \ldots$ are independent with density $f(x)$\pause, then $\{L_1, L_2, \ldots \}$ is a martingale with respect to $\{X_1, X_2, \ldots \}$.
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{Continuous time martingale} \pause
%\framesubtitle{} 
A stochastic process $Y(t)$ is said to be a martingale with respect to the stochastic process $X(t)$ if for all $t$, \pause
\begin{itemize}
    \item $E(|Y(t)|)<\infty$ \pause
    \item $E(Y(t)|\{X(\tau): \tau\leq s\}) = Y(s)$ 
\end{itemize} \pause
Example: If $\widehat{S}(t)$ is the Kaplan-Meier estimate\pause, then under mild technical conditions, $\sqrt{D}(\widehat{S}(t)-S(t))$ is a continuous time martingale.
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{frame}
\frametitle{Martingale convergence theorems}
\framesubtitle{There are many versions} \pause
Let $X_n$ be a martingale satisfying $\sup_{t>0}E(|X|^p<\infty)$ for some $p>1$. \pause Then there exists a random variable $X$ such that \pause
%\framesubtitle{} 
{\LARGE
\begin{displaymath}
    P(\lim_{n \rightarrow \infty} X_n = X) = 1
\end{displaymath} 
} % End size
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{frame}
\frametitle{Martingale Central Limit Theorems}
\framesubtitle{Again there are quite a few versions}  \pause
Under some technical conditions, sums of (normalized) independent martingales converge to a Brownian motion process $B(t)$, for which \pause

\begin{itemize}
%    \item $B(t)$ has continuous sample paths. \pause
    \item $B(0)=0$. \pause
    \item $E(B(t))=0$ for all $t$. \pause
    \item Independent increments: \pause $B(t)-B(u)$ is independent of $B(u)$ for any $0 \leq u \leq t$. \pause
    \item Gaussian process: \pause For any positive integer $n$ and time points $t_1, \ldots, t_n$, the joint distribution of $B(t_1), \ldots, B(t_n)$ is multivariate normal. 
\end{itemize}
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{Doob-Meyer decomposition Theorem} \pause
%\framesubtitle{} 
Any counting process $N_i(t)$ can be decomposed into  \pause
{\LARGE
\begin{displaymath}
    N(t) = \Lambda(t) + M(t), \pause
\end{displaymath} 
} % End size

\noindent
where $M(t)$ is a martingale and $\Lambda(t)$ is a ``predictable" stochastic process. \pause
\vspace{4mm}

``Predictable" has an intense mathematical definition, \pause but the idea is that the distribution of $\Lambda_{n+1}(t)$ depends on the distribution of $\Lambda_1(t), \ldots, \Lambda_n(t)$. 
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{Decomposition for the Proportional Hazards Model}
\framesubtitle{Special case of survival (one event) and right censored data}  \pause
Let $N_i(t)=1$ if unit $i$ failed in $(0,t]$, and zero otherwise. \pause
{\LARGE
\begin{displaymath}
    N_i(t)  = H_i(t) + M_i(t),
\end{displaymath} \pause
} % End size
where $H_i(t) = \int_0^y h_i(s) \, ds$ is the cumulative hazard. 
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Residuals}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{Martingale Residuals}
\framesubtitle{Based on $N_i(t)  = H_i(t) + M_i(t)$}  \pause
\begin{displaymath}
    \widehat{M}_i(t) = N_i(t) - \widehat{H}_i(t)
\end{displaymath} \pause

Evaluated at $t_i$, the \emph{estimated} martingale residual is \pause
\begin{eqnarray*}
    \widehat{M}_i(t_i) & = & \delta_i - \widehat{H}_i(t) \\ \pause
    & = &  \delta_i + e^{\mathbf{x}_i(t)^\top \widehat{\boldsymbol{\beta}}}\log\left(\widehat{S}_0(t_i) \right)
\end{eqnarray*} \pause

\begin{itemize}
    \item Martingale residuals are martingales. \pause
    \item Add to zero. \pause
    \item Large values need investigation. \pause
    \item Plots against $x$ variables can reveal the functional form of the dependence of survival time on $x$.  % \pause
%    \item Starting with a model that has no explanatory variables, plot residuals against variables not in the model. \pause Add a smoothed curve.
\end{itemize}
\end{frame}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{Schoenfeld residuals} \pause
%\framesubtitle{} 
We have already seen
\begin{displaymath}
    \sum_{i=1}^D \left( x_{(i)} - \sum_{j \in R_i}x_j \frac{e^{\widehat{\beta}x_j} }
                                {\sum_{k \in R_j}e^{\widehat{\beta}x_k}} \right) = 0
\end{displaymath} \pause

\begin{itemize}
    \item The terms that add to zero are called the Schoenfeld residuals \pause
    \item There is one set for each explanatory variable. \pause
    \item Unusually large or small values are worthy of investigatoin. \pause
    \item They can be approximately standardized, which helps. \pause
    \item They can be used to form a chi-squared test of $H_0:$ Proportional hazards. \pause (Thereau and Grambsch, Chapter 6).
\end{itemize}
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{Case Deletion Residuals} \pause
%\framesubtitle{} 
\begin{itemize}
    \item Let $\widehat{\boldsymbol{\beta}}_{(i)}$ denote the partial MLE of $\boldsymbol{\beta}$ with case $i$ deleted. \pause
    \item Calculate $\widehat{\boldsymbol{\beta}}_{(i)} - \widehat{\boldsymbol{\beta}}$. \pause
    \item There will be $p$ differences. \pause
    \item These are called \texttt{dfbeta}. \pause
    \item They can be standardized. \pause
    \item The standardized versions are called \texttt{dfbetas}. \pause
    \item They can reveal observations that are overly influential.
\end{itemize}
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{Copyright Information}

This slide show was prepared by  \href{http://www.utstat.toronto.edu/~brunner}{Jerry Brunner},
Department of Statistics, University of Toronto. It is licensed under a 
\href{http://creativecommons.org/licenses/by-sa/3.0/deed.en_US}
     {Creative Commons Attribution - ShareAlike 3.0 Unported License}. Use any part of it as you like and share the result freely. The \LaTeX~source code is available from the course website:
\href{http://www.utstat.toronto.edu/~brunner/oldclass/312s19} {\footnotesize \texttt{http://www.utstat.toronto.edu/$^\sim$brunner/oldclass/312s19}}

\end{frame}

\end{document}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Removed this for simplicity
\begin{frame}
\frametitle{Decomposition for the Proportional Hazards Model}
%\framesubtitle{} 
Let $N_i(t)$ represent the number of events in $(0,t]$ for unit $i$, and let $Y_i(t)$ be the risk process. \pause
{\Large
\begin{eqnarray*}
    N_i(t)  & = & \Lambda_i(t) + M_i(t) \\ \pause
            & = & \int_0^t Y_i(s) h_i(s) \, ds + M_i(t) \\ \pause
            & = & \int_0^t Y_i(s) h_0(s)e^{\mathbf{x}_i(t)^\top \boldsymbol{\beta}} \, ds + M_i(t) \\ \pause
            & = & e^{\mathbf{x}_i(t)^\top \boldsymbol{\beta}}\int_0^t Y_i(s) h_0(s) \, ds + M_i(t) 
\end{eqnarray*}
} % End size
\end{frame}


P. Grambsch and T. Therneau (1994), Proportional hazards tests and diagnostics based on weighted residuals. Biometrika, 81, 515-26.

Therneau, T. M., P. M. Grambsch, and T. R. Fleming. 1990. Martingale based residuals for survival models. Biometrika 77:147–60.