 Methodology
 Open Access
 Published:
A simulation study on estimating biomarker–treatment interaction effects in randomized trials with prognostic variables
Trials volume 19, Article number: 128 (2018)
Abstract
Background
To individualize treatment decisions based on patient characteristics, identification of an interaction between a biomarker and treatment is necessary. Often such potential interactions are analysed using data from randomized clinical trials intended for comparison of two treatments. Tests of interactions are often lacking statistical power and we investigated if and how a consideration of further prognostic variables can improve power and decrease the bias of estimated biomarker–treatment interactions in randomized clinical trials with timetoevent outcomes.
Methods
A simulation study was performed to assess how prognostic factors affect the estimate of the biomarker–treatment interaction for a timetoevent outcome, when different approaches, like ignoring other prognostic factors, including all available covariates or using variable selection strategies, are applied. Different scenarios regarding the proportion of censored observations, the correlation structure between the covariate of interest and further potential prognostic variables, and the strength of the interaction were considered.
Results
The simulation study revealed that in a regression model for estimating a biomarker–treatment interaction, the probability of detecting a biomarker–treatment interaction can be increased by including prognostic variables that are associated with the outcome, and that the interaction estimate is biased when relevant prognostic variables are not considered. However, the probability of a falsepositive finding increases if too many potential predictors are included or if variable selection is performed inadequately.
Conclusions
We recommend undertaking an adequate literature search before data analysis to derive information about potential prognostic variables and to gain power for detecting true interaction effects and prespecifying analyses to avoid selective reporting and increased falsepositive rates.
Background
Treatment individualization, i.e. finding the right treatment with the right dose at the right time for a specific patient based on certain patient characteristics, is one of the great goals in modern medicine [1]. One requirement for treatment individualization based on, e.g. a certain biomarker like a genetic characteristic or a blood parameter, is the existence of a relevant association between the biomarker and the treatment effect [2], often referred to as the biomarker–treatment interaction. Only a small number of trials have been planned to analyse biomarker–treatment interactions [3], but often the association between one or more biomarkers and a treatment effect is evaluated post hoc in data collected in randomized clinical trials intended for overall comparison of treatment groups, like e.g. the detection of the association between the response to cetuximab and the presence or absence of the Kras mutation in the tumours of patients with advanced colorectal cancer [4].
While often the treatment effect is analysed in different subgroups (prespecified or post hoc specified) to identify patients that benefit from one or another treatment [5], it is widely recognized that the comparison of treatment groups in many different subgroups can lead to spurious results [6]. Therefore, it is often recommended to assess the biomarker–treatment interaction in a regression model, which directly allows us to estimate and test for an interaction effect under common model assumptions [7]. Various authors who provide methods for estimating biomarker–treatment interactions stress the importance of the adequate inclusion of prognostic factors in the model [8, 9]. For treatment effect estimation in a randomized clinical trial, the European Medicines Agency’s guideline on ‘Points to consider on adjustment for baseline covariates’ recommends including other prognostic factors, i.e. covariates that are assumed to be associated with the outcome, as covariates in the regression model to increase the precision of the estimate of the treatment effect [10]. Furthermore, it has been shown that the estimate for the treatment effect is biased in a Cox regression model, if relevant prognostic covariates are not included [11]. While defining the model used for effect estimation and hypothesis testing a priori and including all relevant covariates can be considered as best practices [12], adequate information about prognostic factors might not be available for all research questions, especially when molecular information that has not been well studied and for which limited information from prior investigations is available is included in a regression model. Various approaches to determining the covariates that are to be included in a regression model are presented in the literature [13].
The focus of this article is estimating the interaction between one certain prespecified biomarker of major interest and the treatment. A simulation study was performed to evaluate how the presence and inclusion of further prognostic covariates affect the estimate of the biomarker–treatment interaction. Different strategies for model building, such as including only the main effects of treatment, the biomarker and their interaction, additionally including covariates that are significantly associated with the outcome, or using variable selection methods based on Akaike’s information criterion (AIC) [14] are considered. Scenarios with varying proportions of censored observations, different strengths of association of the prognostic covariates and the outcome, different correlations between prognostic covariates and the biomarker of interest, and different numbers of potential prognostic covariates are considered. The different strategies of covariate inclusion are compared in the control of type I error probabilities and the power to reject the null hypothesis of no biomarker–treatment interaction. A special focus was placed on the socalled rule of ten [12, 15]. This is often considered for predictive models, but (to the best of our knowledge) has not been investigated for the number of additional covariates considered in a regression model, when the primary goal was estimation of an interaction effect.
Methods
Assessing the biomarker–treatment interaction
The interaction between a continuous biomarker of major interest B, or a continuous covariate in general, and treatment T, which is assumed to be binary throughout the article (T∈{0;1}), can be assessed by including an interaction term between the biomarker and the treatment in an adequate regression model. This means the product of B and T is included in the regression model as an additional covariate (see e.g. [13]). The Cox regression model [16], also known as the proportional hazards model, is commonly considered in the analysis of survival data in medical research. In the Cox model, the effect of the biomarker B, the treatment T, their interaction T×B and K other covariates described through the matrix X_{ k } on the hazard rate λ(t) is modelled as
where a linear association between a covariate and the loghazard ratio is assumed. In Eq. (1), λ_{0}(t) is the (unspecified) baseline hazard rate, β_{ T } the regression coefficient for treatment T, β_{ B } the coefficient for the biomarker of interest B, β_{T×B} the regression coefficient for their interaction term and β_{ k } the vector of regression coefficients for the K additional covariates, X_{1},…,X_{ K }. When an interaction term is present, the main effects of the treatment T and the biomarker B can be interpreted as the expected treatment difference at a (fictitious) biomarker value of B=0 and the effect of the biomarker B under treatment T=0 conditional on all other covariates. Regression coefficients are estimated by numerical maximization of the partial loglikelihood PL(β). The variancecovariance matrix of the estimated regression coefficients can be derived as the inverse of the observed information matrix \( {I}^{1}(\boldsymbol {\hat {\beta }})\) (see e.g. [16] or [17] for more details).
Strategies for covariate inclusion
In the simulation study, various approaches for including covariates are compared. In all models, the main effects of the treatment and biomarker as well as their interaction term are included. Obviously, the best choice would be to fit the true model to the data, which includes all covariates that are truly associated with the outcome and ignoring those covariates that are not. This model will be estimated using the simulated data, but in practice the true model will not be known and therefore, the model must be chosen based on plausibility and previous knowledge or based on information gathered from the observed data. Therefore, the following models and strategies were investigated. The names are used for the models/strategies in the figures and tables presented in this article:

Main: A model including only the main effects of treatment T and the biomarker B and their interaction T×B, ignoring all other possible prognostic covariates.

True: A model including the main effects of treatment T, the biomarker of interest B and their interaction T×B, as well as all covariates that are truly associated with the outcome, indicating perfect prior knowledge of relevant covariates.

AIC _{ A }: A model that includes the main effects of treatment T and the biomarker B and their interaction T×B and additionally all covariates that were selected in a forward variable selection procedure based on Akaike’s information criterion (AIC) [14] given T, B and T×B are included (a model including T, B and T×B was used as a starting and minimal model). Additional covariates were selected as long as the AIC criterion
$$ \operatorname{AIC} = 2\operatorname{ll}(\hat{\boldsymbol{\beta}})  2p $$(2)was increased, where \(\operatorname {ll}(\hat {\boldsymbol {\beta }})\) is the partial loglikelihood evaluated at the maximum likelihood estimator \(\hat {\boldsymbol {\beta }}\) and p is the number of estimated regression coefficients.

AIC _{ B }: A modelling strategy similar to AIC _{ A } described above, but prognostic factors were selected based on the AIC criterion considering just the main effect of treatment T as a starting model and not including B or T×B in the variable selection process. After prognostic factors were chosen according to the AIC criterion, B and T×B were added to the model to estimate the biomarker–treatment interaction.

Significance: A model that includes the main effects of treatment T, the covariate of interest and their interaction, as well as all covariates that were significantly associated with the outcome in a Cox regression model including only one covariate (often referred to as univariate Cox models in the medical literature). While this strategy is generally not recommended from a statistical point of view [18], it appears to be a quite popular approach in practice.

Full: A model that includes the treatment T, the biomarker B and their interaction T×B as covariates as well as the main effects of all K potential predictors X_{1},…,X_{ K }.
Data generation and simulation settings
Numerous different settings were considered to evaluate the modelling strategies under varying conditions. For each simulation scenario, 500 subjects were generated. The matrix of continuous covariates (covariate of interest B and potential predictors X_{1},…,X_{ K }) was drawn from a multivariate normal distribution using the R package mvtnorm [19]. For each variable, a mean of 0 and a standard deviation of 1 were used. The correlation structure was specified as described below. Since a randomized controlled trial was intended to be simulated, the treatment variable was drawn independently from all other patient characteristics with Pr(T=1)= Pr(T=0)=0.50 for each individual. For all scenarios, β_{ T } and β_{ B } were chosen as β_{ T }= ln(0.75)=−0.288 (i.e. exp(β_{ T })=0.75) and β_{ B }= ln(1.25)=0.223 (i.e. exp(β_{ B })=1.25).
For each scenario, a timeconstant baseline hazard rate of λ_{0}(t)=1 was used. The hazard rate for each individual was calculated according to Eq. 1 considering the patient’s characteristics and the regression coefficients for the specific scenario. Event times were generated from an exponential distribution using each individual’s hazard rate. All aspects of the simulation study including data generation, estimating regression coefficients and summarizing the results were performed with the statistical software R [20].
The following aspects were varied in the simulation study.
Censoring distribution
Administrative censoring after 5 years was assumed for all scenarios. Additionally, censoring times were generated independently of the event times from an exponential distribution. The hazard rate of the censoring distribution was chosen to produce scenarios with

1.
a low proportion of censored observations (between 30% and 40% censored observations corresponding to 300 to 350 observed events)

2.
a high proportion of censoring (between 60% and 70% censored observations corresponding to 150 to 200 observed events).
Strength of interaction
The strength of the interaction effect between the covariate of interest B and treatment T was varied to consider scenarios with no, quantitative or qualitative biomarker–treatment interaction [21] (see also Fig. 1):

1.
Simulation of data under the null hypothesis of no biomarker–treatment interaction: β_{T×B}=0.

