Who needs the Cox model anyway?
SDCC
November 2023

http://bendixcarstensen.com/WntCma.pdf
http://bendixcarstensen.com/WntCma.R
Version 4

Compiled Sunday 19thNovember, 2023, 17:50

from: C:\Bendix\teach\AdvCoh\art\WntCma/WntCma.tex

Bendix Carstensen

Steno Diabetes Center Copenhagen, Herlev, Denmark& Department of Biostatistics, University of Copenhagen
bcar0029@regionh.dk b@bxc.dk
http://BendixCarstensen.com

1

Theory1.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

Time: Response or covariate?

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.2.1

Likelihood for a single rate . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.2.2

Poisson likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

The Cox-likelihood as a prole likelihood . . . . . . . . . . . . . . . . . . . . . .

4

1.3
1.4

2

1

Practical data processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.4.1

Estimation of baseline hazard. . . . . . . . . . . . . . . . . . . . . . . .

6

1.4.2

Estimation of survival function. . . . . . . . . . . . . . . . . . . . . . .

7

Examples

9

2.1

9

Equality of Cox and Poisson modeling: The lung cancer example . . . . . . . . .
2.1.1

Parametric baseline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.1.2

Rates, cumulative rates and survival. . . . . . . . . . . . . . . . . . . .

13. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

Natural spline vs. penalized splines . . . . . . . . . . . . . . . . . . . . .

15

Parametric models
Cox models

Comparison with Cox
2.1.3. . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

Practical splitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

2.2

Stratied models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

2.3

Time-varying coecients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

Parametrizations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

Simplifying code

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

3

So who do need the Cox-model?

25

4Simple code for parametric model

26

4.1

4.2

Cookbook for a parametric survival curve . . . . . . . . . . . . . . . . . . . . . .

26

4.1.1

In summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

4.1.2

Comparison with K-M: . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

4.1.3

The rates: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

Parametric proportional hazards model . . . . . . . . . . . . . . . . . . . . . . .

30

4.2.1

31

Proportional hazards . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

References

35iiChapter 1
Theory1.1 Introduction
First it should be noted that the content of this note by no means is original; John Whitehead
devised a similar machinery already in 1980 [3], using GLIM. Here things are laid out using R.
In the last about 50 years, survival analysis has been virtually synonymous with applicationof the Cox-model. The common view of survival analysis (and teaching of it) from the
Kaplan-Meier-estimator to the Cox-model is based on time as the response variable,
incompletely observed (right-censored). This has automatically lent a certain aura ofcomplexity to concepts such as time-dependent covariates, stratied analysis, delayed entry
and time-varying coecients.
More unfortunate, however, is that the use of this particular technique for survival analysis
has become a dominant tool in epidemiology too, largely restricting models for occurrence
rates to models with only one time scale.
If survival studies is viewed in the light of the demographic tradition, the basic observation
is not one time to event (or censoring) for each individual, but rather many small pieces offollow up from each individual. This makes concepts clearer as modeling of rates rather thantime to response becomes the focus; the basic response is now a 0/1 outcome in each interval,
albeit not independent anymore, but still with a likelihood which is a product across intervals.
In this set up, time(scale) is then correctly viewed as a covariate rather than a response.
From a practical point of view time-dependent covariates will not have any special status
relative to other covariates. Stratied analysis becomes a matter of interaction between timeand a categorical covariate, and time-varying coecients becomes interactions between timeand a continuous covariate. Finally, the modeling tools needed reduces to Poisson regression(and ultimately logistic regression)  standard generalized linear models.

ad absurdum

The Cox-model may actually be viewed as a special case of a Poisson model where thedetail in modeling of the time covariate has been taken, namely with one

parameter per failure time. The main advantage of the demographic view is therefore that
researchers will be forced to explicitly consider which time-scale(s) to use and to what degreeof detail it is relevant to model interactions between time scales and other covariates.
Contrary to this, Poisson modeling of disease rates and follow-up studies in epidemiology
has traditionally (and until 1990 for good computational reasons) been restricted to analysis
of tables where rates have been assumed constant over fairly broad time-spans, typically 5
years, as most methods have been developed in cancer epidemiology, where 5 years is
considered a short age-span. This approach is essentially one where initial tabulation of data

1

2

1.2

Time: Response or covariate?

WmtCmaunnecessarily limits the exibility of modeling (and discards information).
If follow-up time both in survival and cohort studies are considered in small intervals, thesmoothing of rates can be done with standard regression tools in Poisson modeling.
The practicalities of this type of modeling focus requires a splitting of the follow-up in
many small intervals, and hence Poisson modeling of datasets with many records, each
representing a small piece of the follow-up time.
The only remaining advantage of the Cox-model is the ability to easily produce estimates ofsurvival probabilities in (clinical) studies with a well-dened common entry time for all
individuals, and hence with a single timescale. This can however also be produced from amodel using a smooth parametric form for the occurrence rates.

1.2 Time: Response or covariate?
Both, actually.
One common exposition of survival analysis is as analysis of data (X, Z), where we only
observe min(X, Z) and δ = 1{Z < X}. This is an approach which takes the survival time X ,as response variable, albeit not fully observed, limited by the censoring time, Z .

when

The snag here is that apart from being a measure of

X is also a measure of

how long

the person has been at risk,

the person was at risk; if we refer to time since entry to thestudy, namely at times 0 through X . The entry time to the study is implicitly taken to be 0,
so easy to confuse risk time (X − 0) and time scale ((0, X]).
From a life-table (demographic) point of view, we will break the interval (0, X] in, say,i.e. risk time

one-year intervals, each one with a dierent time associated (be that age or time since entry).
In this sense time is a covariate, and only dierences () on a timescale should beconsidered responses. In a life-table, dierences on the timescale (interval length) are
accumulated as risk time whereas the position on the age-scale for these are used as acovariate classifying the table.
Now consider a follow-up (survival) study where the follow-up time for each individual is
divided into small intervals of equal length y , say, and each with an exit status recorded (thiswill be 0 for the vast majority of intervals and only 1 for the last interval for individuals
experiencing an event).

empirical rate (d, y)Each small interval for an individual contributes an observation of what I will term an,, where d is the number of events in the interval (0 or 1), and y is thelength of the interval, i.e. the risk time. This is slightly dierent from the traditionaldenition of an empirical rate as d/y (or

P

d/

P

y ); it is designed to keep the entire

information content in the demographic observation, even if the number of events is 0. This isin order to make it usable as a response variable, showing that in a follow-up study the basic

theoreticalobservation is a (set of ) rate(s).
Therate of event occurrence is dened as a function, usually depending on sometimescale, t:P{event in (t, t + h]| at risk at time t}h↘0hλ(t) = lim

The rate may depend on any number of covariates; incidentally on none at all. Note that inthis formulation time(scale) t has the status of a covariate and h the status of risk time, h is
the dierence between two points on the timescale (in this case t + h and t).

1.2

Time: Response or covariate?

1.2.1

Theory3

Likelihood for a single rate

