Skip Navigation


Biostatistics Advance Access originally published online on April 14, 2005
Biostatistics 2005 6(3):479-485; doi:10.1093/biostatistics/kxi023
This Article
Right arrow Abstract Freely available
Right arrow FREE Full Text (PDF) Freely available
Right arrow All Versions of this Article:
6/3/479    most recent
kxi023v2
kxi023v1
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Add to My Personal Archive
Right arrow Download to citation manager
Right arrowRequest Permissions
Right arrow Disclaimer
Google Scholar
Right arrow Articles by Sousa, I.
Right arrow Articles by Diggle, P. J.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Sousa, I.
Right arrow Articles by Diggle, P. J.
Social Bookmarking
Add to CiteULike   Add to Connotea   Add to Del.icio.us  
What's this?

© The Author 2005. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: journals.permissions@oupjournals.org.

Exploratory analysis of longitudinal trials with staggered intervention times

Inês Sousa*, Amanda G. Chetwynd and Peter J. Diggle

Department of Mathematics and Statistics, Lancaster University, LA1 4YF, Lancaster, UK i.pereirasilvacunhasousa{at}lancaster.ac.uk

* To whom correspondence should be addressed.


   SUMMARY
TOP
 SUMMARY
1. INTRODUCTION
 2. A NONPARAMETRIC ESTIMATOR...
 3. ANALYSIS OF THE...
 4. DISCUSSION
 REFERENCES
 
Longitudinal trials involving surgical interventions commonly have subject-specific intervention times, due to constraints on the availability of surgeons and operating theatres. Moreover, the intervention often effects a discontinuous change in the mean response. We propose a nonparametric estimator for the mean response profile of longitudinal data with staggered intervention times and a discontinuity at the times of intervention, as an exploratory tool to assist the formulation of a suitable parametric model. We use an adaptation of the standard generalized additive model algorithm for estimation, with smoothing constants chosen by a cross-validation criterion. We illustrate the method using longitudinal data from a trial to assess the effect of lung resection surgery in the treatment of emphysema patients.

Keywords: Back-fitting algorithm; Cross-validation; Exploratory analysis; Longitudinal trials; Lung resection surgery; Nonparametric estimator


   1. INTRODUCTION
TOP
 SUMMARY
 1. INTRODUCTION
2. A NONPARAMETRIC ESTIMATOR...
 3. ANALYSIS OF THE...
 4. DISCUSSION
 REFERENCES
 
A common protocol in trials to evaluate the effectiveness of nonacute surgical interventions is the following. At time zero, subjects recruited to the trial are randomized to receive standard medical therapy, or medical therapy plus a surgical procedure. For convenience, we refer to the two arms of the trial as medical and surgical, respectively. Measurements are taken on relevant outcomes for each patient at a sequence of follow-up times, and the objective is to compare the mean response profiles over time for subjects in the medical and surgical treatment arms. In principle, the schedule of follow-up times is prespecified, but in practice the actual follow-up times inevitably vary between subjects. More fundamentally, the following two features are typical, and have motivated us to develop a novel method of exploratory analysis: firstly, surgical intervention times necessarily vary between subjects and the actual times of surgery may, in extreme cases, be toward or even beyond the end of the prespecified period of follow-up; secondly, any method for estimating the mean response profile in the surgical arm needs to allow for a discontinuity at the time of surgery.

An example is provided by a trial to assess the effect of lung resection surgery as an addition to standard medical therapy in the treatment of emphysema patients (Lim et al., 2004Go). The primary outcome variable is forced expiratory volume in one second (FEV1). Forty-eight patients were recruited to the trial and randomized between medical and surgical treatment arms. The trial protocol initially specified that measurements of FEV1 would be taken at 3, 6, and 12 months, where zero denotes the time of randomization. After 12 months, follow-up continued but in an irregular pattern. Surgery times varied between 3 days and 26.4 months, with an average of 7.6 months.

In these circumstances, an intention-to-treat analysis (Pocock, 1983Go) could be very misleading. We therefore favor a model-based analysis as described in Diggle et al. (2002Go, Chapters 4 and 5, henceforth DHLZ) and in Molenberghs et al. (2004)Go, but adapted to take account of the staggered intervention times within the surgical arm.

