mgus.Rnw

<<echo=FALSE>>=
library("survival", quietly=TRUE)
opts_chunk$set(comment=NA, tidy=FALSE, highlight=FALSE, echo=FALSE,
               fig.with=7, fig.height=5.5, fig.path="figures/",
               out.width="\\textwidth", out.height="!", device="pdf",
               cache=FALSE, background="#ffffff",
               warning=FALSE, error=FALSE, prompt=TRUE)
options(contrasts= c("contr.treatment", "contr.poly"),
        show.signif.starts = FALSE, continue=" ", width=60)
par(mar=c(4.1, 4.1, 1.1, 1.1))

crisk <- function(what, horizontal = TRUE, ...) {
    nstate <- length(what)
    connect <- matrix(0, nstate, nstate,
                      dimnames=list(what, what))
    connect[1,-1] <- 1  # an arrow from state 1 to each of the others
    if (horizontal) statefig(c(1, nstate-1),  connect, ...)
    else statefig(matrix(c(1, nstate-1), ncol=1), connect, ...)
}
state3 <- function(what, horizontal=TRUE, ...) {
    if (length(what) != 3) stop("Should be 3 states")
    connect <- matrix(c(0,0,0, 1,0,0, 1,1,0), 3,3,
                      dimnames=list(what, what))
    if (horizontal) statefig(1:2, connect, ...)
    else statefig(matrix(1:2, ncol=1), connect, ...)
}
@ 

\section{Competing risk}
\begin{frame}
{\Large 3. Competing risks}
\end{frame}

\begin{frame}
<<fourfig>>=
oldpar <- par(mfrow=c(2,2), mar=c(.1, .1, .1, .1))
crisk(c("Alive", "Dead"))
crisk(c("A", "D1", "D2", "D3"))

what <- c("Health", "Illness", "Death")
connect <- matrix(0, 3,3, dimnames=list(what, what))
connect[,3] <- 1
statefig(c(1,2), connect)
delta <- c(-.02, .02)
arrows(c(.35, .65), c(.55, .70) + delta, c(.65, .35), c(.7, .55) + delta,
       angle=20, length= .1)

what <- c("0", "1", "2", "...")
connect <- matrix(0, 4,4, dimnames=list(what, what))
connect[1,2] <- connect[2,3] <- connect[3,4] <- 1
statefig(c(1,1,1,1), connect)
par(oldpar)
@ 
\end{frame}

\begin{frame}{Monoclonal Gammopathy of Undetermined Significance (MGUS)}
  \begin{itemize}
    \item Subjects with a dominant clone in their plasma cell population,
      but without malignancy ($\ge 2$\% of plasma cells in the clone).  
    \item Normally found incidentally to other tests.
    \item Should the patient be worried?
    \item About 1\% per year convert to overt malignancy.
    \item Essentially independent of age and sex.
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{Progression of MGUS}
  \begin{itemize}
    \item 1384 subjects with monoclonal gammopathy of undetermined significance
      (MGUS)
    \item R. Kyle, New Engl J Med 346:564-569 (2002)
    \item Questions
      \begin{itemize}
        \item Pattern of death and progression
        \item Relationship to age, sex, hemoglobin, creatinine,
          and amount of protein in the ``spike''
      \end{itemize}
  \end{itemize}
\end{frame}

\begin{frame}[fragile]
<<kyle1, echo=TRUE>>=
mgus2[1:5,]
@

\begin{itemize}
  \item Classic data set
    \begin{itemize}
      \item One row per subject
      \item Separate columns for death and plasma cell malignancy
      \item Follow-up continues after PCM
    \end{itemize}
\end{itemize}
\end{frame}

\begin{frame}
<<mgus1>>= 
oldpar <- par(mfrow=c(1,2))
hist(mgus2$age, nclass=30, main='', xlab="Age at diagnosis")
#with(mgus2, tapply(age, sex, mean))

