Estimating the distribution of heterogeneous treatment effects from treatment responses and from a predictive biomarker in a parallel‐group RCT: A structural model approach

When the objective is to administer the best of two treatments to an individual, it is necessary to know his or her individual treatment effects (ITEs) and the correlation between the potential responses (PRs) Yi1 and Yi0 under treatments 1 and 0. Data that are generated in a parallel‐group design RCT does not allow the ITE to be determined because only two samples from the marginal distributions of these PRs are observed and not the corresponding joint distribution. This is due to the “fundamental problem of causal inference.” Here, we present a counterfactual approach for estimating the joint distribution of two normally distributed responses to two treatments. This joint distribution of the PRs Yi1 and Yi0 can be estimated by assuming a bivariate normal distribution for the PRs and by using a normally distributed baseline biomarker Zi functionally related to the sum Yi1+Yi0 . Such a functional relationship is plausible since a biomarker Zi and the sum Yi1+Yi0 encode for the same information in an RCT, namely the variation between subjects. The estimation of the joint trivariate distribution is subjected to some constraints. These constraints can be framed in the context of linear regressions with regard to the proportions of variances in the responses explained and with regard to the residual variation. This presents new insights on the presence of treatment–biomarker interactions. We applied our approach to example data on exercise and heart rate and extended the approach to survival data.