For an outcome variable Y, let Yij be the jth follow-up measurement of Y on the ith subject, tij the time at which this measurement is taken, and µij = E[Yij]. For a subject randomized to the surgical arm, let si be the time of their surgery. Ignoring for the time being any adjustment for baseline covariates, we assume that for subjects in the medical arm (nonsurgical), µij = µ0(tij), whereas for subjects in the surgical arm,

(1.1)

so that the function {delta}(u) represents the mean effect of the surgical intervention as a function of time since surgery. We assume that µ0(t) is smoothly varying and that {delta}(u) is smoothly varying except for a discontinuity at u = 0 together with the constraint that {delta}(u) = 0 for all u < 0. In this paper, we propose nonparametric estimators for the functions µ0(·) and {delta}(·) which together determine the mean response profiles. As described in DHLZ, the residuals from the estimated mean response profiles can then be used to suggest a suitable parametric model for the covariance structure of the data, while the estimator itself may suggest a suitable parametric model for the mean response profiles. Formal inference can then proceed using likelihood-based methods in conjunction with the resulting parametric model.

In Section 2 of this paper, we propose nonparametric estimators of the two functions of interest, µ0(·) and {delta}(·), as solutions to a penalized least squares criterion. We also suggest a cross-validation method for choosing values for the corresponding smoothing constants. In Section 3, we report the results of applying the method to the lung resection surgery trial data, showing in particular how the nonparametric estimator assists in the formulation of a suitable parametric model. A set of R functions which implement the methods can be downloaded from http://www.maths.lancs.ac.uk/~pereiras.


   2. A NONPARAMETRIC ESTIMATOR FOR THE MEAN RESPONSE PROFILES
TOP
 SUMMARY
 1. INTRODUCTION
 2. A NONPARAMETRIC ESTIMATOR...
3. ANALYSIS OF THE...
 4. DISCUSSION
 REFERENCES
 
We propose to estimate the functions µ0(·) and {delta}(·) in (1.1) by spline smoothing within a generalized additive model (GAM, Hastie and Tibshirani, 1990Go). However, we cannot use the standard GAM algorithm, as implemented for example in the R function ‘gam’ because of the combination of staggered surgical intervention times and the discontinuity in {delta}(·). To do so would require us to express the complete set of mean responses for the data as E[yij] = µ0(x1ij) + {delta}(x2ij), where x1ij and x2ij are the values of explanatory variables x1 and x2 attached to the jth observation in the ith subject, and µ0(·) and {delta}(·) are smoothly varying functions. Setting x1ij = tij meets the requirement for µ0(·). Setting x2ij = 0 for all j and subject i in the medical arm, and x2ij = max( 0, tij – si) for subject i in the surgical arm, gives the correct interpretation for E[yij] but does not allow {delta}(u) to be discontinuous at u = 0, and would result in all preintervention measurements contributing indirectly to the estimation of {delta}(·). We therefore propose an adaptation of the standard GAM algorithm, as follows.

Let m denote the number of subjects and ni the number of follow-up measurements for the ith subject. We estimate µ0(·) and {delta}(·) to minimize a penalized least squares criterion,

(2.1)

An adaptation of the argument in Green and Silverman (1994Go, Chapter 2) shows that the resulting estimators and are cubic splines with knots at the measurement times tij and at the shifted measurement times tij – si postsurgery, respectively. Given values for {lambda}1 and {lambda}2, we use a back-fitting algorithm to obtain the two estimates, alternating between minimizing S(·) with respect to µ0(·) holding {delta}(·) fixed at its current estimate and vice versa. In our experience, the algorithm typically converges within about 10–15 iterations.

To choose the values of the smoothing constants {lambda}1 and {lambda}2, we minimize a cross-validated least squares criterion in the spirit of Rice and Silverman (1991)Go. Write {lambda} = ({lambda}1, {lambda}2) and let denote the estimate of µij obtained by omitting all observations from subject i, with {lambda} held fixed. Then, the cross-validated least squares criterion is