mfit1 <- survfit(Surv(futime, death) ~ sex, data=mgus2)
plot(mfit1, col=c(1,2), xscale=12, mark.time=FALSE, lwd=2,
     xlab="Years post diagnosis", ylab="Survival")
legend(200, .8, c("female", "male"), col=1:2, lwd=2, bty='n')
par(oldpar)
@
\end{frame}


\begin{frame} {Common analysis}
<<kyle2, echo=FALSE>>=
etime <- with(mgus2, pmin(futime, ptime))
event <- with(mgus2, ifelse(pstat==0, 2*death, 1))
event <- factor(event, 0:2, c("censor", "PCM", "death"))
kdeath <- survfit(Surv(etime, event=='death') ~ sex, data=mgus2)
kprog  <- survfit(Surv(etime, event=='PCM') ~ sex, data=mgus2)
plot(kdeath, col=1:2, mark.time=FALSE, xscale=12, 
     xlab="Years")
lines(kprog, col=1:2, mark.time=FALSE, xscale=12, lty=2)
legend(132, .7, c("F", "M"), col=1:2, lty=1, bty='n')
text(240, .92, "PCM, censoring at death")
text(200, .36, "Death, censoring at PCM", adj=0)
@ 
\end{frame}

\begin{frame}
<<kyle2b>>=
plot(kdeath, col=1:2, mark.time=FALSE, xscale=12, 
     xlab="Years")
lines(kprog, col=1:2, mark.time=FALSE, xscale=12, lty=2)
legend(132, .7, c("F", "M"), col=1:2, lty=1, bty='n')
text(240, .92, "PCM, censoring at death")
lines(mfit1, col=1:2, lwd=2)
@ 
\end{frame}

\begin{frame}{Competing Risks}
  \begin{itemize}
    \item For this illustration we are only interested in the first 
      event for each subject.
    \item Formally we are treating progression to a plasma cell malignancy (PCM)
      as an \emph{absorbing state}, i.e., one
      that subjects never exit.

    \item The event variable was created as a factor.
      The first level of the factor must be censoring, which is the status code
      for those whose follow-up terminated without reaching either endpoint.
      Codes for the remaining states can be in any order. The labels are
      arbitrary.
    \item A simple print of the \code{mfit2} object shows the order in
      which the curves will be displayed.  This information was used to
      choose the line types and colors for the curves.
    \item Curves start at 0.
  \end{itemize}
\end{frame}


\begin{frame}{Competing Risk (Aalen-Johansen)}
<<mgus2, echo=FALSE>>=
etime <- with(mgus2, ifelse(pstat==0, futime, ptime))
event <- with(mgus2, ifelse(pstat==0, 2*death, 1))
event <- factor(event, 0:2, labels=c("censor", "pcm", "death"))

mfit2 <- survfit(Surv(etime, event) ~ sex, data=mgus2)

oldpar <- par(mar=c(5.1, 4.1, 1,.1))
layout(matrix(1:2,1,2), widths=2:1)
plot(mfit2, col=c(1,2,1,2), lty=c(1,1,2,2),
     mark.time=FALSE, lwd=2,  xscale=12,
     xlab="Years post diagnosis", ylab="Probability in State")
legend(170, .5, c("death w/o PCM:f", "death w/o PCM:m", 
                  "PCM:female", "PCM:male"), 
       col=c(1,2,1,2), lty=c(2,2,1,1), lwd=2, bty='n')

crisk(c("MGUS", "PCM", "Death"))
par(oldpar)
layout(1)
@ 
\end{frame}

\begin{frame}[fragile]
<<mgus2b, echo=TRUE>>=
mfit2 <- survfit(Surv(etime, event) ~ sex, data=mgus2)
print(mfit2, rmean=240, scale=12)
@ 
\end{frame}

