On model pre-specification in confirmatory randomized studies.

David Edwards

June 10, 1998

Statistics Department, Novo Nordisk A/S, DK-2880 Bagsvaerd, Denmark. tlf: +45 44 42 61 35 fax: +45 44 42 10 65 e-mail: DEd@novo.dk

On model pre-specification in confirmatory randomized studies.

Summary

Typically, the primary purpose of confirmatory randomized trials, such as drug trials sponsored by the pharmaceutical industry, is to determine whether there is a treatment effect, and if there is, to estimate the size of the effect. For such studies it is accepted practice to pre-specify the statistical model to be used in the primary analysis. The reason for this is a concern that if the model were to be chosen on the basis of the data, the model most favourable to the sponsor might be chosen, with consequent inflation of the type I error. The purpose of this article is to show that, in a sense, this concern is needless. It is shown that if the model is chosen in a blinded fashion and randomisation-based tests for no treatment effect are used, then the type I error is controlled. A similar technique to derive unbiased estimates of treatment effect is also described. This approach may be of value when there is uncertainty as to the correct model when the study is being planned.

1 Introduction

In settings where confirmatory inference is required, it is sometimes assumed that the method of statistical analysis must be fully specified in advance. This is the case, for example, with comparative efficacy trials sponsored by the pharmaceutical industry, whose primary purpose is to determine whether a drug is efficacious and if it is, to estimate the size of the effect. Due to the competitive and expensive nature of drug development, the company sponsoring the trial and analyzing the data often has strong vested interests in the outcome. If the method of analysis were selected on the basis of the data, there would be a risk that the most favourable method would be selected, thereby giving rise to biased conclusions. For this reason, regulatory guidelines^{1,2,3} stipulate that the trial protocol must contain a full description of the method of analysis to be used, including^{1,3} a specification of the statistical model. (These requirements are not as rigid as they might seem, since the guidelines^{1,3} also allow the planned analysis to be modified prior to unblinding: see Section 6 below). Consequently it is currently standard practice for such trials to specify a model in the study protocol and use this in the analysis.

Clearly, when the appropriate model can be known in advance, this practice is unproblematic. Usually results from similar trials are available, and will provide a good basis for model choice. However, sometimes there may be uncertainty about the model. This can be due to a variety of factors: for example, although similar trials may be have performed, their results may not be available, or the trials may differ in ways that necessitate modifications to the analysis.

When there is uncertainty about which model to use, choosing one in advance may lead to suboptimal inference. In principle, this could lead to the unnecessary exposure of patients to suboptimal treatments, and so could have ethical and economic ramifications. This presents a dilemma: on the one hand, it is desirable that the model be pre-specified to prevent bias, but on the other hand this may lead to suboptimal inference.

To resolve this dilemma, I propose the following procedure. Instead of pre-specifying a particular model, a family of possible models and a model selection criterion are pre-specified. These models describe the distribution of the response, possibly incorporating some covariates, but without including a treatment term. When the study is complete, but before the blind is broken, a model is selected using the criterion. After the blind is broken, a randomisation test of no treatment effect is made using test statistics appropriate to the selected model. A similar technique is used to derive point and interval estimates of the size of the effect. I show in Section 2 that this procedure is valid, in the sense that the type I error does not exceed its nominal level. In Section 3, some related issues are discussed. Section 4 addresses computational aspects, and Section 5 presents two examples. Section 6 discusses some operational aspects particular to drug trials, and the final section relates the results to recent work on model uncertainty.

2 Main Results

Consider a randomised study with N patients, and let the covariates (prerandomisation measurements) be denoted by an N × q matrix X, the randomised treatment allocation as Z, and the response as Y. Let Ω be the set of possible treatment allocations. I follow the usual convention of distinguishing random variables from their realizations by using upper case for the former and lower case for the latter.

The hypothesis of no treatment effect, denoted H, states that the response of the i-th unit, Y_{i}, would have been the same whatever treatment had been allocated to this unit (and whatever treatments had been allocated to the other units).

Consider a family of parametric models F for the distribution of Y given X. Note that they do not include any treatment effect, and so are consistent with H. Suppose that each model m in F can be augmented to a model m′ by adding a treatment term, and denote the model family thereby generated as F′. For example, F could be the family of linear regression models of the form y_{i} = b_{0} + ∑_{j∈S_{m}} b_{j}x_{ij} + ϵ, where the value of the j-th covariate for the i-th observation is written x_{ij}, S_{m} is some subset of the set {1, … q}, b_{j} (j = 0, … q) are unknown parameters, and ϵ is an error term, normally distributed with zero mean and unknown variance. The corresponding model in F′ would add a treatment term.

For a given model m in F, consider a test for m versus m′, using a test statistic t(m, x, y, z). The randomisation test using this statistic is based on the distribution of t(m, x, y, Z): that is, X and Y are considered to be fixed, and Z random. Under H, for a level α test,

Pr{t(m, x, y, Z) > t^{α}(m, x, y)} ≤ α,

where t^{α}(m, x, y) is the 1 − α point of the distribution of t(m, x, y, Z). This can be re-written as

∑_{z} Pr(Z = z)I{t(m, x, y, z) > t^{α}(m, x, y)} ≤ α,  (1)

