class2.Rnw

\documentclass[11pt]{article}
\usepackage{amsmath}
\title{Analysis of the NAFLD data set}
\author{Terry Therneau}

\addtolength{\textwidth}{1in}
\addtolength{\oddsidemargin}{-.5in}
\setlength{\evensidemargin}{\oddsidemargin}

\newcommand{\code}[1]{\texttt{#1}}

\begin{document}
<<echo=FALSE>>=
opts_chunk$set(comment=NA, tidy=FALSE, highlight=FALSE, echo=TRUE,
               fig.with=5, fig.height=4, fig.path="figures/",
               out.width="0.9\\textwidth", out.height="!", device="pdf",
               cache=FALSE, background="#FFFFFF", mypar=TRUE,
               warning=FALSE, error=FALSE)
knit_hooks$set(mypar = function(before, options, envir) {
    # see Xie, Dynamic documents with R and knitr, 10.1.1 and 10.1.4
    if (before) { 
        mar=c(4.1,4.1, .2, .2)  # we will float figures, so no titles
    } else NULL
})

library(survival, quietly=TRUE)
library(splines, quietly=TRUE)
library(mstate,  quietly=TRUE)
library(timereg, quietly=TRUE)
options(show.signif.stars=FALSE)  # display statistical intelligence
load('nafld.rda')

overlap <- with(nafld3, sum(event=="nafld" & days > 0)) #case and control
@ 
\maketitle

\section{NAFLD}
\subsection{Background}
The NAFLD data set is used as an example in the course; this file works
through those examples in more detail and shows alternate computations.
The primary reference for the NAFLD study is Allen \cite{Allen2018}.
The incidence of non-alcoholic fatty liver disease (NAFLD) has been rising
rapidly in the last decade and it is now one of the main drivers of
hepatology practice \cite{Tapper2018}.
It is essentially the presence of excess fat in the
liver, and parallels the ongoing obesity epidemic.
Approximately 20-25\% of NAFLD patients will develop the inflammatory state
of non-alcoholic steatohepatitis (NASH), leading to faster fibrosis 
progression to end-stage disease.
NAFLD can be accurately diagnosed by MRI methods, 
but NASH diagnosis currently requires a biopsy.

\subsection{Data}
Population-based epidemiological research can be
conducted in Olmsted County, Minnesota, because
medical care is effectively self-contained within the
community and there are only a few providers. 
All medical information from each of these providers for each 
individual resident of Olmsted County
is linked to a single identification number allowing
linkage across health care systems. The diagnoses
recorded in billing data from each institution are
indexed so that all individuals with any specific diagnosis
can be identified across systems and the corresponding patient 
records can easily be retrieved from each site.
The system and infrastructure that collects, collates,
and indexes the data is the Rochester Epidemiology Project REP
\cite{Rocca2012, Sauver2012}, which has been
funded by the National Institutes of Health since
1966.

Using REP resources, the current study constructed a community cohort of
all adult NAFLD subjects from 1997 to 2014  along with 4 potential controls
for each case.
Because NAFLD is often a disease of exclusion, a NAFLD diagnosis followed
shortly by the diagnosis of another liver disease is considered a false
positive.  For this study we only consider ``confirmed'' NAFLD, i.e., if
someone were diagnosed on 2001-06-20, the index date for confirmed NAFLD 
would be 2002-06-20, assuming that another liver diagnosis, death, or
incomplete follow-up did not intervene.  Control subjects are matched on
age and sex, and for each their follow-up also commences on the 
``confirmed NAFLD'' date.

After the case-control set was assembled, we removed any subjects with less
than 7 days of follow-up.  These subjects add little information, and it
prevents a particular confusion that can occur with a multi-day medical visit
where two results from the same encounter have different dates.
To protect patient confidentiality all time intervals are in days since
the index date; none of the dates from the original data were retained.
Subject age is their integer age at the index date, and the subject
identifier is an arbitrary integer. 
As a final protection, we include only a 90\% random sample of the data.
As a consequence analyses results will not exactly match the
original paper.

The data consists of 3 R data sets.
\begin{itemize}
  \item \code{nafld1} has one observation for each subject containing
    baseline covariates.
  \item \code{nafld2} has selected post-index covariate values.  
    A given subject can have anywhere from 0 to several hundred rows.
  \item \code{nafld3} has subject outcomes, one per line.  Outcomes such
    as hyperlipidemia may appear on multiple dates, each corresponding to
    a patient visit.
\end{itemize}
The \code{nafld2} data set has been restricted to time points after the
index date, but \code{nafld3} needs events from both before and after to
set a subject's starting state in the multi-state models,
e.g., hypertension before or after NAFLD.

\subsection{Cases and controls}
A recurring and fundamental mistake in survival analysis is immortal time
bias, variations of which occur whenever either covariate values or
patient selection at a particular time are allowed to ``peek into the
future''.
When selecting the age/sex matched controls for a given case at age $a$,
it was important to \emph{not} exclude subjects who will become NAFLD at
some later age $>a$ as potential matches.
(Convincing an investigator of this fact can sometime be very hard.)
In this data set \Sexpr{overlap}
of the NAFLD cases were also selected as controls.
There are at least three ways to analyse them.
\begin{enumerate}
  \item One approach is to enter such a subject into the 
    study twice,
    from the first index until study end, copy 1 of the subject 
    is a control, and from NAFLD to
    study end copy 2 of the subject is a case. 
    This is the best way to approach the usual ``table 1'' of a paper, which
    lists characteristics of all the cases and controls.
    For regression models the approach requires
    using a robust variance since the two copies are obviously not independent.
  \item Treat NAFLD status as a time dependent covariate.  The hazard ratio is
    now the comparison of someone with current NAFLD to others 
    who currently do not have the condition.  
  \item Treat them as a single control. This reduces our sample size of
    NAFLD subjects however.  
    Note that this approach is statistically justified since a 
    ``new'' NAFLD is accepted or rejected based on something in their past,
    i.e., prior use as a control. 
\end{enumerate}
Many of the issues in this choice are philosophical and I haven't 
digested them all.

\subsection{Time dependent covariates and endpoints}
Our primary data set will use NAFLD as a time dependent covariate.
The \code{tmerge} routine allows for easy creation of time dependent
covariates and endpoints.
The routine's first three arguments are always two data sets and a
subject identifier variable, then other arguments that perform actions.
On the first call, \code{data1} contains one row per subject and any
baseline variables which do not change over time, such as enrollment age or sex,
and \code{data2} is used to set the initial time line for the subject.
Subsequent calls build on this data set, inserting new rows as necessary for
time dependent covariates (created by tdc calls) and events (created by event
calls).  
The \code{tmerge} function automatically takes care of an important 
distinction with tied times, namely that if subject Jones has a death on
day 400 and subject Smith converts from control to NAFLD on day 400,
then the time dependent covariate changes \emph{after} that event occurs.
The new NAFLD value of 1 for Smith will apply to deaths that happen after
day 400 but not to the one on day 400.
Short variable names like \code{diab} for diabetes were chosen
below to help printouts fit the page width.

<<data1, cache=FALSE>>=
data1 <- tmerge(nafld1[,1:8], nafld1, id=id, death = event(futime, status)) 
data1 <- tmerge(data1, subset(nafld3, event=="nafld"), id, 
                nafld= tdc(days))
data1 <- tmerge(data1, subset(nafld3, event=="diabetes"), id,
                diab= tdc(days), e1= event(days))
data1 <- tmerge(data1, subset(nafld3, event=="htn"), id, 
                htn= tdc(days), e2= event(days))
data1 <- tmerge(data1, subset(nafld3, event=="dyslipidemia"), id, 
                dyslip = tdc(days), e3= event(days))
data1$age1 <- with(data1, age + tstart/365.25)
data1$age2 <- with(data1, age + tstop/365.25) 
attr(data1, 'tcount')
@ 

<<echo=FALSE>>=
tc <- attr(data1, 'tcount')   # shorter name for use in Sexpr below
icount <- table(table(nafld3$id)) #number with 1, 2, ... intervals
ncount <- sum(nafld3$event=="nafld")
@ 

The \code{tcount} attribute tells us a lot about the creation process.
Each addition of a new endpoint or covariate to the data generates one 
row in the table.  Column labels are explained by figure \ref{fig:timeline}.
\begin{itemize}
  \item There are \Sexpr{tc[1,7]} last fu/death additions, 
    which by definition fall at the
    trailing end of a subject's observation interval: they define the
    interval.
  \item There are \Sexpr{tc[2,2]} nafld splits that fall after the end
    of follow-up (`late').
    These are subjects whose first NAFLD fell within a year of the end of
    their time line, and the one year delay for ``confirmed'' pushed them
    over the end.  (The time value in the \code{nafld3} data set is 1 year
    after the actual notice of NAFLD; no other observations have this
    offset added).  The time dependent covariate \code{nafld} never turns
    from 0 to 1 for these subjects. 
    (Why were these subject not removed earlier by my ``at least 7 days of
    follow-up rule?  They are all controls for someone else and so appear
    in the data at a younger age than their NAFLD date.)
  \item \Sexpr{tc[3,1]} of the diabetes diagnoses are before entry, i.e., 
    these are the prevalent cases.  One diagnosis occurred on the day of
    entry (``leading''), and will not be counted as a post-enrollment endpoint,
    all the other fall somewhere between study entry and last follow-up.
  \item Conversely, \Sexpr{tc[5,7]} subjects were diagnosed with hypertension
    at their final visit (``trailing'').  These will be counted as an
    occurrence of a hypertension event (\code{e2}), but the time dependent 
    covariate \code{htn} will never become 1.  
\end{itemize}

Such a detailed look at data set construction may seem overmuch.   
Our experience is that issues with covariate and event timing
occur in nearly all data sets, large or small.  The 13 NAFLD cases ``after
last follow-up'' were for instance both a surprise and a puzzle to us;
but we have learned through experience that it is best not to proceed until
such puzzles are understood.  (This particular one was benign.)
If, for instance, some condition is noted at autopsy, do we want the related
time dependent covariate to change before or after the death event?
Some sort of decision has to be made, and it is better to look and understand
than to blindly accept some programming default.

\subsection{Simple counts}
<<table1, echo=FALSE>>=
tcount <- function(data) {
    first <- !duplicated(data$id)
    data2 <- subset(data, first, c(age, male, bmi, diab, htn, 
                                   dyslip, death))

    temp1 <- sapply(data2, mean, na.rm=TRUE)
    temp2 <- sapply(data2, function(x) diff(quantile(x, c(.25, .75), na.rm=TRUE)))
    n <- nrow(data2)
    count <- function(x) 100*x
    c(n, temp1[1], temp2[1], count(temp1[2]),
      temp1[3], temp2[3], count(temp1[4:7]))
    }
                
tab1 <- rbind(tcount(data1), tcount(subset(data1, nafld==0)),
              tcount(subset(data1, nafld==1)))
dimnames(tab1) <- list(c("All subjects", "Control", "Case"),
                       c("N", "age(mean)", "age (IQR)", "male %",
                         "bmi(mean)", "bmi (IQR)",
                         "diabetes %", "htn %", "dyslipidemia %", "deaths %"))
round(t(tab1), 1)
@ 
The above table gives counts as of the index date for the control subjects
and as of the NAFLD date for the NAFLD subjects.
Subjects who are both cases and controls are counted in both columns.
NAFLD subjects have nearly a 2x higher rate of metabolic
comorbidities at diagnosis and a somewhat higher mean BMI.
The age and sex distributions are matched by the study design.
(Although four potential controls were selected for each case, about 1/3 have
fewer after subject exclusions for prior liver disease; otherwise the
age and sex distributions would match exactly.)


\section{Time to death}
\subsection{Overall survival}
<<sfig1, cache=FALSE>>=
data1$nafld0 <- with(data1, nafld[match(id, id)])  # nafld status at entry
sfit1 <- survfit(Surv(tstart, tstop, death) ~ male + nafld0, data1)
sfit2 <- survfit(Surv(age1, age2, death) ~ male + nafld0, data1,
                 start.time=20)
@ 

% float the figure, and don't show all the plotting code
<<sfig1b, echo=FALSE, fig.cap='Survival curves for NAFLD vs. control subjects, on age scale and on time since entry scale.'>>= 
oldpar <- par(mfrow=c(1,2), mar=c(5,5,1,1))
plot(sfit1, xscale=365.25, fun="event", lwd=2, col=c(1,1,2,2), lty=c(2,1,2,1),
     xlab="Years from index date", ylab="Deaths")
legend(1, .35, c("Control, female","NAFLD, female", "Control, male",
                 "NAFLD, male"), col=c(1,1,2,2), lty=c(2,1,2,1), lwd=2, bty='n')
plot(sfit2, fun="event", lwd=2, col=c(1,1,2,2), lty=c(2,1,2,1),
     xlab="Age", ylab="Deaths")
par(oldpar)
@   

Figure \ref{fig:sfig1b} shows survival curves for the data on both the age
and time since index scales. 
The curves for female NAFLD subjects and male control subjects almost overlay,
and this will be reprised in the Cox model coefficients further below
for male gender and for NAFLD; they each increase the death rate by about
1.5 fold.
Survival curves with a time dependent covariate are problematic: if you think
at all deeply about them it becomes difficult to say exactly what they
represent.  The curves shown in figure \ref{fig:sfig1b} use the NAFLD status
at first entry for each subject and thus correspond to case 3 of our
earlier list.

\subsection{Initial fits}
Now for the analysis.
A first question is what time scale to use.
There are three obvious choices.
\begin{itemize}
 \item Fits on the ``time since entry'' scale are the obvious approach 
   for clinical trial data, where due to a long list of eligibility requirements
   the subjects are similar at their entry date. \\
   A Cox model compares each subject who has an event to all others who ``could
have had'' the event and who share the same baseline risk.
In the NAFLD data set, however, time since discovery of NAFLD is not nearly
as compelling due to the unknown, and possibly quite long, duration of the
condition prior to detection. 
For the controls the time since ``your identifier was picked out of a hat''
is completely odd. \\
If this scale is used, then age and sex need to also be in
the model, and the model for them needs to be \emph{right}.
The yearly death rate ranges from .03 to 500 per thousand for this population
and age range; any small lack of fit in the age*sex modeling 
can dominate all other covariates.

    \item Use age as the time scale.  This naturally matches subjects on age.
The model can go further and stratify on sex, so that proportional hazards is
sidestepped for that effect as well.  
One downside is that we do not get a coefficient estimate
for age, but that is minor since we already know that age is related
to the outcomes.

\item A third option is the use a matched analysis.  Each case/control set
is a separate stratum with follow-up time used within the stratum.  This also
matches on age and sex, but raises some technical and programming issues if
we want further summaries than the hazard ratio.
\end{itemize}

Start with a simple fit using follow-up time.
The results show that, as expected, the age effect is dominating.
There is a statistically significant non-linear effect of age on the
log hazard;
however, figure \ref{fig:cfit1} shows that the effect is very close to
linear.  
The simple tests for non proportional hazards are not remarkable,
but figure \ref{fig:zp1} shows a possibly interesting early effect for
age and a strong quadratic effect for NAFLD.  
For the latter the simple proportional hazards (PH) test used by 
\code{cox.zph}, namely a linear fit plus a test of non-zero slope,
was clearly insufficient.

<<cfit1, fig.cap="Estimated age effect from the model on fu-time scale">>=
cfit1 <- coxph(Surv(tstart, tstop, death) ~ age + male + nafld, 
               data=data1)
cfit1

zp1 <- cox.zph(cfit1, transform="identity")
zp1

cfit1b <-  coxph(Surv(tstart, tstop, death) ~ pspline(age,4) + male + nafld, 
               data=data1)
cfit1b
termplot(cfit1b, term=1, se=TRUE)
@ 

\begin{figure}
<<zp1b, echo=FALSE>>=
oldpar <- par(mfrow=c(2,2), mar=c(5,5,1,1))
plot(zp1[1], resid=FALSE, xlab="Time since enrollment")
abline(h=coef(cfit1)[1], col=2, lty=2)
plot(zp1[2], resid=FALSE, xlab="Time since enrollment")
abline(h=coef(cfit1)[2], col=2, lty=2)
plot(zp1[3], resid=FALSE, xlab="Time since enrollment")
abline(h=coef(cfit1)[3], col=2, lty=2)
par(oldpar)
@ 
\caption{Estimated $\beta(t)$ functions for three covariates, in a
model using time since enrollment.}
\label{fig:zp1}
\end{figure}

A fit on age scale is more revealing about the NAFLD effect.
The test for PH is now highly significant. 
Separate \code{cox.zph} tests on the male (\code{cfit2m})
and female (\code{cfit2f}) fits leads to 
figure \ref{fig:zp2},  which clearly shows the interaction.  
For males the effect of
NAFLD is a classic proportional hazard effect, increasing their
risk by 1.6 fold across age, while for females the effect is
much larger at younger ages.
Interestingly, the overall time-averaged NAFLD effect for females,
which is what the coefficient from \code{cfit2f} estimates, is about the same
size as the effect for males.

<<cfit2>>=
cfit2 <- coxph(Surv(age1, age2, death) ~ nafld + male, data1)
cfit2

cox.zph(cfit2)

cfit2m <- coxph(Surv(age1, age2, death) ~ nafld, data1,
                subset= (male==1))
cfit2f <- coxph(Surv(age1, age2, death) ~ nafld, data1,
                subset= (male==0))
round(exp(c(female = cfit2f$coef, male = cfit2m$coef)),2)
@ 

\begin{figure}
<<zp2>>=
oldpar <- par(mfrow=c(1,2), mar=c(5,5,1,1))
plot(cox.zph(cfit2m, transform= 'identity'), resid=FALSE, xlab='Age',
     ylim=c(-.5, 2.5))
abline(h = coef(cfit2m), lty=2, col=2)
plot(cox.zph(cfit2f, transform= 'identity'), resid=FALSE, xlab='Age',
     ylim=c(-.5, 2.5))
abline(h = coef(cfit2f), lty=2, col=2)
par(oldpar)
@ 
\caption{Estimated NAFLD effect for males (left panel) and females (right).
The estimated overall coefficients for each separate fit is in red.}
\label{fig:zp2}
\end{figure}

To better understand what is going on, make a simple table of the
absolute hazards of death over fixed age intervals.  
Such tables of observed and expected are quite common in epidemiology where
they are called ``person-years'' analyses.
We will pick age epochs that have similar numbers of deaths.
The \code{pyears} function creates tables with the number of events and the
total amount of follow-up time in each cell.

<<rate1, fig.cap="Death rates by sex and age group.  Dashed lines are NAFLD and solid are control">>=
dtime <- with(data1, age2[death==1])
round(quantile(dtime, 1:9/10),1)  # use 50, 65, 75, 85
agecut <- tcut(data1$age1, c(0, 50, 65, 75, 85, 200),
               c("<=50", "50-65", "65-75", "75-85", ">85"))
pyfit <- pyears(Surv(age2-age1, death) ~ agecut + nafld + male, data1,
                scale=1000)
rate <- pyfit$event/pyfit$pyears

# reshape this for printing and plotting.  A 5x4 matrix is handier than 5x2x2
r2 <- matrix(rate, ncol=4) 
dimnames(r2) <- list(rownames(rate), c("Female control", "Female NAFLD", 
                  "Male control", "Male NAFLD"))
round(r2,1)

matplot(1:5, r2, log='y', pch="ffmm", col=c(2,2,1,1), lty=c(1,2,1,2),
        type= 'b', xlab="", xaxt='n', ylab="Death rate per 1000")
axis(1, 1:5, levels(agecut))
@  

Figure \ref{fig:rate1} shows yearly death rates per 1000 (absolute) and
figure \ref{fig:zp2} shows the ratios for NAFLD versus controls.
We see that NAFLD levies a particularly heavy burden on younger females,
in terms of relative rates, 
while the hazard ratio for males is fairly constant.
Interestingly, the NAFLD effect for females is  close to constant on
the $\sqrt{\mbox{rate}}$ scale.  (I don't yet have a good use for that fact.)

\subsection{Summarizing a non-proportional effect}
A formal way to understand non-proportionality is to fit an extended
Cox model that has a time dependent coefficient, e.g.,
$$
\lambda(t) = \lambda_0(t) \exp(\beta_1(t)x_1 + \beta_2 x_2 + \ldots)
$$
A first order approximation to this is plotted as part of the 
\code{cox.zph} function. 

\subsubsection{Using survSplit}
A more refined version can be based on a simple slight of hand.
Assume for instance that $\beta(t) = a + bt$ for some covariate.
Then $\beta(t)x = ax  + b(tx)$, and we can fit this with an ordinary Cox
model along with a modified data set containing the variables $x$ and
$tx$. 
Because the underlying code of a Cox model only evaluates its arguments
at the event times, the modified data would have one row per subject/time
combination, i.e., if subject Smith was at risk for 8 of the events then
he/she would have 8 rows in the data set.
Normally such a fine division of time is not necessary. 
Below, the survSplit routine is used to expand the original data set into 1 year
intervals of age, which is not continuous time but close enough to 
continuous over our 60 year scale. 
The code below shows an unimportant time effect for males but an
substantial drop off for females.  
Time dependent covariates should always depend on the \emph{past}, so our
constructed variable below uses \code{age1}, the age at the start of each
interval.

<<slight>>=
timedata <- survSplit(Surv(age1, age2, death) ~ ., data1, cut=20:99)
timedata$nafld2 <- with(timedata, nafld *age1)  # the dummy variable
tfit1 <- coxph(Surv(age1, age2, death) ~ nafld + nafld2, 
                  timedata, subset= (male==0))
tfit2 <- coxph(Surv(age1, age2, death) ~ nafld + nafld2,
                  timedata, subset= (male==1))
round(summary(tfit1)$coef,4)
round(summary(tfit2)$coef,4)

# plot, not shown
if (FALSE) {
    ages <- seq(20, 100, 5)
    tdata <- data.frame(nafld=1, nafld2 = ages)
    matplot(ages, exp(cbind(predict(tfit1, tdata), predict(tfit2, tdata))),
            col=2:1, pch='fm',  log='y', xlab="Age", ylab="Hazard ratio")
    abline(0,0, lty=2)
}
@ 

The time dependent trick gives a simple and interpretable picture of
$\beta(t)$, but with a major problem: we picked the form for $\beta(t)$
by looking at the \code{cox.zph} plot. 
As with a linear model that is chosen after looking at the scatterplot,
any inference from such a fit is suspect.  
Many authors suggest using $\beta(t) = a + b\log(t)$ based on the fact that
covaratiate effects have often been observed to fade with time;
however fitting an arbitrary fixed shape without even looking at the data
is an even worse course choosing a fit based on looking:  
we might conclude that nothing is going on simply because of a bad choice
for $\beta(t)$.
A more sophisticated approach is to fit a spline function.

<<slight2>>=
# The anova function below will not compare models fit to two different
#  size data sets.  tfit3a will have exactly the same coefficients and
#  likelihood as cfit2f above
tfit3a<- coxph(Surv(age1, age2, death) ~ nafld, timedata, subset=(male==0))
tfit3b <- coxph(Surv(age1, age2, death) ~ nafld + nafld:ns(age1,3), timedata,
               subset= (male==0))
anova(tfit3a, tfit3b)

tdata <- data.frame(age1=20:90, nafld=1)
yhat <- predict(tfit3b, tdata, se=TRUE)
yhat <- yhat$fit + outer(yhat$se.fit, c(0, -1.96, 1.96),'*')
matplot(20:90, yhat,  type='l', lty=c(1,2,2), col=1, ylim=c(-.2,5),
     xlab="Age", ylab="Estimated relative NAFLD risk, females")
abline(h=0,lty=3)
@ 

One can, of course, start by creating the full data set in \code{survSplit}
by setting the cuts value to a vector containing all of the unique death
times in data1.  
The same analysis can be done using the built in \code{tt} operator in
the coxph function.  As shown below, the log-likelihood hardly changes
between the two fits.
(Because the ns function chooses data based knots, the coefficients of the
two fits are not comparable.)

<<tt>>=
tfit3c <- coxph(Surv(age1, age2, death) ~ nafld + tt(nafld), data1,
                subset= (male==0), 
                tt = function(x, time, ...) x*ns(time, df=3))
c("survSplit loglik"= 2*diff(tfit3b$loglik), 
  "tt loglik"= 2*diff(tfit3c$loglik))
@ 

The advantage of the \code{tt} approach is that it is less work ---- you do
not need to explicitly create an expanded data set.
A primary disadvantage is that because the expanded data set is
internal to the function, follow-up analyses of the fit such as
predictions and residuals are not available. 
For large data sets the \code{tt} approach can be slower as well, since
it uses all unique event times to split the data.


\subsubsection{Timreg package}
A more formal approach is to use the
\code{timereg} package, which directly implements the extended Cox model
\begin{equation}
  \lambda(t) = \lambda_0(t)\exp(\beta_1(t)x_1 + \beta_2(t) x_2 + \ldots)
  \label{eq:time}
\end{equation}
where $\beta(t)$ are time dependent coefficients.

The \code{timecox} function currently does not handle the subset argument.
It can also have issues if the risk sets get too small, so we truncate the
estimation at age 95.  
<<timereg, fig.cap='Cumulative regression coefficient for females and males'>>=
fdata <- subset(data1, male==0) 
tcox1 <- timecox(Surv(age1, age2, death) ~ nafld, fdata,
                 max.time= 95, start.time= 30, id= fdata$id)
mdata <- subset(data1, male==1)
tcox2 <- timecox(Surv(age1, age2, death) ~ nafld, mdata,
                 max.time= 95, start.time= 30, id= mdata$id)
temp <- matrix(c(tcox1$pval.testBeq0[2], tcox1$pval.testBeqC[2],
                 tcox2$pval.testBeq0[2], tcox2$pval.testBeqC[2]), 2, 2,
               dimnames=list(c("Overall", "Constant"), c("Female", "Male")))
temp

plot(tcox1$cum[,"time"], tcox1$cum[,"nafld"], xlab="Age", ylab="Cumulative")
points(tcox2$cum[,"time"], tcox2$cum[,"nafld"], col=2)
@ 

The resulting fit shows tests for overall significance and for a
constant hazard ratio over time, e.g., proportional hazards.
For both male and female the NAFLD effect is significant, while for 
females it is strongly non-proportional.
Figure \ref{fig:timereg} shows the resulting cumulative regression
plots for males and females, the female one deviates strongly from
linear.

\subsection{Other summaries}

\section{Multi-state analysis}
\subsection{Data}
\begin{figure}
<<nfig1, echo=FALSE>>=
states <- c("No comorbidity", "1 comorbidity", "2 comorbidities", 
            "3 comorbitities", "Death")
cmat <- matrix(0, 5,5)
cmat[,5] <- 1
cmat[1,2] <- cmat[2,3] <- cmat[3,4] <- 1
cmat[1,3] <- cmat[2,4] <- 1.6
cmat[1,4] <- 1.6
dimnames(cmat) <- list(states, states)
statefig(cbind(4,1), cmat)
@   
\caption{Schema of the multi-state model.}
\label{statefig}
\end{figure}

NAFLD is closely linked to obesity and metabolic syndrome,
and Allen considered
the occurrence of three dysmetabolic comorbidities --- hypertension,
dyslipidemia, and diabetes.  
Analysis was then based on the 5 state model shown in figure \ref{statefig}.
Using `0MC', '1MC', etc. as shorthand names for the number of metabolic
comorbidities, we can define both current state and outcome
variables for each subject. 
Note one major difference between current state and outcome, namely that the
first is a time-dependent variable.
If a piece of a subject's time
were divided into multiple sub-intervals, time dependent covariates
repeat in each interval.  An endpoint, however, does not repeat.
The 2MC endpoint is not like a recurrent infection; the subject does
not get new instances of it.

<<mstate>>=
data1$cstate <- with(data1, diab + htn + dyslip)  # td covariate

# verify the e1, e2, e3 only happen once per person (we count on this)
check <- with(data1, pmax(tapply(e1, id, sum), tapply(e2, id, sum),
                          tapply(e3, id, sum)))
if (any(check >1)) cat("Duplicates exist!\n")

tcount <- with(data1, e1 + e2 + e3)
temp2 <- with(data1, ifelse(death, 4, 
              ifelse(tcount ==0, 0, cstate + tcount)))
data1$endpoint <- factor(temp2, 0:4, 
         c("censored", "1mc", "2mc", "3mc", "death"))
data1$cstate <- factor(data1$cstate, 0:3,
                c("0mc", "1mc", "2mc", "3mc"))       
with(data1, table(cstate, endpoint))
@ 

For any given row of data set \code{data1} a subject can begin
that time period in any state except death,
and end that time period either with ``no transition at this time''
which is coded as censored, or end it with a transition to a higher state. 
From the table above direct jumps of 2 or 3 comorbid states are
uncommon, but they do occur when multiple new conditions are
detected at a single visit.

\subsection{Estimated hazard rates}
The hazard rates reprise values shown in table 4 of the paper.
(Not exactly, since this data set is a subset).
These models were fit
on age scale, assuming proportional hazards for the NAFLD effect.
The small number of subjects who jump 2 states are considered under
their current state, so \code{nfit2a} below is a model for the
risk of progression to a more severe metabolic state given you are 
currently in state 0mc. 
This is essentially a fit under the constraint that the 0:1, 0:2, and 0:3
hazards have common coefficients, and similarly for the 1:2 and 1:3
transitions.
Such a reduction is necessary for state pairs that have only a small number of
transitions,
and is essentially equivalent to leaving selected interaction terms out of
a model. 
The rule of thumb that one should have 10--20 events per covariate is
operative here.

<<mfit1>>=
nfit1 <- coxph(Surv(age1, age2, death) ~ male + 
                   strata(cstate)/nafld,  data= data1)
nfit2a <- coxph(Surv(age1, age2, endpoint %in% c("1mc", "2mc", "3mc")) ~ 
                      male + nafld,
               data=data1, subset= (cstate=="0mc"))
nfit2b <- coxph(Surv(age1, age2, endpoint %in% c("2mc", "3mc")) ~ 
                    male  + nafld,
               data=data1, subset= (cstate== "1mc"))
nfit2c <- coxph(Surv(age1, age2, endpoint=="3mc") ~ male + nafld,
               data=data1, subset= (cstate=="2mc"))
@ 

<<mfit1b, echo=FALSE>>= 
table4 <- matrix(NA, 4, 5, dimnames= list(from=paste(0:3, "mc"),
                to = c("0mc", "1-3mc", "2-3mc", "3mc", "Death")))
table4[, "Death"] <- coef(nfit1)[-1]
table4[1,2] <- coef(nfit2a)[-1]
table4[2,3] <- coef(nfit2b)[-1]
table4[3,4] <- coef(nfit2c)[-1]
print(round(exp(table4), 2), na.print="")
@ 

The above table largely agrees with the values found on figure 4 of the paper.
Given that the effect of NAFLD on death was quite different for males and
females, a more careful analysis would look at males and females 
separately.

\subsection{Time in state}
The set of hazard rates is a useful summary of the model on a micro level,
but does not give a firm feel for the overall effect of NAFLD on a subject.
One way to look at this is with the $p(t)$ function, which gives the 
probability of being in each state, for any given time.

In this data set NAFLD is a time dependent covariate.  Exactly how to compute
predicted curves when there is a time dependent covariate is a hard question
and has bedeviled the statistical community for many years.  
My own thinking divides the approach into three groups according to whether
the original fit uses time dependent covariates or only baseline ones, and
whether the prediction uses time dependent covariates or not.  Some further
thoughts can be found in the ``Adjusted Survival Curves'' vignette 
in the survival package \cite{Therneau2016}.

\begin{description}
  \item [baseline] 
    Fit the model using baseline covariates as predictors and predict for
    a given subject using that subject's baseline values.  This is the classic
    method and also the simplest to grasp.  A sub-category of the approach
    is the landmark method: which is to apply this approach at some fixed
    future time point, using the covariate values that will apply at that
    instant.
  \item [partial time dependent] Fit the model using 
    time dependent covariates, and then do prediction
     for a given subject using only their baseline values.  This is essentially
     the approach of Simon and Makuch \cite{Simon84}.  The approach is
     commonly used but below I discuss some problems.
   \item [time dependent]
     Fit the model using time dependent covariates, and do prediction
     for a subject using time dependent covariates.  This is the most
     elegant but can be very difficult.
\end{description}

To better understand the issue with method 2 look at a simpler data
set and only the death endpoint.
Primary biliary cirrhosis is a progressive auto-immune liver disease. 
Chronic low level inflammation slowly but inexorably generate scar 
tissue leading eventually to cirrhosis and death.
The \code{pbcseq} data set in the \code{survival} library contains ongoing
biochemical measurements for the subjects.
Fit a simple time dependent model using the 5 variable found in the
PBC risk score \cite{Dickson89}.

<<pbc>>=
first <- pbcseq[ !duplicated(pbcseq$id),]  # first row for each
pdata <- tmerge(first[, c("id", "trt", "age", "sex")], first, id,
                death = event(futime, 1*(status==2)))
pdata <- tmerge(pdata, first, id, bili = tdc(day, bili),
                edema = tdc(day, edema), albumin= tdc(day, albumin),
                protime = tdc(day, protime))
pfit <- coxph(Surv(tstart, tstop, death) ~ age + log(bili) + edema +
              log(protime) + log(albumin), pdata)
pfit
@ 

Consider an 'ordinary' subject, who has no
edema (the most prevalent category) and median values for all the other
covariates at baseline.
For that subject obtain the predicted survival curve from the fitted model
using these baseline values, i.e., a curve of type 2 above. 
At the same time get
the overall KM of the data set and plot them together.
The result is shown in figure \ref{fig:pbc2}.
Why is the predicted curve so high?  This 'average' subject's prediction
lies completely outside the confidence bounds for the overall curve.
The answer is that this is the predicted curve for someone who starts with
average values for their liver function tests and then those values
\emph{never change}.
Such a subject, whose disease does not progress, is predicted to do quite
well; that is, of course, if such a patient were ever to exist.
The problem with method 2 is that it is all too easy to create predicted
curves for uninteresting, and in fact impossible, subjects.  

<<pbc2, fig.cap='Overlaid curves for tertiles of baseline risk score along with the time-dependent prediction'>>=
tdata <- with(first, data.frame(age=median(age), bili=median(bili), edema=0,
                    protime= median(protime), albumin= median(albumin)))
pfit1 <- survfit(Surv(futime, status==2) ~ 1, data=first)
pfit2 <- survfit(pfit, newdata=tdata)
plot(pfit1, col=1, lty=1, lwd=2, xscale=365.25, xmax=3652.5, xaxs='r',
     xlab="Years", ylab="Survival")
lines(pfit2, col=2, lty=1, lwd=2, conf.int=FALSE)
legend(730, .5, c("Overall KM", "Predicted from the model"), col=1:2, lty=1,
       lwd=2)
@ 
  
There are certainly cases where type 2 curves make sense.  For instance,
assume that we were modeling cardiac risk with smoking status as a time
dependent covariate. 
A prediction of ``your future, if you don't stop'' could be quite useful
in counseling a patient.
It is not the public health question of expected future mortality for a
cohort of current smokers, however, because the future trajectory of the
cohort as a whole will include some subjects who later quit.
For NAFLD, which is estimated to affect 1 in 4 people \cite{Tapper2018}, 
long term prediction
for a population who never get NAFLD is hard to interpret simply because
such a population is a rather nebulous concept.  

We can estimate NAFLD survival under method 1 using non-parametric 
Aalen-Johansen estimates.
For this method a control at entry remains a control from that time forward,
ignoring any change to NAFLD status -- after all, the patient for whom you want
a prediction may do the same.  NAFLD subjects at entry remain NAFLD subjects.
An underlying idea here is the notion of a direct estimate.
For some client sitting across the table from me, the oracle, I want to
create a pair of curves comparing their future to that of some other
group of comparator subjects.
The direct method is to recruit two such cohorts and follow each forward 
in time; done in some parallel universe where I can quickly complete that
followup and then return to the present time to inform my client.
In this thought experiment no subjects in either group are dropped, even
though they might gain or lose disease markers in the future.

<<ag1>>=
agfit1 <- survfit(Surv(age1, age2, endpoint) ~ male + nafld0, data1,
                  id= id, istate= cstate, start.time=40)
dim(agfit1)
# plot(agfit1[1:2,-4], lty=1:2, col=rep(1:5, each=2))
rtime <- summary(agfit1, rmean=100)$table[,3]
rtime <- matrix(rtime, nrow(agfit1), ncol(agfit1))
dimnames(rtime) <- list(c("female, control", "female, NAFLD", "male, control",
                          "male, nafld"), agfit1$states)
round(rtime[c(2,1,3,4),c(5,1,2,3,4)], 1)
@  

<<ag1b, echo=FALSE, fig.cap="Expected time in state for each group">>=
if (TRUE) {
barplot(height= t(rtime[c(2,1,3,4), c(5,1,2,3,4)]), space=c(.05, .5),
        names.arg=c("F NAFLD", "F Control", "M Control", 
                    "M NAFLD"),
        legend= TRUE, args.legend=list(x=11, y=24),
        horiz=FALSE, las=1, beside=TRUE,
        col= c(rainbow(5)[4:1], 'black'), ylab="Expected years in state")
} else {
oldpar <- par(mar=c(5.1, 7.1, 2, 2))
barplot(height= t(rtime[c(2,1,3,4), c(5,1,2,3)]), space=c(.1),
        names.arg=c("F NAFLD", "F Control", "M Control", 
                    "M NAFLD"),
        legend= TRUE, horiz=TRUE, las=1, xlim=c(0,60),
        col= rainbow(5)[4:1], xlab="Expected years in state")
par(oldpar)
}
@ 

A problem with this is the sheer number of curves: there are 4 NAFLD/sex groups
times 5 states = 20 curves.
Even a simpler plot that has only females and omits death, commented out in
the above code, is too much.  
We instead summarize the data by the expected time in state, which is shown
in the final table above and more succinctly in figure \ref{fig:ag1b}.
Comparing the middle two groups in the figure, women have a longer expected 
lifetime than men (less time in the death state) and spend somewhat fewer years
of that longer lifetime in the 3MC state.
Comparing control to NAFLD (left or right pair), NAFLD subjects
have a shorter expected lifetime, but the biggest effect is that far more of
those remaining years are spent with 2 or 3 of the metabolic comorbidies.
Recognize that a significant portion of this difference is due to the
fact that a 40 year old NAFLD subject is more likely than their matched
control to already be in the 2MC or 3MC state. 

Now repeat the curves using fitted Cox models. 
I will first do this ``by hand'' because it is illustrative. 
(Later it will be demonstrated using the mstate package.)
The key idea is that we want to replicate the underlying steps in
the Aalen-Johansen estimate, which builds off of a nstate by nstate
matrix at each event time whose $jk$ element are $\lambda_{jk}(t)$,
the hazard at that time.  Transitions that cannot occur do not need to
be filled in, nor does the diagonal of the matrix.
Each of these transitions will come from a separate Cox model.
In those prior fits we fit all the death transitions at once using
a strata by covariate interaction; the survival curve portion of the
package isn't quite smart enough to deal with that so we redo them
as separate fits (yielding the same coefficients).
The result is given numerically below and in figure \ref{fig:mstack2}.
Since it used a time dependent covariate in the model, but predicts
for a subject who is always NAFLD or non-NAFLD, this corresponds
to method 2 given earlier.

Figures \ref{fig:ag1b} and \ref{fig:mstack2} differ in another important way,
chosen to illustrate the difference.  
The former tries to address the question the future of an ``average'' NAFLD
subject who is 40 years old.  Many of these will already have comorbidities,
and at a higher rate than a matched reference subject without NAFLD.
The second figure asks the question for a given pair of subjects, both with
0 comorbidities; it is also based on the fitted model.  
The two computations answer different questions.


<<coxcurve1>>=
sname <- c("0mc", "1mc", "2mc", "3mc", "death")
lambda <- matrix(vector("list", 25), 5,5,
                 dimnames = list(from=sname, to=sname))
tdata <- expand.grid(nafld=0:1, male=0:1)

# survival curves
for (i in 1:4) {
    tfit <- coxph(Surv(age1, age2, death) ~ nafld + male, data1,
                  subset=(cstate==sname[i]))
    lambda[i,5] <- list(survfit(tfit, newdata=tdata, start.time=40))
}

# state transitions
lambda[1, 2:4] <- list(survfit(nfit2a, newdata=tdata, start.time=40))
lambda[2, 3:4] <- list(survfit(nfit2b, newdata=tdata, start.time=40))
lambda[3, 4]   <- list(survfit(nfit2c, newdata=tdata, start.time=40))

# fix a zero divide bug, found when working up this example, which
#  occurs at age 102 (beyond the range I care about)
temp <- lambda[4,5][[1]]
temp$cumhaz <- ifelse(is.finite(temp$cumhaz), temp$cumhaz, 0)
lambda[4,5] <- list(temp)

# Now proceed  (this next line isn't fast)
agfit2 <- survfit(lambda, start.time=40, p0=c(1,0,0,0,0))
rtime2 <- summary(agfit2, rmean=100)$table[,3]
rtime2 <- matrix(rtime2, nrow(agfit2), ncol(agfit2))
dimnames(rtime2) <- list(c("female, control", "female, NAFLD", "male, control",
                          "male, nafld"), agfit2$states)
round(rtime2[c(2,1,3,4),],1)
@ 

\begin{figure}
<<mstack2>>=
if (FALSE) {  # individual bars
    barplot(height= t(rtime2[c(2,1,3,4),]), space=c(.05, .5),
        names.arg=c("F NAFLD", "F Control", "M Control", 
                    "M NAFLD"),
        legend= TRUE, args.legend=list(x=11, y=28),
        horiz=FALSE, las=1, beside=TRUE,
        col= c(rainbow(5)[4:1], 'black'), 
        ylab="Expected years in state, starting in 0mc at age 40")
} else { #stacked bars
    oldpar <- par(mar=c(5.1, 7.1, 2, 2))
    barplot(height= t(rtime2[c(2,1,3,4), 1:4]), space=c(.1),
        names.arg=c("F NAFLD", "F Control", "M Control", 
                    "M NAFLD"),
        legend= TRUE, horiz=TRUE, las=1, xlim=c(0,60),
        col= rainbow(5)[4:1], 
        xlab="Expected years in state, starting in 0mc at age 40")
    par(oldpar)
}
@
\caption{Predicted average number of years in each state, for a hypothetical
  subject starting with 0 metabolic comorbidities at age 40.}
\label{fig:mstack2}
\end{figure}

Now repeat the modeling
using the \code{mstate} package.
The first step is to create a transition matrix, which guides all of
the further routines.
<<mstate1>>=
sname <- c("0mc", "1mc", "2mc", "3mc", "death")
transitions <- matrix(NA, 5, 5, dimnames= list(from=sname, to=sname))
transitions[1,2:4] <- 1:3
transitions[2,3:4] <- 4:5
transitions[3,4] <- 6
transitions[1:4, 5] <- 7:10
transitions
@ 

From this point forward all of the transitions are referred to by the
numeric value found in this matrix. 
Transitions marked as NA are impossible.
The second step is to create a special data set.  
The mstate package fits all of the models at once, unlike our code
further above which fit separate models \code{nfit1}, \code{nfit2a}, etc.
for the different transitions. 
To do this it requires a ``stacked'' data set:
the first set of rows are the ones that would be used for a Cox model
of 0mc$\rightarrow$1mc transitions, then the data needed for 
1mc$\rightarrow$2mc transitions, etc.
for all of the 10 possible transitions given in our matrix.
Each of the individual subsets is marked with a (from, to, transition) triplet
along with a new variable \code{status} which is 1 if the
transition in question occurred. 
A special function \code{msprep} can create such a data set from a particular
kind of input data set, but that approach does not allow for 
time dependent covariates.

<<mstate2>>=
stacker <- function(s1, s2, trans) {
    temp <- subset(data1, as.numeric(cstate) == s1)
    data.frame(from=s1, to=s2, trans=trans, 
               status = 1*(as.numeric(temp$endpoint) == s2),
               temp)
}
ntrans <- max(transitions, na.rm=TRUE)
temp <- vector("list", ntrans)
for (i in 1:5) {
    for (j in 1:5) {
        k <- transitions[i,j]
        if (!is.na(k)) 
            temp[[k]] <- stacker(i, j, k)
    }
}
msdata <- do.call(rbind, temp)
attr(msdata, "trans") <- transitions
class(msdata) <- c("msdata", "data.frame")
events(msdata)$Frequencies
@ 

The result of the \code{events} function can be used as a check for
a proper data set creation.
It is also a good idea to select 2 or 3 subjects and print out their data.
We can now fit the models for all of the transitions.
Because we decided to lump transitions from 0MC to 1, 2, or 3MC together as
a single rate, we need to create a variable that lumps them.

<<mstate3>>=
msdata$trans2 <- c(1,1,1,2,2,6:10)[msdata$trans]
mcox <- coxph(Surv(age1, age2, status) ~ strata(trans2)/(male + nafld), 
              msdata)
mcox
@ 

The coefficients of this model are the same as individual fits.  The
advantage is that we can get all of the predictions at once, using
a data frame that has all of the states.  This prediction matrix needs
a set of rows for each of the transitions and needs to be of the
\code{msdata} class.

<<mstate4>>=
temp <- msdata[match(1:10, msdata$trans), c('trans', 'male', 'nafld', 'trans2')]
pdata <- expand.grid(nafld=0:1, male=0:1, trans=1:10,
                     trans= c(1,1,1,2,2,6:10))
attr(pdata, "trans") <- transitions
class(pdata) <- c("msdata", "data.frame")

# there is a problem with my mhaz call, not yet addressed: the next line fails
#mhaz <- msfit(mcox, newdata=pdata, trans=transitions)
@ 

\section{Using a detailed time line}

\begin{figure}
<<timeline, echo=FALSE>>=
frame()
plot(c(0, 75), c(0,2), xlab="Time", ylab="", yaxt='n', type='n')
segments(c(10, 26, 34), c(1,1,1), c(18, 31, 65), c(1,1,1), lwd=2)
arrows(c(5, 15, 21, 26, 31, 70), rep(c(1.4, .6), 3), c(5, 15, 21, 26, 31, 70),
       rep(c(1.05, .95), 3), angle=20, length=.1)
text(c(5, 15, 21, 26, 31, 70), rep(c(1.5, .5), 3),
     c("early", "within", "gap", "leading", "trailing", "late"))
@ 
 \caption{Time Line for a hypothetical subject who moved to the region at
age 10, left for college at 18, returned at 26 for a residency, left again
for further training at ages 31--32, returned at 33, and the retired to a
warmer climate at age 70.}
\label{fig:timeline}
\end{figure}

A fourth data set \code{nafld0} is special to the REP.  A given
subject may be in the Olmsted County population for only a short period,
since birth, or even for disjoint intervals;
figure \ref{fig:timeline} gives an example. Hence the definition of
the at risk period is a bit more complicated than a single follow-up time.
This data set contains begin/end intervals for each subject.  
We have eliminated time line information before the index date, 
so each subject's first interval starts at or after 0.  

Let us repeat the exercise using the full REP time line data.
Formally, this is the correct way to proceed since it accounts for the
detail information about who is at risk.
The tmerge calls are similar to the prior set.

<<data2>>=
data2 <- tmerge(nafld1[,1:8], nafld0, id=id, tstart=begin, tstop=end) 
data2 <- tmerge(data2, subset(nafld3, event=="death"), id=id, 
                              death= event(days))
data2 <- tmerge(data2, subset(nafld3, event=="nafld"), id, 
                nafld= tdc(days))
data2 <- tmerge(data2, subset(nafld3, event=="diabetes"), id,
                diab= tdc(days), e1= event(days))
data2 <- tmerge(data2, subset(nafld3, event=="htn"), id, 
                htn= tdc(days), e2= event(days))
data2 <- tmerge(data2, subset(nafld3, event=="dyslipidemia"), id, 
                dyslip = tdc(days), e3= event(days))
data2$age1 <- with(data2, age + tstart/365.25)
data2$age2 <- with(data2, age + tstop/365.25)
                   
attr(data2, 'tcount')
@ 

<<echo=FALSE>>=
tc <- attr(data1, 'tcount')   # shorter name for use in Sexpr below
excess <- with(data1, sum(age2-age1))/ with(data2, sum(age2-age1))
excess.pct <- 100*(excess -1)
@ 
The overall number of person-years of follow-up in this data is the
sum of (age2-age1).
It turns out that \code{data1} has an excess total follow-up of of about 
\Sexpr{round(excess.pct,1)}\%, so we do not expect the answers to change much.
The tcount matrix has a few new features.


After adding in the state variables in the same way as before, the 
estimated state transition effects end up as below.
<<state2>>=
data2$cstate <- with(data2, diab + htn + dyslip)  # td covariate

# verify the e1, e2, e3 only happen once per person (we count on this)
check <- with(data2, pmax(tapply(e1, id, sum), tapply(e2, id, sum),
                          tapply(e3, id, sum)))
if (any(check >1)) cat("Duplicates exist!\n")

tcount <- with(data2, e1 + e2 + e3)
temp2 <- with(data2, ifelse(death, 4, 
              ifelse(tcount ==0, 0, cstate + tcount)))
data2$endpoint <- factor(temp2, 0:4, 
         c("censored", "1mc", "2mc", "3mc", "death"))
data2$cstate <- factor(data2$cstate, 0:3,
                c("0mc", "1mc", "2mc", "3mc"))       

nfit3 <- coxph(Surv(age1, age2, death) ~ male + 
                   strata(cstate)/nafld,  data= data2)
nfit4a <- coxph(Surv(age1, age2, endpoint %in% c("1mc", "2mc", "3mc")) ~ 
                      male + nafld,
               data=data2, subset= (cstate=="0mc"))
nfit4b <- coxph(Surv(age1, age2, endpoint %in% c("2mc", "3mc")) ~ 
                    male  + nafld,
               data=data2, subset= (cstate== "1mc"))
nfit4c <- coxph(Surv(age1, age2, endpoint=="3mc") ~ male + nafld,
               data=data2, subset= (cstate=="2mc"))

# print out a table comparing these models to prior ones
temp <- cbind(c(coef(nfit1), coef(nfit2a), coef(nfit2b), coef(nfit2c)),
              c(coef(nfit3), coef(nfit4a), coef(nfit4b), coef(nfit4c)))
temp <- temp[c(2:5, 7,9,11), ]  #only the death coefs
dimnames(temp) <- list(c("0mc:death", "1mc:death", "2mc:death", "3mc:death",
                         "0mc:1mc", "1mc:2mc", "2mc:3mc"),
                       c("Simple", "Timeline"))
round(exp(temp), 2)
@ 

We see that adding in the careful time line data makes almost no change in
the Cox model fits, at least
for this data set.  We will not pursue the timeline data further.

\bibliography{refer}
\bibliographystyle{plain}
\end{document}