A Multi-parameter regression model for interval censored survival data

by   Defen Peng, et al.
University of Limerick

We develop flexible multi-parameter regression survival models for interval censored survival data arising in longitudinal prospective studies and longitudinal randomised controlled clinical trials. A multi-parameter Weibull regression survival model, which is wholly parametric, and has non-proportional hazards, is the main focus of the paper. We describe the basic model, develop the interval-censored likelihood and extend the model to include gamma frailty and a dispersion model. We evaluate the models by means of a simulation study and a detailed re-analysis of data from the Signal Tandmobiel^ study. The results demonstrate that the multi-parameter regression model with frailty is computationally efficient and provides an excellent fit to the data.



There are no comments yet.


page 1

page 2

page 3

page 4


Semiparametric analysis of clustered interval-censored survival data using Soft Bayesian Additive Regression Trees (SBART)

Popular parametric and semiparametric hazards regression models for clus...

An Ensemble Method for Interval-Censored Time-to-Event Data

Interval-censored data analysis is important in biomedical statistics fo...

Multi-Parameter Regression Survival Modelling with Random Effects

We consider a parametric modelling approach for survival data where cova...

Multi-Parameter Regression Survival Modelling: An Alternative to Proportional Hazards

It is standard practice for covariates to enter a parametric model throu...

Parmsurv: a SAS Macro for Flexible Parametric Survival Analysis with Long-Term Predictions

Health economic evaluations often require predictions of survival rates ...

Proportional hazards model with partly interval censoring and its penalized likelihood estimation

This paper considers the problem of semi-parametric proportional hazards...

Joint Modeling of Longitudinal and Survival Data with Censored Single-index Varying Coefficient Models

In medical and biological research, longitudinal data and survival data ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Interval censored survival data can arise in longitudinal epidemiological studies where the response variable

is binary. Typically, at baseline, , patients start in an initial state, e.g., for the th patient (say) and later, as follow-up proceeds at scheduled inspection times, the event of interest may occur at a time whence, , where . This leads naturally to the use of “time to event” survival modelling in order to determine the effect of selected risk factors measured at baseline on the time to the event of interest. The interval censoring arises because patients are not monitored continuously but rather a finite schedule of follow-up examinations at times , . Thus, if, for the th subject, the event occurs between the th and the th follow-up examination the binary indicators become for and for and then the time to event is .

As a motivating study we consider the Tandmobiel study which is a longitudinal prospective oral health study conducted in Flanders (Belgium) from 1996 to 2001. A cohort of 4430 randomly sampled school children who attended the first year of the primary school at the beginning of the study were examined annually (6 times). The response was time to the emergence of the permanent upper left first premolars (tooth 24 in European dental notation). When emergence occurs between annual follow-ups the exact time of emergence is not known and the time to the emergence is said to be interval censored.

Modelling such data presents a variety of challenges. MacKenzie and Peng (2013)

developed interval censoring methods for parametric models, some of which were non-PH and they compared the use of standard right censored likelihoods based on mid-points with interval censored likelihoods. They showed that the use of mid-points led to artificially precise estimators in PH models when analysing time to loss of vision in a longitudinal trial of Age-Related Macular Degeneration (ARMD)

(Hart, 2002). Survival data arising in longitudinal vision studies, were analysed earlier, by MacKenzie (1999) and later by Al-tawarah and MacKenzie (2002). See also Finkelstein (1986) for an important early paper in the field, Huang and Wellner (1997) for a more theoretical review, and the books by Sun (2006) and Bogaerts et al. (2017) for comprehensive treatments of the subject. Nearly all of these papers employ models where covariates enter through a single parameter. In this paper such models are designated as single parameter regression (SPR) models.

In contrast, the concept of multi-parameter regression (MPR) survival modelling was developed in Burke and MacKenzie (2017). In MPR survival models the scale and shape parameters are modelled simultaneously by means of two separate linear predictors: these models are parametric and intrinsically more flexible than classical proportional hazards (PH) survival models. In the 2017 paper, MPR models were investigated in the context of right censored survival data from the Northern Ireland Lung Cancer Study (Wilkinson, 1995).

In this paper, we extend MPR models for interval censored survival data arising in longitudinal studies and introduce a MPR model with gamma frailty and a dispersion model to re-analyse data from the Signal Tandmobiel study (Bogaerts et al., 2002; Gómez et al., 2009).

2 MPR modelling framework

We envisage the class of two-parameter parametric survival models supporting scale and shape parameters. Within that class we model the scale and shape parameters simultaneously by means of two separate linear predictors which may involve the same set, or different sets, of covariates. In this paper we focus on the Weibull MPR model in order to illustrate the MPR approach to analysing interval-censored data. This model has proved useful in other contexts and has the advantage of directly extending a standard proportional hazards model.

2.1 Weibull MPR survival model

The Weibull MPR model of Burke and MacKenzie (2017) can be defined through its hazard function


which arises from a basic (no-covariate) Weibull hazard function, , (Lawless, 2003; Marshall and Olkin, 2007) allowing both its scale, , and shape, , parameters to depend on covariates via the MPR specification


where and

are the scale and shape covariate vectors, and

and are the associated regression coefficient vectors. Clearly the hazard function, (2.1), depends on covariates but we avoid for notational convenience; similarly we avoid and in (2.2). Note that when the model reduces to a standard PH regression survival model, i.e., an SPR model. Although the model is non-PH, when an coefficient, say , is zero, the corresponding covariate is PH, whence the hazard ratio is (Burke and MacKenzie, 2017). Thus, the Weibull MPR model usefully supports a mix of PH and non-PH effects, including crossing hazards and, through the coefficients, provides covariate-specific tests of proportionality of hazards.

2.2 Frailty extension

We extend the Weibull MPR model of Burke and MacKenzie (2017) to incorporate a multiplicative frailty term via the conditional hazard

where is defined in (2.1), and

is the frailty term, which we will assume follows a gamma distribution with density

where , such that

and the frailty variance is