To derive the likelihood for a rate we must compute the probability of an observation from asingle person as a function of the rate λ. Suppose a person is alive from time te (entry) to tx(exit) and that the person's status at tx is d, where d = 0 means alive and d = 1 means dead.
If we choose, say, two time points, t1 , t2 between te and tx , standard use of conditional
probability (formally, repeated use of Bayes' formula) givesP{d at tx | entry at te } = P{survive (te , t1 ] | alive at te } ×
P{survive (t1 , t2 ] | alive at t1 } ×
P{survive (t2 , tx ] | alive at t2 } ×
P{d at tx | alive just before tx }(1.1)

By choosing more intermediate time points we can make the intervals arbitrarily small, so wejust need to derive the probability of surviving a small piece of time, as a function of themortality rate.
For a start assume that the mortality is constant over time λ(t) = λ. From the denition of
a rate we have (conditional on being alive at t):

P{death during (t, t + h]} ≈ λh
⇒ P{survive (t, t + h]} ≈ 1 − λh(1.2)

where the approximation gets better the smaller h is. Suppose we have survival for a timespan y = tx − te and that this is subdivided in N intervals, each of length h = y/N , then thesurvival probability for the entire span from te to tx is the product of probabilities of surviving
each of the small intervals, conditional on being alive at the beginning of each interval:

N
λyP{survive te to tx } ≈ (1 − λh) = 1 −NNnFrom mathematics it is known that (1 + x/n) → exp(x) as n → ∞ (some dene exp(x) thisway). So if we divide the time span y in successively smaller pieces we will have that N → ∞,and hence that:
P{survive te to tx } ≈

λy
1−
N

N
→ exp(−λy),

N →∞

(1.3)

Therefore the contribution to the likelihood from a person observed for a time span of lengthy is exp(−λy), and the contribution to the log-likelihood is therefore −λy .
If we observe a person dying at the end of the last interval, the contribution to thelikelihood from the last interval will be the probability surviving till just before the end of theinterval, multiplied by the probability of dying in the last tiny instant (of length ϵ, say) of theinterval (the last term in (1.1)). This probability is by (1.2) λϵ, and hence the log-likelihood
contribution from this last instant is log(λϵ) = log(λ) + log(ϵ).
The total likelihood for one person is the product of all these terms from the follow-up
intervals (i) for the person; and the log-likelihood (ℓ) is therefore:ℓ(λ) = −λ

X
i

=

X
i

yi +

X

di log(λ) +iXdi log(ϵ)

i Xdi log(λ) − λyi +di log(ϵ)
i

4

1.3

The Cox-likelihood as a prole likelihood

WmtCmathwhere yi is observation time (risk time or person-years) in the iindicator (0/1) for the ithinterval and di the deathinterval (which is 0 for all intervals, except possibly for the last).

We want to estimate λ, so terms that does not involve λ can be ignoredsuch as the lastterm. Very convenient, we then do not have to decide precisely how small ϵ is.
Thus we have established that the contribution to the log-likelihood from a single person's
follow-up is the sum of a number of terms of the form di log(λ) − λyi where di is the eventth

indicator (0/1) and yi the length of the iinterval.

If we subdivide the follow-up of a person across several records, each representing yi offollow-up, and keep track of the deaths in each interval, di (which is 0 for all intervals, except
possibly for the last if the person dies), then each record (di , yi ) will represent a contribution
to the log-likelihood of di log(λ) − λyi . This will be exploited in the practical modeling ofrates.
By the assumption that the rate λ is constant over time, the log-likelihood contribution
from a person with d deaths (0 or 1) at the end of a follow-up period of y is d log(λ) − λy .
But once we have subdivided follow-up we do not need the assumption of a constant rate
across the intervals, we can allow separate rates (λi s) in dierent intervals, relaxing theassumption to only constant rates within each of the very small intervals. These rates can ofcourse not be estimated based on the observation from a single person. A model for the λi s
must be specied in terms of covariates associated with each interval. Among these will
normally be the value of one or more time scales at the beginning of each interval.

1.2.2

Poisson likelihood

The Poisson distribution with mean µ is a distribution on the non-negative integers(x = 0, 1, 2, 3, . . .), where:

P{X = x} =µx exp(−µ)
x!

Thus, the log-likelihood from observation of a Poisson variate x with mean µ is

x log(µ) − µ − log(x!). The last term does not depend on the parameter µ so it can beomitted when maximizing the log-likelihood over values of µ.
Changing the notation, the log-likelihood contribution from a Poisson-variate d with mean
λy is d log(λy) − λy = d log(λ) − λy + d log(y), which is the same as the likelihood for d events
during y follow-up time with rate λ, except for the term d log(y). But this term does notdepend on the parameter λ and therefore can be ignored when maximizing the log-likelihood.
This means that maximum-likelihood estimation for observations from follow up, (d, y), canbe done using a function that can maximize the likelihood for independent Poisson variables.
We just need to pretend that the ds are independent Poisson variates with means λy . But
recall that the ds are neither independent nor Poisson distributed.
In R we have the functions glm and gam to do this, depending on how we will model thecovariate eects.

1.3 The Cox-likelihood as a prole likelihood
The Cox model [2] species the intensity (rate) as a function of time (t) and the covariates
through the linear predictor ηi = β1 x1i + · · · + βp xpi as:λ(t, xi ) = λ0 (t) exp(ηi )1.3

The Cox-likelihood as a prole likelihood

Theory5leaving the baseline hazard λ0 unspecied.
Cox devised thepredictor

partial(log-)likelihood for the parameters β = (β1 , . . . , βp ) in the linearXℓ(β) =

log

death timeseηdeathP
ηii∈Rt e

where Rt is the risk set at time t, i.e. the set of individuals at risk at time t.
Suppose the time-scale has been divided into small time intervals with at most one death ineach, and that we in addition to the regression parameters describing the eect of covariates
use one parameter per time interval to describe the eect of time (i.e. the chosen timescale).
Thus the model with constant rates in each small interval is:
log λ(t, xi ) = log λ0 (t) + β1 x1i + · · · + βp xpi = αt + ηiusing αt = log λ0 (t) Assume w.l.o.g. the y for these empirical rates are 1. The log-likelihood
contributions that contain information on a specic time-scale parameter αt , relating to aparticular time t, will be contributions from the empirical rate (d, y) = (1, 1) with the deathat time t, and the empirical rates (d, y) = (0, 1) from all other individuals at risk at time t.
Note that there is exactly one contribution from each individual at risk to this part of the
log-likelihood:ℓt (αt , β) =XXdi (αt + ηi ) − eαt +ηi = αt + ηdeath − eαt
eηi
i∈Rt

i∈Rt

where ηdeath is the linear predictor for the individual that died at t. For those intervals on the

1time-scale where no deaths occur the estimate of the αt will be −∞ , and so these intervals
will not contribute to the log-likelihood.
The derivative w.r.t. αt is:

Dαt ℓ(αt , β) = 1 − eαtXeηi = 0⇒αt = Pec

1i∈Rt ei∈Rt

ηiαIf this estimate of e t is fed back into the log-likelihood for αt , we get the(with αt proled out):


log



1
P

i∈Rt eηi+ ηdeath − 1 = logeηdeathP
ηii∈Rt eprole likelihood

−1which is the same as the contribution from time t to Cox's partial likelihood.
Thus we may estimate the regression parameters from the Cox model by standard
Poisson-regression software by splitting the data nely and specifying the model as having one
rate parameter per time interval.
The Cox model could therefore have been formulated as model with a baseline rate modeledby a timescale parameter for each time recorded. This is an exchangeable model for thebaseline rate parameters, thus using neither the ordering nor the absolute scaling of the times.
The results will be the same, also for the s.e.s, since everything is derived . This is
illustrated in section 2.1, where fully parametric alternatives to the Cox model is describedtoo.

1 This is because the term α +ηtdeath vanishes if all di = 0, and the last term is maximal if eαt

= 0 ⇔ αt = −∞

6

1.4

Practical data processing

WmtCma

1.4 Practical data processing
Implementation of the Poisson-approach in practice requires that follow-up for each individualis split in small pieces of follow-up along one or more time scales. The relevant time-varying
covariates should be computed for each interval and xed covariates should be carried over toall intervals for a given individual.
Presently there are (at least the following) tools for this in:

Stata: The function stsplit is part of standard Stata, it is a descendant of stlexis written
by Michael Hills & David Clayton.

SAS: A macro %Lexis, available at http://BendixCarstensen/Lexis, written by Bendix
Carstensen. Another macro is by Klaus Rostgaard [5]https://sourceforge.net/p/pyrsstep/wiki/Home/.
R:

Function survSplit from the survival package does the job. The Epi package has afunction splitLexis that does this for Lexis objects [4, 1], and in the popEpi package

there is a faster data.table based version, splitMulti, which also has a more friendly
syntax.

These tools expand a traditional survival dataset with one record per individual to one withseveral records per individual, one record per follow-up interval. In the following we shallrestrict attention to the Lexis tools in R.

time scale

The split data makes a clear distinction between
interval andrisk time

which is the length of each

which is the value of the timescale at (the beginning of ) each interval,be that time since entry, current age, calendar time, etc.
In traditional Poisson modeling the log-risk time is used as oset and the time is used ascovariate. Thus Poisson modeling of follow-up data makes a clear distinction between risktime as the response variable and time scale(s) as covariate(s). A recent addition to the Epi

2package is the family poisreg , which uses the more intuitive specication of the response asa two-column vector of events and person-years.

1.4.1

Estimation of baseline hazard

Once data has been split in little pieces of follow-up time, the eect of any timescale can be
estimated using parametric regression tools such as splines. This will directly produce
estimated baseline rates by using standard prediction machinery for generalized linear models
with a given set of covariates.
Suppose h(t) is a parametric function which is parametrized linearly by the parameters inτ , h(t) = w′ τ (w and τ are column vectors). The model can be formulated as:
 

τ
′
′
′
′
log λ(t, x) = h(t) + x γ = w τ + x γ = (w x)γStandard prediction machinery can be used to produce estimates of log-rates with standard
errors for a set of values of t (and hence w ), and some chosen values of the variables in x.
This is a standard tool in any statistical package able to t generalized linear models. Rate

2 This means that in a glm or gam model you can specify family=poisreg

1.4

Practical data processing

Theory7estimates with condence intervals are then derived by taking the exponential function of theestimates for the log-rates with condence intervals.
In the Epi package this is handled by the ci.pred function that produces predicted rates
for a specied set of prediction points.

1.4.2

Estimation of survival function

In studies where entry time 0 is meaningful the survival function is a simple, albeit non-linear
function of the rates:

 Z t
S(t) = exp − λ(s) ds
0

so in order to estimate this from a parametric model for the log-rates we need to derive theintegral, i.e. a cumulative sum of predictions. If we want standard errors for this we musthave not only standard errors for the λs, but the entire the variance-covariance matrix of
estimated values of λ.
From a generalized linear model we can easily extract estimates for log

λ(t)



at any set ofpoints. This is just a linear function of the parameters, and so the variance-covariance matrixof these can be computed from the variance-covariance matrix of the parameters.
A Taylor approximation of the variance-covariance matrix for λ(t) can be obtained fromthis by using the derivative of the function that maps log

λ(t)



to λ(t). This is thecoordinate-wise exponential function, so the matrix required is the diagonal matrix withentries log


λ(t) .

Finally the cumulative sum is obtained by multiplying with a matrix with 1s on and below
the diagonal, so this matrix just needs to be pre and post-multiplied in order to produce thevariance-covariance of the cumulative hazard at the prespecied points.
In technical terms we let fˆ(ti ) be estimates for the log-rates for a certain set of covariate
values (x) at points ti , i = 1, . . . , I , derived by:fˆ(ti ) = B β̂where β = (τ, γ) is the parameter vector in the model (of length p, say), including theparameters that describe the baseline hazard. Here, fˆ(ti ) is I × 1, B is I × p and β̂ is p × 1.
Now, let the estimated variance-covariance matrix of β be Σ, a p × p matrix. Then the′variance-covariance of fˆ(ti ) is BΣB which is I × I . The transformation to the rates is thecoordinate-wise exponential function so the derivative of this (evaluated at the m.l.e.) is the
diagonal matrix with entries exp

fˆ(ti ) , so the variance-covariance matrix of the rates at the

points ti isˆˆdiag(ef (ti ) ) B Σ B′ diag(ef (ti ) )′Finally, the transformation to the cumulative hazard (assuming that all intervals have lengthℓ) is by an I × I matrix of the form:


1
 1

L=ℓ×
 1
 1
1

0
1
1
1
1

0
0
1
1
1

0
0
0
1
1


0
0 

0 

0 
1

8

1.4

Practical data processing

WmtCmaso the (approximate) variance-covariance matrix for the cumulative hazard is:ˆˆ

L diag(ef (ti ) ) B Σ B′ diag(ef (ti ) )′ L′
This formula for the variance of the cumulative hazard does not guarantee that the derived
condence intervals' lower endpoints are larger than 0. This can be xed by computing
condence intervals for the log-cumulative hazard using the delta-rule, and back-transforming
to the rate scale.
These calculations are implemented in the Epi package function ci.cum, which requires (atleast) 3 objects as arguments: 1) a model object representing a multiplicative model for
occurrence rates, 2) a prediction data frame which will produce (log) rate-estimates from themodel at a set of equidistant times since some origin, and 3) a scalar representing the distance
between the prediction times (in the units in which the person-years was supplied to themodel). The function also has a facility for computing the condence limits on the
log-cumulative hazard scale and back transforming to ensure positive lower condence bounds
for the integrated hazard.
Once we have estimated the cumulative hazard function as a function of time we cantransform it to the survival function by the exponentialthis is implemented in the function