INTRODUCTION
Numerous examples in the medical literature demonstrate the existence of different treatments for the same disease with similar or different modes of action. While in some indications, one treatment is uniformly superior to the other, there are also cases, where a percentage of subjects benefited more from the one specific treatment compared to alternatives. This is called a heterogeneous treatment effect or subject-treatment interaction. The identification of these effects within respective subjects allows selecting the appropriate individual treatment. These effects play a key role in the concept of personalized medicine.
Starting with a randomized controlled trial (RCT) that allocates two treatments denoted by 1 (for the new treatment) and 0 (for the standard treatment or placebo), let be the potential response (PR) to the treatment of subject (with ∈ {0, 1}). However, we cannot simultaneously observe both responses 1 and 0 on subject due to the "fundamental problem of causal inference" (Holland, 1986, p. 947). Formally, the observed response is given by where ∈ {0; 1} denotes the treatment allocated to subject i. This equation holds true due to the consistency assumption, which binds the potential to the observed response. Nonetheless, it is of interest to know the joint distribution of the PRs from which the "individual treatment effect" (ITE) 1 − 0 is derived and to determine the subject's probability to benefit from which treatment best. This is a decision aid for a physician to make the best treatment choice within two (or several) treatment options. But due to the "fundamental problem of causal inference" (Holland, 1986, p. 947), only the average of the expected difference of the PRs can be estimated when no further information is available. The term is commonly known as "average treatment effect" (ATE) and can be estimated from the observed marginals of the PRs. If the ITE, 1 − 0 , is not constant within the trial population, an interaction between treatments and the subjects (subject-treatment interaction) is present. It is possible to approximate the joint distribution of the PRs in a replicated randomized crossover design. This design is commonly used for evaluating individual bioequivalence. It allows separate estimates of between-subject variation, subjecttreatment interactions (in the bioequivalence context called "subject-formulation interaction") and within-subject variation (Senn, 2001). Unreplicated crossover trials can also be used to determine the correlation between responses under two treatments. A positive correlation reflects treatments with similar "modes of action," a negative correlation reflects treatments with different "modes of actions" (Cleophas, 2000).
In situations where a crossover design is not feasible (e.g., eradication therapies in infectious diseases and cancer), a parallelgroup RCT remains the only option to compare treatment effects. Here, subject effects (between-subject variation) can be approximated by biomarkers (known as "prognostic biomarkers") and subject-treatment interactions can be approximated by treatment-biomarker interactions (known as "predictive biomarkers") as proposed by Senn (2001). Predictive biomarker guide treatment options and are labeled "companion diagnostic (cDx)." The most prominent example is the biomarker Her-2/neu. If a breast cancer patient is treated with trastuzumab depends on the value measured (Hudis, 2007).
A treatment-biomarker interaction allows assessing if a biomarker is predictive or not. Here, two kinds of interactions are distinguished: quantitative and qualitative interaction. A quantitative interaction is observed if the treatment effects differ with respect to their size but not with respect to their sign over the range of biomarker values. A qualitative interaction is observed if the treatment effects differ with respect to their size as well as their sign over the range of biomarker values. From a treatment point of view, only qualitative interactions are of interest since there are subgroups of patients who benefit from treatment 1 as well as subgroups of patients who benefit from treatment 0. In case of a quantitative interaction, all patients will benefit from the same treatment but in different extent.
These effects can be seen by looking at the conditional treatment effect. It is the average treatment difference given a specific biomarker value. In general, it is estimated by the difference of two regression functions which fit the biomarker value to the specific treatment response. The goal of this paper is to sharpen this concept by introducing the correlation between the two individual potential responses in order to predict the treatment effect for an individual with a given biomarker value. For this aim, we derive prediction intervals. The consequence of not considering the dependence between the responses 1 and 0 is shown in the bottom left plot of Figure 1. In this plot, the prediction intervals for a positive correlation between 1 and 0 ( 10 = 0.8) (dashed lines) are much closer around the conditional treatment effects compared to the prediction intervals where independence (in this case thus a correlation 10 = 0.0) is assumed (dotted lines). The prediction intervals are wider for a negative correlation ( 10 = −0.5) (dashed-dotted lines) than those based on a zero correlation and much wider than those based on a positive correlation.
Homogeneous treatment effects are present for 10 = 1. As a consequence, the treatment effect is additive and is constant in the trial population. Heterogeneous treatment effects are present if 10 ≠ 1. Here, treatment effects differ among the patients. Knowing the joint distribution of the responses 1 and 0 formalizes the aims of "personalized medicine": Applying the best F I G U R E 1 In the top left, the outcome for two different treatments is depicted, is the value of the biomarker where both treatments are equally effective. The top right figure shows the ATE conditional on the biomarker together with the confidence and the prediction interval under the assumption of 10 = 0. The bottom left figure shows prediction intervals for different values of 10 . The bottom right figure shows the probailities ( 1 > 0 ) for different values of 10 F I G U R E 2 Plot of the heart rate observed after treatment under exercise 1, 2, and 3 versus baseline heart rate 1 with estimated linear relationships of two (or more) treatment alternatives to a patient based on his or her baseline biomarker value = . We assume a linear dependence between the potential responses and the biomarker. An example of Schwenke (1990 , Table 1) demonstrates our ideas. The RCT randomized 24 patients equally to three different exercise programs to improve the heart rate. After 8 weeks of training, each participant's heart rate was recorded after a 6 min run. Schwenke also provides data on the resting heart rate at baseline, which we consider as the biomarker of interest. Lower heart rates after the exercise are considered as more favorable. Figure 2 shows the data. It indicates a linear relationship between biomarker and outcome as well as a quantitative and qualitative interaction. This example is discussed in more detail in Section 4.
The article is organized as follows. In Section 2, the model for reconstructing the joint distribution of the PRs is presented. We assume normal marginal distributions of the responses under treatment 1 and 0, respectively. Further, the assumptions on the biomarker for reconstructing the joint distribution are presented. In Section 3, estimators for reconstructing the joint distribution and corresponding variability are derived using the maximum likelihood theory. Conditions for the existence of variation in ITEs are presented. In Section 4, a medical data example is presented which is analyzed by the developed methodology for reconstructing joint distributions. In Section 5, the limitations of the presented methodology are discussed and an outlook for further issues of research is given.

