A quantile cure model with partially functional covariate effects

2.1. Model assumptions

Let η * denote the cure status of an individual. An individual with η * = 1 will eventually experience the event of interest with event time T * if there is no censoring, while an individual with η * = 0 is cured with survival time equal to ∞ , indicating that the individual will never experience the event of interest. Hence, the failure time of an individual can be written as T ~ * = T * η * + ∞ ( 1 − η * ) . In practice, there may exist a censoring time C * and a finite administrative censoring time C ~ such that T ~ * can be righted censored by C = min { C * , C ~ } . Let Δ = I ( T ~ * ≤ C ) and T ~ = min { T ~ * , C } . An individual with Δ = 1 is known to be uncured and η * = 1 . However, the cure status of an individual with Δ = 0 is unknown. Let X be a set of covariates that may have effects on T * , η * , and C , and assume that C is independent of η * and T * given X . Let observed data { ( T ~ i , Δ i , X i ) , i = 1 , … , n } be n independent and identically distributed ( i . i . d . ) copies of ( T ~ , Δ , X ) .

Assume that the uncured probability is associated with a p -vector of covariates X

P ( η * = 1 | X ) = π ( X ′ γ )

(1)

For the distribution of T * , let τ be a quantile level of interest and Q T * ( τ | X ) = inf { t : P ( T * ≤ t | X ) ≥ τ } denote the τ th quantile of T * given X . Assume that at the specific τ , Q T * ( τ | X ) is modeled by

Q T * ( τ | X ) = exp ( X ′ β ( τ ) )

(2)

Let γ 0 be the true value of γ . Since Δ i = 1 implies η i * = 1 and C i > T ~ i * , using the conditional independence between C and ( η * , T ~ * ) given X and following Chen et al., ¹² we have the following equation E { Δ i / G ( T ~ i | X i ) | X i } = E [ I ( η i * = 1 ) E { I ( C i > T ~ i * ) / G ( T ~ i * | X i ) | T ~ i * , X i } | X i ] = π ( X i ′ γ 0 ) , where G ( t | X ) is the survival function of censoring time C that may depend on X . If G ^ ( ⋅ ) is a consistent estimator of G ( ⋅ ) , we have the following estimating equation for γ

1 n ∑ i = 1 n X i ( Δ i G ^ ( T ~ i | X i ) − π ( X i ′ γ ) ) = 0

(3)

To estimate β ( τ ) , consider E ( Δ i I { T ~ i ≤ t } / G ( T ~ i | X i ) | X i ) = E ( E { I { T ~ i ≤ t , η i * = 1 , C i > T ~ i * } / G ( T ~ i | X i ) | T ~ i * , η i * , X i } | X i ) = π ( X i ′ γ 0 ) P ( T i * ≤ t | η i * = 1 , X i ) . Hence, E ( Δ i I { T ~ i ≤ Q T * ( τ | X i ) } / G ( T ~ i | X i ) ) = π ( X i ′ γ 0 ) τ implies that

1 n ∑ i = 1 n X i ( Δ i I { log T ~ i ≤ X i ′ β ( τ ) } G ^ ( T ~ i | X i ) − π ( X i ′ γ ~ ) τ ) = 0

(4)

S n ( β ( τ ) ) = 1 n ∑ i = 1 n X i Δ i G ^ ( T ~ i | X i ) ( I { log T ~ i ≤ X i ′ β ( τ ) } − τ ) = 0

(5)

h n ( β ( τ ) ) = ∑ i = 1 n Δ i G ^ ( T ~ i | X i ) | log T ~ i − β ( τ ) ′ X i | + | M + β ( τ ) ′ ∑ i = 1 n X i Δ i G ^ ( T ~ i | X i ) | + | M − β ( τ ) ′ ∑ i = 1 n 2 X i Δ i τ G ^ ( T ~ i | X i ) |

(6)

We suggest estimating G ( ⋅ ) using the Kaplan–Meier estimator if it does not depend on X . If G ( ⋅ ) depends on X , we may model C with some survival models, for example, Cox’s PH model, the AFT model, or a nonparametric estimation method ¹⁴ to take the effect of X into account in the estimate G ^ ( t | X ) . When modeling survival data with a cure model, the support of C is required to cover the support of T * to meet the sufficient follow-up requirement so that the cure probability can be properly estimated from the data.^15,16 Hence, consistency and the weak convergence property of G ^ ( ⋅ ) hold on the support of T * ¹⁷ because the estimate of G ^ ( ⋅ ) will be a proper survival function, and it is well defined for a time within the support of T * .

