We present the frailtyHL package for fitting semi-parametric frailty models using h-likelihood. This package allows lognormal or gamma frailties for random-effect distribution, and it fits shared or multilevel frailty models for correlated survival data. Functions are provided to format and summarize the frailtyHL results. The estimates of fixed effects and frailty parameters and their standard errors are calculated. We illustrate the use of our package with three well-known data sets and compare our results with various alternative R-procedures.
Frailty models with a non-parametric baseline hazard have been widely adopted for the analysis of survival data (Hougaard 2000; Duchateau and P. Janssen 2008). The frailtyHL package (Ha, M. Noh and Y. Lee 2012) implements h-likelihood procedures (Ha, Y. Lee and J. K. Song 2001; Ha and Y. Lee 2003, 2005) for frailty models. The package fits Cox proportional hazards models (PHMs) with random effects, or frailties. The lognormal or gamma distribution can be adopted as the frailty distribution, corresponding to the normal or log-gamma distribution for the log frailties. H-likelihood obviates the need for marginalization over the frailty distribution, thereby providing a statistically efficient procedure for various random-effect models (Lee, J. A. Nelder and Y. Pawitan 2006).
The main function is frailtyHL()
. For instance,
> frailtyHL(Surv(time, status) ~ x + (1|id),
+ RandDist = "Normal",
+ mord = 0, dord = 1,
+ Maxiter = 200, convergence = 10^-6,
+ varfixed = FALSE, varinit = 0.1)
fits a lognormal frailty model. The first argument is a formula object,
with the response on the left of a ~
operator, and the terms for the
fixed and random effects on the right. The response is a survival object
as returned by the Surv
function (Therneau 2011). Here, time
and
status
denote survival time and censoring indicator taking value 1 or
0 for uncensored or censored observations, respectively; x
denotes a
fixed covariate and id
denotes the subject identifier. The expression
(1|id)
specifies a random intercept model ((x|id)
would specify a
random intercept, random slope model). The parameters mord
and dord
are the orders of Laplace approximations to fit the mean parameters
(mord
= 0 or 1) and the dispersion parameters (dord
= 1 or 2),
respectively. The Maxiter
parameter specifies the maximum number of
iterations and convergence
specifies the tolerance of the convergence
criterion. If varfixed
is specified as TRUE
(or FALSE
), the value
of one or more of the variance terms for the frailties is fixed (or
estimated) with starting value (e.g. 0.1) given by varinit
.
Previously, frailty models have been implemented in several R functions
such as the coxph()
function in the
survival package
(Therneau 2010) and the coxme()
function in the
coxme package
(Therneau 2011), based on penalized partial likelihood (PPL), the
phmm()
function in the
phmm package
(Donohue and R. Xu 2012), based on a Monte Carlo EM (MCEM) method, and the
frailtyPenal()
function in the
frailtypack package
(Gonzalez, V. Rondeau Y. Mazroui, A. Mauguen, and A. Diakité 2012), based on penalized marginal
likelihood. The phmm package fits one-component frailty models,
although it does allow multivariate frailties. The coxme()
function
can fit the multi-component model as shown in Example 2. Results from
frailtyHL are compared with those from survival, coxme and phmm.
Recently, the h-likelihood procedures of Lee and J. A. Nelder (1996) for fitting
hierarchical generalized linear models (HGLMs) have been implemented
using the hglm()
function in the
hglm package
(Alam, L. Ronnegard, and X. Shen 2010), the HGLMfit()
function in the
HGLMMM package
(Molas 2010) and the dhglmfit()
function in the
dhglm package
(Noh and Y. Lee 2011). The frailtyHL package is for fitting frailty models
with a non-parametric baseline hazard, providing estimates of fixed
effects, random effects, and variance components as well as their
standard errors. In addition, it provides a statistical test for the
variance components of frailties and also three AIC criteria for the
model selection.
The frailtyHL package makes it possible to
fit models with log-normal and gamma frailty distributions and
estimate variance components when the frailty structure is shared or nested.
For illustration, we present two models below and show how to fit these
models using frailtyHL()
in the Examples section.
Suppose that data consist of right censored time-to-event observations
from
Given the common unobserved frailty for the
where
The model represented by (1) is known as a shared or one-component frailty model (Hougaard 2000).
We can fit a multi-component model as below (Ha, Y. Lee and G. MacKenzie 2007):
For example, the CGD data (Fleming and D. P. Harrington 1991) have a multilevel structure in which patients nested within centers have recurrent infection times.
Later, we analyze these data using model (2) with
where
The h-likelihood
where
The functional form of
where
Lee and J. A. Nelder (1996, 2001) have proposed the use of the
Laplace approximation based on the h-likelihood when the marginal
likelihood,
The h-likelihood procedures use the Breslow method for handling tied
event times, while the PPL procedures allow the Efron method. For
estimating
Criterion | Literature | ||
Method | |||
HL(0,1) | Ha and Y. Lee (2003) | ||
HL(0,2) | Ha and Y. Lee (2003) | ||
HL(1,1) | Ha, M. Noh and Y. Lee (2012) | ||
HL(1,2) | Ha, M. Noh and Y. Lee (2012) | ||
coxph() |
Therneau (2010) | ||
for LN FM | |||
coxph() |
Therneau (2010) | ||
for gamma FM | |||
coxme() |
Therneau (2011) | ||
for LN FM | |||
phmm() |
Donohue and R. Xu (2012) | ||
for LN FM |
Table 1 shows historical developments of estimating criteria for frailty
models. The frailtyHL()
function provides estimators based on the
h-likelihood. As the orders in mord
and dord
increase, the bias of
estimator is reduced, but the calculation becomes computationally
intensive due to the extra terms. We recommend the use of HL(1,1) for
the lognormal frailty and of HL(1,2) for the gamma frailty. However,
HL(0,1) for the lognormal frailty and HL(0,2) for the gamma frailty
often perform well if
Based on the PPL methods, the coxph()
and coxme()
functions,
respectively, implement REMPL and MPL estimators for coxph()
function the maximum
likelihood (ML) estimators, maximizing the marginal likelihood coxph()
and coxme()
functions; Therneau (2010) recommended the
Efron method. For the lognormal frailty the ML estimator maximizing phmm()
function, but care must be taken to ensure
that the MCEM algorithm has converged (Donohue and R. Xu 2012). However, the
ML estimator can be biased when the number of nuisance parameters
increases with the sample size (Ha, M. Noh and Y. Lee 2010).
Furthermore, for the lognormal frailty the coxph()
function uses the
existing codes in linear mixed models so that it misses the
Suppose that HL(0,1) is used. The fitting algorithm is as follows:
Take
Given
Repeat Step 2 until the maximum absolute difference between the
previous and current estimates for
After convergence, we compute the estimates of the standard errors of
To demonstrate differences of various estimation methods in small
cluster size time
) is the time until infection since the insertion of the
catheter. The survival times for the same patient are likely to be
related because of a shared frailty describing the common patient’s
effect. The catheter is later removed if infection occurs and can be
removed for other reasons, which we regard as censoring; about 24% of
the data were censored.
We fit frailty models with two covariates, the sex (1 = male; 2 =
female) and age of each patient, using the functions (frailtyHL()
,
coxph()
, coxme()
and phmm()
) in the four packages. The results
are summarized in Table 2.
Sex | Age | Patient | |
Method | |||
lognormal model | |||
HL(0,1) | -1.380 | 0.005 | 0.535 |
(0.431) | (0.012) | (0.338) | |
HL(1,1) | -1.414 | 0.005 | 0.545 |
(0.432) | (0.012) | (0.340) | |
coxph() |
-1.388 | 0.005 | 0.551 |
(Breslow) | (0.441) | (0.012) | ( – ) |
coxph() |
-1.411 | 0.005 | 0.569 |
(Efron) | (0.445) | (0.013) | ( – ) |
coxme() |
-1.332 | 0.005 | 0.440 |
(Breslow) | (0.414) | (0.012) | ( – ) |
coxme() |
-1.355 | 0.004 | 0.456 |
(Efron) | (0.417) | (0.012) | ( – ) |
phmm() |
-1.329 | 0.004 | 0.378 |
(0.452) | (0.012) | ( – ) | |
gamma model | |||
HL(0,2) | -1.691 | 0.007 | 0.561 |
(0.483) | (0.013) | (0.280) | |
HL(1,2) | -1.730 | 0.007 | 0.570 |
(0.485) | (0.013) | (0.281) | |
coxph() |
-1.557 | 0.005 | 0.398 |
(Breslow) | (0.456) | (0.012) | ( – ) |
coxph() |
-1.587 | 0.005 | 0.412 |
(Efron) | (0.461) | (0.012) | ( – ) |
In PPL procedures (coxph()
and coxme()
), the Breslow method provides
slightly smaller estimate for fraityHL()
and coxph()
) give
larger estimates for phmm()
) and MPL (coxme()
)
procedures. However, both ML and MPL estimates from phmm()
and
coxme()
are somewhat different when cluster size is small,
coxph()
uses the
ML procedure, but it still gives smaller estimate for
The current version of the frailtyHL package allows multi-component (multilevel) frailties. Allowance for correlation (Ha, R. Sylvester, C. Legrand and G. MacKenzie 2011) between random effects, for example correlations between random center effect and random treatment effect, is currently in progress. Other developments include dispersion frailty models based on double HGLMs (Lee and J. A. Nelder 2006), which allow fixed and random effects in dispersion parts (i.e. variance of random effects) of the model as well as the hazard functions, and joint modelling (Ha, T. Park and Y. Lee 2003) of HGLMs and frailty models.
The data set presented by Mantel, N. R. Bohidar and J. L. Ciminera (1977) is based on a
tumorigenesis study of 50 (time
) is the time to
development of tumor, measured in weeks. Death before occurrence of
tumor yields a right-censored observation; forty rats developed a tumor,
leading to censoring of about 73%. The survival times for rats in a
given litter may be correlated due to a random effect representing
shared genetic or environmental effects.
We fit frailty models with one covariate, rx (1 = drug; 0 = placebo),
using frailtyHL()
. Below, we present the code and results for the
lognormal frailty model with HL(1,1). The output from the R code shows
that the effect of rx is significant (t-value = 2.823 with p-value =
0.005). That is, the rx group has a significantly higher risk than in
control group. Here, the variance estimate of the frailty is
Note that although we report the SE of varfixed = TRUE
and varinit = 0
in the frailtyHL()
: see below.
> library(survival)
> data(rats)
> frailtyHL(Surv(time,status) ~ rx + (1|litter),
+ data = rats,
+ varfixed = TRUE, varinit = 0)
iteration :
4
convergence :
4.801639e-09
[1] "converged"
[1] "Results from the Cox model"
[1] "Number of data : "
[1] 150
[1] "Number of event : "
[1] 40
[1] "Model for conditional hazard : "
Surv(time, status) ~ rx + (1 | litter)
[1] "Method : HL(0,1)"
[1] "Estimates from the mean model"
Estimate Std. Error t-value p-value
rx 0.8982 0.3174 2.83 0.004655
[1] "Estimates from the dispersion model"
Estimate Std. Error
litter "0" "NULL"
-2h0 -2*hp -2*p_b,v(hp)
[1,] 363.69 363.69 364.15
cAIC mAIC rAIC
[1,] 365.69 365.69 364.15
> frailtyHL(Surv(time,status) ~ rx + (1|litter),
+ data = rats, RandDist = "Normal",
+ mord = 1, dord = 1)
iteration :
87
convergence :
9.97616e-07
[1] "converged"
[1] "Results from the log-normal frailty model"
[1] "Number of data : "
[1] 150
[1] "Number of event : "
[1] 40
[1] "Model for conditional hazard : "
Surv(time, status) ~ rx + (1 | litter)
[1] "Method : HL(1,1)"
[1] "Estimates from the mean model"
Estimate Std. Error t-value p-value
rx 0.9107 0.3226 2.823 0.004754
[1] "Estimates from the dispersion model"
Estimate Std. Error
litter 0.4272 0.4232
-2h0 -2*hp -2*p_v(hp) -2*p_b,v(hp)
[1,] 335.97 397.36 362.14 362.56
cAIC mAIC rAIC
[1,] 362.22 366.14 364.56
The code and results for the gamma frailty model with HL(1,2) are
presented below. The output shows that these results are similar to
those of the lognormal frailty, particularly for estimation of
> frailtyHL(Surv(time,status) ~ rx + (1|litter),
+ data = rats, RandDist = "Gamma",
+ mord = 1, dord = 2)
iteration :
170
convergence :
9.567765e-07
[1] "converged"
[1] "Results from the gamma frailty model"
[1] "Number of data : "
[1] 150
[1] "Number of event : "
[1] 40
[1] "Model for conditional hazard : "
Surv(time, status) ~ rx + (1 | litter)
[1] "Method : HL(1,2)"
[1] "Estimates from the mean model"
Estimate Std. Error t-value p-value
rx 0.9126 0.3236 2.82 0.004806
[1] "Estimates from the dispersion model"
Estimate Std. Error
litter 0.5757 0.5977
-2h0 -2*hp -2*p_v(hp) -2*s_v(hp)
[1,] 331.60 413.85 365.35 361.71
-2*p_b,v(hp) -2*s_b,v(hp)
365.77 362.12
cAIC mAIC rAIC
[1,] 365.30 365.71 364.12
For the selection of a model among nested or non-nested models such as
lognormal and gamma frailty models, we may use three Akaike information
criteria (AIC)
(Lee, J. A. Nelder and Y. Pawitan 2006; Ha, Y. Lee and G. MacKenzie 2007; Donohue, R. Overholser, R. Xu, and F. Vaida 2011)
based on conditional likelihood, marginal likelihood and restricted
likelihood, respectively, defined by
In the data set, in the Cox model cAIC=365.69, mAIC=365.69 and
rAIC=364.15, and in the lognormal frailty model cAIC=362.22, mAIC=366.14
and rAIC=364.56, and in the gamma frailty model cAIC=365.30, mAIC=365.71
and rAIC=364.12. The likelihood tests based upon the REMPL showed the
absence of frailty effect (
The CGD data set presented by Fleming and D. P. Harrington (1991) consists of a
placebo-controlled randomized trial of gamma interferon (rIFN-g) in the
treatment of chronic granulomatous disease (CGD). 128 patients (id
)
from 13 centers (tstop
-tstart
) are the recurrent infection times
of each patient from the different centers. Censoring occurred at the
last observation for all patients, except one, who experienced a serious
infection on the day he left the study; in the CGD study about 63% of
the data were censored. The recurrent infection times for a given
patient are likely to be correlated. However, each patient belongs to
one of the 13 centers; hence, the correlation may also be attributed to
a center effect.
Ignoring important random components may render invalid many of the
traditional statistical analysis techniques. We fit a multilevel
lognormal frailty model with two frailties and a single covariate, treat
(rIFN-g, placebo), using frailtyHL()
. Here, the two frailties are
random center and patient effects. The code and results using HL(1,1)
are provided below. The output shows that the effect of treatment is
significant (t-value = -3.476 with p-value < 0.001), indicating that
rIFN-g significantly reduces the rate of serious infection in CGD
patients. The estimate of variance of patient frailty
(
> library(survival)
> data(cgd)
> frailtyHL(Surv(tstop-tstart,status) ~ treat +
+ (1|center) + (1|id),
+ data = cgd,
+ RandDist = "Normal", mord = 1, dord = 1)
iteration :
157
convergence :
9.336249e-07
[1] "converged"
[1] "Results from the log-normal frailty model"
[1] "Number of data : "
[1] 203
[1] "Number of event : "
[1] 76
[1] "Model for conditional hazard : "
Surv(tstop-tstart,status)~treat+(1|center)+(1|id)
[1] "Method : HL(1,1)"
[1] "Estimates from the mean model"
Estimate Std. Error t-value p-value
treatrIFN-g -1.184 0.3407 -3.476 0.0005085
[1] "Estimates from the dispersion model"
Estimate Std. Error
center 0.02986 0.1572
id 1.00235 0.5089
-2h0 -2*hp -2*p_v(hp) -2*p_b,v(hp)
[1,] 603.30 853.66 692.63 692.95
cAIC mAIC rAIC
[1,] 684.92 698.63 696.95
For testing the absence of a random component (
Cox model without frailty
model without patient effect
model without center effect
multilevel model above requiring both patient and center effects
The deviance difference between M3 and M4 is
Let us choose a final model using information criteria. For M1 we have cAIC=708.68, mAIC=708.68 and rAIC=707.48; for M2 cAIC=702.96, mAIC=706.88 and rAIC=705.66; for M3 cAIC=684.84, mAIC=696.68 and rAIC=694.99; for M4 cAIC=684.92, mAIC=698.63 and rAIC=696.95. All of the three criteria choose M3 in the CGD data set.
Using the examples in the previous section, we compare the outcomes from
frailtyHL and other packages. We consider the three functions (
coxph()
, coxme()
and phmm()
) for the lognormal frailty model and
the coxph()
function for the gamma frailty model.
Example 1: Rat data The codes of coxph()
, coxme()
and phmm()
for fitting lognormal frailty model are, respectively, as follows:
>coxph(Surv(time, status) ~ rx +
+ frailty(litter, dist = "gauss"),
+ method = "breslow", rats)
> coxme(Surv(time, status) ~ rx + (1|litter),
+ ties = "breslow", rats)
> phmm(Surv(time, status) ~ rx + (1|litter),
+ data = rats, Gbs = 2000, Gbsvar = 3000,
+ VARSTART = 1, NINIT = 10,
+ MAXSTEP = 200, CONVERG=90)
Table 3 summarizes the estimation results. Even though cluster size
coxme()
and phmm()
are somewhat
different in Table 2 when
Next, the code of coxph()
for fitting the gamma frailty model is
below:
> coxph(Surv(time, status) ~ rx +
+ frailty(litter, dist = "gamma"),
+ method = "breslow", rats)
The results of frailtyHL()
(HL(0,2), HL(1,2)) and coxph()
for gamma
frailty are also presented in Table 3. For the estimation of frailtyHL()
and coxph()
are similar, but for
frailtyHL()
(coxph()
(
Rx | Litter | |
Method | ||
lognormal model | ||
HL(0,1) | 0.906 (0.323) | 0.427 (0.423) |
HL(1,1) | 0.911 (0.323) | 0.427 (0.423) |
coxph() |
0.905 (0.322) | 0.395 ( – ) |
(Breslow) | ||
coxph() |
0.913 (0.323) | 0.412 ( – ) |
(Efron) | ||
coxme() |
0.905 (0.322) | 0.406 ( – ) |
(Breslow) | ||
coxme() |
0.913 (0.323) | 0.426 ( – ) |
(Efron) | ||
phmm() |
0.920 (0.326) | 0.449 ( – ) |
gamma model | ||
HL(0,2) | 0.908 (0.324) | 0.575 (0.598) |
HL(1,2) | 0.913 (0.324) | 0.576 (0.598) |
coxph() |
0.906 (0.323) | 0.474 ( – ) |
(Breslow) | ||
coxph() |
0.914 (0.323) | 0.499 ( – ) |
(Efron) |
Example 2: CGD data The code of the coxme()
function for fitting
multilevel lognormal frailty model is as follows:
> coxme(Surv(tstop - tstart,status) ~
+ treat + (1|center) + (1|id),
+ ties = "breslow", cgd)
The results of frailtyHL()
(HL(0,1), HL(1,1)) and coxme()
are
summarized in Table 4. The results from HL and PPL methods for frailty
parameters become more similar because the cluster sizes (the number of
patients from different centers) are somewhat large, ranging from 4 to
26.
Treat | Center | Patient | |
Method | |||
lognormal model | |||
HL(0,1) | -1.074 | 0.026 | 0.982 |
(0.335) | (0.153) | (0.501) | |
HL(1,1) | -1.184 | 0.030 | 1.002 |
(0.341) | (0.157) | (0.509) | |
coxme() |
-1.074 | 0.033 | 0.939 |
(Breslow) | (0.333) | ( – ) | ( – ) |
coxme() |
-1.074 | 0.032 | 0.947 |
(Efron) | (0.333) | ( – ) | ( – ) |
The h-likelihood method offers novel possibilities to fit various models
with random effects. The frailtyHL package for frailty models
eliminates the nuisance parameters
This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (No. 2009-0088978, No. 2010-0021165 and No. 2010-0011372).
frailtyHL, survival, coxme, phmm, frailtypack, hglm, HGLMMM, dhglm
Agriculture, ClinicalTrials, Econometrics, MixedModels, Spatial, Survival
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
Ha, et al., "frailtyHL: A Package for Fitting Frailty Models with H-likelihood", The R Journal, 2012
BibTeX citation
@article{RJ-2012-010, author = {Ha, Il Do and Noh, Maengseok and Lee, Youngjo}, title = {frailtyHL: A Package for Fitting Frailty Models with H-likelihood}, journal = {The R Journal}, year = {2012}, note = {https://rjournal.github.io/}, volume = {4}, issue = {2}, issn = {2073-4859}, pages = {28-37} }