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 followup 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 followup examinations at times , . Thus, if, for the th subject, the event occurs between the th and the th followup 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 followups 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 nonPH and they compared the use of standard right censored likelihoods based on midpoints with interval censored likelihoods. They showed that the use of midpoints led to artificially precise estimators in PH models when analysing time to loss of vision in a longitudinal trial of AgeRelated Macular Degeneration (ARMD)
(Hart, 2002). Survival data arising in longitudinal vision studies, were analysed earlier, by MacKenzie (1999) and later by Altawarah 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 multiparameter 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).
2 MPR modelling framework
We envisage the class of twoparameter 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 intervalcensored 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
(2.1) 
which arises from a basic (nocovariate) Weibull hazard function, , (Lawless, 2003; Marshall and Olkin, 2007) allowing both its scale, , and shape, , parameters to depend on covariates via the MPR specification
(2.2) 
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 nonPH, 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 nonPH effects, including crossing hazards and, through the coefficients, provides covariatespecific 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 personspecific 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
(2.3) 
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 (nonfrailty) 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 personspecific via another regression model, i.e.,
(2.4) 
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 meandispersion 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 followow examinations is too rigid as many subjects fail to respect their exact reexamination appointment dates. Accordingly it is usual to allow the intervals to be personspecific 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 intevalcensoring (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 rightcensored. In total, there are patients of whom are right censored or withdrawn at specific times, leaving patients who are intervalcensored. Note that the intervalcensored subjects play the same role as “failures” in the rightcensored setting and often rightcensoring occurs at times completely unrelated to the scheduled followup examinations, e.g., an early withdrawal from the study. Thus, this notational setup is advantageous when it is important to distinguish between intervalcensored and rightcensored observations. See MacKenzie and Peng (2013) for more details on this approach.
Here, however, we rewrite the likelihood above as
(3.5) 
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
(3.6) 
where , and are the cumulative hazard functions for the th individual evaluated at the endpoints 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 loglikelihood 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 Weibullgamma 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 multiparameter 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.
Regression  
Model  Baseline  Frailty  Yes  No 
PH  PH  No  
PHF  PH  Yes  
PHDM  PH  Yes  
MPR  MPR  No  —  
MPRF  MPR  Yes  
MPRDM  MPR  Yes  — 
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 nonPH covariate, and, in a frailty dispersion model, indicates that the frailty variance differs in the subgroups 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 Waldbased significance tests (which account only for the variance of estimates, and not covariance – which is important in this context) might render a particular covariate nonsignificant 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 intervalcensored 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 noninformative 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 . Letwhere 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 rightcensoring and inspection tend to reduce performance, but, even in the worst case of 30% rightcensoring 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.
Estimates  
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 
Estimates  
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       
Estimates  
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 
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 premolars 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 bestfitting 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 nonparametric 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 singlefactor 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 
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 nonPH 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 personspecific 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 
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 nonPH process for time to emergence of tooth24. Here, both of these bestfitting models are nonPH, and the nonPHness 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 bestfitting model by both AIC and BIC measures, and we can confirm that the fit is nearperfect 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 fullcovariate 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 nonPH effect due to its appearance in the shape component.
MPRF(III)  MPRF(IV)  
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)  — 
MPRF(III)R  MPRF(IV)R  
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)  —  — 
PH(IV)  MPR(IV)  
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)  — 
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 hazardbased interpretations.
(a) PH(IV) Hazard Functions  (b) MPRF(IV)R Hazard Functions 
(c) PH(IV) Hazard Ratios  (d) MPRF(IV)R Hazard Ratios 
Dmf  Sex  PH(IV)  MPRF(IV)R 
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] 
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 nondmf) in boys and girls. We can see that the effect of dmf, which is highly timedependent 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 nondmf 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 timedependent 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 nondmf  (d) Sex hazard ratio in dmf 
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 rightcensored 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 splinebased 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 multiparameter 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 meanscale model (citing earlier work by Pan and MacKenzie (2003) on meancovariance 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 timedependent nature of the dmf hazard ratios is quite interesting, and suggests that the dmf group becomes more like the nondmf 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 (nondmf 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 intervalcensored 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.
Acknowledgements
This research was supported by the Irish Research Council (New Foundations Award).
References
 Altawarah and MacKenzie (2002) Altawarah, 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 IntervalCensored 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 intervalcensored data. Statistics in Medicine, 21:3775 –3787.
 Burke and MacKenzie (2017) Burke, K. and MacKenzie, G. (2017). Multiparameter 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 134. 36, (2).
 Finkelstein (1986) Finkelstein, D. M. (1986). A proportional hazards model for intervalcensored 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 intervalcensored 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 agerelated 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. SpringerVerlag, 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 intervalcensored 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 intervalcensored 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= 9783319045788.
 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). Intervalcensored 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 meancovariance 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 IntervalCensored 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 199192. MD Thesis, The Queen’s University of Belfast; Vol 1, pp 1204.
 Zhang (2009) Zhang, Z. (2009). A class of transformed regression models for interval censoring. Statistical Modelling, 9,4:321 –343.
Comments
There are no comments yet.