. This is the classical frailty model in which the random effect, , measures additional person-specific heterogeneity not accounted for by the covariates. (Vaupel et al., 1979)

Since is an unobserved variable, the marginal distribution, obtained by integrating over , has survivor function given by


where the subscript “” indicates marginal, and is the cumulative hazard function associated with (2.1). Note that, in the absence of frailty, i.e., , we have that restoring the familiar (non-frailty) relationship between a survivor function and its cumulative hazard. See Hougaard (1995, 2000) and Duchateau and Janssen (2008) for more details on frailty models.

2.3 Dispersion model

The MPR frailty model can be extended further with advantage. Usually the frailty variance, , is a constant, but we can allow the frailty variance to be person-specific via another regression model, i.e.,


where and are vectors of covariates and their coefficients respectively. This dispersion model (DM) allows one to investigate the structure of the dispersion and provides a convenient framework for testing the homogeneity of frailty variances among covariates (e.g., between sexes). In addition, when the frailty variance is unstructured, the model reduces to the Weibull MPR frailty model of Section 2.2.

The concept of modelling the structure of the dispersion can be traced back to joint mean-dispersion modelling (Smyth, 1989; Lee and Nelder, 2001; Pan and MacKenzie, 2003), but its adoption in the frailty paradigm, in survival analysis, is more recent (Lynch and MacKenzie, 2014). Furthermore, the combination of frailty dispersion modelling in combination with an underlying MPR model is novel in the literature. It will be apparent that the frailty dispersion model which introduces a third regression (2.4) is an entirely natural development in the MPR paradigm.

3 Estimation

3.1 Likelihood Functions

In most longitudinal studies the idea of a fixed schedule of follow-ow examinations is too rigid as many subjects fail to respect their exact re-examination appointment dates. Accordingly it is usual to allow the intervals to be person-specific such that . In general, is close to the scheduled and of course . This notation, whilst accurate, is rather cumbersome and it is convenient to abbreviate it to in the equations which follow.

Then a general likelihood for inteval-censoring (IC) data is given by

with subject specific intervals , where time . Here denotes an interval censored observation and denotes a right censored observation with censoring time . In this IC setting, the th subject either “fails” in interval , or is right-censored. In total, there are patients of whom are right censored or withdrawn at specific times, leaving patients who are interval-censored. Note that the interval-censored subjects play the same role as “failures” in the right-censored setting and often right-censoring occurs at times completely unrelated to the scheduled follow-up examinations, e.g., an early withdrawal from the study. Thus, this notational setup is advantageous when it is important to distinguish between interval-censored and right-censored observations. See MacKenzie and Peng (2013) for more details on this approach.

Here, however, we re-write the likelihood above as


for notational convenience, in which we define the th right censored observation as lying in the interval . Accordingly, there are now intervals, with . This representation is the most commonly occurring form in the IC literature and we use it below.

The IC likelihood for the Weibull MPR DM model (i.e., the most general model of Section 2) is


where , and are the cumulative hazard functions for the th individual evaluated at the end-points of (where, for notational convenience, we avoid expressions such as and ), and from (2.2), from (2.4), and , , and are the scale, shape, and dispersion covariate vectors respectively (which may or may not contain covariates in common).

3.2 Score Functions

We now let so that the log-likelihood function can be written as . The score functions are then given by

all of which have a similar functional form, differing only with respect to the “weight” functions, which are given by

with , , and analogously defined by replacing with .

Although, in the above, we intend that for the purpose of the current paper (i.e., that of a Weibull MPR model), we have written the above score functions in a generic form so that can be replaced by any cumulative hazard function. If the underlying MPR model had another (positive) shape parameter, say, , modelled as , then we would gain an additional score function where , i.e., this score function has the same structure as that of and , but with a different . On the other hand, if the frailty distribution was changed, the factor in all of the score functions would change (and not only through changing), and, of course, the form of would also change. Thus, although we focus on a Weibull-gamma frailty model, the above is easily adapted to a wide range of MPRDM models for IC data.

4 Model selection

For the purpose of selecting among models within the MPRDM class, standard information criteria may be used, namely, the Akaike Information Criterion, AIC , and the Bayesian Information Criterion, BIC where is the maximum likelihood estimator and . There are two levels of model selection, both of which can be handled by these information criteria, namely: (a) the overall model type, and (b) the covariate set for each regression component within a given model type.

4.1 Model types

The MPRDM modelling framework introduced in Section 2 2 is quite general, containing a range of new and existing regression model types. A natural hierarchy of model types emerges as follows: the underlying model may be PH ( regression) or MPR ( and regressions), and the frailty component may be absent, present, or present with a regression component. The six model types are summarised in Table 1. Note that models PH and PHF are single parameter regression (SPR) models (only one parameter depends on covariates) while all other types are multi-parameter regression (MPR) models – in particular, the PHDM model has a PH baseline component (i.e., SPR), but the overall marginal model is MPR since depends on covariates. We have found models without a regression to be less useful and, so, these are not considered here.

Model Baseline Frailty Yes No
“Model” is the name of the model; “Baseline” is the baseline covariate structure such that “PH” is a Proportional Hazards structure where only the scale parameter, , depends on covariates whereas “MPR” is a Multi-Parameter Regression structure where the shape parameter, , also depends on covariates (see “Regression” columns); “Frailty” indicates the presence of a frailty term; “Regression” highlights the regression components via the distributional parameters which depend on covariates (“Yes”) and which do not depend on covariates (“No”).
Table 1: Model types

4.2 Covariates

Given a particular model type from Table 1, we will generally wish to select from a set of candidate covariates, say , to appear in the model (note: is used for the intercept term). In the most general MPRDM model, this amounts to the selection of scale covariates, , shape covariates, , and frailty dispersion covariates, , where the subsets may or may not overlap. While the union, , is of interest as these covariates affect survival in some way, so too are the , , and vectors themselves as these characterise the nature of specific covariate effects, e.g, in an MPR model without frailty, implies that is a non-PH covariate, and, in a frailty dispersion model, indicates that the frailty variance differs in the sub-groups defined by .

