Weighted censored quantile regression

Section 4. Performance analysis

• Previous
• Next

We conduct extensive simulation studies to compare the performance between our proposed EL based weighted censored quantile regression estimator and the standard censored quantile regression estimator. For our simulation, we use the models discussed in Tang and Leng (2012).

The simulation models used to generate the logarithmic event time $\left(T_r\right)$ and logarithmic censoring time $\left(C_r\right)$ for the $r^{\text{th}}$ $\left(r=1,2,\dots ,N\right)$ subject are given in Table 4.1 under four Cases (i)-(iv).

In Cases (i) and (ii), event times and censoring times are generated from the homoscedastic models and in Cases (iii) and (iv), we considered heteroscedastic models to examine the efficiency gain of our proposed method over the standard censored quantile regression. We set the parameter values as ${\theta }^{\top }=\left(0,-1,0.2\right),{\pi }^{\top }=\left(0.3,-0.1,0.1\right)$ and selected ${\gamma }^{\top }$ to maintain approximately 30% of the censoring proportion in each case. We generated explanatory variables from zero mean bivariate normal distribution with covariance,

$$\Sigma =\left[\begin{array}{cc}1& {\sigma }_{x_1,x_2}\\ {\sigma }_{x_1,x_2}& 1\end{array}\right].$$

We considered two different ways to compute the EL based probability weights. In numerical study -I, we compute $P_i’\mathrm{s}$ based on the auxiliary information related to the failure time, $T_i,$ whereas in numerical study -II, $P_i’\mathrm{s}$ are computed using the observed survival time, $Y_i=\min \left(T_i,C_i\right).$ In numerical study -II, we employ the synthetic variable approach (Koul et al., 1981; Qin and Jing, 2001; Li and Wang, 2003) to compute the EL based data driven probability weights.

4.1  Numerical study -I

To compute $P_i’\mathrm{s},$ first we need to have a known population parameter, $\theta ,$ or its estimate. We considered a linear relation between $T$ and $X=\left(X_1,X_2\right)$ with slopes $({\theta }_1$ and ${\theta }_2)$ and intercept $\left({\theta }_0\right)$ as the auxiliary information. We estimated $\theta$ using the standard linear regression (least square) based on a large, finite population with size, $N=$10,000. We need to generate censoring times as well to compute the event indicator, ${\delta }_i=\mathbb{I}\left(T_i\le C_i\right)$ and survival time, $Y_i=\min \left(T_i,C_i\right)$ to estimate the censored quantile regression parameters. To fit the weighted censored quantile regression model given in (1.2), we generated another $n$ observations ${\left\{y_i,x_i\right\}}_{i=1}^n$ with $n\ll N,$ using the same models given in Table 4.1. We considered the sample sizes, $n=$100 and 200 and three quantiles, $\tau =$ 0.25, 0.50, 0.75. For our proposed method, we estimated $P_i’\mathrm{s}$ using the estimating function, $g\left(t_i,x_i;\theta \right)$ defined based on the normal equations of the linear least squares method as,

$$g_i\left(z_i;\theta \right)=g\left(t_i,x_i;\theta \right)=x_i\left(t_i-x_i^{\top }\widehat{\theta }\right),i=1,2,\dots ,n.(4.1)$$

For a given quantile, $\tau ,$ the true value of the censored quantile regression parameters ${\beta }_0\left(\tau \right)$ are estimated from the population of size, $N=$10,000. In general, under a linear model assumption, the true value of the censored quantile regression slope parameters are the same as the $\theta$ (i.e., ${\beta }_1\left(\tau \right)={\theta }_1,{\beta }_2\left(\tau \right)={\theta }_2).$ But for the intercept, it is ${\beta }_0\left(\tau \right)={\theta }_0+F^{-1}\left(\tau \right),$ where $F$ is the error distribution. We conducted 1,000 simulations and computed mean bias, standard error (SE) and 95% coverage probability (CP) of the model parameter estimates for different sample sizes using 250 bootstrap samples. We compared the performance of our proposed method (CQR-EL1) with the standard censored quantile regression (CQR) model. We present the simulation results in Tables 4.2 to 4.5 respectively for Cases (i)-(iv) with ${\sigma }_{x_1,x_2}=0.$

