Regression Strategies – Understanding Interaction with %ia% in rms and Three-Way Interactions

interactionregressionregression-strategiesrmssplines

In this very illustrative post on evaluating added value of predictors by Frank Harrell, he codes a logistic regression model as such:

lrm(sigdz ~ rcs(age,4) * sex + rcs(choleste,4) + rcs(age,4) %ia%
         rcs(choleste,4), data=acath)

The "%ia%" expression is new to me. He justifies its use:

The nonlinear interaction between age and cholesterol is a restricted
one such that terms that are nonlinear in both predictors are
excluded. This is to save degrees of freedom.

Although I do not quite understand what's behind the explanation, it does make sense to me that we can make a less-demanding interaction, since coding continuous-continuous interactions of spline variables does seem to make model heavier.

Q1: What does %ia% mean, in coding terms? Couldn't find much refs on the expression.

Q2: If I wanted to code a three-way interaction between two continuous predictors modelled with rcs() and a categorical one, how could I code in this df-sparing %ia% manner? (in the case of the post example, involving cholesterol, age and sex, for example)

Know coding Q's are better suited to SO, but this has a more stats background.

Thanks

Best Answer

This would be better posted in https://stackoverflow.com or even better at https://discourse.datamethods.org/t/rms-discussions but here goes:

Suppose you had two predictors $a, b$ that are modeled as quadratic effects using the rms regular polynomial function. The model would be specified as y ~ pol(a,2) * pol(b, 2). With full interactions you'd have these terms in the model: $a, b, a^{2}, b^{2}, ab, a^{2}b, ab^{2}, a^{2}b^{2}$. If you use the restricted interaction operator in the R rms package you'd drop $a^{2}b^{2}$, the doubly nonlinear term. %ia% does not extend to three-way interactions.

Related Solutions

Solved – Interaction of risk factor with age in survival model with age as time scale

Bendix Carstensen has a nice presentation in which he shows how this can be accomplished by splitting individual follow up time into small pieces and use a Poisson model for the rates thereby using the time scale (in your case age) as a covariate rather than part of the outcome. If you use R, the Epi package contains the tools you need.

Another option is to extend the Cox model to include time-dependent coefficients. See this vignette for the survival package in R. In that vignette it is also explained why it is not appropriate to include a covariate * ~~age~~ follow-up time interaction directly in the Cox model.

EDIT: Removed the bit about baseline age. That would of course not make sense with age as the time scale. Typed a bit too fast...

EDIT 2: Yes, the age * covariate interaction is an example of a time dependent coefficient and one of the ways for assessing the proportional hazards assumption.

As for the paper, it is possible that the critisism could apply. I can't find any info that the authors allowed for time-dependent coefficients. Including the age * covariate interaction is certainly possible, and R prints a warning message but estimates the coefficients. See the example below which is very inspired by the vignette mentioned above.

require(survival)
data(veteran)

## Use age as the time scale
## age * karno in the model
fit1 <- coxph(Surv(age, age + time, status == 1) ~ trt + celltype + karno + age:karno, data = veteran)
Warning message:
In coxph(Surv(age, age + time, status == 1) ~ trt + celltype + age *  :
  a variable appears on both the left and right sides of the formula

## karno as a time-dependent coefficient, simple interaction similar to fit1
fit2 <- coxph(Surv(age, age + time, status == 1) ~ trt + celltype + karno + tt(karno),
          tt = function(x, t, ...) x * t, data = veteran)

## Coefficients from fit1
> fit1
Call:
coxph(formula = Surv(age, age + time, status == 1) ~ trt + celltype + 
karno + age:karno, data = veteran)

                       coef exp(coef)  se(coef)     z       p
trt                2.59e-01  1.30e+00  2.05e-01  1.26  0.2072
celltypesmallcell  8.32e-01  2.30e+00  2.73e-01  3.05  0.0023
celltypeadeno      1.17e+00  3.21e+00  2.95e-01  3.95 7.8e-05
celltypelarge      3.95e-01  1.48e+00  2.83e-01  1.39  0.1633
karno             -3.33e-02  9.67e-01  1.07e-02 -3.10  0.0020
karno:age          4.75e-05  1.00e+00  1.77e-04  0.27  0.7890

Likelihood ratio test=62  on 6 df, p=1.75e-11
n= 137, number of events= 128

## Coefficients from fit2

> fit2
Call:
coxph(formula = Surv(age, age + time, status == 1) ~ trt + celltype + 
karno + tt(karno), data = veteran, tt = function(x, t, ...) x * 
t)

                       coef exp(coef)  se(coef)     z       p
trt                1.70e-01  1.19e+00  2.04e-01  0.84 0.40333
celltypesmallcell  9.08e-01  2.48e+00  2.75e-01  3.30 0.00097
celltypeadeno      1.22e+00  3.40e+00  3.00e-01  4.08 4.5e-05
celltypelarge      4.04e-01  1.50e+00  2.83e-01  1.42 0.15422
karno             -4.69e-02  9.54e-01  8.33e-03 -5.63 1.8e-08
tt(karno)          1.15e-04  1.00e+00  4.71e-05  2.45 0.01443

Likelihood ratio test=68.2  on 6 df, p=9.45e-13
n= 137, number of events= 128

Note the difference in the estimated coefficients. Reading the vignette again I see that section 4.2 deals with follow-up time * covariate interaction and not necessarily age * covariate interaction. But I think the example illustrates that results can be misleading when not properly allowing for time-dependent coefficients.

Hope this helps!

Cox PH Model – Interpreting Coefficient Changes Involving Interactions

That is exceptionally well laid out. As I detail in RMS it is not so fruitful to emphasize individual coefficient interpretation except for the highest-order term. That is, the interaction term can be interpreted (especially when one of the interacting variables is binary) but lower order terms (often called "main effects" but this is a bit misleading) are heavily depending on centering and indicator variable coding. Think instead of contrasts, which are most easily calculated in terms of differences in predicted values (here, differences in log relative hazard, so that the anti-log of the contrast is a hazard ratio).

What I'm about to describe doesn't work as well if any of the effects is nonlinear. In that case the choice of intervals is arbitrary and can even be misleading if there is a flat or non-monotonic piece of one of the covariate effects. Supposing effects are linear (which has to be the case for binary predictors), suppose you want to estimate the effect of x1 where x1 interacts with x2. Pick a series of values of x2. For each x2 contrast the relative log hazard when x1 is set to a vs. when x1 is set to b. For binary x1, a=0, b=1. For continuous x1 you might let a=first quartile of x1, b=upper quartile of x1.

In the R rms package the contrast function (long name contrast.rms) makes contrasts over a series of x2 values very easy to do no matter how complex the model. When getting differences in predicted values one can easily compute the standard error of each difference, as contrast does.

As shown in the RMS course notes, for the continuous-continuous x1 x2 case it works well to graph a color image plot where the 'temperature' of the the color represents the log hazard or survival probability.

If you evaluate any of the contrasts manually you'll see that the x1 effect is a simple function of x2. And as you referenced, large negative interaction terms can lead to reversals of x1 effects as x2 changes.

Note that the univariable estimates you provided are not helpful.

Best Answer

Related Solutions

Solved – Interaction of risk factor with age in survival model with age as time scale

Cox PH Model – Interpreting Coefficient Changes Involving Interactions

Related Question