Background
Heterogeneity of treatment effect (HTE), better called differential treatment effect, is variation in a measure of treatment effect on a scale for which it is mathematically possible that such variation be absent even if the treatment has a nonzero effect. The most commonly used scales for measuring HTE are
 the original scale for continuous response variables that are properly transformed, in which case the treatment effect is the (adjusted) difference in means
 odds ratios, for binary or ordinal outcomes
 hazard ratios, for timetoevent outcomes
Even though the more complex hazard ratio seems to be well accepted as a summary measure of treatment effect in a timetoevent randomized clinical trial (RCT), there is still a good deal of resistence to odds ratios (OR) from some clinical researchers. This resistence is difficult to understand, although it is clear that ORs are more difficult to understand than risk ratios (RR) or absolute risk reduction (ARR; risk difference). The reasons for choosing ORs are:
 ORs come directly from logistic models, and logistic models are as likely as any model to fit patterns leading to binary responses. This is primarily because the logistic model places no restrictions on the regression coefficients.
 ORs are capable of being constant over a range of baseline risk all the way from 0 to 1.0.
 In the multitude of forest plots present in journal articles depicting RCT results, the constancy of ORs over patient types is impressive.
 Unlike RRs, ORs are invariant to the choice of the “event” vs. the “no event”. If you interchange event and nonevent you would get the reciprocal of the original OR, but the RR would change arbitrarily.
RR and ARR are not capable of being constant. For example, a risk factor that doubles the risk of an outcome can only apply to a patient having a baseline risk of 0.5 or less, otherwise the risk with the risk factor would exceed 1.0. A treatment that lowers the risk of an outcome of 0.04 cannot apply to a patient with a baseline risk of 0.02.
The purposes of this article are to show how
 a rigorous test for HTE in the frequentist domain is done for an actual large RCT with a binary endpoint
 the difficulty in estimating and testing interaction effects
 the OR varies over patient types if one allows all the baseline variables to interact with treatment, whether the interactions are “significant” or not
 RR varies over patient types
 ARR varies over patient types
Difficulty of the Task
Establishing the existence of differential treatment effects in RCTs is difficult because RCTs are typically sized just large enough to detect an overall average treatment effect. In the ideal situation for showing evidence of HTE, a potentially interacting binary covariate yields optimum precision and power when the covariate has equal sample sizes for its two values. Even then, the precision of the interaction effect (a double difference) is four times as bad as the precision of the overall treatment effect. And as Andrew Gelman has discussed, the power is as if the sample size is divided by 16 because one typically wants to detect a differential effect that is only half as large as the overall detectable treatment effect. A related example will be shown below.
It is important to note that assessing treatment effect in an isolated subgroup defined by a categorical covariate does not establish HTE and results in unreliable estimates. Differential treatment effect must be demonstrated.
The GUSTOI Study and Basic Outcome Model
The GUSTOI Study was a 40,000 patient international trial comparing four treatment arms for 30day mortality following acute myocardial infarction: streptokinase (SK) with subcutaneous heparin, SK with intravenous heparin, accelerated dosing of tissue plasminogen activator (tPA) with intravenous heparin, and SK combined with tPA and intravenous heparin. The covariateadjusted analysis that forms the basis for this article was done by Steyerberg, who replaced a limited number of missing covariate values with singly imputed values. We will be primarily comparing the accelerated tPA arm (n=10,320) with the two combined SK arms that did not involve tPA (n=20,162).
require(Hmisc)
mu < markupSpecs$html # in Hmisc
knitrSet(lang='blogdown')
load(url('http://hbiostat.org/data/gusto.rda'))
gusto < upData(gusto, keep=Cs(day30, tx, age, Killip, sysbp, pulse, pmi, miloc))
## Input object size: 5241552 bytes; 29 variables 40830 observations
## Kept variables day30,tx,age,Killip,sysbp,pulse,pmi,miloc
## New object size: 1476752 bytes; 8 variables 40830 observations
html(describe(gusto), scroll=TRUE)
8 Variables 40830 Observations
day30
n  missing  distinct  Info  Sum  Mean  Gmd 

40830  0  2  0.195  2851  0.06983  0.1299 
Killip: Killip Class
n  missing  distinct 

40830  0  4 
Value I II III IV Frequency 34825 5141 551 313 Proportion 0.853 0.126 0.013 0.008
age
n  missing  distinct  Info  Mean  Gmd  .05  .10  .25  .50  .75  .90  .95 

40830  0  5575  1  60.9  13.62  40.89  44.69  52.09  61.52  69.86  76.25  79.45 
pulse: Heart Rate beats/min
n  missing  distinct  Info  Mean  Gmd  .05  .10  .25  .50  .75  .90  .95 