From Tables 4.2-4.5, we see that our proposed estimator has approximately zero bias. A comparison of SE of CQR-EL1 with CQR indicates that the SE of CQR-EL1 reduces remarkably for all the parameters irrespective of any quantile. For example, we consider the scenario of $n=$100 and $\tau =$ 0.25 for comparison purposes throughout this paper. From Table 4.2, for CQR, SE of ${\widehat{\beta }}_1$ is 0.1533 and for CQR-EL1, SE of ${\widehat{\beta }}_1$ is reduced to 0.1159. When the sample size is increased to 200, SE of ${\widehat{\beta }}_1$ of our proposed method further is reduced to 0.0746. If we compare the CP of our proposed method with the nominal level of 95%, CQR-EL1 provides approximately 95% coverage and becomes more stable when the sample size increases. Similar conclusions can be reached for Case (ii) (results are in Table 4.3) even though we considered heavy tailed distribution for the failure time compared to Case (i). For example, SE of ${\widehat{\beta }}_1$ using CQR is 0.1946, whereas it is only 0.1555 for the CQR-EL1 based estimate. We also observed that SE is comparatively high in Case (ii) compared to Case (i).

In Cases (iii) and (iv), the error depends on the covariates. Simulation results for these Cases (Tables 4.4 and 4.5) are almost similar to the cases where error is independent of covariates. For example, in Case (iii) (Table 4.4), SE of ${\widehat{\beta }}_1$ is 0.0566 and 0.0424 for CQR and CQR-EL1 respectively. Similarly, in Case (iv) (Table 4.5), SE of ${\widehat{\beta }}_1$ is 0.0667 and 0.0514 for CQR and CQR-EL1 respectively. Here, we could also see a slight increase in the SE of estimates for Case (iv) because of the heavy tailed distribution assumption for the failure time compared to Case (iii).

4.2  Numerical study -II

In most of the survival data with random right censoring, the observed data are the triplet $\left\{Y=\min \left(T,C\right),X,\delta \right\}.$ We consider a linear relationship between the survival time $\left(Y\right)$ and the covariates as the auxiliary information. Here we cannot use the EL estimating function, $g(\cdot )$ defined in (4.1) because of the censoring. There are other methods available in the literature which take care of the right censoring in the linear regression.

Koul et al. (1981) introduced a synthetic data approach by transforming the survival time, $Y_r$ to a synthetic variable, ${\tilde{Y}}_r$ as

$${\tilde{Y}}_r=\frac{{\delta }_r{Y}_r}{1-G\left(Y_r\right)};r=1,2,\dots ,N,(4.2)$$

where ${\delta }_r$ is the censoring indicator and $G(\cdot )$ is the distribution of the censoring time. $E\left(\tilde{Y}|X\right)=E\left(Y|X\right)$ if $C$ is independent of both $X$ and $Y.$ When $G(\cdot )$ is unknown, we can replace it with its Kaplan-Meier estimator. The estimator of $G(\cdot )$ using the Kaplan-Meier (Kaplan and Meier, 1958) estimator is

$$1-{\widehat{G}}_N\left(t\right)={\displaystyle \prod _{r=1}^N{\left(\frac{N-r}{N-r+1}\right)}^{\mathbb{I}\left(Y_{\left(r\right)}\le t,{\delta }_{\left(r\right)}=0\right)}\text{},}(4.3)$$

where $Y_{\left(r\right)}’\mathrm{s}$ are ordered and the corresponding censoring indicator is ${\delta }_{\left(r\right)}.$ We can estimate $\theta$ as

$$\tilde{\theta }={\left(X^{\top }X\right)}^{-1}{X}^{\top }{\tilde{Y}}_r.(4.4)$$

