Test for interaction - Columbia University

Report
Are you looking for the right
interactions?
Part 2: Statistically testing for interaction effects with
dichotomous outcome variables
Melanie M. Wall
Departments of Psychiatry and Biostatistics
New York State Psychiatric Institute and Mailman School of Public Health
Columbia University
[email protected]
Joint Presentation with Sharon Schwartz (Part 1)
Department of Epidemiology
Mailman School of Public Health
Columbia University
1
Data from Brown and Harris (1978) – 2X2X2 Table
Vulnerability
Lack of
Intimacy
Exposure
Stress Event
No
Yes
Risk of
Depression
No
No
191
2
0.010
P00
0.102
RD = 0.092
(0.027,0.157)
P01
0.032
P10
0.316
RD=0.284
(0.170,0.397)
P11
Yes
Effect of Stress given No Vulnerability ->
Yes
No
Yes
Effect of Stress given Vulnerability ->
Outcome (Depression)
79
9
OR = 10.9 (2.3, 51.5)
RR = 9.9 (2.2, 44.7)
60
2
52
24
OR = 13.8 (3.1,61.4)
RR = 9.8 (2.4, 39.8)
OR = Odds Ratio (95% Confidence Interval) <-compare to 1
RR = Risk Ratio (95% Confidence Interval) <-compare to 1
RD = Risk Difference (95% Confidence Interval) <-compare to 0
2
Does Vulnerability Modify the Effect
of Stress on Depression?
• On the multiplicative Odds Ratio scale, is 10.9 sig different from 13.8?
– Test whether the ratio of the odds ratios
(i.e. 13.8/10.9 = 1.27) is significantly different from 1.
• On the multiplicative Risk Ratio scale, is 9.9 sig different from 9.8?
– Test whether the ratio of the risk ratios
(i.e. 9.8/9.9 = 0.99) is significantly different from 1.
• On the additive Risk Difference scale, is 0.092 sig different from 0.284?
– Test whether the difference in the risk differences
(i.e. 0.28-0.09 = 0.19) is significantly different from 0.
Rothman calls this difference in the risk differences the “interaction
contrast (IC)”
IC = (P11 - P10) – (P01 - P00)
3
Comparing stress effects across vulnerability groups
Different conclusions on multiplicative vs additive scale
Comparing Risk Ratios (RR)
Comparing Risk Differences (RD)
60
Risk Ratio of Stress Event on Depression
10
20
30
40
50
OR=
OR =
RR =
9.9
9.8
0
13.8
0
10.9
RR=
no
Lack of Intimacy
yes
no
Lack of Intimacy
yes
Risk Difference of Stress Event on Depression
0.0
0.1
0.2
0.3
60
Odds Ratio of Stress Event on Depression
10
20
30
40
50
0.4
Comparing Odds Ratios (OR)
RD =
0.284
RD=
0.092
no
Lack of Intimacy
yes
95% confidence intervals for Odds Ratios overlap
-> no statistically significant multiplicative interaction OR scale
95% confidence intervals for Risk Ratios overlap
-> no statistically significant multiplicative interaction RR scale
95% confidence intervals for Risk Differences do not overlap
-> statistically significant additive interaction
In general, it is possible to reach
different conclusions on the two
different multiplicative scales
“distributional interaction”
(Campbell, Gatto, Schwartz 2005)
4
Modeling Probabilities
1.0
Binomial modeling with logit, log, or linear link
Binomial model
logit link
0.8
linear link
Probability of Depression
0.4
0.6
log link
Lack Intimacy
0.0
0.2
No Lack of Intimacy
no
Stress Event
yes
5
Test for multiplicative interaction on the OR scaleLogistic Regression with a cross-product
IN SAS:
proc logistic data = brownharris descending;
model depressn = stressevent
lack_intimacy
stressevent*lack_intimacy;
oddsratio stressevent / at(lack_intimacy = 0 1);
oddsratio lack_intimacy / at(stressevent = 0 1);
run;
Analysis of Maximum Likelihood Estimates
Parameter
Intercept
stressevent
lack_intimacy
stresseve*lack_intim
DF
1
1
1
1
Standard
Error
0.7108
0.7931
1.0109
1.0984
Estimate
-4.5591
2.3869
1.1579
0.2411
Wald
Chi-Square
41.1409
9.0576
1.3120
0.0482
Pr > ChiSq
<.0001
0.0026
0.2520
0.8262
exp(.2411) = 1.27 =
Ratio of Odds ratios =13.846/10.880
Not significantly different from 1
Wald Confidence Interval for Odds Ratios
Label
Estimate
stressevent at lack_intimacy=0
10.880
stressevent at lack_intimacy=1
13.846
lack_intimacy at stressevent=0
3.183
lack_intimacy at stressevent=1
4.051
95% Confidence
2.299
3.122
0.439
1.745
Limits
51.486
61.408
23.086
9.405
“multiplicative
interaction” on
OR scale is not
significant
6
Test for interaction: Are the lines Parallel?
Probability scale
0
Log Odds scale
Lack Intimacy
0.4
Lack Intimacy
No Lack of Intimacy
Test for whether lines are
parallel on probability scale is
same as H0: IC = 0.
Need to construct a statistical
test for IC = P11-P10-P01+P00
P11
-6
no
Stressful event
yes
P01
P10
0.0
-5
Cross product term in logistic
regression is magnitude of deviation
of these lines from being parallel…
p-value = 0.8262 -> cannot reject that
lines on logit scale are parallel
Thus, no statistically significant
multiplicative interaction on OR scale
0.1
Log Odds of Depression
-4
-3
Probability of Depression
0.2
-2
0.3
-1
No Lack of Intimacy
P00
no
Stress event
yes
7
The Problem with Comparing Statistical Significance
of Effects Across Groups
• Don’t fall into the trap of concluding there must be effect
modification because one association was statistically significant
while the other one was not.
• In other words, just because a significant effect is found in one
group and not in the other, does NOT mean the effects are
necessarily different in the two groups (regardless of whether
we use OR, RR, or RD).
• Remember, statistical significance is not only a function of the
effect (OR, RR, or RD) but also the sample size and the baseline
risk. Both of these can differ across groups.
• McKee and Vilhjalmsson (1986) point out that Brown and Harris
(1978) wrongfully applied this logic to conclude there was
statistical evidence of effect modification (fortunately there
conclusion was correct despite an incorrect statistical test )
8
Testing for additive interaction on the probability scale
Strategy #1: Use linear binomial regression with a cross-product
Risk = b0 + b1 * STRESS + b2 * LACKINT + b3*STRESS*LACKINT
NOTE: b3 = IC
link=identity dist=binomial tells SAS to do linear binomial
regression. lrci outputs likelihood ratio (profile likelihood)
confidence intervals.
IN SAS:
proc genmod data = individual descending;
model depressn = stressevent lack_intimacy stressevent*lack_intimacy/
link = identity dist = binomial lrci;
estimate 'RD of stressevent when intimacy = 0' stressevent 1;
estimate 'RD of stressevent when intimacy = 1' stressevent 1 stressevent*lack_intimacy 1;
run;
Analysis Of Maximum Likelihood Parameter Estimates
Standard
Parameter
DF
Estimate
Error
Intercept
1
0.0104
0.0073
stressevent
1
0.0919
0.0331
lack_intimacy
1
0.0219
0.0236
stresseve*lack_intim
1
0.1916
0.0667
Likelihood Ratio
95% Confidence
Limits
0.0017
0.0317
0.0368
0.1675
-0.0139
0.0870
0.0588
0.3219
Wald
Chi-Square
2.02
7.70
0.86
8.26
Pr>ChiSq
0.1551
0.0055
0.3534
0.0040
Contrast Estimate Results
Label
RD of stressevent when intimacy = 0
RD of stressevent when intimacy = 1
Mean
Estimate
0.0919
0.2835
Mean
Confidence Limits
0.0270
0.1568
0.1701
0.3969
Standard
Error
0.0331
0.0578
Interaction is statistically significant “additive interaction”.
Reject H0: IC = 0, i.e. Reject parallel lines on probability 9scale
Different strategies for statistically testing
additive interactions on the probability scale
The IC is the Difference of Risk Differences. IC = (P11 - P10) – (P01 - P00) = P11-P10-P01+P00
From Cheung (2007) “Now that many commercially available statistical packages have the capacity to fit log binomial and linear
binomial regression models, ‘there is no longer any good justification for fitting logistic regression models and estimating odds
ratios’ when the odds ratio is not of scientific interest” Inside quote from Spiegelman and Herzmark (2005).
1.
Fit a linear binomial regression Risk = b0 + b1 * EXPO + b2 * VULN + b3*EXPO*VULN. The b3 = IC and so a
test for coefficient b3 is a test for IC. Can be implemented directly in PROC GENMOD. PROS: Contrast of
interest is directly estimated and tested and covariates easily included CONS: Linear model for probabilities
can be greater than 1 and less than 0 and thus maximum likelihood estimation can be a problem. Waldtype confidence intervals can have poor coverage (Storer et al 1983), better to use profile likelihood
confidence intervals.
2.
Fit a logistic regression log(Risk/(1-Risk)) = b0 + b1 * EXPO + b2 * VULN + b3*EXPO*VULN, then backtransform parameters to the probability scale to calculate IC. Can be implemented directly in PROC
NLMIXED. PROS: logistic model more computationally stable since smooth decrease/increase to 0 and 1.
CONS: back-transforming can be tricky for estimator and standard errors particularly in presence of
covariates. Covariate adjusted probabilities are obtained from average marginal predictions in the fitted
logistic regression model (Greenland 2004). Homogeneity of covariate effects on odds ratio scale is not the
same as homogeneity on risk difference scale and this may imply misspecification (Kalilani and Atashili
2006; Skrondal 2003).
3.
Instead of IC, use IC ratio. Divide the IC by P00 and get a contrast of risk ratios:
IC Ratio = P11/P00 -P10/P00 -P01/P00+P00/P00 = RR(11) – RR(10) – RR(01) + 1 called the
Relative Excess Risk due to Interaction (RERI).
Many papers on inference for RERI
10
Test for additive interaction on the probability scale
Strategy #2: Use logistic regression and back-transform estimates
to form contrasts on the probability scale
PROC NLMIXED DATA=individual;
***logistic regression model is;
odds = exp(b0 +b1*stressevent + b2*lack_intimacy + b3*stressevent*lack_intimacy);
pi = odds/(1+odds);
MODEL depressn~BINARY(pi);
estimate 'p00' exp(b0)/(1+exp(b0));
estimate 'p10' exp(b0+b1)/(1+exp(b0+b1));
estimate 'p01' exp(b0+b2)/(1+exp(b0+b2));
estimate 'p11' exp(b0+b1+b2+b3)/(1+exp(b0+b1+ b2+b3));
estimate 'p11-p10' exp(b0+b1+b2+b3)/(1+exp(b0+b1+ b2+b3))- exp(b0+b1)/(1+exp(b0+b1));
estimate 'p01-p00' exp(b0+b2)/(1+exp(b0+b2)) - exp(b0)/(1+exp(b0));
estimate 'IC= interaction contrast = p11-p10 - p01 + p00'
exp(b0+b1+b2+b3)/(1+exp(b0+b1+ b2+b3)) - exp(b0+b1)/(1+exp(b0+b1)) - exp(b0+b2)/(1+exp(b0+b2)) +
exp(b0)/(1+exp(b0));
estimate 'ICR= RERI using RR = p11/p00 - p10/p00 - p01/p00 + 1'
exp(b0+b1+b2+b3)/(1+exp(b0+b1+ b2+b3))/ (exp(b0)/(1+exp(b0)))
- exp(b0+b1)/(1+exp(b0+b1))/ (exp(b0)/(1+exp(b0)))
- exp(b0+b2)/(1+exp(b0+b2)) / (exp(b0)/(1+exp(b0))) + 1;
estimate 'ICR= RERI using OR' exp(b1+b2+b3) - exp(b1) - exp(b2) +1;
RUN;
11
Strategy #2 Output from NLMIXED
Parameter Estimates
Parameter Estimate
b0
b1
b2
b3
-4.5591
2.3869
1.1579
0.2411
Standard
Error
0.7108
0.7931
1.0109
1.0984
Label
Estimate
p00
0.01036
p10
0.1023
p01
0.03226
p11
0.3158
p11-p10
0.2135
p01-p00
0.02190
IC =p11-p10-p01+p00 0.1916
RERI using RR
RERI using OR
18.4915
31.0138
DF
419
419
419
419
t Value
Pr > |t|
Alpha
-6.41
<.0001
0.05
3.01
0.0028
0.05
1.15
0.2527
0.05
0.22
0.8264
0.05
Additional Estimates
Lower
Upper
-5.9563
0.8280
-0.8291
-1.9180
-3.1620
3.9458
3.1450
2.4002
Standard
Error
0.00728
0.03230
0.02244
0.05332
0.06234
0.02359
0.06666
DF
419
419
419
419
419
419
419
t Value
1.42
3.17
1.44
5.92
3.43
0.93
2.87
Pr > |t|
0.1559
0.0017
0.1513
<.0001
0.0007
0.3539
0.0042
Lower
-0.00397
0.03878
-0.01185
0.2110
0.09098
-0.02448
0.06060
Upper
0.0246
0.1658
0.0763
0.4206
0.3361
0.0682
0.3226
13.8661
24.3583
419
419
1.33
1.27
0.1831
0.2036
-8.7644
-16.8659
45.7473
78.8936
Gradient
-0.00002
-0.00003
2.705E-6
-0.00001
IC estimator same as
strategy #1, but
slightly different s.e.,
p-value, 95% conf
interval
•
The IC estimator is same as before (slide 9) but slightly different s.e., p-value and 95% confidence interval
– still conclude there is a significant additive interaction.
•
Results for RERI (using RR and OR) indicate that there is NOT a significant additive interaction. This
conflicts with the conclusion that the IC is highly significant. The cause of the discrepancy is related to
estimation of standard errors and confidence intervals. Literature indicates Wald-type confidence intervals
perform poorly for RERI (Hosmer and Lemeshow 1992; Assman et al 1996).
•
Proc NLMIXED uses Delta method to obtain standard errors of back-transformed parameters and Waldtype confidence intervals, i.e. (estimate) +- 1.96*(standard error) . Possible to obtain profile likelihood
confidence intervals using a separate macro (Richardson and Kaufman 2009) or PROC NLP (nonlinear
programming) (Kuss et al 2010). Also possible to bootstrap (Assman et al 1996 and Nie et al 2010) or
incorporate prior information (Chu et al 2011)
12
Conclusion
• The appropriate scale on which to assess interaction effects with
dichotomous outcomes has been a controversial topic in epidemiology for
years, but awareness of this controversy is not yet wide spread enough.
• This would not be a problem if the status quo for examining effect
modification (i.e. testing interaction effects in logistic regression) was
actually the “RIGHT” thing to do, but, persuasive arguments have been
made from the sufficient cause framework that the additive probability
scale (not the multiplicative odds ratio scale) should be used to assess the
presence of synergistic effects (Darroch 1997, Rothman and Greenland
1998, Schwartz 2006, Vanderwheel and Robins 2007,2008)
• There are now straightforward ways within existing software to estimate
and test the statistical significance of additive interaction effects.
• Additional work is needed getting the word out that effect modification
should not (just) be looked at using Odds Ratios.
13
Appendix Material
14
Reading in the Brown and Harris data into SAS
data a;
input stressevent lack_intimacy depressn count;
***When entering the counts, it is necessary to subtract the cases
from the denominator to get the non-cases, e.g. 9/88 becomes 9 cases and 79 (88-9) noncases;
cards;
0 0 0 191
0012
0 1 0 60
0112
1 0 0 79
1019
1 1 0 52
1 1 1 24
;
****This data step creates 419 records corresponding to the counts above, basically one record for each of
individuals in the study;
data individual; set a;
do i = 1 to count;
output;
end;
drop i;
run;
15
Appendix Appendix Appendix Appendix Appendix Appendix Appendix Appendix
SAS (Proc GENMOD) code for fitting
logit, log, and linear binomial models
***logistic binomial regression;
proc genmod data = individual descending;
model depressn = stressevent lack_intimacy stressevent*lack_intimacy/
link = logit dist = binomial lrci type3;
estimate 'OR of stressevent when intimacy = 0' stressevent 1/exp;
estimate 'OR of stressevent when intimacy = 1' stressevent 1 stressevent*lack_intimacy 1/exp;
estimate 'OR interaction contrast' stressevent*lack_intimacy 1/exp;
run;
***log binomial regression;
proc genmod data = individual descending;
model depressn = stressevent lack_intimacy stressevent*lack_intimacy/
link = log dist = binomial lrci type3;
estimate 'RR of stressevent when intimacy = 0' stressevent 1/exp;
estimate 'RR of stressevent when intimacy = 1' stressevent 1 stressevent*lack_intimacy 1/exp;
estimate 'RR interaction contrast' stressevent*lack_intimacy 1/exp;
run;
***linear binomial regression;
proc genmod data = individual descending;
model depressn = stressevent lack_intimacy stressevent*lack_intimacy/
link = identity dist = binomial lrci type3;
estimate 'RD of stressevent when intimacy = 0' stressevent 1;
estimate 'RD of stressevent when intimacy = 1' stressevent 1 stressevent*lack_intimacy 1;
estimate 'RD interaction contrast' stressevent*lack_intimacy 1;
run;
16
Appendix Appendix Appendix Appendix Appendix Appendix Appendix Appendix
Partial output from SAS Proc GENMOD code on
previous slide
Logistic binomial regression
Analysis Of Maximum Likelihood Parameter Estimates
Parameter
Intercept
stressevent
lack_intimacy
stresseve*lack_intim
DF
Estimate
Standard
Error
1
1
1
1
-4.5591
2.3869
1.1579
0.2411
0.7108
0.7931
1.0109
1.0984
Likelihood Ratio
95% Confidence
Limits
-6.3576
1.0038
-0.9791
-2.0342
-3.4207
4.2824
3.2951
2.5246
Wald
Chi-Square
Pr > ChiSq
41.14
9.06
1.31
0.05
<.0001
0.0026
0.2520
0.8263
Wald
Chi-Square
Pr > ChiSq
42.20
8.82
1.32
0.00
<.0001
0.0030
0.2510
0.9938
Wald
Chi-Square
Pr > ChiSq
2.02
7.70
0.86
8.26
0.1551
0.0055
0.3534
0.0040
Log binomial regression
Analysis Of Maximum Likelihood Parameter Estimates
Parameter
Intercept
stressevent
lack_intimacy
stresseve*lack_intim
DF
Estimate
Standard
Error
1
1
1
1
-4.5695
2.2894
1.1356
-0.0081
0.7034
0.7711
0.9893
1.0521
Likelihood Ratio
95% Confidence
Limits
-6.3593
0.9596
-0.9675
-2.2075
-3.4529
4.1567
3.2388
2.1996
Linear binomial regression
Analysis Of Maximum Likelihood Parameter Estimates
Parameter
Intercept
stressevent
lack_intimacy
stresseve*lack_intim
DF
Estimate
Standard
Error
1
1
1
1
0.0104
0.0919
0.0219
0.1916
0.0073
0.0331
0.0236
0.0667
Likelihood Ratio
95% Confidence
Limits
0.0017
0.0368
-0.0139
0.0588
0.0317
0.1675
0.0870
0.3219
17
Appendix Appendix Appendix Appendix Appendix Appendix Appendix Appendix
Data from Brown and Harris (1978) – 2X2X2 Table
Exposure
Stress Event
No
Vulnerability
Lack of
Intimacy
No
Yes
Effect of Vulnerability given No Stress->
Yes
No
Yes
Effect of Vulnerability given Stress ->
Outcome (Depression)
No
191
60
Yes
2
2
OR = 3.2 (0.4,23.1)
X2 = 1.46 p-value = 0.22
79
52
9
24
OR = 4.1 (1.7,9.4)
X2 = 11.57 p-value < 0.01
Risk of Depression
0.010
0.032
RD = 0.022
(-0.024,0.068)
0.102
0.316
RD=0.214
(0.091,0.336)
Brown and Harris Concluded that since the Chi-square test in the No Stress group was NOT
significant while the Chi-square test in the Stress group was significant that this was statistical
evidence of an interaction effect. This was an incorrect use of the Chi-square tests. Need to
compare the effects themselves as described.
18
Appendix Appendix Appendix Appendix Appendix Appendix Appendix Appendix
References
•
Andersson T, Alfredsson L, Kallberg H, Zdravkovic S (2005) Calculating measures of biological interaction, European Journal of Epidemiology 20: 575-579.
•
Assmann SF, Hosmer DW, Lemeshow S, et al. Confidence intervals for measures of interaction. Epidemiology 1996;7: 286–90.
•
Brown and Harris (1978) Social origins of depression: A study of psychiatric disorder in women, New York: Free Press.
•
Campbell UB, Gatto NM, Schwartz S (2005) Distributional interaction: Interpretational problems when using incidence odds ratios to assess interaction.
Epidemiologic Perspectives and Innovations 2:1
•
Cheung YB (2007) A Modified least-squares regression approach to the estimation of risk difference
•
Chu, Nie, Cole. Estimating the relative excess risk due to interaction: a Bayesian approach. Epidemiology. 2011; 22:242-248
•
Darroch J: Biologic synergism and parallelism. American Journal of Epidemiology 1997, 145:661-668.
•
Hosmer DW, Lemeshow S. Confidence interval estimation of interaction. Epidemiology 1992;3:452–6.
•
Kalilani L and Atashili J (2006). Measuring additive interaction using odds ratios. Epidemiologic Perspectives and Innovations. 3:5.
•
Knol MJ van der Tweel I, Grobbee DE Numans M, Geerlings M (2007) Estimating interaction on an additive scale between continuous determinants in a logistic
regression model.
•
Kuss, Schmidt-Pokrzywniak and Stang (2010) Confidence Intervals for the Interaction Contrast Ratio Epidemiology, 21(2), 273-4.
•
Lundberg M, Fredlund P, Hallquist J, Diderichsen F. (1996). A SAS Program Calculating Three Measures of Interaction with Confidence Intervals. Epidemiology, Vol. 7,
No. 6 (Nov., 1996), pp. 655-656Published
•
Greenland S. Model-based estimation of relative risks and other epidemiologic measures in studies of common outcomes and in case-control studies. Am J
Epidemiol. 2004;160(4): 301–305.
•
McKee and Vilhjalmsson (1986) Life stress, vulnerability and depression: A methodological critique of Brown et al., Sociologyu 20: 589.
•
Nie L, Chu H, Li F, Cole S (2010) Relative Excess Risk Due to Interaction Resampling based confidence intervals, 21 (4), 552-556.
•
Richardson B and Kaufman JS (2009) Estimation of the Relative Excess Risk Due to Interaction and Associated Confidence Bounds, American Journal of Epidemiology,
169(6) 756-760.
•
Rothman KJ, Greenland S: Modern Epidemiology. Second edition. Philadelphia, Lippincott-Raven Publishers; 1998.
•
Schwarz S (2006) Modern Epidemiologic Approaches to Interaction: Applications to the study of Genetic Interactions.. In Genes, Behavior, and the Social
Environment: Moving Beyond the Nature/Nurture Debate http://www.nap.edu/openbook.php?record_id=11693&page=310
•
Spiegelman, D. and Hertzmark, E Easy SAS Calculations for Risk or Prevalence Ratios and Differences, American Journal of Epidemiology, 2005, 162, 199-205.
•
Skrondal A (2003) Interaction as Departure from Additivity in Case-Control Studies: A Cautionary Note, American Journal of Epidemiology, 158 (3) 251-258.
•
Storer BE Wacholder and Breslow Maximum likelihood fitting of general risk models to stratified data Appl Stat 1983: 32(2): 172-181
•
VANDERWEELE, T. J. & ROBINS, J. M. (2007). The identification of synergism in the sufficient-component cause framework. Epidemiol. 18, 329–39.
•
VANDERWEELE, T. J. & ROBINS, J.M. (2008). Empirical and counterfactual conditions for sufficient cause interactions. Biometrika 95, 49–61.
19

similar documents