The asymptotic properties of γ ~ have been investigated in Chen et al. ¹² and its asymptotic variance can be estimated by applying the plug-in method. We established the consistency and the asymptotic normality of β ~ ( τ ) . The details of the asymptotic properties and their proof are relegated to Section A in the Supplemental materials. Following Peng and Fine, ¹³ the plug-in estimator for the asymptotic variances can be derived. However, according to our simulation experience, the variance estimate for β ~ ( τ ) using the plug-in method is not reliable and can sometimes cause a smaller empirical coverage rate than 95% for the Wald-type 95% confidence interval based on the variance estimate. This is also reported and discussed by Choi et al. ¹⁸ Consequently, we suggest applying the weighted bootstrap approach ¹⁹ for the variance estimation. The details of this approach are provided in Section B in the Supplemental materials.

It is noteworthy that our approach is fundamentally different from that of Wu and Yin. ³ Wu and Yin ³ developed a two-stage estimation method, where γ is estimated in Stage 1 and then β ( τ ) is estimated in Stage 2 using the multiple imputation and the results from Stage 1. The local Kaplan–Meier estimators for the latency survival function are repeatedly used in both stages. However, the locally nonparametric estimator suffers from the curse of dimensionality. Moreover, the estimation in the latency part of their method requires consistent estimation results in Stage 1. Thus, a misspecification in the incidence part can have a nonignorable impact on the quantile regression estimation in Stage 2. In contrast, in our proposed approach, the estimation in one part does not depend on the other part. This can be seen from the fact that the estimating equation (3) does not depend on the latency part and the estimating equation (5) does not depend on the incidence part, which makes the estimation doubly robust to the misspecification of either part of the model. This advantage holds for the quantile cure model only because P ( log ( T i * ) ≤ X i ′ β ( τ ) | X i ) is equal to the common τ for all i . A disadvantage of the proposed method is that we need to estimate the censoring distribution. If the censoring mechanism depends on a high-dimensional set of covariates in a way that cannot be sufficiently modeled by the Cox PH or the AFT models, and a nonparametric estimation has to be considered, then the proposed method may also suffer from the curse of dimensionality. However, this concern is at a much lesser degree than that in the method of Wu and Yin ³ because the censoring time C is usually considered separately from the event time of interest, and it is not uncommon to assume a simpler model for the association between C and a limited number of covariates. On the other hand, Wu and Yin ³ suggested modeling F T * ( t | X ) nonparametrically. The function F T * ( t | X ) is for the event of interest, and if X is high-dimensional, the concern of the curse of dimensionality will inevitably arise. Moreover, this approach basically considered two models for T * , one is the conditional quantile model and the other is the nonparametric estimation method. There is no guarantee that the two models produce coherent estimates, which can in turn affect the final estimates adversely.

3.1. Fully functional covariate effects

If the linearity assumption (2) holds for all covariates over a range of quantiles in [ τ L , τ U ] , where 0 < τ L < τ U < 1 , β ~ ( τ ) can be estimated from (5). We consider β ( τ ) as a right-continuous step function of τ at artificially selected jump points I K ( n ) = { τ L = τ 0 < τ 1 < ⋯ < τ K ( n ) = τ U } and estimate β ( τ k ) , k = 0 , 1 , … , K ( n ) , in a pointwise manner using the estimating equation (5) and computation algorithm (6). One advantage of our approach is that the estimation of β ( τ k ) at τ k does not require to be carried out sequentially from lower quantiles. We set the mesh size ‖ I K ( n ) ‖ = max k = 1 , … , K ( n ) | τ k − τ k − 1 | as order of o ( n − 1 / 2 ) . Under regularity conditions, we can prove that the estimator β ~ ( τ ) on [ τ L , τ U ] is uniformly consistent on [ τ L , τ U ] and converges to a Gaussian process at a root- n rate. The details of the asymptotic properties and their proof are provided in Section C in the Supplemental materials.

3.2. Partially functional covariate effects