ci.surv.

Chapter 2
Examples

2.1 Equality of Cox and Poisson modeling: The lungcancer example
In this section we use the lung cancer example data from the survival package to illustratethat the results from a Cox model are identical to results from a particular Poissonmodelalbeit quite an absurd one. We also illustrate two ways to use a parametricallysmoothed version of the linear predictor in the Poisson model to obtain a sane estimate forthe baseline hazard.

> library(Epi)> library(popEpi)> library(survival)> library(mgcv)> data(lung)> str(lung)'data.frame':
228 obs. of 10 variables:$ inst: num 3 3 3 5 1 12 7 11 1 7 ...
$ time: num 306 455 1010 210 883 ...
$ status : num 2 2 1 2 2 1 2 2 2 2 ...
$ age: num 74 68 56 57 60 74 68 71 53 61 ...
$ sex: num 1 1 1 1 1 1 2 2 1 1 ...
$ ph.ecog : num 1 0 0 1 0 1 2 2 1 2 ...
$ ph.karno : num 90 90 90 90 100 50 70 60 70 70 ...
$ pat.karno: num 100 90 90 60 90 80 60 80 80 70 ...
$ meal.cal : num 1175 1225 NA 1150 NA ...
$ wt.loss : num NA 15 15 11 0 0 10 1 16 34 ...
> lung[1:10, ]
1
2
3
4
5
6
7
8
9
10

inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
3 306
2 74 1
1
90
100
1175
NA
3 455
2 68 1
0
90
90
1225
15
3 10101 56 1
0
90
90
NA
15
5 210
2 57 1
1
90
60
1150
11
1 883
2 60 1
0
100
90
NA
0
12 1022
1 74
1
1
50
80
513
0
7 310
2 68 2
2
70
60
384
10
11 361
2 71 2
2
60
80
538
1
1 218
2 53 1
1
70
80
825
16
7 166
2 61 1
2
70
70
271
34

9

10

2.1

Equality of Cox and Poisson modeling: The lung cancer example

WmtCmaHow many deaths are there?

> table(lung$status)
1 2
63 165
Convert sex to a factor:> lung$sex <- factor(lung$sex, labels = c("M", "F"))How many distinct event times?

> addmargins(table(table(lung$time)))1 2 3 Sum146 382 186
To avoid tied event times we add a small random quantity to each time:> set.seed(1952)> lung$time <- lung$time + runif(lung$time, -2, 2)> addmargins(table(table(lung$time)))1 Sum
228 228First we t a traditional Cox-model for the Mayo clinic data

> m0.cox <- coxph(Surv(time, status == 2) ~ age + sex, data = lung)> summary(m0.cox)Call:coxph(formula = Surv(time, status == 2) ~ age + sex, data = lung)n= 228, number of events= 165
coef exp(coef) se(coef)z Pr(>|z|)age 0.016985 1.017130 0.009222 1.842 0.06550 .
sexF -0.518435 0.595452 0.167496 -3.095 0.00197 **--Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
age
sexF

exp(coef) exp(-coef) lower .95 upper .95
1.0171
0.98320.99891.03570.59551.67940.42880.8268

Concordance= 0.603 (se = 0.025 )Likelihood ratio test= 14.3 on 2 df, p=8e-04
Wald test= 13.64 on 2 df, p=0.001
Score (logrank) test = 13.9 on 2 df, p=0.001
Now create a Lexis object from the dataset lung, representing the follow-up time and events:> Lung <- Lexis(exit = list(tfe = time),+
exit.status = factor(status,+
labels = c("Alive", "Dead")),+data = lung)
NOTE: entry.status has been set to "Alive" for all.
NOTE: entry is assumed to be 0 on the tfe timescale.

2.1

Equality of Cox and Poisson modeling: The lung cancer example

Examples11> summary(Lung)
Transitions:To
From
Alive Dead
Alive
63 165

Records: Events: Risk time: Persons:228
165 69618.73
228

Split data in small intervals, dened by all recorded events and censoring times> Lung.s <- splitMulti(Lung, tfe = c(0, sort(unique(Lung$time))))> summary(Lung.s)
Transitions:To
From
Alive Dead Records: Events: Risk time: Persons:
Alive 25941 165
26106
165 69618.73
228> Lung.s[1:10, 1:10]
lex.id tfe lex.dur lex.Cst lex.Xst inst
time status age sex
1 0.006.78 Alive Alive
3 307.804
2 74 M
1 6.783.25 Alive Alive
3 307.804
2 74 M
1 10.04
0.15
Alive
Alive
3 307.804
2 74 M
1 10.19
0.38
Alive
Alive
3 307.804
2 74 M
1 10.57
1.83
Alive
Alive
3 307.804
2 74 M
1 12.40
0.63
Alive
Alive
3 307.804
2 74 M
1 13.03
1.62
Alive
Alive
3 307.804
2 74 M
1 14.65
1.01
Alive
Alive
3 307.804
2 74 M
1 15.66
8.94
Alive
Alive
3 307.804
2 74 M
1 24.59
5.02
Alive
Alive
3 307.804
2 74 M
Now t the Cox model to the Lexis data set as well as the time-split Lexis data set; note thecode is exactly the same, only the data= argument diers:> mL.cox <- coxph(Surv(tfe, tfe+lex.dur, lex.Xst == "Dead") ~ age + sex,+
eps = 10^-11, iter.max = 25, data = Lung)> mLs.cox <- coxph(Surv(tfe, tfe+lex.dur, lex.Xst == "Dead") ~ age + sex,+
eps = 10^-11, iter.max = 25, data = Lung.s)> round(rbind(ci.exp( m0.cox),+
ci.exp( mL.cox),+
ci.exp(mLs.cox))[c(1,3,5,2,4,6),], 6)exp(Est.)
2.5%
97.5%
age 1.017130 0.998911 1.035681
age 1.017130 0.998911 1.035681
age 1.017130 0.998911 1.035681
sexF 0.595452 0.428818 0.826837
sexF 0.595452 0.428818 0.826837
sexF 0.595452 0.428818 0.826837
We see we get the same results from the three dierent sets of data  they contain exactly
the same amount of information about the rates.
Now we t the corresponding Poisson model with factor modeling of the time scale  notethat we use the poisreg family where we enter events and person-years as a 2 column matrix:> nlevels(factor(Lung.s$tfe))[1] 228122.1