(2.2)

As explained in Rice and Silverman (1991)Go, minimizing (2.2) leads to consistent estimation of the bandwidth values {lambda} which minimize the mean square error, whatever the within-subject correlation structure. Moreover, evaluation of (2.2) does not require explicit computation of all m leave-one-out estimates of the mean response. For fixed {lambda} and any t, the estimates of µ0(t) and {delta}(t) at the current iteration within the back-fitting algorithm can be written as

and

where, for all t, wij(t) = 0 if subject i is in the medical arm, and if subject i is in the surgical arm and tij < si, i.e. the corresponding Yij is a presurgery measurement. It follows that for any subject i in the medical arm,

where For a subject i in the surgical arm, we need to distinguish between measurements observed before (tij < si) and after (tij ≥ si) the surgical intervention. Then, for tij < si,

and for tij ≥ si,

where


   3. ANALYSIS OF THE LUNG RESECTION TRIAL DATA
TOP
 SUMMARY
 1. INTRODUCTION
 2. A NONPARAMETRIC ESTIMATOR...
 3. ANALYSIS OF THE...
4. DISCUSSION
 REFERENCES
 
We now consider the data from the lung resection surgery trial. One of the subjects randomized to the surgical arm declined to proceed with surgery, while seven randomized to the medical arm subsequently underwent surgery. Within our analysis framework, the subject who declined with surgery presents no special difficulties because the results are unaffected by whether we regard them as right censored prior to surgery or as being transferred to the medical arm. For the subjects who crossed over to surgery, the situation is less straightforward because the decision to undertake surgery was influenced by their good general clinical condition, potentially leading to selection bias affecting our results. We considered three scenarios with respect to these patients: retain them in the medical arm but censor each at the time of their surgery, retain and analyze them as part of the medical arm (intention-to-treat), and consider them as being transferred into the surgical arm. For the analysis reported here, we used the first of these three strategies. In principle, it would be wiser to analyze the data in such a way as to allow for the possibility of informative noncompliance (Smith and Diggle, 1999Go). However, for these data, our results were not materially affected by the three scenarios we used to deal with the noncompliant subjects, because their times of surgery mostly came toward the end of the follow-up period; for example, only one of the seven received surgery during the first 12 months of follow-up. For further discussion of ways of dealing with deviations from protocol in randomized clinical trials, see for example Matthews (2000Go, Chapter 9) or Senn (2002Go, Chapter 9).

For the minimization of (2.2), we used an initial grid-search followed by the Nelder–Mead simplex algorithm (Nelder and Mead, 1965Go) as implemented by the R function optim( ), using four widely separated starting values and checking that convergence was to the same solution in each case. In the initial grid-search, we found no evidence of secondary local minima in the CV({lambda}) surface.

Using this value of {lambda}, we then estimated the functions µ0(·) and {delta}(·) using the back-fitting algorithm. Figure 2 shows the resulting estimates. The form of is somewhat complicated, perhaps reflecting chance fluctuations in the general health of subjects within this rather small sample. In contrast, is approximately linear, suggesting that a linear form, {delta}(t) = {alpha} + ßt, may be adequate in a parametric analysis. This has the advantage that the parameters {alpha} and ß have natural interpretations as the immediate average benefit of surgery (assuming, as appears to be the case, that {alpha} > 0) and the rate over time at which this initial benefit is lost (similarly assuming ß < 0), respectively. Of course, extrapolation beyond, or even as far as, the maximum observed follow-up time would be unwarranted.



View larger version (15K):
[in this window]
[in a new window]
  		Fig. 2. Nonparametric estimates and, together with pointwise 95% envelopes (solid lines) from 500 simulated data sets (see text for details of simulation model), and parametric estimate for {delta}(t) (dashed line).

 
To obtain interval estimates of µ0(t) and {delta}(t), we proceed as follows. Firstly, we construct residuals for each subject, , where the are defined by (1.1) but using the nonparametric estimates of µ0(·) and {delta}(·). Secondly, we construct the empirical variogram of these residuals (DHLZ, Section 3.4) and use nonlinear ordinary least squares to identify a suitable parametric model for the covariance structure of the residuals. Figure 1 shows the empirical and fitted variograms, the latter taking the form