Qin and Jing (2001) and Li and Wang (2003) independently provided the estimating function to compute the EL based data driven probabilities as

$$g_i\left(z_i;\tilde{\theta }\right)=g\left(y_i,x_i,{\delta }_i;\tilde{\theta }\right)=x_i\left({\tilde{y}}_i-x_i^{\top }\tilde{\theta }\right),i=1,2,\dots ,n.(4.5)$$

We can compute the ${\tilde{y}}_i$ and ${\widehat{G}}_n\left(t\right)$ using the sample analogues of (4.2) and (4.3) respectively.

To compute $P_i’\mathrm{s},$ we consider a linear relation between $Y$ and $X=\left(X_1,X_2\right)$ with slopes $({\theta }_1$ and ${\theta }_2)$ and intercept $\left({\theta }_0\right).$ We estimate $\theta$ using (4.4) based on a large, finite population with size, $N=$10,000. To fit the weighted censored quantile regression model given in (1.2), we generate another $n$ observations ${\left\{y_i,x_i\right\}}_{i=1}^n$ with $n\ll N$ using the same models given in Table 4.1. For our proposed method, we estimate $P_i’\mathrm{s}$ using the estimating function, $g\left(y_i,x_i,{\delta }_i;\tilde{\theta }\right)$ given in (4.5).

Similar to numerical study -I, we present the results based on 1,000 simulations and report the bias, standard error (SE) and empirical coverage probability (CP) for the nominal level of 95% based on 250 bootstrap samples. We provide the summary of the simulation results for this study in Tables 4.6-4.9.

Similar to the population information related to $T$ (numerical study -I), conclusions are almost similar for uncorrelated covariates. From Tables 4.6-4.9 we see that our proposed method (CQR-EL2) provides unbiased estimates irrespective of any sample size and quantile. If we consider the coverage probability, both CQR and CQR-EL2 provide approximately 95% coverage. For any quantile, there is a reduction in the standard error of CQR-EL2 parameter estimates compared to CQR parameter estimates. If we consider Case (i) as a basic model, CQR-EL2 with Case (ii) has reasonably higher SE along with CQR because of the heavy tailed distribution of the observed survival time. When the error depended on the covariates (Cases (iii) & (iv)), the SE of CQR-EL2 reduced considerably.

We also conducted large number of simulations with correlated covariates with ${\sigma }_{x_1,x_2}=$ 0.5 as well as constructed weights based on simple relationship with one covariate only for both numerical studies. The results of these simulations are not provided here to save the space. The conclusions arrived are almost similar to the uncorrelated covariate cases.

In numerical study -I, we noticed that there is a slight reduction in SE of ${\widehat{\beta }}_2$ using heteroscedastic models for CQR-EL1. But use of the estimating function, $g\left(y_i,x_{1i},{\delta }_i;\tilde{\theta }\right)$ (CQR-EL2), does not reduce the SE of ${\widehat{\beta }}_2$ under heteroscedastic models. Since we utilized only partial population information in relation to $X_1,$ the standard error of ${\widehat{\beta }}_0$ and ${\widehat{\beta }}_1$ reduced for CQR-EL2 compared to CQR. The standard error of ${\widehat{\beta }}_2$ was not changed.

Our simulation studies reveal that auxiliary information greatly enhances the efficiency of estimation, if the population information related to both $X_1$ and $X_2$ is available. If the population information is only related to $X_1,$ the efficiency gain is limited to ${\beta }_0$ and ${\beta }_1$ only. However, under heteroscedastic models, the efficiency of estimating ${\beta }_2$ slightly improved in numerical study -I, but not in numerical study -II.

