Original Research

A conditional frailty model for bivariate interval-truncated failure time data: an application to a study on siblings diagnosed with schizophrenia

© 2016 Sutradhar et al; licensee Herbert Publications Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

---

Introduction

Interval truncated failure time data arise when only those individuals with a failure time that lies within a certain interval are observed. When a failure time falls outside this truncation interval, no information regarding the subject is available. One may view this as a screening procedure in which all observations outside the interval known as the truncation region are removed. It is the concept of "who is included" that distinguishes truncation from censoring. In the case of censoring, we have partial information on all observations in a sample in that subjects are known to have failure times greater (or less) than some value. However in the case of truncation, individuals with failure times outside the truncation interval are not considered for inclusion in the study. This important difference between truncation and censoring is reflected in the construction of the data likelihood.

Approaches for dealing with censoring in the presence of interval truncation were first discussed by Turnbull [1]. He proposed a method for obtaining the nonparametric maximum likelihood estimate of the distribution function under interval truncation that also accommodates left, right, and interval censoring. Turnbull described a self-consistent approach to estimation that turns out to be a version of the EM algorithm [2]. Frydman noted that the method proposed by Turnbull was not applicable in the case of both truncation and censoring, and suggested appropriate modifications [3]. Alioum and Commenges provided a further detailed correction of Turnbull's approach and offered an extension to regression analysis [4].

We are interested developing methodology for examining interval-truncated failure times in a bivariate setting. That is, when available information consists of paired failure times (e.g., we consider families of size two), where, due to the observation scheme, paired failure times are known only if they both fall within their respective truncation intervals. Such data arises in settings where families are observed through family disease registries, for example, or where data are grouped by family by design. When paired data are collected subject to truncation, it is frequently of interest to investigate associations among failure times within families. The existing literature on methods for analyzing interval-truncated failure times in a bivariate setting is very limited and mostly refers to the case when only one component of the bivariate vector is subject to truncation. Exceptions to this assumption include the studies by van der Laan and Gurler [5,6], in which nonparametric estimation of the bivariate survival function is discussed where both components of the bivariate vector are randomly left or right-truncated, but not interval-truncated. This paper proposes a model and estimation approach for examining interval-truncated bivariate failure time data. Straight forward extensions are possible to deal with the multivariate setting.

Methods

Notation
Suppose a pair of individuals is included in a study if both members' event times are observed to fall in the calendar time interval [L,R]. Let Y_ij represent the year of birth for the j^th member in the ith family and let X_ij denote the year of event for that individual. Let the age at the year of the event for the jth subject of the ith pair be denoted by T_ij=X_ij-Y_ij, and let B_ij=[L_ij,R_ij]=[L-Y_ij,R-Y_ij] be the corresponding truncation interval in terms of age, where L_ij is the age at calendar time L and R_ij is the age at calendar time R. Although we observe T_ij if it belongs to its corresponding truncation interval B_ij, most information regarding familial association is available if two or more event times within a family meet truncation requirements.

For the i^th family with two members or siblings j=1,2, under the assumption of independent truncation, the joint failure times are sampled from a conditional distribution.

F( t i1 , t i2 | T i1 ∈ B i1 , T i2 ∈ B i2 )=P( T i1 ≤ t i1 , T i2 ≤ t i2 | T i1 ∈ B i1 , T i2 ∈ B i2 ). (1) (1)

The likelihood for bivariate truncated failure times with observed data (t_ij, B_ij=[l_ij,r_ji]; j=1,2, i=1,...n) can be expressed simply in terms of the bivariate survivor function:

L= ∏ i=1 n f( t i1 t i2 | T i1 ∈ B i1 , T i2 ∈ B i2 )= ∏ i=1 n { 1 P( T i1 ∈ B i1 , T i2 ∈ B i2 × ∂ 2 ∂ t i1 ∂ t i2 S * ( t i1 , t i2 ) }, (2) (2)

where

S * ( t i1 , t i2 )=P( T i1 ≥ t i1 , T i2 ≥ t i2 )

and

P( T i1 ∈ B i1 , T i2 ∈ B i2 )= S * ( l i1 , l i2 )− S * ( l i1 , r i2 )− S * ( r i1 , l i2 )+ S * ( r i1 , r i2 ),

is the probability of inclusion in the truncation rectangle Bi=B_i1XB_i2.

Model formulation: joint distributions using a frailty model
A conditional formulation is a common approach adopted for modeling bivariate failure time data. Joint distributions are routinely formulated by assuming responses are conditionally independent given a common scalar random effect. Integrating the conditional joint distribution with respect to the unobserved random effect provides the joint distribution of the failure times. Suppose (T₁,T₂) are bivariate failure times for a particular pair. The joint survivor function under a conditional model can be written as

S * ( t 1 , t 2 )=P( T 1 > t 1 , T 2 > t 2 )= ∫ ​ S 1 ( t 1 |α) S 2 ( t 2 |α )dG( α;ϕ ), (3) (3)

and the marginal survivor function for T_j is

S j * (t)=P( T j >t )= ∫ ​ S j (t|α)dG( α;ϕ ),j=1,2. (4) (4)

Paired failure times T₁ and T₂ are independent given the pairspecific random effect α, which has distribution G(.) indexed by parameter φ.

Various models have been proposed in the literature for the distribution of the random effects: one-parameter gamma distribution [8], the log-normal distribution [13], the positive stable distribution [9], and the inverse gamma distribution [10]. We assume the random effects arise from a one-parameter gamma distribution with mean 1 and variance φ^-1. Furthermore, consider a frailty model in which random effects act multiplicatively on the hazard rate for each member of a family. Conditional on the random effect, the survivor function for T_j is

S j ( t|α )= e −α Λ 0j * (t)

where ^*_0j (t) is the cumulative baseline hazard function for T_j (j=1,2) under the conditional model. These assumptions provide closed form expressions for the joint and marginal survivor function, respectively.

S * ( t 1 , t 2 )= ∫ ​ exp{− αΛ 01 * ( t 1 )}exp{− αΛ 02 * ( t 2 )} ϕ ϕ α ϕ−1 e −ϕα Γ(ϕ) dα= ( 1+ Λ 01 * ( t 1 ) ϕ + Λ 02 * ( t 2 ) ϕ ) −ϕ (5) (5)

and

S j * ( t )= ∫ ​ exp{ − αΛ 0j * ( t ) } ϕ ϕ α ϕ−1 e −ϕα Γ( ϕ ) da= ( 1+ Λ 0j * ( t ) ϕ ) −ϕ ,j=1,2, (6) (6)

Note that the marginal survivor function S_j*(t), under the conditional formulation, depends on both the parameter φ and the parameters of the conditional baseline hazard function Λ_0j*(t). Moreover, as ϕ→∞, S_j*(t) approaches e^{-Λ0j* (t)} for j=1,2 and S*(t₁,t₂) approaches.

e −( Λ 01 * ( t )+( Λ 02 * ( t ))

Thus, as the variance of the random effect distribution tends to zero, S* (t₁,t₂ )= S₁* (t₁)S₂* (t₂), indicating that T₁ and T₂ are independent.

It is important to interpret the parameters of the model appropriately. Parameter φ measures both the lack of fit of the conditional hazard function, as well as the association among paired failure times. Large values of φ could indicate that standard survival models may be inadequate by failing to sufficiently account for variation in the data; this can arise simply from poor fit of a model. That is, large values of φ may be an indication of unaccounted for dependence, overdispersion, or other forms of model misspecification including violations of the proportional hazards assumption. Moreover, as the variance approaches zero, failure times among individuals within a pair approach independence.

Proposed estimation approach: an EM algorithm
Methodological complications arise when both failure times in a pair are observed subject to interval truncation. The inclusion of shared random effects among individuals within a pair, in addition to the presence of interval truncation for each member in the pair, results in a non-trivial likelihood function. Direct maximization of this likelihood function requires integrating the conditional joint distribution with respect to the unobserved random effect, which is not straightforward particularly if strong distributional assumptions are to be avoided. In addition, if a piecewise constant model is used, direct maximization becomes more challenging as the number of pieces (and hence parameters) increases. Instead of trying to directly maximize over a high dimensional space, the EM algorithm offers an alternative method of estimation that is sometimes slow to converge but avoids occasional divergence problems. We develop an EM algorithm for bivariate interval-truncated failure times under a conditional (random effect) formulation. We use the concept of "ghosts", which was first discussed by Turnbull [1] with the aim of obtaining a nonparametric estimate of a truncated distribution in the univariate setting. Here the notion of "ghost pairs" is used for the analysis of interval-truncated failure times in the bivariate case.

Let t_i=(t_i1,t_i2) denote the vector of observed failure times for the subjects of the ith pair and let B_i=(B_i1XB_i2)=[L_i1,R_i1]X[L_i2, R_i2] be the corresponding truncation rectangle. The observed data is represented as X=(X₁,....,,X_n) with X_i=(t_i,B_i) being the observed data for the ith pair (i=1,...,n). To define the complete data, we focus on the complement of the bivariate truncation rectangle B_i^c. Figure 1 provides an illustration of the data for a particular pair of individuals. The truncation region is represented by the center rectangle and the point t_i is observed because both T_i1∈B_i1 and T_i2∈B_i2.

Due to truncation, the pair giving the observations t_i=(t_i1,t_i2) may be considered a remnant of an unknown number of pairs from the same "birth cohort". Pairs are said to be in the same birth cohort if they have the same dates of birth. We define "ghost pairs" corresponding to pair i as pairs of individuals in the same birth cohort who were not observed because their failure times did not fall in the truncation region B_i. Thus these "ghost pairs" would have been observed had there been no truncation.

To define the complete data, let J_i be the number (unknown) of ghost pairs for the i^th pair of observed event times. Suppose u_i=(u_i1,u_i2) for j=1,...,J_i are the corresponding unknown failure times for these ghost pairs. These are the paired times that fall outside the truncation rectangle B_i for the i^th observed pair, and are hence considered to arise from the unobserved ghosts. Times within an observed pair or ghost pair are considered independent given their corresponding random effect, arising from the distribution G(.) indexed by
parameter φ. Here the complete data is given by Y=(Y₁....,Y_n) with Y_i=(X_i,α_i,J_i,{u_j,α_j;j=1,...,J_i}) being the complete data for the ith pair. Notice here we consider the times of the events for the ghost pairs, as well as the random effects for the observed pair and corresponding ghost pairs, as part of the complete data. The reason for this is that the logarithm of the resulting complete data likelihood is relatively convenient to work with and, as will be seen, maximization is particularly straight forward.

The complete data likelihood and corresponding loglikelihood for the ith pair of individuals are constructed as

L ci ( θ; Y i )=f( t i1 | α i )f( t i2 | α i )g( α i )× ∏ j=1 J i {f( u j1 | α j )f( u j2 | α j )g( α j )} (7) (7)

and

l ci ( θ; Y i )=log f( t i1 | α i )+log f( t i2 | α i )+logg( α i )+ ∑ j=1 J i {logf( u j1 | α j )+logf( u j2 | α j )+logg( α j )}, (8) (8)

respectively, where

f( t|α )=α λ 0 * ( t )exp{−α ∫ 0 t λ 0 * ( w )dw}. (9) (9)

We further assume a gamma distribution for the random effect distribution with p.d.f. given as

g( α )= ϕ ϕ Γ(ϕ) α ϕ−1 e −ϕα ,

where E(α)=1 and Var (α)=φ^-1. Based on these assumptions, the complete data log-likelihood can now be re-expressed -the parameters of this function consist of the parameters of the baseline hazard functions and the parameter of the random effects distribution. These parameters will be estimated via an Expectation-Maximization algorithm, where a weakly parametric approach will be taken by assuming a piecewise constant model for the baseline hazard function. The mathematical details of this algorithm are provided in the Appendix. We show in the maximization step that the parameter estimates of the baseline hazard function have a closed simple form and the parameter estimate of the random effects distribution can be easily obtained using standard Newton-Raphson procedures.

Evaluating performance of proposed methodology via simulations
We conducted simulations for investigating the performance of our proposed EM algorithm. We generated random effects α₁,...,α_n from a gamma distribution with p.d.f.

g( α )= ϕ ϕ α ϕ−1 e −ϕα Γ(ϕ) ,

where E(α)=1 and Var(α)=φ^-1. Bivariate failure times (T₁,T₂) were simulated from the joint distribution

S * ( t 1 , t 2 )= ∫ ​ S( t 1 |α )S( t 2 |α )dG(α),

with

S( t|α )= e −α (ψt) γ ,

where the cumulative baseline hazard function Λ₀* (t)=(ψt)^γ had a Weibull form with shape parameter γ and scale parameter ψ^-1. We generated n=1000 pairs of failure times. The median of the marginal distribution

S * (t)= ∫ ​ S(t|α)dG(α)= [1+ ( ψt ) γ /ϕ] −ϕ

was specified at n=40. Note that T₁ and T₂ were taken to have the same marginal distribution. The true value of the parameter φ was assumed to be 1/0.25 and 1/0.75 (log φ=1.386 and 0.288, respectively), and the true value of the shape parameter γ was taken to be 1.0 and 1.2. Each failure time was truncated over the same interval B=[L,R]. The left endpoint L was taken to be 0.0 and the right endpoint R was calculated at the 95th and 90th percentiles of the true marginal distribution. In summary, we generated data under 8 scenarios, each of which was simulated M=100 times.

The value of ψ was expressed in terms of φ and γ. It was obtained from the marginal survival function as follows:

0.5= S * ( η )= [1+ ( ψη ) γ ϕ ] −ϕ ,

which gave

ψ= [ ( 0.5 −1/ϕ −1)ϕ η γ ] 1/γ .

The truncation percentiles were calculated in a similar manner. They depended on both parameters φ and γ:

R= [ ( (1−p) −1/ϕ −1)ϕ ψ γ ] 1/γ = [ ( (1−p) −1/ϕ −1) ( 0.5 −1/ϕ −1) ] 1/γ *η,

where P(T≤R)=p for p=0.95 and 0.90. Upon generating uniform random variables U~ U(0,1) and V~ U(0,1), the bivariate failure times were derived from

S(t|α)= e −α (ψt) γ

as

t 1 = [ − log(1−u) α ψ γ ] 1/γ

and

t 2 = [ − log(1−v) α ψ γ ] 1/γ .

To accommodate for bivariate interval truncation, (t₁,t₂) was retained only if it belonged to rectangle B₁XB₂, that is if both times satisfied their truncation requirements. This procedure was continued until n=1000 pairs were generated from the truncation region.

The generated data of paired failure times were fit assuming a gamma distribution for the random effects and a cumulative baseline hazard function of the following forms: (i) 1-piece or exponential model (true model when γ=1) (ii) 2-piece model, and a (iii) 4-piece model.

The cut points were determined based on the percentiles of the true marginal distribution. That is,

q= [ ( (1−p) −1/ϕ −1)ϕ ψ γ ] 1/γ = [ ( (1−p) −1/ϕ −1) ( 0.5 −1/ϕ −1) ] 1/γ *η,

where P(T≤q)=p. For a 2-piece model, the cut point was computed for p=0.50, and for a 4-piece model, the cut points were computed for p=0.25,0.50 and 0.75. Parameter estimates under the bivariate truncation conditions were obtained using the proposed EM algorithm with ghost pairs discussed above. Note that the the crude rates of events during each of the piecewise constant intervals from our simulated data were taken to be the initial values of the piecewise constant rate parameters for the EM algorithm, and the tolerance level for convergence was specified at 1e-04.

Application of proposed methodology to data on siblings with schizophrenia
L.S Penrose (1945) discussed a survey on the cases of familial mental illness, conducted in 1944 by the Division of Psychiatric Research for the Ontario Department of Health. An investigation was made on all families in which two or more members had been admitted to a mental hospital in Ontario. The records date back to 1926, and thus a period of 18 years was covered. Of the families that were admitted to a mental hospital, only those families in which two siblings were admitted for a diagnosis of schizophrenia were included in our study. The age at diagnosis of schizophrenia was taken to be the age at first admission to a mental hospital for schizophrenia.

The primary objective was to determine whether there was an association within pairs in diagnosis times among schizophrenics. Estimating the cumulative hazard function and the marginal survival function of the event times were also of interest. Various forms for the baseline hazard function were examined, including: (i) Weibull with shape parameter γ and scale parameter ψ^-1 i.e.,

Λ 0 * ( t )= ( ψt ) γ

and (ii) piecewise constant models with two, three, and four pieces i.e.,

Λ 0 * ( t )= ∑ r=1 R ρ r w r ( t ),

where ρ_r is the rth unknown rate parameter corresponding to the r^th piecewise constant interval, and wr(t) is the length of the intersection of interval (0,t) and the rth piecewise constant interval. When aiming for a flexible weakly parametric model using piecewise constant rates, the selection of cut points should be considered carefully. The choice of cut points a₁,...,a_R-1) are typically determined either i) a priori (for example, if there is prior evidence indicating the event rate may begin to change after a specific time) or ii) based on the percentiles from the estimated marginal survivor function of the data. It may also be useful to begin with a small number of cut points (say 2) and then examine the effect on parameter estimates as the number of cut points increases. In our data, the breakpoints for the piecewise constant models were selected based on the estimated marginal survivor function under the Weibull fit. The two-piece model had a breakpoint at t=20; the three-piece model had breakpoints at t=20 and 40; and the four-piece model had breakpoints at t=16,28, and 48. Upon constructing the complete data log-likelihood, the parameters of the model were estimated using the EM algorithm with ghost pairs under the conditional formulation, as proposed above.

