\SweaveOpts{results=verbatim,keep.source=TRUE,include=FALSE,eps=FALSE,prefix.string=MS}
\documentclass[a4paper,tipsiness,twoside,12pt]{report}

% ----------------------------------------------------------------------
% General information for the title page and the page headings
\newcommand{\Title}{EPIC multistate models}
\newcommand{\Tit}{EPIC-MS}
\newcommand{\Version}{Version 3}
\newcommand{\Dates}{April 2018}
\newcommand{\Where}{SDCC: BxC, IARC: P Ferrari, H Lennon}
\newcommand{\Homepage}{\url{http://bendixcarstensen.com/} }
\newcommand{\Faculty}{\begin{tabular}{rl}
Bendix Carstensen
  & Steno Diabetes Center Copenhagen, Gentofte, Denmark\\
  & {\small \& Department of Biostatistics,
               University of Copenhagen} \\
  & \texttt{b@bxc.dk}\\
  & \url{http://BendixCarstensen.com} \\[1em]
                      \end{tabular}}

\input{/home/bendix/util/tex/topreport}
\renewcommand{\rwpre}{MS}

%----------------------------------------------------------------------
% Here comes the substance part(s)
%
\chapter{Multimorbidity and mortality in EPIC}

\section{3-state modeling}

We want to model occurrence of DM, CVD and Cancer in the EPIC cohorts,
and assess to which extent these depend on baseline characteristics.

Fist we visualize the rates we are concerned about:
<<>>=
states <- c("Well","DM","CVD","Ca","DM+CVD","DM+Ca","CVD+Ca","DM+CVD+Ca","Dead")
ns <- length(states)
TM <- matrix(NA,ns,ns)
rownames(TM) <- colnames(TM) <- states
TM["Well",c("DM","CVD","Ca","Dead")] <- 1
TM[c("DM","Ca"),c("DM+Ca","Dead")] <- 1
TM[c("DM","CVD"),c("DM+CVD","Dead")] <- 1
TM[c("Ca","CVD"),c("CVD+Ca","Dead")] <- 1
TM[c("DM+CVD","CVD+Ca","DM+Ca"),c("DM+CVD+Ca","Dead")] <- 1
TM[c("DM+CVD+Ca"),c("Dead")] <- 1
TM
@ %
With the transition matrix we can set up the graph of transitions:
<<boxes,fig=TRUE,height=12,width=12>>=
library(Epi)
# x11(h=16,w=24,p=24)
boxes.matrix(TM,boxpos=list(x=c(10,35,35,35,65,65,65,90,50),
                            y=c(40,60,40,20,60,40,20,40,90)-10),
             hmult=2,wmult=1.2,col.arr="gray",
             col.border=c(rep("transparent",8),"black"))
par(new=TRUE)
boxes.matrix(TM[-9,-9],boxpos=list(x=c(10,35,35,35,65,65,65,90),
                                   y=c(40,60,40,20,60,40,20,40)-10),
             hmult=2,wmult=1.2,lwd.arr=3,col.bg=gray(0.9),
             col.border=gray(0.9))
@ %
\insfig{boxes}{1.0}{States of interest in the MM study. Note that the
  NE arrows correspond to T2D occurrence, the E arrows to CVD
  occurrence and the SE arrows to Ca occurrence. Also note that we
  have not separated persons with CVD before Ca from persons with Ca
  before CVD; they are all in the state CVD+Ca. Implicitly this
  assumes that the T2D occurrence rates are the same regardless of
  which of the two conditions occurred first.}
There are two sets of analyses of immediate interest, one addressing
the disease rates --- the black rates in figure \ref{fig:boxes} (there
are 12 of these), and
one addressing the mortality rates --- the gray lines (there are 8 of
these rates).

\section{Setting up data for analysis}

The base dataset contains the relevant dates for each person:
\begin{itemize}
\item \texttt{doBth} --- date of birth
\item \texttt{doExm} --- date of examination (entry to EPIC)
\item \texttt{doDM} --- date of DM
\item \texttt{doCVD} --- date of CVD
\item \texttt{doCa} --- date of Cancer
\item \texttt{doDth} --- date of death
\item \texttt{doX} --- date of exit, \textit{i.e.} censoring, last
  date of follow-up for \emph{all} events.
\end{itemize}
Suppose the data set is called \texttt{EPIC}, the we set up a Lexis object by
first using death:
<<eval=FALSE>>=
L0 <- Lexis( entry = list(per=doExm,
                          age=doExm-doBth,
                          tfx=0),
              exit = list(per=doX),
       exit.status = factor(!is.na(doDth),labels=c("Well","Dead")),
              data = EPIC)
@ %
In order to get the follow-up split between the various other states we
cut follow-up in different states:
<<eval=FALSE>>=
LD <- cutLexis( subset(L0,!is.na(doDM) & 
                          (doDM<doCVD|is.na(doCVD)) & 
                          (doDM<doCa |is.na(doCa )) ),
                cut=doDM,
                new.state="DM",
                pre="Well" )
LV <- cutLexis( subset(L0,!is.na(doCVD) & 
                          (doCVD<doDM|is.na(doDM)) & 
                          (doCVD<doCa|is.na(doCa)) ),
                cut=doDM,
                new.state="CVD",
                pre="Well" )
LC <- cutLexis( subset(L0,!is.na(doCa) & 
                          (doCa<doDM |is.na(doDM )) & 
                          (doCa<doCVD|is.na(doCVD)) ),
                cut=doDM,
                new.state="Ca",
                pre="Well" )
L0 <- subset( L0, is.na(doDM) & is.na(doCVD) & is.na(doCa) )
@ %
The three data sets \texttt{LD}, \texttt{LV} and \texttt{LC} now has
the follow-up from all person with at least one data of diagnosis,
subdivided at the first diagnosis, and \texttt{L0} that of the
remaining persons (\textit{i.e} persons without any disease event).

The three should now be further subdivided according to the next
event, for example for those with diabetes as first event, into those
with no further event (\texttt{LD0}), those with a CVD event
(\texttt{LDV}) and those with a cancer event (\texttt{LDC}):
<<eval=FALSE>>=
LD0 <-           subset(LD, is.na(doCVD) & is.na(doCa) )
LDV <- cutLexis( subset(LD,!is.na(doCVD) &
                            doCVD<doCa|is.na(doCa)),
                 cut=doCVD,
                 new.state="DM+CVD" )
LDC <- cutLexis( subset(LD,!is.na(doCa) &
                            doCa<doCVD|is.na(doCVD)),
                 cut=doCa,
                 new.state="DM+Ca" )
@ % 
Similar subdivisions is then made for persons with CVD, resp. cancer as
first event, and finally those with a third event are considered from
the 6 subsets with two initial events. 

Finally all data sets (Lexis objects) are collected in a single object,
\texttt{LE}, say, and using this we can analyze rates (as for example
outlined in \cite{Iacobelli.2013}) by a simple \texttt{glm} with
Poisson error:
<<eval=FALSE>>=
glm( (lex.Xst=="Ca") ~ Ns(age) + sex,
     family = poisson, 
     offset = log(lex.dur),
       data = subset(LE,lex.Cst=="Well") )
@ %

\section{4 states}

Now suppose that we instead have 4 disease states of interest, DM, Ca,
CHD and Str. The we must use these 4 disease states and all possible combinations of
these as possible states:
<<>>=
dn1 <- c("DM","Ca","CHD","Str")
dn2 <- outer(dn1,dn1,paste,sep="+")
dn2 <- dn2[upper.tri(dn2)]
dn3 <- NULL
for( i in 4:1 )
dn3  <- c(dn3,paste( dn1[-i], collapse="+" ))
dn4 <- paste( dn1, collapse="+" )
dn1
dn2
dn3
dn4
TM <- matrix(NA,15,15)
rownames(TM) <- colnames(TM) <- c(dn1,dn2,dn3,dn4)
TM[1,c(5,6,8)] <-
TM[2,c(5,7,9)] <-
TM[3,c(6,7,10)] <-
TM[4,c(8,9,10)] <-
TM[5,c(11,12)] <-
TM[6,c(11,13)] <-
TM[7,c(11,14)] <-
TM[8,c(12,13)] <-
TM[9,c(12,14)] <-
TM[10,c(13,14)] <-
TM[11,15] <-
TM[12,15] <-
TM[13,15] <-
TM[14,15] <- 1
@ %
We also need an initial ``Well'' state (without recorded disease), from
which only the 4 states with one disease present can be reached:
<<>>=
TM <- rbind( c(NA,rep(1,4),rep(NA,11)),
             cbind( NA, TM ) )
rownames(TM)[1] <-
colnames(TM)[1] <- "Well"
@ %
Finally we should expand with a death state which can be reached from
any other state:
<<>>=
TM <- rbind( cbind( TM, 1 ), NA )
rownames(TM)[nrow(TM)] <-
colnames(TM)[ncol(TM)] <- "Dead"
TM
@ %
Drawing the transitions it is relevant to color the arrows to get a
proper overview --- transition rates are colored by a color
corresponding to the disease \emph{added} to the current
state. Technically, we draw the same diagram twice; first with all
arrows in gray, and subsequently (on top of the first one) with the
colored incidence arrows, so that we have the mortality arrows
underneath.
<<tbox4,fig=TRUE,height=10,width=12>>=
clr <- c("limegreen","blue","brown","orange",gray(0.7))
# Order colors according to the ordering convention of transitions: 
o.clr <- c(1:4,
           
           2:4,
         1,3,4,
         1,2,4,
           1:3,
         
           3,4,
           2,4,
           1,4,
           2,3,
           1,3,
           1,2,
             
             4,
             3,
             2,
             1)
# Position of the boxes
bp <- list(x=c(3,
               rep(20,4),
               c(50,50,43,57,50,50)-5,
               rep(73,4),92,
               95),
           y=c(c(50,
                 seq(78,22,,4),
                 seq(97,03,,5)[c(1,2,3,3,4,5)],
                 seq(78,22,,4),50)*0.95,
               95))
boxes.matrix( TM, boxpos=bp,
              hmult=2, wmult=1.3, eq.wd=FALSE,
              col.bg=gray(0.95),col.border=gray(0.95),
              col.arr=clr[5] )
par(new=TRUE)
boxes.matrix( TM[-17,-17], boxpos=lapply(bp,function(x)x[-length(x)]),
              hmult=2, wmult=1.3, eq.wd=FALSE,
              col.bg=gray(0.9),col.border=gray(0.9),
              col.arr=clr[o.clr], lwd.arr=4 )
@ %
\insfig{tbox4}{1.0}{Expanded model with both stroke and CHD as possible
diagnoses of interest. Green arrows are DM occurrence, blue cancer
occurrence, orange stroke occurrence, maroon arrows CHD occurrence,
while gray arrows represent deaths.}

In principle it should be possible to model all 32 incidence rates (8
of each kind) and all 16 mortality rates shown in figure
\ref{fig:tbox4} separately, but hardly feasible to report, let alone
summarize in a sensible form. The relevant approach would be to fit
separate models for each \emph{color} of rates, using the originating
state for each rate as covariate. The rates would in addition depend
on sex, age, calendar time, country and population/EPIC status.

The big task is setting up the follow-up; fairly simple for the EPIC
data --- just as outlined for the 3-event-type situation above. But
this is hardly doable for the population data for want of data --- it
is not realistic to get population rates corresponding to all arrows
in figure \ref{fig:tbox4}; population rates are at the best available
as sex-, age- and period-specific rates; that is weighted averages of
arrows of the same color. 

Thus assuming that only EPIC data can be used for the modeling illustrated
in figure \ref{fig:tbox4}, each set of rates can be related to crude
population rates of DM, cancer, stroke and CHD.

Likewise, mortality rates \emph{from} each of the 16 states can be
related to the overall mortality rates. 

Note that cause-specific mortality rates are not discussed at all
here; only mortality rates \emph{from} different disease states. In
the vein of relative survival it may be deliberated whether additive
rather than multiplicative models might be relevant for description of
the differences between mortalities in the EPIC cohort (gray lines).

\input{read}
\input{ana}
\input{comments1}

\renewcommand{\bibname}{References}
\bibliographystyle{unsrt}
\bibliography{%
/home/bendix/art/bibtex/BxC,%
/home/bendix/art/bibtex/Stat,%
/home/bendix/art/bibtex/DMCa,%
/home/bendix/art/bibtex/Diabetes%
}
\addcontentsline{toc}{section}{\bibname}

\end{document}