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. But it is only a set of wrappers fro glm, gam and coxph.

Usage

glm.Lexis( Lx,         # Lexis object	
           resp,         # 'to' states	
        formula,         # ~ model	
           xpos,         # 'from' states
           link = "log", # link function
          scale = 1,     # scaling of PY
        verbose = TRUE,  # report what is done?
             ... )
  gam.Lexis( Lx,         # Lexis object	
           resp,         # 'to' states	
        formula,         # ~ model	
           xpos,         # 'from' states
           link = "log", # link function
          scale = 1,     # scaling of PY
        verbose = TRUE,  # report what is done?
             ... )
coxph.Lexis( Lx, # Lexis object	
           resp, # 'to' states	
        formula, # timescale ~ model	
           xpos, # 'from' states
        verbose = TRUE,  # report what is done?
             ... )

Arguments

A Lexis object representing cohort follow-up.

resp

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.

formula

Model formula describing the model for the intensity. 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 as the l.h.s.

xpos

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. If missing (that is not supplied), the entire Lexis object is used in the modeling.

link

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?

…

Arguments passed on to the methods.

Value

glm returns a glm object, gam returns a gam object and coxph returns a coxph object. The returned objects all have an extra attribute, Lexis; a list with names Exposure and Events; character vectors of state names from which transitions are modeled and that are considered events, respectively. The glm object also has a scale element in the list, an scalar indicating the scaling of lex.dur before modeling. The coxph object also has a Timescale element in the list, a character indicating the underlying timescale variable.

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 returns the estimated rates in units of the lex.dur in the Lexis object, scaled by scale.

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.

The functions are slightly experimental so far. Argument names and ordering may change in the future. The update methods do not always work.

Examples

Run this code

# NOT RUN {
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") )
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
mt <- glm.Lexis( sL, resp="Dead",
                     ~ 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 takes quite some time so we leave it out
# }
# NOT RUN {
mt.gam <- gam.Lexis( sL, "Dead", ~ sex + lex.Cst + s(age), scale=1000 )
        
# }
# NOT RUN {
# Fit a Cox model with age as baseline time scale and insulin as time-dependent
mt.cox <- coxph.Lexis( sL, "Dead", age ~ sex + lex.Cst )

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