The fully functional covariate effects in Section 3.1 assume that the effect β ( τ ) in (2) holds for all covariates in X on τ ∈ [ τ L , τ U ] . However, in practice, the quantile effects of some covariates may be invariant over the range [ τ L , τ U ] . To make the estimation more efficient, it is important to take this information into account in the estimation procedure. Let X = ( X 1 ′ , X 2 ′ ) ′ , where X 1 is a p 1 -dimensional covariate vector with quantile-varying effect β ( τ ) and X 2 is a p 2 -dimensional covariate vector with quantile-invariant effects α , where p 1 + p 2 = p . The quantile regression model in the latency part assumes the partially functional covariate effects

Q T * ( τ | X ) = exp ( X 1 ′ β ( τ ) + X 2 ′ α )

(7)

Similar to Section 2, we assume that the observed data { ( T ~ i , Δ i , X 1 i , X 2 i ) , i = 1 , … , n } are i . i . d . copies of ( T ~ , Δ , X 1 , X 2 ) . The estimating equation of γ in the incidence part is still (3). The estimating equation for β ( τ ) for any τ ∈ [ τ L , τ U ] can be obtained from (5) as

U 1 n ( β ( τ ) ; α ) = 1 n ∑ i = 1 n X 1 i Δ i G ^ ( T ~ i | X i ) ( I { log T ~ i ≤ X 1 i ′ β ( τ ) + X 2 i ′ α } − τ ) = 0

(8)

U 2 n ( α ; β ) = 1 n ∑ i = 1 n ∫ τ L τ U X 2 i Δ i G ^ ( T ~ i | X i ) ( I { log T ~ i ≤ X 1 i ′ β ( u ) + X 2 i ′ α } − u ) d u = 0

(9)

h 1 n ( β ( τ ) ; α ) = ∑ i = 1 n Δ i G ^ ( T ~ i | X i ) | log T ~ i − X 2 i ′ α − β ( τ ) ′ X 1 i | + | M + β ( τ ) ′ ∑ i = 1 n X 1 i Δ i G ^ ( T ~ i | X i ) | + | M − β ( τ ) ′ ∑ i = 1 n 2 X 1 i Δ i τ G ^ ( T ~ i | X i ) |

U 2 n ( α ; β ) = 1 n ∑ i = 1 n ∑ k = 0 K ( n ) − 1 X i Δ i G ^ ( T ~ i | X i ) I { log T ~ i ≤ X 1 i ′ β ( τ k ) + X 2 i ′ α } ( τ k + 1 − τ k ) − 1 2 n ∑ i = 1 n X 2 i Δ i G ^ ( T ~ i | X i ) ( τ U 2 − τ L 2 )

(10)

h 2 n ( α ; β ) = ∑ i = 1 n ∑ k = 0 K ( n ) − 1 Δ i ( τ k + 1 − τ k ) G ^ ( T ~ i | X i ) | log T ~ i − X 1 i ′ β ( τ k ) − X 2 i ′ α | + | M + α ′ ∑ i = 1 n X 2 i Δ i G ^ ( T ~ i | X i ) ( τ U − τ L ) | + | M − α ′ ∑ i = 1 n X 2 i Δ i G ^ ( T ~ i | X i ) ( τ U 2 − τ L 2 ) |

The choice of τ U in the interval [ τ L , τ U ] is very crucial, since it may cause the nonidentifiability issue and affect the estimation of α and β ( τ ) if it is too wide. Wu and Yin ⁵ proposed an adaptive approach to select the largest τ U such that τ U < 1 and the identifiability of the parameters for τ ∈ [ τ L , τ U ] are ensured, and this approach can be employed here.

