Fit intensity models to follow-up data in Lexis objects
Modeling intensities based on Lexis objects, exploiting the structure of the
Lexis objects where the events and risk time have predefined
representations. This allows a simpler syntax than the
traditional explicit modeling using glm, gam
and coxph. Requires that lex.Cst and lex.Xst
are defined as factors.
But it is just a set of wrappers fro
glm, gam and coxph.
glm.Lexis( Lx, # Lexis object
formula, # ~ model
from = preceding(Lx,to), # 'from' states
to = absorbing(Lx) , # 'to' states
paired = FALSE, # only the pairwise
link = "log", # link function
scale = 1, # scaling of PY
verbose = TRUE, # report what is done?
... ) # further arguments to glm
gam.Lexis( Lx, # Lexis object
formula, # ~ model
from = preceding(Lx,to), # 'from' states
to = absorbing(Lx) , # 'to' states
paired = FALSE, # only the pairwise
link = "log", # link function
scale = 1, # scaling of PY
verbose = TRUE, # report what is done?
... ) # further arguments to gam
coxph.Lexis( Lx, # Lexis object
formula, # timescale ~ model
from = preceding(Lx,to), # 'from' states
to = absorbing(Lx) , # 'to' states
paired = FALSE, # only the pairwise
verbose = TRUE, # report what is done?
... ) # further arguments to coxphLx |
A |
formula |
Model formula describing the model for the
intensity(-ies). For |
from |
Character vector of states from which transitions
are considered. May also be an integer vector in which case the
reference will be to the position of levels of
|
to |
Character vector of states to which a transition is
considered an event. May also be an integer vector in which case the
reference will be to the position of levels of |
paired |
Logical. Should the states mentioned in |
link |
Character; name of the link function used, allowed values
are |
scale |
Scalar. |
verbose |
Print information on the states modeled? |
... |
Further arguments passed on to |
The glm and gam models are fitted using the family
poisreg which is a bit faster than the traditional
poisson family. The response variable for this family is a
two-column vector of events and person-time respectively, so the
predictions, for example using ci.pred does not require
lex.dur (and would ignore this) as variable in the
newdata. ci.pred will return the estimated rates in
units of the lex.dur in the Lexis object, scaled by
scale, which has a default value of 1.
The default is to model all transitions into any absorbing state by
the same model (how wise is that??). If only from is given,
to is set to all states reachable from from, which may
be a really goofy model and if so a warning is issued. If only
to is given, from is set to the collection of states
from which to can be reached directly — see
preceding and its cousins. This convention means that if
you have a Lexis object representing a simple survival
analysis, with states, say, "alive" and "dead", you can dispense with
the from and to arguments.
Occasionally you only want to model a subset of the possible
transitions from states in from to states in to, in
which case you specify from and to as character vectors
of the same length and set paired=TRUE. Then only transitions
from[i] to to[i], i=1,2,... will be modeled.
There is no working update functions for these objects (yet).
Strictly speaking, it is a bit counter-intuitive to have the time-scale
on the l.h.s. of the formula for the coxph since the time scale
is also a predictor of the occurrence rate. On the other hand, calling
coxph directly would also entail having the name of the time
scale in the Surv object on the l.h.s. of the formula. So the
inconsistency is merely carried over from coxph.
glm.Lexis returns a glm object, which is
also of class glm.lex,
gam.Lexis returns a gam object, which is
also of class gam.lex, and
coxph.Lexis returns a coxph object, which is
also of class coxph.lex. These extra class attributes are meant
to facilitate the (still pending) implementation of an update function.
The returned objects all have an extra attribute, Lexis which
is a list with entries
data, the name of the Lexis object modeled (note that it
is not the object, only the name of it, which may not be portable);
trans, a character vector of transitions modeled;
formula, the model formula; and
scale, the scaling applied to lex.dur before modeling.
Only the glm and gam objects have the scale element
in the list; a scalar indicating the scaling of lex.dur before
modeling. Note that the formula component of the Lexis
attribute of a coxph object is a
two-sided formula with the baseline time scale as the l.h.s.
Bendix Carstensen, http://bendixcarstensen.com.
library( Epi )
library( survival )
data( DMlate )
# Lexis object of total follow-up
mL <- Lexis( entry = list(age=dodm-dobth,per=dodm),
exit = list(per=dox),
exit.status = factor(!is.na(dodth),labels=c("Alive","Dead")),
data = DMlate )
# Cut follow-up at start of insulin use
cL <- cutLexis( mL, cut = mL$doins,
timescale = "per",
new.state = "Ins",
precursor.states = "Alive" )
# Split follow-up on age-axis
system.time( sL <- splitLexis( cL, breaks=0:25*4, time.scale="age") )
# ( consider splitMulti from the popEpi package )
summary( sL )
# glm models for rates based on the time-split dataset by insulin and sex
# Proportional hazards model with insulin as time-dependent variable
# - uses the defaul of modeling all transitions from both transient
# states ("Alive" and "Ins") to the absorbing state ("Dead").
mt <- glm.Lexis( sL, ~ sex + lex.Cst + Ns(age,knots=c(15,3:8*10)) )
# prediction of mortality rates from "Alive" with and without PH assumption
nA <- data.frame( age=40:70, sex="M", lex.Cst="Alive" )
nI <- data.frame( age=40:70, sex="M", lex.Cst="Ins" )
matshade( nA$age, cbind( ci.pred(mt,nA),
ci.pred(mt,nI) )*1000, plot=TRUE,
lwd=3, lty=1, log="y", col=c("black","blue","red"),
xlab="Age", ylab="Mortality per 1000 PY" )
# gam models may take some time to run so we leave it out
## Not run:
mt.gam <- gam.Lexis( sL, ~ sex + lex.Cst + s(age), to="Dead",
scale=1000 )
## End(Not run)
# Fit a Cox model for mortality with age as baseline time scale and
# insulin (lex.Cst) as time-dependent covariate
mt.cox <- coxph.Lexis( sL, age ~ sex + lex.Cst, c("Alive","Ins"), "Dead" )
# Pretty much the same results for regression paramters as the glm:
ci.exp( mt , subset="ex" )
# ci.exp( mt.gam, subset="ex" )
ci.exp( mt.cox, subset="ex" )Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.