Equality of Cox and Poisson modeling: The lung cancer example

WmtCma> system.time(+ mLs.pois.fc <- glm(cbind(lex.Xst == "Dead", lex.dur) ~ 0 + factor(tfe) + age + sex,+family = poisreg, data = Lung.s))user system elapsed
15.39
0.47 15.86> length(coef(mLs.pois.fc))[1] 230> rbind(ci.exp(mLs.cox),+
ci.exp(mLs.pois.fc, subset = c("age", "sex")))[c(1,3,2,4),]exp(Est.)
2.5%
97.5%
age 1.0171302 0.9989114 1.0356813
age 1.0171302 0.9989114 1.0356813
sexF 0.5954519 0.4288185 0.8268369
sexF 0.5954519 0.4288185 0.8268369
In accordance with the mathematical derivations above, we see that the estimates of theregression coecients are exactly the same from the Cox model and the Poisson model. The
latter has an extra 228 parameters estimated, which is what causes the very long estimation
time.

2.1.1

Parametric baseline

To get a more realistic model for the baseline rate we now dene knots for a spline basis andt the model with natural splines for the baseline eect of tfe. These knots are basicallytaken out of thin air:> t.kn <- c(0, 25, 100, 500, 1000)> system.time(+ mLs.pois.sp <- glm(cbind(lex.Xst == "Dead", lex.dur)+~ Ns(tfe, knots = t.kn) + age + sex,+family = poisreg,+
data = Lung.s))user system elapsed
0.09
0.04
0.12
Finally we t the model with a penalized spline model for the eect of tfe using gam from the

mgcv package:> system.time(+ mLs.pois.ps <- gam(cbind(lex.Xst == "Dead", lex.dur)+~ s(tfe) + age + sex,+family = poisreg,+
data = Lung.s))user system elapsed
0.66
0.39
1.05
> summary(mLs.pois.ps)
Family: poisson
Link function: log
Formula:cbind(lex.Xst == "Dead", lex.dur) ~ s(tfe) + age + sex

2.1

Equality of Cox and Poisson modeling: The lung cancer example

Examples13

Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)(Intercept) -6.94486
0.59468 -11.678 < 2e-16 ***age
0.016290.00920 1.771 0.07663 .
sexF-0.50706
0.16731 -3.031 0.00244 **--Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df Chi.sq p-values(tfe) 2.131 2.687 17.74 0.000649 ***--Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 9.74e-06 Deviance explained = 1.7%
UBRE = -0.93042 Scale est. = 1
n = 26106
We see that the eective d.f. for the time scale eect (tfe) is about 2, so some indication thatthe arbitrary spline may be over-modeling data.
Finally we make an overall comparison of estimates of age and sex eects from the dierent
approaches:> ests <+ rbind(ci.exp(m0.cox),+
ci.exp(mLs.cox),+
ci.exp(mLs.pois.fc, subset = c("age", "sex")),+ci.exp(mLs.pois.sp, subset = c("age", "sex")),+ci.exp(mLs.pois.ps, subset = c("age", "sex")))> cmp <- cbind(ests[c(1, 3, 5, 7, 9) , ],+ests[c(1, 3, 5, 7, 9)+1, ])> rownames(cmp) <+ c("Cox", "Cox-split", "Poisson-factor", "Poisson-spline", "Poisson-Pspline")> colnames(cmp)[c(1, 4)] <- c("age", "sex")> round(cmp, 7)
age
2.5%
97.5%
sex
2.5%
97.5%
Cox
1.017130 0.9989114 1.035681 0.5954519 0.4288185 0.8268369
Cox-split
1.017130 0.9989114 1.035681 0.5954519 0.4288185 0.8268369
Poisson-factor 1.017130 0.9989114 1.035681 0.5954519 0.4288185 0.8268369
Poisson-spline 1.016214 0.9980575 1.034701 0.5994793 0.4317352 0.8323977
Poisson-Pspline 1.016423 0.9982591 1.034917 0.6022668 0.4338818 0.8360001
So even if the factor model, and by that token also the Cox-model, seem pretty far fetched intheir (lack of ) assumptions, there is minimal dierence to the regression parameter estimatesfrom the models with more realistic assumptions for the baseline rates.

2.1.2

Rates, cumulative rates and survival

Parametric models
Now we compute the estimated rates and cumulative rates over 10-day periods for 60 year oldmen, and then the survival function at these points.
In order to get the predictions from the spline model we specify a prediction data frame,
where we predict rates at equidistant points, using ci.pred for the rates. Note that since weused the poisreg family, the predicted rates are by denition per one unit of lex.dur, whichin our case is days, so we multiply by 365.25 to get rates per 1 PY.

14

2.1

Equality of Cox and Poisson modeling: The lung cancer example

WmtCmaFurther, when we compute the survival function (and by that token also the cumulative
rates), we must supply the interval length (distance between values of tfe in the prediction

data frame) by intl= in the units of lex.dur. If we do not supply it, it will be taken as thedierence between the two rst elements of the rst column in the prediction data frame:> # midpoints of 10-day intervals and other covariates> nd <- data.frame(tfe = seq(5, 995, 10), age = 60, sex = "M")> #> # the rates from the spline model (events/year)> lambda <- ci.pred(mLs.pois.sp, nd) * 365.25> # the survival function> survP <- ci.surv(mLs.pois.sp, nd)
NOTE: interval length chosen from as tfe[2] - tfe[1]> # same same for the penalized spline model:> lambdap <- ci.pred(mLs.pois.ps, nd) * 365.25> survPp <- ci.surv(mLs.pois.ps, nd)
NOTE: interval length chosen from as tfe[2] - tfe[1]So now we have the incidence rates per 1 PY as well as cumulative incidence rates and thecorresponding survival function(s) based both on natural splines and a penalized likelihood
via gam.

Cox models
The Breslow-estimator of the survival curve from the corresponding Cox-model for a maleaged 60 is obtained from the m0.cox object:> sf <- survfit(m0.cox, newdata = data.frame(sex = "M", age = 60))
We can extract the baseline rates from the corresponding Poisson model as well. Since

lex.dur is supplied in units of days to mLs.pois.fc, the predicted rates from using ci.exp
on the model will be in events per day, hence we rescale to events per year. We extract thetimes from the names of the parameters:> (nc <- length(coef(mLs.pois.fc)))[1] 230> br <- ci.exp(mLs.pois.fc,+
ctr.mat = cbind(diag(nc - 2), 60, 0)) * 365.25> bt <- as.numeric(gsub("factor\\(tfe)", "", names(coef(mLs.pois.fc))[1:(nc - 2)]))> head(cbind(bt, br))bt exp(Est.)
2.5%
97.5%
[1,] 0.00000 0.2653605 0.03724572 1.890585[2,] 6.78194 0.5547910 0.07787287 3.952507[3,] 10.03619 11.8535350 1.66384425 84.446782[4,] 10.18946 4.8591849 0.68210629 34.615834[5,] 10.56603 1.0050413 0.14108059 7.159794[6,] 12.39700 2.9190141 0.40975739 20.794362
Now we have the predicted rates in intervals between the times observed; the Poisson model

end

left

tted implicitly assumes that event rates are constant within intervals between times. Since
the deaths occur at theof the intervals, and intervals are named by theirplotting of the rates must use typewhere the curverst

endpoint,= "s", which creates steps between successive points

moves horizontally, then vertically.

2.1

Equality of Cox and Poisson modeling: The lung cancer example

Examples15