Results

Results from simulations
Results from the simulation study under the true value γ=1 and γ=1.2 are provided in Tables 1 and 2, respectively. The mean and simulation-based standard error are given for the estimates of the model parameters. Interest lied in examining how the average bias and standard error for φ and ρ vary with respect to (i) the number of pieces (ii) the amount of truncation, and (iii) the true value of the variance of the random effect distribution. Recall that as the variance φ^-1 approaches zero, there was no association among failure times within a pair.

From Table 1, when increasing the number of pieces in the model of the baseline hazard function, there was no notable effect on the bias of the estimates, however the standard error increased. The same occurred as the degree of truncation increased. Furthermore, a larger value for the variance of the random effect distribution resulted in an increase in both the bias and the standard error of the estimate for φ. Table 2 showed a strong decreasing trend in the bias of the estimate for φ as the number of pieces increased. Furthermore, a rise in the degree of truncation resulted in a notable increase in the standard error of the estimate for φ. Note that the magnitude of the bias of the estimates for all parameters in Tables 1 and 2 were in fact very small, irrespective of any change in conditions, and the standardized bias also indicated that the bias was negligible throughout.

Results from application to data on siblings with schizophrenia
The data was comprised of 173 pairs of siblings, in which each member of a pair was hospitalized with schizophrenia between the years 1926 and 1944. The age at first hospitalization was considered as the age of diagnosis (the event time). This observation was left and right truncated between the age at 1926 and the age at 1944, respectively. The ages at first hospitalization for schizophrenia for the first-born and secondborn siblings are represented by T₁ and T₂, respectively. The ages at 1926 for the first-born and second-born siblings are denoted by L₁ and L₂, respectively. And the ages at 1944 for the first-born and second-born siblings are represented by R₁ and R₂, respectively. Note that in this particular case, the length of the truncation interval is the same across individuals, and thus the area of the truncation square is the same across pairs. Data for the first 20 pairs of siblings can be found in Table 3. The data is further illustrated in Figure 2. The x-axis provides the age of the first-born sibling and the y-axis represents the age of the second-born sibling. The truncation intervals for individuals in a pair create a truncation square, and their event times are given by the corresponding point, which is only observed because it falls within their truncation square.

