Epi (version 2.56)

mod.Lexis: Fit intensity models to follow-up data in Lexis objects

Description

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.

Usage

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 coxph

Value

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.

Arguments

Lx

A Lexis object representing cohort follow-up.

formula

Model formula describing the model for the intensity(-ies). For glm and gam, the formula should be one-sided; for coxph the formula should be two-sided and have the name of the time-scale used for baseline hazard as the l.h.s.

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 lex.Cst. Defaults to the collection of transient states immediately preceding the absorbing states.

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 lex.Xst. Defaults to the set of absorbing states.

paired

Logical. Should the states mentioned in to, rep. from be taken as pairs, indicating the only transitions modeled. If FALSE all transitions from any of the states in from to any states in to are modeled.

link

Character; name of the link function used, allowed values are 'log' (the default), 'identity' and 'sqrt', see the family poisreg.

scale

Scalar. lex.dur is divided by this number before analysis, so that you can get resulting rates on a scale of your wish.

verbose

Print information on the states modeled?

...

Further arguments passed on to glm, glm or coxph

Author

Bendix Carstensen, http://bendixcarstensen.com.

Details

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.

See Also

Lexis, cutLexis, mcutLexis, addCov.Lexis, absorbing, transient

Examples

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
if (FALSE) {
mt.gam <- gam.Lexis( sL, ~ sex + lex.Cst + s(age), to="Dead",
                     scale=1000 )
        }

# 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" )