where I() is the indicator function and the summation is over all z in Ω. See, for example, Rosenbaum^{4}, Section 2.4.1.

Similarly, the attained significance level of the observed statistic t^{∗} is

Pr{t(m, x, y, Z) > t^{∗}} = ∑_{z} Pr(Z = z)I{t(m, x, y, z) > t^{∗}}.

If this is small, either H does not hold, or an unlikely event has occurred.

Suppose now that the model is not pre-specified, but is stochastic. Then the type I error is

Pr{t(M, x, y, Z) > t^{α}(M, x, y)}

where this is evaluated under H. This can be written as

∑_{z} ∑_{m} Pr(M = m, Z = z)I{t(m, x, y, z) > t^{α}(m, x, y)}  (2)

where the summation is over all m in F.

If M and Z are independent under H then it follows directly from (1) that the type I error is at most α. To satisfy this condition, a model selection algorithm that uses only the response Y and the covariates X may be employed. Write this as m = msa(x, y). To see that the condition is fulfilled, note that Z is independent of Y under H, and that since Z is random, Z is independent of X. So Z is independent of (X, Y) under H and the result follows.

In principle model selection need not be performed by an automatic algorithm, but could incorporate more complex, subjective or ad hoc judgements by the analyst. But it would be difficult to ensure that M and Z are independent under H, since this would require that the analyst does not use any post-treatment information (apart from Y) that could concievably be related to treatment. That is, information of events that are, as it were, ”causally downstream” from treatment. If the analyst is aware of such events, he or she may not be able to suppress this knowledge. (To see this, try to run around the house three times without thinking of the word ”gnu”). So it would be difficult to fulfill this requirement rigorously, except by using a pre-specified algorithm.

Point and interval estimates closely related to the test just described can also be derived. First I describe these briefly for the case when the model is given, and then discuss the effect of letting the model be stochastic.

Suppose for the sake of simplicity that there are two treatment groups, so that Z is coded 0 for the first group and 1 for the second, and that Y is a metric variable, so an additive model for Y is appropriate. Consider the hypothesis H_{τ} that the effect of the second treatment is additive and of magnitude τ. In other words that the response of the i-th unit, Y_{i}, would have been τ larger if the second treatment had been allocated to this unit than if the first treatment had been allocated.

Consider testing the null hypothesis H_{τ}. This is equivalent to testing that the adjusted response Ỹ = Y − τZ satisfies the hypothesis of no treatment effect. Thus given a pre-specified model m in F, and hence a test statistic, a randomisation test of H_{τ} can be performed. Under H_{τ}, when the true value (τ) is tested, it is rejected with probability no greater than α. So under H_{τ},

Pr{t(m, x, ỹ, Z) > t^{α}(m, x, ỹ)} ≤ α.

This can be re-written as

∑_{z} Pr(Z = z)I{t(m, x, ỹ, z) > t^{α}(m, x, ỹ)} ≤ α.  (3)

This test is inverted in the usual way to derive a 1 − α confidence region for τ, this being defined as the set of τ for which the test does not lead to rejection at level α. It can be shown^{4} that if the test is order-preserving, then the confidence region takes the form of an interval, and line-searching methods can be used to find the endpoints.

Now consider point estimation of τ. Often an estimator τ̃ = τ̃(m, x, y, z) derived from model-based theory is available and has good properties in the randomisation framework. Suppose that this is unbiased in the randomisation distribution when m is pre-specified, that is,

∑_{z} Pr(Z = z)τ̃(m, x, y, z) = τ

under H_{τ}. If the model is stochastic, then the expectation becomes

∑_{z} ∑_{m} Pr(M = m, Z = z)τ̃(m, x, y, z)

where this is evaluated under H_{τ}. If M and Z are independent under H_{τ}, then the estimator is unbiased. If the model is selected using m = msa(x, y), then M and Z are independent under H, but not necessarily under H_{τ}, so the estimator may be biased. But if the model is selected using the adjusted reponses, that is, m = msa(x, y − τz), then M and Z are independent under H_{τ}. It follows that the bias can be removed if a pair (m̂, τ̂) can be found that jointly satisfy

τ̂ = τ̃(m̂, x, y, z)

and

m̂ = msa(x, y − τ̂z).

An obvious approach is to use an iterative process, the ith step of which calculates τ̂_{i} = τ̃(m̂_{i−1}, x, y, z) and then m̂_{i} = msa(x, y − τ̂_{i}z). A convenient starting point would set m̂_{0} = msa(x, y), the model used in the test of H. The iteration converges when m̂_{i} = m̂_{i−1}.

This method of bias correction appears to be novel. It would be useful to clarify when such an estimate (m̂, τ̂) can be found, but this is not attempted here. Note that even if τ̃ is biased when m is pre-specified, the method can still be used to eliminate model selection bias. Note also that when Y is nonmetric it is not appropriate to adjust the responses directly, but the method may still be applicable. For example, for logistic regression models the step m = msa(x, y − τz) would be performed by modifying the models in F to include the term τz as an offset, that is to say, as a variable in the linear predictor whose coefficient is unity. This is illustrated in Section 5.2 below.

