%\VignetteIndexEntry{Recurrent Events}
%\VignetteEngine{R.rsp::tex}
%\VignetteKeyword{R}
%\VignetteKeyword{package}
%\VignetteKeyword{vignette}
%\VignetteKeyword{LaTeX}
\PassOptionsToPackage{usenames}{xcolor}
\PassOptionsToPackage{hidelinks,colorlinks,linkcolor={blue!50!black},urlcolor={blue!50!black},citecolor={blue!50!black}}{hyperref}
\documentclass[a4paper]{tufte-handout}
\usepackage{etex}
\usepackage{etoolbox}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{mathtools}
\usepackage[strict=true]{csquotes}
\usepackage{xcolor}
\usepackage{natbib}
\usepackage{amsmath,amssymb,textcomp}
\usepackage{bm}
\usepackage{graphicx}
\usepackage{wrapfig}
\usepackage{rotating}
\usepackage{capt-of}
\usepackage{listings}
\usepackage{booktabs}
\usepackage{multicol}
\usepackage{longtable}
\usepackage{zlmtt}
\usepackage{float}
\usepackage{fancyvrb}
\usepackage{verbatim}
\usepackage[english]{babel}
\usepackage{hyperref}
\usepackage{environ}
\NewEnviron{mnote}{\marginnote{\BODY}}
\NewEnviron{snote}{\sidenote{\BODY}}
\usepackage{xspace}
\usepackage{url}
\usepackage{amsthm}
\newtheorem{thm}{Theorem.}
\newtheorem{prop}{Proposition.}
\theoremstyle{definition}
\newtheorem{defn}[thm]{Definition.}
\newtheorem{exa}[thm]{Example}
\newcommand{\R}{\ensuremath{\mathbb{R}}} 
\newcommand{\Real}{\ensuremath{\mathbb{R}}} 
\newcommand{\Complex}{\ensuremath{\mathbb{C}}} 
\newcommand{\Z}{\ensuremath{\mathbb{Z}}} 
\newcommand{\N}{\ensuremath{\mathbb{N}}} 
\newcommand{\norm}[1]{\ensuremath{\left\Vert#1\right\Vert}} 
\newcommand{\var}{\ensuremath{\mathbb{V}\text{ar}}} 
\newcommand{\cov}{\ensuremath{\mathbb{C}\text{ov}}} 
\newcommand{\cor}{\ensuremath{\mathbb{C}\text{or}}} 
\newcommand{\E}{\ensuremath{\mathbb{E}}} 
\newcommand{\pr}{\ensuremath{\mathbb{P}}} 
\newcommand{\from}{\leftarrow} 
\newcommand{\To}[2]{\ensuremath{#1\!\to\!#2}}
\newcommand{\From}[2]{\ensuremath{#1\!\from\!#2}}
\newcommand{\chain}[3]{\ensuremath{#1\!\to\!#2\to\!#3}}
\newcommand{\ichain}[3]{\ensuremath{#1\!\from\!#2\from\!#3}}
\newcommand{\fork}[3]{\ensuremath{#1\!\from\!#2\to\!#3}}
\newcommand{\ifork}[3]{\ensuremath{#1\!\to\!#2\from\!#3}}
\newcommand{\pa}{\text{pa}}
\newcommand{\abs}[1]{\ensuremath{\left\vert#1\right\vert}} 
\newcommand{\ipr}[1]{\langle#1\rangle} 
\newcommand{\set}[1]{\left{#1\right}} 
\newcommand{\seq}[1]{\left<#1\right>} 
\renewcommand{\subset}{\subseteq}  
\renewcommand{\supset}{\supseteq} 
\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\argmin}{arg\,min}
\DeclareMathOperator{\sgn}{sgn}
\DeclareMathOperator{\mvec}{vec}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator{\im}{im}
\DeclareMathOperator{\diag}{diag}
\DeclareMathOperator{\bias}{bias}
\DeclareMathOperator{\rank}{rank}
\DeclareMathOperator{\logit}{logit}
\DeclareMathOperator{\expit}{expit}
\makeatletter
 \def\mathcolor#1#{\@mathcolor{#1}}
 \def\@mathcolor#1#2#3{%
   \protect\leavevmode
   \begingroup
     \color#1{#2}#3%
  \endgroup
 }
\newcommand{\Dto}{\overset{\mathcal{D}}{\longrightarrow}}
\newcommand{\Pto}{\overset{P}{\longrightarrow}}
\newcommand{\Wto}{\overset{W}{\longrightarrow}}
\newcommand{\VV}{\bm{\Omega}_\theta}
\newcommand{\independenT}[2]{\mathrel{\setbox0\hbox{$#1#2$}\copy0\kern-\wd0\mkern4mu\box0}} 
\newcommand{\indep}{\protect\mathpalette{\protect\independenT}{\perp}}
\DefineVerbatimEnvironment{verbatim}{Verbatim}{xleftmargin=0.5em,xrightmargin=0em,numbers=none,frame=none,fontsize=\footnotesize,formatcom={\color[rgb]{0.2,0.2,0.2}}}
\setlength{\parindent}{0pt} % Kills annoying indents. 
\let\iint\relax 
\let\iiint\relax

\lstset{basicstyle=\ttfamily\small,keywordstyle=\color{black},commentstyle=\color{gray}\ttfamily\itshape,stringstyle=\color[rgb]{0,0,0.5},columns=fullflexible,alsoletter=.,texcl=true,escapeinside={*@}{@*)},escapebegin=\lst@commentstyle\,,breaklines=true,breakatwhitespace=false,numbers=left,numberstyle=\ttfamily\tiny\color{gray},stepnumber=1,numbersep=10pt,backgroundcolor=\color{white},tabsize=4,showspaces=false,showstringspaces=false,xleftmargin=.23in,frame=lines
,rulesepcolor=\color[rgb]{0.85,0.85,0.85},basewidth={0.5em,0.42em},language=r,label=
,caption= ,captionpos=b}
\lstset{basicstyle=\ttfamily\small,keywordstyle=\color{black},commentstyle=\color{gray}\ttfamily\itshape,stringstyle=\color[rgb]{0,0,0.5},columns=fullflexible,alsoletter=.,texcl=true,escapeinside={*@}{@*)},escapebegin=\lst@commentstyle\,,breaklines=true,breakatwhitespace=false,numbers=left,numberstyle=\ttfamily\tiny\color{gray},stepnumber=1,numbersep=10pt,backgroundcolor=\color{white},tabsize=4,showspaces=false,showstringspaces=false,xleftmargin=.23in,frame=lines
,rulesepcolor=\color[rgb]{0.85,0.85,0.85},basewidth={0.5em,0.42em},language=python,label=
,caption= ,captionpos=b}

\setlength{\parindent}{0pt} % Kills annoying indents.
\newcommand{\n}{}
\author{Klaus Holst \& Thomas Scheike}
\date{\today}
\title{Recurrent Events}
\hypersetup{
 pdfauthor={Klaus Holst \& Thomas Scheike},
 pdftitle={Recurrent Events},
 pdfkeywords={},
 pdfsubject={},
 pdfcreator={Emacs 26.1 (Org mode 9.1.14)}, 
 pdflang={English}}
\begin{document}

\maketitle
\noindent\rule{\textwidth}{0.5pt}

\section*{Overview}
\label{sec:org01ebcd2}

For recurrent events data it is often of interest to compute basis descriptive
quantities as a first go at getting some basic understanding of the phenonmenon
studied. We here demonstrate how one can compute 

\begin{itemize}
\item the marginal mean
\item the variance
\item the probability of exceeding k events
\end{itemize}

In addition several tools can be used for simulating recurrent events and
bivariate recurrent events data, in the case with a possible terminating event. 

We start by simulating some recurrent events data with two type of events with
cumulative hazards 

\begin{itemize}
\item \(\Lambda_1(t)\)
\item \(\Lambda_2(t)\)
\item \(\Lambda_D(t)\)
\end{itemize}

where we consider types 1 and 4 and with a rate of the terminal event given by 
\(\Lambda_D(t)\). We let the events be independent, but could also specify a random effects
structure to generate dependence. 

When simulating data we can impose various random-effects structures to generate dependence  
\begin{itemize}
\item We can draw normally distributed random effects \(Z_1,Z_2,Z_d\) were the variance (var.z) 
and correlation can be specified (cor.mat) (dependence=2).  
Then the intensities are 
\begin{itemize}
\item \(\exp(Z_1) \lambda_1(t)\)
\item \(\exp(Z_2) \lambda_2(t)\)
\item \(\exp(Z_3) \lambda_D(t)\)
\end{itemize}
\item We can one gamma distributed random effects \(Z\).  Then the intensities are (dependence=1)
\begin{itemize}
\item \(Z \lambda_1(t)\)
\item \(Z \lambda_2(t)\)
\item \(Z \lambda_D(t)\)
\end{itemize}
\item We can draw gamma distributed random effects \(Z_1,Z_2,Z_d\) were the sum-structure can be speicifed via a matrix
cor.mat. Then we compute \(\tilde Z_j = \sum_k  Z_k^{cor.mat(j,k)}\) for \(j=1,2,3\)  (dependence=3)
Then the intensities are
\begin{itemize}
\item \(\tilde Z_1 \lambda_1(t)\)
\item \(\tilde Z_2 \lambda_2(t)\)
\item \(\tilde Z_3 \lambda_D(t)\)
\end{itemize}
\item The intensities can be independent (dependence=0)
\end{itemize}

We return to how to run the different set-ups later and start by simulating independent processes. 


\subsection*{Utility functions}
\label{sec:org9be05a8}

We here mention two utility functions 

\begin{itemize}
\item tie.breaker for breaking ties among jump-times which is expected in the functions below.
\item count.history that counts the number of jumps previous for each subject that is \(N_1(t-)\) and \(N_2(t-)\).
\end{itemize}


\subsection*{Marginal Mean}
\label{sec:orga8b27b9}

We start by estimating the marginal mean \(E(N_1(t \wedge D))\) where \(D\) is the timing of the terminal event. 

This is based on a  rate model  for 

\begin{itemize}
\item the type 1 events
\item the terminal event
\end{itemize}

and is defined as \(\mu_1(t)=E(N_1^*(t))\) 
\begin{align}
   \int_0^t S(u) d R_1(u)	
\end{align}
where \(S(t)=P(D \geq t)\) and \(dR_1(t) = E(dN_1^*(t) | D \geq t)\) 

and can therefore be estimated by a 

\begin{itemize}
\item Kaplan-Meier estimator, \(\hat S(u)\)
\item Nelson-Aalen estimator for \(R_1(t)\)
\end{itemize}

\begin{align}
  \hat R_1(t) & =   \sum_i \int_0^t  \frac{1}{Y_\bullet (s)}  dN_{1i}(s)
\end{align}
where \(Y_{\bullet}(t)= \sum_i Y_i(t)\) such that the estimator is 
\begin{align}
  \hat \mu_1(t) & =    \int_0^t \hat S(u) d\hat R_1(u).
\end{align}

Cook \& Lawless (1997), and developed further in Gosh \& Lin (2000). 

The variance can be estimated based on the asymptotic expansion 
of \(\hat \mu_1(t) - \mu_1(t)\)
\begin{align*}
  & \sum_i \int_0^t \frac{S(s)}{\pi(s)} dM_{i1}  - \mu_1(t) \int_0^t  \frac{1}{\pi(s)} dM_i^d +  \int_0^t \frac{\mu_1(s) }{\pi(s)} dM_i^d,
\end{align*}

with mean-zero processes 

\begin{itemize}
\item \(M_i^d(t) = N_i^D(t)- \int_0^t Y_i(s) d \Lambda^D(s)\),
\item \(M_{i1}(t) = N_{i1}(t) - \int_0^t Y_{i}(s) dR_1(s)\).
\end{itemize}

as in Gosh \& Lin (2000)

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
library(mets)
set.seed(1000) # to control output in simulatins for p-values below.

data(base1cumhaz)
data(base4cumhaz)
data(drcumhaz)
ddr <- drcumhaz
base1 <- base1cumhaz
base4 <- base4cumhaz
rr <- simRecurrent(1000,base1,death.cumhaz=ddr)
rr$x <- rnorm(nrow(rr)) 
rr$strata <- floor((rr$id-0.01)/500)
dlist(rr,.~id| id %in% c(1,7,9))
\end{lstlisting}

\begin{verbatim}
id: 1
  entry time  status rr dtime fdeath death start stop  x     strata
1 0     133.1 0      1  133.1 1      1     0     133.1 1.185 0     
------------------------------------------------------------ 
id: 7
     entry  time   status rr dtime fdeath death start  stop   x       strata
7       0.0  813.3 1      1  1729  1      0        0.0  813.3  1.5495 0     
1004  813.3 1288.4 1      1  1729  1      0      813.3 1288.4  1.0535 0     
1658 1288.4 1315.4 1      1  1729  1      0     1288.4 1315.4  1.5330 0     
2150 1315.4 1449.4 1      1  1729  1      0     1315.4 1449.4  0.8944 0     
2539 1449.4 1726.1 1      1  1729  1      0     1449.4 1726.1 -0.1931 0     
2851 1726.1 1729.4 0      1  1729  1      1     1726.1 1729.4  0.4081 0     
------------------------------------------------------------ 
id: 9
     entry  time   status rr dtime fdeath death start  stop   x       strata
9       0.0  433.5 1      1  5110  0      0        0.0  433.5 -0.4660 0     
1006  433.5 2451.1 1      1  5110  0      0      433.5 2451.1  1.0647 0     
1659 2451.1 3629.7 1      1  5110  0      0     2451.1 3629.7 -0.2506 0     
2151 3629.7 3644.7 1      1  5110  0      0     3629.7 3644.7 -0.6748 0     
2540 3644.7 3695.8 1      1  5110  0      0     3644.7 3695.8  0.6510 0     
2852 3695.8 3890.7 1      1  5110  0      0     3695.8 3890.7 -0.2033 0     
3112 3890.7 5110.0 0      1  5110  0      0     3890.7 5110.0 -1.6981 0
\end{verbatim}

The status variable keeps track of the recurrent evnts and their type, and death the timing of 
death.


\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label=rec1,caption= ,captionpos=b}
\begin{lstlisting}
#  to fit non-parametric models with just a baseline 
xr <- phreg(Surv(entry,time,status)~cluster(id),data=rr)
dr <- phreg(Surv(entry,time,death)~cluster(id),data=rr)
par(mfrow=c(1,3))
bplot(dr,se=TRUE)
title(main="death")
bplot(xr,se=TRUE)
# robust standard errors 
rxr <-   robust.phreg(xr,fixbeta=1)
bplot(rxr,se=TRUE,robust=TRUE,add=TRUE,col=4)

# marginal mean of expected number of recurrent events 
out <- recurrentMarginal(xr,dr)
bplot(out,se=TRUE,ylab="marginal mean",col=2)
\end{lstlisting}

\begin{marginfigure}
\begin{center}
\includegraphics[width=\textwidth]{rec1.jpg}
\end{center}


\captionof{figure}{Marginal mean for number of type 1 events, rate for death (panel (a)), rate for type 1 among survivors (panel (b)), and marginal mean (panel (c)).}
\label{fig:rec}
\end{marginfigure}

We can do the same with strata 

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label=rec2,caption= ,captionpos=b}
\begin{lstlisting}
xr <- phreg(Surv(entry,time,status)~strata(strata)+cluster(id),data=rr)
dr <- phreg(Surv(entry,time,death)~strata(strata)+cluster(id),data=rr)
par(mfrow=c(1,3))
bplot(dr,se=TRUE)
title(main="death")
bplot(xr,se=TRUE)
rxr <-   robust.phreg(xr,fixbeta=1)
bplot(rxr,se=TRUE,robust=TRUE,add=TRUE,col=1:2)

out <- recurrentMarginal(xr,dr)
bplot(out,se=TRUE,ylab="marginal mean",col=1:2)
\end{lstlisting}

\begin{marginfigure}
\begin{center}
\includegraphics[width=\textwidth]{rec2.jpg}
\end{center}


\captionof{figure}{Recurrent events}
\label{fig:rec2}
\end{marginfigure}


Furhter, if we adjust for covariates for the two rates we can still do
predictions of marginal mean, what can be plotted is the baseline marginal mean, 
that is for the covariates equal to 0 for both models.

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label=rec3,caption= ,captionpos=b}
\begin{lstlisting}
# cox case
xr <- phreg(Surv(entry,time,status)~x+cluster(id),data=rr)
dr <- phreg(Surv(entry,time,death)~x+cluster(id),data=rr)
par(mfrow=c(1,3))
bplot(dr,se=TRUE)
title(main="death")
bplot(xr,se=TRUE)
rxr <- robust.phreg(xr)
bplot(rxr,se=TRUE,robust=TRUE,add=TRUE,col=1:2)

out <- recurrentMarginal(xr,dr)
bplot(out,se=TRUE,ylab="marginal mean",col=1:2)

# predictions witout se's 
outX <- recmarg(xr,dr,Xr=1,Xd=1)
bplot(outX,add=TRUE,col=3)

\end{lstlisting}

\begin{marginfigure}
\begin{center}
\includegraphics[width=\textwidth]{rec3.jpg}
\end{center}


\captionof{figure}{Recurrent events with cox models for rates.}
\label{fig:rec3}
\end{marginfigure}






\subsection*{Other marginal properties}
\label{sec:org2fcaab7}



\begin{itemize}
\item \(P(N_1^*(t) \ge k)\) 
\begin{itemize}
\item cumulative incidence of \(T_{k} = \inf \{ t: N_1^*(t)=k \}\) with competing \(D\).
\end{itemize}
\end{itemize}

We note also that \(N_1^*(t)^2\) can be written as
\begin{align*}
   \sum_{k=0}^K  \int_0^t I(D > s) I(N_1^*(s-)=k) f(k) dN_1^*(s)
\end{align*}
with \(f(k)=(k+1)^2 - k^2\), such that its mean can be written as 
\begin{align*}
	\sum_{k=0}^K \int_0^t S(s) f(k) P(N_1^*(s-)= k  | D  \geq s) E( dN_1^*(s)  | N_1^*(s-)=k, D> s) 
\end{align*}
and estimated by
\begin{align*}
\hat \mu_{1,2}(t) & = 
	\sum_{k=0}^K \int_0^t \hat S(s) f(k) 
	\frac{Y_{1\bullet}^k(s)}{Y_\bullet (s)} \frac{1}{Y_{1\bullet}^k(s)} d N_{1\bullet}^k(s)= \sum_{i=1}^n \int_0^t \hat S(s) f(N_{i1}(s-)) \frac{1}{Y_\bullet (s)} d N_{i1}(s),
\end{align*}
Compared to "product-limit" estimator for \(E( (N_1^*(t))^2 )\) 
\begin{align}
  \hat \mu_{1,2}(t) & =    \sum_{k=0}^K k^2 ( \hat F_{k}(t) - \hat F_{k+1}(t) ).
\end{align}


Probabilty  of exceeding "k"

Note also that  \(I(N_1^*(t) \geq k)\) is 
\begin{align*}
	\int_0^t I(D > s) I(N_1^*(s-)=k-1) dN_1^*(s),
\end{align*}
suggesting that its mean can be computed as
\begin{align*}
\int_0^t S(s) P(N_1^*(s-)= k-1  | D  \geq s) E( dN_1^*(s)  | N_1^*(s-)=k-1, D> s) 
\end{align*}
and estimated by 
\begin{align*}
\tilde F_k(t) = \int_0^t \hat S(s)  \frac{Y_{1\bullet}^{k-1}(s)}{Y_\bullet (s)} 
          	\frac{1}{Y_{1\bullet}^{k-1}(s)} d N_{1\bullet}^{k-1}(s).
\end{align*}


\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label=rec4,caption= ,captionpos=b}
\begin{lstlisting}
cor.mat <- corM <- rbind(c(1.0, 0.6, 0.9), c(0.6, 1.0, 0.5), c(0.9, 0.5, 1.0))
rr <- simRecurrent(1000,base1,cumhaz2=base4,death.cumhaz=ddr)
rr <-  count.history(rr)
dtable(rr,~death+status)

oo <- prob.exceedRecurrent(rr,1)
bplot(oo)
\end{lstlisting}

\begin{marginfigure}
\begin{center}
\includegraphics[width=\textwidth]{rec4.jpg}
\end{center}


\captionof{figure}{Recurrent events: probability of exceeding k events}
\label{fig:rec4}
\end{marginfigure}


\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label=rec4,caption= ,captionpos=b}
\begin{lstlisting}
cor.mat <- corM <- rbind(c(1.0, 0.6, 0.9), c(0.6, 1.0, 0.5), c(0.9, 0.5, 1.0))
rr <- simRecurrent(1000,base1,cumhaz2=base4,death.cumhaz=ddr)
rr <-  count.history(rr)
dtable(rr,~death+status)

oo <- prob.exceedRecurrent(rr,1)
bplot(oo)
\end{lstlisting}

\begin{marginfigure}
\begin{center}
\includegraphics[width=\textwidth]{rec4MV.jpg}
\end{center}


\captionof{figure}{Recurrent events: probability of exceeding k events}
\label{fig:rec4MV}
\end{marginfigure}


Mean and variance of number of recurrent events 
\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label=rec4MV,caption= ,captionpos=b}
\begin{lstlisting}
par(mfrow=c(1,2))
with(oo,plot(time,mu,col=2,type="l"))
#
with(oo,plot(time,varN,type="l"))
\end{lstlisting}




\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label=rec4Bi,caption= ,captionpos=b}
\begin{lstlisting}
# Bivariate probability of exceeding 
oo <- prob.exceedBiRecurrent(rr,1,2,exceed1=c(1,5,10),exceed2=c(1,2,3))
with(oo, matplot(time,pe1e2,type="s"))
nc <- ncol(oo$pe1e2)
legend("topleft",legend=colnames(oo$pe1e2),lty=1:nc,col=1:nc)
\end{lstlisting}


\begin{marginfigure}
\begin{center}
\includegraphics[width=\textwidth]{rec4Bi.jpg}
\end{center}



\captionof{figure}{Recurrent events: probability of exceeding k events}
\label{fig:rec4Bi}
\end{marginfigure}


\subsection*{Dependence between events: Covariance}
\label{sec:org0e76a6d}

Covariance among two types of events 
\begin{align}
\rho(t) &  = \frac{ E(N_1^*(t) N_2^*(t) )  - \mu_1(t) \mu_2(t) }{ \mbox{sd}(N_1^*(t)) \mbox{sd}(N_2^*(t)) }
\end{align}
where  
\begin{itemize}
\item \(E(N_1^*(t) N_2^*(t))\).
\end{itemize}

\begin{align*}
  E(N_1^*(t) N_2^*(t)) &   = E( \int_0^t N_1^*(s-) dN_2^*(s) ) + E( \int_0^t N_2^*(s-) dN_1^*(s) ) 
\end{align*}

Recall that \(N_1^*(t \wedge D)\)  and \(N_2^*(t \wedge D)\). 


\begin{align*}
E(\int_0^t N_1^*(s-) dN_2^*(s) ) & = \sum_k E( \int_0^t k I(N_1^*(s-)=k) I(D \geq s)  dN_2^*(s) ) 
\end{align*}
\begin{align*}
= \sum_k \int_0^t S(s) k P(N_1^*(s-)= k  | D  \geq s) E( dN_2^*(s)  | N_1^*(s-)=k, D \geq s) 
\end{align*}
estimated by 
\begin{align*}
  & \sum_k \int_0^t \hat S(s) k \frac{Y_1^k(s)}{Y_\bullet (s)} \frac{1}{Y_1^k(s)} d \tilde N_{2,k}(s),
\end{align*}
\begin{itemize}
\item \(Y_j^k(t) = \sum Y_i(t) I( N_{ji}^*(s-)=k)\) for \(j=1,2\),
\item \(\tilde N_{j,k}(t) = \sum_i \int_0^t I(N_{ij^o}(s-)=k) dN_{ij}(s)\)
\end{itemize}

Estimate of \$ E(N\(_{\text{1}}^{\text{*}}\)(t) N\(_{\text{2}}^{\text{*}}\)(t))\$ 
\begin{align*}
  \sum_k \int_0^t \hat S(s) k \frac{Y_1^k(s)}{Y_\bullet (s)} \frac{1}{Y_1^k(s)} d \tilde N_{2,k}(s) +
   \sum_k \int_0^t \hat S(s) k \frac{Y_2^k(s)}{Y_\bullet (s)} \frac{1}{Y_2^k(s)} d \tilde N_{1,k}(s).
\end{align*}


\begin{itemize}
\item Without terminating event covariance is useful nonpar measure
\item With terminating event dependence generated by terminating event.
\item In reality what is of interest would be
independence among survivors 
\begin{itemize}
\item if \(N_1\) not predicitive for \(N_2\)
\end{itemize}
\begin{align}
   E( dN_2^*(t)  | N_1^*(t-)=k, D \geq t) =  E( dN_2^*(t)  | D \geq t) 
\end{align}
\begin{itemize}
\item if \(N_2\) not predicitive for \(N_1\)
\end{itemize}
\begin{align}
   E( dN_1^*(t)  | N_2^*(t-)=k, D \geq t) =  E( dN_1^*(t)  | D \geq t) 
\end{align}
\end{itemize}


If the two processes are independent among survivors then 
\begin{align}
 E( dN_2^*(t)  | N_1^*(t-)=k, D \geq t) =  E( dN_2^*(t)  | D \geq t) 
\end{align}
so 
\begin{align*}
  E( \int_0^t N_1^*(s-) dN_2^*(s) )  & =  \int_0^t S(s) E(N_1^*(s-) | D  \geq s) E( dN_2^*(s)  | D \geq s) 
\end{align*}
and 
\begin{align*}
  \int_0^t \hat S(s) \{  \sum_k k \frac{Y_1^k(s)}{Y_\bullet (s)} \} \frac{1}{Y_\bullet (s)} dN_{2\bullet}(s),
\end{align*}
where \(N_{j\bullet}(t) = \sum_i \int_0^t dN_{j,i}(s)\).

Under the independence \(E(N_1^*(t) N_2^*(t))\)  is estimated 
\begin{align*}
  \int_0^t \hat S(s) \{  \sum_k k \frac{Y_1^k(s)}{Y_\bullet (s)} \} \frac{1}{Y_\bullet (s)} dN_{2\bullet}(s) 
    + \int_0^t \hat S(s) \{  \sum_k k \frac{Y_2^k(s)}{Y_\bullet (s)} \} \frac{1}{Y_\bullet (s)} dN_{1\bullet}(s).
\end{align*}


Both estimators, \(\hat E(N_1^*(t) N_2^*(t))\) and \(\hat E_I(N_1^*(t) N_2^*(t))\),  as well as 
\(\hat E(N_1^*(t))\) and \(\hat E(N_2^*(t))\), 
have asymptotic expansions that can be written as a sum of iid processes, similarly to the arguments 
of Ghosh \& Lin 2000, \(\sum_i \Psi_i(t)\).  

We can thus estimate the standard errors and of the estimators and their difference  
\(\hat E(N_1^*(t) N_2^*(t))- \hat E_I(N_1^*(t) N_2^*(t))\). 

Terms for 
\begin{itemize}
\item N1 -> N2 : \(E( \int_0^t N_1^*(s-) dN_2^*(s) )\)
\item N2 -> N1 : \(E( \int_0^t N_2^*(s-) dN_1^*(s) )\)
\end{itemize}


\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label=rec5,caption= ,captionpos=b}
\begin{lstlisting}
rr$strata <- 1
dtable(rr,~death+status)

covrp <- covarianceRecurrent(rr,1,2,status="status",death="death",
			start="entry",stop="time",id="id",names.count="Count")
par(mfrow=c(1,3)) 
plot(covrp)

# with strata, each strata in matrix column, provides basis for fast Bootstrap
covrpS <- covarianceRecurrentS(rr,1,2,status="status",death="death",
	start="entry",stop="time",strata="strata",id="id",names.count="Count")
\end{lstlisting}

\begin{marginfigure}
\begin{center}
\includegraphics[width=\textwidth]{rec5.jpg}
\end{center}

\captionof{figure}{Covariance between events}
\label{fig:rec5}
\end{marginfigure}


\subsection*{Bootstrap standard errors for terms}
\label{sec:orgc567a4f}

First fitting the model again to get our estimates of interst, and then 
computing them for some specific time-points  

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
times <- seq(500,5000,500)

coo1 <- covarianceRecurrent(rr,1,2,status="status",start="entry",stop="time")
#
mug <- Cpred(cbind(coo1$time,coo1$EN1N2),times)[,2]
mui <- Cpred(cbind(coo1$time,coo1$EIN1N2),times)[,2]
mu2.1 <- Cpred(cbind(coo1$time,coo1$mu2.1),times)[,2]
mu2.i <- Cpred(cbind(coo1$time,coo1$mu2.i),times)[,2]
mu1.2 <- Cpred(cbind(coo1$time,coo1$mu1.2),times)[,2]
mu1.i <- Cpred(cbind(coo1$time,coo1$mu1.i),times)[,2]
cbind(mu2.1,mu2.i)
cbind(mu1.2,mu1.i)
\end{lstlisting}

\begin{verbatim}
           mu2.1      mu2.i
 [1,] 0.04101096 0.03656491
 [2,] 0.09303668 0.08572694
 [3,] 0.22613687 0.21906324
 [4,] 0.35727148 0.34562539
 [5,] 0.60258982 0.59071900
 [6,] 0.80089841 0.79020220
 [7,] 1.03031183 1.03424672
 [8,] 1.16860632 1.16686717
 [9,] 1.25782175 1.25105963
[10,] 1.38716306 1.40250244
           mu1.2      mu1.i
 [1,] 0.03501045 0.03259566
 [2,] 0.08803686 0.08526834
 [3,] 0.16709531 0.16634828
 [4,] 0.27720710 0.29485672
 [5,] 0.38034407 0.41985665
 [6,] 0.53057410 0.56459585
 [7,] 0.69387628 0.72234676
 [8,] 0.87226707 0.88771625
 [9,] 0.96949736 0.99728527
[10,] 1.06074066 1.06854228
\end{verbatim}

To get the bootstrap standard errors there is a quick memory demanding function 
(with S for speed and strata)  BootcovariancerecurrenceS
and slow function that goes through the loops in R Bootcovariancerecurrence. 

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
bt1 <- BootcovariancerecurrenceS(rr,1,2,status="status",start="entry",stop="time",K=100,times=times)
#bt1 <- BootcovariancerecurrenceS(rr,1,2,status="status",start="entry",stop="time",K=K,times=times)

output <- list(bt1=bt1,mug=mug,mui=mui,
bse.mug=bt1$se.mug,bse.mui=bt1$se.mui,
dmugi=mug-mui,
bse.dmugi=apply(bt1$EN1N2-bt1$EIN1N2,1,sd),
mu2.1 = mu2.1 , mu2.i = mu2.i , dmu2.i=mu2.1-mu2.i,
mu1.2 = mu1.2 , mu1.i = mu1.i , dmu1.i=mu1.2-mu1.i,
bse.mu2.1=apply(bt1$mu2.i,1,sd), bse.mu2.1=apply(bt1$mu2.1,1,sd),
bse.dmu2.i=apply(bt1$mu2.1-bt1$mu2.i,1,sd),
bse.mu1.2=apply(bt1$mu1.2,1,sd), bse.mu1.i=apply(bt1$mu1.i,1,sd),
bse.dmu1.i=apply(bt1$mu1.2-bt1$mu1.i,1,sd)
)
\end{lstlisting}

We then look at the test for overall dependence  in the different time-points. 
We here have no suggestion of dependence.

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
tt  <- output$dmugi/output$bse.dmugi
cbind(times,2*(1-pnorm(abs(tt))))
\end{lstlisting}

\begin{verbatim}
      times          
 [1,]   500 0.3572253
 [2,]  1000 0.4577012
 [3,]  1500 0.7136132
 [4,]  2000 0.7956959
 [5,]  2500 0.3837459
 [6,]  3000 0.5134406
 [7,]  3500 0.4209237
 [8,]  4000 0.7632914
 [9,]  4500 0.6836682
[10,]  5000 0.6598813
\end{verbatim}

We can also take out the specific components for whether \(N_1\) is predictive for 
\(N_2\) and vice versa.  We here have no suggestion of dependence.

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label= ,caption= ,captionpos=b}
\begin{lstlisting}
t21  <- output$dmu1.i/output$bse.dmu1.i 
t12  <- output$dmu2.i/output$bse.dmu2.i 
cbind(times,2*(1-pnorm(abs(t21))),2*(1-pnorm(abs(t12))))
\end{lstlisting}

\begin{verbatim}
      times                     
 [1,]   500 0.71706002 0.3918872
 [2,]  1000 0.81454942 0.3202626
 [3,]  1500 0.95715638 0.6006314
 [4,]  2000 0.21300406 0.4942293
 [5,]  2500 0.02182129 0.6086128
 [6,]  3000 0.11688970 0.6805457
 [7,]  3500 0.25587816 0.8965495
 [8,]  4000 0.63373150 0.9578608
 [9,]  4500 0.41743073 0.8548733
[10,]  5000 0.83041113 0.6805618
\end{verbatim}


We finally plot the boostrap samples 

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label=rec6,caption= ,captionpos=b}
\begin{lstlisting}
par(mfrow=c(1,2))
matplot(bt1$time,bt1$EN1N2,type="l",lwd=0.3)
matplot(bt1$time,bt1$EIN1N2,type="l",lwd=0.3)
\end{lstlisting}

\begin{marginfigure}
\begin{center}
\includegraphics[width=\textwidth]{rec6.jpg}
\end{center}

\captionof{figure}{Bootstrap samples}
\label{fig:rec6}
\end{marginfigure}


\subsection*{Looking at other simulations with dependence}
\label{sec:orga99167a}


Using the normally distributed random effects we plot 4 different settings. We have variance \(0.5\) for all
random effects and change the correlation. We let the correlation between the random effect associated with
\(N_1\) and \(N_2\) be denoted \(\rho_{12}\) and the correlation between the random effects 
associated between \(N_j\) and \(D\) the terminal event be denoted as \(\rho_{j3}\), and organize all correlation
in a vector \(\rho=(\rho_{12},\rho_{13},\rho_{23})\).

\begin{itemize}
\item Scenario I \(\rho=(0,0.0,0.0)\) Independence among all efects.
\item Scenario II \(\rho=(0,0.5,0.5)\) Independence among survivors but dependence on terminal event
\item Scenario III \(\rho=(0.5,0.5,0.5)\) Positive dependence among survivors and dependence on terminal event
\item Scenario IV \(\rho=(-0.4,0.5,0.5)\) Negative dependence among survivors and positive dependence on terminal event
\end{itemize}

\lstset{numbers=left,numberstyle=\ttfamily\tiny\color{gray},language=r,label=rec7,caption= ,captionpos=b}
\begin{lstlisting}
par(mfrow=c(2,2))

data(base1cumhaz)
data(base4cumhaz)
data(drcumhaz)
dr <- drcumhaz
base1 <- base1cumhaz
base4 <- base4cumhaz

var.z <- c(0.5,0.5,0.5)
# death related to  both causes in same way 
cor.mat <- corM <- rbind(c(1.0, 0.0, 0.0), 
			 c(0.0, 1.0, 0.0), 
			 c(0.0, 0.0, 1.0))
rr <- simRecurrentII(3000,base1,base4,death.cumhaz=dr,var.z=var.z,cor.mat=cor.mat,dependence=2)
rr <- count.history(rr,types=1:2)
cor(attr(rr,"z"))
coo <- covarianceRecurrent(rr,1,2,status="status",start="entry",stop="time")
par(mfrow=c(2,2))
with(coo, {
      plot(time, EN1N2, type = "l", lwd = 2,lty=1,ylab="",xlab="time (a)")
      lines(time, EN1EN2, col = 2, lwd = 2,lty=2)
      lines(time, EIN1N2, col = 3, lwd = 2,lty=3)
  })
  legend("topleft", c("E(N1N2)", "E(N1) E(N2) ", "E_I(N1 N2)-independence"),lty = 1:3, col = 1:3)
  title(main ="Scenario I")


var.z <- c(0.5,0.5,0.5)
# death related to  both causes in same way 
cor.mat <- corM <- rbind(c(1.0, 0.0, 0.5), 
			 c(0.0, 1.0, 0.5), 
			 c(0.5, 0.5, 1.0))
rr <- simRecurrentII(3000,base1,base4,death.cumhaz=dr,
		      var.z=var.z,cor.mat=cor.mat,dependence=2)
rr <- count.history(rr,types=1:2)
coo <- covarianceRecurrent(rr,1,2,status="status",start="entry",stop="time")
with(coo, {
      plot(time, EN1N2, type = "l", lwd = 2,lty=1,ylab="",xlab="time (b)")
      lines(time, EN1EN2, col = 2, lwd = 2,lty=2)
      lines(time, EIN1N2, col = 3, lwd = 2,lty=3)
  })
  legend("topleft", c("E(N1N2)", "E(N1) E(N2) ", "E_I(N1 N2)-independence"),lty = 1:3, col = 1:3)
  title(main ="Scenario II")

var.z <- c(0.5,0.5,0.5)
# positive dependence for N1 and N2 all related in same way
cor.mat <- corM <- rbind(c(1.0, 0.5, 0.5), 
			 c(0.5, 1.0, 0.5), 
			 c(0.5, 0.5, 1.0))
rr <- simRecurrentII(3000,base1,base4,death.cumhaz=dr,
		      var.z=var.z,cor.mat=cor.mat,dependence=2)
rr <- count.history(rr,types=1:2)
coo <- covarianceRecurrent(rr,1,2,status="status",start="entry",stop="time")
with(coo, {
      plot(time, EN1N2, type = "l", lwd = 2,lty=1,ylab="",xlab="time (d)")
      lines(time, EN1EN2, col = 2, lwd = 2,lty=2)
      lines(time, EIN1N2, col = 3, lwd = 2,lty=3)
  })
  legend("topleft", c("E(N1N2)", "E(N1) E(N2) ", "E_I(N1 N2)-independence"),lty = 1:3, col = 1:3)
  title(main ="Scenario III")


var.z <- c(0.5,0.5,0.5)
# negative dependence for N1 and N2 all related in same way
cor.mat <- corM <- rbind(c(1.0, -0.4, 0.5), 
			 c(-0.4, 1.0, 0.5), 
			 c(0.5, 0.5, 1.0))
rr <- simRecurrentII(3000,base1,base4,death.cumhaz=dr,
		      var.z=var.z,cor.mat=cor.mat,dependence=2)
rr <- count.history(rr,types=1:2)
coo <- covarianceRecurrent(rr,1,2,status="status",start="entry",stop="time")
with(coo, {
      plot(time, EN1N2, type = "l", lwd = 2,lty=1,ylab="",xlab="time (d)")
      lines(time, EN1EN2, col = 2, lwd = 2,lty=2)
      lines(time, EIN1N2, col = 3, lwd = 2,lty=3)
})
legend("topleft", c("E(N1N2)", "E(N1) E(N2) ", "E_I(N1 N2)-independence"),lty = 1:3, col = 1:3)
title(main ="Scenario IV")
\end{lstlisting}

\begin{figure}
\begin{center}
\includegraphics[width=\textwidth]{rec7.jpg}
\end{center}

\captionof{figure}{Bootstrap samples}
\label{fig:rec7}
\end{figure}
\end{document}