2.
Quantitative biomarker–treatment interaction with a difference in the magnitude of the treatment effect between individuals with a low value of B and individuals with a large value of B: β_{T×B}= ln(1.1)=0.095, leading to a hazard ratio between the treatment groups (T=1 vs. T=0) of about 0.6 for a given value of B=−2 and a hazard ratio of about 0.9 for B=2.

3.
Qualitative biomarker–treatment interaction indicating an expected lower risk for an event from treatment T=1 for patients with a small value of B and a lower risk under treatment T=0 for patients with a large value of B: β_{T×B}= ln(1.33)=0.285, providing a hazard ratio between the treatment groups smaller than 1 for B<1 and a hazard ratio larger than 1 for B>1 (dotted line in Fig. 1).
Number of potential prognostic variables to be included in the model Three settings for the number K of potential candidate predictors that can be included in the regression model were considered:

1.
K=12: Here 12 additional prognostic covariates are considered, so the rule of ten is fulfilled under both censoring distributions for most simulation runs, as 150 to 200 events are expected in the settings with a high amount of censoring and up to 15 regression coefficients are to be estimated (12 prognostic variables plus the main effects of treatment T and the covariate of interest B and their interaction T×B).

2.
K=24: Here 24 additional prognostic covariates are considered, so the rule of ten will be violated for most scenarios with high censoring.

3.
K=36: Here 36 additional prognostic covariates are considered. Again, the rule of ten will be violated under high censoring.
Correlation structure between prognostic variables and covariate of interest Three different correlation structures between the covariate of interest B and the potential prognostic variables X_{1},…,X_{ K } were considered:

1.
Firstly, a scenario with a biomarker of interest B that is independent of the potential prognostic variables, and independence between all the prognostic variables was investigated, with
$$\Sigma_{1} = \left(\begin{array}{ccccc} 1 & 0 & \cdots & \cdots & 0 \\ 0 & 1 & 0 & \cdots & 0 \\ \vdots & & \ddots & &\vdots \\ 0 & \cdots & 0 & 1 & 0 \\ 0 & \cdots & \cdots & 0 & 1 \\ \end{array} \right). $$ 
2.
As a second setting, the correlation coefficients between B and all other covariates X_{1},…,X_{ K }, as well as between each pair of covariates X_{ i },X_{ j } with i≠j was set to r=0.5, indicating a moderate correlation between all variables:
$$\Sigma_{2} = \left(\begin{array}{ccccc} 1 & 0.5 & \cdots & \cdots & 0.5 \\ 0.5 & 1 & 0.5 & \cdots & 0.5 \\ \vdots & & \ddots & &\vdots \\ 0.5 & \cdots & 0.5 & 1 & 0.5 \\ 0.5 & \cdots & \cdots & 0.5 & 1 \\ \end{array} \right). $$ 
3.
A block correlation structure between the covariates was considered, with a high correlation of r=0.7 between the biomarker B and a set of variables as well as between those variables, a moderate correlation of r=0.4 for another set and a correlation of r=0.1 or r=0 for the other variables:
$$\Sigma_{3} = \left(\begin{array}{ccccccccccccc} 1 & 0.7 & 0.7 & 0.7 & 0.4 & 0.4 & 0.4 & 0.4 & 0.1 & 0.1 & 0.1 & 0 & 0 \\ 0.7 & 1 & 0.7 & 0.7 & 0.4 & 0.4 & 0.4 & 0.4 & 0.1 & 0.1 & 0.1 & 0 & 0 \\ 0.7 & 0.7 & 1 & 0.7 & 0.4 & 0.4 & 0.4 & 0.4 & 0.1 & 0.1 & 0.1 & 0 & 0 \\ 0.7 & 0.7 & 0.7 & 1 & 0.4 & 0.4 & 0.4 & 0.4 & 0.1 & 0.1 & 0.1 & 0 & 0 \\ 0.4 & 0.4 & 0.4 & 0.4 & 1 & 0.4 & 0.4 & 0.4 & 0.1 & 0.1 & 0.1 & 0 & 0 \\ 0.4 & 0.4 & 0.4 & 0.4 & 0.4 & 1 & 0.4 & 0.4 & 0.1 & 0.1 & 0.1 & 0 & 0 \\ 0.4 & 0.4 & 0.4 & 0.4 & 0.4 & 0.4 & 1 & 0.4 & 0.1 & 0.1 & 0.1 & 0 & 0 \\ 0.4 & 0.4 & 0.4 & 0.4 & 0.4 & 0.4 & 0.4 & 1 & 0.1 & 0.1 & 0.1 & 0 & 0 \\ 0.1 & 0.1 & 0.1 & 0.1 & 0.1 & 0.1 & 0.1 & 0.1 & 1 & 0.1 & 0.1 & 0 & 0 \\ 0.1 & 0.1 & 0.1 & 0.1 & 0.1 & 0.1 & 0.1 & 0.1 & 0.1 & 1 & 0.1 & 0 & 0 \\ 0.1 & 0.1 & 0.1 & 0.1 & 0.1 & 0.1 & 0.1 & 0.1 & 0.1 & 0.1 & 1 & 0 & 0 \\ 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1 & 0 \\ 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0 & 1 \end{array}\right). $$
For the scenarios with K=24 or K=36 potential predictors, the correlation matrices were adapted accordingly.
Strength of association between prognostic variables and outcome
For the strength of association between the potential prognostic variables X_{1},…,X_{ K } and the outcome, two different settings were chosen:

1.
For all covariates X_{1},…,X_{ K }, the same regression coefficient was chosen:
$$\boldsymbol{\beta_{k}} = \boldsymbol{\beta_{eq}} = (\ln(1.1), \ldots, \ln(1.1))^{T} = (0.095, \ldots, 0.095)^{T}. $$ 
2.
Varying strengths of association between the potential predictors and the risk for an event were considered. The vector of regression coefficients was chosen to be
$$\boldsymbol{\beta_{k}} = \boldsymbol{\beta_{v}} = \left(\begin{array}{cc} \ln(1.2) \\ \ln(1.1)\\ \ln(1)\\ \ln(1.2)\\ \ln(1.1)\\ \ln(1)\\ \vdots\\ \ln(1.2)\\ \ln(1.1)\\ \ln(1) \end{array} \right) = \left(\begin{array}{cc} 0.182 \\ 0.095\\ 0\\ 0.182\\ 0.095\\ 0\\ \vdots\\ 0.182\\ 0.095\\ 0 \end{array} \right). $$
As all combinations of the different settings described above were considered in the simulation study, a total of 2 censoring distributions × 3 strengths of interaction between biomarker B and treatment T × 3 numbers of potential prognostic variables × 3 different correlation structures × 2 settings for association between the potential prognostics variables and the outcome = 108 settings were considered in the simulation. For each of these settings, 1000 simulation runs were performed.
Analysis and presentation of results
In each simulation run, all of the methods or strategies described in “Strategies for covariate inclusion” section were fitted or applied, respectively. Estimation of the regression coefficients from the Cox regression models was performed with the function coxph in the survival library [22] of the statistical software R [20]. For the variable selection based on the AIC criterion, the function stepAIC in the library MASS [23] was applied.
For each model in each simulation run, the estimated regression coefficient for the biomarker–treatment interaction term \(\hat {\beta }_{T\times B}\) and its estimated variance as well as the p value of the Wald test for the null hypothesisv H_{0}: β_{T×B}=0 was saved. Additionally, a 95% confidence interval for β_{T×B} was estimated as
where ϕ_{0.975} denotes the 97.5% quantile of the standard normal distribution and \(\widehat {\operatorname {var}}(\hat {\beta }_{T\times B})\) is the estimated variance of the interaction coefficient obtained in the corresponding simulation run for the respective modelling approach. If the algorithm for numerical maximization of the partial loglikelihood did not converge, this information was saved. All results presented in ‘Results’ rely on only estimations for which the numerical optimization algorithm converged. The number of runs for which no result was returned is presented.
For each model and strategy, the confidence interval coverage, i.e. the fraction of simulation runs in which the estimated confidence interval for the biomarker–treatment interaction covered the true value, was derived. The proportion of simulation runs in which the null hypothesis was rejected and a statistically significant biomarker–treatment interaction was detected for the conventional significance level of 5%, i.e. the power of the statistical test if H_{0} were false or the probability of a type I error if H_{0} were true (β_{T×B}=0), was determined [24].
Results
The observed proportions of rejected null hypotheses are summarized in Table 1. Results are presented stratified for different values of K, strength of interaction and proportion of censored observations, but were aggregated over different values of β_{ k } and Σ. In Tables 2, 3 and 4, the observed proportions of simulation runs with rejected null hypotheses are shown separately for the scenarios with K=12 (Table 2), K=24 (Table 3) and K=36 (Table 4), for scenarios with no true biomarker–treatment interaction (top), true quantitative biomarker–treatment interaction (middle) and true qualitative biomarker–treatment interaction (bottom). An observed type I error probability of 7% was considered to be acceptable. For scenarios with no interaction (β_{ T×B }=0), observed type I error proportions larger than 7% are in italics. For scenarios with data generated under H_{1} (quantitative interaction and qualitative interaction), the proportions of rejected null hypotheses are in bold if the type I error probability for the approach at the given scenario was not larger than 7%.
The mean numbers of included additional covariates are given for each method or strategy for sets of scenarios stratified for β_{ k } and amount of censoring in the bottom rows of Tables 2, 3 and 4 and for each of the 108 simulated scenarios in Additional file 7: Table S1 (for K=12), Additional file 8: Table S2 (for K=24) and Additional file 9: Table S3 (for K=36).
The distributions of the obtained estimates are illustrated in Fig. 2 for one exemplary set of scenarios. The observed distributions of the regression coefficient estimates for the biomarker–treatment interaction \(\hat {\beta }_{T\times B}\) are displayed as box plots for the scenarios with Σ=Σ_{ 3 }, β_{ k }=β_{ v } and low (a) or high number of censored observations (b). In the top rows, scenarios with no true biomarker–treatment interaction are shown, and in the bottom rows, results for data simulated with true qualitative biomarker–treatment interactions are presented. Scenarios with different numbers of (potential) prognostic variables (K=12, K=24 and K=36) are shown in separate columns. Distributions of estimated regression coefficients are illustrated for all scenarios with no true interaction (under H_{0}) or with true qualitative interaction in Additional file 1: Figure S1, Additional file 2: Figure S2, Additional file 3: Figure S3, Additional file 4: Figure S4, Additional file 5: Figure S5 and Additional file 6: Figure S6. In each figure, the true value of the interaction regression coefficient is illustrated by the horizontal red line. Additionally, the confidence interval coverage for each modelling strategy (triangles and blue lines) and the probability of rejection of the null hypothesis of no biomarker–treatment interaction, i.e. the estimated probability for a type I error in the first row and the observed power in the second row, are illustrated (circles and green lines).
The type I error probabilities for the biomarker–treatment interaction term, which are presented in the lines indicated with no interaction (Table 1) and in the upper parts of Tables 2, 3 and 4 for the scenarios with no interaction, were acceptable for almost all methods and strategies, when K=12 further (potential) prognostic variables were considered. Only for strategy AIC _{ A } an unacceptably high probability of type I errors (defined as larger than 7%) was observed for one setting (Table 2). For scenarios with K=24 (potential) prognostic variables, increased type I error probabilities were observed for each method for at least one scenario, except for AIC _{ B }. For AIC _{ A }, type I error probabilities above 7% were observed for six of the 12 settings (Table 3) and for scenarios with a high proportion of censored observations (60% to 70%) when scenarios with different β_{ k } and Σ were aggregated (Table 1). When K = 36 potential predictors were considered, an increased type I error probability was observed for AIC _{ A } for all scenarios. For main, significance and full, elevated false positive rates were obtained for three to five scenarios with a high proportion of censored observations. For the true model, only two scenarios with a high proportion of censored observations led to rejection of the null hypothesis in more than 7% of the observed simulation runs (Σ_{ 2 }, β_{ eq } and Σ_{ 1 }, β_{ eq }). For all other scenarios, the observed type I error probabilities were between 5% and 7%. For the strategy AIC _{ B }, all observed type I error probabilities were between 5% and 7%.
For the main model, regression coefficients for the biomarker–treatment effect were underestimated when a true biomarker–treatment interaction was present (Fig. 2), with the largest bias observed for scenarios with Σ=Σ_{ 2 } (see second rows of Additional file 1: Figure S1A and Figure S1B, Additional file 2: Figure S2A and Figure S2B, Additional file 3: Figure S3A and Figure S3B, Additional file 4: Figure S4A and Figure S4B, Additional file 5: Figure S5A and Figure S5B, and Additional file 6: Figure S6A and Figure S6B). This also led to a loss of power, which was reduced as compared to the true model for most of the scenarios (Tables 1, 2, 3 and 4, green dots in Additional file 1: Figure S1, Additional file 2: Figure S2, Additional file 3: Figure S3, Additional file 4: Figure S4, Additional file 5: Figure S5 and Additional file 6: Figure S6). Generally, the highest power was observed for the true model. The power for AIC _{ A } cannot be interpreted adequately for most of the scenarios due to its increased type I error probabilities. The full model is identical to the true model for β_{ k }=β_{ eq }, as all covariates are truly associated with the outcome. For β_{ k }=β_{ v }, the power of the full model was similar to the power of the true model for K=12 and K=24 in our simulation runs, but was slightly lower for simulations with K=36. The strategy AIC _{ B }, which appears to have an adequate false positive rate, showed (slightly) lower power than the true model for (almost) all of the scenarios. A slightly decreased power was also observed for the strategy including all covariates that were significantly associated with the outcome (significance). The type I error probability was acceptable for most scenarios with a small or moderate number of potential predictors (K=12 and K=24), but an increased type I error probability was observed for scenarios with many potential predictors (K=36).
Coverage was adequate for most of the models and strategies. For main, the coverage was reduced for some scenarios due to biased estimates. For AIC _{ A }, the coverage was under 93% for 52 of the 108 scenarios (48.1%), indicating standard errors for the regression coefficient of interest were underestimated following the variable selection procedure.
In the last rows of Tables 2, 3 and 4, the mean numbers of additionally included covariates are summarized for each method/strategy stratified for the amount of censoring and β_{ k } (which determines the number of truly prognostic variables). It was observed that for our settings, the procedure including variables that were significantly associated with the outcome in univariate Cox models selected more variables than the AICbased methods, and that slightly more variables were chosen with AIC _{ B } than with AIC _{ A }. For scenarios with β_{ k }=β_{ eq }, the true and full models were identical by definition. More detailed information on the numbers of covariates included are given in Additional file 7: Table S1, Additional file 8: Table S2 and Additional file 9: Table S3.
The optimization algorithm for numerical maximization of the partial loglikelihood of the Cox regression model for estimating the regression coefficients did not converge for some simulation runs. The problem especially occurred for AIC _{ A }. Over all 108,000 simulation runs (108 scenarios × 1,000 runs per scenario), the estimation algorithm did not converge 11 times (0.010%) for main, twice (0.002%) for true, 895 times (0.829%) for AIC _{ A }, 27 times (0.025%) for AIC _{ B }, three times (0.003%) for significance and no times (0%) for full.
Discussion
The ultimate goal in individualized or tailored medicine is to find the best treatment for each individual based on the patient’s characteristics like age, sex, comorbidities, disease history and molecular and genetic information, which are often referred to as biomarkers. The existence and detection of a biomarker–treatment interaction can be considered as a requirement for such treatment individualization [2], and consequently an interaction between the biomarker of interest and treatment has to be established in a first step, e.g. by finding statistically significant and clinically relevant interactions based on data from (multiple) randomized clinical trials. Decision rules for treatment selection based on the characteristics of a certain patient have to be investigated and established afterwards, also considering the benefits and costs of the application of a certain treatment strategy for a given patient.
To detect relevant associations and interactions, it is well known that splitting a quantitative variable into different categories, leading to a comparison of treatment effects between different subgroups, will result in a loss of information and will consequently decrease the probability of detecting a true biomarker–treatment interaction [25]. So, using all the quantitative information is recommended for analysis of biomarker–treatment interactions [7]. To estimate a treatment effect in a randomized clinical trial, the inclusion of relevant prognostic variables is recommended [10] to increase the precision of the estimate and consequently the probability of detecting real group differences. For this article, we performed a simulation study to investigate whether the probability of detecting a biomarker–treatment interaction in data derived from a randomized clinical trial can be improved by including further potentially prognostic variables in a Cox regression model for timetoevent data. Different settings for the strength of interaction between the biomarker and the treatment, the correlation between the biomarker of interest and other potential predictors, the strength of association between the predictors and outcome, the number of (potential) further predictors, and the number of events and censored observations were considered. When a biomarker–treatment interaction is assessed using data from a randomized clinical trial, obviously the best choice is to include in the final model all covariates truly associated with the outcome, which was covered by the true model in our simulation study. As this true model often is not known in practice, especially in investigations including molecular or genetic information, more flexible approaches might be needed. So, we also investigated strategies using datadriven variable selection procedures based on AIC [14] or on the results of Cox regression models with single covariates.
In our simulation study, we observed that including the correct prognostic variables leads to an increased probability of detecting a true biomarker–treatment interaction and reduced bias of the estimated interaction effect, with the magnitude of improvement depending on the strength of association between the prognostic variables and the outcome and between the prognostic variables and the biomarker of interest. In contrast, including too many variables per event can lead to the opposite effect and increased probabilities of false positives. This problem is well known for multiple regression models [15, 26]. Our results support the rule of ten, which was proposed for predictive modelling [27], since the type I error probability was increased for the biomarker of interest, even for the true model, when a large number of covariates was considered. The simulation study also revealed that ignoring relevant prognostic factors leads to biased estimates for the biomarker–treatment interaction effect, which has been described for estimating the group effect from a randomized clinical trial using a Cox regression model [11]. Generally, the datadriven selection of prognostic variables by an inclusion procedure based on the AIC after including the main effects of the biomarker of interest, the treatment and their interaction in the model increases the type I error probabilities and reduces the confidence interval coverage. This was not observed in a strategy that selected the relevant prognostic variables in a first step and added the biomarker main effect and the biomarker–treatment interaction afterwards (called AIC _{ B } in our article). In our simulated scenarios, the strategy including all covariates that were found to be significantly associated with the outcome performed similarly to that approach. Automated variable selection procedures are criticized in the literature for various reasons (see e.g. [28]). Based on the results of our simulation study, we strongly discourage using an automated variable selection procedure to choose additional prognostic variables after including the biomarker–treatment interaction of interest, as this may lead to unreliable results.
An obvious limitation of our study is that we observed only a moderate number of different scenarios with three correlation structures, three strengths of interaction between the biomarker and treatment, two strengths/structures of association between the additional prognostic variables and treatment, two censoring distributions, three numbers of (potential) prognostic variables, and a fixed number of 500 observations, due to limited time and space. All these aspects influenced the results and other settings may have led to different findings and consequently recommendations. In particular, the number of observed events, which is more important than the total sample size for a timetoevent outcome, was varied only by choosing two different censoring proportions, but it has a major impact on the power of the interaction test. We also investigated only a small number of strategies for inclusion or selection of further covariates based on the AIC and significant associations with the outcome. Other strategies (like backward selection), other criteria (like the Bayesian information criterion [29]) or other procedures for variable selection (like the least absolute shrinkage and selection operator [30]) were not considered. Furthermore, we considered only normally distributed biomarkers and linear associations and interactions in our simulations and fitted Cox regression models assuming linear associations and timeconstant effects to our data. Recently introduced methods for estimating nonlinear interactions, like local partial likelihood estimation [31], multivariable fractional polynomials for interaction [8] or the modified covariate approach [9], were not investigated.
It has to be considered that in our scenario, only one prespecified biomarker of interest is assessed. It was identified as being of interest e.g. in an observational study or was found to be relevant for a similar kind of disease. If more than one biomarker is investigated, multiplicity issues arise that have to be adequately considered [32]. When an analysis is an additional analysis to a standard group comparison for a randomized clinical trial, it can only be exploratory in nature. Nevertheless, the method used for statistical analysis should be specified a priori to generate reliable results and avoid problems of datadredging and selective reporting, and consequently generating unreliable results and increased falsepositive rates [33]. Further algorithms or strategies should be used in sensitivity analyses to assess the stability of the observed results. If the investigation of a biomarker–treatment interaction is of major importance for a clinical trial, this should be considered in the design stage and consequently in the sample size calculation.
Conclusions
Based on the results of our simulation study, we recommend considering prognostic covariates in regression models when estimating biomarker–treatment interactions, as the power for detecting true interactions can be increased. However, including too many variables can lead to unreliable results. The choice of variables included should be based on prior information and subject knowledge. Automatic variable selection procedures have to be handled with care.
References
Hamburg MA, Collins FS. The path to personalized medicine. N Engl J Med. 2010; 363(4):301–4.
Chen JJ, Lu TP, Chen YC, Lin WJ. Predictive biomarkers for treatment selection: statistical considerations. Biomark Med. 2015; 9(11):1121–35.
Rothwell PM. Subgroup analysis in randomised controlled trials: importance, indications, and interpretation. Lancet. 2005; 365(9454):176–86.
Karapetis CS, KhambataFord S, Jonker DJ, O’Callaghan CJ, Tu D, Tebbutt NC, et al.Kras mutations and benefit from cetuximab in advanced colorectal cancer. N Engl J Med. 2008; 359(17):1757–65.
Assmann SF, Pocock SJ, Enos LE, Kasten LE. Subgroup analysis and other (mis)uses of baseline data in clinical trials. Lancet. 2000; 355(9209):1064–9.
Naggara O, Raymond J, Guilbert F, Roy D, Weill A, Altman DG. Analysis by categorizing or dichotomizing continuous variables is inadvisable: an example from the natural history of unruptured aneurysms. Am J Neuroradiol. 2011; 32(3):437–40.
Royston P, Sauerbrei W. Interactions between treatment and continuous covariates: a step toward individualizing therapy. J Clin Oncol. 2008; 26(9):1397–9.
Royston P, Sauerbrei W. A new approach to modelling interactions between treatment and continuous covariates in clinical trials by using fractional polynomials. Stat Med. 2004; 23(16):2509–25.
Tian L, Alizadeh AA, Gentles AJ, Tibshirani R. A simple method for estimating interactions between a treatment and a large number of covariates. J Am Stat Assoc. 2014; 109(508):1517–32.
Committee for Proprietary Medicinal Products. Points to consider on adjustment for baseline covariates. Stat Med. 2004; 23(5):701.
Langner I, Bender R, LenzTönjes R, Küchenhoff H, Blettner M. Bias of maximumlikelihood estimates in logistic and Cox regression models: a comparative simulation study. Sonderforschungsbereich 386. LudwigMaximiliansUniversität München; 2003.
Harrell F. Regression modeling strategies: with applications to linear models, logistic and ordinal regression, and survival analysis. New York: Springer; 2001.
Royston P, Sauerbrei W. Multivariable modelbuilding: a pragmatic approach to regression analysis based on fractional polynomials for modelling continuous variables, Vol. 777. Chichester: Wiley; 2008.
Akaike H. A new look at the statistical model identification. IEEE Trans Autom Control. 1974; 19(6):716–23.
Babyak MA. What you see may not be what you get: a brief, nontechnical introduction to overfitting in regressiontype models. Psychosom Med. 2004; 66(3):411–21.
Cox DR. Regression models and life tables (with discussion). J Royal Stat Soc. 1972; 34:187–220.
Therneau TM, Grambsch PM. Modeling survival data: extending the Cox model. New York: Springer; 2013.
Vach W. Regression models as a tool in medical research. Boca Raton: CRC Press; 2012.
Genz A, Bretz F, Miwa X, Tetsuhisa abd Mi, Leisch F, Scheipl F, Hothorn T. Mvtnorm: multivariate normal and T distributions. R package version 1.05. 2016. http://CRAN.Rproject.org/package=mvtnorm.
R Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2016. https://www.Rproject.org/.
Polley MYC, Freidlin B, Korn EL, Conley BA, Abrams JS, McShane LM. Statistical and practical considerations for clinical evaluation of predictive biomarkers. J Natl Cancer Inst. 2013; 105(22):1677–83.
Therneau T. A package for survival analysis in S. Version 2.38. 2015. http://CRAN.Rproject.org/package=survival.
Venables WN, Ripley BD. Modern applied statistics with S, 4th ed. New York: Springer; 2002. http://www.stats.ox.ac.uk/pub/MASS4.
Burton A, Altman DG, Royston P, Holder RL. The design of simulation studies in medical statistics. Stat Med. 2006; 25(24):4279–92.
Royston P, Altman DG, Sauerbrei W. Dichotomizing continuous predictors in multiple regression: a bad idea. Stat Med. 2006; 25(1):127–41.
Concato J, Peduzzi P, Holford TR, Feinstein AR. Importance of events per independent variable in proportional hazards analysis. I. Background, goals, and general strategy. J Clin Epidemiol. 1995; 48(12):1495–501.
Peduzzi P, Concato J, Feinstein AR, Holford TR. Importance of events per independent variable in proportional hazards regression analysis. II. Accuracy and precision of regression estimates. J Clin Epidemiol. 1995; 48(12):1503–10.
Sainani KL. Multivariate regression: the pitfalls of automated variable selection. PM&R. 2013; 5(9):791–4.
Schwarz G, et al. Estimating the dimension of a model. Annals Stat. 1978; 6(2):461–4.
Tibshirani R. Regression shrinkage and selection via the lasso. J R Stat Soc Ser B Methodol. 1996; 58(1):267–288.
Liu Y, Jiang W, Chen BE. Testing for treatmentbiomarker interaction based on local partiallikelihood. Stat Med. 2015; 34(27):3516–30.
European Medicines Agency. Guideline on multiplicity issues in clinical trials. 2017. http://www.ema.europa.eu/ema/index.jsp?curl=pages/regulation/general/general_content_001220.jsp&mid=.
Ioannidis JP. Why most published research findings are false. PLoS Med. 2005; 2(8):124.
Acknowledgments
Not applicable
Funding
This work was supported by the German Research Foundation (DFG) and the Technische Universität München within the funding programme Open Access Publishing.
Availability of data and materials
No patient data were used when writing this article. The R code for generating the data sets used in the simulation study, for applying the approaches and strategies described, and for analysing the results obtained in the simulation study can be obtained from the first author upon reasonable request.
Author information
Authors and Affiliations
Contributions
BH designed and implemented the simulation study and drafted the manuscript. KU critically reviewed the manuscript for intellectual content. Both authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable
Consent for publication
Not applicable
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional files
Additional file 1
Figure S1. Distribution of \(\hat {\beta }_{T\times B}\) for scenarios with K=12, β_{ k } = β_{ eq }, and low censoring (A) or high censoring (B) for no biomarker–treatment interaction (β_{T×B}= ln(1.0)=0, top rows) or qualitative biomarker–treatment interaction (β_{T×B}= ln(1.33)=0.285, bottom rows). Results for different correlation structures are shown in separate columns. The dashed red lines indicate the true value of β_{T×B}, the blue triangles represent the observed confidence interval coverages, the green dots the observed probability for a type I error (A) or estimated power (B). (PDF 20 kb)
Additional file 2
Figure S2. Distribution of \(\hat {\beta }_{T\times B}\) for scenarios with K=12, β_{ k } = β_{ v }, and low censoring (A) or high censoring (B) for no biomarker–treatment interaction (β_{T×B}= ln(1.0)=0, top rows) or qualitative biomarker–treatment interaction (β_{T×B}= ln(1.33)=0.285, bottom rows). Results for different correlation structures are shown in separate columns. The dashed red lines indicate the true value of β_{T×B}, the blue triangles represent the observed confidence interval coverages, the green dots the observed probability for a type I error (A) or estimated power (B). (PDF 20 kb)
Additional file 3
Figure S3. Distribution of \(\hat {\beta }_{T\times B}\) for scenarios with K=24, β_{ k } = β_{ eq }, and low censoring (A) or high censoring (B) for no biomarker–treatment interaction (β_{T×B}= ln(1.0)=0, top rows) or qualitative biomarker–treatment interaction (β_{T×B}= ln(1.33)=0.285, bottom rows). Results for different correlation structures are shown in separate columns. The dashed red lines indicate the true value of β_{T×B}, the blue triangles represent the observed confidence interval coverages, the green dots the observed probability for a type I error (A) or estimated power (B). (PDF 20 kb)
Additional file 4
Figure S4. Distribution of \(\hat {\beta }_{T\times B}\) for scenarios with K=24, β_{ k } = β_{ v }, and low censoring (A) or high censoring (B) for no biomarker–treatment interaction (β_{T×B}= ln(1.0)=0, top rows) or qualitative biomarker–treatment interaction (β_{T×B}= ln(1.33)=0.285, bottom rows). Results for different correlation structures are shown in separate columns. The dashed red lines indicate the true value of β_{T×B}, the blue triangles represent the observed confidence interval coverages, the green dots the observed probability for a type I error (A) or estimated power (B). (PDF 20 kb)
Additional file 5
Figure S5. Distribution of \(\hat {\beta }_{T\times B}\) for scenarios with K=36, β_{ k } = β_{ eq }, and low censoring (A) or high censoring (B) for no biomarker–treatment interaction (β_{T×B}= ln(1.0)=0, top rows) or qualitative biomarker–treatment interaction (β_{T×B}= ln(1.33)=0.285, bottom rows). Results for different correlation structures are shown in separate columns. The dashed red lines indicate the true value of β_{T×B}, the blue triangles represent the observed confidence interval coverages, the green dots the observed probability for a type I error (A) or estimated power (B). (PDF 20 kb)
Additional file 6
Figure S6. Distribution of \(\hat {\beta }_{T\times B}\) for scenarios with K=36, β_{ k } = β_{ v }, and low censoring (A) or high censoring (B) for no biomarker–treatment interaction (β_{T×B}= ln(1.0)=0, top rows) or qualitative biomarker–treatment interaction (β_{T×B}= ln(1.33)=0.285, bottom rows). Results for different correlation structures are shown in separate columns. The dashed red lines indicate the true value of β_{T×B}, the blue triangles represent the observed confidence interval coverages, the green dots the observed probability for a type I error (A) or estimated power (B). (PDF 20 kb)
Additional file 7
Table S1. Mean number of additionally included prognostic variables for all scenarios with K=12. (PDF 68 kb)
Additional file 8
Table S2. Mean number of additionally included prognostic variables for all scenarios with K=24. (PDF 68 kb)
Additional file 9
Table S3. Mean number of additionally included prognostic variables for all scenarios with K=36. (PDF 68 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Haller, B., Ulm, K. A simulation study on estimating biomarker–treatment interaction effects in randomized trials with prognostic variables. Trials 19, 128 (2018). https://doi.org/10.1186/s1306301824910
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1306301824910
Keywords
 Biomarker–treatment interaction
 Randomized trial
 Stratified medicine
 Predictive covariates
 Variable selection