Calculation of interval estimates when the model is stochastic appears difficult. Consider first testing the point hypothesis H_{τ}. To extend equation (3) to the stochastic setting requires that the model be chosen using the adjusted reponses m = msa(x, y−τz). Inverting this test to derive a confidence region would mean that a model should be selected for each τ tested. This would require much computation, and from an interpretative point of view it would be strange that the model could vary over the region. It seems preferable to derive a confidence interval conditional on a specific model, for example the model chosen in the test for no treatment effect, or if this test is rejected, the model chosen when estimating τ.

In summary, the following analysis procedure suggests itself:

1. Before breaking the blind, select a model m_{0} = msa(x, y).

2. After breaking the blind, test H using a test statistic based on m_{0}. It has been shown that the type I error is controlled.

3. If H is rejected at step 2, estimate τ. A technique to perform unbiased estimation has been described: this entails selection of an appropriate model, say m_{1}.

4. Calculate a confidence interval for τ, based on m_{0} if H was accepted, otherwise m_{1}.

3 Some related issues

In this section some issues related to the procedure described are discussed. After some general remarks about randomisation inference, two questions are addressed. Firstly, would the procedure be valid if model-based, rather than randomisation-based inference were used? Secondly, it has been shown that the procedure is valid, but is it useful, and if so, when?

Randomisation-based inference was introduced by Fisher^{5} and its merits and demerits have been discussed by several authors^{6,7}. Presumably a consensus view is that the randomisation framework is elegant and powerful, but is of somewhat limited utility, and may be computationally prohibitive. It is most useful for making cautious comparisons between randomised groups; but that is precisely the present objective. Computational aspects are described briefly below. A general limitation is that the space of possible randomisations (Ω) must be sufficiently rich — if it is not, then conditioning on Y may be too severe.

One point at issue is whether the null hypothesis H is too strong, since it stipulates that there should be no effect of treatment for any patient in the study. For some purposes, for example in a safety context, this seems reasonable enough: for example, we would say that a drug causes a severe adverse event even if it only does so in a single patient. But in other contexts we might be more interested in the weaker null hypothesis that the average effect is zero^{8}. Under this weaker null hypothesis the size of a randomisation test may exceed its nominal level: however, Romano^{9} showed that in the two-sample case, when there is equal allocation to treatment groups, randomisation tests have asymptotic validity (i.e., the size tends to the nominal level) under the weaker null hypothesis. It is easy to see that the randomisation tests described in the previous section inherit this property.

It is instructive to compare model-based and randomisation-based tests, for a given m and test statistic t(m, x, y, z). The former are based on the distribution of t(m, x, Y, z), and the latter on the distribution of t(m, x, y, Z). That is, model-based tests are based on a set of hypothetical repetitions of the experiment in which the covariates and the treatment allocation are held fixed, and the outcome measure Y varies according to the models posited. They are constructed such that under m the type I error is controlled, and under the alternative m′ the power is in some sense maximized. In contrast, randomisation tests are based on hypothetical repetitions of the experiment in which the covariates and the response are held fixed, and the treatment allocation varies. Proponents of randomisation tests emphasize that since randomisation is under the control of the experimenter, the possible allocations z and their probabilities are known: whereas model-based tests ”depend on unverifiable assumptions about distributions, independence and the like”^{7}.

Hoeffding^{10} studied the asymptotic power of randomisation tests and found that for many tests, randomisation tests have the same asymptotic power as the corresponding model-based tests. Thus when m′ is true, there is generally no loss of asymptotic efficiency in using the randomisation-based test for m versus m′ compared with the model-based test.

The key difference appears when m′ does not hold, since the size of the model-based test is only controlled under m, whereas the size of the randomisation test is retained under the more general hypothesis H. The behaviour of model-based tests when the model is misspecified has been much studied for particular model families, but only a few general results are available^{11}. There has recently been a focus in several areas on studying the effect of omitting covariates: see, for example, Keiding, Andersen and Klein^{12} in the context of survival analysis. Gail, Tan and Piantadosi^{13} study tests for no treatment effect in randomised trials. In particular, they examine the behaviour of score tests under various generalized linear models and survival models when one relevant covariate is omitted from the regression. They show, using simulation, that in many cases supranominal size is obtained: for example, for a poisson regression model, a test with a nominal size of 0.05 had an actual size of 0.23 (op. cit., Table 1). However, some models, including the normal linear model and the logistic regression model, are robust to omitted covariates. Note that pre-specifying the model can be regarded as a special case of the model selection algorithm of Section 2, where msa(x, y) = m, ∀x, y. As the poisson regression example shows, the modelbased test does not always retain its size when the model is false. It follows that the procedure cannot in general be valid when model-based tests are used, even approximately. But it may be true that for robust model families (coupled, perhaps, with specific model selection algorithms) the procedure could be valid in some approximate sense: this would be an interesting research topic but is not pursued further here.

When the model is not fixed beforehand, but is chosen from the data using m = msa(x, y), then the overall size of the randomisation test of H is retained. Ensuring that model choice is independent of treatment allocation under the null hypothesis prevents that the most favourable model be chosen, for example a model corresponding to a small p-value in a test for no treatment effect. Similarly, ensuring the independence of model choice and treatment allocation under an alternative H_{τ} prevents bias in the estimation of τ: this lies behind the method of bias correction described above. In both cases, the purpose of selecting the model on the basis of the data is to adopt one that approximates the data well, in order to capture the efficiency of the model-based inference that would obtain if the model were true.