40830  0  166  0.999  75.41  19.55  50  55  62  73  86  98  107 
sysbp: Systolic Blood Pressure mmHg
n  missing  distinct  Info  Mean  Gmd  .05  .10  .25  .50  .75  .90  .95 

40830  0  209  0.999  129  26.56  92  100  112  130  144  160  170 
miloc: MI Location
n  missing  distinct 

40830  0  3 
Value Inferior Other Anterior Frequency 23495 1435 15900 Proportion 0.575 0.035 0.389
pmi: Previous MI
n  missing  distinct 

40830  0  2 
Value no yes Frequency 34104 6726 Proportion 0.835 0.165
tx: Tx in 3 groups
n  missing  distinct 

40830  0  3 
Value SK+tPA SK tPA Frequency 10320 20162 10348 Proportion 0.253 0.494 0.253
Base Risk Model
A simplified covariateadjusted risk model as developed by Steyerberg is fitted below. For simplicity when later interacting baseline variables with treatment, the important age by Killip class interaction is omitted. Pulse rate is modeled using a linear spline with a knot at 50 beats/minute.
require(rms)
options(prType='html')
dd < datadist(gusto); options(datadist='dd')
f < lrm(day30 ~ tx + age + Killip + pmin(sysbp, 120) +
lsp(pulse, 50) + pmi + miloc, data=gusto,
eps=0.005, maxit=30)
f
lrm(formula = day30 ~ tx + age + Killip + pmin(sysbp, 120) + lsp(pulse, 50) + pmi + miloc, data = gusto, eps = 0.005, maxit = 30)
Model Likelihood Ratio Test 
Discrimination Indexes 
Rank Discrim. Indexes 


Obs 40830  LR χ^{2} 4143.09  R^{2} 0.243  C 0.820 
0 37979  d.f. 12  g 1.418  D_{xy} 0.641 
1 2851  Pr(>χ^{2}) <0.0001  g_{r} 4.128  γ 0.641 
max ∂log L/∂β 3×10^{9}  g_{p} 0.081  τ_{a} 0.083  
Brier 0.055 
β  S.E.  Wald Z  Pr(>Z)  

Intercept  3.4692  0.7095  4.89  <0.0001 
tx=SK  0.0742  0.0513  1.45  0.1483 
tx=tPA  0.1358  0.0609  2.23  0.0258 
age  0.0796  0.0021  37.12  <0.0001 
Killip=II  0.6011  0.0511  11.76  <0.0001 
Killip=III  1.2092  0.1053  11.49  <0.0001 
Killip=IV  1.9292  0.1391  13.87  <0.0001 
sysbp  0.0391  0.0017  23.40  <0.0001 
pulse  0.0210  0.0141  1.48  0.1376 
pulse’  0.0406  0.0144  2.82  0.0048 
pmi=yes  0.4714  0.0486  9.70  <0.0001 
miloc=Other  0.3144  0.1162  2.70  0.0068 
miloc=Anterior  0.5441  0.0444  12.27  <0.0001 
SemiSaturated Model
The only way that ORs for treatment can vary in a logistic model is for there to be one or more interactions between baseline covariates and treatment. Here we allow for all interactions with treatment, which will allow ORs, RRs, and ARRs to vary as the data dictate, assuming that threeway interactions are not needed, e.g., the infarct locationspecific treatment effect does not depend on age.
g < lrm(day30 ~ tx * (age + Killip + pmin(sysbp, 120) +
lsp(pulse, 50) + pmi + miloc), data=gusto,
eps=0.005, maxit=60)
print(g, coefs=FALSE)
lrm(formula = day30 ~ tx * (age + Killip + pmin(sysbp, 120) + lsp(pulse, 50) + pmi + miloc), data = gusto, eps = 0.005, maxit = 60)
Model Likelihood Ratio Test 
Discrimination Indexes 
Rank Discrim. Indexes 