Table 4 provides the results from our analysis. The estimates and their standard errors are given for the parameters of the Weibull, 2-piece, 3-piece, and 4-piece models, along with the corresponding observed log-likelihood at the maximum. Once these estimates were obtained, the 25th, 50th, and 75th percentiles under each model were computed by solving

[1+ϕ Λ 0 * ( t q )] −1/ϕ =1−q

for t_q, where q=0.25,0.50,0.75. Furthermore, the natural log of the odds ratio and its standard error were computed at t=21.0 under each model. The focus on the age of 21 was determined by our clinical collaborators, who define this age to be an indicator of early versus late schizophrenia onset; it is also believed that genetic factors are strongly associated with onset of schizophrenia prior to the age of 21 [14]. Note that the natural log of the odds ratio at time t was calculated as:

log( OR( t ) )=log[ P( T 2 >t| T 1 >t)/P( T 2 ≤t| T 1 >t) P( T 2 >t| T 1 ≤t)/P( T 2 ≤t| T 1 ≤t) ],

which was simply a function of marginal and joint survival functions. The estimate of the variance of the log(OR(t)) estimator was computed using a sandwich estimator. The results of the odds ratio under all models led to similar conclusions. Specifically, under a Weibull baseline hazard model, the 95% confidence interval for the log(OR(t=21)) was (0.080, 3.780), indicating a strong presence of association among diagnosis times for schizophrenia within pairs of siblings.