Clearly, the choice of model selection algorithm is critical. Many different algorithms have been proposed to choose the ”best” model on the basis of the data, and it is not the purpose of this paper to recommend any particular one. However, a general approach that appears appropriate is to select the model that minimizes either the Akaike information criterion^{14} (AIC), given as −2 ln(Q) − 2p, or the Bayesian information criterion^{15} (BIC), given as −2 ln(Q) − ln(N)p. Here Q is the likelihood, N the number of observations and p the number of free parameters in the model. The AIC is derived from information-theoretic arguments, whereas the justification for using the BIC is that it corresponds, asymptotically, to choosing the model with the largest posterior probability.

There is no guarantee of a gain in efficiency for the proposed procedure. A typical situation where a loss of efficiency will occur is the following. Suppose that, by chance, a covariate thought to be potentially prognostic is confounded with the treatment allocation. Then if there is a treatment effect, the covariate will be predictive of the response, and so will tend to be included in the model. This will lead to loss of power for the test of no treatment effect. But from another point of view, the alleged treatment effect could in fact be due to the covariate, so it is not clear what conclusion should be drawn: the study seems inherently flawed.

A related issue concerns the use of model selection algorithms. It is well known (for an overview, see Chatfield^{16}) that such procedures, when used with linear normal models, have a tendency to select a model with a smaller residual variance than the true model: in other words, the residual variance is under-estimated. This appears to contradict the statement that the procedure controls the type I error, since if the residual variance is underestimated, the variance of the estimated effect will also be under-estimated. So it would seem that a test for no treatment effect based on this model will have inflated type I error. However, this argument presupposes that the magnitude of the estimated effect remains constant, which is not the case: it is also reduced. It is important to distinguish between model search among models without a treatment term, and among models including a treatment term. The latter is not ”safe” with respect to testing for no treatment effect, whereas the former is, in the randomisation framework at least.

A general point is that it appears that different models are required for different purposes. A model chosen to perform a valid test of H weighs the evidence in favour of H and thus will tend to give a biased estimate of τ: conversely, a model chosen to give an unbiased estimate of τ will tend to have supranominal size if used to test H.

The procedure rests on an assumption that the treatment effect is small. For example, when Y is metric and Z binary, we may assume models for Y given (X = x, Z = z) of the form f(x) + τz + ϵ, where f(x) is a systematic component and ϵ a random component. The distribution of Y given X = x is then a mixture of these for z = 0 and 1. If τ is sufficiently small, then this mixture will be close enough to the conditional distributions to be informative about the form of f(x) and/or the distribution of ϵ. This assumption is part of the rationale of clinical trials. As Piantadosi^{17} remarks,

... clinical trials are suited to situations where it is important to learn about treatment effects or differences that are about the same magnitude as the random error and bias that invariably characterize medical studies.

A common situation where the approach is not likely to be useful can be described as follows. Suppose that the model family F contains a full model, i.e. a model containing all the other models in the family as submodels, say m_{f}. Suppose also that the number of free parameters in m_{f} is small, compared to the number of patients in the study. If some m′_{0} in F′ is a good approximation to the data, then m′_{f} will also be a good approximation, and the loss of efficiency in using m′_{f} instead of m′_{0} will be minimal. There will be little to gain in selecting a model from F from the data. If we are willing to assume that m′_{f} is true, if then model-based inference using m′_{f} will be optimal or near-optimal. (Clearly, if the model family is robust to omitted covariates, then this assumption is less perilous.) Otherwise, randomisationbased inference using m′_{f} can be adopted.

4 Computational Aspects

In principle, it is straightforward to perform randomisation inference. The study design determines the set of possible treatment allocations Ω, and for each z in Ω, Pr(Z = z). By calculating t(m, x, y, z) for each z in Ω, the distribution of t(m, x, y, Z) under H can be generated, and thereby, for example, the critical point t^{α}(m, x, y). However, in practice this is often infeasible. For moderately sized studies, the number of allocations can be astronomical. Moreover, for models requiring iterative parameter estimation, a separate iterative loop may be needed to calculate each t(m, x, y, z). A powerful general approach is to use Monte Carlo sampling: in other words, to sample a number of random z’s from Ω, calculate t(m, x, y, z) for each, and estimate the p-value as the proportion exceeding t(m, x, y, z_{obs}). This can compute significance levels to any desired accuracy.

An overview of various approaches, and a good bibliography, to these computational problems can be found in Good^{18}. See also Tritchler^{19} and Garthwaite^{20} for algorithms for the calculation of confidence intervals.

Gail et al.^{13} propose a simple method of performing randomisation tests for generalized linear models and certain survival models, that is closely tied to the current framework. This consists of calculating residuals under a model without a treatment term (i.e. an m in F), and then comparing the group means of these residuals using randomisation tests. They show that under m′ this test has the same efficiency as the score test for m versus m′. This reduces the task to that of performing a randomisation t-test (in the two-sample case), or randomisation F-test (when there are more than two treatments). They also derive an expression for the asymptotic variance of the test statistic, enabling a simple test based on a standardized deviate to be constructed.

In the examples below, the program MIM (Edwards^{21}) was used to perform exact randomisation t-tests, to calculate the associated confidence intervals, and to perform the model selection.

5 Two examples

