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