Introduction

A new joint model

stjm

Application

Future work

References

Flexible joint modelling of longitudinal and

survival data: The stjm command
17th Stata UK Users Group Meeting
Michael J. Crowther1 ,
Keith R. Abrams1 and Paul C. Lambert1,2
1 Centre

for Biostatistics and Genetic Epidemiology

Department of Health Sciences
University of Leicester, UK.

2 Department

of Medical Epidemiology and Biostatistics

Karolinska Institutet
Stockholm, Sweden.
mjc76@le.ac.uk

Michael J. Crowther

Stata UK Meeting

1 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Outline

Outline

I
I
I
I
I

Introduction to joint modelling

A new joint model
stjm
Example application - Primary Biliary Cirrhosis (PBC)
Future work

Michael J. Crowther

Stata UK Meeting

2 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Background

Background

Longitudinal response data affected by informative

dropout
Inclusion of time-varying covariates in survival analyses

Michael J. Crowther

Stata UK Meeting

3 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Background

Background

Longitudinal response data affected by informative

dropout
Inclusion of time-varying covariates in survival analyses

Michael J. Crowther

Stata UK Meeting

3 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Background

Background

Longitudinal response data affected by informative

dropout
Inclusion of time-varying covariates in survival analyses

Approaches:
I Latent class approach (Proust-Lima and Taylor, 2009)
I Shared parameter models - dependence through shared
random effects (Wulfsohn and Tsiatis, 1997)

Michael J. Crowther

Stata UK Meeting

3 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Background

Background

Longitudinal response data affected by informative

dropout
Inclusion of time-varying covariates in survival analyses

Approaches:
I Latent class approach (Proust-Lima and Taylor, 2009)
I Shared parameter models - dependence through shared
random effects (Wulfsohn and Tsiatis, 1997)

Michael J. Crowther

Stata UK Meeting

3 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Example dataset

Example dataset

Primary Biliary Cirrhosis (PBC) dataset - 312 patients

with 1945 repeated measurements of serum bilirubin
(Murtagh et al., 1994).
158 randomised to receive D-penicillamine and 154 to
placebo
Interested in the treatment effect after adjusting for the
longitudinal biomarker

Michael J. Crowther

Stata UK Meeting

4 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Data structure

Data structure
Longitudinal Data

Survival Data

yij = {yi(tij); j = 1,,mi}

Ti = min(Si, Ci), di = Si Ci

eg. Linear mixed

effects model

eg. Proportional
hazards model

Subject-specific
random effects
bi ~ MVN(0,)

Michael J. Crowther

Stata UK Meeting

5 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Data structure

Data structure
. list id logb drug _t0 _t _d if id==3 | id==20, noobs sepby(id)
id

logb

drug

_t0

_t

_d

3
3
3
3

.3364722
.0953102
.4054651
.5877866

D-penicil
D-penicil
D-penicil
D-penicil

0
.48187494
.99660498
2.0342789

.48187494
.99660498
2.0342789
2.7707808

0
0
0
1

20
20
20
20

1.629241
2.727853
2.406945
3.465736

placebo
placebo
placebo
placebo

0
.49556455
.91446722
3.6797721

.49556455
.91446722
3.6797721
3.7126274

0
0
0
1

Michael J. Crowther

Stata UK Meeting

6 / 27

Introduction

A new joint model

Application

stjm

Future work

References

Data structure

Data structure
KaplanMeier survival estimates

1.0

Survival probability

0.8

0.6

0.4

0.2

0.0
0

10

15

Followup time
95% Confidence Interval
95% Confidence Interval

Michael J. Crowther

Stata UK Meeting

Log(serum bilirubin) < 1

Log(serum bilirubin) > 1

7 / 27

Introduction

A new joint model

Application

stjm

Future work

References

Data structure

stjmgraph logb, panel(id)

15

Longitudinal response
0
2
2

Longitudinal response
0
2

Event

Censored

10
5
Time before censoring

15

10
5
Time before event

Figure: Longitudinal trajectories. Adjusted timescale.

Michael J. Crowther

Stata UK Meeting

8 / 27

Introduction

A new joint model

Application

stjm

Future work

References

Model specification

Longitudinal submodel
I

Linear mixed effects model (Laird and Ware, 1982):

yi (tij ) = Wi (tij ) + eij ,

eij N(0, e2 )

Wi (tij ) = xi0 (tij ) + zi0 (tij )bi + ui

.
I

xtmixed logb time drug || id:

time, cov(unstr)

Increase flexibility through the use of fixed/random

fractional polynomials of time (Royston and Altman,
1994).