In general, the basic parameters of survival models (including the frailty variance) are rarely orthogonal, i.e., estimates of these parameters will be correlated. When covariates enter these parameters, this correlation propagates to the regression coefficients. In particular, if the covariate appears in all regression components simultaneously, then the estimates of its corresponding , , and coefficients tend to be quite correlated. It is important to emphasize that this correlation does not lead to convergence issues in model fitting, nor does it imply that a covariate must only appear in one regression component within the model. However, it does have implications for variable selection, e.g., individual Wald-based significance tests (which account only for the variance of estimates, and not covariance – which is important in this context) might render a particular covariate non-significant in all regression coefficients, when, in fact, the overall effect is significant. Covariate selection in MPR models was developed in Burke and MacKenzie (2017) who suggested the use of stepwise procedures in which covariate additions/deletions are carried out for each regression component separately as well as simultaneously, e.g., in an MPR model the covariate could be added to first (but not to ), then to (but not to ), and finally to and simultaneously.

5 Simulation study

We conducted a simulation study to assess the estimation properties of the Weibull MPR model for interval-censored data. Failure times were generated from the Weibull regression model (with or without frailty) with two covariates: , a binary covariate where mimicking the treatment effect for example, and a continuous baseline covariate distributed as .

In addition, we constructed trajectories for each individual in the study by constructing intervals such that and , where and

are independent continuous variables with uniform distribution in the interval

. Zhang (2009) used this approach which, by construction, defines intervals which are non-informative about the survival time distribution, . Furthermore, it can be shown that (proof omitted). In this simulation study, we set , i.e., , so that the average inspection length is proportional to the average survival time, , where we use .

In the simulations, the proportion of right (random) censoring was controlled by using an exponential distribution where the estimate of the controlling parameter,

, was obtained from the “-function” approach of MacKenzie and Peng (2013). Suppose the independent censoring times follow an exponential distribution with density . Let

where is the survival function and is the censoring proportion required. Then, ensures that, on average, proportion of censored individuals in each simulation equals . We set in this simulation study. The entire simulation was conducted in the R software package. (R Core Team, 2018)

Three sample sizes were used . Each scenario was replicated 5000 times, and, for each replicate, the model was fitted using the likelihood function given in (3.6). Note that the true coefficient values were set as .

5.1 Results

Table 2

shows the median estimate, standard error, and average relative percentage bias (

) for each parameter arising from our simulation study. We see that the estimators from the MPR IC likelihood have very little bias, and that both the bias and standard error tend to reduce with increasing sample size. As expected, larger right-censoring and inspection tend to reduce performance, but, even in the worst case of 30% right-censoring and average inspection length of 0.5 , the results are quite good.

Table 3 and Table 4 show the results of the simulation from the Weibull MPR with frailty, and Weibull MPR with dispersion model, , respectively. For the smallest sample size () the frailty estimates can be somewhat biased (and even for for the more complicated dispersion model case), but the bias is much reduced for . Of course, we expect that more complicated models tend to require larger sample sizes, but, overall, our simulation suggests that the model parameters are recoverable for reasonable sample sizes.

