%\VignetteIndexEntry{Multistate example from Crowther and Lambert - with multiple timescales}
\SweaveOpts{results=verbatim,keep.source=TRUE,include=FALSE,eps=FALSE}
\documentclass[a4paper,dvipsnames,twoside,12pt]{article}
<<echo=FALSE>>=
options( width=90,
#        prompt=" ", continue=" ", # Absence of prompts makes it easier for
                                   # students to copy from the final pdf document
         SweaveHooks=list( fig=function()
         par(mar=c(3,3,1,1),mgp=c(3,1,0)/1.6,las=1,bty="n") ) )
@ %
% ----------------------------------------------------------------------
% General information for the title page and the page headings
\newcommand{\Title}{Multistate example from\\ Crowther \& Lambert\\---
with multiple timescales}
\newcommand{\Tit}{Crowther \& Lambert}
\newcommand{\Version}{Version 5}
\newcommand{\Dates}{February 2018}
\newcommand{\Where}{SDCC}
\newcommand{\Homepage}{\url{http://bendixcarstensen.com/AdvCoh} }
\newcommand{\Faculty}{\begin{tabular}{rl}
Bendix Carstensen
  & Steno Diabetes Center Copenhagen, Gentofte, Denmark\\
  & {\small \& Department of Biostatistics,
               University of Copenhagen} \\
  & \mailto{bcar0029@regionh.dk} \quad \mailto{b@bxc.dk} \\
  & \url{http://BendixCarstensen.com} \\[1em]
                      \end{tabular}}

\input{toparticle}
\renewcommand{\rwpre}{./BrCaMS}

\section{Introduction}

This is a re-do (and extension) of (parts of) the example from the
short-titled paper by Crowther \& Lambert \cite{Crowther.2017}. The
data provided by the authors are available as the data set
\texttt{BrCa} in the \texttt{Epi} package in a slightly modified form,
where dates of relapse, metastasis and death are only non-\texttt{NA} for
those that actually do see the events.

First we load the relevant packages and then the example data from the
\texttt{Epi} package:
<<read>>=
library( Epi )
library( popEpi )
load( file="./BrCa.rda" )
# data( BrCa )
head( BrCa )
@ %

\subsection{Setting up a \texttt{Lexis} object for the follow-up}

Now we are in a position to set up the survival data as a Lexis
object. The age and date of entry are only given as integral years, so in
order to make the data credible we add a random number between 0 and 1
to mimic a real age and date at entry. We define the time scale
\texttt{tfd} (\texttt{t}ime \texttt{f}rom \texttt{d}iagnosis) as time
since entry into the study:
<<defLex>>=
set.seed( 1952 )
Lbc <- Lexis( entry = list( tfd = 0,
                              A = age  + runif(nrow(BrCa)),
                              P = year + runif(nrow(BrCa)) ),
               exit = list( tfd = tox ),
        exit.status = xst,
                 id = pid,
               data = BrCa )
summary( Lbc )
names( Lbc )
@ %
Now we want to cut the follow up at the times of relapse (including
metastasis), but keep track of whether a person died with or without
relapse, so we set \texttt{split.states} to true, and since time since
relapse is presumably of interest too we ask for that time scale to be
defined as well (using the argument \texttt{new.scale}):
<<cut>>=
Rbc <- cutLexis( Lbc, 
                 cut = pmin( Lbc$tor, Lbc$tom, na.rm=TRUE ),
           timescale = "tfd",
    precursor.states = "Alive",
           new.state = "Rel",
        split.states = TRUE,
           new.scale = "tfr" )
summary( Rbc, timeScale = TRUE )
@ % $
From the summary we see that the transitions to death are to different
states, depending on whether a relapse had occurred or not (this is
the result of \texttt{split.states}), this will eventually allow us to
assess the cumulative risk of relapse. Moreover \texttt{new.scale}
ensured that a new time scale, \texttt{tfr}, time from relapse has
been added to the \texttt{Lexis} object.

We can illustrate the transitions by a plot that gives a convenient
overview of transitions: 
<<box-r,fig=TRUE, height=5>>=
boxes( Rbc, boxpos=list(x=c(15,15,85,85),
                        y=c(85,15,85,15)),
            show.BE=TRUE, scale.R=100, )
@ % 
\insfig{box-r}{0.7}{Transitions in the correctly set up multistate
  model for the breast cancer survival dataset. Numbers \emph{in} the
  boxes are person-years and (at the bottom) the number of persons
  starting resp. ending their follow-up in each state. Numbers on the
  arrows are the number of transitions and transition rates per 100
  PY (by the \textrm{\tt scale.R} argument).}

\section{Modeling rates}

In line with Crowther and Lambert we now model the transition
rates. To this end we first split the data in smaller chunks of length
1 month --- with some 20,000 PY we would expect to have some
250,000 records:
<<>>=
system.time(
Sbc <- splitLexis( Rbc, breaks=seq(0,100,1/12), "tfd" ) )
summary( Sbc )
@ %
In the \texttt{popEpi} package is a similar function with more elegant
syntax and somewhat faster particularly for large data sets:
<<>>=
system.time(
Sbc <- splitMulti( Rbc, tfd=seq(0,100,1/12) ) )
summary( Sbc )
@ %

\subsection{Stacking?}

We could model all 3 rates jointly by stacking the data --- the function
\texttt{stack.Lexis} would do this, and create variables
\texttt{lex.Tr} (transition type) and \texttt{lex.Fail} (event indicator):
<<>>=
Stbc <- stack( Sbc )
round( ftable( xtabs( cbind(lex.Fail,lex.dur) ~ lex.Tr + lex.Xst, 
                      data=Stbc ),
               row.vars=c(3,1),
       1 ) )
@ %
However, stacking data is needed only when all transitions are to be
modeled jointly, or more specifically, when more than one transition
\emph{out} of a given state are modeled jointly. This type of modeling
is rarely wanted, since rates of different types of events (in this
case relapse and death) are unlikely to depend on the same variable in
the same way.

It is much more likely that different mortality rates depend on
covariates in the same way --- in this case that mortality from
``Alive'' and from ``Rel'' depend on time since entry and on the
clinical parameters the same way. Additionally we may take time since
relapse into account.

In such an instance, the original \texttt{Lexis} object where the
total follow-up time is represented exactly once in \texttt{lex.dur},
will suffice as database for the analysis, because at most \emph{one}
transition out of each state is considered. So we shall leave aside
the stacking, and model the three rates separately.

\subsection{Initial model by C \& L}

The initial approach is basically to model each of the transitions
separately; here we use natural splines with 4 knots placed at the
quantiles of the transition times (we refer to the transitions as
\texttt{ad} (\texttt{a}live to \texttt{d}ead),
\texttt{ar} (\texttt{a}live to \texttt{r}elapse),
\texttt{rd} (\texttt{r}elapse to \texttt{d}ead). For the sake of
completeness we also compute knots on the scale of time since relapse,
as well as for the (fixed) difference between \texttt{tfd} and
\texttt{tfr} (the time \emph{at} relapse --- note that we do not
construct a separate variable for this): 
<<>>=
( kd.ad <- with( subset( Sbc, lex.Cst=="Alive" & lex.Xst=="Dead"),
                 quantile( tfd+lex.dur, probs=(1:4-0.5)/4) ) )
( kd.ar <- with( subset( Sbc, lex.Cst=="Alive" & lex.Xst=="Rel"),
                 quantile( tfd+lex.dur, probs=(1:4-0.5)/4) ) )
( kd.rd <- with( subset( Sbc, lex.Cst=="Rel" & lex.Xst=="Dead(Rel)"),
                 quantile( tfd+lex.dur, probs=(1:4-0.5)/4) ) )
( kr.rd <- with( subset( Sbc, lex.Cst=="Rel" & lex.Xst=="Dead(Rel)"),
                 quantile( tfr+lex.dur, probs=(1:4-0.5)/4) ) )
( ka.rd <- with( subset( Sbc, lex.Cst=="Rel" & lex.Xst=="Dead(Rel)"),
                 quantile( tfd-tfr, probs=(1:4-0.5)/4) ) )
@ %
With these vectors of knots in place we can fit models for the three
rates --- note the similarity of the modeling code for the different
models and the immediate readability of what is being modeled;
\texttt{lex.Cst} is used to define the risk set (using
\texttt{subset}) and \texttt{lex.Xst} to define the event type:
<<>>=
m.ad <- glm( (lex.Xst=="Dead") ~ Ns( tfd, knots=kd.ad ),
             offset = log( lex.dur ),
             family = poisson,
               data = subset( Sbc, lex.Cst=="Alive" ) )
m.ar <- glm( (lex.Xst=="Rel") ~ Ns( tfd, knots=kd.ar ),
             offset = log( lex.dur ),
             family = poisson,
               data = subset( Sbc, lex.Cst=="Alive" ) )
m.rd <- glm( (lex.Xst=="Dead(Rel)") ~ Ns( tfd, knots=kd.rd ),
             offset = log( lex.dur ),
             family = poisson,
               data = subset( Sbc, lex.Cst=="Rel" ) )
x.rd <- update( m.rd, . ~ . + Ns( tfr, knots=kr.rd) )
r.rd <- update( x.rd, . ~ . - Ns( tfd, knots=kd.rd) )
anova( m.rd, x.rd, r.rd, test="Chisq" )
@ %
We see that the mortality rates in relapse depends strongly on the
time since relapse, a deviance reduction of 77 on 3 df! Ditching the
effect of \texttt{tfd} is clearly neither a feasible option with a
deviance difference of 198 on 3 df. We shall deal with this extension
later.

First we turn to the transition rates as function of time since
diagnosis. Note that since the \texttt{lex.dur} is in units of PY,
setting the value of it (as a covariate) to 100, means that we get
the rates in units of 100 PY --- basically rates in \% per year:
<<rate-pr>>=
nd <- data.frame( tfd = seq(0,15,0.1),
              lex.dur = 100 )
ad.rate <- ci.pred( m.ad, nd )
ar.rate <- ci.pred( m.ar, nd )
rd.rate <- ci.pred( m.rd, nd )
@ %
We can plot the three sets of estimated rates in the same graph:
<<pr-pl,fig=TRUE>>=
clr <- rainbow(3) # ; yl <- c(0.03,60)
matplot( nd$tfd, cbind( ad.rate,
                        ar.rate,
                        rd.rate ),
         type="l", lty=1, lwd=c(3,1,1), col=rep(clr,each=3), las=1,
         log="y", xlab="Time since diagnosis (years)",
                  ylab="Rate per 100 PY" )
text( par("usr")[2]*0.95, (10^par("usr"))[3]*1.4^(1:3), 
      c("A->D","A->R","R->D"), col=clr, adj=1, font=2 )
matlines( nd$tfd, ci.ratio( rd.rate, ad.rate ), 
          lty=1, lwd=c(3,1,1), col=gray(0.6) )
abline( h=1, col=gray(0.6) )
@ %
\insfig{pr-pl}{0.8}{Transition rates as function of time since
  diagnosis, the gray line is the mortality rate-ratio between persons
  with and without relapse --- it seems as if the earlier the relapse,
  the higher the impact on mortality.}

From the graph in figure \ref{fig:pr-pl} we see that the occurrence of
relapse almost doubles over the first two years and then decreases. We
also observe that the mortality RR between persons with relapse and
those without decreases from extremely high to about 5, a combination
of decreasing mortality among persons with relapse and an increasing
mortality among persons without relapse.

\section{The two time scales --- and their difference}

We noted that the model \texttt{x.rd} above with effects of both time
since diagnosis and time since relapse represented a substantial
improvement over the models with only one of these time-scales.

We could expand this model further with an effect of time \emph{at}
relapse, $\mathtt{tfd}-\mathtt{tfr}$:
<<>>=
xx.rd <- update( x.rd, . ~ . + Ns( tfd-tfr, knots=ka.rd) )
anova( m.rd, x.rd, xx.rd, test="Chisq" )
@ %
We see there is a formally statistically significant effect of time at
relapse, but the deviance change is much smaller than for the two
timescales.

What we are doing here is adding interactions between timescales,
popularly known as ``testing for non-proportionality''. Adding time
since relapse as a time scale is one extension of the model with
proportional mortality rates between persons with and without relapse,
by letting the HR dpend on time since relapse. A further extension is
to add an effect of the difference of the two is yet another
interaction term.

The tests are however not particularly relevant; a considerably large
dataset as the current may yield statistical significance where no
clinically relevant significant effects are present. Therefore,
testing of proportionality must necessarily be supported by dispays of
the \emph{shape} of the interactions.

We can show how the addition of time since relapse and time at relapse
affects the estimated mortality by showing mortality after relapse as
a function of time since diagnosis for different times of relapse ---
by showing curves starting at the times of relapse.
<<rd-int,fig=TRUE>>=
nd <- data.frame( expand.grid( tfd=c(NA,seq(0,15,0.1)), 
                               tad=c(0,0.5,1,2,3,5,8) ), 
                  lex.dur=100 )
nd <- subset( transform( nd, tfr = tfd - tad ), tfr>=0 | is.na(tfr) )
head( nd )
matplot( nd$tfd, cbind( ci.pred( x.rd, nd )[,1],
                        ci.pred(xx.rd, nd )[,1] ),
         type="l", lty=c("solid","22"), lend="butt", 
         lwd=3, col=clr[3], las=1,
         log="y", xlab="Time since diagnosis (years)",
                  ylab="Mortality rate per 100 PY" )
matlines( seq(0,15,0.1), rd.rate,
          type="l", lwd=c(3,1,1), lty=1, col=gray(0.3) )
@ %
\insfig{rd-int}{0.8}{Estimated mortality among women in relapse. The
  blue lines represent mortality for women relapsed at 0, 0.5, 1, 2,
  3, 5, 8 years after diagnosis. The broken lines are predictions from
  the model where the time \emph{at} relapse is modeled too. The gray
  line is from the model where only time since diagnosis is included
  (``proportional hazards model''), corresponding to the blue line in
  figure \ref{fig:pr-pl}.}

From figure \ref{fig:rd-int} we see that the simple model completely
misses to describe the initial increase in mortality, and that the
model without the time \emph{at} relapse overestimates the mortality
among women with early relapse.

\section{Including covariates}

Following the example in the paper, we include the available
covariates in the models:
<<>>=
c.ar <- update( m.ar, . ~ . + age + size + nodes + pr.tr + hormon )
c.ad <- update( m.ad, . ~ . + age + size + nodes + pr.tr + hormon )
c.rd <- update( m.rd, . ~ . + age + size + nodes + pr.tr + hormon )
cx.rd<- update(xx.rd, . ~ . + age + size + nodes + pr.tr + hormon )
@ %

\subsection{Testing for interaction with time}

Further, we can now include terms allowing for interaction between
covariates and time since diagnosis (often termed
``non-proportionality'' in the vein of never foregoing an opportunity
to invent yet another term for a well-known concept). It is not
entirely clear from the models shown in the paper how the
non-proportionality is taken into account, but here we have used the
product of the variable with log-time + 0.5 years. In total we have 4
models and 5 variables that we can test for interaction with
\texttt{tfd}, so we set up an array to hold the p-values for the
tests.
<<int-test,eval=FALSE>>=
int.test <- NArray( list( model=c("c.ar","c.ad","c.rd","cx.rd"),
                            var=c("age","size","nodes","pr.tr","hormon"),
                           what=c("d.f.","Dev","P") ) )
str( int.test )

int.test[1,1,]<-as.numeric(anova( c.ar,update( c.ar,.~.+log(tfd+0.5):age   ),test="Chisq")[2,3:5])
int.test[1,2,]<-as.numeric(anova( c.ar,update( c.ar,.~.+log(tfd+0.5):size  ),test="Chisq")[2,3:5])
int.test[1,3,]<-as.numeric(anova( c.ar,update( c.ar,.~.+log(tfd+0.5):nodes ),test="Chisq")[2,3:5])
int.test[1,4,]<-as.numeric(anova( c.ar,update( c.ar,.~.+log(tfd+0.5):pr.tr ),test="Chisq")[2,3:5])
int.test[1,5,]<-as.numeric(anova( c.ar,update( c.ar,.~.+log(tfd+0.5):hormon),test="Chisq")[2,3:5])

int.test[2,1,]<-as.numeric(anova( c.ad,update( c.ad,.~.+log(tfd+0.5):age   ),test="Chisq")[2,3:5])
int.test[2,2,]<-as.numeric(anova( c.ad,update( c.ad,.~.+log(tfd+0.5):size  ),test="Chisq")[2,3:5])
int.test[2,3,]<-as.numeric(anova( c.ad,update( c.ad,.~.+log(tfd+0.5):nodes ),test="Chisq")[2,3:5])
int.test[2,4,]<-as.numeric(anova( c.ad,update( c.ad,.~.+log(tfd+0.5):pr.tr ),test="Chisq")[2,3:5])
int.test[2,5,]<-as.numeric(anova( c.ad,update( c.ad,.~.+log(tfd+0.5):hormon),test="Chisq")[2,3:5])

int.test[3,1,]<-as.numeric(anova( c.rd,update( c.rd,.~.+log(tfd+0.5):age   ),test="Chisq")[2,3:5])
int.test[3,2,]<-as.numeric(anova( c.rd,update( c.rd,.~.+log(tfd+0.5):size  ),test="Chisq")[2,3:5])
int.test[3,3,]<-as.numeric(anova( c.rd,update( c.rd,.~.+log(tfd+0.5):nodes ),test="Chisq")[2,3:5])
int.test[3,4,]<-as.numeric(anova( c.rd,update( c.rd,.~.+log(tfd+0.5):pr.tr ),test="Chisq")[2,3:5])
int.test[3,5,]<-as.numeric(anova( c.rd,update( c.rd,.~.+log(tfd+0.5):hormon),test="Chisq")[2,3:5])

int.test[4,1,]<-as.numeric(anova(cx.rd,update(cx.rd,.~.+log(tfd+0.5):age   ),test="Chisq")[2,3:5])
int.test[4,2,]<-as.numeric(anova(cx.rd,update(cx.rd,.~.+log(tfd+0.5):size  ),test="Chisq")[2,3:5])
int.test[4,3,]<-as.numeric(anova(cx.rd,update(cx.rd,.~.+log(tfd+0.5):nodes ),test="Chisq")[2,3:5])
int.test[4,4,]<-as.numeric(anova(cx.rd,update(cx.rd,.~.+log(tfd+0.5):pr.tr ),test="Chisq")[2,3:5])
int.test[4,5,]<-as.numeric(anova(cx.rd,update(cx.rd,.~.+log(tfd+0.5):hormon),test="Chisq")[2,3:5])
save( int.test, file="int-test.Rda")
@ %
<<>>=
load( file="int-test.Rda")
round( int.test[,,2], 2 )
round( int.test[,,3], 4 )
round( int.test[,,1], 0 )
@ %
Thus it seems that there are interactions between time from diagnosis
and progesterone for all transition rates, and that relapse rates
additionally have interactions between time from diagnosis and size
and hormone therapy. The p-values would of course have looked slightly
differently if some other parametric shape of the interactions were
chosen. This is merely a reflection of the fact that there is no
well-defined concept of test for proportionality; as in all cases of
interaction with at least one quantitative variable involved the test
for interaction is always a test versus some pre-specified alternative
in the form of a specific \emph{shape} of the interaction.

\subsection{The interaction models (non-proportionality)}

It is bad practice to make interaction tests without showing how the
interactions look; however this is not a trivial task with three
different interactions, but if you do not bother to show the shape and
size of estimated interactions, then you should refrain from
interaction tests in the first place.

So we include the identified interactions in the models for the rates.
Note that we also for the sake of notational convenience also include
a void update of the model for mortality after relapse where we take
time since relapse into account:
<<>>=
 i.ar <- update(  c.ar, . ~ . + log(tfd+0.5):size 
                              + log(tfd+0.5):pr.tr
                              + log(tfd+0.5):hormon )
 i.ad <- c.ad
 i.rd <- update(  c.rd, . ~ . + log(tfd+0.5):pr.tr )
ix.rd <- update( xx.rd, . ~ . + log(tfd+0.5):pr.tr )
round( ci.lin( i.ad ), 4 )
round( ci.lin( i.ar ), 4 )
round( ci.lin( i.rd ), 4 )
round( ci.lin( cx.rd ), 4 )
@ %
Note that we have one aliased parameter (\texttt{NA} for \texttt{z}
and \texttt{P}) in the model with effects of the two timescales
(\texttt{tfd}, \texttt{tfr}) and their difference. This is because the
natural spline parametrization include the linear effects of the
variables modeled.

In the following we shall use reference values for each of the
covariates, and show mortality rates as function of time since
diagnosis for select values of the interaction variables:

For each of the three covariates with interactions we construct a
prediction frame with varying levels of the interaction variables:
<<>>=
nd.size <- data.frame( tfd = rep( c(NA,seq(0,15,0.1)), 3 ),
              lex.dur = 100,
                  age = 45,
                 size = rep( levels(Lbc$size), each=152 ),
                nodes = 5,
                pr.tr = 3,
               hormon = levels(Lbc$hormon)[1] )
nd.pr <- data.frame( tfd = rep( c(NA,seq(0,15,0.1)), 6 ),
              lex.dur = 100,
                  age = 45,
                 size = levels(Lbc$size)[2],
                nodes = 5,
                pr.tr = rep( 0:5, each=152 ),
               hormon = levels(Lbc$hormon)[1] )
nd.hormon <- data.frame( tfd = rep( c(NA,seq(0,15,0.1)), 2 ),
              lex.dur = 100,
                  age = 45,
                 size = levels(Lbc$size)[2],
                nodes = 5,
                pr.tr = 3,
               hormon = rep( levels(Lbc$hormon), each=152 ) )
@ % $
For each of these prediction frames we can plot the three estimated
transition rates as we did for the overall rates (or rather the rates
estimated using only the \texttt{tfd} variable as covariate). Moreover
we will plot the estimated rates both from the interaction models
(\texttt{i.}) and the main-effects models (\texttt{c.}):
<<int-size,fig=TRUE>>=
clr <- rainbow(3) ; yl <- c(0.03,60)
ad.c.rate <- ci.pred( c.ad, nd.size ) ; ad.i.rate <- ci.pred( i.ad, nd.size )
ar.c.rate <- ci.pred( c.ar, nd.size ) ; ar.i.rate <- ci.pred( i.ar, nd.size )
rd.c.rate <- ci.pred( c.rd, nd.size ) ; rd.i.rate <- ci.pred( i.rd, nd.size )
matplot( nd.size$tfd, cbind( ad.c.rate, ad.i.rate,
                             ar.c.rate, ar.i.rate,
                             rd.c.rate, rd.i.rate ),
         type="l", lty=rep(c("22","solid"),each=3), 
                   lwd=c(2,0,0), 
                   col=rep(clr,each=6), las=1, lend="butt",
         log="y", xlab="Time since diagnosis (years)",
         ylim=yl, ylab="Rate per 100 PY" )
text( par("usr")[2]*0.95, (10^par("usr"))[3]*1.4^(1:3), 
      c("A->D","A->R","R->D"), col=clr, adj=1, font=2 )
@ % $
\insfig{int-size}{0.8}{Transition rates as function of time since
  diagnosis; the broken lines are from the main effects models and the
  full lines from the interaction model with
  \texttt{\tt age}=54, 
  \texttt{\tt nodes}=5, 
  \texttt{\tt pr.tr}=3, 
  \texttt{\tt hormon}=no and where 
  \texttt{\tt size} assumes the values $<20$ mm, 20--50 mm and $>50$
  mm (only the Alive$\rightarrow$Rel transition).
  Thus the test of interaction is the comparison of the sets of
  parallel broken lines with the non-parallel full lines.} 

<<int-pr,fig=TRUE>>=
ad.c.rate <- ci.pred( c.ad, nd.pr ) ; ad.i.rate <- ci.pred( i.ad, nd.pr )
ar.c.rate <- ci.pred( c.ar, nd.pr ) ; ar.i.rate <- ci.pred( i.ar, nd.pr )
rd.c.rate <- ci.pred( c.rd, nd.pr ) ; rd.i.rate <- ci.pred( i.rd, nd.pr )
matplot( nd.pr$tfd, cbind( ad.c.rate, ad.i.rate,
                           ar.c.rate, ar.i.rate,
                           rd.c.rate, rd.i.rate ),
         type="l", lty=rep(c("22","solid"),each=3), 
                   lwd=c(2,0,0), 
                   col=rep(clr,each=6), las=1, lend="butt",
         log="y", xlab="Time since diagnosis (years)",
         ylim=yl, ylab="Rate per 100 PY" )
text( par("usr")[2]*0.95, (10^par("usr"))[3]*1.4^(1:3), 
      c("A->D","A->R","R->D"), col=clr, adj=1, font=2 )
@ % $
\insfig{int-pr}{0.8}{Transition rates as function of time since
  diagnosis, the broken lines are from the main effects models and the
  full lines from the interaction model with
  \texttt{\tt age}=54, 
  \texttt{\tt size}=20--50 mm, 
  \texttt{\tt nodes}=5, 
  \texttt{\tt hormon}=no and where 
  \texttt{\tt pr.tr} assumes the values 0--6.
  Thus the test of interaction is the comparison of the sets of
  parallel broken lines with the non-parallel full lines --- no
  interaction for the Alive$\rightarrow$Dead transition.}

<<int-hormon,fig=TRUE>>=
ad.c.rate <- ci.pred( c.ad, nd.hormon ) ; ad.i.rate <- ci.pred( i.ad, nd.hormon )
ar.c.rate <- ci.pred( c.ar, nd.hormon ) ; ar.i.rate <- ci.pred( i.ar, nd.hormon )
rd.c.rate <- ci.pred( c.rd, nd.hormon ) ; rd.i.rate <- ci.pred( i.rd, nd.hormon )
matplot( nd.hormon$tfd, cbind( ad.c.rate, ad.i.rate,
                               ar.c.rate, ar.i.rate,
                               rd.c.rate, rd.i.rate ),
         type="l", lty=rep(c("22","solid"),each=3), 
                   lwd=c(2,0,0), 
                   col=rep(clr,each=6), las=1, lend="butt",
         log="y", xlab="Time since diagnosis (years)",
         ylim=yl, ylab="Rate per 100 PY" )
text( par("usr")[2]*0.95, (10^par("usr"))[3]*1.4^(1:3), 
      c("A->D","A->R","R->D"), col=clr, adj=1, font=2 )
@ % $
\insfig{int-hormon}{0.8}{Transition rates as function of time since
  diagnosis, the broken lines are from the main effects models and the
  full lines from the interaction model with
  \texttt{\tt age}=54, 
  \texttt{\tt size}=20--50 mm, 
  \texttt{\tt nodes}=5, 
  \texttt{\tt pr.tr}=3 and where
  \texttt{\tt hormon} assumes the values no and yes. Thus the test of
  interaction is the comparison of the sets of parallel broken lines
  with the non-parallel full lines.}

The general picture from the figures
\ref{fig:int-size},
\ref{fig:int-pr} and
\ref{fig:int-hormon} is the the major interactions are with the relapse
rates, where it seems that the interactions mainly reveal that the major
effects are early, and are possibly even reversed later. If
exploration of interactions were a major concern we might have used 

\section{Predicting state occupancy}

As done in the SiM paper \cite{Crowther.2017} we predict state
occupancy for a patient aged 54, with a transformed progesterone level
of 3, and no hormone therapy (?), for different tumour groups and node
numbers 0, 10 and 20. We shall also compute the expected time alive,
so the calculations will be made for node numbers 0, 5, 10, 15 and 20
for this purpose.

\subsection{Initial cohort}

To this end we construct a Lexis object from \texttt{Rbc}; the main
thing here is to maintain the \texttt{Lexis}-specific attributes which
will be used in the simulation process. And all the time scale
variables too, even if \texttt{A} and \texttt{P} will not be used in
the simulation (because they are not in any of the models) --- the
latter is a feature (or bug) in \texttt{simLexis}; the function will
refer to all timescales in the object even if they are not in the
models and hence not explicitly used in the calculations:
<<>>=
names( Rbc )
Lini <- Rbc[NULL,c("tfd","A","P","tfr",
                   "lex.Cst","lex.Xst","lex.dur","lex.id",
                   "age","size","nodes","pr.tr","hormon")]
pr.nodes <- seq(0,20,5)
npr <- nlevels(Rbc$size) * length(pr.nodes)
Lini[1:npr,"tfd"] <- 0
Lini[1:npr,"tfr"] <- NA
Lini[1:npr,"lex.Cst"] <- "Alive"
Lini[1:npr,"age"] <- 54
Lini[1:npr,"size"] <- rep( levels(Rbc$size), length(pr.nodes) )
Lini[1:npr,"nodes"] <- rep( pr.nodes, each=nlevels(Rbc$size) )
Lini[1:npr,"pr.tr"] <- 3
Lini[1:npr,"hormon"] <- "no"
Lini
str( Lini )
@ %

\subsection{Transition rates}

In order to simulate a number of persons initiating follow-up
(=diagnosed with breast cancer) with these covariate patterns
according to our model, we must also define the transition objects
(that is, specify models for the three transition rates) --- we make
one designed to mimic the models used in the SiM paper
\cite{Crowther.2017} and one using the better fitting model for death
after relapse:
<<>>=
TR  <- list( Alive = list( Dead = i.ad,
                            Rel = i.ar ),
               Rel = list( "Dead(Rel)" = i.rd ) )
TRx <- list( Alive = list( Dead = i.ad,
                            Rel = i.ar ), 
               Rel = list( "Dead(Rel)" = ix.rd ) )
lapply( TR, names )
lapply( TR, lapply, class )
@ %

\subsection{Simulation of a cohort}

With this in place we can simulate:
<<simL,eval=FALSE>>=
sL  <- simLexis( Tr=TR , init=Lini, N=2000, t.range=16 )
sLx <- simLexis( Tr=TRx, init=Lini, N=2000, t.range=16 )
save( sL, sLx, file="sL.Rda" )
@ % $
We asked for simulation of 2000 persons with each of the 15 covariates
patterns in \texttt{Lini}, a total of 30,000 persons:
<<loadsim>>=
load( file="sL.Rda" )
summary( sLx )
@ %

\subsection{State occupancy probabilities}

We can now devise the state probabilities by using \texttt{nState} and
\texttt{pState} --- here we just use an arbitrary subset to get the
object structure:
<<nStstr>>=
nn <- nState( sLx[1:1000,], at=seq(0,16,0.1), from=0, time.scale="tfd" )
pp <- pState( nn, perm=c(1,2,4,3) )
str( pp )
@ %
However this is not what we want; we want the calculation for the 15
different combinations of node and size; so we devise these levels too:
<<defNArr>>=  
( tt <- with( sLx, table( nodes, size ) ) )
prX <- prA <- NArray( c( dimnames( tt ), dimnames( pp ) ) )
str( prA )
@ %
So now we have two arrays to hold the state occupancy probabilities
for all combinations of nodes, size and time from diagnosis; thus we
need a loop over the 15 subsets to devise the relevant probabilities
and put them in the arrays:
<<fillNArr,eval=FALSE>>=
for( nn in dimnames(prA)[[1]] )
for( ss in dimnames(prA)[[2]] )
   {
prA[nn,ss,,] <- pState( nState( subset( sL , nodes==as.numeric(nn) &
                                              size==ss ),
                                at = seq(0,16,0.1), 
                              from = 0, 
                        time.scale = "tfd" ),
                        perm = c(1,2,4,3) )
prX[nn,ss,,] <- pState( nState( subset( sLx, nodes==as.numeric(nn) &
                                              size==ss ),
                                at = seq(0,16,0.1), 
                              from = 0, 
                        time.scale = "tfd" ),
                        perm = c(1,2,4,3) )
   }
save( prA, prX, file="pr.Rda" )
@ %
With this array of probabilities we can now plot the state occupancy
probabilities as a function of time:
<<states,fig=TRUE,width=7,height=8.5>>=
load( file="pr.Rda" )
clr <- col2rgb( c("forestgreen","maroon") )
clr <- cbind( clr, clr[,2:1]*0.6 + matrix(255,3,2)*0.4 )
clr <- rgb( t(clr), max=255 ) 
par( mfrow=c(3,3), mar=c(1,1.5,1,1), mgp=c(3,1,0)/1.6, oma=c(2,2,2,2) )
nnn <- dimnames(prA)[[1]]
sss <- dimnames(prA)[[2]]
for( nn in nnn[c(1,3,5)] ) # only nodes as in the SiM paper
for( ss in sss )
   {
   plot.pState( prX[nn,ss,,], col=clr, xlim=c(0,15), ylab="", xlab="" )
   lines( as.numeric(dimnames(prX)[[3]]), prX[nn,ss,,  2], lwd=3, lty=1, col="black" ) 
matlines( as.numeric(dimnames(prA)[[3]]), prA[nn,ss,,1:3], lwd=1, lty=1, col="white" ) 
   axis( side=2, at=0:10/10, labels=NA, tcl=-0.4 )
   axis( side=4, at=0:10/10, labels=NA, tcl=-0.4 )
   axis( side=2, at=0:50/50, labels=NA, tcl=-0.2 )
   axis( side=4, at=0:50/50, labels=NA, tcl=-0.2 )
   }
mtext( paste( "Size"  ,sss), side=3, at=c(1,3,5)/6, outer=TRUE, line=0, cex=0.66, las=0 )
mtext( paste( "Nodes=",nnn[c(1,3,5)]), side=4, at=c(5,3,1)/6, outer=TRUE, line=0, cex=0.66, las=0 )
mtext( "Time since diagnosis (years)", side=1,      outer=TRUE, line=1, cex=0.66, las=0 )
mtext( "Probability", side=2, outer=TRUE, line=1, cex=0.66, las=0 )
@ %
\insfig{states}{1.0}{Probabilities of being alive without relapse
  (green), with relapse (purple), dead after relapse (light purple),
  and dead without relapse (light green) The black line is the
  estimated survival curve. Computed from the model with effects of
  time since diagnosis as well as since relapse. The white lines
  indicates what would have been obtained with the model with only
  time since diagnosis, that is plots corresponding to those in the
  SiM paper \cite{Crowther.2017}.}

\clearpage

\section{Years lived with and without relapse}

We have the estimated probabilities from the simulation in the arrays
\texttt{prA}, respectively \texttt{prX}. If we want to compute the
years lived during the first 15 years, we want the integral under the
curves. To this end we need a function that does the
triangulation of the area. Here we compute the area under the
curves up til 15 years past diagnosis; first based on the naive
models, then on the models taking time since relapse into account: 
<<>>=
cA <- apply( prA[,,1:151,1:3], c(1,2,4),
             function(x) (sum(x[-1])+sum(x[-length(x)]))/2 * 1/10 )
cA[,,3] <- cA[,,2] - cA[,,1]
dimnames( cA )[[3]] <- c("noRel","Total","Rel")
cA <- cA[,,c(1,3,2)]
round( ftable( cA, row.vars=c(3,2) ), 2 )

cX <- apply( prX[,,1:151,1:3], c(1,2,4),
             function(x) (sum(x[-1])+sum(x[-length(x)]))/20 )
cX[,,3] <- cX[,,2] - cX[,,1]
dimnames( cX )[[3]] <- c("noRel","Total","Rel")
cX <- cX[,,c(1,3,2)]
round( ftable( cX, row.vars=c(3,2) ), 2 )
@ %
Thus it is clear that both the number of nodes and the tumour size
influences the expected lifetime during the first 15 years, although
they primarily influence the relapse-free years lived; the years lived
with relapse is not that much affected.

Note that if we had a simulation-based sample of the probabilities as
outlines above, we would be able to put confidence limits on the
entries in this table as well.

The numbers in the tables above correspond to points at 15 years on
the curves of ``length of stay'' in the SiM paper, so we could have
generated these curves by using the cumulative sums instead, and the
differences and ratios would then have been operations inside the
resulting arrays.

Again, confidence intervals would be easiest to compute by using
simulated datasets from many bootstrap samples, which are not
implemented yet.

\section{Metastases}

A further state, ``metastases'' is recorded too. We included these
among the relapses --- relapse without metastases is at time
\texttt{tor}, whereas metastases is at \texttt{tom}, regardless of
previous relapse.

If we are willing to dispense with subdividing the deaths by the state
from which they occurred we can split the original follow-up (in the
\texttt{Lexis} object \texttt{Lbc}) in one go, using the
\texttt{mcutLexis} function. Note that this requires that relapse
dates recorded as equal to the metastasis dates be coded as
\texttt{NA} thus treating relapse and metastasis as separate events
(that can not occur at the same time). This is what we did when
grooming the data initially, so we can cut the original \texttt{Lexis}
object:
<<>>=
mbc <- mcutLexis( Lbc,
            timescale = "tfd", 
                   wh = c("tor","tom"),
     precursor.states = "Alive",
           new.states = c("Rel","Met"),
           seq.states = TRUE,
           new.scales = c("tfr","tfm") )
summary( mbc, timeScale = TRUE )
mbc <- Relevel( mbc, list( 1, 3, Met=4:5, 2 ) )
summary( mbc )
subset( mbc, lex.id %in% (1328+0:2) )[,1:10]
@ %
The lack of subdivision of deaths by state immediately preceding death
can of course be remedied ``by hand'':
<<>>=
xbc <- transform( mbc, lex.Xst = factor( ifelse( lex.Xst=="Dead" & 
                                                 lex.Cst!="Alive",
                                                 paste( "D(",lex.Cst,")",sep=""),
                                                 as.character(lex.Xst) ) ) )
xbc <- Relevel( xbc )
levels( xbc )
xbc <- Relevel( xbc, c(1:3,5,4) )
levels( xbc )
summary( xbc )
@ % 
<<box-rmx,fig=TRUE,width=8>>=
boxes( xbc, boxpos=list(x=c(15,40,15,85,85,85),
                        y=c(85,50,15,15,85,50)),
            show.BE=TRUE, scale.R=100, wmult=1.1 )
@ % 
\insfig{box-rmx}{1.0}{Transitions when metastases are taken into account.}

We could now model all 6 transitions, exploring the possible effects
of time since entry to the relapse and metastasis states as well as
possible interactions. We might even model mortality rates from
relapse and metastasis with some common parameters.

Eventually we would have specified some model for each of the
transitions, and we could repeat the exercise from above, simulating
state occupancies and time spent in different states.

So far this is left as an exercise to the reader\ldots

\renewcommand{\bibname}{References}
\bibliographystyle{plain}
\begin{thebibliography}{1}
\bibitem{Crowther.2017}
M.~J. Crowther and P.~C. Lambert.
\newblock {{P}arametric multistate survival models: {F}lexible modelling
  allowing transition-specific distributions with application to estimating
  clinically useful measures of effect differences}.
\newblock {\em Stat Med}, 36(29):4719--4742, Dec 2017.
\end{thebibliography}
\addcontentsline{toc}{section}{\bibname}

\end{document}

\clearpage

\section{What is still missing}
 
The arrays \texttt{prA} and \texttt{prX} contain the probabilities of
being in each of the four states (well, cumulated over states) as a
function of time. Additionally, there are two more dimensions to the
arrays corresponding to $5\times 3$ combinations of two covariates
(nodes and size) whereas other covariates (age, progesterone and
hormone therapy) are fixed.

If we wanted some sort of uncertainty associated with the estimates
we cold either simulate using repeat samples from the ``posterior''
distribution of the model parameters, or we could do a bootstrap of
the original sample, re-estimating the models.

In terms of the simulated cohort, we would instead end up with, say
1000 cohorts, each of 100 people, and a corresponding extra dimension
of 1000 on the arrays of probabilities. The could then be used for
computation of confidence intervals for \emph{any} type of measure we
were to derive from the simulated cohorts.

Essentially measures of uncertainty would be referring to quantiles of
the simulated probabilities (well, empirical fractions) from each of
the samples of say 100, persons. Since each sample is devised to
represent a probability we should take the sampling uncertainty into
account when devising probabilities --- that is not just use the
empirical fractions but replace them by a sample from the posterior
distribution of the probability given the empirical fraction.

If we use a flat prior for the probability, the posterior distribution
of the probability given an observed fraction of $x/n$ is Beta with shape
$(x+1,n-x+1)$. Thus a simple deterministic jitter of the array of
probabilities applied before computing the confidence limits. However,
this does not take the time-dependence of the probabilities into account.

To be continued \ldots

\subsection{Technical note on \texttt{simLexis} implementation}

The transition objects are large and clumsy, and may even contain the same
models more than once. It would be better to only have the contents as
the \emph{names} of the transition models, and inside use \texttt{get}
to construct the objects currently used:
<<eval=FALSE>>=
Tr <- lapply( Tr, lapply, get )
@ %
This will also make it easier to use bootstrapped data for evaluation
of uncertainty. For a given bootstrap sample of data we would make updated model
objects with names appended with some string, so that the input for
each cycle of the simulation loop over bootstrap samples of data would
be using an input transition object of the form:
<<eval=FALSE>>=
bootTr <- lapply( Tr, lapply, function(x) paste("BOOT",x,sep="") )
@ %
Generation of the model objects with these names would be using only
the unique elements, avoiding fitting the same model more than once:
<<eval=FALSE>>=
unique.models <- unique( unlist( Tr ) )
for( m in unique.models )
   {
   assign( paste("BOOT",m,sep=""),
           update( get(m), data=boot.Lexis(data) ) )
   }
@ %

\end{document}