We considered methods to check the fits of the proposed models for all τ ∈ ( 0 , 1 ) . There are a few methods proposed in the literature to check the fit of a cure model. We adapted the Cox–Snell residuals proposed for the mixture cure model by Peng and Taylor. ²⁰ Let S u ( t | X ) = P ( T * > t | X ) . The residual for subject i based on an estimated cure model is defined as r i = − log { 1 − π ( X i ′ γ ^ ) + π ( X i ′ γ ^ ) S ^ u ( T ~ i | X i ) } . It is proved that if the estimated cure model fits the data well for all τ ∈ ( 0 , 1 ) , the residual can be viewed as a random number from a mixed-type distribution with a unit exponential distribution as the continuous component on ( 0 , − log { 1 − π ( X i ′ γ ^ ) } ) and a probability mass 1 − π ( X i ′ γ ^ ) at − log { 1 − π ( X i ′ γ ^ ) } . The distribution of the censored residuals ( r i , Δ i ) , i = 1 , … , n , can be compared with the unit exponential distribution to determine the fit of the estimated cure model and this is usually done by plotting the estimated cumulative hazard function of the censored residuals and comparing it with a 45 ∘ line. This approach, however, cannot be applied directly to the proposed model as the proposed method does not estimate S u ( t | X ) . We proposed the following approach to approximate the residuals: Consider R grid points τ k = k / ( R + 1 ) , k = 1 , … , R between ( 0 , 1 ) . We use the proposed method to estimate γ ^ and all β ^ ( τ k ) , k = 1 , … , R , and then estimate S u ( T ~ i | X i ) by τ k if exp { X i ′ β ^ ( τ k ) } ≤ T ~ i < exp { X i ′ β ^ ( τ k + 1 ) } . The corresponding Cox–Snell residual for the subject i is obtained as r i = − log { 1 − π ( X i ′ γ ^ ) + π ( X i ′ γ ^ ) τ k } . If the model contains an estimated quantile-invariant effect α ^ , then τ k is determined by exp { X 1 i ′ β ^ ( τ k ) + X 2 i ′ α ^ } ≤ T ~ i < exp { X 1 i ′ β ^ ( τ k + 1 ) + X 2 i ′ α ^ } . We suggest that R = 99 .

An interesting question in practice is how to choose between X 1 and X 2 for a set of covariates. We suggest that a model with all covariates in X 1 should be fit first, and then the plot of each estimate in β ^ ( τ ) against τ is visually examined and followed by a hypothesis test to examine the significance of the quantile-varying effect in an interval identified in the visual inspection. If the plot of an estimate does not show a substantial change in the interval and the test shows no significant quantile-varying effect, the corresponding covariate is moved from X 1 to X 2 . This approach provides investigators with a method to select a parsimonious quantile model in the latency part and efficiently summarizes the effect of a covariate via one parameter instead of a quantile-varying effect when it is unnecessary.

To test whether or not a covariate has a quantile-varying effect in a prespecified interval, we follow the approach proposed by Peng and Huang ⁴ to conduct hypothesis testing. Suppose that we are interested in testing the significance of the quantile-varying effect β ^ 1 j ( τ ) in β ^ ( τ ) for the j th covariate in X 1 over the interval [ τ L , τ U ] . The hypotheses are

H 0 : β 1 j ( τ ) = c j for all τ ∈ [ τ L , τ U ] versus H a : β 1 j ( τ ) ≠ c j for some τ ∈ [ τ L , τ U ]

(11)

T = n 1 / 2 ∫ τ L τ U ( β ^ 1 j ( u ) − c ^ j ) Ξ ( u ) d u

(12)

4.1. Covariate effects on quantiles

We consider X = ( 1 , X 1 , X 2 ) ′ , where X 1 and X 2 are respectively generated independently of the Bernoulli distribution with success probability 0.5 and uniform distribution on ( 0 , 1 ) . Three link functions for π ( X ′ γ ) in (1) are considered in data generation:

(i)

the logit link π ( X ′ γ ) = exp ( X ′ γ ) / [ 1 + exp ( X ′ γ ) ] ,

(ii)

the probit link π ( X ′ γ ) = Φ ( X ′ γ ) ,

(iii)

the complementary log–log link π ( X ′ γ ) = 1 − exp [ − exp ( X ′ γ ) ] ,