Intercept 2.04 0.11 2.55 2.04 0.08 2.94
0.51 0.16 1.15 0.25 0.11 -0.71
0.30 0.16 0.82 -0.11 0.11 5.08
Intercept 2.01 0.06 0.62 2.01 0.05 1.05
0.51 0.10 1.51 0.25 0.07 1.05
0.30 0.10 0.27 -0.10 0.07 3.22
Intercept 2.00 0.05 0.32 2.01 0.04 0.73
0.50 0.07 0.84 0.25 0.05 0.10
0.30 0.07 0.47 -0.10 0.05 0.14
Intercept 2.04 0.13 2.84 2.04 0.09 3.14
0.51 0.20 1.32 0.25 0.13 0.31
0.30 0.20 1.00 -0.10 0.13 1.86
Intercept 2.01 0.08 0.97 2.02 0.06 1.33
0.51 0.12 1.21 0.25 0.08 0.66
0.30 0.13 0.86 -0.10 0.08 0.96
Intercept 2.01 0.06 0.51 2.01 0.04 0.78
0.50 0.09 0.44 0.25 0.06 0.05
0.30 0.09 0.03 -0.10 0.06 0.79
Intercept 2.03 0.11 1.99 2.03 0.09 2.20
0.51 0.17 1.08 0.25 0.12 0.48
0.30 0.17 1.29 -0.10 0.13 -1.10
Intercept 2.01 0.07 1.02 2.02 0.05 1.19
0.50 0.11 0.58 0.25 0.08 0.45
0.30 0.11 1.29 -0.10 0.08 2.17
Intercept 2.01 0.05 0.41 2.01 0.04 0.59
0.50 0.08 0.71 0.25 0.06 0.57
0.30 0.07 -0.82 -0.10 0.06 0.01
Intercept 2.03 0.13 2.21 2.05 0.10 3.79
0.51 0.22 2.83 0.25 0.14 1.16
0.30 0.21 0.14 -0.10 0.14 0.32
Intercept 2.00 0.08 0.27 2.03 0.06 2.27
0.51 0.13 1.32 0.25 0.09 1.32
0.30 0.13 0.94 -0.10 0.09 -2.83
Intercept 2.00 0.06 0.09 2.02 0.04 1.67
0.50 0.09 0.11 0.25 0.06 0.86
0.30 0.09 0.33 -0.10 0.06 -1.55
The true parameters: .    
Table 2: Simulation: ML Estimates for the MPR Weibull model with various sample sizes and censoring rates – without frailty
Intercept 1.99 0.22 -0.80 2.02 0.12 1.49 0.48 0.40 5.87
0.51 0.23 1.08 0.25 0.12 1.12 - - -
0.30 0.22 -0.38 -0.10 0.12 2.47 - - -
Intercept 2.00 0.14 0.28 2.01 0.07 1.03 0.50 0.24 0.81
0.50 0.14 0.26 0.25 0.07 0.44 - - -
0.30 0.14 0.99 -0.10 0.07 -0.46 - - -
Intercept 2.00 0.10 0.00 2.00 0.05 0.35 0.50 0.17 0.95
0.50 0.10 0.13 0.25 0.05 0.37 - - -
0.30 0.10 1.06 -0.10 0.05 -0.11 - - -
Intercept 2.02 0.28 1.44 2.03 0.14 2.26 0.48 0.53 4.53
0.51 0.27 2.53 0.26 0.14 2.12 - - -
0.30 0.26 1.06 -0.10 0.13 -4.50 - - -
Intercept 2.01 0.18 0.76 2.02 0.09 1.24 0.50 0.32 0.76
0.50 0.17 0.68 0.25 0.09 1.23 - - -
0.30 0.18 -1.43 -0.10 0.09 3.46 - - -
Intercept 2.01 0.12 0.66 2.01 0.06 0.97 0.50 0.23 -0.18
0.51 0.12 1.45 0.25 0.06 0.32 - - -
0.30 0.13 1.39 -0.10 0.06 0.59 - - -
Intercept 1.98 0.25 -1.55 2.01 0.14 0.88 0.48 0.45 7.05
0.50 0.24 -0.15 0.25 0.13 1.80 - - -
0.30 0.22 -0.06 -0.10 0.13 4.62 - - -
Intercept 1.99 0.16 -0.95 2.00 0.09 0.19 0.49 0.27 3.69
0.50 0.15 0.15 0.25 0.08 -0.24 - - -
0.30 0.15 -0.79 -0.10 0.08 0.46 - - -
Intercept 2.00 0.11 -0.16 2.00 0.06 0.02 0.49 0.19 1.68
0.50 0.11 0.52 0.25 0.06 0.70 - - -
0.30 0.10 0.39 -0.10 0.06 -0.33 - - -
Intercept 2.01 0.32 0.41 2.05 0.16 3.65 0.48 0.60 5.41
0.50 0.29 0.96 0.25 0.15 -0.71 - - -
0.31 0.29 3.71 -0.10 0.15 0.55 - - -
Intercept 1.99 0.19 -0.72 2.03 0.10 2.29 0.50 0.36 0.66
0.50 0.18 0.60 0.25 0.09 0.11 - - -
0.30 0.17 1.55 -0.10 0.09 -4.23 - - -
Intercept 2.00 0.13 0.00 2.02 0.07 1.71 0.50 0.25 -0.42
0.50 0.12 0.49 0.25 0.07 0.94 - - -
0.30 0.12 1.52 -0.10 0.07 -1.61 - - -
The true parameters: .    
Table 3: Simulation: ML Estimates for the MPR Weibull model with various numbers of sample sizes and censoring rates – with frailty variance
Intercept 2.01 0.29 0.57 2.03 0.14 2.27 0.46 0.57 12.54
0.48 0.44 -3.53 0.24 0.20 -2.51 0.15 0.78 -0.67
0.32 0.32 5.77 -0.09 0.15 -6.89 -0.18 0.56 -11.71
Intercept 2.00 0.17 -0.11 2.01 0.09 0.62 0.48 0.33 5.47
0.51 0.28 2.04 0.25 0.13 1.44 0.17 0.47 10.98
0.30 0.24 1.39 -0.10 0.11 0.88 -0.20 0.41 -0.68
Intercept 2.00 0.12 0.01 2.00 0.06 0.31 0.49 0.24 3.07
0.50 0.20 0.99 0.25 0.09 0.64 0.15 0.32 2.68
0.31 0.18 2.12 -0.10 0.08 -1.31 -0.20 0.29 -2.08
Intercept 2.11 0.35 7.59 2.06 0.16 4.25 0.48 0.72 5.93
0.50 0.56 0.81 0.25 0.22 1.51 0.13 1.02 -11.36
0.31 0.43 4.02 -0.10 0.18 -2.31 -0.19 0.75 -7.38
Intercept 2.01 0.21 0.50 2.01 0.10 1.04 0.47 0.47 8.34
0.52 0.34 3.76 0.25 0.14 0.40 0.19 0.65 24.02
0.31 0.26 3.92 -0.10 0.11 -2.16 -0.18 0.49 -8.27
Intercept 2.00 0.16 0.17 2.01 0.07 0.45 0.49 0.34 4.34
0.51 0.24 2.36 0.26 0.10 2.18 0.17 0.44 13.61
0.31 0.20 2.51 -0.10 0.09 -0.95 -0.20 0.37 0.24
Intercept 1.99 0.32 -0.69 2.02 0.18 1.11 0.44 0.66 19.10
0.49 0.50 -2.86 0.25 0.24 -0.94 0.19 0.90 29.52
0.31 0.32 1.84 -0.10 0.16 -2.32 -0.18 0.56 -11.54
Intercept 1.98 0.19 -1.66 2.00 0.11 0.10 0.47 0.38 8.38
0.50 0.31 -0.10 0.25 0.15 0.88 0.17 0.52 13.91
0.31 0.23 4.95 -0.09 0.11 -5.28 -0.19 0.39 -4.39
Intercept 2.00 0.14 -0.33 2.00 0.07 0.07 0.49 0.26 4.14
0.50 0.22 -0.81 0.25 0.11 1.43 0.16 0.36 3.93
0.30 0.19 1.49 -0.10 0.09 -1.22 -0.20 0.31 -1.01
Intercept 2.12 0.39 8.46 2.11 0.19 7.63 0.49 0.79 3.67
0.50 0.62 0.90 0.24 0.26 -2.37 0.13 1.09 -13.02
0.34 0.45 14.57 -0.09 0.20 -12.54 -0.15 0.75 -24.64
Intercept 2.01 0.24 0.88 2.04 0.12 2.77 0.48 0.52 7.28
0.52 0.38 3.62 0.25 0.17 1.39 0.17 0.70 12.32
0.31 0.28 2.48 -0.10 0.13 -4.14 -0.18 0.50 -11.03
Intercept 2.01 0.17 10.65 2.03 0.08 2.48 0.49 0.36 1.58
0.50 0.27 -0.80 0.25 0.12 -0.42 0.14 0.48 -3.85
0.31 0.22 2.51 -0.09 0.10 -6.28 -0.20 0.39 -0.80
The true parameters: .    
Table 4: Simulation: ML Estimates for the MPR Weibull model with various numbers of sample sizes and censoring rates – with frailty dispersion

