PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Stat Med. Author manuscript; available in PMC 2010 March 30.
Published in final edited form as:
Stat Med. 2009 March 30; 28(7): 1054–1066.
doi:  10.1002/sim.3529
PMCID: PMC2715957
NIHMSID: NIHMS112743

Regression analysis of mean quality-adjusted survival time based on pseudo-observations

SUMMARY

Regression models for the mean quality-adjusted survival time are specified from hazard functions of transitions between two states and the mean quality-adjusted survival time may be a complex function of covariates. We discuss a regression model for the mean quality-adjusted survival (QAS) time based on pseudo-observations, which has the advantage of directly modeling the effect of covariates in the QAS time. Both Monte Carlo simulations and a real data set are studied.

Keywords: quality of life, survival analysis, multistate models, pseudo-observations, jackknife, regression model

1. INTRODUCTION

In many clinical studies the primary outcome is the time to some event, usually survival, of patients. Of interest is often the comparison of treatments or the modeling of factors possibly affecting outcome. Recently it has become of interest to consider the quality of the survival experience as well as its length. Quality of life may depend on physical, functional, mental and social well-being in addition to survival and disease-specific responses.

Quality-adjusted survival (QAS) time is an approach that attempts to account for both the length and quality of a patients survival experience in a single variable. It is an extension of the Q-TWiST method, proposed by Gelber et al. [1].

In order to define the QAS time, we assume that the health history of a given patient can be described by a process {V(t), t≥0}, where V(t) may assume any of K + 1 health states in which a patient may belong. We denote the space state by Γ = {0,1,…,K}. States 1,2,…,K are transient health states and state 0 is an absorbing state, typically death. Notice that the usual survival time is given by T = inf{t: V(t) = 0}. We must also define the function Q that maps the state space into a pre-specified set of real numbers, know as utility coefficients, ranging in the interval [0,1]. The utility coefficient associated with the absorbing state is zero. The perfect health state coefficient is taken to be one and the utility coefficients of the other health states reflect the health experience of patients in that state. A utility coefficient of 12 for a health state, for example, would imply that a year of survival in this state is the equivalent of 6 months in the perfect health state.

The QAS time up to τ is then defined as

U=0TτQ{V(t)}dt
(1)

The QAS time can also be written as

U=q1Y1++qKYK
(2)

where Yi is the sojourn time in state i and qi = Q(i) is the utility coefficient associated with the health state i.

Modeling of the QAS function, defined as SQ(t) = P[U >t], is difficult when the data are right censored. The censoring in this setting is informative as pointed out by Gelber et al. [1], since patients with poor quality of life tend to accumulate QAS slower than the ones with very good quality of life. Thus, the Kaplan–Meier estimator of SQ is a biased estimator. Zhao and Tsiatis [2] have proposed an efficient estimator for the distribution of the QAS time.

Approximately unbiased estimators of the mean QAS restricted to τ,μτ=0τufU(u)du, for censored quality of life data are available. This includes estimators proposed by Huang and Louis [3] and by Zhao and Tsiatis [4]. These models do not account for covariates. Recently, Tunes-da-Silva et al. [5] derived an expression for the covariate adjusted mean QAS based on multistate models. The intensity of transitions between health states in this model was modeled by a proportional hazards model so that the mean QAS can be estimated for a given vector of covariates. While this model is useful for evaluating the effects of covariates on the transition intensities from one state to the next, the effect of a covariate on the mean survival itself is hard to determine since it is a highly non-linear effect based on all the transition intensities.

In this paper we present a technique to model directly the effects of covariates on the restricted mean QAS time. The technique is based on using a generalized estimating equation (GEE) model on the pseudo-values of the restricted mean QAS time. The pseudo-values are obtained as differences of the complete sample and leave-one-out estimator of the mean QAS time. This approach has been used successfully to model multistate model probabilities (see Andersen et al. [6], Andersen and Klein [7]); to model restricted mean survival (see Andersen et al. [8]); or to model competing risks data [9]. An alternative approach to direct modeling of the mean QAS up to time τ can be obtained by modifying the approach of Bang and Tsiatis [10] to model median costs with censored data. The approach is based on an estimating equation derived using the inverse probability weighting technique, given by