To illustrate the approach I consider two simulated examples. Table 1 displays a set of data consisting of 20 observations with 20 variables. The first 17 variables, x_{1}, … x_{17}, have been taken from a set of data whose origin is not germane to the present topic. These will be used as covariates. The variable Z is a binary factor that has been generated using a pseudo-random number generator. In the examples it is assumed that the numbers in each group are fixed by design, so that Ω consists of all treatment allocations with 10 in each group. The two remaining variables, y_{1} and y_{2}, are response variables, and are described in Sections 5.1 and 5.2 respectively.

5.1 Example 1: multiple regression

This example is constructed to illustrate a situation in which a large number of potentially prognostic variables are available. In Table 1, y_{1} has been generated using the expression

y_{1} = 0.5 + 0.2x_{1} + 0.4x_{3} + 0.9z + ϵ  (4)

where ϵ ∼ N(0, 0.5). Suppose that the data stem from a study in which y_{1} is the primary outcome, and whose purpose is to determine whether or not the randomised treatment Z has an effect on this outcome. Suppose also that the 17 covariates are considered potential prognostic factors.

 x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 Z   y1      y2
 4  1  1  3  2  1  1  5  4   2   5   4   1   6   4   1   5  0  2.303 26.3818
 3  3  6  3  3  3  4  6  5   5   5   4   3   5   2   1   6  1  4.877 15.1875
 2  1  3  3  2  1  3  6  4   2   6   4   1   6   3   6   6  1  3.161  9.6196
 2  2  6  2  2  2  4  5  5   4   1   3   3   5   2   1   6  1  4.602 12.4505
 4  4  4  4  4  3  4  5  4   3   6   3   2   4   2   1   5  0  3.752 30.2431
 2  1  4  2  2  2  1  6  4   2   5   3   1   6   3   1   6  1  3.296  9.2725
 2  2  1  4  2  3  2  5  4   3   6   5   1   5   5   5   3  1  1.815 11.1357
 3  3  4  5  2  2  4  4  4   3   3   4   2   5   5   1   4  1  3.758 16.9926
 3  2  4  5  3  3  3  5  4   3   6   3   2   5   4   2   4  0  2.822 11.7300
 3  3  3  4  3  3  4  3  4   3   4   5   3   5   3   4   4  0  2.682 18.2954
 3  2  5  5  4  5  2  5  5   2   5   5   2   5   2   1   4  0  2.783 14.9717
 4  6  6  5  3  2  5  4  6   3   6   5   3   6   2   1   4  0  4.271 28.1513
 4  1  4  4  3  2  2  4  4   3   2   4   2   5   3   1   4  1  3.732 27.6856
 2  2  5  5  2  2  2  3  4   2   4   5   3   5   2   1   5  0  2.710  6.5885
 3  3  6  4  3  3  5  5  4   4   5   4   3   5   4   1   5  0  3.531 14.8071
 3  1  4  5  4  2  2  6  4   2   5   3   1   6   5   1   5  1  4.295 16.3655
 3  2  4  3  3  3  4  5  6   4   6   5   3   5   2   5   4  1  3.621 16.7027
 4  2  3  4  4  4  2  4  4   3   4   2   3   5   5   4   5  0  2.975 24.4056
 3  2  4  4  2  3  4  5  4   3   6   5   2   5   6   1   6  1  3.668 15.4263
 4  2  2  5  2  2  1  6  5   1   6   5   2   6   4   5   6  0  2.317 24.9944

Table 1: Data used in Examples 1 and 2.

The model family F consists of standard multiple regression models of the form

y_{i} = b_{0} + ∑_{j∈S_{m}} b_{j}x_{ij} + ϵ_{i}  (5)

where S_{m} is some subset of the set {1, … 17}, and the ϵ_{i} are normally distributed with unknown variance.

A model selection algorithm that attempts to find the model with minimum BIC value is used^{21}. This starts with the null model, and successively adds terms, choosing the term which gives the greatest reduction in BIC. It stops when no further terms reduce the BIC. (This does not guarantee to find the model with minimum BIC.)

Table 2 compares the results of model-based inference for three models: the null model, the full model and the true model. The tests for no treatment effect based on the null and full models are marginally non-significant at the conventional 5% level. The confidence interval for the treatment effect under the full model is very wide, due to over-fitting (in all 19 degrees of freedom, but 18 parameters).

In the first step of the procedure — choice of model in F on which to

Model  p  τ̂ S.E.  95% C.I.
null .0633 .668 0.337 [−0.041, 1.377]
full .0870 .968 0.133 [−3.319, 5.255]
true .0000 .944 0.159 [0.607, 1.281]

Table 2: Comparison of three pre-specified models in Example 1.

base the test of no treatment effect — the algorithm selected the model m_{0} = {3, 5, 6, 10, 17}.

The test of the null hypothesis was carried out using the method proposed by Gail et al.^{13}, namely by calculating the residuals under m_{0} and comparing these between treatment groups using a randomised t-test. This gave p = .0025. This tallies closely with the corresponding model-based test for m_{0} versus m′_{0}, for which p = 0.0021.