\begin{frame}{Footnote: Alternate model}
<<cralt>>=
par(mfrow=c(1,2), mar=c(5.1, 4.1, 1,1))
crisk(c("MGUS", "PCM", "Death"))
state3(c("MGUS", "PCM", "Death"))
@ 
\end{frame}

\begin{frame}[fragile]
  \begin{itemize}
    \item Left figure: ever PCM and death without PCM
    \item Right figure: currently in PCM, ever dead
    \item Same status variable, different data set
  \end{itemize}

Raw data
<<>>=
mgus2[80:83, c("id", "age", "sex", "ptime", "pstat", "futime","death")]
@ 
\end{frame}

\begin{frame}[fragile]
CR data: add two variables
<<>>=
m2 <- mgus2
m2$etime <-etime; m2$event <- event
m2[80:83, c("id", "age", "sex", "ptime", "pstat", "futime","death", "etime",
            "event")]
@ 

Figure 2: add lines
<<>>=
ptemp <- with(mgus2, ifelse(ptime==futime & pstat==1, ptime-.1, ptime))
mgus3 <- tmerge(mgus2, mgus2,  id=id, death=event(futime, death),
                  pcm = event(ptemp, pstat))
mgus3 <- tmerge(mgus3, mgus3, id, enum=cumtdc(tstart))
mgus3$event <- factor(mgus3$pcm + 2*mgus3$death, 0:2, 
                      c("censor", "PCM", "death"))
subset(mgus3, id %in% 80:83, c(id, age, sex, tstart, tstop, event))
@ 
\end{frame}

\begin{frame}
  \begin{itemize}
    \item The first version of the data set generated errors.
    \item 9 subjects have PCM and death declared at the same time.
    \item Treated as PCM in the CR analysis
    \item For the multi-state model we need to be explicit
    \item Push progressions back by .1 month when there is a tie
      \pause
    \item Many data sets have case like this.
  \end{itemize}
\end{frame}

\begin{frame}
<<mguspair>>=
par(mfrow=c(1,2), mar=c(5.1, 4.1, 1,1))
plot(mfit2, xscale=12, col=1:2, lwd=2, lty=c(1,1,2,2))
pfit3 <- survfit(Surv(tstart, tstop, event) ~ sex, mgus3, id=id)
plot(pfit3, xscale=12, col=1:2, lwd=2, lty=c(1,1,2,2))
@ 
\end{frame}

\begin{frame}[fragile]
<<kyle10, echo=FALSE>>=
plot(pfit3[,1], mark.time=FALSE, lwd=2, col=1:2,
     xscale=30.5, xlab="Years post diagnosis", ylab="Alive with PCM")
@ 
\end{frame}

\begin{frame}[fragile]{Look at multiple factors}
<<kyle4, echo=TRUE, fig.show="hide">>=
oldage <- (mgus2$age > 70)
bigm   <- mgus2$mspike > 1
pfit2 <- survfit(Surv(etime, event) ~ oldage + sex + 
                 bigm, data= mgus2)

plot(pfit2, col=1:8, lty=rep(1:2, each=8),
     xscale=12, xlab="Years since MGUS")
@ 
\end{frame}


\begin{frame}{Tangle of yarn plot}
  \myfig{kyle4-1}
\end{frame}


\begin{frame}{Hazard models}
  \begin{itemize}
     \item Model the individual rates $\lambda_{jk}(t)$ from state $j$
       to state $k$
     \item Separate model for each transition
	\begin{itemize}
	  \item which covariates for each
	  \item which, if any, coefficients are shared
	  \item what time scale for each transition (baseline hazard)
	  \item do any transitions share a baseline
	\end{itemize}
      \item A given arrow only depends on the starting box and the transitions
        \begin{itemize}
          \item At risk = in the starting box
          \item Event = transition of \emph{this} type, all others are treated
            as censored
        \end{itemize}
  \end{itemize}
\end{frame}

\begin{frame}[fragile]
<<cmugus1, echo=TRUE>>=
cfit1 <- coxph(Surv(etime, event=="death") ~ sex + age +mspike,
               data= mgus2)