Natural spline vs. penalized splines
First we just compare the two smooth curves:> par(mfrow = c(1, 2), mar = c(3, 3, 1, 1), mgp = c(3, 1, 0)/1.6,+
bty = "n", las = 1, lend = "butt")> matshade(nd$tfe, cbind(lambda, lambdap), plot = TRUE,+
col = c("blue", "red"), lwd = 3, lty = c("solid", "21"),+xlim = c(0, 900), xaxs = "i", ylim = c(1/5, 20), log = "y",+xlab = "Days since diagnosis",+ylab = "Mortality rate per 1 year")> matshade(nd$tfe-5, cbind(survP[, -4], survPp[, -4]), plot = TRUE,+
col = c("blue", "red"), lwd = 3, , lty = c("solid", "21"),+xlim = c(0, 900), xaxs = "i", yaxs = "i", ylim = 0:1,+
xlab = "Days since diagnosis",+ylab = "Survival probability")

1.0
20.0

0.8

5.0
Survival probability

Mortality rate per 1 year

10.0

2.0

1.0

0.5

0.6

0.4

0.2

0.2
0.0
0

200

400
600
Days since diagnosis

800

0

200

400
600
Days since diagnosis

800

Left panel: Estimated mortality rates by Poisson models; blue is the natural splinewith pre-chosen knots, the red is the gam model with penalization. Right panel: The resulting
survival curves. Shaded areas indicate 95% condence intervals.
./lung-rtsurv-sm
Figure 2.1:

Comparison with Cox
The we make the same plots overlaid with the results from the Cox model:

16

2.1

Equality of Cox and Poisson modeling: The lung cancer example

WmtCma> par(mfrow = c(1, 2), mar = c(3, 3, 1, 1), mgp = c(3, 1, 0)/1.6,+
bty = "n", las = 1, lend = "butt")> plot(NA, xlim = c(0, 900), xaxs = "i", ylim = c(1/5, 20), log = "y",+xlab = "Days since diagnosis",+ylab = "Mortality rate per 1 year")> lines(bt, br[, 1], type = "s", col = gray(0.6))> matshade(nd$tfe, cbind(lambda, lambdap), # plot = TRUE,+
col = c("blue", "red"), lwd = 3, lty = c("solid", "21"))> matshade(nd$tfe, cbind(survP[, -4], survPp[, -4]), plot = TRUE,+
col = c("blue", "red"), lwd = 3, , lty = c("solid", "21"),+xlim = c(0, 900), xaxs = "i", yaxs = "i", ylim = 0:1,+
xlab = "Days since diagnosis",+ylab = "Survival probability")> lines(sf, lwd = 1, lty = c(1, 1))> lines(sf, lwd = 2, conf.int = FALSE)
1.0
20.0

0.8

5.0
Survival probability

Mortality rate per 1 year

10.0

2.0

1.0

0.5

0.6

0.4

0.2

0.2
0.0
0

200

400
600
Days since diagnosis

800

0

200

400
600
Days since diagnosis

800

Left panel: Estimated mortality rates by Poisson models; blue is the natural splinewith pre-chosen knots, the red is the gam model with penalization, and the thin gray line indicate the estimated baseline hazard from the Poisson (Cox) model with one parameter per
event/censoring time. Right panel: The resulting survival curves, over-laid in black with theBreslow-estimator of the survival curve. Shaded areas and thin lines indicate 95% condence
intervals.
./lung-rtsurv-cp
Figure 2.2:

From gure 2.1 we see that there is only slight dierence between the two parametricapproaches; the penalized splines (red broken curve) smooths a bit more hat the natural
splines with arbitrarily chosen knots. When transformed to the survival scale, the two
approaches are practically indistinguishable.

2.1

Equality of Cox and Poisson modeling: The lung cancer example

Examples17

Figure 2.1 has the Cox-model estimates overlaid; strictly speaking the baseline hazard is
not really part of the Cox-model, the underlying hazard comes from the corresponding
Poisson model, the survival curve is the Breslow estimator. There is a faint indication thatthe parametric curves produces slightly narrower condence bands for the survival
probabilities than the Breslow-estimator.

2.1.3

Practical splitting

In practical applications the splitting of time need not be at the times of events andcensorings; this was only done above to demonstrate the connection between the Cox modeland the Poisson model.
The assumption behind the Poisson approach is essentially only the assumption that amodel with constant rates in each small interval gives an adequate description of data. So inpractice we would split data in small equidistant intervals. In the lung cancer dataset there
are 165 deaths and the total observation period is some 1000 days, some 2.8 years, so we split
the follow-up in intervals of 20 days:> sL <- splitMulti(Lung, tfe = seq(0, 1100, 20))> summary(Lung)
Transitions:To
From
Alive Dead Records: Events: Risk time: Persons:Alive
63 165228
165 69618.73
228> summary(sL)Transitions:To
From
Alive Dead Records: Events: Risk time: Persons:
Alive 3443 1653608
16569618.73
228
so we have much fewer records but the same number of events and person-time.
We can then compare to the estimates from the parametric model mLs.pois.sp, if we
instead use the equidistantly cut dataset:> mLs.pois.se <- update(mLs.pois.sp, data = sL)> ee <- cbind(ci.exp(mLs.pois.sp),+
ci.exp(mLs.pois.se),+ci.exp(mLs.pois.sp)/+
ci.exp(mLs.pois.se))> colnames(ee)[c(1,4,7)] <- "expEst"> round(ee, 3)
expEst 2.5% 97.5% expEst 2.5% 97.5% expEst 2.5% 97.5%(Intercept)
0.001 0.000 0.002 0.001 0.000 0.002 0.869 0.774 0.976
Ns(tfe, knots = t.kn)1 2.773 1.037 7.415 2.639 1.103 6.312 1.051 0.940 1.175
Ns(tfe, knots = t.kn)2 2.817 0.902 8.797 2.454 0.880 6.841 1.148 1.025 1.286
Ns(tfe, knots = t.kn)3 3.836 0.536 27.479 2.894 0.586 14.296 1.325 0.914 1.922
Ns(tfe, knots = t.kn)4 3.253 0.705 15.010 3.329 0.792 13.999 0.977 0.891 1.072
age
1.016 0.998 1.035 1.016 0.998 1.035 1.000 1.000 1.000
sexF
0.599 0.432 0.832 0.599 0.432 0.832 1.000 1.000 1.000
We see only minor dierences in the estimated values of the regression parameters, while itappears that the spline parameters are somewhat dierent. This does however not translate
any relevant dierences in the estimated curves:

18

2.2

Stratied models

WmtCma> plot(NA, xlim = c(0, 900), xaxs = "i", ylim = c(1/5, 20), log = "y",+xlab = "Days since diagnosis",+ylab = "Mortality rate per 1 year")> matshade(nd$tfe, cbind(ci.pred(mLs.pois.sp, nd),+ci.pred(mLs.pois.se, nd)) * 365.25,+
lwd = 2, col = c('blue', 'black'), log = "y" )20.0

Mortality rate per 1 year

10.0

5.0

2.0

1.0

0.5

0.2
0

200

400
600
Days since diagnosis

800

Comparing the same model tted to data split at all 228 recorded event and censoringtimes (26,106 records) (blue), and tted to a data set only cut every 20 days (3,591 records)(black).
./lung-spcmp
Figure 2.3:

Thus from gure 2.3 it appears that the splitting of the follow-up time in 20-day intervals issucient to render the estimation of the baseline hazard reliable.

2.2 Stratied models
A stratied Cox-model is a model where the underlying hazard is allowed to dier between
strata, i.e. between levels of a categorical variable.
Thus this is merely an interaction between time and a categorical variable. If a spline basis
has been chosen as model for the time variable a model with separate baseline hazards foreach level of a factor F this is easily modelled by saying:
For illustration we use the lung cancer example again, the Lexis object where time is split
every 20 days:

2.2

Stratied models

Examples19

