In this section, we illustrate the application of the methods, mentioned in the previous section, by using data from the individual-patient-data metaanalysis of patients with advanced colorectal cancer (see Section 2.2.4). The data were analyzed before by Burzykowski, Molenberghs, and Buyse (2004). The goal of the analysis is to evaluate tumor response as a surrogate for overall survival (OS). In particular, we consider two forms of the response: the four-category (complete response, CR; partial response, PR; stable disease, SD; progressive disease, PD) and the binary (complete/partial versus stable disease/progression) ones.

Using SAS

Four-Category Tumor Response

We first conduct the analysis of the four-category response using the Plackett- copula model with the marginal models (6.3) and (6.4); the latter is specified by assuming the Weibull hazard functions. Note that, as the response, we use the best response observed across the entire observation period for a patient.

We fit the model with the help of the SAS macro (see Section 12.5.4). For this model, the estimated value of the copula parameter в is equal to 6.52 (95% CI: [5.79, 7.25]). It indicates a substantial association between response and OS at the individual-patient level: the odds for surviving beyond time t given a higher (for instance, complete) response are about 6.5 times higher than the odds of surviving beyond time t given any lower response (partial response, stable disease, or progression). Note that the estimated value differs slightly from the value of 6.78 reported by Burzykowski, Molenberghs, and Buyse (2004). This is due to the fact that these authors used the marginal model (6.2) with some of the terms щ_{0},1,..., щ_{0},к_{-}1 constrained to 0 for trials with missing response categories.

It should be noted that в involves a comparison of survival times of patients classified according to tumor response. It is well known that such a comparison is likely to be length-biased, because a response to treatment is not observed instantaneously. As a result, patients who enjoy long survival times are more likely to be responders than non-responders, and therefore the survival of responders is likely to be biased upward compared with that of nonresponders. One method of correcting for length bias in such a comparison is a landmark analysis (Anderson et al., 1983). We do not attempt to conduct such an analysis here; it was considered by Burzykowski, Molenberghs, and Buyse (2004).

Figure 6.1 presents the estimated trial-specific treatment effects, i.e., logarithms of hazard ratios on OS and log-odds ratios for tumor response; each trial is represented by a circle of the size proportional to the trial sample size. In the dataset, higher response categories (PD, SD, PR, CR) were coded with higher scores (1, 2, 3, 4). Hence, according to the parameterization used in model (6.3), a is the logarithm of the ratio of the odds of observing a lower tumor-response category for the experimental treatment to the odds for the control treatment. On the other hand, following (6.4), Д is the logarithm of the ratio of the hazard of death for the experimental treatment to the hazard for the control treatment. Thus, Figure 6.1 indicates that the larger experimental- treatment-induced reduction of the odds of a “worse” response (OR<1), the larger the corresponding reduction in the hazard of death (HR<1). However,

FIGURE 6.1

Advanced Colorectal Cancer. Triad-level association between copula-model- based treatment effects on four-category tumor response and OS (both axes on a log scale). The circle surfaces are proportional to trial size.

the association between the effects appears weak. A simple linear regression model, fitted without any adjustment for the estimation error present in the estimated treatment effects, yields the following regression equation:

with the standard errors of the intercept and slope estimated to be equal to 0.044 and 0.043, respectively. The corresponding value of Rriai(r) is equal to 0.125 (95% CI: [-0.108, 0.358]); the negative lower limit of the CI is due to the use of the normal approximation to the distribution of (R?t_{rial}(_{r}))^{2} (see Section 5.3.1.1). The value is based on the following estimate of the variance- covariance matrix D of the (random) treatment effects:

for which the condition number (the ratio of the largest to the smallest eigenvalue) is equal to 21.4. This value is small enough to regard the obtained estimate as numerically stable.

After adjusting for the estimation error by using model (5.13)-(5.14) fitted with the PROC MIXED code described in Section 5.3.1, R2_{rial}(_{r}) is estimated to be equal to 0.979 (95% CI: [-10.3,12.3]). The value is based on the following estimate of the variance-covariance matrix D of the (random) treatment effects:

for which the condition number is equal to 9216.7. This value is very large and indicates that the obtained estimate of R_{r}i_{a}i(_{r}) is numerically unstable. This is also seen from the CI, which is very wide due to a huge estimate of the standard error of R_{r}i_{a}i(_{r}), equal to 5.77. Thus, the analysis after accounting for the fact that we analyze estimated, and not the true random, treatment effects, is actually non-informative.

As an alternative, one could consider an analysis, in which the estimated treatment effects would be appropriately weighted to reflect the precision of their estimation. An obvious issue is the choice of the weights. A simple choice is to use the sample size. However, the precision of the estimation of the treatment effect on OS depends rather on the number of deaths observed in the trial, not on the total sample size. On the other hand, the precision of the estimation of the treatment effect on tumor response depends on the observed numbers of different types of responses. Hence, it is not clear to what extent the analysis weighted by the trial sample size can correct for the estimation error and provide a correct estimate of R_{r}i_{a}i(_{r}).

If the weighted estimate of R_{rial}(_{r}) is desired, it can be obtained by using any procedure providing the coefficient of determination for a simple linear regression model, e.g., PROC REG. However, an important issue is the estimation of the precision of the obtained estimate. Toward this aim, PROC MIXED can be used. In particular, we need data in the “long” format, with two records per trial, one providing the estimate of effect for OS, and the other one for tumor response. An illustration of a few first records in a such dataset is:

trial effect ssize endp

1 -0.23164 306 MAIN

1 -0.50734 306 SURR

2 -0.27657 249 MAIN

2 -0.52413 249 SURR

3 0.03843 164 MAIN

3 -0.14893 164 SURR

Then we fit a (weighted) linear model with the help of the following SAS syntax:

PROC MIXED data=ests;

CLASS trial endp;

WEIGHT ssize;

MODEL effect=endp / noint;

REPEATED endp / subject=trial type=unr;

RUN;

Note the use of the WEIGHT statement, which applies the sample sizes, contained in the variable ssize, as weights. Moreover, the type=unr option in the REPEATED statement provides directly an estimate of the correlation coefficient between the treatment effects, i?triai(r), and its estimated variance. Subsequently, we can obtain the estimate of R2_{rial}(_{r}) as R?2_{rial}(_{r}), with an estimate of its variance computed by using the delta method, i.e.,

^{Var}(^{R}t^{2}rial(r)) ~ ^{4 X R}t^{2}rial(r) ^{X Var(} ^triagr) ^

The aforementioned approach is applied to the treatment-effect estimates obtained from the copula model by using a dedicated SAS macro (see Section 12.5.4). It yields an estimate of R2_{rial}(_{r}) equal to 0.144 (95% CI: [-0.106,0.394]). The underlying weighted linear regression model is

with the standard errors of the intercept and slope estimated to be equal to 0.039 and 0.048, respectively. This is the line indicated in Figure 6.1 as “predicted.” The “weighted” estimates are based on the following estimate of the variance-covariance matrix D of the (random) treatment effects:

for which the condition number is equal to 17.3.

The “weighted” results are quite similar to those obtained based on the unweighted, “naive” linear regression model (6.6). It is worth mentioning that the use of PROC REG yields the same estimates of R2_{rial}(_{r}) and linear-regression coefficients, but with somewhat larger standard errors for the intercept (0. 04) and slope (0.049) of equation (6.9). This is due to the fact that the PROC MIXED estimates are based on the (restricted) likelihood based on a normal distribution, whereas the PROC REG estimates are obtained by using the weighted least-squares method.

If we focus on the trial-level surrogacy, we can conduct an analysis based solely on the marginal models (6.2) and (6.4), with the latter specified without assuming any particular form of the baseline hazards. Toward this aim, model (6.2) can be fitted by using PROC LOGISTIC and (6.4) by using PROC PHREG. In particular, the call to PROC LOGISTIC can be as follows:

PROC LOGISTIC data=dataset;

CLASS trial treat(ref=first) / param=glm;

MODEL response=trial treat(trial);

RUN;

Variable trial is the trial identifier, treat is the treatment indicator (with 0 for the control group and 1 for the experimental one), and response is the response. Of interest are the trial-specific treatment-effect estimates, provided by the nested effect treat(trial).

On the other hand, (6.4) can be fitted by using PROC PHREG, with a call as follows:

PROC PHREG data=dataset;

CLASS trial treat(ref=first) / param=glm; MODEL surv*survind(0)=treat(trial);

STRATA trial;

RUN;

Variable surv is the observed survival time, while survind is the event indicator (with 0 for censored times and 1 for event times).

After fitting the models, the variance-covariance of the estimated coefficients of the models can be estimated while taking into account the association between the response and OS. Toward this aim, a dedicated SAS macro can be used (see Section 12.5.4).

Figure 6.2 presents the estimated log-odds ratios and log-hazard ratios. As was the case in the copula-based estimates (see Figure 6.1), the association between the treatment effects is only moderate. A simple linear regression model, fitted without any adjustment for the estimation error present in the estimated treatment effects, yields the value of Rj-iaii» equal to 0.145 (95% CI: [-0.105, 0.395]), somewhat smaller than the estimate obtained from the copula model with Weibull marginal hazards. The value is based on the following estimate of matrix D:

for which the condition number is equal to 25.3. The resulting regression equation is

ln(HR_{OS}) = -0.023 + 0.082 x ln(OR_{PFS}), (6.12)

with the standard errors of the intercept and slope estimated to be equal to 0.052 and 0.144, respectively. The regression equation is very similar to (6.6) obtained for the model based on the Plackett copula.

Analysis adjusted for the estimation error leads to an estimate of variance- covariance matrix D which is singular. The analysis weighted by the trial sample size, conducted with the help of the dedicated SAS macro (see Section 12.5), yields the following estimate of the variance-covariance matrix D of the (random) treatment effects:

for which the condition number is equal to 22.0. The corresponding estimate of

FIGURE 6.2

Advanced Colorectal Cancer. Triad-level association between the marginal- model-based treatment effects on OS and four-category tumor response (both axes on a log scale). The circle surfaces are proportional to trial size.

Rtriai(r) i^{s е}Ч^{иа}1 to 0.144 (95% CI: [-0.111,0.399]), with the linear regression equation

and standard errors of the intercept and slope estimated equal to 0.036 and 0.043, respectively. The line corresponding to equation (6.14) is indicated in Figure 6.2 as “predicted.”