Table 3 shows the progress of the bias correction process. Starting with m_{0}, the process converged after four iterations. The selected model was m_{1} = {1, 2, 3, 12, 17}, corresponding to an estimate of τ̂ = 0.9768. From Section 2 it is known that this estimate is unbiased. Table 4 shows estimates and both model-based and randomisation confidence intervals using m′_{0} and m′_{1}. It can be seen that m′_{0} under-estimates the true effect (0.9), and confidence intervals do not cover the true effect. Under m′_{1} a better estimate and confidence interval are calculated.

In summary, the analysis rejected the hypothesis of no treatment effect very strongly (p = 0.0021). The treatment effect was estimated as τ̂ =

\begin{matrix}Model & τ̂ \\ m_{0} = \{3, 5, 6, 10, 17\} & 0.5474 \\ \{1, 3, 6, 7, 8, 12\} & 0.7736 \\ \{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17\} & 0.9681 \\ \{1, 2, 3, 12, 17\} & 0.9768 \\ m_{1} = \{1, 2, 3, 12, 17\} & 0.9768 \end{matrix}

Table 3: Results of bias correction process in Example 1.

Model  τ̂ S.E. 95% C.I.  95% C.I.
m′_{0} 0.5474 0.143 [0.238, 0.856] [0.239, 0.859]
m′_{1} 0.9768 0.127 [0.701, 1.252] [0.694, 1.262]

Table 4: Comparison of the two selected models. The first confidence interval is model-based, the second is randomisation-based.

0.9768, with a 95% confidence interval of [0.694, 1.262]. If the procedure had used model-based rather than randomisation-based inference, then the conclusions would have remained essentially unchanged.

Some more general comments on the inclusion of covariates are in place here. It is quite common for many potential covariates to be recorded in clinical trials. These include the obvious demographic variables, information on previous medical history plus batteries of biochemical tests taken prior to randomisation. A sceptical attitude to these is appropriate. As Piantadosi^{17} remarks: ”Many times prognostic factors are simply convenient to use rather than having been collected because of a priori interest or biological plausibility. Such factors should probably be discarded from further consideration.” As seen in Section 3, the uncritical inclusion of covariates runs the risk of obscuring treatment effect through confounding with treatment allocation. See also Senn^{22} for a helpful discussion.

It is interesting to note in passing that a related method was proposed by Tukey^{7,23} for accounting for many covariates in clinical trials. This used composite scores constructed by regressing the response on each covariate, using the data from the combined treatment arms, i.e. margining over treatment. The score was then used as a regressor in testing and estimating treatment effects. This is related to the approach of Gail et al.^{13} described above in Section 4.

5.2 Example 2: Additive or Multiplicative Effect.

This example is contructed to illustrate use of the procedure when there is uncertainty as to whether a treatment effect is additive or multiplicative. In Table 1, y_{2} has been generated using the expression

y_{2} = exp(1.2 + 0.5x_{1} + 0.009x_{2} + 0.08z) + ϵ  (6)

f() covariates BIC
exp  x_{1}, x_{2} 76.43
exp  x_{1} 81.72
exp  x_{2} 792.26
identity  x_{1}, x_{2} 115.29
identity  x_{1} 117.36
identity  x_{2} 810.51

Table 5: Model choice for Example 2.

where ϵ ∼ N(0, 2.0). Suppose these data stem from a randomised confirmatory study to examine the effect on the response. It is believed that x_{1} and x_{2} are predictive of the response, but it is unclear whether they (and the treatment) act additively or multiplicatively. Models of the form

y_{2} = f(b_{0} + b_{1}x_{1} + b_{2}x_{2} + τz) + ϵ  (7)

are considered, where f() can be either the exponential function or the identity, and either b_{1} or b_{2} can be zero. The error term ϵ is normally distributed with unknown variance. This specifies the model family F′, the family F being defined by setting τ = 0. Both belong to the family of generalized linear models. The minimum BIC criterion is used to select models from F.

Two models that might have been prespecified are the additive and multiplicative models including both covariates. The model-based test for no treatment effect based on the former gave p = .1042, while the latter gave p = .0410, so the study conclusions are sensitive to the choice of model.

Applying the procedure of Section 2, the first step is to choose a model in F. Table 5 displays the results.

The multiplicative model with both covariates is selected. The randomisation test of no treatment effect based on the residuals from this model, gave p = .0458, consistent with the model-based test.

The estimate of treatment effect under the model was τ̂ = 0.1038, with an estimated standard error of 0.0480. The question of whether this estimate is subject to model selection bias involves a subtlety that is of interest in its own right. Observe that for the multiplicative models, the additive alternative H_{τ} is not applicable. Instead, a multiplicative hypothesis must be used: that the effect of the second treatment is multiplicative and of magnitude e^{τ}. Call this H_{τ}^{∗}. Then for the estimate to be unbiased it is required that its expectation is τ under H_{τ}^{∗}. To examine this, the model selection process is repeated, this time including an offset term 0.1038z in each model, and only considering the multiplicative models. The selected model still has the smallest BIC. So the conclusion is that given that the model is multiplicative, the estimate is not subject to model selection bias.

6 Operational Issues

In this section I discuss issues arising when considering how to operationalize the approach in the context of drug trials.

Power calculations are made difficult when a model is selected from the data: it would be natural to take a conservative approach based on a default model.