S(β)=i=1nΔiK^(Tiτ)Zi(UiβTZi)=0
(3)

where Ui is the QAS up to time τ for the ith patient, Δi is the indicator that individual i dies at Ti or is alive at time τ and [K with circumflex](·) is the Kaplan–Meier estimator of the censoring distribution. Estimates of the standard errors of [beta] that maximizes (3) can be found by a bootstrap technique.

To illustrate our approach we consider two examples taken from analysis of the bone marrow transplantation (BMT) recovery process. The first example involves a set of 614 chronic myeloid leukemia (CML) patients treated at the Hammersmith hospital in London and first reported in Klein et al. [11]. Here patients right after transplantation are in complete remission of their CML. At a future time they relapse and their utility of quality of life decreases. To induce a second remission patients are given an infusion of donor leukocytes (a DLI) to induce a second remission using the so-called graft-versus-leukemia effect (see Horowitz et al. [12]). Most patients then have a durable second remission. Of course patients may die from any of the health states as depicted in Figure 1. Two sets of utilities will be used to model a patient’s QAS. One utility assigns a weight of one to the two remission states and zero to all other states and gives an expected QAS analogous to the expected TWiST of Gelber et al. [1]. A second set of utilities was obtained by polling 10 transplant physicians at the Center for International Blood and Marrow Transplantation (CIBMTR). The mean utilities were: First remission, q = 1; First relapse, no DLI, q =0.76; First Relapse after DLI, q =0.70; Second Remission, q =0.89; and Second Relapse, q =0.61. Of interest is how the mean QAS is related to the choice of donor (HLA identical sibling or unrelated donor), age at transplant, disease status at transplant and graft-versus-host prophylaxis.

Figure 1
CML data process.

The second example is based on a study of 517 BMT patients performed to evaluate grading systems for acute graft-versus-host disease (aGVHD) reported in Cahn et al. [13]. Patients were assessed for aGVHD weekly for the first 14 weeks after transplant. The severity of the aGVHD was assessed weekly using the Glucksberg grade. This grade is based on the severity and involvement of three organ systems (skin, gastrointestinal track and liver) and ranges from zero (no aGVHD) to four. Our panel of CIBMTR physicians assigned utilities of heath to the grades of aGVHD as follows: for grade 0, the utility coefficient is q = 1; for grade 1, q =0.94; for grade II, q =0.82; for grade III, q =0.63; and for grade IV, q =0.43. Of interest is the mean QAS in the first 14 weeks post BMT as a function of donor type.

The remainder of the paper is as follows. In Section 2 we review univariate estimators of the mean QAS time. In Section 3 we use these estimators in the pseudo-value approach to regression modeling. In Section 4 we present results of a Monte Carlo study of the performance of the regression estimators. In Section 5 we revisit the two examples and in the final section we present some concluding remarks.

2. ESTIMATORS FOR THE MEAN QAS TIME

In this section we briefly review the estimation of the mean QAS time truncated at time τ. The pseudo-value approach discussed in the next section requires an (approximate unbiased) estimator of this quantity and a number of estimators are available [3,4]. To allow for multiple entries into a health state and because we can use standard statistical software to compute the estimator we will base our pseudo-values on the estimator of Huang and Louis [3].

To describe the estimator we define a mode as the description of the state and the number of times that state was visited and Θ as the set of all modes. For example, in our first example (see Figure 1) possible modes are ‘1st remission’, ‘relapse for first time’, ‘DLI following relapse’, ‘2nd remission’, ‘2nd relapse’, and ‘dead’. Associated with each potential mode, J[set membership]Θ is the time XJE, when a patient first reaches mode J and a time XJL, when the patient leaves mode J. Note that XJE and XJL are measured from time 0 and XJLXJE is the sojourn time in mode J.

The mean QAS time can be written as

μ=JΘqJE(XJLXJE)

where qJ is the utility associated with mode J.

It can be shown that

μ=JΘI=E,LqJI0GJI(u)du

where