Obs 40830  LR χ^{2} 4159.64  R^{2} 0.244  C 0.821 
0 37979  d.f. 32  g 1.428  D_{xy} 0.641 
1 2851  Pr(>χ^{2}) <0.0001  g_{r} 4.170  γ 0.642 
max ∂log L/∂β 2×10^{10}  g_{p} 0.082  τ_{a} 0.083  
Brier 0.055 
The likelihood ratio test is the goldstandard frequentist method for testing for added value of the treatment interactions, i.e., whether the treatment ORs are constant over covariate settings. By using a chunk test to test whether any of the baseline variables interact with treatment, this method has a perfect multiplicity adjustment in the frequentist sense.
lrtest(f, g)
Model 1: day30 ~ tx + age + Killip + pmin(sysbp, 120) + lsp(pulse, 50) +
pmi + miloc
Model 2: day30 ~ tx * (age + Killip + pmin(sysbp, 120) + lsp(pulse, 50) +
pmi + miloc)
L.R. Chisq d.f. P
16.552157 20.000000 0.681833
We see no evidence to suggest HTE for tPA or SK effects. We can also use AIC to assess whether allowing for interactions will likely result in better patientspecific outcome predictions.
c(AIC(f), AIC(g))
[1] 16558.31 16581.76
The simpler model has a lower (better) AIC, indicating that treatment interactions are not strong enough to overcome the overfitting that allowing for them would entail. In other words, the interactions are likely more noise than signal.
We can also compare the two models by computing the relative explained variation on the risk scale as detailed here.
rexv < var(predict(f, type='fitted')) /
var(predict(g, type='fitted'))
rexv
[1] 0.9959584
From this we see that even without penalizing for overfitting, all the treatment interactions account for only 0.004 of the predictive information. The nointeraction logitadditive model that assumes constancy of treatment ORs is at least 0.996 adequate on a scale from 0 to 1.
Back to the Difficulty of the Task
Covariatespecific treatment effects are combinations of main effects and interaction effects, and one needs evidence that the interaction effect is nonzero before covariatespecific relative efficacy is interesting. As an example of how difficult it is to estimate differential treatment effects (here, double differences on the logit scale, or ratios of odds ratios), let’s temporarily refit the logistic model including only a single treatment interaction—with location of infarct. Consider just inferior vs. anterior infarcts, which are the dominant categories and are not too imbalanced.
We estimate the double difference in log odds for tPA vs. SK in anterior infarcts minus tPA vs. SK in inferior infarcts. Then we get the frequencyweighted overall treatment effect in these two groups, which is similar to dropping the interaction from the model.
i < lrm(day30 ~ tx * miloc + age + Killip + pmin(sysbp, 120) +
lsp(pulse, 50) + pmi, data=gusto,
eps=0.005, maxit=30)
contrast(i, list(tx='tPA', miloc='Inferior'),
list(tx='SK', miloc='Inferior'),
list(tx='tPA', miloc='Anterior'),
list(tx='SK', miloc='Anterior'))
Contrast S.E. Lower Upper Z Pr(>z)
11 0.1546134 0.1081337 0.05732489 0.3665516 1.43 0.1528
Confidence intervals are 0.95 individual intervals
w < with(gusto, c(sum(miloc=='Inferior'), sum(miloc=='Anterior')))
contrast(i, list(tx='tPA', miloc=c('Inferior', 'Anterior')),
list(tx='SK', miloc=c('Inferior', 'Anterior')),
type='average', weights=w)
Contrast S.E. Lower Upper Z Pr(>z) var
1 0.1838079 0.05588118 0.293333 0.07428276 3.29 0.001 0.003122706
Confidence intervals are 0.95 individual intervals
The first contrast (double difference; differential treatment effect) has a standard error that is twice that of the second contrast (treatment effect averaged over the two infarct locations), and only the latter overall treatment effect has evidence for being nonzero (i.e., a benefit of tPA).
Penalization
Bayesian modeling of HTE uses hierarchical random effects with shrinkage of interaction effects towards zero, e.g., shrinkage of “subgroup” effects towards the grand mean treatment effect. Borrowing information in this fashion is an optimal way to get covariatespecific treatment effect estimates. The frequentist counterpart that handles mixtures of continuous and categorical baseline variables is penalized maximum likelihood estimation. We will leave the main effects unpenalized, favoring an additive (in the log odds) model with many fewer parameters but allowing the completely unpenalized model to also have a chance at winning, and find the penalty for interaction terms that optimizes the effective AIC, i.e., the AIC computed using the effective degrees of freedom. Note that the pentrace
function uses the likelihood ratio χ^{2} scale for AIC unlike the traditional AIC formula, i.e., here bigger AIC is better.
h < update(g, x=TRUE, y=TRUE)
p < pentrace(h,
list(simple=0,
interaction=1000*
c(1, 5, 10, 15, 20, 30, 40, 100, 200, 500, 1000)),
maxit=30)
p
Best penalty:
simple interaction df
0 1e+06 12.03211
simple interaction df aic bic aic.c
0 1000 20.33720 4111.002 3935.753 4110.981
0 5000 15.75213 4115.700 3979.961 4115.687
0 10000 14.32236 4117.036 3993.618 4117.025
0 15000 13.69510 4117.613 3999.600 4117.603
0 20000 13.33778 4117.937 4003.003 4117.928
0 30000 12.94300 4118.289 4006.757 4118.280
0 40000 12.72881 4118.477 4008.791 4118.469
0 100000 12.30905 4118.837 4012.768 4118.829
0 200000 12.15779 4118.963 4014.197 4118.955
0 500000 12.06394 4119.040 4015.083 4119.032
0 1000000 12.03211 4119.066 4015.383 4119.058
h < update(g, penalty=list(simple=0, interaction=30000),
maxit=30, eps=0.005)
effective.df(h)
Original and Effective Degrees of Freedom
Original Penalized
All 32 12.94
Simple Terms 11 11.00
Interaction or Nonlinear 21 1.94
Nonlinear 3 1.05
Interaction 20 0.94
Nonlinear Interaction 2 0.05
The following plots show first the main effect coefficients (with interaction effects penalized vs. not penalized), and second the interaction effect coefficients. Note that main effect coefficients can change just because of penalizing interactions, just as removing interactions from models can greatly change main effects.
ia < grep('\\*', names(coef(g)))
x < coef(g)[c(1,ia)]; y < coef(h)[c(1,ia)]
r < range(c(x, y))
plot(x, y,
xlab='Unpenalized Coefficients',
ylab='Penalized Coefficients',
main='Main Effects', xlim=r, ylim=r)
abline(a=0, b=1, col=gray(.85))
x < coef(g)[ia]; y < coef(h)[ia]
r < range(c(x, y))
plot(x, y,
xlab='Unpenalized Coefficients',
ylab='Penalized Coefficients',
main='Interaction Parameters', xlim=r, ylim=r)
abline(a=0, b=1, col=gray(.85))
Of the penalties studied, the best penalty by both AIC and corrected AIC is 1e+06 which effectively makes the model estimate 12.03 parameters instead of the original 12 parameters. The optimum penalty is really infinity, consistent the the likelihood ratio test for presence of any interactions with treatment. To give the benefit of the doubt we will use a penalty of 30000 to achieve an effective d.f. for interactions of 0.94.
Using the Penalized SemiSaturated Model to Estimate Distributions of Various Effects
Even though there is no statistical evidence for HTE (i.e., interactions with treatment on an appropriate scale), we will still allow for all interactions, to give the benefit of the doubt. But we are penalizing the interaction parameters closer to what will crossvalidate.
Let’s start with what I believe is the best way to communicate results to an individual patient: estimating the absolute risk of mortality under two treatment options, without hinting how the patient should summarize the two estimates (absolute vs. absolute difference, etc.). The baseline covariate values for one patient are listed below.
d < data.frame(tx=c('SK', 'tPA'), Killip='II', age=55, pulse=40, sysbp=100,
miloc='Inferior', pmi='no')
p < predict(h, d, type='fitted')
names(p) < c('SK', 'tPA')
p
SK tPA
0.04043122 0.03302851
Now turn to estimating treatment effects for all the patients in the RCT. To estimate patientspecific treatment effects on various scales, we estimate risks of dying within 30 days for each patient with their own covariate settings (baseline covariate values). This is done first with treatment set to SK, no matter which treatment the patient actually received. Then predictions are done setting treatment to tPA.
d < gusto
d$tx < 'SK'
psk < predict(h, d, type='fitted')
d$tx < 'tPA'
ptpa < predict(h, d, type='fitted')
k < c(5,10,15,20,50,100,200,400,800,1600,3200,20000)
pl < c(0, quantile(psk, 0.99))
j < pmax(psk, ptpa) <= pl
ggfreqScatter(psk[j], ptpa[j], cuts=k,
xlab='P(deathSK)', ylab='P(deathtPA)')
Here is the distribution of estimated ARR across all patient types.
hist(psk  ptpa, xlim=c(.05, .1), nclass=200, main=NULL)
And the distribution of estimated RR.
hist(ptpa / psk, nclass=200, xlim=c(.75, 1), main=NULL)
And finally the distribution of estimated ORs.
or < exp(qlogis(ptpa)  qlogis(psk))
hist(or, nclass=200, xlim=c(.75, 1), main=NULL)
The following scatterplot shows the relationship between estimated RRs and estimate ORs.
ggfreqScatter(or, ptpa / psk, cuts=k, xlab='OR', ylab='RR')
See that the variation in absolute or relative risk reduction is large, but adjusted odds ratios for tPA vs. SK have a very narrow distribution. If we relied on the formal interaction test, this distribution would have a single value.
To see how the various measures depend on baseline risk, see the following graphs.
xl < 'Risk Under Control Treatment (SK)'
j < psk <= pl
ggfreqScatter(psk[j], psk[j]  ptpa[j], cuts=k, xlab=xl, ylab='ARR')
ggfreqScatter(psk[j], ptpa[j] / psk[j], cuts=k, xlab=xl, ylab='RR')
ggfreqScatter(psk[j], or[j], cuts=k, xlab=xl, ylab='OR')
Discussion
Add your comments, suggestions, and criticisms on datamethods.org