General setting
In accordance with Cox and Reid (2000, p. 20), we extend the potential responses as follows: with as constant mean parameters and ( 1 , 0 ) ′ as independent and identically distributed from a bivariate normal distribution with ( ) = 0. We denote the potential outcome of individual dependent on treatment , ∈ 0, 1 by . Substituting 1 and 0 in (1) yields where 0 can be interpreted as global mean, ( 1 − 0 ) as treatment effect, 0 as subject effect and ( 1 − 0 ) as subjecttreatment interaction. This can be interpreted as a structural model of heterogeneous treatment effects considering subjectspecific effects.
In an RCT based on the parallel-group design, we can approximate the potential responses 1 and 0 in (2) by a covariate and a linear relationship leading to 1 = 1 + 1 + 1 with 1 and 0 being identically and independently distributed additive errors. This in turn yields the following structural model of heterogeneous treatment effects approximating subject-specific effects by a covariate : We assume that the variables 1 , 0 , and follow a trivariate normal distribution given by where is the mean of , is the mean of , 2 is the variance of , and 2 is the variance of . The correlation between and is denoted by . The trivariate normal distribution implies that the ATE is given by Further, it also implies that the joint distribution of and is bivariate normal. The conditional distribution of given is hold. Generally, we cannot measure both PRs in a parallel-group RCT: we can only estimate the quantities , and 2 | for both treatment groups. Thus, we can calculate the expected value of the response conditional on the biomarker , [ | ], by a linear relationship described by the intercept and the slope . Based on the linear relationship, the expected difference [ 1 − 0 | ], can be estimated. If the difference [ 1 − 0 | ] is constant for every value of the biomarker then treatment-biomarker additivity will be present. If not, an interaction between the treatments and the biomarker is present.
The expected conditional ITE [ 1 − 0 | ] could of course be estimated via the standard linear model approach by including a biomarker-treatment interaction term and predicting the outcome under both treatments. However, this approach is not sufficient, as the variance of the ITEs is dependent on the correlation between the PRs. As this correlation is unknown, we propose a different approach which allows us estimating this parameter. Futhermore, having an estimate of this correlation, our approach allows obtaining individual response probabilities (and prediction intervals).
In order to determine the variance of the ITE the dependence structure of the two PRs has to be known. As a consequence, it is necessary to make some assumptions about the joint distribution of the two PRs. In the case of normally distributed PRs, Gadbury and Iyer (2000) and Gadbury, Iyer, and Allison (2001) rely on the trivariate normal distribution and derive bounds for the correlation coefficient 10 . The correlation coefficient 10 can be bounded due to the fact that the matrix Σ of the trivariate normal distribution has to be positive definite. Knowing the correlation coefficient 10 allows, in consequence, the bounding of the variance of the ITEs (unconditional or conditional on ) as outlined by Gadbury and Iyer (2000) and Gadbury et al. (2001). These bounds usually cover a wide range of 10 and are thus of limited practical value as already noticed by Lord (1955a).
The article examines a trivariate normal distribution as a model to evaluate the variation of ITEs when the PRs 1 and 0 and the biomarker are generated in a parallel-group RCT (Gadbury & Iyer, 2000;Gadbury et al., 2001;Lord, 1955a). With a point estimate of the correlation between the PRs, their joint distribution can be reconstructed and a point estimator for the variation of the ITEs can be derived. This allows the assessment of whether or not subject-treatment interactions are present. Additionally, it should be noted that knowledge about the correlation parameter 10 can be used for planning RCTs with less subjects if this correlation is positive. A more rapid and efficient drug development could be achieved. There is medical evidence that the assumption of independence is not reasonable when comparing treatments. In a series of articles, Cleophas et al. (Cleophas, 1996a(Cleophas, , 1996bCleophas and de Vogel, 1998;Cleophas, 2000) showed by using results of crossover RCTs that a positive correlation 10 between the responses 1 and 0 is observed if the compared treatments 1 and 0 have a similar pharmacological mode of action whereas a negative correlation 10 between the responses 1 and 0 is observed if the compared treatments 1 and 0 each have different pharmacological modes of action.

Reconstruction
We start by making assumptions about the PRs with regard to the mechanism of how the data are generated based on Cheng, Small, Tan, and Ten Have (2009, p. 21). First the PR of subject is assumed to be independent of the treatment allocation other than that of subject . This is commonly referred to as "stable unit treatment value assumption (SUTVA)." The second point would be that a subject enrolled in the RCT is an independent and identically distributed random sample from a well-defined population. The third assumption is that the allocation and are independent as guaranteed by random allocation. Aside from these standard assumptions we relate the biomarker to the ITE. Assume for the moment that we would know the joint distribution of the normally distributed marginals of the responses 1 and 0 as shown in Figure 3a. The diagonal line shows where the responses to the compared treatments are exactly equivalent. As described above, we make the reasonable assumption that the joint distribution of 1 and 0 follows a bivariate normal distribution.
Rotating Figure 3a by 45 • in a clock-wise manner yields Tukey's well-known "sum-difference plot," that is this plot is obtained by graphing the differences 1 − 0 on the sums 1 + 0 . Note that the y-axis shows the treatment effects for each F I G U R E 3 Marginal distributions 1 and 0 and "sum-difference plot." A linear relation between 1 + 0 and is assumed with only small measurement errors. The different gray tones denote the quantiles of . It can be seen that lighter values benefit more likely from treatment 0 while darker values benefit more likely from treatment 1 single subject. In our hypothetical example, the sum-difference plot corresponding to Figure 3a is shown in the plot of Figure 3b. The diagonal line in Figure 3a is now a horizontal line. Variation alongside the x-axis quantifies between-subject variation whereas variation along the y-axis quantifies within-subject variation. Since there is no independence between the sums 1 + 0 and differences 1 − 0 , subject-treatment interactions are present. The between-subject variation and the subject-treatment interactions are approximated in an RCT based on the parallel-group design by a biomarker . Since the biomarker and the sum 1 + 0 encode for the same information, namely the between-subject variation and the subject-treatment interactions, it is reasonable to assume a functional relationship between the biomarker and the sum 1 + 0 given by where and are constants and the error term is assumed to be independent of , formally ⟂ , independent of ( 1 + 0 ), formally ⟂ ( 1 + 0 ) and identically and independently distributed (iid) as | 1 + 0 ∼ N[0, 2 ] where 2 denotes the variance of . If the random variables ( 1 , 0 , ) ′ have a normal distribution then the linear combination in (9) will be normally distributed with mean value + ( 1 + 0 ) and variance 2 ( 2 1 + 2 0 + 2 10 1 0 ) + 2 . The use of the reconstruction variable leads to the following model of heterogeneous treatment effects: When looking at any model of heterogeneous treatment effects shown in (5) and (10), we can see the following restrictions.
• If 1 = 0 , then no interaction between the biomarker and the treatment is present. This might be due to the fact that there are indeed no heterogeneous treatment effects and thus no treatment-covariate interaction are present or due to the fact that the covariate simply cannot "recover" any heterogeneous treatment effects although indeed present. Alternatively, in case of 0 = 1 ≠ 0, a prognostic biomarker is present.
• Besides, if without loss of generality 1 = − 0 holds, then no between-subject variability is present and a "pure" subjecttreatment interaction is present. However, only knowledge about 0 is necessary since the in case of = 0 the biomarker effect is quantified by 0 or in case of = 1 the biomarker effect is quantified by − 0 .
The proof can be found in the Supporting Information. The relation between and 1 + 0 can be checked in the marginal distribution plots of and . If a linear relation between and both of them can be assumed, then the relation between and 1 + 0 is also linear. Now, the trivariate normal distribution is completely described by the reconstruction parameters = ( 1 , 0 , 2 1 , 2 0 , 10 , , , 2 ) ′ . The information about the parameter 10 , which is necessary for reconstructing the joint distribution of 1 and 0 , is "contained" in the correlations between and . These correlations are observable in parallel-group RCTs.

MAXIMUM LIKELIHOOD ESTIMATION
In the following, we propose maximum likelihood estimators of the parameters for model (11) which are based on data generated in a parallel-group RCT by ( , , ) with two treatment groups. We start with deriving maximum likelihood estimators of the parameters ( , 2 , , , 2 | ). Further, conditions where the estimation of the joint distribution of the PRs is possible are shown. For each subject either response 1 or 0 is observable. This results in a missing value problem with respect to the PRs 1 or 0 per subject which was first recognized by Lord (1955a) and Anderson (1957).
We start by introducing the data structure of an RCT based on the parallel-group design and corresponding notation. Let the treatment indicator be independently distributed from , ⟂ . Additionally, let = 1 + 0 be the sample size with 1 subjects allocated to treatment 1 and 0 to treatment 0. Ordering by the treatment indicator we get the following data structure (adapted from (Anderson, 1957, p. 202)) with missing data Then the maximum likelihood estimators for and 2 denoted bŷand̂2 are given bŷ Further, the maximum likelihood estimators for , , and 2 | denoted bŷ,̂, and̂2 | are given bŷ =̄−̂̄, Before proceeding to the estimation of the reconstruction parameters = ( , 2 , 10 , , , 2 ) ′ , we provide marginal param- Substituting the parameters of by the corresponding estimators for = ( , 2 , , , 2 | ) ′ provided in (12) Similar versions of the estimateŝand̂were derived by Lord (1955a), Lord (1955b) and Anderson (1957). For reconstruction of the joint distribution of the PRs, the reconstruction parameters = ( , 2 , 10 , , , 2 ) ′ have to be estimated.
The parameters for and 2 are already given in (17)-(18). The parameters ( 10 , , , 2 ) ′ are given by Substituting the parameters by the estimators for provided in (12)-(16) gives the estimators for denoted by The asymptotic variances (and covariances) of can be derived by using the multivariate delta method (Greene, 2003, chapter D.2.7) and can be found in the Supporting Information.
Besides, some restrictions appear in the parameter space and are given in terms of by ( ) or in terms of by ) .
The proof for deriving these restrictions can be found in Laubender (2014). In practice, checking the regression lines for both outcomes 0 and 1 on the biomarker will give an impression, whether there is an interaction effect. In order to understand the constraints 2 0 − 2 1 < 0 and 2 0| − 2 1| < 0, assume without loss of generality that we want to develop a cDx for a new treatment 1. A high response is more favorable and the variance 2 1 is higher than the variance 2 0 indicating the presence of subject-treatment interactions (see, e.g., Cox and Reid (2000, p. 21)). In this case, it is essential to see how the comparison treatment 0 performs under that cDx so that a treatment-biomarker interaction can be established. However, as outlined in Section 1, it is not sufficient to only consider the mean of the ITEs but the variance of the ITEs (either unconditional or conditional on ) should also be taken into account. The restriction of different sized residual variances 2 0| and 2 1| < 0 is a consequence of the structural model of heterogeneous treatment effects shown in (5) and (10).
When looking at the ITEs conditional on , three types of (qualitative) interactions can be distinguished and examples are shown in Figure 4. The scatter plots in Figure 4 show simulated responses of treatment 1 (gray crosses) and 0 (black circles) stratified by the cDx with corresponding regression lines superimposed. The first interaction is shown in the top of Figure 4 and is in accordance with the constraints 2 0 − 2 1 < 0 and 2 0| − 2 1| < 0. For high values of , we see that the conditional mean of treatment 1 is higher than that of treatment 0, and higher responses under treatment 1 compared to treatment 0 can be reached.
The second interaction is shown in the middle of Figure 4. In this case 2 0 − 2 1 > 0 and 2 0| − 2 1| < 0 hold so that the constraints are not fulfilled. For high values of we can now see that the conditional mean of treatment 0 is higher than under treatment 1, but under treatment 1 more higher responses are achieved under treatment 1 than under treatment 0. Thus, the cDx does not capture the variation of the responses 1 as accurate as the variation of the responses 0 . In this case the strict focusing on mean effects is misleading.
The third interaction is shown in the bottom of Figure 4. In this case 2 0 − 2 1 < 0 and 2 0| − 2 1| > 0 hold so that the constraints are not fulfilled. For very high values of , we see that the conditional mean of treatment 0 is higher than under treatment 1 and that higher responses under treatment 1 can be reached than under treatment 0. However, the lower the value of the cDx becomes the higher the values of response under treatment 0 are even though the mean of 1 conditional on is higher than the mean of 0 conditional on .
These constraints also imply the following well-known situations of linear regression modeling where no estimation of the reconstruction parameters and thus the joint distribution of the PRs 1 and 0 is possible: First, if 2 1 = 2 0 = 0 hold, then an uninformative biomarker is present. Nonetheless, there might be subject-treatment interactions present which cannot be modeled by the uninformative biomarker . Second, if 1 = 0 ≠ 0 and simultaneously 2 1| = 2 0| hold, then subject-treatment additivity is usually assumed in this situation. Third, if 1 = − 0 or 0 = − 1 and simultaneously 2 1| = 2 0| hold then subject-treatment interactions are present without any simultaneous subject effects. This finding is in accordance with the statement that subject-treatment interactions "cannot be estimated separately from variation among the units" (Cox & Reid, 2000, p. 20). An interaction without main effect for the biomarker is present, from the point of view of the linear regression model.