6 Data analysis

We analysed a subset of data from the Signal Tandmobielstudy. Following other authors (e.g., Gómez et al. (2009)), time to emergence, , was measured as “child’s age minus five years” (i.e., age 5) since emergence of pre-molars does not occur before age five; more generally, of course, this threshold parameter could be estimated. Two covariates were analysed, namely, sex, where 0 = boy  (52%) and 1 = girl  (48%), and dmf, where 0 (57%) indicates that the primary predecessor tooth was sound and 1 (43%) indicates that it was decayed, missing due to caries, or filled. Gómez et al. (2009) excluded 44 (1%) school children in whom the dmf status was unknown thus leaving 4386 children for our analysis. It should be noted that a more extensive data set from the Signal Tandmobiel study has been analysed by Bogaerts et al. (2002), Lesaffre and Komárek (2005) and Komárek and Lesaffre (2009).

We investigated four covariate structures: (I) sex only, (II) dmf only, (III) sex and dmf together, and (IV) sex and dmf together along with their interaction term. Each of these covariate structures were included in the six different model types designated by: PH, PHF, PHDM, MPR, MPRF, MPRDM (see Table1). Once the best-fitting model type is identified (using the mean AIC averaged over covariate structures as a guide), covariate selection proceeds within the best type, and the final model selected, may, if MPR, have different covariates in the regressions.

Thus, there are 24 initial models in total defined by the combination of four covariate structures in each of six model types. These models were fitted to the data using a specially written R programme (R Core Team, 2018) which called the routine nlm to maximize the likelihood function (3.6) and compute the observed information matrix. Furthermore, we computed the non-parametric maximum likelihood estimator of the survivor function (Turnbull, 1976) using the R package interval (Fay and Shaw, 2010).

6.1 Results

A summary of the 24 initial models fitted to the data is given in Table 5. Firstly, looking at the four regression structures (I) - (IV), we see that, for all model types, the dmf models, (II), have lower AICs and BICs than the sex models, (I), suggesting that the status of the primary predecessor tooth (sound versus decayed/missing/filled) has a greater effect than the sex of the child. This is not to say that sex is unimportant, as the simultaneous inclusion of both sex and dmf, (III), yields a greater reduction in AIC and BIC relative to the two single-factor models. On the other hand, interestingly, the addition of the interaction term, (IV), reduces the AICs, but increases the BICs, suggesting that there may be a weak interaction effect.

Model Covariates AIC BIC dAIC dBIC
PH(I) sex -5562.1 3 11130.1 11149.3 180.0 147.3
PH(II) dmf -5559.2 3 11124.5 11143.7 174.4 141.7
PH(III) sexdmf -5523.9 4 11055.7 11081.3 105.6 79.3
PH(IV) sexdmf -5520.2 5 11050.3 11082.3 100.2 80.3
mean: 11090.2 11114.1 140.0 112.2
PHF(I) sex -5540.9 4 11089.7 11115.3 139.6 113.3
PHF(II) dmf -5526.6 4 11061.2 11086.7 111.0 84.8
PHF(III) sexdmf -5488.2 5 10986.4 11018.3 36.3 16.4
PHF(IV) sexdmf -5485.1 6 10982.2 11020.5 32.0 18.5
mean: 11029.9 11060.2 79.7 58.3
PHDM(I) sex -5540.8 5 11091.7 11123.6 141.5 121.6
PHDM(II) dmf -5516.3 5 11042.5 11074.5 92.4 72.5
PHDM(III) sexdmf -5475.2 7 10964.4 11009.1 14.3 7.1
PHDM(IV) sexdmf -5472.6 9 10963.3 11020.8 13.2 18.8
mean: 11015.5 11057.0 65.3 55.0
MPR(I) sex -5560.8 4 11129.7 11155.2 179.6 153.3
MPR(II) dmf -5538.3 4 11084.6 11110.2 134.5 108.2
MPR(III) sexdmf -5501.7 6 11015.4 11053.7 65.3 51.8
MPR(IV) sexdmf -5493.7 8 11003.4 11054.4 53.2 52.5
mean: 11058.3 11093.4 108.1 91.4
MPRF(I) sex -5540.7 5 11091.4 11123.4 141.3 121.4
MPRF(II) dmf -5511.3 5 11032.6 11064.5 82.5 62.6
MPRF(III) sexdmf -5471.6 7 10957.2 11001.9 7.1 0.0
MPRF(IV) sexdmf -5466.1 9 10950.1 11007.6 0.0 5.7
mean: 11007.8 11049.4 57.7 47.4
MPRDM(I) sex -5540.7 6 11093.3 11131.7 143.2 129.7
MPRDM(II) dmf -5511.2 6 11034.4 11072.7 84.3 70.8
MPRDM(III) sexdmf -5469.8 9 10957.6 11015.1 7.5 13.2
MPRDM(IV) sexdmf -5465.6 12 10955.2 11031.9 5.1 29.9
mean: 11010.2 11062.8 60.0 60.9
“sex” and “dmf” indicate a single-factor model in which one of sex or dmf appears; “sexdmf” indicates a model with both sex and dmf; “sexdmf” indicates a models with both sex and dmf along with the interaction between these two; note that all models contain an intercept terms in the scale and shape; is the log-likelihood value; is the number of parameters in the model; “AIC” and “BIC” are the Akaike Information Criteria and Bayesian Information Criteria respectively; where represents the lowest AIC among the models shown in this table which is that of model MPRF(IV) whose AIC is ; where corresponds to model MPRF(III) whose BIC is 11001.9. Note that, for each of the six model types, the mean AIC, BIC, dAIC, and dBIC values are shown to facilitate quick comparison of the basic model types.
Table 5: Summary of initial models fitted