It has been assumed throughout that the null hypothesis is of no treatment difference, so that the burden of proof lies in demonstrating that there is a difference. Equivalence studies^{1,22}, however, are designed to demonstrate the equivalence of treatments, so for these the burden of proof is reversed. The null hypothesis is inequivalence, namely that |τ| > τ_{0}, for some specific τ_{0}. Equivalence is inferred when the two-sided confidence interval for τ falls within [−τ_{0}, τ_{0}]. Operationally this is equivalent to using two simultaneous one-sided tests, namely (i) of τ = −τ_{0} versus τ > −τ_{0} and (ii) of τ = τ_{0} versus τ < τ_{0}. These tests can be performed using the present approach, and any gain in efficiency obtained thereby would be in the sponsors interest, since it would increase power to reject inequivalence.

Some regulatory guidelines^{1,3} allow for a blind review of the data after the study has been conducted. Since this review is blind, it is taken to be unbiased in a broad sense. That is, decisions concerning, for example, classification of protocol deviations cannot be intentionally biased. In general, however, the lack of intentional bias does not guarantee that decisions taken do not in fact bias treatment comparisons. The results given in Section 2 can be regarded as clarifying how the primary analysis may be modified in the blind review, without biasing the test for no treatment effect in the strict sense. It should be noted that the term ”blind review” is sometimes taken to mean that the randomisation codes (say, A and B) for each patient are known to the statistician, but not their meaning - for example, that A is placebo and B active. Blinding in this sense is clearly inadequate for the present purpose.

In Section 2 it was emphasized that the selection process should not involve any post-treatment information that could concievably be related to treatment. It is useful to consider more closely which data-analytic actions are legitimate when selecting a model. It would be invalid to modify or redefine the outcome measure Y (say, by categorizing a continuous variable), or to exclude patients from the analysis, since these actions would redefine the primary comparison. Needless to say, exclusion of patients would also violate the intention-to-treat principle. Redefining covariates, and modelling variance heterogeneity by means of weights, are legitimate if they help attribute variation in Y to factors not related to treatment. There may, however, be grey areas. For example, there may be a covariate which explains much variation in Y, but which for some obscure reason is unrecorded for some patients. In this case one approach could be to impute the missing values, so as to avoid excluding patients: presumably this would not bias the treatment comparison, though this would seem difficult to address formally.

It is helpful to cite a recent guideline^{1} on this topic:

This pre-analysis review, blinded to treatment, should cover decisions concerning, for example, the exclusion of subjects or data from the analysis sets; possible transformations may be checked, and outliers defined; important covariates identified in other recent research may be added to the model; the use of parametric or non-parametric methods may be reconsidered. ... When the blinding is compromised by the possibility that treatment induced effects may be apparent in the data, special care will be needed for the blind review.

In a similar vein, Lewis^{24} discusses issues related to the pre-specification of analysis. He emphasizes its benefits, including ”it avoids accusations of bias arising from the selection of the most favourable analysis”, and goes on to state that

There is certainly no need for the statistician to feel restricted during analysis. If the state of knowledge has changed since the protocol was written in a way which clearly affects the analysis, then it is appropriate for the analysis plan to be changed. A written explanation should then be provided; if the reason for the change is sound, this is generally easy to write and easy for others to accept. Sometimes preliminary examination of the data leads to modified methods of analysis. For example, earlier distributional assumptions may be found to be wrong (preferably during a preliminary blinded analysis of the data), and this provides an acceptable justification for change which can be satisfactorily explained in writing.

Although there is a difference of emphasis, and it is not clear whether selecting a model on the basis of the data would be interpreted by the regulatory authorities as a change in ”the state of knowledge”, these views appear broadly compatible with the approach described in the present paper.

To sum up: the natural way of operationalizing the approach in the clinical trials context would be to pre-specify the model family and the criterion for selecting a model at an early stage, for example in the protocol or in the analysis plan, and then select the model for the test of no treatment effect in the blind review. This would make model choice objective and so subject to regulatory control if need be. If a significant treatment effect is found, then this should be estimated using the bias correction technique decribed in Section 2. This will involve selection of an appropriate model, this time after unblinding; but again, the choice will be based on objective criteria and so demonstrably free of bias. In this process, if unexpected features are discovered which lead to deviations from the plan, then these should be carefully argued and the preferred analysis defended.

7 Discussion

The proposed procedure can be seen as an extension to the work of Gail et al^{13}, who advocated the use of randomisation inference using a model-based test statistic, on the grounds that this captures the efficiency of the modelbased inference when the model is true, yet retains its nominal size if the model is mis-specified. A similar view was expressed by Brillinger et al.^{25}: ”For crucial comparisons ... the appropriate role for the classical kind of parametric analysis would seem to be confined to assistance in the selection of the test statistics to be used ... in a randomisation analysis.”

I now turn to a broader perspective. The dilemma described in Section 1 is related to a familiar dilemma in applied statistics, viz. that most statistical methods based on models require that the model be fixed in advance, but that choosing a model without looking at the data would often lead to the use of inappropriate models, giving rise to suboptimal, and perhaps incorrect, conclusions.

A simple example of the problems that can arise when a model is chosen from the data is linear regression, with a response Y and a possible explanatory variable X, where models of the type E(Y |x) = α + βx are being considered. A common procedure is to calculate β̂, the least squares estimate for β, and test whether this is significantly different from zero. If it is, then the model E(Y |x) = α + βx is adopted. It is intuitively clear that