DATA EXAMPLE
A classic example of an RCT is taken from Schwenke (1990, table 1). The RCT randomized 24 patients equally to three types of exercise programs. The clinically relevant endpoint is the heart rate observed after treatment. A lower value is more favorable. The biomarker is the baseline heart rate (Schwenke, 1990, p. 444). Similar to Schwenke (1990, Figure 1), the scatter plot in Figure 2 shows the data. It shows the linear relationships between the heart rate observed after treatment and the baseline heart rate separately for each exercise group. Following Schwenke (1990, table 3, "Without Bonferroni Adjustment"), the contrast exercise 1 compared to exercise 3 (contrast A), and exercise 2 compared to exercise 3 (contrast B) are considered so that in the following exercise 3 is regarded as reference treatment ( = 0).
The crossing of the both lines of exercises 1 and 3 indicates a qualitative interaction for contrast A whereas the nearly noncrossing lines of exercise 2 and 3 pronounces a quantitative interaction for contrast B. A first check shows that the restrictions on the parameter space Θ ⋆ are fulfilled. It can be seen from Table 1 that for both contrasts A and B̂2 0 −̂2 1 < 0 and̂2 0| −̂2 1| < 0 hold. Table 1 also shows the estimates of the reconstruction parameters with asymptotic 95% confidence intervals. Exercise 1 is on average lower than exercise 3 with an ATEΔ = −5.75 and the corresponding 95% confidence interval includes 0. In contrast,