We now discuss the merits of the various model types. First we notice in Table 5 that generally the MPR models outperform their simpler PH counterparts in terms of AIC and BIC (the only exceptions being some MPRDM versus PHDM comparisons). Visually, this improvement in model fit is clear by comparing Figures 1 (a) and (b) to Figures 1 (c) and (d). This shows the additional value gained by modelling the shape, and highlights that one or more of the covariates have non-PH effects.

The addition of frailty to PH or MPR models improves the fit in all cases. This suggests that additional heterogeneity (e.g., via unobserved covariates) exists within the data. That this is so for the MPR models is noteworthy, since the MPR model already explains variation beyond the PH model by means of its person-specific shape regression – but, it is clear, that there is additional heterogeneity present in these data.

  (a) PH(IV) – Boys   (b) PH(IV) – Girls
  (c) MPR(IV) – Boys   (d) MPR(IV) – Girls
  (e) MPRF(IV)R – Boys   (f) MPRF(IV)R – Girls
Figure 1: Non-parametric Turnbull survivor curves (solid) with model-based curves (dash) overlayed for PH(IV) (panels (a) and (b)), MPR(IV) (panels (c) and (d), and MPRF(IV)R (panels (e) and (f)). The curves for boys are on the left, while the curves for boys are to the right. Within each plot we show the sound and decayed predecessor tooth groups (i.e., non-dmf and dmf respectively).

It is noteworthy that the inclusion of dispersion models (DM) is supported for the PH model cases, but not for the MPR cases. From (2.3), we see that plays a role in describing the shape of the marginal distribution (albeit that is the primary shape parameter). Thus, the support for DM on top of a basic PH model essentially highlights the benefit of modelling shape in addition to scale as advocated by Burke and MacKenzie (2017). In particular, the comparison of PHDM against MPRF allows us to assess the value of modelling shape via the frailty variance, , or the shape parameter, . The latter models outperform the former in terms of AIC and BIC suggesting that there is more value in modelling the shape than the frailty variance, at least, in these data. This is perhaps not surprising since is the main shape parameter as mentioned.

Overall, we find that the best model in terms of AIC is MPRF(IV). This is the MPRF model with and regression components both containing sex, dmf, and the sexdmf interaction. The best model in terms of BIC is the MPRF(III) model, obtained by omitting the interaction term from the regression components of MPRF(IV). Gómez et al. (2009) did not consider a non-PH process for time to emergence of tooth24. Here, both of these best-fitting models are non-PH, and the non-PH-ness arises in two ways: (a) the shape parameter, , depends on covariates, and, (b) through the presence of frailty. The coefficients for these two models are shown in Table 6 where we see that, in MPRF(III), the sex effect appears to be statistically significant in neither the scale nor shape. A naive approach to variable selection might treat the scale and shape regression components completely separately, thereby removing sex entirely from the model. Such a removal would bring us back to MPRF(II) which, from Table 5, has much higher AIC and BIC values. This highlights the more involved nature of variable selection within MPR models, i.e., we cannot simply rely on separate Wald tests for (potentially highly) correlated scale and shape effects (see Burke and MacKenzie (2017) for further details). In this case, it is clear that the sex effect is required, but eliminating it from the shape yielded model MPRF(III)R which has improved AIC and BIC; note that eliminating sex from the scale instead produces similar AIC and BIC values (not shown). Similarly, careful reduction of MPRF(IV) yielded MPRF(IV)R which has the lowest AIC and BIC of all models considered. Both of these reduced models are shown in Table 6 and, of course, the reduction produces more interpretable models.

MPRF(IV)R is the best-fitting model by both AIC and BIC measures, and we can confirm that the fit is near-perfect via Figures 1 (e) and (f). Before we interpret this model however, we first highlight PH(IV) shown in Table 6 which is the simplest full-covariate model fitted. The scale coefficients of this model suggest that tooth24 emerges first in girls with dmf, then girls without dmf and boys with dmf (whose emergence times are close), and, lastly, boys without dmf. On the other hand, model MPR(IV), also shown in Table 6, directly extends this model and improves the fit which suggests that in fact something more complex arises. Returning to MPRF(IV)R, because this has a reduced shape covariate structure (compared with MPR(IV)), we can see that dmf has a highly non-PH effect due to its appearance in the shape component.

Scale Shape Frailty Scale Shape Frailty
intercept -12.97 (0.42) 1.98 (0.04) -0.45 (0.15) -13.68 (0.54) 2.03 (0.04) -0.48 (0.16)
girl 0.19 (0.30) 0.02 (0.03) 1.52(0.59) -0.07 (0.05)
dmf 2.73(0.38) -0.19(0.03) 3.36(0.57) -0.23(0.05)
girldmf -1.08 (0.76) 0.06 (0.06)
Scale Shape Frailty Scale Shape Frailty
intercept -13.05 (0.42) 1.99 (0.03) -0.46 (0.15) -13.22 (0.43) 1.99 (0.03) -0.46 (0.15)
girl 0.47(0.06) 0.62(0.08)
dmf 2.65(0.37) -0.18(0.03) 2.93(0.39) -0.19(0.03)
girldmf -0.33 (0.11)
Scale Shape Frailty Scale Shape Frailty
intercept -9.95 (0.16) 1.68 (0.02) -11.82 (0.39) 1.86 (0.03)
girl 0.43(0.05) 1.64(0.50) -0.11(0.04)
dmf 0.45(0.06) 3.08(0.48) -0.26(0.05)
girldmf -0.21(0.08) -1.08 (0.63) 0.07 (0.06)
“girldmf” represents the interaction between being a girl and having dmf (primary predecessor tooth decayed / missing / filled); MPRF(III)R is a reduced form of MPRF(III) and, similarly, MPRF(IV)R is a reduced form of MPRF(IV); AIC, BIC, dAIC and dBIC are as described in Table 5. Note that regression coefficients with asterisk are statistically significant at 5% level according to the Wald test (but we do not highlight statistically significant intercepts).
Table 6: Selection of fitted models including additional reduced models

Figures 2 (a) and (b) show the hazard functions for the four groups defined by the sexdmf interaction for PH(IV) and MPRF(IV)R, while Figures 2 (c) and (d) display hazard ratios relative to the boys without dmf group. We can see that the MPRF(IV)R model permits quite different hazard shapes for each group, whereas, for the PH(IV) model the shapes are constrained to be the same; in both cases the hazards increase with time, indicating the inevitability of the emergence of tooth24 later in time. Overall, both models agree in terms of the highest and lowest emergence hazards which correspond, respectively, to girls with dmf and boys without dmf groups. On the other hand, within the MPRF(IV)R model, it seems that boys with dmf have a higher hazard than girls without dmf earlier in time which reverses later in time – at about time point 5 (i.e., when the children are aged 10); this crossing effect cannot be handled by the simpler PH(IV) model which forces these two groups to be equal. Table 7 shows the estimated median emergence times for the four groups under the PH(IV) and MPRF(IV)R models; the results are in line with the hazard-based interpretations.

  (a) PH(IV) Hazard Functions   (b) MPRF(IV)R Hazard Functions
  (c) PH(IV) Hazard Ratios   (d) MPRF(IV)R Hazard Ratios
Figure 2: Hazard functions (panels (a) and (b)) and hazard ratios (panels (c) and (d)) for boys without dmf (solid), girls without dmf (dash), boys with dmf (dot), and girls with dmf (dash-dot) for the PH(IV) (panels (a) and (c)) and MPR(IV)R (panels (b) and (d)) models respectively. Note that, for the hazard ratios, we are using boys without dmf (the lowest hazard group) as the reference group for the four groups formed via the sexdmf interaction.
Yes Girls 5.27  [5.19 5.34] 5.10  [5.01 5.20]
Yes Boys 5.49  [5.40 5.57] 5.35  [5.26 5.44]
No Girls 5.51  [5.43 5.58] 5.49  [5.41 5.57]
No Boys 5.97  [5.88 6.05] 5.98  [5.89 6.06]
Recall that we have modelled years so that adding 5 to the above emergence times gives the emergence age.
Table 7:

Estimated median emergence times with 95% confidence interval

While the nature of the covariate effects can be determined by examining Figure 2 it is instructive to consider the effects of dmf and sex separately (albeit they do interact). Thus, Figures 3 (a) and (b) present the dmf hazard ratios (i.e., dmf versus non-dmf) in boys and girls. We can see that the effect of dmf, which is highly time-dependent and greater in boys, is to increase the hazard of emergence – although, later in time, the strength of this effect reduces. Figures 3 (c) and (d) present the sex hazard ratios (i.e., girls versus boys) in non-dmf and dmf groups. We can see that the girls have a greater hazard of emergence than boys; the effect is reduced when dmf is present. Compared with the dmf hazard ratios, we can see that the sex effect is weaker (as was apparent from Table 5) and is much less time-dependent due to the lack sex in the shape regression (in fact the MPRF(IV)R sex hazard ratios are much closer to their PH(IV) counterparts).

  (a) Dmf hazard ratio in boys   (b) Dmf hazard ratio in girls
  (c) Sex hazard ratio in non-dmf   (d) Sex hazard ratio in dmf
Figure 3: Ratio of dmf to non-dmf marginal hazards in boys (panel (a)) and girls (panel (b)), and ratio of girls to boys marginal hazards in non-dmf (panel (c)) and dmf (panel (d)) groups. These marginal hazard ratios are shown for both the PH(IV) (dash) and MPR(IV)R (solid) models; a reference line at one is also shown.

7 Discussion

In this paper we investigated the utility of MPR models in the context of interval censored data, by way of simulation and a practical application, which, to the best of our knowledge, is the first time such an investigation has been carried out. In particular, we have found that the parameters can be estimated with reasonable precision even in relatively small samples of interval censored data (albeit the most complicated model structures work best in larger samples). Moreover, we found that the MPR extension, and the additional extensions of frailty and dispersion modelling, to be fruitful in the context of the Signal Tandmobiel study. Thus, as we might expect, the utility of the MPR Weibull model, and indeed the MPR framework in general, extends beyond right-censored data and the specific lung cancer application considered in Burke and MacKenzie (2017)

Our analysis of the tooth data considered a variety of additional model structures not previously explored in the existing literature (Bogaerts et al., 2002; Gómez et al., 2009; Lesaffre and Komárek, 2005; Komárek and Lesaffre, 2009). It is noteworthy that Lesaffre and Komárek (2005), who developed a spline-based AFT model, suggest that “parametric methods do not offer enough flexibility to correctly model survival data”. In contrast, we have demonstrated the appropriateness of relatively simpler parametric models for these data which achieve flexibility through a combination of multi-parameter regression and frailty modelling. Interestingly, Lesaffre and Komárek (2005) and Komárek and Lesaffre (2009) briefly considered a dispersion model extension of their AFT model which they refer to as a mean-scale model (citing earlier work by Pan and MacKenzie (2003) on mean-covariance modelling in longitudinal studies). While an investigation of MPR approaches was not their focus, they nonetheless found, like us, that an MPR extension yielded a superior fit to the data.

It is fortunate that these data are sufficiently extensive to permit the investigation of such models. We find that to emergence of the permanent upper left first premolars depends on sex and dmf. In particular, emergence times are significantly earlier in children whose predecessor tooth was decayed, missing, or filled (dmf), and that emergence times are earlier in girls. However, we have also found that frailty effects are supported within the data, i.e., there may be further unmeasured features at play. The time-dependent nature of the dmf hazard ratios is quite interesting, and suggests that the dmf group becomes more like the non-dmf group later in time. Recall that the dmf group is a mixture of individuals with decayed, missing, and filled teeth. With this in mind, one might speculate that filled teeth are more similar to sound teeth (non-dmf teeth) with larger emergence times, while missing teeth are perhaps quite different with shorter emergence times. This, at least, would be a frailty interpretation of time varying effects, i.e., there is a mixture of groups, some of which “fail” earlier. However, the time variation is much greater than that supported by frailty alone due to the presence of the significant dmf shape effect within the MPRF(IV)R model.

The extension of the MPR framework to interval-censored survival data with frailty permits the examination of a variety of potential data structures. An appealing aspect of this approach is that the breadth of models supported exist within a reasonably straightforward parametric setup which is not computationally intensive. Such a framework provides a practical and useful adjunct to existing methods which may reveal new insights.


This research was supported by the Irish Research Council (New Foundations Award).


  • Al-tawarah and MacKenzie (2002) Al-tawarah, Y. and MacKenzie, G. (2002). A logistic ph regression model for interval censored survival data. Proceedings of the 17th International Workshop on Statistical Modelling. Chania, Crete Eds. M. Stasinopoulis and Goito Touloumi., 17:105–113.
  • Bogaerts et al. (2017) Bogaerts, K., Komarek, A., and Lesaffre, E. (2017). Survival Analysis with Interval-Censored Data: A Practical Approach with Examples in R, SAS, and BUGS. Chapman and Hall/CRC.
  • Bogaerts et al. (2002) Bogaerts, K., Leroy, R., Lesaffre, E., and Declerck, D. (2002). Modelling tooth emergence data based on multivariate interval-censored data. Statistics in Medicine, 21:3775 –3787.
  • Burke and MacKenzie (2017) Burke, K. and MacKenzie, G. (2017). Multi-parameter regression survival modelling - an alternative to proportional hazards. Biometrics, 73:361–716.
  • Duchateau and Janssen (2008) Duchateau, L. and Janssen, P. (2008). The Frailty Model. Springer.
  • Fay and Shaw (2010) Fay, M. P. and Shaw, P. A. (2010). Exact and Asymptotic Weighted Logrank Tests for Interval Censored Data: The interval R package. Journal of Statistical Software. pp 1-34. 36, (2).
  • Finkelstein (1986) Finkelstein, D. M. (1986). A proportional hazards model for interval-censored failure time data. Biometrics, 42:845 –854.
  • Gómez et al. (2009) Gómez, G., Calle, M. L., Oller, R., and Langohr, K. (2009). Tutorial on methods for interval-censored data and their implementation in r. Statistical Modelling, 9,(4):259 –297.
  • Hart (2002) Hart, P. M. et al.. (2002). The subfoveal radiotherapy study: Visual outcome at 12 and 24 months in a randomised controlled trial of external beam radiotherapy for subfoveal choroidal neovasularisation of age-related maculuar degeneration. Archives: 120., 8:1029 –1038.
  • Hougaard (1995) Hougaard, P. (1995). Frailty models for survival data. Lifetime Data Analysis, 1, 3:255 –273.
  • Hougaard (2000) Hougaard, P. (2000). Analysis of multivariate survival Data. Springer-Verlag, New York.
  • Huang and Wellner (1997) Huang, J. and Wellner, J. (1997). Interval censored survival data: A review of recent progress. Proceedings of the First Seattle Symposium in Biostatistics, Springer, New York, NY, 123:123 –169.
  • Komárek and Lesaffre (2009) Komárek, A. and Lesaffre, E. (2009).

    The regression analysis of correlated interval-censored data: illustration using accelerated failure time models with flexible distributional assumptions.

    Statistical Modelling, 9,4:299 –319.
  • Lawless (2003) Lawless, J. (2003). Statistical Models and Methods for Lifetime Data. Wiley, second edition.
  • Lee and Nelder (2001) Lee, Y. and Nelder, J. (2001). Hierarchical generalised linear models: A synthesis of generalised linear models, random effect models and structured dispersions. Biometriks, 88, 4:987 –1006.
  • Lesaffre and Komárek (2005) Lesaffre, E. and Komárek, A. (2005). An overview of methods for interval-censored data with an emphasis on applications in dentistry. Statistical Methods in Medical Research, 14:539 –552.
  • Lynch and MacKenzie (2014) Lynch, J. and MacKenzie, G. (2014). Frailty models with structural dispersion. In MacKenzie, G. and Peng, D., editors, Statistical Modelling in Biostatistics and Bioinformatics, pages 35 –44. Springer, Munich Germany. ISBN= 978-3-319-04578-8.
  • MacKenzie (1999) MacKenzie, G. (1999). Survival analysis for longitudinal data. Proceedings of the 14th International Workshop on Statistical Modelling. Graz, Austria., 14:259–264.
  • MacKenzie and Peng (2013) MacKenzie, G. and Peng, D. (2013). Interval-censored parametric regression survival models and the analysis of longitudinal trials. Statistics in Medicine., 32:2804–2822.
  • Marshall and Olkin (2007) Marshall, A. W. and Olkin, I. (2007). Life Distributions: Structure of Nonparametric, Semiparametric, and Parametric Families. Springer.
  • Pan and MacKenzie (2003) Pan, J. and MacKenzie, G. (2003). On modelling mean-covariance structures in longitudinal studies. Biometrika, 90, 1:239–244.
  • R Core Team (2018) R Core Team (2018). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
  • Smyth (1989) Smyth, G. (1989). Generalized linear models with varying dispersion. JRSS, Series B (Methodological), 51,1:47 –60.
  • Sun (2006) Sun, J. (2006). Statistical Analysis of Interval-Censored Failure Time Data. Springer, New York, 1st edition.
  • Turnbull (1976) Turnbull, B. (1976). The empirical distribution function with arbitrarily grouped, censored and truncated data. J. Royal Statist. Soc. B, 38, 3:290–295.
  • Vaupel et al. (1979) Vaupel, J., Manton, K., and Stallard, E. (1979). The impact of heterogeneity in individual frailty on the dynamics of mortality. Demography, 16:439–454.
  • Wilkinson (1995) Wilkinson, P. (1995). Lung Cancer in Northern Ireland 1991-92. MD Thesis, The Queen’s University of Belfast; Vol 1, pp 1-204.
  • Zhang (2009) Zhang, Z. (2009). A class of transformed regression models for interval censoring. Statistical Modelling, 9,4:321 –343.