Michael J. Crowther

Stata UK Meeting

9 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Model specification

Survival submodel
I

Many choices proposed including standard parametric

and, of course, Cox proportional hazards models
hi (Ti ) = h0 (Ti )exp(Wi (Ti ) + vi )
where vi is a set of baseline covariates and, for example;
Wi (Ti ) = 0i + 1i Ti + ui

is termed the association parameter; in this case we

assume the association is based on the current value of
the biomarker

Michael J. Crowther

Stata UK Meeting

10 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Model specification

We can now write down the full joint likelihood:

"Z
!
#
mi
n

Y
Y
f (yi (tij )|bi , ) f (bi |)f (Ti , di |bi , ) dbi
i=1
where

j=1



yi (tij ) Wi (tij )
f (yi (tij )|bi , ) = (2e2 )1/2 exp
,
2e2
(
)
b 0 V 1 bi
f (bi |) = (2|V |)1/2 exp i
,
2

and
Ti

 Z
f (Ti , di |bi , ) = [h0 (Ti )exp(Wi (t) + vi )]di exp


h0 (u)exp(Wi (u) + vi )du

Michael J. Crowther

Stata UK Meeting

11 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Model specification

We can now write down the full joint likelihood:

"Z
!
#
mi
n

Y
Y
f (yi (tij )|bi , ) f (bi |)f (Ti , di |bi , ) dbi
i=1
where

j=1



yi (tij ) Wi (tij )
f (yi (tij )|bi , ) = (2e2 )1/2 exp
,
2e2
(
)
b 0 V 1 bi
f (bi |) = (2|V |)1/2 exp i
,
2

and
Ti

 Z
f (Ti , di |bi , ) = [h0 (Ti )exp(Wi (t) + vi )]di exp


h0 (u)exp(Wi (u) + vi )du

Michael J. Crowther

Stata UK Meeting

11 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Numerical integration

Gauss-Hermite quadrature

Numerical method to approximate analytically intractable

integrals (Pinheiro and Bates, 1995)
Z

x 2

f (x)dx

m
X

wq f (xq )

q=1

Can be extended to multivariate integrals i.e. multiple

random effects

Michael J. Crowther

Stata UK Meeting

12 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Numerical integration

Slope
0

Node locations

0
Intercept

Figure: 2-dimensional 15 point node locations.

Michael J. Crowther

Stata UK Meeting

13 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Numerical integration

Slope
0

Node locations

0
Intercept

Figure: Adjusted 2-dimensional 15 point node locations.

Michael J. Crowther

Stata UK Meeting

14 / 27

Introduction

A new joint model

stjm

Application

Future work

References

A new joint model

Survival submodel

Flexible parametric survival model (Royston and Parmar,

2002; Lambert and Royston, 2009)
Modelled on the log cumulative hazard scale using
restricted cubic splines (Durrleman and Simon, 1989)
log{H0 (t)} = s{log(t)|, k0 }

Can evaluate the likelihood directly

Michael J. Crowther

Stata UK Meeting

15 / 27

Introduction

A new joint model

Application

stjm

Future work

References

Motivation

Marginal survival 1 degree of freedom

Survival probability

.8

.6

.4

.2
KaplanMeier curve
Marginal survival
0
0

2
3
Followup time

Figure: Predicted marginal survival function from joint model with 1 degree of
freedom, overlaid on the Kaplan-Meier survival curve.
Michael J. Crowther

Stata UK Meeting

16 / 27

Introduction

A new joint model

Application

stjm

Future work

References

Motivation

Marginal survival 5 degrees of freedom

Survival probability

.8

.6

.4

.2
KaplanMeier curve
Marginal survival
0
0

2
3
Followup time

Figure: Predicted marginal survival function from joint model with 5 degrees of
freedom, overlaid on the Kaplan-Meier survival curve.
Michael J. Crowther

Stata UK Meeting

17 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Motivation

Survival model linear predictor:

log{H(Ti |bi , vi )} = i = s{log(Ti )|, k0 } + Wi (Ti ) + vi

Michael J. Crowther

Stata UK Meeting

18 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Motivation

Survival model linear predictor:

log{H(Ti |bi , vi )} = i = s{log(Ti )|, k0 } + Wi (Ti ) + vi

Transform to the hazard and survival scales:



1 ds{log(Ti )|, k0 }
dW (Ti )
h(Ti |bi , vi ) =
+
exp(i )
Ti
d log(Ti )
dTi
S(Ti |bi , vi ) = exp{exp(i )}
(Crowther et al., 2011)