F I G U R E 4 Scatter plots of simulated responses
and a biomarker for treatment 1 (gray crosses) and for treatment 0 (black circles) with regression lines superimposed. All three plots show a qualitative interaction. Only the top plot shows data which fulfill the constraints 2 0 − 2 1 < 0 and 2 0| − 2 1| < 0 T A B L E 1 Estimates of selected parameters with corresponding asymptotic 95% confidence intervals (CIs). The lower and upper limits of the confidence intervals for some estimates are outside their range which is due to the fact that asymptotic estimators instead of exact estimators are used exercise 2 is on average higher than exercise 3 withΔ = 7.20, here the corresponding 95% confidence interval does not include 0. However, the ATE does not inform whether or not there are patients which will more likely profit from exercise 1 (or 2) or more likely from exercise 3. In order to get this information, it is necessary to know the correlation 10 for the contrasts A and B. The correlationŝ1 0 are positive for both contrasts. The resulting variation in ITEs is described bŷΔ and we can see that the corresponding 95% confidence intervals do not include 0 for both contrasts. This indicates the presence of subjecttreatment interactions. Based on the estimates of the reconstruction parameters , both bivariate normal densities are drawn for both contrasts in Figure 5 (top row). From this reconstructed joint distribution of 1 and 0 , we can quantify the subgroup of patients which will more likely benefit from exercise 1 (or 2) compared to exercise 3 and vice versa: In the upper wedge is the proportion of patients who will benefit from exercise 1 (or 2) whereas in the lower wedge is the proportion of patients who will benefit from exercise 3. This subgroup can be quantified by the probability that 1 > 0 hold, that is the probability that a higher (unfavorable) response is observed under exercise 1 (or 2) than under exercise 3. This probability is estimated by Φ[Δ∕̂Δ], where Φ[⋅] denotes the cumulative distribution function of the standard normal distribution. Note that this probability heavily depends on the correlation coefficient 10 . For contrast A, Φ[Δ∕̂Δ] is 0.35, that is on average 35 out of 100 randomly selected patients from a population will reach a higher unfavorable response under exercise 1 than under exercise 3 and thus should be treated by exercise 3. Vice versa, on average 65 out of 100 randomly selected patients will benefit from exercise 3 compared to exercise 1. For contrast B Φ[Δ∕̂Δ] = 0.77 is observed. A similar reasoning holds for contrast B. Although the ATE is in this case in clear favor for exercise 3, there are nonetheless 23 out of 100 randomly selected patients which will be harmed by exercise 3.
For clinical decision-making, it is of interest to not only have estimates with "integrated out" likeΔ and Φ[Δ∕̂Δ] but also to have estimates conditional on = . Such estimates are shown in Figure 5 (bottom row) where the treatment effects conditional on are shown with 95% confidence and prediction intervals. Here, probabilites that 1 > 0 conditional on holds are shown with asymptotic 95% confidence intervals. Details on these estimators can be found in Laubender (2014). For contrast A, it can be seen that there is a clear distinction between those patients benefitting from exercise 1 or from exercise 3 illustrating the qualitative nature of the treatment-biomarker interaction whereas for contrast B exercise 3 should be recommended for patients with higher values of the baseline heart rate.

