*! version 0.1.2 2017-04-20

Title

Measures of interaction contrast (biological interaction) - ic, ici and icp

Summary

It has become more common to investigate not only single factors effect
on a outcome but also to look at the interaction between factors, as is
facilitated by the last decades massive increase in computer resources to
gather and analyses large databases. Several approaches haves been promoted.

For the analyzes of binary (or count) data, a fairly well established
position holds, that loglinear models estimating measures of relative risk type
represent a convenient choice, whereas interactions are often best interpreted
if estimated on a linear/additive scale.

This code implement the procedure described in
Hosmer & Lemeshow (1992) and
Alfredsson et all (2005) purely in Stata.

These two sources however both uses odds ratios as approximations for
relational risks. This code on the other hand opens up for using proper
estimates of relational risks.

Authors and support

Author: Niels Henrik Bruun,
Section for General Practice,
Dept. Of Public Health,
Aarhus University

Coauthor: Morten Fenger-Grøn,
Research Unit for General Practice,
Aarhus University

Coauthor: Anders Prior,
Research Unit for General Practice,
Aarhus University

Support: nhbr@ph.au.dk

Syntax

ic varlist (min=3)
[, ic_options]

ici matrix_estimates matrix_variance
[macro_with_3_labels]

icp [icp_options] :
regression_model_to_estimate_relative_risks

ic_options

referenceA The reference value for the first variable in the interaction.
Default is 0.

referenceB The reference value for the second variable in the interaction.
Default is 0.

exposedA The exposed value for the first variable in the interaction.
Default is 1.

exposedB The exposed value for the second variable in the interaction.
Default is 1.

show|SHOW Show the underlying regression estimating the relative risks.
Default is no show.
The option show doesn't remove the dummy variables used in the regression. This way the regression is reproducible.

rrby Indicate which estimation method to use for estimating
adjusted RRs, see Cummings (2009).

The following methods are implemented:

* or the classical method using logistic regression to approximate RRs with ORs as described in Alfredsson et all (2005) (default)