qJI={qJifI=EqJifI=L

and GJI(t)=P(XJI>t),I=E,L.

This expression motivates the ‘event-marginal estimator’, given by

μ^τ=JΘ,I=E,LqJI[0,τ]G^J(t)dt
(4)

where ĜJ(t) is the Kaplan–Meier estimator of GJ(t). Note that here and in the sequel we restrict the mean and the integrals to the range (0, τ) to ensure that all estimators are finite.

3. PSEUDO-VALUE REGRESSION FOR MEAN QAS

In this section we present a regression model for the mean QAS time using the pseudo-value regression approach in Andersen et al. [6]. In this approach we let Y1, Y2,…,Yn be independent and identically distributed random variables and θ a parameter of the form

θ=E[f(Yi)]

Suppose that we have an (approximately) unbiased estimator [theta w/ hat] for this parameter. Let Zi, i = 1,…,n, be independent and identically distributed sample of covariates. Define the conditional expectation by

θi=E[f(Yi)|Zi]

The ith pseudo-observation is

θ^i=nθ^(n1)θ^i

where [theta w/ hat]i is the ‘leave-one-out’ estimator for θ based on Yj, ji. Note that if there is no censoring one can show that [theta w/ hat]i = f(Yi).

The pseudo-observation can be shown to have the same expectation as θi [6]. A regression model for θ can be based on a generalized linear model for θi, namely,

g(θi)=βTZi

where g(·) is a link function and β is a set of regression coefficients. Standard GEE packages can be used to fit these models.

Our goal is to apply this approach to a regression model for the restricted mean quality of life, μτ. An identity link function is most appropriate so our model is

μτ|Z=βTZ

where we take the first element of Z to be one to allow for an intercept term.

Estimates of β are based on the pseudo-values observations, [theta w/ hat]i derived using [mu]τ[3]. The GEE is

U(β)=i=1nUi(β)=i=1n((βTZi)β)(θ^iβTZi)
(5)

As shown in Andersen et al. [6], consistent estimators of β are obtained as the solution to (5) and the variance of β can be obtained from a standard sandwich estimator

Σ=I1(β^)V(U(β^))I1(β^)

where

I(β)=i(βTZiβ)T(βTZiβ)

and

V(U(β))=iUi(β)UiT(β)

Hence once we compute the pseudo-values, estimation of β and Σ can be carried out using standard statistical software such as SAS’s Proc Genmod or R’s GEE.

4. MONTE CARLO STUDY

To study the performance of this approach we performed a small Monte Carlo study. The study was focused on the effect of censoring on the regression estimates and on how the choice of estimator for the restricted mean QAS affects the estimation of model parameters.

In the study we generated data from a three state forward transition model (see Figure 2). We considered models with two covariates, one, Z1, Bernoulli distributed with p=12 and the other, Z2, uniformly distributed on [−1, 1]. Transition intensities were assumed to be of the form λkl=exp{βklTZ} where Z ={1,Z1,Z2}. In the study we used: β12=(0.4,1, − 0.8),β10 = (0.2,−0.5,0.5), β23 = (0.4,1,−0.8), β20 = (0.2,−0.5,0.5)and β30 = (0.2,−1,0.8). The utility coefficients used in the simulation are q1 = 1,q2 = 0.6 and q3=0.8.

Figure 2
Simulated process.

As noted earlier it is hard to determine the regression model for the mean QAS based on regression models for the transition intensities. To obtain ‘true’ values for this model we use an approach suggested by Andersen et al. [8]. We first create a set of 402 Z vectors with Z1,i = 1 for the first half and 0 for the second half. We took Z2,i =−1 + (i − 1)/100 for the first half of the sample and Z2,i+201 = −1 + (i − 1)/100 for the second half. Using the theoretical expressions for the mean QAS (see Tunes-da Silva et al. [5]) we compute the value of μL for each of the Z vectors. Here we take τ= 12. We then fit a linear regression model to these values of μτ used as ‘true’ values in the tables.

We simulated data with samples sizes of n = 100 and n=400. We generated the censoring in two different ways so that we could have a better idea of the effect of censoring on estimators properties. In the first simulation study, exponential censoring was used to obtain data with 0, 10 and 20 per cent censoring [5] and the same censoring distribution was used to all observations. In the second simulation study, we generated exponential censoring as well, but with different hazard rates for the two groups (so that we would have the same proportion of censored observations in each group).We examined two approximately unbiased estimators of mean QAS, the Huang and Louis [3] estimator discussed in Section 2 and an estimator of Zhao and Tsiatis [4], and applied the pseudo-values technique for the regression model. We also considered the inverse probability-weighted estimating equation approach (3).

A total of 1000 replicates were used in each scenario of our simulation study. Table I examines the bias of the regression estimates and also report SDsim, the sample standard deviation of the 1000 estimates of β and SEest = square root of the mean of the 1000 variance estimates from the GEE model fitting. We did not compute the variance estimatives of parameters for the inverse probability-weighted estimating equation because no estimator for that has been proposed. Here we see that there is less variability when using the Huang and Louis method than the Zhao and Tsiatis method. The sandwich estimator of the variance seems to perform reasonably well at least for large samples. The results of the pseudo-observation regression seem to be quite robust to the choice of estimation technique with the Huang and Louis approach a bit better with more censoring. The method based on the inverse probability weighting gives comparable results to the pseudo-observation approach. However, given that the linear model is only approximately true here it is a bit hard to decide between the models. A simulation study was also conducted including a normally distributed covariate and the conclusions are very similar to the ones presented here (results not shown). Figure 3 examines how well the models predict the true mean QAS. An approximate value of true mean based on theoretical calculations was computed and the bias is given by the predicted value for a given estimator minus the true value. The lines are the mean of 1000 estimated bias based on the Huang and Louis estimator, Zhao and Tsiatis estimator and also for the estimator computed using the inverse probability-weighted estimating equation for the situation with 20 per cent of censoring.

Figure 3
Simulation results for 20 per cent of censoring.
Table I
Simulation results: estimated average of 1000 parameter estimates; SDsim square root of empirical variance of 1000 simulation replications; and SEest square root of average estimated variance based on the generalized linear model.

5. EXAMPLES

We now examine the use of these techniques in the two examples. The first example dealt with CML patients given the potential for a second remission using a DLI (see Figure 1). Here there were 614 patients initially in remission after transplant (State 1) of which 246 relapse (State 2) and 202 died. Of the relapsed patients 124 had a DLI (State 3) and 92 died. Of those given a DLI, 77 had a second remission (State 4) of which 5 died and 8 were still in the second relapse state at the end of the study. Covariates were age (mean 34 range 4–59), donor type (243 HLA identical siblings and 371 matched unrelated), disease status (early 453, intermediate 138 and advanced 23) and GVHD prophylaxis (CSA + MTX 252, CSA +Other 68, ATG 208, and T-depletion 86). Using the pseudo-regression approach on the mean QAS restricted to 5 years we fit three models. In the first we take qi = 1 for all health states so the mean QAS reduces to the mean time to death. We also fit the model with qi = 1 for the two remission states and q = 0 for all other states. This gives an estimator analogous to the mean TWiST. Finally, we fit a model with the consensus utilities: q1 = 1; q2 = 0.76; q3 =0.70; q4 = 0.89 and q5 =0.61. The results are in Table II. Note all times are in months. Here we see that there a significant drop of about 8.8 months (95 per cent) confidence interval (−16.5 to −0.5 months) in overall mean survival for patients with an unrelated transplant. However, for the mean time spent in the remission states (TWiST), there is no significance difference between donor types. For the mean QAS using the consensus weights, the effect of donor is more pronounced but it does not quite obtain significance at the 5 per cent level.

Table II
Regression models for 612 CML patients given a DLI, identity link function.

For this model we used an identity link function so that the effect of covariates is additive on a baseline mean quality-adjusted life time. The values β(Z1Z2) can be interpreted as the excess mean QAL for a patient with covariate Z1 as compared with a patient with covariate Z2. While any link function can be used, a link function that leads to a nice interpretation is the log link function, g(x)=log(X). Here we have that μ(Z) = exp(α+βZ)and we can interpret eα as a baseline mean QAS and exp{β(Z1Z2)} as the relative excess mean QAL for a patient with covariate Z1 as compared with a patient with covariate Z2. Table III shows the results of fitting the three models using the log link function. As one can see when this table is compared with Table II the basic conclusions remain the same, namely the mean quality-adjusted life is better for early disease patients and patients given CSA+MTX. In Table II we see that patient’s with intermediate disease have about 9 months shorter mean TWiST than an early stage patient, while in Table III we see they have a mean TWiST that is 74 per cent of that for an early disease patient, for example.

Table III
Regression models for 612 CML patients given a DLI, log link function.

The second example is based on study of 517 patients given a BMT [13]. Patients had their grade of acute GVHD assessed weekly for the first 14 weeks of the study. Patients could have no GVHD or grade I–IV GVHD. For simplicity, we assume that a patient’s GVHD status changes only on the first day of the week and that patients remain in that GVHD state until the first day of the corresponding week or until they die. The health utility for each state was assessed by our panel of transplant physicians as q = 1 for grade 0, q=0.94 for grade I, q=0.82 for grade II, q = 0.63 for grade III and q=0.43 for grade IV Covariates of interest are the patient’s age at transplantation, donor type (matched related, matched unrelated, mismatched related, mismatched unrelated), graft type (bone marrow or blood stem cells) and disease. In Table IV we report the results of a Cox model analysis of the risk factors for 14-month survival. Here we ignore aGVHD and we censor patient at 14 weeks. The results show that 14-week survival is related to age and donor. Table IV also shows the regression model for the mean QAS up to 14 weeks. We computed parameter estimatives using both the pseudo-values approach and the inverse probability-weighted estimating equation. Since no variance estimator for parameter estimatives is available, we propose here the use of bootstrap to compute confidence intervals. We used the non-parametric bootstrap-t method [14] based on 800 bootstrap samples. Here we see that the Cox model provides similar results to the quality-adjusted approach. Donor type and age seem to be the most important variables concerning the QAS time as well as survival. We see a significant drop out in the mean QAS time for patients with unrelated transplant (both matched and mismatched).

Table IV
Results for mortality and Quality-adjusted survival restricted to 14 weeks post transplantation.

DISCUSSION

We have presented a flexible method for modeling directly the effects of covariates on the expected QAS time. The few methods currently available for analyzing mean QAS are based on models for the effects of covariates on the transition intensities. The mean QAS time is a complex function of these transition intensities and it is extremely difficult to understand how some covariates influence the expected QAS time. In this paper we have presented an alternative to these techniques based on pseudo-observations. This approach can be applied with any (approximately) unbiased univariate estimator of mean quality-adjusted time. Hence we directly model the effect of covariates on the mean QAS time. This method allows the use of standard generalized linear models software to estimate parameters and their standard errors once the pseudo-values are computed. All that is required is a single computation of the pseudo-values. In our simulation study we found that the approach is relatively robust to the choice of the underlying estimator of the expected QAS time. The estimator proposed by Huang and Louis [3], which is based on the difference of Kaplan–Meier curves, is a simple estimator that performed well in our simulations studies. The Zhao and Tsiatis [4] estimator depends on the censoring distribution and it may be more effected by censoring than Huang and Louis estimator. The method based on the inverse probability weighting method seems to give comparable results to our method. This method as opposed to the pseudo-value method requires the development of special software. The variance of the regression estimates is obtained by bootstrapping, which can be quite computer intensive. The pseudo-observation technique once the values are computed uses standard GEE software for estimation of the regression coefficients and to compute the sandwich estimator of the standard error. This seems to be an advantage for our approach.

Using the results presented here an investigator can obtain a direct regression model for the mean QAL time of an individual. By modeling each of the transition intensities an investigator, using methods in Tunes-da-Silva [5, 15], can obtain a regression model for the mean QAL of an individual. Both methods have advantages and disadvantages. The first approach allows the user to have a model for the aggregate effect of a covariate on mean QAL which is easy to interpret. This simple model can be used to make predictions for patients of their mean QAL based on their characteristics at baseline. It does this modeling at the expense of having models for each of the transition intensities that may be interest to the clinical investigator trying to determine the mechanism of good mean QAL. This is a strength of the second model. Of course when risk factors have quite different effects for different transitions this is hard to interpret. Both approaches require additional programming to implement, the first a routine to compute pseudo-observations and the second to synthesize the Cox models. The direct modeling approach may be a bit more robust since it does not require a Markov assumption on the states transition intensities, but if the goal is to predict mean QAL both methods perform well.

ACKNOWLEDGEMENTS

Professor Tunes-da-Silva’s work was supported by Fundação de Amparo à Pesquisa do Estado de São Paulo—FAPESP, Brazil, grant number 2007/02823-3. Professor Klein’s work was supported by a grant from the National Cancer Institute (R01 CA54706-14). We thank the referees for their helpful comments and for suggesting the Bang and Tsiatis estimator.

REFERENCES

1. Gelber RD, Gelman RS, Goldhrisch A. A quality-of-life-oriented endpoint for comparing therapies. Biometrics. 1989;45(3):781–795. [PubMed]
2. Zhao H, Tsiatis AA. Efficient estimation of the distribution of quality-adjusted survival time. Biometrics. 1999;55:1101–1107. [PubMed]
3. Huang YJ, Louis TA. Expressing estimators of expected quality-adjusted survival as functions of Nelson–Aalen estimators. Lifetime Data Analysis. 1999;5(3):199–212. [PubMed]
4. Zhao H, Tsiatis AA. Estimating mean quality-adjusted lifetime with censored data. Sankhya: The Indian Journal of Statistics, Series B. 2000;62:175–188.
5. Tunes-da-Silva G, Sen PK, Pedroso-de-Lima AC. Estimation of the mean quality-adjusted survival using a multistate model for the sojourn times. Journal of Statistical Planning and Inference. 2008;138(8):2267–2282.
6. Andersen PK, Klein JP, Rosthoj S. Generalized linear models for correlated pseudo-observations, with applications to multi-state models. Biometrika. 2003;90:15–27.
7. Andersen PK, Klein JP. Regression analysis for multistate models based on a pseudo-value approach, with applications to bone marrow transplantation studies. Scandinavian Journal of Statistics. 2007;34:3–16.
8. Andersen PK, Hansen MG, Klein JP. Regression analysis of restricted mean survival time based on pseudo-observations. Lifetime Data Analysis. 2004;10:335–350. [PubMed]
9. Klein JP, Andersen PK. Regression modeling of competing risks data based on pseudo-values of the cumulative incidence function. Biometrics. 2005;61:223–229. [PubMed]
10. Bang H, Tsiatis AA. Median regression with censored cost data. Biometrics. 2002;58:643–649. [PubMed]
11. Klein JP, Szydlo RM, Craddock C, Goldman JM. Estimation of current leukaemia-free survival following donor lymphocyte infusion therapy for patients with leukaemia who relapse after allografting: application of a multistate model. Statistics in Medicine. 2000;19:3005–3016. [PubMed]
12. Horowitz MM, Gale RP, Sondel PM, Goldman JM, Kersey J, Kolb HJ, Rimm AA, Ringden O, Rozman C, Speck B, Truitt RL, Zwaan FE, Bortin MM. Graft-versus-leukemia reactions after bone marrow transplantation. Blood. 1990;75(3):555–556. [PubMed]
13. Cahn JY, Klein JP, Lee SJ, Milpied NL, Blaise D, Antin JH, Leblond V, Ifrah N, Jouet JP, Loberiza F, Ringden O, Barrett AJ, Horowitz MM, Socie G. Prospective evaluation of 2 acute graft-versus-host (GVHD) grading systems: a joint Société Française de Greffe de Moelle et Thérapie Cellulaire (SFGM-TC), Dana Farber Cancer Institute (DFCI) and International Bone Marrow Transplant Registry (IBMTR) prospective study. Blood. 2005;106(4):1495–1500. [PubMed]
14. Efron B, Tibshirani RJ. An Introduction to the Bootstrap. Boca Raton, FL: CRC Press; 1993.
15. Tunes-da-Silva G, Sen PK, Pedroso-de-Lima AC. A semi-Markov multistate model for estimation of the mean quality-adjusted survival for non-progressive processes. Lifetime Data Analysis [PubMed]