with ( γ 0 , γ 1 , γ 2 ) = ( 1 , − 1 , 1 ) . These three link functions respectively produce 72%, 80%, and 88% as the average of uncure rates. For the latency part, assume log T * = b 0 + b 1 X 1 + α X 2 + ( 1 + X 1 ) ϵ with ( b 0 , b 1 , α ) = ( 1 , − 1 , 0.5 ) . It implies log Q T * ( τ | X ) = β 0 ( τ ) + β 1 ( τ ) X 1 + β 2 ( τ ) X 2 . Two types of distributions for ϵ are considered: the standard normal distribution N ( 0 , 1 ) with P ( ϵ ≤ u ) = Φ ( u ) and the extreme value distribution with P ( ϵ ≤ u ) = F ( u ) = 1 − exp ( − exp ( u ) ) . When ϵ follows the standard normal distribution, β 0 ( τ ) = Φ − 1 ( τ ) + b 0 , β 1 ( τ ) = Φ − 1 ( τ ) + b 1 , β 2 ( τ ) = α . When ϵ follows the extreme value distribution, β 0 ( τ ) = F − 1 ( τ ) + b 0 , β 1 ( τ ) = F − 1 ( τ ) + b 1 , β 2 ( τ ) = α . The censoring time C follows a mixture distribution with 80% from the exponential distribution with a rate of 0.025 and right-truncated at time 80 and 20% from a degenerated distribution at 80 if ϵ follows the standard normal distribution or a mixture distribution with 80% from the exponential distribution with a rate of 0.025 and right-truncated at time 60 and 20% from a degenerated distribution at 60 if ϵ follows the extreme value distribution. The resulting censoring rate among the uncured subjects is about 12% if ϵ follows the standard normal distribution and 8% if ϵ follows the extreme value distribution. The overall average censoring rates are respectively 36%, 29%, and 21% for the logit, probit, and complementary log–log links in the incidence parts if ϵ follows the standard normal distribution, and 34%, 25%, and 18% if ϵ follows the extreme value distribution. Additional simulation studies with higher censoring rates are provided in Section E in the Supplemental materials.

We generate 500 Monte Carlo datasets with a sample size of n = 300 for each data configuration. For each dataset, the estimation method in Section 2 is used to estimate γ and β ( τ ) with the logit link used in the incidence part (an R program was developed for the proposed estimation methods and tests, and it is available in the Supplemental materials). Consequently, our approach correctly specifies the incidence part in case (i), but misspecifies it in cases (ii) and (iii). The standard error estimates SEEs for β ~ ( τ ) are based on the weighted bootstrap method with B = 200 , while those for γ ~ are based on the plug-in type standard error estimator of Chen et al. ¹² The 95% confidence intervals for β ( τ ) are based on the 2.5th and 97.5th percentiles from the weighted bootstrap estimates. The simulation results for ϵ following the standard normal distribution and extreme value distribution are respectively summarized in Tables 1 and 2. Both tables include the average estimates (EST), the sample standard deviations (SD), the averages of the SEEs, and the coverage rates of 95% confidence intervals (CP), and the mean square errors (MSEs) from the proposed method. We also fit the data using the method proposed by Wu and Yin ³ (referred to as the WY method) for comparison, but only report ESTs, SDs, and MSEs to save time. The biquadratic kernel function with 10 imputation times is used in the WY method. The relative efficiency is measured by the ratios of MSEs of the proposed method to those of the WY method.

OPEN IN VIEWER

Table 1.

Simulation results of the parameter estimates from the proposed method and the WY method for the data generated under three link functions in the incidence part and the standard normal distribution for ϵ in the latency part.

		Proposed method					WY method
	True value	EST	SD	SEE	CP	MSE	EST	SD	MSE	Ratio
