Implementation of Additive Interactions in Cox proportional hazards model?

[Scott Pilla]

I am performing analysis with stset data in long format using the cox proportional hazards model (stcox). I would like to test for additive (rather than multiplicative) interactions between 2 binary covariates. The paper below publishes and implementation of this in SAS. Does anyone have a method for implementing this in stata? Thank you for any help.

http://www.sciencedirect.com/science...4727970600264X

[Scott Pilla]

The stata team figured this out -- their message to me is copied below:

Dear Scott,

Thank you for raising this question and providing the reference.

I have read the paper you sent us. The method proposed by this paper has
not been implemented by the -stcox- command in Stata at present. I have
not seen any user-written commands for this either. However, the method
to test the additive effect of interaction in this paper is based on
the Cox PH model output. The author first fitted a cox PH model, then
wrote some commands to calculate the estimation and confidence interval
of "RERI" based on formula (4). In fact, we could easily calculate
"RERI" and perform test using the commands -nlcom- or -testnl- in Stata.

We could first fit the Cox PH model by -stcox-, then obtain an
estimation and test results for HERI by -nlcom- or -testnl-. For
illustration purpose, I have written a piece of code with my synthetic
data, which includes two covariates: gene and environment. My code and
its output are as follows.

Below I use the factor variable notation to include the interaction.
Type
help factor variable
to read more about the Stata factor variable notation.

. list in 1/5

+----------------------------+
| failur~e gene enviro~t |
|----------------------------|
1. | 2 1 1 |
2. | 3 0 1 |
3. | 4 1 1 |
4. | 5 1 1 |
5. | 6 0 1 |
+----------------------------+

. stset failuretime

failure event: (assumed to fail at time=failuretime)
obs. time interval: (0, failuretime]
exit on or before: failure

------------------------------------------------------------------------------
18 total observations
0 exclusions
------------------------------------------------------------------------------
18 observations remaining, representing
18 failures in single-record/single-failure data
97 total analysis time at risk and under observation
at risk from t = 0
earliest observed entry t = 0
last observed exit t = 11

. stcox i.gene##i.environment, nohr//fit the Cox model

failure _d: 1 (meaning all fail)
analysis time _t: failuretime

Iteration 0: log likelihood = -37.565721
Iteration 1: log likelihood = -30.75554
Iteration 2: log likelihood = -30.375598
Iteration 3: log likelihood = -30.324598
Iteration 4: log likelihood = -30.324419
Iteration 5: log likelihood = -30.324419
Refining estimates:
Iteration 0: log likelihood = -30.324419

Cox regression -- Breslow method for ties

No. of subjects = 18 Number of obs = 18
No. of failures = 18
Time at risk = 97
LR chi2(3) = 14.48
Log likelihood = -30.324419 Prob > chi2 = 0.0023

----------------------------------------------------------------------------------
_t | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-----------------+----------------------------------------------------------------
1.gene | 2.572562 1.310805 1.96 0.050 .0034316 5.141692
1.environment | .3904787 .7351701 0.53 0.595 -1.050428 1.831386
|
gene#environment |
1 1 | .4060818 1.100093 0.37 0.712 -1.75006 2.562224
----------------------------------------------------------------------------------


. stcox, coeflegend //display legend for all the coefficients to test

Cox regression -- Breslow method for ties

No. of subjects = 18 Number of obs = 18
No. of failures = 18
Time at risk = 97
LR chi2(3) = 14.48
Log likelihood = -30.324419 Prob > chi2 = 0.0023

----------------------------------------------------------------------------------
_t | Coef. Legend
-----------------+----------------------------------------------------------------
1.gene | 2.572562 _b[1.gene]
1.environment | .3904787 _b[1.environment]
|
gene#environment |
1 1 | .4060818 _b[1.gene#1.environment]
----------------------------------------------------------------------------------



. nlcom (RERI:exp(_b[1.gene]+_b[1.environment]+_b[1.gene#1.environment])-exp(_b[1.gene])-exp(_b[1.environment])+1)

RERI: exp(_b[1.gene]+_b[1.environment]+_b[1.gene#1.environment])-exp(_b[1.gene])-exp(_b[1.environment])+1

------------------------------------------------------------------------------
_t | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
RERI | 15.47599 23.00612 0.67 0.501 -29.61518 60.56715
------------------------------------------------------------------------------

. testnl exp(_b[1.gene]+_b[1.environment]+_b[1.gene#1.environment])-exp(_b[1.gene])-exp(_b[1.environment])+1=0

(1) exp(_b[1.gene]+_b[1.environment]+_b[1.gene#1.environment])-exp(_b[1.gene])-exp(_b[1.environment])+1 = 0

chi2(1) = 0.45
Prob > chi2 = 0.5011


Please note that -nlcom- and -testnl- both adopt the delta method, so
that we don't need to compute the variance manually as the paper shows.
From the output of both -nlcom- and -testnl-, we can see the test
result. From the output of -nlcom-, we can obtain the estimation and
confidence interval.

Note that you don't need to write a program to perform this test in
Stata. You just need to write one line to implement the formula (4)
from the paper.

Please let us know if you have any further questions.

[Carlo Lazzaro]

Thanks Scott for sharing this with the list.