F I G U R E 5
The left column of plots refers to contrast A whereas the right column of plots refers to contrast B. Top row: reconstructed joint distribution of heart rate after treatment 1 and 0 where the corresponding values of 1 and 0 are added by vertical and horizontal dotted lines. Further, the line of no difference between the two treatments is indicated by the diagonal; middle row: plot of the ITEs conditional on the baseline heart rate (solid line) with 95% confidence intervals (dashed lines) and 95% prediction intervals (dotted lines) and with reference line (dashed-dotted line) of no treatment difference between exercises 1 and 3 and exercises 2 and 3; bottom row: corresponding conditional response probabilities [ 1 > 0 | ] (solid line) with reference line (dashed-dotted line) of no treatment effect between exercises 1 and 3 and exercises 2 and 3

Conclusions
The knowledge of ITEs is helpful to select the best of two treatments for an individual subject. To improve the assessment of uncertainty of ITE estimates it is necessary to know the counterfactual correlation between the PRs 1 and 0 . In a parallelgroup RCT information of the marginal distributions of the single PR components is available, but no information of their joint distribution due to the "fundamental problem of causal inference" (Holland, 1986, p. 947).
We present a counterfactual approach for estimating the joint distribution of two normally distributed responses to two treatments. We assume a bivariate normal joint distribution for the PR. Additionally, we assume the presence of a normally distributed baseline biomarker which is functionally related to the sum 1 + 0 . Such a functional relationship is plausible since a biomarker and the sum 1 + 0 encode for the same information in an RCT : the between-subject variation. The biomarker as an indicator for the sum 1 + 0 , which creates constraints for the estimation of the bivariate PR distribution. To perform successfully a maximum likelihood estimation of the reconstruction parameters , the constraints 2 0 − 2 1 < 0 and 2 0| − 2 1| < 0 have to be fulfilled. The quantity 2 can be interpreted as proportion of the variance of explained by the biomarker . Similarly, the quantity 2 | can be regarded as unexplained or residual variance. Lord (1955a) formulated this problem in the early 1950s and described the estimation of parameters in the given specific incomplete data setting. He used a trivariate normal distribution. Our model combines this idea and ideas proposed by Gadbury and Iyer (2000). This paper derives bounds of the correlation 10 . For the first time, we provide a point estimate of the correlation between the individual potential responses 10 together with a confidence interval.
Based on the more developed theory, our approach facilitates a more-informed assessment of a biomarker's relevance for treatment selection than the classical approach of estimating an interaction effect between marker and treatment. As Huang, Gilbert, and Janes (2012) have demonstrated, a strong interaction coefficient is important for a biomarker to have value for treatment selection but is not useful for summarizing its predictive performance . The predictive performance depends on other coefficients in the risk model as well as the functional form of the model. Therefore, the interaction coefficient is not directly comparable between biomarkers (and models).
We Most of the methodological literature on treatment-selection markers discusses the issue in the context of a randomized trial. It also stresses that the statistical interaction between marker value and treatment assigned is the primary measure of marker performance (Buyse, 2007;Freidlin & Simon, 2005;Sargent, Conley, Allegra, & Collette, 2005;Simon, 2008;Simon, Paik, & Hayes, 2009). However, a strong interaction is important but not sufficient for adequate marker performance (Janes, 2011). Specifically, two markers can have the same interaction but very different performance. Huang et al. (2012) present an example where two biomarkers have the same numerical interaction coefficient estimates but show different capacity in terms of classifying a subject according to treatment effectiveness. Therefore, we see the response function as important information to be communicated to clinicians (see Figure 5).