Logit
γ 0	1.000	1.026	0.405	0.374	0.946	0.164	1.119	0.468	0.232	0.707
γ 1	−1.000	−1.064	0.376	0.358	0.956	0.145	−1.138	0.424	0.199	0.730
γ 2	1.000	0.987	0.623	0.605	0.964	0.388	1.062	0.719	0.520	0.746
β 0 ( 0.25 )	0.326	0.334	0.257	0.268	0.954	0.066	0.343	0.300	0.090	0.731
β 1 ( 0.25 )	−1.674	−1.715	0.293	0.334	0.978	0.087	−1.700	0.366	0.134	0.650
β 2 ( 0.25 )	0.500	0.497	0.418	0.453	0.968	0.174	0.499	0.500	0.249	0.699
β 0 ( 0.5 )	1.000	1.008	0.236	0.255	0.958	0.056	1.030	0.273	0.075	0.739
β 1 ( 0.5 )	−1.000	−1.073	0.307	0.311	0.950	0.099	−1.052	0.347	0.123	0.809
β 2 ( 0.5 )	0.500	0.494	0.401	0.431	0.958	0.161	0.479	0.462	0.213	0.754
β 0 ( 0.75 )	1.674	1.694	0.275	0.303	0.954	0.076	1.732	0.307	0.098	0.774
β 1 ( 0.75 )	−0.326	−0.443	0.350	0.380	0.930	0.136	−0.435	0.379	0.155	0.876
β 2 ( 0.75 )	0.500	0.468	0.469	0.513	0.962	0.220	0.445	0.530	0.283	0.777
Probit
γ 0	1.000	1.860	0.712	0.875	0.862	1.246	2.094	0.953	2.103	0.592
γ 1	−1.000	−1.955	0.738	0.908	0.858	1.456	−2.160	0.940	2.228	0.653
γ 2	1.000	1.741	0.910	0.862	0.938	1.375	1.823	0.984	1.642	0.838
β 0 ( 0.25 )	0.326	0.329	0.243	0.250	0.958	0.059	0.338	0.282	0.079	0.746
β 1 ( 0.25 )	−1.674	−1.719	0.282	0.317	0.972	0.081	−1.709	0.361	0.131	0.622
β 2 ( 0.25 )	0.500	0.505	0.392	0.425	0.974	0.153	0.490	0.473	0.223	0.688
β 0 ( 0.5 )	1.000	1.010	0.223	0.237	0.944	0.050	1.020	0.272	0.074	0.671
β 1 ( 0.5 )	−1.000	−1.068	0.287	0.295	0.946	0.087	−1.042	0.353	0.126	0.688
β 2 ( 0.5 )	0.500	0.494	0.374	0.405	0.950	0.139	0.474	0.452	0.204	0.683
β 0 ( 0.75 )	1.674	1.697	0.266	0.281	0.956	0.071	1.727	0.306	0.096	0.742
β 1 ( 0.75 )	−0.326	−0.436	0.333	0.359	0.928	0.123	−0.416	0.402	0.169	0.728
β 2 ( 0.75 )	0.500	0.460	0.448	0.480	0.958	0.202	0.427	0.505	0.260	0.776
Complementary log–log
γ 0	1.000	2.954	1.091	3.164	0.778	5.006	3.428	1.877	9.408	0.532
γ 1	−1.000	−2.516	1.138	3.251	0.996	3.590	−2.974	1.851	7.316	0.491
γ 2	1.000	2.060	1.611	1.563	0.990	3.713	2.329	2.149	6.373	0.583
β 0 ( 0.25 )	0.326	0.328	0.235	0.239	0.950	0.055	0.364	0.276	0.078	0.707
β 1 ( 0.25 )	−1.674	−1.708	0.272	0.294	0.960	0.075	−1.673	0.331	0.109	0.690
β 2 ( 0.25 )	0.500	0.492	0.376	0.409	0.966	0.141	0.431	0.448	0.205	0.689
β 0 ( 0.5 )	1.000	0.996	0.215	0.227	0.954	0.046	1.021	0.241	0.058	0.793
β 1 ( 0.5 )	−1.000	−1.053	0.264	0.276	0.944	0.072	−1.039	0.309	0.097	0.746
β 2 ( 0.5 )	0.500	0.495	0.365	0.391	0.958	0.133	0.454	0.414	0.173	0.771
β 0 ( 0.75 )	1.674	1.677	0.245	0.263	0.956	0.060	1.714	0.284	0.082	0.731
β 1 ( 0.75 )	−0.326	−0.417	0.313	0.334	0.950	0.106	−0.392	0.346	0.124	0.854
β 2 ( 0.75 )	0.500	0.467	0.425	0.452	0.960	0.181	0.409	0.491	0.248	0.731

WY: Wu and Yin; EST: estimate; SEE: standard error estimate; CP: coverage rates of 95% confidence interval; MSE: mean square error.

OPEN IN VIEWER

Table 2.

Simulation results of the parameter estimates from the proposed method and the WY method for the data generated under three link functions in the incidence part and the extreme value distribution for ϵ in the latency part.

		Proposed method					WY method
	True value	EST	SD	SEE	CP	MSE	EST	SD	MSE	Ratio