The results in Table 4 were based on the joint modeling of paired failure times. For comparison purposes, however, similar measures may be obtained based on independence assumptions among event times within a pair. With a naive matched-pairs analysis under a binary response (age of onset ≤21 or age of onset >21), the estimate of the log(OR(t=21)) is 0.788 with standard error 0.381, and the 95% confidence interval for the log(OR(t=21)) was (0.035, 1.525). Furthermore, McNemar's test of marginal association for such matched-pair data gave a chi-square test statistic of 3.781 with a resulting p-value of 0.052. These naive results agreed with our more formal analysis, indicating a notable presence of association among diagnosis times for schizophrenia within pairs of siblings.

An illustration of the performance of the piecewise constant models can be seen from the plot of the estimated marginal survivor functions given in Figure 3. The 95% confidence bands are also provided under the Weibull fit at years 20, 40, and 60. The results from the piecewise models agree with that from the Weibull model. However, during the first two decades in which very few events are observed to occur, the four-piece model provides a more suitable fit than the two-piece or three-piece model.

Although the piecewise model fit could be carried out using direct maximization, the EM algorithm was used to avoid divergence problems often arising when there is a large number of pieces. In the case of a 2-piece model, the computing time using a standard maximization function was faster than implementing an EM approach. However, if one wanted to fit an 8-piece or 10-piece model, the EM algorithm was attractive since it did not require obtaining the derivative functions, and thus avoided possible problems with convergence.