Note that the value of the auxiliary parameter value plays a big role in the efficiency of the weighted censored quantile regression parameter estimates. If the estimate of $\theta$ based present study data and previous study (or known $\theta$ value) are very close, then all weights will be close to $1/n$ and solutions to (1.1) and (1.2) remain the same. If data on previous studies are not available, we can make of the data available in the present study to estimate the value of $\theta .$ In this case, if dimensions of $\theta$ and estimating equation $g\left(z,\theta \right)$ are same, then all weights will be equal to $1/n$ and solutions to (1.1) and (1.2) remain same. However, if the dimensions of $g\left(z,\theta \right)$ is greater than that of $\theta ,$ the weights $p\left(\widehat{\theta }\right)$ is no longer equal to $1/n$ and this scheme provides an efficiency gain over the conventional QR estimates (Tang and Leng, 2012).

4.3  Case example

The North Central Cancer Treatment Group (NCCTG) was initiated by a group of physicians from the north central region of the United States of America and the Mayo Clinic in Rochester , Minnesota . This study was conducted by NCCTG to determine whether the conclusions from the patient-completed questionnaire and those already obtained by the patient’s physician were independent or not (Loprinzi, Laurie, Wieand, Krook, Novotny, Kugler, Bartel, Law, Bateman and Klatt, 1994). They used the performance scores (ECOG and Karnofsky) to assess the patient’s daily activities. The dataset is available in the “survival” package of R software with readings of 228 patients. Because of the incompleteness of some of the variables, we had to limit the dataset to 167 observations. For the illustration of our proposed method, we changed our focus to identify the effect of following covariates over the observed survival time at different quantiles. We considered “age”, patient’s age in years; “sex”, (Male = 1 Female = 2); “ph.ecog”, ECOG performance score measured by physician (0 = good 5 = dead); “meal.cal”, calories consumed at meals and “wt.loss”, weight loss in the last six months as the covariates. After removing the incomplete patient readings, the available ECOG scores were 0,1 and 2 only. We defined two dummy categorical variables for “ph.ecog” as follows.

$$\begin{array}{l}\begin{array}{l}\text{ecog}1=\{\begin{array}{ll}1,\hfill & \text{if}\text{ph}\text{.ecog}=1\hfill \\ \text{0,}\hfill & \text{otherwise}\hfill \end{array}\hfill \\ \text{ecog2}=\{\begin{array}{ll}1,\hfill & \text{if}\text{ph}\text{.ecog}=2\hfill \\ \text{0,}\hfill & \text{otherwise}\text{.}\hfill \end{array}\hfill \end{array}\hfill \end{array}$$

To demonstrate the usefulness of our proposed method, we randomly selected a part (100 observations) of the complete data (167 observations) by considering it to be the data available from the previous study. We assumed that there exists a linear relation between the logarithm of the observed survival time and all the continuous explanatory variables (age, meal.cal and wt.loss) as the available auxiliary information. We estimated the $\theta =\left({\theta }_0,{\theta }_{\text{age}},{\theta }_{\text{meal}},{\theta }_{\text{wt}}\right)$ by the least square method based on 100 observations where the response is the synthetic variable defined by (4.2). Then we computed the EL based data driven probability weights for the present study data points (67 observations). After computing the weights, we estimated the weighted censored quantile regression parameters using Peng and Huang (2008) method with all the covariates. For the present study data, the censoring proportion is 0.283. Interestingly, we estimated the regression parameters using CQR up to the ${86}^{\text{th}}$ quantile, where as we could estimate to the ${90}^{\text{th}}$ quantile using CQR-EL2. Along with the estimates for the quantiles, $\tau =$ 0.25, 0.50, 0.75, we report standard error (SE) and 95% confidence limits using 250 bootstrap samples as well in Table 4.10.

From Table 4.10, we see that the standard error of the estimates of all the continuous variable parameters and the intercept reduced considerably because we considered the auxiliary information related to them. For the remaining variables, a reduction of standard error can also be seen, even though we did not consider any auxiliary information related to them. In the censored quantile regression with the EL based data driven probability weights, we see narrower 95% confidence limits for all the variables compared to those using the standard censored quantile regression.

• Previous
• Article main page
• Next