(3.1)

with {tau}2 = 0.016, {sigma}2 = 0.062, and {phi} = 43.018.



View larger version (18K):
[in this window]
[in a new window]
  		Fig. 1. Empirical (solid line) and fitted (dashed line) variograms of residuals from the nonparametric estimates of the mean response.

 
Finally, we simulated replicate data sets from a Gaussian linear model (DHLZ, Chapters 4 and 5) with observation times equal to those for the actual data, mean response profiles derived from the nonparametric estimates of µ0(·) and {delta}(·), and covariance structure derived from the fitted variogram (3.1). Figure 2 includes pointwise 95% envelopes from these simulated data sets. These suggest that much of the fine structure of is attributable to statistical error, but more particularly that a linear form for {delta}(·) should suffice. The second conclusion is important because, as noted earlier, it enables a straightforward interpretation of the effect of the surgical intervention. The function µ0(·) is not of direct interest, and there is no interpretive advantage to its having a simple parametric form. Note also that the widths of the simulation envelopes increase with time, as the data are progressively censored.


   4. DISCUSSION
TOP
 SUMMARY
 1. INTRODUCTION
 2. A NONPARAMETRIC ESTIMATOR...
 3. ANALYSIS OF THE...
 4. DISCUSSION
REFERENCES
 
The specific contribution of this paper has been to propose nonparametric estimators for the mean response profiles in longitudinal studies involving an intervention for which two conditions apply: the intervention may effect a discontinuous change in the mean response; intervention times are subject specific. We believe that both conditions commonly hold in studies involving surgical interventions, which are necessarily scheduled according to the availability of surgeons and operating theatres. For the estimation, we have developed an adaptation of the standard back-fitting algorithm used in generalized additive modeling to obtain smoothing spline estimates of the functions defining the mean response, with smoothing parameters selected by a cross-validation criterion which allows for arbitrary covariance structure in the sequence of longitudinal measurements from a single subject. To assess the precision of our estimates, we have used a simple Monte Carlo method, combining our nonparametric estimates of the mean response with a parametric, correlated Gaussian model for residual variation about the mean.

Our use of this method has been as an exploratory tool, to help in the formulation of a parametric model for the data which we could then use for formal inference. In the case of the lung surgery data used in the paper, we subsequently fitted a parametric model in which the mean response in the absence of surgery, which was not of particular scientific interest, was modeled as a quartic polynomial in time while the surgical intervention effect, which was the primary focus of interest, was modeled linearly. This gave us an efficient analysis which could be reported in terms of two physically interpretable parameters: the average immediate benefit of surgery and the rate at which, on average, this benefit reduced over time. Specifically, we estimated the immediate benefit to be a 0.259l increase in mean FEV1, with associated 95% confidence interval (0.179, 0.339), and the mean rate of change postsurgery as – 0.005l per month, with 95% confidence interval ( – 0.009, – 0.001). Note that the parametric estimate of {delta}(·) lies within the nonparametric confidence limits, as it should, but not conversely because the more efficient parametric analysis produces narrower confidence limits.

Extending the method to allow an adjustment for baseline covariates is straightforward. We simply add an intercept term to the right-hand side of (1.1), where xi is the vector of covariate values for subject i and include optimization with respect to ß as an additional step in each iteration of the back-fitting algorithm which, for ß, reduces to ordinary least squares with the current estimates of µ0(·) and {delta}(·) treated as offsets.

The method may seem somewhat elaborate, given its exploratory purpose. However, in practice it is considerably easier to implement than more ad hoc approaches, such as shifting individual subjects' sequences of intervention times to a common origin and smoothing by simple averaging in time-bands, because the fitting process is automated. Furthermore, the spline estimator with smoothing parameters selected by cross-validation avoids arbitrariness in the choice of time-bands within which to average, and adapts itself to variation in the intensity of observation times by applying more smoothing where the data are sparse in time. Indeed, it was our experience of encountering these difficulties when using an ad hoc method which led us to develop the approach described in the present paper.


   ACKNOWLEDGMENTS
 