Logit
γ 0	1.000	1.025	0.333	0.346	0.960	0.111	1.094	0.422	0.187	0.595
γ 1	−1.000	−1.031	0.303	0.326	0.964	0.093	−1.074	0.376	0.147	0.633
γ 2	1.000	0.991	0.543	0.554	0.976	0.294	1.015	0.632	0.398	0.738
β 0 ( 0.25 )	−0.246	−0.256	0.374	0.387	0.952	0.139	−0.267	0.455	0.558	0.250
β 1 ( 0.25 )	−2.246	−2.218	0.468	0.484	0.962	0.220	−2.225	0.552	0.305	0.720
β 2 ( 0.25 )	0.500	0.487	0.631	0.650	0.948	0.397	0.512	0.768	0.589	0.675
β 0 ( 0.5 )	0.633	0.636	0.273	0.287	0.966	0.074	0.640	0.336	0.242	0.307
β 1 ( 0.5 )	−1.367	−1.363	0.321	0.353	0.960	0.103	−1.379	0.400	0.160	0.646
β 2 ( 0.5 )	0.500	0.470	0.474	0.481	0.958	0.225	0.473	0.556	0.310	0.727
β 0 ( 0.75 )	1.327	1.334	0.246	0.265	0.962	0.060	1.343	0.307	0.204	0.296
β 1 ( 0.75 )	−0.673	−0.676	0.304	0.333	0.962	0.092	−0.682	0.370	0.136	0.677
β 2 ( 0.75 )	0.500	0.475	0.436	0.452	0.958	0.190	0.478	0.526	0.276	0.689
Probit
γ 0	1.000	1.835	0.534	0.562	0.770	0.982	1.996	0.741	1.541	0.638
γ 1	−1.000	−1.907	0.551	0.578	0.766	1.125	−2.072	0.716	1.661	0.677
γ 2	1.000	1.796	0.770	0.758	0.866	1.225	1.874	0.888	1.552	0.789
β 0 ( 0.25 )	−0.246	−0.262	0.364	0.364	0.944	0.132	−0.260	0.428	0.525	0.252
β 1 ( 0.25 )	−2.246	−2.224	0.442	0.458	0.954	0.195	−2.215	0.527	0.278	0.701
β 2 ( 0.25 )	0.500	0.510	0.621	0.612	0.940	0.385	0.497	0.720	0.517	0.745
β 0 ( 0.5 )	0.633	0.629	0.259	0.267	0.940	0.067	0.633	0.303	0.227	0.296
β 1 ( 0.5 )	−1.367	−1.361	0.315	0.333	0.942	0.099	−1.353	0.383	0.147	0.677
β 2 ( 0.5 )	0.500	0.489	0.445	0.451	0.960	0.198	0.478	0.496	0.246	0.807
β 0 ( 0.75 )	1.327	1.329	0.239	0.244	0.962	0.057	1.343	0.279	0.187	0.304
β 1 ( 0.75 )	−0.673	−0.679	0.290	0.313	0.954	0.084	−0.676	0.362	0.131	0.643
β 2 ( 0.75 )	0.500	0.487	0.423	0.420	0.960	0.178	0.459	0.465	0.217	0.821
Complementary log–log
γ 0	1.000	2.981	0.913	1.834	0.688	4.756	3.424	1.475	8.045	0.591
γ 1	−1.000	−2.540	0.922	1.898	0.982	3.221	−3.003	1.461	6.139	0.525
γ 2	1.000	2.203	1.159	1.193	0.938	2.787	2.357	1.127	3.110	0.896
β 0 ( 0.25 )	−0.246	−0.276	0.343	0.349	0.948	0.118	−0.264	0.400	0.160	0.739
β 1 ( 0.25 )	−2.246	−2.217	0.403	0.429	0.960	0.163	−2.207	0.479	0.231	0.705
β 2 ( 0.25 )	0.500	0.512	0.571	0.593	0.950	0.326	0.514	0.658	0.433	0.754
β 0 ( 0.5 )	0.633	0.625	0.240	0.257	0.964	0.058	0.632	0.281	0.079	0.733
β 1 ( 0.5 )	−1.367	−1.340	0.292	0.314	0.962	0.086	−1.342	0.364	0.133	0.644
β 2 ( 0.5 )	0.500	0.479	0.415	0.438	0.960	0.172	0.471	0.467	0.218	0.789
β 0 ( 0.75 )	1.327	1.323	0.231	0.233	0.952	0.053	1.347	0.245	0.060	0.880
β 1 ( 0.75 )	−0.673	−0.663	0.270	0.288	0.958	0.073	−0.673	0.311	0.097	0.754
β 2 ( 0.75 )	0.500	0.485	0.406	0.404	0.950	0.165	0.456	0.429	0.185	0.891