Limitations and issues for future research
The main limitation of our approach consists in its distributional and structural assumptions: linear relationship between biomarker and PRs and the distributional assumptions, especially that the joint distribution of the PRs follows a bivariate normal PR distribution. We explore copula models as alternative to the bivariate normal PR distribution. They allow marginal normal distributions of the biomarker and the PRs 1 and 0 as well as a bivariate normal distributions of ( 1 , ) and of ( 0 , ) to be normal. Conditional copula models may have the potential to work out such distributions (Veraverbeke, Omelka, & Gijbels, 2011). It is important to stress that these observable and assessable relationships-even made normally distributed by a transformation-do not necessarily imply that the PRs follow a joint normal distribution.
It is well-known that univariate normality of 1 , 0 , and and bivariate normality of 1 and and of 0 and does not necessarily imply trivariate normality, however it is usually sufficient to evaluate these assumptions as outlined in standard text books on multivariate statistics: "in practice, …, the presence of joint nonnormality is likely to be detected quite often by methods directed at studying the marginal normality of the observations on each variable" (Gnanadesikan, 1977, p. 168) and "many types of nonnormality are often reflected in the marginal distributions and scatterplots. Moreover, for most practical work, one-dimensional and two-dimensional investigations are ordinarily sufficient. Fortunately, pathological data sets that are normal in lower-dimensional representations but nonnormal in higher dimensions are not frequently encountered in practice" (Johnson & Wichern, 1992, p. 153).
A similar critical aspect consists in the distributional assumptions of the reconstruction variable created by its linear relationship with the sum 1 + 0 . Here, it is an issue of research to explore the potential of proper scoring rules (see, e.g., Gneiting and Raftery (2007)) to assess the correct functional relationship between biomarker and the PRs 1 and 0 . It is also of interest to study models where the correlation between the PRs 1 and 0 depends on a biomarker value. This aspect is not studied so far and is an issue for future research. Finally, the model is developed for univariate biomarkers. It is an issue for future research to generalize to the use of multiple biomarkers simultaneously. Attempts have been made by Kaiser and Gadbury (2013). It is also of interest to apply the basic idea of this paper to other common clinical outcomes like binary endpoints and count data, especially to remove restrictions as introduced by Huang et al. (2012) to make the counterfactual model identifiable.
The extension of our work to multivariate measurements per subject, as in longitudinal data, are of interest. Another potential extension is the analysis of time-to-event data in the described context. The concept described herein could also be applied when 1 and 0 are (potentially right-censored) survival times that follow a log-normal distribution. However, for right-censored survival times, values would have to be imputed. This is a topic for future work. Log-normally distributed event times and handling right censoring with the EM-algorithm allows to translate the event data problem in the setting of our approach. The interesting question is, how many real data settings do fit in this structural straightjacket.

Summary
The presented work follows a line of ideas which started with Lord's paper on the calibration of scores, Anderson's likelihood approach to Lord's problem, and the work of Gadbury and colleagues to infer information on the correlation between the individual potential responses 01 . We complete this line of thinking by providing a structural model with a maximum likelihood estimation of 01 and related measures of uncertainty. Having presented the formal machinery, we demonstrate its application for individualized treatment decisions and derive decision tools for a straightforward interpretation of the formal outcomes. We stress the difference between our counterfactual and the usual approach, which combines exchangeability arguments and regression models with interaction effects between treatment and biomarker. Our approach exploits strong formal assumptions: trivariate normal distribution as well as linear relationships between biomarker and outcome. We present an example, which fit this formal setting and discuss strategies to assess these assumptions. The presented structural model offers approaches to causal sensitivity analyses of clinical trial data. There are options to generalize our model by using copula and monotone relationships between treatment outcome and biomarker. It is also of interest to explore Bayesian approaches to handle issues of non-identifiability, which easily raise for more general potential response models or when modeling non-Gaussian types of outcomes. We see our work as a proof of principle and as encouragement to explore the relevance of more general potential response models in clinical research.