This work was supported by the Portuguese Foundation of Science and Technology through the award of a studentship to Inês Sousa (SFRH/BD/10266/2002), and by the UK Engineering and Physical Sciences Research Council through the award of a Senior Fellowship to Peter Diggle (grant number GR/S48059/01). We also thank Eric Lim for making available to us the lung resection trial data.


   REFERENCES
TOP
 SUMMARY
 1. INTRODUCTION
 2. A NONPARAMETRIC ESTIMATOR...
 3. ANALYSIS OF THE...
 4. DISCUSSION
 REFERENCES
 

    DIGGLE, P. J., HEAGERTY, P., LIANG, K. Y. and ZEGER, S. L. (2002). Analysis of Longitudinal Data, 2nd edition. Oxford: Oxford University Press .

    GREEN, P. J. and SILVERMAN, B. W. (1994). Non Parametric Regression and Generalised Linear Models—A Roughness Penalty Approach. London: Chapman and Hall .

    HASTIE, T. J. and TIBSHIRANI, R. J. (1990). Generalized Additive Models. London: Chapman and Hall .

    LIM, E., ALI, A., CARTWRIGHT, N., SOUSA, I., CHETWYND, A., POLKEY, M., GEDDES, D., PEPPER, J., DIGGLE, P. and GOLDSTRAW, P. (2004). Effect and duration of lung volume reduction surgery on pulmonary function: mid term results of the Brompton trial. Preprint.

    MATTHEWS, J. N. S. (2000). An Introduction to Randomized Controlled Clinical Trials. London: Arnold.

    MOLENBERGHS, G., THIJS, H., JANSEN, I., BEUNCKENS, C., KENWARD, M., MALLINCKRODT, C. and CARROLL, R. (2004). Analyzing incomplete longitudinal clinical trial data. Biostatistics 5, 445–464.[Abstract]

    NELDER, J. A. and MEAD, R. (1965). A simplex algorithm for function minimization. Computer Journal 7, 308–313.

    POCOCK, S. (1983). Clinical Trials, A Practical Approach. Chichester: Wiley.

    RICE, J. A. and SILVERMAN, B. W. (1991). Estimating the mean and covariance structure nonparametrically when the data are curves. Journal of the Royal Statistical Society, Series B 53, 233–243.

    SENN, S. (2002). Cross-over Trials in Clinical Research, 2nd edition. Statistics in Practice. Chichester, UK: John Wiley & Sons.

    SMITH, D. M. and DIGGLE, P. J. (1999). Compliance in an anti-hypertension trial: a latent process model for binary longitudinal data. Statistics in Medicine 18, 357–370.

    Received October 25, 2004; revised February 25, 2005; accepted for publication February 25, 2005.


    Add to CiteULike CiteULike   Add to Connotea Connotea   Add to Del.icio.us Del.icio.us    What's this?



    This Article
    Right arrow Abstract Freely available
    Right arrow FREE Full Text (PDF) Freely available
    Right arrow All Versions of this Article:
    6/3/479    most recent
    kxi023v2
    kxi023v1
    Right arrow Alert me when this article is cited
    Right arrow Alert me if a correction is posted
    Services
    Right arrow Email this article to a friend
    Right arrow Similar articles in this journal
    Right arrow Similar articles in PubMed
    Right arrow Alert me to new issues of the journal
    Right arrow Add to My Personal Archive
    Right arrow Download to citation manager
    Right arrowRequest Permissions
    Right arrow Disclaimer
    Google Scholar
    Right arrow Articles by Sousa, I.
    Right arrow Articles by Diggle, P. J.
    Right arrow Search for Related Content
    PubMed
    Right arrow PubMed Citation
    Right arrow Articles by Sousa, I.
    Right arrow Articles by Diggle, P. J.
    Social Bookmarking
    Add to CiteULike   Add to Connotea   Add to Del.icio.us  
    What's this?