Discussion

This paper was concerned with inference for bivariate failure time data that were interval-truncated. A conditional (random effect) formulation was discussed for modeling the marginal and joint distribution of the failure times, where particular interest lied in measuring the association within pairs. Under the conditional model, the marginal distributions were determined by both the variance parameter of the random effect distribution and the parameters of the conditional hazard function. To minimize the need for strong distributional assumptions and facilitate a likelihood analysis based on parametric models, the baseline hazard functions were assumed to be of a piecewise constant form.

The literature offers limited approaches in estimation when both failure times in a pair are observed subject to interval truncation. This paper developed an EM approach for estimation based on a conditional formulation. The idea of "ghosts" was introduced in the bivariate setting so that the complete data log-likelihood could be easily maximized. The formulation was convenient to implement since the E-step was relatively straightforward, and it also provided closed form expressions for the parameter estimates of the piecewise constant conditional hazard function in the M-step. Results from the simulation studies indicated that all model parameters could be estimated with negligible bias in finite samples. The data simulated for Table 1 was based on an exponential baseline hazard function (as the shape parameter γ=1), which is why more pieces were not needed to improve model fit. On the other hand, as the baseline hazard function was not constant for the data simulated in Table 2, we see that increasing the number of pieces provided a better approximation to the true marginal survivor function. In the motivating schizophrenia data set, our proposed model and method of estimation revealed a strong presence of dependence among schizophrenia diagnosis times within pairs of siblings. In addition, since the rate of onset for schizophrenia was quite low in the early years, the four-piece model was more sensitive to this shift and was able to detect this shift, thus providing different results than the two-piece and three-piece models. Moreover, as the estimate of the shape parameter under the Weibull distribution was much greater than 1.0, it was reasonable that increasing the number of pieces provided a better fit, as seen by the plots of the estimated marginal survivor functions in Figure 3.

Although we considered the scenario where there was no censoring, our proposed methods could be naturally extended to handle the case where the event times were interval-censored within their respective truncation intervals. We also considered the scenario in which there were paired failure times; our conditional approach could be easily extended to handle a multivariate setting with more than 2 individuals within a group.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

Authors' contributions	RS	RJC
Research concept and design	√	√
Collection and/or assembly of data	√	√
Data analysis and interpretation	√	√
Writing the article	√	√
Critical revision of the article	√	√
Final approval of article	√	√
Statistical analysis	√	√

Acknowledgement and funding

This work was supported by the Institute for Clinical Evaluative Sciences. Richard Cook is a Canada Research Chair in Statistical Methods for Health Research. We would also like to thank Dr. Janice Husted from the University of Waterloo, Canada, and Dr. Anne Bassett from the University of Toronto, Canada, for their valuable input and guidance regarding the data.

Publication history

EICs: Jimmy Efird, East Carolina University, USA.
Max K. Bulsara, University of Notre Dame, Australia.
Received: 17-Nov-2015 Final Revised: 31-Dec-2015
Accepted: 13-Jan-2016 Published: 29-Jan-2016