* poisson estimating RRs by robust poisson regression (see Method 4 in Cummings (2009)

* binomial using binomial regression with initiating estimates from a poisson regression for the RR estimates (see Method 2 in Cummings (2009)

icp_options

referenceA The reference value for the first variable in the interaction.
Default is 0.

referenceB The reference value for the second variable in the interaction.
Default is 0.

exposedA The exposed value for the first variable in the interaction. Default is 1.

exposedB The exposed value for the second variable in the interaction.
Default is 1.

show|SHOW Show the underlying regression estimating the relative risks.
Default is no show.
The option show doesn't remove the dummy variables used in the regression.

Description

Tool for calculating 3 different measures of interaction contrast
(biological interaction) on an additive scale: relative excess risk due to
interaction [RERI], attributable proportion [AP] and synergy index [S].

Corresponding 95% confidence intervals aand two-tailed tests for no
interaction (RERI=0, AP=0, S=1, respectively) are given as well.

Assume that we want to analyze the effects of 2 dichotome factors, A and
B, on a dichotome outcome Y.
Define Pij as the average risk (the incidence proportion) when A = i and B = j,
i = 0,1 and j = 0,1.

And define the relative risk as RRij = Pij / P00.

The biological interaction measures is defined as

* RERI = RR11 - RR10 - RR01 - 1

* AP = RERI / R11

* S = (RR11 - 1) / (RR10 + RR01 - 2)

VanderWeele (2014) Rothman (2008).

The estimation/approximation of the RRs are described in
Hosmer & Lemeshow (1992) Alfredsson et all (2005),
Cummings (2009) and VanderWeele (2014).

Understanding biological interaction

Below are some (sometimes slightly modified) quotes to help understand
biological interaction.

The term statistical interaction is intended to denote the
interdependence between the effects of two or more factors within the confines
of a given model of risk. The evaluation of interaction depends on the model
chosen. (Rothman (1980))

Biologic interaction may be defined as the interdependent operation of
two or more causes to produce disease. (Rothman (1980))

Although we may not be able to estimate the additive interaction, i.e.
p11 - p10 - p01 + p00, directly, we can still proceed as follows.

If we divide by p00 we obtain the following: RERI = RR11 - RR10 - RR01 + 1"
VanderWeele (2014)

"However, again, only the direction, rather than the magnitude, of RERI is
needed to draw conclusions about the public health relevance of interaction.
VanderWeele (2014)

Another measure of additive interaction that is sometimes used is called
the attributable proportion and is defined as: AP = RERI / RR11 and essentially
measures the proportion of the risk in the doubly exposed group that is due to
the interaction itself. VanderWeele (2014)

"The logistic regression model and the Cox regression model are probably
the most commonly used statistical models in epidemiologic analysis to day.

Because these models are exponential they are inherently multiplicative and
become additive only after logarithmic transformation.

Thus, absence of an interaction term in such a model implies a multiplicative
relation between the disease rates and the presence of an interaction term
implies departure from multiplicativity, rather than from additivity.

Therefore, the interaction term, in one of these models, has no direct
relevance for the issue of whether or not biological interaction is present.

However, the presence of biological interaction can still be asessed from the
results of a logistic regression model or a Cox regression model, but this
requires that the model is defined in a special way and that the analysis is
done adequately." Alfredsson et all (2005)

"In this study we showed that calculating measures of interaction on an
additive scale using preventive factors can give inconsistent results.

Researchers should therefore be aware to not use preventive factors to
calculate these measures unless they have been recoded."
Knol et all (2011)

ic

Given a variable list staring with outcome, factor A, factor B end
possibly some adjustment variables for the regression the ic estimates the RR
for A, B, the interaction between A and B, RERI, AP and S as well as their 95%
confidence intervals.

Note that the interaction term is automatically created
based on variable 2 and variable 3 (outcome is variable number 1) in the
regression.

ici

The immediate form of ic. Requires a 1 by 3 or 3 by 1 matrix
containing the logs of the RRs and their 3 by 3 covariance matrix as input.
A local containing labels like "labelA labelB labelAB" is optional.

icp

This is the advanced prefix form letting the user define regression
settings to use to get estimates of log relative risks to use interaction
contrast.

Note that the interaction term is automatically created
based on variable 2 and variable 3 (outcome is variable number 1) in the
regression.

Also note that icp only accepts logistic, logit, binreg, nbreg,
poisson or stcox (stcox needs a stset first) as regression inputs.

The key point for the regression is that it must produce estimates or at least
approximations to relative risks.

Examples

Code for generating example data
The essential output from icp, ici and ic
Prefix icp and how it works
ici and how it works
ic and how it works
Using ic on the rheumatoid arthritis dataset
Using icp with stcox

Code for generating example data

Below is code to generate datasets to use in the examples.

To reproduce the datasets just copy the whole code blocks into the
command window and press enter.

The oral cancer dataset

To get the example data run the following code:

. clear all
. capture cls // Cls only work in Stata 13
. input alcohol smoking cancer count
0 0 0 20
0 0 1 3
0 1 0 18
0 1 1 8
1 0 0 12
1 0 1 6
1 1 0 166
1 1 1 225
end
* Expand the dataset (You can't use weights in ic):
. expand count
. drop count
* And add labels and label values:
. label variable alcohol "Alcohol use"
. label variable smoking "Smoker"
. label variable cancer "Oral cancer"
. label define no_yes 0 "No" 1 "Yes"
. label define cc 0 "Controls" 1 "Cases"
. label values alcohol no_yes
. label values smoking no_yes
. label values cancer cc
* And hence the dataset is (Table 1 transposed, more or less):
. table smoking cancer alcohol

This is the example dataset from Hosmer & Lemeshow (1992).
The data concerns smoking and alcohol use in relation to oral cancer among male
veterans under age 60.

The rheumatoid arthritis dataset

To get the example data run the following code:

. clear all
. capture cls // Cls only work in Stata 13
. input se sex smoke ra count
0 0 0 0 73
0 0 0 1 63
0 0 1 0 43
0 0 1 1 44
0 1 0 0 13
0 1 0 1 16
0 1 1 0 17
0 1 1 1 16
1 0 0 0 75
1 0 0 1 98
1 0 1 0 46
1 0 1 1 79
1 1 0 0 29
1 1 0 1 29
1 1 1 0 11
1 1 1 1 28
2 0 0 0 25
2 0 0 1 40
2 0 1 0 11
2 0 1 1 47
2 1 0 0 6
2 1 0 1 9
2 1 1 0 3
2 1 1 1 16
end
* Expand the dataset (You can't use weights in ic):
. expand count
. drop count
* Add labels and label values:
. label variable se "SE genes"
. label variable sex "Sex"
. label variable smoke "Smoking"
. label variable ra "Rheumatoid arthritis"
. label define se 0 "No SE" 1 "Single SE" 2 "Double SE"
. label values se se
. label define sex 0 "Women" 1 "Men"
. label values sex sex
. label define smoke 0 "Never smokers" 1 "Current Smokers"
. label values smoke smoke
. label define ra 0 "No RA" 1 "RA"
. label values ra ra
* And hence the table 4 can be reproduced by
. egen col = group(se ra), label
. egen row = group(sex smoke), label
. table row col

Table 4 in Padyukov et all (2004). Relative risk of
developing rheumatoid arthritis among subjects exposed to different combinations
of current cigarette smoking habits and SE genes and attributable proportion
due to interaction between cigarette smoking and SE genes (single, double,
or any), by sex.

Back to examples

The essential output from icp, ici and ic

The output (example below from ic with no labels set, compare with
table 5 in Hosmer & Lemeshow (1992)) starts with the 3
relative risks used in the calculation of RERI, AP and S. Then comes the
estimates of RERI, AP and S.

Column 1 are the labels for the estimates, column 2 contains the estimates,
column 3 contains the P-values for testing away the measure at hand versus the dual
alternative (for S it is the P-value for testing S = 1).

The 2 final columns (4 and 5) are lower and upper bound in the 95%
confidence interval.

To reproduce the result below run the code from example 1
to generate the oral cancer dataset.

Then run the following command:

. ic cancer alcohol smoking, rrby(or)

And the output is:

+------------------------------------------------------------------------------+
| Summary measures | Estimates P-value Lower bound Upper bound |
|-----------------------------+------------------------------------------------|
| alcohol_NOT_smoking | 3.3333 0.1303 0.7006 15.8593 |
| smoking_NOT_alcohol | 2.9630 0.1480 0.6800 12.9097 |
| alcohol_AND_smoking | 9.0361 0.0005 2.6413 30.9131 |
|-----------------------------+------------------------------------------------|
| RERI | 3.7398 0.1887 -1.8367 9.3164 |
| AP | 0.4139 0.0957 -0.0731 0.9008 |
| S | 1.8705 0.2483 0.6460 5.4156 |
+------------------------------------------------------------------------------+
Interaction exists if RERI != 0 or AP != 0 or S != 1

Back to examples

Prefix icp and how it works

Using the oral cancer dataset what is done is shown below. In the example
the log odds are estimated using a logistic regression with option coef.

* Generate dummy variables for the regression:
. generate marker = inlist(alcohol, 0, 1) & inlist(smoking, 0, 1)
. generate ind11 = (alcohol == 1 & smoking == 1) if marker
. replace ind11 = 0 if !missing(alcohol, smoking) & !marker
. generate ind10 = (alcohol == 1 & smoking == 0) if marker
. replace ind10 = 0 if !missing(alcohol, smoking) & !marker
. generate ind01 = (alcohol == 0 & smoking == 1) if marker
. replace ind01 = 0 if !missing(alcohol, smoking) & !marker
* The values 0 and 1 above can be changed by options in icp
* Run the regression on the dummy variables:
. logistic cancer ind10 ind01 ind11, coef
* Get the log relative risk estimates from the regression:
. matrix b = r(table)
. matrix b = b[1, 1..3]
* Get the variance from the regression:
. quietly vce
. matrix V = r(V)
. matrix V = V[1..3, 1..3]
* The data used ici is:
. matlist (b', V)
* Call ici with estimates:
. ici b V
* And clean up:
. drop ind10 ind01 ind11 marker
. matrix drop b V

This is much like the description in Hosmer & Lemeshow (1992)
and Alfredsson et all (2005).
With the exception that their spreadsheet is replaced with the Stata command
ici which described further down. Instead of these commands the prefix
icp could have been used:

. icp: logistic cancer alcohol smoking, coef

So in other words icp takes a regression command with at least 3
variables and lets the first be the outcome, the next 2 variables be the factors
used as base for the dummy variables and the rest of the variables are used for
adjusting.

And if the underlying regression is wanted then use the option show:

. icp, show: logistic cancer alcohol smoking, coef

Note that the option show doesn't remove the dummy variables used in the
regression.

If the calculation of the dummy variables should be based on other values
than 0 and/or 1 then use the options referenceA, referenceB, exposedA and/or
exposedB.

So

. icp, referenceA(2) referenceB(3) exposedA(4) exposedB(5): logistic cancer alcohol smoking, coef

would be the same as the following dummy generating commands:

. generate marker = inlist(alcohol, 2, 4) & inlist(smoking, 3, 5)
. generate ind11 = (alcohol == 4 & smoking == 5) if marker
. replace ind11 = 0 if !missing(alcohol, smoking) & !marker
. generate ind10 = (alcohol == 4 & smoking == 3) if marker
. replace ind10 = 0 if !missing(alcohol, smoking) & !marker
. generate ind01 = (alcohol == 2 & smoking == 5) if marker
. replace ind01 = 0 if !missing(alcohol, smoking) & !marker

Please note that this might lead to a nonsense regression.

Of course it is possible to chose your own regression after icp, eg:

. icp: poisson cancer alcohol smoking, irr vce(robust)

or

. icp: binreg cancer alcohol smoking, rr

However icp only accepts logistic, logit, binreg, nbreg, poisson or
stcox (stcox needs a stset first) as regression inputs. The key point for the
regression is that it must produce estimates or at least approximations to
relative risks.

Back to examples

ici and how it works

ici can be used in interactively much the same way as the spreadsheet
in Alfredsson et all (2005).

It takes 3 estimates of the log relational risk (or approximating odds ratio) as
first argument in a row or column vector (matrix), second argument is the
3 by 3 variance matrix and finally (and optionally) labels to use in the output.

So if 2 matrices and a macro for labels are defined as:

. matrix b = 1.2039728, 1.0861898, 2.2012326
. matrix V = .63333333, .38333333, .38333333
. matrix V = V \ .38333333, .56388889, .38333333
. matrix V = V \ .38333333, .38333333, .39380187

. local labels "alcohol_NOT_smoking smoking_NOT_alcohol alcohol_AND_smoking"

Note that matrix b can b either a row or a column vector.

Also one can enter a matrix as a string argument right away.

So the call of ici could be:

. ici b V

or

. ici "1.2039728, 1.0861898, 2.2012326" V "`labels'"

So the parameter label is optional.

Back to examples

ic and how it works

ic is probably how most would like to use the code. Based on
Cummings (2009) it offers some standard ways of estimating
relative risks using ic. This example is based on the oral cancer
dataset.

The command below perform the calculations based on odds ratio estimates
, option rrby(or), and show, option show, the underlying regression output:

. ic cancer alcohol smoking, show rrby(or)

This is the classical method (and the default).

Here the poisson regression is used as a base for the relative risk
estimates, option rrby(poisson). Also the controls and the cases are reversed
for smoking, options referenceB(1) exposedB(0):

. ic cancer alcohol smoking, referenceB(1) exposedB(0) rrby(poisson)

According to eg Cummings (2009) the poisson regression
makes robust estimates of the relative risks. However it tends to give too wide
confidence intervals.

Finally there is a stable estimation, option rrby(binomial), of the
relative risks using estimates from a similar poisson regression as initiation
for the binomial regression.

. ic cancer alcohol smoking, rrby(binomial)

The binomial regression should give better estimates of the confidence
intervals. But the algorithm isn't as stable. Hence we initiate with point
estimates from a similar poisson regression.

Note that if alcohol had eg 3 levels (eg 0 "No" 1 "Moderate" 2 "Extreme").
Then the interaction estimates for "No" (reference) versus "Moderate" (exposed)
would be the traditional usage of ic:

. ic cancer alcohol smoking

And then the interaction estimates for "No" (reference, default = 0)
versus "Extreme" (exposed) would be:

. ic cancer alcohol smoking, exposedA(2)

Back to examples

Using ic on the rheumatoid arthritis dataset

In this example we look at a typical use of ic. First we need to run
the code to get the rheumatoid arthritis dataset.

Then we first need to generate a dichotome variable se2 of se to reproduce the
result from Padyukov et all (2004):

. recode se (0=0 "No SE") (1 2 = 1 "SE"), generate(se2)

Note that we also want to adjust for sex in the calculations.So by
running the following command one gets:

. ic ra smoke se2 i.sex, rrby(or)
+------------------------------------------------------------------------------+
| Summary measures | Estimates P-value Lower bound Upper bound |
|-----------------------------+------------------------------------------------|
| smoke_NOT_se2 | 1.0843 0.7366 0.6766 1.7377 |
| se2_NOT_smoke | 1.4159 0.0724 0.9689 2.0693 |
| smoke_AND_se2 | 2.5999 0.0000 1.7206 3.9287 |
|-----------------------------+------------------------------------------------|
| RERI | 1.0997 0.0141 0.2217 1.9777 |
| AP | 0.4230 0.0028 0.1461 0.6999 |
| S | 3.1985 0.1310 0.7072 14.4658 |
+------------------------------------------------------------------------------+
Interaction exists if RERI != 0 or AP != 0 or S != 1

The factor smoke_AND_se2 is significantly different from 1
(P-value = 0.0000), so there is a (statistical) multplicative interaction.
And also by looking at RERI and AP there is a significant biological interaction
(P-values are 0.0141 0.0028 respectively).

Now one could be interested in whether this result were the same if there
were a single SE gene present (exposed == 1):

. ic ra smoke se i.sex, exposedB(1) rrby(or)
+------------------------------------------------------------------------------+
| Summary measures | Estimates P-value Lower bound Upper bound |
|-----------------------------+------------------------------------------------|
| smoke_NOT_se | 0.6815 0.0751 0.4468 1.0396 |
| se_NOT_smoke | 0.8339 0.2979 0.5923 1.1739 |
| smoke_AND_se | 1.2833 0.2112 0.8680 1.8975 |
|-----------------------------+------------------------------------------------|
| RERI | 0.7679 0.0044 0.2394 1.2965 |
| AP | 0.5984 0.0004 0.2643 0.9325 |
| S | . . . . |
+------------------------------------------------------------------------------+
Interaction exists if RERI != 0 or AP != 0 or S != 1

Here we note that values for S could not be calculated. This is due to estimates
being to big to exponentiate. However here is only a biological interaction
present (P-values: smoke_AND_se=0.2112, RERI=0.0044, AP=0.0004).

Finally one could be interested in whether this result was the same if
there was a double SE gene present (exposed == 2):

. ic ra smoke se i.sex, exposedB(1) rrby(or)
+------------------------------------------------------------------------------+
| Summary measures | Estimates P-value Lower bound Upper bound |
|-----------------------------+------------------------------------------------|
| smoke_NOT_se | 0.7865 0.2336 0.5297 1.1677 |
| se_NOT_smoke | 1.2505 0.3613 0.7738 2.0208 |
| smoke_AND_se | 3.5467 0.0000 1.9412 6.4801 |
|-----------------------------+------------------------------------------------|
| RERI | 2.5097 0.0230 0.3467 4.6727 |
| AP | 0.7076 0.0000 0.4579 0.9573 |
| S | 68.8603 0.6664 0.0000 1.5604e+10 |
+------------------------------------------------------------------------------+
Interaction exists if RERI != 0 or AP != 0 or S != 1

Once again the estimates leading to S have become to big. So they are to
be ignored. Here there are both a multplicative and a biological interaction
present (P-values: smoke_AND_se=0.0000, RERI=0.0230, AP=0.0000)

Back to examples

Using icp with stcox

Cox regression is treated a bit different with icp. First the example
data

. use http://www.stata-press.com/data/r12/drugtr.dta, clear
. egen age_grp = cut(age), group(2) label

Now icp can be used:

. icp, show: stcox drug age_grp, nolog cformat(%7.3f)

failure _d: died
analysis time _t: studytime

Cox regression -- Breslow method for ties

No. of subjects = 48 Number of obs = 48
No. of failures = 31
Time at risk = 744
LR chi2(3) = 27.98
Log likelihood = -85.923498 Prob > chi2 = 0.0000

-----------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
------------------+----------------------------------------------------------------
_drug_NOT_age_grp | 0.107 0.069 -3.46 0.001 0.030 0.379
_age_grp_NOT_drug | 1.967 0.927 1.44 0.151 0.782 4.952
_drug_AND_age_grp | 0.254 0.137 -2.55 0.011 0.089 0.729
-----------------------------------------------------------------------------------

+------------------------------------------------------------------------------+
| Summary measures | Estimates P-value Lower bound Upper bound |
|-----------------------------+------------------------------------------------|
| drug_NOT_age_grp | 0.1072 0.0005 0.0303 0.3795 |
| age_grp_NOT_drug | 1.9675 0.1507 0.7818 4.9516 |
| drug_AND_age_grp | 0.2540 0.0108 0.0885 0.7286 |
|-----------------------------+------------------------------------------------|
| RERI | -0.8207 0.3619 -2.5850 0.9437 |
| AP | -3.2315 0.3434 -9.9159 3.4529 |
| S | . . . . |
+------------------------------------------------------------------------------+
Interaction exists if RERI != 0 or AP != 0 or S != 1

Back to examples

Stored results

icp store the following in r():

Macros

r(labels) The labels used in the output from icp

Matrices

r(est) A column vector (3 x 1) containing the log estimates

r(var) A (3 x 3) matrix containing the variances

ici and ic store the following in r():

Macros

r(labels) The labels used in the output from icp

Matrices

r(est) A column vector (3 x 1) containing the log estimates of the relative risks: RR10, RR01 and RR11.

r(var) A (3 x 3) matrix containing the variances of the relative risks: RR10, RR01 and RR11

r(output) A (5 x 5) matrix containing the data for the result table

References

Andersson T, Alfredsson L, Kallberg H, Zdravkovic S, Ahlbom A.
Calculating measures of biological interaction
Eur J Epidemiol. 2005;20(7):575-9.

Cummings, P., Methods for estimating adjusted risk ratios, Stata Journal",
Stata Press 9-2 2009, 175-196(22),
Link to pdf

Hosmer, David W., Lemeshow, Stanley,
Confidence Interval Estimation of Interaction,
Epidemiology, Vol. 3, No. 5 (Sep., 1992), pp. 452-456

Knol, Mirjam J., VanderWeele, Tyler J., Groenwold, Rolf H. H., Klungel, Olaf H.,
Rovers, Maroeska M., Grobbee, Diederick E.
Estimating measures of interaction on an additive scale for preventive exposures
Eur J Epidemiol (2011) 26:433438 DOI 10.1007/s10654-011-9554-9

Padyukov L, Silva C, Stolt P, Alfredsson L, Klareskog L.
A gene-environment interaction between smoking and shared epitope genes in
HLA-DR provides a high risk of seropositive rheumatoid arthritis.
Arthritis Rheum 2004; 50: 3085-3092.

Rothman, K. J., Greenland, S., and Walker, A. M. (1980),
Concepts of interaction;
American Journal of Epidemiology, 112:467470

Rothman Kenneth J., Lash Timothy L., Greenland Sander,
Modern Epidemiology 3rd ed., 2008 Lippincott Williams & Wilkins, p. 74-83

VanderWeele T.J., Knol, M.J., A tutorial on interaction.
Epidemiologic Methods, Published Ahead of Print May 27, 2014,
DOI 10.1515/em-2013-0005,
Link to pdf