NOTE: entry.status has been set to "Alive" for all.
NOTE: entry is assumed to be 0 on the tfe timescale.
> summary(sL)Transitions:To
From
Alive Dead Records: Events: Risk time: Persons:
Alive 3432 165
3597
165
69593
228
For the modeling of the baseline rate (timescale tfe) we dene the knots and t a natural
spline, one with main eect of sex, the other with an interaction. Note that there is norequirement that the time-part of the interaction is parametrized in the same way as the maineect. The model m3 below uses a simpler time-eect in the interaction:> kn <- c(0, 50, 150, 450)
> m1 <- glm(cbind(lex.Xst=="Dead", lex.dur)+~ Ns(tfe, knots = kn) + sex + age,+family = poisreg,+
data = sL)> m2 <- glm(cbind(lex.Xst=="Dead", lex.dur)+~ Ns(tfe, knots=kn) * sex + age,+family = poisreg,+
data = sL )> m3 <- update(m1, . ~ . + Ns(tfe, knots = c(0, 50, 200)):sex )> anova(m3, m1, m2, test = "Chisq")
Analysis of Deviance Table
Model 1: cbind(lex.Xst == "Dead", lex.dur) ~ Ns(tfe, knots = kn) + sex +age + sex:Ns(tfe, knots = c(0, 50, 200))
Model 2: cbind(lex.Xst == "Dead", lex.dur) ~ Ns(tfe, knots = kn) + sex +ageModel 3: cbind(lex.Xst == "Dead", lex.dur) ~ Ns(tfe, knots = kn) * sex +age
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1
3588
1293.2
23591
1298.6 -3 -5.3226 0.1496
33588
1294.2 34.32920.2280
There is no signicant interaction here, but the test statistic would actually have beensignicant if the interaction were only on two degrees of freedom. Hence we may want toinspect the two tted baseline rates, as well as their ratio:> par(mfrow = c(1, 2), mar = c(3, 3, 1, 3), mgp = c(3, 1, 0), las = 1)
> nm <- data.frame(tfe = seq(0, 1000, 10),+age = 65,+
sex = factor("M", levels = c("M", "F")))> nf <- data.frame(tfe = seq(0, 1000, 10),+age = 65,+
sex = factor("F", levels = c("M", "F")))> plot(NA, xlim = c(0,900), xaxs = "i", ylim = c(1/100, 5), log = "y",+xlab = "Days since diagnosis",+ylab = "Mortality rate per 1 year" )> matshade(nm$tfe, cbind(ci.pred(m2, nm) * 365.25,+
ci.pred(m2, nf) * 365.25,+
ci.exp (m2, list(nm,nf)) / 20),

20

2.2

Stratied models

WmtCma

5.00

5.00

2.00

2.00

1.00

1.00

0.50

0.20
0.10

2
1.5

0.05

1

Mortality rate per 1 year

Mortality rate per 1 year+
lwd = 2, col = c('blue', 'red', 'black'))> abline(h = 1/20, lty = 3)
> axis(side = 4, at = c(2, 5, 10, 15, 20) / 200, labels = c(2, 5, 10, 15, 20) / 10 )> axis(side = 4, at = c(2:9)/200, labels = NA, tcl = -0.3 )> plot(NA, xlim = c(0,900), xaxs = "i", ylim = c(1/100,5), log = "y",+xlab = "Days since diagnosis",+ylab = "Mortality rate per 1 year" )> matshade(nm$tfe, cbind( ci.pred(m3,nm)*365.25,+
ci.pred(m3,nf)*365.25,+
ci.exp (m3,list(nm,nf))/20 ),+lwd = 2, col = c('blue','red','black') )> abline(h = 1 / 20, lty = 3)
> axis(side = 4, at = c(2, 5, 10, 15, 20) / 200,+
labels = c(2, 5, 10, 15, 20) / 10)> axis(side = 4, at = c(2:9) / 200, labels = NA, tcl = -0.3)

0.50

0.20
0.10

2
1.5

0.05

1

0.5
0.02

0.5
0.02

0.01

0.2
0

200

400

600

800

0.01

0.2
0

200

400

600

800

Baseline rates for 65 year old men (blue) resp. women (red), and the rate-ratio
between these (black). The leftmost panel uses the same set of knots for the main eect and theinteraction, the rightmost a more parsimonious interaction specication.
./strat-prcmp
Days since diagnosis

Days since diagnosis

Figure 2.4:

From gure 2.4, it is clear that although there is no formal interaction (p-values are
15-20%), there is a clear tendency that the mortality among men is higher during the rstyear or so after diagnosis.
For illustration we repeat the same exercise with the gam machinery. Note that theinteraction specication s(tfe,bymust be specied:= sex) does not contain the main eect of sex, so this2.3

Time-varying coecients

Examples21> p1 <- gam(cbind(lex.Xst == "Dead", lex.dur) ~ s(tfe)+ sex + age,+family = poisreg,+
data = sL)> p2 <- gam(cbind(lex.Xst == "Dead", lex.dur) ~ s(tfe, by = sex) + sex + age,+family = poisreg,+
data = sL)> anova(p1, p2, test = "Chisq")
Analysis of Deviance Table
Model 1: cbind(lex.Xst == "Dead", lex.dur) ~ s(tfe) + sex + age
Model 2: cbind(lex.Xst == "Dead", lex.dur) ~ s(tfe, by = sex) + sex +age
Resid. Df Resid. Dev
Df Deviance Pr(>Chi)
13591.41299.523590.91298.0 0.52822 1.4688 0.1092> par(mar = c(3,3,1,3), mgp = c(3,1,0), las = 1 )> plot(NA, xlim = c(0,900), xaxs = "i", ylim = c(1/100,5), log = "y",+xlab = "Days since diagnosis",+ylab = "Mortality rate per 1 year" )> matshade(nm$tfe, cbind(ci.pred(p2, nm) * 365.25,+
ci.pred(p2, nf) * 365.25,+
ci.exp (p2, ctr.mat = list(nm, nf)) / 20 ),+lwd = 2, col = c('blue', 'red', 'black') )> abline(h = 1 / 20, lty = 3)
> axis(side = 4, at = c(2, 5, 10, 15, 20) / 200,+
labels = c(2, 5, 10, 15, 20) / 10)> axis(side = 4, at = c(2:9) / 200, labels = NA, tcl = -0.3)
From gure 2.5 we see the same overall tendency, but substantially more smoothed. But withthis type of analysis we have a more rm evidence that male mortality actually is higher inthe rst year or so.
This is one of the main advantages of a fully parametric approach to modeling of rates: It is
possible to show how the dierent baseline rates from a stratied model looks. The baseline
rates are not readily available from a Cox model.

2.3 Time-varying coecients
When it is suspected that eects of a given variable is not constant, one may allow thecoecient of a variable to vary by time:

λi (t) = λ0 (t) exp(β(t)xi + · · · )When we think of time as a covariate, this corresponds to an interaction between the covariate
and time, which is restricted by letting the x-eect be linear for any xed value of time.
The substantial reason for this particular choice of this particular form of interaction isslightly opaque. Given that one variable (time) is meticulously modelled it seems strange toinsist on a conditionally linear eect of x. It seems to be more productive to explore moreparsimonious parametrizations of interactions that were more directly addressing biologicallymeaningful deviations from the log-linear additivity of the eects.
There is however a tradition in epidemiological analysis of trends in rates to summarize
calendar time trends separately in each age-group by computing the average trend within

22

2.3

Time-varying coecients

WmtCma5.00

2.00

Mortality rate per 1 year

1.00
0.50

0.20
0.10

2
1.5

0.05

1
0.5

0.02
0.01

0.2
0

200

400

600

800

Estimated rates and rate-ratio by the gam tting machinery.
Days since diagnosis

Figure 2.5:./strat-pnsh

each age class. The continuous time version of this is precisely a varying coecients model
where the eect of calendar time is taken as linear at each age.
The simplest parametric form of such an interaction is to have separate eects of x, say, atdierent times. This would correspond to adding an interaction between x and some groupingof time. Again this approach can be taken

ad absurdumwith increasingly ne groupings oftime until we end up with the Cox-model formulation of the problem.
But when the main eect of time is modelled by a spline or any other smooth function,
implemented as columns of the model matrix in the Poisson regression model, we canestimate time-varying coecients by adding the same columns multiplied by x to the model
matrix. The coecients of these will then be the ones that determine the (time-varying) eectof the covariate x.
We use the same dataset as before, but now we have the interaction with the quantitative
variable age.

> summary(sL)Transitions:To
From
Alive Dead Records: Events: Risk time: Persons:
Alive 3432 165
3597
165
69593
228
> kn <- c(0, 50, 150, 450)
> m1 <- glm(cbind(lex.Xst == "Dead", lex.dur)+~ Ns(tfe, knots=kn) + age + sex,+family = poisreg,+
data = sL)2.3

Time-varying coecients

Examples23> mv <- glm(cbind(lex.Xst == "Dead", lex.dur)+~ Ns(tfe, knots=kn) + sex + Ns(tfe, knots=kn, i=T):age,+family = poisreg,+
data = sL)> anova(m1, mv, test = "Chisq")
Analysis of Deviance Table
Model 1: cbind(lex.Xst == "Dead", lex.dur) ~ Ns(tfe, knots = kn) + age +sexModel 2: cbind(lex.Xst == "Dead", lex.dur) ~ Ns(tfe, knots = kn) + sex +Ns(tfe, knots = kn, i = T):age
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1
3591
1298.623588
1280.3 318.245 0.0003914 ***--Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Here we see that there actually is a massive interaction  the age-eects
considerably by time.

does

vary