WY: Wu and Yin; EST: estimate; SEE: standard error estimate; CP: coverage rates of 95% confidence interval; MSE: mean square error.

The results in both tables show that when the link function in the incidence part is correctly specified, that is, in the logit link case, the estimates of the proposed method are approximately unbiased, SEE is close to SD, and the coverage probabilities are close to the nominal level. The estimates from the proposed method have smaller SD and MSE than those from the WY method, and the MSE ratios are about 70% to 88%, indicating that the proposed method gains substantial efficiency over the WY method in this case. When the link function in the incidence part is misspecified by both methods, that is, in the probit and complementary log–log link function cases, the estimates of β ( τ ) from the proposed method are still robust and have smaller MSEs than those from the WY method. It is noteworthy that the biases in β ^ 1 ( 0.75 ) and β ^ 2 ( 0.75 ) tend to be larger in both methods for the cases of standard normally distributed ϵ than those of extreme value distributed ϵ . It may be related to a higher censoring rate in the former than that in the latter such that the former provides less information for estimating β ( τ ) when τ is large. In both tables, most of the estimates of β ( τ ) from the proposed method have smaller biases than those from the WY method, indicating that the proposed method tends to be more robust to the misspecification of the incidence part.

It is not unexpected that the estimates of γ do not perform well under the probit and the complementary log–log links because of the misspecified model for the incidence part.

We further investigated the performance of the test (12) proposed in Section 3.3 to test the significance of quantile-varying effects. In the settings considered above, β 1 ( τ ) is a quantile-varying effect while β 2 ( τ ) is quantile-invariant. Thus, we can see the power and the size of the test when applying the test to the simulated data to test β 1 ( τ ) and β 2 ( τ ) respectively. When using the test, we set Ξ ( u ) = 2 I ( u < ( τ L + τ U ) / 2 ) / ( τ U − τ L ) , where [ τ L , τ U ] = [ 0.25 , 0.75 ] . The results are presented in Table 3. They show that the power of the test is good and the size of the test is slightly conservative. The results are insensitive to the error distribution choices and the specification of the link function in the incidence part.

OPEN IN VIEWER

Table 3.

Simulation results for the size and power of the homogeneous test of the proposed quantile regression model.

	Logit		Probit		Cloglog
Error distribution	Power	Size	Power	Size	Power	Size
Standard normal	0.866	0.042	0.918	0.028	0.962	0.030
Extreme value	0.880	0.040	0.918	0.052	0.956	0.040

4.2. Partially functional covariate effects

We further investigate the performance of the estimates of the partially functional covariate effects proposed in Section 3.2 using the same simulated data in Section 4.1. Under both settings, it is known that β 2 ( τ ) is constant, and we will treat it as α in (7) and use the method in Section 3.2 to estimate β 0 ( τ ) , β 1 ( τ ) , and α on the interval [ τ L , τ U ] = [ 0.25 , 0.75 ] . We will compare the estimates with the fully functional covariate effect estimates from the method in Section 3.1. Since the estimates for γ are unchanged, we are only interested in the estimates of β 0 ( τ ) , β 1 ( τ ) and α . Table 4 shows the biases, SDs, SEEs, CPs, and MSEs of the estimates from the method in Section 3.2. The gain in efficiency by taking the homogeneous effect information of X 2 into account is demonstrated in the table by the ratios of the MSEs of β ^ 0 ( τ ) and β ^ 1 ( τ ) to those of β ~ 0 ( τ ) and β ~ 1 ( τ ) , and the ratios of the MSEs of α ^ to the MSEs of α ~ , where α ~ = ∫ 0.25 0.75 β ~ 2 ( τ ) d τ / 0.5 estimated via the fully functional covariate effects model.

The results in the table show that the performance of the proposed estimation method is satisfactory: the biases are small, the SDs are similar to the SEEs, and the coverage probabilities are close to the nominal level. All the ratios in the estimation of β 0 ( τ ) and β 1 ( τ ) are between 0.502 and 1.03, indicating that taking the constant effect into account produces efficiency gain in estimating parameters in the latency part for most of the parameters. It is interesting to see that all the MSE ratios for α are between 0.79 and 0.90. It indicates that a substantial efficiency gain is achieved by incorporating the information of quantile-invariant effects in the estimating procedure when compared to a naive method by averaging the estimate β ~ 2 ( τ ) over [ τ L , τ U ] .