E(β̂|β̂ significantly different from 0) ≠ β

so that β̂ is no longer unbiased. This example, and many others, are described in a recent discussion paper by Chatfield^{16}, who writes: ”Perhaps the main message of this paper is that, when a model is formulated and fitted to the same data, inferences made from it will be biased and overoptimistic when they ignore the data-analytic actions which preceded the inference. Statisticians must stop pretending that model uncertainty does not exist and find ways of coping with it.”

Although model selection techniques may be the most flagrant examples of data-driven modelling, they are by no means the only ones. For example, model diagnostic procedures, and techniques for identifying outliers, constitute data-driven model adjustments which lead to bias in the same way as model selection techniques.

Confirmatory clinical trials is one of the few areas in applied statistics for which it is normal practice to pre-specify models, rather than choose them on the basis of the data. Perhaps it is fair to describe the situation in the following way. In confirmatory clinical trials, a model is pre-specified, and it is necessary to pretend that it is true. In most other statistical applications, the choice of model is data-driven, but it is necessary to pretend that it is not.

Chatfield remarks that ”despite the limited progress described above, the overall impression is that the frequentist approach to statistical inference does not adapt easily to cope with model uncertainty.” The present paper has described a way of coping with the problem in the frequentist framework, albeit in the restricted context of randomised studies. Although it does not remove or reduce model uncertainty, it prevents model uncertainty having an adverse impact on inference about the treatment effect.

References

1. ‘Statistical Considerations in the Design of Clinical Trials’. International Conference on Harmonisation of Technical Requirements for Registration of Pharmaceuticals for Human Use (1998).

2. Food and Drug Administration, ‘Guideline for the Format and Content of the Clinical and Statistical Sections of an Application’. FDA, US Department of Health and Human Services, Rockville, MA, USA (1988).

3. CPMP Working Party on Efficacy of Medicinal Products, ‘Biostatistical Methodology in Clinical Trials in Applications for Marketing Authorizations for Medicinal Products’, Statistics in Medicine, 14, 1659-1682 (1995).

4. Rosenbaum, P.R. Observational Studies, Springer-Verlag, New York, 1995.

5. Fisher, R.A. The Design of Experiments, Hafner, New York.

6. Kempthorne, O. ‘Some Aspects of Experimental Inference’, Journal of the American Statistical Association, 61, 11-34 (1966).

7. Tukey, J.W. ‘Tightening the Clinical Trial’, Controlled Clinical Trials, 14, 266-285 (1993).

8. Gail, M.H., Mark, S.D., Carroll, R.J., Green, S.B. and Pee, D. ‘On Design Considerations and randomisation-based inference for Community Intervention Trials’, Statistics in Medicine, 15, 1069-1092 (1996).

9. Romano, J.P. ’On the Behavior of Randomisation Tests without a Group Invariance Assumption’, Journal of the American Statistical Association, 85, 411, 696-692 (1990).

10. Hoeffding, W. ’The Large-Sample Power of Tests Based on Permutations of Observations’, Annals of Mathematical Statistics, 23, 169-192 (1952).

11. Neuhaus, J.M. ‘Misspecification’, Encyclopedia of Biostatistics, Volume 4, 2654-2657. Wiley, Chicester (1998).

12. Keiding, N., Andersen, P.K., Klein, J.P. ‘The role of frailty models and accelerated failure time models in describing heterogeneity due to omitted covariates’, Statistics in Medicine, 16, 215-224 (1997).

13. Gail, M.H., Tan, W.Y. and Piantadosi, S. ‘Tests for No Treatment Effect in randomised Clinical Trials’, Biometrika, 75, 1, 57-64 (1988).

14. Akaike, H., ‘Fitting Autoregressive Models for Prediction’, Annals of the Institute of Statistical Mathematics, 21, 243-247 (1969).

15. Schwarz, G. ‘Estimating the dimension of a model’, Annals of Statistics, 6, 461-4 (1978).

16. Chatfield, C. ‘Model Uncertainty, Data Mining and Statistical Inference (with discussion)’, Journal of the Royal Statistical Society, A, 158, 3, 419-466 (1995).

17. Piantadosi, S. Clinical Trials: A Methodological Perspective, Wiley, New York, 1997.

18. Good, P. Permutation Tests, Springer-Verlag, New York, 1994.

19. Trichler, D. ‘On Inverting Permutation Tests’, Journal of the American Statistical Association, 79, 385, 200-207 (1984).

20. Garthwaite, P.H. ‘Confidence Intervals from Randomisation Tests’, Biometrics, 52, 1387-1393 (1996).

21. Edwards, D. Introduction to Graphical Modelling, Springer-Verlag, New York, 1995.

22. Senn, S. Statistical Issues in Drug Development, Wiley, Chichester, 1997.

23. Tukey, J.W. ‘Use of Many Covariates in Clinical Trials’, International Statistical Review, 59, 2, 123-137 (1991).

24. Lewis, J. A. ‘Statistical Issues in the Regulation of Medicines’, Statistics in Medicine, 14, 127-136 (1995).

25. Brillinger, D.R., Jones, L.V. and Tukey, J.W. ‘The Role of Statistics in Weather Resource Management’, Report of the Statistical Task Force to the Weather Modification Advisory Board (1978).