The parameters of interest are those from the second Ns term in the model, but of course
taken out as a curve. But it is also possible to extract the age-eect as a dierence between
two predictions, namely the rate-ratio between two persons, say 5 years apart in age,computed at range of times (tfe):> nx <- data.frame(tfe = seq(0, 1000, 10), age = 70, sex = "M")> nr <- data.frame(tfe = seq(0, 1000, 10), age = 65, sex = "M")> aRR <- ci.exp(mv, list(nx, nr)) # computes the ratio of predictions between nx and nr> matshade(nx$tfe / 365.25 * 12, aRR, plot = TRUE,+
lwd = 3,+log = "y", xlab = "Time since diagnosis (months)",+ylab = "RR per 5 years of age at diagnosis")> abline( h=1 )From the gure 2.6 we see that the age at diagnosis matters a lot for the mortality the rstfew months after diagnosis, but after about 3 months there is no eect.

Parametrizations
Note that when we use prediction data frames to tease out the eects, the particularparametrization does not matter, so we could have used a simple expression for the r.h.s. ofthe model formula:> ~ Ns( tfe, knots=kn ) * age + sex
and we would have obtained the same results.

Simplifying code
Since we are using a Lexis object as data base for the analysis, we already have specied thetime-structure in data, so we can shorten the code even further using the glm.Lexis function
for tting data. This function is really is designed to simplify analysis of rates in multistatemodels; basically it is just a wrapper using the poisreg family. It will by default analyze alltransition to any absorbing state, which in this case is "Dead", so the analysis could be done
by:

24

2.3

Time-varying coecients

WmtCma4.0
3.5

RR per 5 years of age at diagnosis

3.0
2.5

2.0

1.5

1.0

0

5

10
15
20
25
Time since diagnosis (months)

30

RR of death for the lung cancer patients, per 5 years of age at diagnosis. Resultsfrom a "varying-coecients" model  interaction between two continuous variables, where theeect of age is constrained to be linear at any time since diagnosis.
./time-var-Aeff
Figure 2.6:

> mL <- glm.Lexis(sL, formula = ~ Ns(tfe, knots = kn) * age + sex)
stats::glm Poisson analysis of Lexis object sL with log link:
Rates for the transition:
Alive->Dead> c(deviance(mv), deviance(mL))[1] 1280.321 1280.321
You want to look up the help page:> ?glm.Lexis

Chapter 3
So who do need the Cox-model?
Since everything which is possible using the Cox-model can be done using the Poissonmodeling of split data, there is no loss only substantial gain of capability by switching toPoisson modeling.
The Cox-model is computationally vastly more ecient, and it is easier to produce asurvival curve by standard software, which is relevant in most clinical studies. A drawback is
the overly detailed modeling of survival curves that may lead to over-interpretation of little
humps and notches on an estimated curve.
When stratication or time-dependent variables are involved, the facilities in the standard
Cox-analysis programs limits the ways in which the desired interactions can be modeled, andmoreover distracts the user from realizing that other interactions between covariates may beof interest.
Thus it seems that the Cox model is useful in the following cases: Clinical follow-up studies with only one relevant timescale and the focus on the eect ofother covariates than time.

 Studies analyzed on computing equipment pre-1985.
In other settings it seems preferable to split time and use the Poisson approach, because: it claries the distinction between events and (risk) time as response variables andtime(scales) as covariates  reected in the poisreg family in the Epi package.

 it enables smoothing of the eect of timescales using standard regression tools. In
particular it allows more credible estimates of survival functions in the simple case with
only time since entry as timescale.

 it enables sensible modeling of interactions between timescales and other variables (andbetween timescales), using standard regression tools.
Moreover, as the necessary computing power and software is available, the computational
problems encountered previously are now non-existent. Extraction of the relevant functions ofmodel parameters has been facilitated by the introduction of the possibility of supplying pairsof prediction data frames to extract rate-ratios (see the help page for ci.lin in the Epi
package).
That said, the user-interface to the Poisson modeling is slightly more complex than thatoered by standard packages for the Cox-model. This is because the Poisson approachrequires an explicit specication of a model for the eect of the timescale.

25

Chapter 4
Simple code for parametric model
The code for this chapter is here: http://bendixcarstensen.com/simple.R

4.1 Cookbook for a parametric survival curve
We use the lung dataset for illustration; categorical variables should always be declared asfactors:> lung$sex <- factor(lung$sex, labels = c("M", "W"))First declare the follow-up as a Lexis data frame:> Lx <- Lexis(exit = list(tfe = time),+
exit.status = factor(status,+
labels = c("Alive", "Dead")),+data = lung)
NOTE: entry.status has been set to "Alive" for all.
NOTE: entry is assumed to be 0 on the tfe timescale.
> summary(Lx)Transitions:To
From
Alive Dead
Alive
63 165

Records: Events: Risk time: Persons:
228
165
69593
228

current

Note that the () time variable is called tfethis will be the time variable to use in aparametric model.

Then split data in small intervals (so small that an assumption of constant rates in each isreasonable)in this case we use 20 days:> sL <- splitLexis(Lx, seq(0, 1100, 20))
Then specify a model for the eect of time using the gam function from the mgcv packagethrough the gam.Lexis function:> pmod <- gam.Lexis(sL, ~ s(tfe))
mgcv::gam Poisson analysis of Lexis object sL with log link:
Rates for the transition:
Alive->Dead

26

4.1

Cookbook for a parametric survival curve

Simple code for parametric model

27

This estimates how the mortality rate depends on tfe.
The Kaplan-Meier estimator provides estimates of the survival function at the observed
event times; a parametric model at any time. Thus, when predicting from a parametric model
we must specify the times at which we want the survival function calculated. This is in theform of a data frame of equidistant times (intervals of 10 days, say, this is unrelated to theinterval lengths used to split the follow-up for modeling):> nd <- data.frame(tfe = seq(0, 1000, 10))The survival function can then be computed (with condence intervals) at each of the times
supplied in nd (new data):> Sf <- ci.surv(pmod, nd)
NOTE: interval length chosen from as tfe[2] - tfe[1]> head(Sf)Estimate
2.5%
97.5%[1,] 1.0000000 1.0000000 1.0000000[2,] 0.9852730 0.9893000 0.9797459[3,] 0.9704037 0.9782776 0.9597349[4,] 0.9553965 0.9669351 0.9399581[5,] 0.9402563 0.9552769 0.9204052[6,] 0.9249886 0.9433092 0.9010654
We can then plot the estimated survival curve:> matshade(nd$tfe, Sf, plot = TRUE,+
lwd = 2, ylim = c(0, 1), yaxs = "i",+xlab = "Time since lung cancer (days)",+ylab = "Survival probability")

4.1.1

In summary

what you need to do in order to get a parametric survival curve is:> # Lexis object:> Lx <- Lexis(exit = list(tfe = time),+
exit.status = factor(status,+
labels = c("Alive", "Dead")),+data = lung)
NOTE: entry.status has been set to "Alive" for all.
NOTE: entry is assumed to be 0 on the tfe timescale.
> # Split follow up in small intervals> sL <- splitLexis(Lx, seq(0, 1100, 20))
> # Model rates in these intervals> pM <- gam.Lexis(sL, ~ s(tfe))
mgcv::gam Poisson analysis of Lexis object sL with log link:
Rates for the transition:
Alive->Dead> # Survival function times> nd <- data.frame(tfe = seq(0, 1000, 10))
> # Survival function estimate at these times> Sf <- ci.surv(pM, nd)
NOTE: interval length chosen from as tfe[2] - tfe[1]

28

4.1

Cookbook for a parametric survival curve

WmtCma> # Survival function plot
> matshade(nd$tfe, Sf, plot = TRUE,+
lwd = 2, ylim = c(0, 1), yaxs = "i",+xlab = "Time since lung cancer (days)",+ylab = "Survival probability")

1.0

Survival probability

0.8

0.6

0.4

0.2

0.0
0

200

400
600
Time since lung cancer (days)

800

1000

Parametric survival curve based on a gam model for the rates, compared with theK-M estimator (step function).
./simple-psurv
Figure 4.1:

4.1.2

Comparison with K-M:

For illustration we can plot the Kaplan-Meier curve for comparison with the parametriccurve:> matshade(nd$tfe, Sf, plot = TRUE,+
lwd = 2, ylim = c(0, 1), yaxs = "i",+xlab = "Time since lung cancer (days)",+ylab = "Survival probability")> km <- survfit(Surv(time, status) ~ 1, data = lung)> lines(km, lty = 1)

4.1

Cookbook for a parametric survival curve

4.1.3

Simple code for parametric model

29

The rates:

The parametric model provides estimates of rates (in this case in units of events / day), so we
multiply by 365.25 to get rates per 1 year:> Rf <- ci.pred(pM, nd) * 365.25> head(Rf)
Estimate
2.5%
97.5%
1 0.5419047 0.3929219 0.7473767
2 0.5554210 0.4091894 0.7539112
3 0.5692682 0.4257944 0.7610862
4 0.5834469 0.4426640 0.7690038
5 0.5979540 0.4597171 0.7777586
6 0.6127833 0.4768741 0.7874267. . . and we can plot the mortality rates as a function of time:

Mortality rate (per 1 PY)> matshade(nd$tfe, Rf, lwd = 2, log = "y", plot = TRUE,+
xlab = "Time since lung cancer (days)",+ylab = "Mortality rate (per 1 PY)")

2.0

1.0

0.5

0

Figure 4.2:

200

400
600
Time since lung cancer (days)

800

Mortality rates from a parametric model based on a gam model.

1000./simple-rates

30

4.2

Parametric proportional hazards model

scaledNote that when we show the rates it is necessary to consider themeasuredthey are

WmtCmaunits

in which they arequantities as opposed to survival probabilities. Note here that wehave rate estimates at more than 0.5 per 1 PY, something like 50% per yearcompatible
with the 1-year survival of approximately 50%.

4.2 Parametric proportional hazards model
The parametric model can easily be expanded with other (xed) covariates:> fM <- gam.Lexis(sL, ~ s(tfe) + age + sex)
mgcv::gam Poisson analysis of Lexis object sL with log link:
Rates for the transition:
Alive->Dead> round(ci.exp(fM), 3)
exp(Est.) 2.5% 97.5%(Intercept)
0.001 0.000 0.003
age
1.016 0.998 1.035
sexW
0.603 0.435 0.837
s(tfe).1
0.747 0.370 1.507
s(tfe).2
1.360 0.537 3.446s(tfe).3
1.208 0.774 1.884s(tfe).4
0.862 0.523 1.422s(tfe).50.850 0.575 1.255s(tfe).6
0.853 0.564 1.291
s(tfe).7
1.171 0.809 1.695s(tfe).8
1.787 0.417 7.663
s(tfe).9
1.293 0.849 1.969
This is a proportional hazards model; the hazards (rates, occurrence rates, mortality rates) as
function of tfe are proportional between, say, women and menthe W/M hazard ratio is thesame at all times, namely 0.60.
The model fM corresponds to a Cox-model; except that a Cox model does not provide
estimates of rates, but only of the hazard ratios:> cM <- coxph(Surv(time, status) ~ age + sex, data = lung)> round(ci.exp(cM), 3)
exp(Est.) 2.5% 97.5%
age
1.017 0.999 1.036
sexW
0.599 0.431 0.831
we see that the estimates of covariate eects are identical for all practical purposes.
There is a machinery to produce estimated survival curves from a Cox-model using theso-called Breslow-estimator, say for 60 year old men and women:> Sc <- survfit(cM, data.frame(age = 60, sex = c("M", "W")))> plot (Sc, col = c("blue", "red"), conf.int = TRUE)> lines(Sc, col = c("blue", "red"), lwd = 2)
From gure 4.3 you observe that the location of the jumps is the same for the two curves,
namely at each of the death times in the data set lung.
We can estimate the same curves from the parametric proportional hazards model, andoverlay these with the step-curves from the Cox model:

4.2

Parametric proportional hazards model

Simple code for parametric model

31

1.0

0.8

0.6

0.4

0.2

0.0
0

200

400

600

800

1000

Estimated survival curves for 60 year old men (blue) and women (red). Breslowestimator from a Cox model.
./simple-plCox
Figure 4.3:

> Mp <- ci.surv(fM, transform(nd, age = 60, sex = "M"))NOTE: interval length chosen from as tfe[2] - tfe[1]> Wp <- ci.surv(fM, transform(nd, age = 60, sex = "W"))NOTE: interval length chosen fromas tfe[2] - tfe[1]> matshade(nd$tfe, cbind(Mp, Wp), col = c("blue","red"),+xlab = "Time since lung cancer (days)",+ylab = "Survival probability",+
lwd = 2, ylim = c(0, 1), yaxs = "i", plot = TRUE)> lines(Sc, conf.int = TRUE)

4.2.1

Proportional hazards

can be tested using the cox.zph function that tests the absence of interactions between
covariates and time (non-proportionality):> cox.zph(cM)

32

4.2

Parametric proportional hazards model

WmtCma1.0

Survival probability

0.8

0.6

0.4

0.2

0.0
0

200

400
600
Time since lung cancer (days)

800

1000

Estimated survival curves for 60 year old men (blue) and women (red) from a modelwith smooth main eects of time since lung cancer and age. Breslow-estimators from a Coxmodel overlaid in black.
./simple-plPH
Figure 4.4:

chisq dfpage
0.209 1 0.65
sex
2.608 1 0.11
GLOBAL 2.771 2 0.25
The parametric counterpart to this is to compare the smooth parametric models with andwithout the interaction; for sex we see that the p-value is pretty much the same:> mM <- gam.Lexis(sL, ~ s(tfe, bs = "cr") + sex + age)mgcv::gam Poisson analysis of Lexis object sL with log link:
Rates for the transition:
Alive->Dead> iM <- gam.Lexis(sL, ~ s(tfe, bs = "cr", by = sex) + sex + age)mgcv::gam Poisson analysis of Lexis object sL with log link:
Rates for the transition:
Alive->Dead> anova(mM, iM, test = "Chisq")4.2

Parametric proportional hazards model

Simple code for parametric model

33

Analysis of Deviance Table
Model 1: cbind(trt(Lx$lex.Cst, Lx$lex.Xst) %in% trnam, Lx$lex.dur) ~ s(tfe,
bs = "cr") + sex + age
Model 2: cbind(trt(Lx$lex.Cst, Lx$lex.Xst) %in% trnam, Lx$lex.dur) ~ s(tfe,
bs = "cr", by = sex) + sex + age
Resid. Df Resid. Dev
Df Deviance Pr(>Chi)
13591.51299.523590.91298.0 0.54129 1.4919 0.1102
But testing for an interaction without showing the shape of it is a bit odd; so we rst plot theestimated rates for 60 year old men and women together, and beside that the estimated M/Wrate-ratio as functions of time since lung cancer:

5.0
Mortality rate ration between M and W

Mortality rate at 60 years at diagnosis

2.0

2.0

1.0

0.5

1.0

0.5

0.2

0.2

0.1
0

Figure 4.5:

200
400
600
800
Time since lung cancer (days)

1000

0

200
400
600
800
Time since lung cancer (days)

1000

Left: Mortality rates for men (blue) and women (red); right: M/W rate ratio.

./simple-2rates> par(mfrow = c(1,2))> # rates from M and W> nd <- data.frame(tfe = seq(0, 1000, 10), age = 60)> Mp <- ci.pred(iM, mutate(nd, sex = "M"))> Wp <- ci.pred(iM, mutate(nd, sex = "W"))> matshade(nd$tfe, cbind(Mp, Wp) * 365, plot = TRUE,+
col = c("blue","red"),+xlab = "Time since lung cancer (days)",+ylab = "Mortality rate at 60 years at diagnosis",+
lwd = 2, log = "y")> # rate-ratio between M and W> RR <- ci.exp(iM, list(mutate(nd, age = 60, sex = "M"),+mutate(nd, age = 60, sex = "W")))

34

4.2

Parametric proportional hazards model

WmtCma> matshade(nd$tfe, RR, plot = TRUE,+
xlab = "Time since lung cancer (days)",+ylab = "Mortality rate ration between M and W",+
lwd = 2, log = "y")> abline(h = 1)
It is pretty clear that there is some kind of interaction even if the test produces a p-value ofsome 10% (both in the parametric and the Cox model). The mortality rates in both sexes
increase by time, but more so among women, so the M/W rate ratio is clearly decreasing.
Neither of these two features of the mortality rates are available from the Cox-model, even if
they are essential for judging the interaction (non-proportionality).

References

Journal of Statistical Software[1] Bendix Carstensen and Martyn Plummer. Using Lexis objects for multi-state models in R.
, 38(6):118, 1 2011.

[2] DR Cox. Regression and life-tables (with discussion).
1972.

Statistics

J. Roy. Statist. Soc B, 34:187220,[3] Whitehead J. Fitting Cox's regression model to survival data using GLIM.
, 29(3):268275, 1980.

AppliedJournal of Statistical Software[4] Martyn Plummer and Bendix Carstensen. Lexis: An R class for epidemiological studies
with long-term follow-up.

, 38(5):112, 1 2011.

Epidemiol Perspect Innov[5] K. Rostgaard. Methods for stratication of person-time and events - a prerequisite for
Poisson regression and SIR estimation.

, 5:7, Nov 2008.

35