round(summary(cfit1, scale=c(1, 10, 1))$coef, 2)

cfit2 <- coxph(Surv(etime, event=="pcm") ~ sex + age + mspike,
               data= mgus2)
round(summary(cfit2, scale=c(1, 10, 1))$coef, 2)

quantile(mgus2$mspike, na.rm=TRUE)
@ 
\end{frame}

\begin{frame}[fragile]
  \frametitle{Simple event rates}
<<mgu2sb, echo=TRUE>>=
pfit <- pyears(Surv(ptime, pstat) ~ sex, scale=12*100, 
                 data=mgus2)
pfit$event/pfit$pyears
@ 
\begin{itemize}
  \item Overall rate is 1\% per year for males, 1.1\% for females
  \item In 15 years we would expect 15--16\% to progress \\
    \emph{if there were no deaths}
\end{itemize}
\end{frame}


\begin{frame}{Progression to PCM}
<<mgus2d>>=
plot(mfit2[,1],  col=c(1,2), xmax=15*12, 
     mark.time=FALSE, lwd=2,  xscale=12,
     xlab="Years post diagnosis", ylab="Progression to PCM")
legend(96, .03, c("Female", "Male"), lty=1, col=1:2, lwd=2,
       bty='n')
@ 
\end{frame}

<<expect, echo=FALSE>>=
# mean years of life remaining, used in lines below, but not shown
efit1 <- survexp(~ sex, data=mgus2, times=365.25*(0:150)/2,
                 rmap=list(sex=sex, age=age*365.25, year=as.Date("1980-01-01")),
                 ratetable=survexp.mn)
lifeleft <- colSums(efit1$surv)/2  # area under the curve
@ 

\begin{frame}{Lesson 1}
  \begin{itemize}
    \item Any given rate (arrow) is modeled using that endpoint and
      ignoring all others (treat them as censored).  Individual rates
      are local.
    \item The probability of being in any one state (box) depends on \emph{all}
      the rates.  States are global.
      \begin{itemize}
        \item Sex has no effect on progression to PCM.
        \item But females have a higher 15 year risk (9.5 vs 7.5) \pause
        \item Mean life remaining for this age distribution is 
          \Sexpr{round(lifeleft[1],1)} and \Sexpr{round(lifeleft[2],1)}
          years, respectively.
      \end{itemize}
  \end{itemize}
\end{frame}

\begin{frame}
  \begin{itemize}
    \item Underlying biology question: does gender affect the rate of
      progression?
    \item Public health question: does gender affect the lifetime risk?
    \item They have different answers.
  \end{itemize}
\end{frame}

\begin{frame}
<<mgus3>>=
pcmbad <- survfit(Surv(etime, pstat) ~ sex, data=mgus2)
plot(pcmbad, lwd=2, fun="event", conf.int=FALSE,
     xscale=12, xmax=20*12, col=1:2, lty=2,
     xlab="Years after MGUS diagnosis", ylab="Fraction with PCM")
lines(mfit2[,1], lty=1, lwd=2, col=1:2)
legend(0, .18, c("PCM, incorrect curve", "PCM, competing risk"),
       col=1, lwd=2, lty=c(2,1), bty='n')
@ 
\end{frame}


\begin{frame}[fragile]
  \frametitle{Consequences of the rates}
  \begin{itemize}
    \item Can we predict multivariate outcome from a set of Cox models? Yes.
    \begin{itemize}
      \item As with all Cox models, prediction must be
        for some particular person(s).
      \item Decision: male/female by age 60/80, with mspike= 1.2
      \item Create a 4 observation data set \code{newdata}
       \pause
     \item Individual hazard, from each Cox model alone, are building
          blocks
          \begin{itemize}
            \item Aalen-Johansen estimator, with $\lambda(t | z)$ as entries.
              \item Requires \emph{all} the fits at once
              \item Result is 8 curves: predicted lifetime risk of PCM and of
                death before PCM, each for m/f  at age 60/80.
          \end{itemize}
        \item Individual ``predicted survival from Cox model'' curves are
          useless.
    \end{itemize}
  \end{itemize} 