Michael J. Crowther

Stata UK Meeting

18 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Syntax

Implementation in Stata
stjm longdepvar [varlist ], panel(varname ) df(#)
[nodes(#) ...]

Michael J. Crowther

Stata UK Meeting

19 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Syntax

Implementation in Stata
stjm longdepvar [varlist ], panel(varname ) df(#)
[nodes(#) ...]
I Longitudinal submodel:
I
I
I

ffracpoly(numlist ) - Fixed FPs of time

rfracpoly(numlist ) - Random FPs of time
[varlist ] - Baseline covariates

Michael J. Crowther

Stata UK Meeting

19 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Syntax

Implementation in Stata
stjm longdepvar [varlist ], panel(varname ) df(#)
[nodes(#) ...]
I Longitudinal submodel:
I
I
I

ffracpoly(numlist ) - Fixed FPs of time

rfracpoly(numlist ) - Random FPs of time
[varlist ] - Baseline covariates

Survival submodel:
I
I

df(#)/knots(numlist ) - Baseline cum. hazard

survcov(varlist ) - Baseline covariates

Michael J. Crowther

Stata UK Meeting

19 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Syntax

Implementation in Stata
stjm longdepvar [varlist ], panel(varname ) df(#)
[nodes(#) ...]
I Longitudinal submodel:
I
I
I

Survival submodel:
I
I

ffracpoly(numlist ) - Fixed FPs of time

rfracpoly(numlist ) - Random FPs of time
[varlist ] - Baseline covariates
df(#)/knots(numlist ) - Baseline cum. hazard
survcov(varlist ) - Baseline covariates

Association:
I
I
I

nocurrent - Current value is the default

derivassoc(1) - 1st derivative
sepintassoc/sepassoc(numlist) - Random
coefficient, e.g. random intercept

Michael J. Crowther

Stata UK Meeting

19 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Predictions

Predictions
predict newvarname, option

Michael J. Crowther

Stata UK Meeting

20 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Predictions

Predictions
predict newvarname, option
I Longitudinal:
I
I
I
I

xb/fitted - Fitted values

residuals - Subject level residuals
rstandard - Standardised residuals
reffects/reses - Empirical Bayes predictions of
random effects

Michael J. Crowther

Stata UK Meeting

20 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Predictions

Predictions
predict newvarname, option
I Longitudinal:
I
I
I
I

xb/fitted - Fitted values

residuals - Subject level residuals
rstandard - Standardised residuals
reffects/reses - Empirical Bayes predictions of
random effects

Survival:
I
I
I
I
I

hazard - Hazard function

survival - Survival function
cumhazard - Cumulative hazard function
martingale - Martingale residuals
stjmcondsurv - Conditional survival

Michael J. Crowther

Stata UK Meeting

20 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Predictions

Predictions
predict newvarname, option
I Longitudinal:
I
I
I
I

xb/fitted - Fitted values

residuals - Subject level residuals
rstandard - Standardised residuals
reffects/reses - Empirical Bayes predictions of
random effects

Survival:
I
I
I
I
I

hazard - Hazard function

survival - Survival function
cumhazard - Cumulative hazard function
martingale - Martingale residuals
stjmcondsurv - Conditional survival

Predictions can be evaluated at measurement/survival times,

or user specified times.
Michael J. Crowther

Stata UK Meeting

20 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Predictions

Random intercept with fixed slope, current value

association, one degree of freedom:
.

stjm logb drug, ffracp(1) nodes(15) df(1) survcov(drug)

Random intercept and slope with fixed time powers 2 and

3, association based on 1st derivative, 3 degrees of
freedom:
.

stjm logb drug, rfracp(1) ffracp(2 3) survcov(drug) ///

nodes(15) df(3) nocurrent derivassoc(1)

Michael J. Crowther

Stata UK Meeting

21 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Application

Application to PBC dataset

. stjm logb trt, panel(id) nodes(15) rfracp(1) df(1) survcov(trt)
-> gen double timevar_1 = X^(1)
(where X = _t0)
Obtaining initial values:
Fitting full model:
Joint model estimates
Patient variable: id
Log-likelihood = -1952.7411
Coef.

Number of obs.
=
Number of patients =
Std. Err.

P>|z|

1945
312

[95% Conf. Interval]

Longitud.:
timevar_1
trt
_cons

.1636943
-.1705512
.682979

.0042372
.0349558
.0348529

38.63
-4.88
19.60

0.000
0.000
0.000

.1553895
-.2390633
.6146686

.171999
-.1020391
.7512894

Survival:
trt
_rcs1
_cons

-.0209648
.8489682
-3.330624

.177859
.0820932
.2450019

-0.12
10.34
-13.59

0.906
0.000
0.000

-.369562
.6880684
-3.810819

.3276325
1.009868
-2.850429

Association:
current

1.010613

.0836087

12.09

0.000

.8467429

1.174483

Michael J. Crowther

Stata UK Meeting

22 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Application

Random effects table

Random effects Parameters

Estimate

Std. Err.

[95% Conf. Interval]

id: Unstructured
sd(timevar_1)
sd(intercept)
corr(timevar_1,intercept)

.1588024
.8564879
.5405518

.0055511
.0173008
.0208844

.1482869
.8232413
.4983468

.1700637
.8910771
.580201

sd(Residual)

.3687027

.0067102

.3557827

.3820918

Longitudinal submodel: Linear mixed effects model

Survival submodel: Flexible parametric model with 1 degree of freedom
Integration method: Gauss-Hermite quadrature using 15 quadrature points

Michael J. Crowther

Stata UK Meeting

23 / 27

Introduction

A new joint model

Application

stjm

Future work

References

Predictions

Individual level predictions

stjmcondsurv, panel(id) id(16) fu(20)

Predicted conditional survival

Patient 16

3.5

1.0

1.0

3.0
0.8

1.0
0.4
0.5
0.0

0.2

0.5

2.0
0.6
1.5
1.0
0.4
0.5
0.0

Survival probability

0.6
1.5

0.8

2.5
Survival probability

2.0

log(serum bilirubin)

2.5
log(serum bilirubin)

Patient 70

3.5

3.0

0.2

0.5

1.0

0.0
0

10
15
Followup time

20

1.0

0.0
0

10
15
Followup time

20

Longitudinal response
Predicted conditional survival
95% pointwise confidence interval

Michael J. Crowther

Stata UK Meeting

24 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Future work

Future work
I

Survival submodels:
I
I
I
I

I
I
I
I
I

Weibull PH model
Gompertz PH model
2-component mixture Weibull PH model
Mixture Weibull-exponential PH model

Extension to competing risks

Extension to include a cure proportion
Longitudinal categorical responses
EM algorithm
Adaptive GH quadrature

Michael J. Crowther

Stata UK Meeting

25 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Future work

Command acknowlegdments

I
I
I
I
I

rcsgen - Paul Lambert

stpm2 - Paul Lambert
ghquadm - Bill Sribney
esttab - Ben Jann
fracgen - Patrick Royston

Michael J. Crowther

Stata UK Meeting

26 / 27

Introduction

A new joint model

stjm

Application

Future work

References

Future work

References I
M. J. Crowther, K. R. Abrams, and P. C Lambert. Flexible parametric joint modelling of longitudinal and survival
data. Submitted, 2011.
S. Durrleman and R. Simon. Flexible regression models with cubic splines. Statistics in Medicine, 8(5):551561,
1989.
N. M. Laird and J. H. Ware. Random-effects models for longitudinal data. Biometrics, 38(4):963974, 1982.
P. C Lambert and P. Royston. Further development of flexible parametric models for survival analysis. The Stata
Journal, 9:265290, 2009.
P. Murtagh, E. Dickson, M. Van Dam, G. Malincho, and P. Grambsch. Primary biliary cirrhosis: Prediction of
short-term survival based on repeated patient visits. Hepatology, 20:126134, 1994.
Jos
e C. Pinheiro and Douglas M. Bates. Approximations to the log-likelihood function in the nonlinear mixed-effects
model. Journal of Computational and Graphical Statistics, 4(1):pp. 1235, 1995. ISSN 10618600.
C. Proust-Lima and J. M. G. Taylor. Development and validation of a dynamic prognostic tool for prostate cancer
recurrence using repeated measures of posttreatment psa: a joint modeling approach. Biostatistics, 10(3):
535549, 2009.
P. Royston and D. G. Altman. Regression using fractional polynomials of continuous covariates: Parsimonious
parametric modelling. Journal of the Royal Statistical Society. Series C (Applied Statistics), 43(3):429467,
1994.
P. Royston and M. K. B. Parmar. Flexible parametric proportional-hazards and proportional-odds models for
censored survival data, with application to prognostic modelling and estimation of treatment effects. Stat Med,
21(15):21752197, 2002.
M. S. Wulfsohn and A. A. Tsiatis. A joint model for survival and longitudinal data measured with error.
Biometrics, 53(1):330339, 1997.

Michael J. Crowther

Stata UK Meeting

27 / 27