\end{frame}

\begin{frame}
\begin{align*}
  H&= \left( \begin{array}{ccc}
     \lambda_{11}(t) & \lambda_{12}(t) &\lambda_{13}(t)  \\
     \lambda_{21}(t) & \lambda_{22}(t) &\lambda_{23}(t)  \\
     \lambda_{31}(t) & \lambda_{32}(t) &\lambda_{33}(t)  \\
       \end{array} \right) \\
  &= \left( \begin{array}{ccc}
       * & \lambda_{12}(t) &\lambda_{13}(t)  \\
       0 &  *  & 0  \\
       0 &  0 &  *  \\
       \end{array} \right) \\
  \end{align*}
\end{frame}

\begin{frame}
<<mfit4, echo=FALSE>>=
cfit3 <- coxph(Surv(etime, event=="death") ~ sex + age + mspike, mgus2)
cfit4 <- coxph(Surv(etime, event=="pcm") ~ sex + age + mspike, mgus2)

newdata <- expand.grid(sex=c("F", "M"), age=c(60, 80), mspike=1.2)

tstate <- c("entry", "pcm", "death")
temp <- matrix(list(), 3,3,
               dimnames=list(from=tstate, to=tstate))
temp[1,2] <- list(survfit(cfit4, newdata, std.err=FALSE))
temp[1,3] <- list(survfit(cfit3, newdata, std.err=FALSE))
csurv  <- survfit(temp, p0 =c(entry=1, PCM=0, death=0))
plot(csurv[,2], xmax=25*12, xscale=12, 
     xlab="Years after MGUS diagnosis", ylab="PCM",
     col=1:2, lty=c(1,1,2,2), lwd=2)
abline(0, .01/12, col=4, lty=3)
legend(10, .14, outer(c("female", "male   "), 
                     c("diagnosis at age 60", "diagnosis at age 80"), 
                      paste, sep=", "),
       col=1:2, lty=c(1,1,2,2), bty='n', lwd=2)
@ 
\end{frame}

\begin{frame}{Prediction at 25 years}
\begin{center} \begin{tabular}{cccc}
    & female & male& delta \\ \hline
   diagnosis at age 60 & .142 & .118 & .024\\
   diagnosis at age 80 & .086 & .060 & .026 
   \end{tabular} \end{center}

\begin{itemize}
  \item 2.5\% increase for females
  \item Not a constant wrt to other covariates, even though the 
    HR for sex is constant within each each of death and PCM.\\
    PH for components $\ne$ PH for the composite.
  \item A good summary is a population average prediction = mean prediction
    over the all the other covariates.
    \begin{itemize}
      \item For each combination of age and mspike in the data set
      \item Compute the CI curves, tabulate the difference at age 90
      \item Bootstrap
    \end{itemize}
\end{itemize}
\end{frame}


\begin{frame}{mstate package}
  \begin{itemize}
    \item Create a \emph{stacked} data set
      \begin{itemize}
        \item 1384 obs for the MGUS to PCM transition
        \item 1384 obs for the MGUS to death transition
        \item 115 obs for the PCM to death transtion (optional)
        \item Add \code{from} and \code{to} as covariates
        \item Each obs has \code{status} =1 if \emph{this} transition occured
      \end{itemize}
    \item Create a 3x3 transition matrix
    \item Fit all the models at once \\
      \code{coxph(Surv(time, status) ~ (age + sex)* strata(from, to), ...)}
    \item The \code{mfit} command will create the AJ curves
      \begin{itemize}
        \item includes variance/covariace
      \end{itemize}
  \end{itemize}
\end{frame}