How to best model interaction effect of two continuous predictor variables?

Question

Consider the following problem:

In a logistic regression model, we believe that two continuous predictor variables $X_1$ and $X_2$ impact the probability of event. It is hypothesized that the effect of $X_2$ on event depends on the value of $X_1$. More specifically, as level of $X_1$ increases, the effect of an unit increase in $X_2$ (at some fixed reference level of $X_2$) on log event odd first increases then plateaus out.

How may we model such interaction effect? I have considered the following three approaches, but none seem satisfactory. What would you suggest? Any comments on the three approaches would be welcomed in the comment section as well!

1. Include the interaction term $X_1\cdot X_2$ in the model. However, I have trouble getting meaningful interpretation of its parameter estimate.
2. Discretize $X_1$ into three indicator variables, $X_{1l}$, $X_{1m}$ and $X_{1h}$. $X_{1l}$ is the indicator for $X_1$ falling into its lowest 33-rd percentile. Include $X_{1l} \cdot spline{X_2}_1$, $X_{1m} \cdot spline{X_2}_2$ and $X_{1h} \cdot spline{X_2}_3$ where $spline{X_2}_i$ is natural cubic spline effect on $X_2$ for each $i$. -- The parameter for this model seems easier to interpret, however, the decision to split into 3-tiles seem arbitrary.
3. (A variation on 2) Split the modeling population into three pieces - the first one containing $X_1$ with lowest 33-tile values, second one the middle 33-tile, and the third one the highest 33-tile values. Fit a model on each of the segments using a cubic spline effect on $X_2$ together with the other predictor variables.

Answer 1

Any solution that involves choosing arbitrary cutpoints, other than for computing knot locations for spline functions, is to be avoided. The direct use of cutpoints as you suggest will result in discontinuities and lack of fit. It is far better to spend time learning to interpret interaction effects.

You didn't state the frequency of $Y=0$ and $Y=1$ which will be crucial in deciding what to do next. When sufficiently large frequencies you can avoid linearity (in the logit) assumptions and use things like restricted cubic splines for main effects and tensor spline products of them for interaction effects. This allows you to create a 3-D surface as well as a series of curves for showing the estimates, and you can form specific contrasts to estimate any effect you want (e.g., effect of change $X_{1}$ from $a$ to $b$ when $X_{2}=c$. [Follow-up: you have a very large absolute number of events and non-events so you can fit as complex a model as you desire. Imbalance in the number of events and non-events is not a significant problem.]

My course notes at http://biostat.mc.vanderbilt.edu/rms under Materials go into great detail. See the long example in the logistic regression chapter where interaction surfaces for age and cholesterol are derived.

It is important to note that percentiles are population characteristics and not characteristics of individuals. Percentiles are useful in picking spline knot locations because in restricted cubic splines the locations don't matter much and we want to typically put knots where data are dense. Using percentiles to directly code predictors ignores the fact that predictors affect individual outcomes in a physical way regardless of the distribution of covariates in the particular sample you've drawn.

Answer 2

Note: This might not serve as an answer, yet I want to point out some points that should be recognized for such a problem.

The first thing I recognized is that this problem seems like a classical problem of model selection. So the easiest way is to screen using something like AIC and BIC. A more advanced way of doing this is to use methods from model selections.

In modern texts, by modern I mean after E.P.Box, we usually do not separate some effect but talk about the regression function as a whole. A naive example is that we use some k-cutoff Fourier basis consisting of only k base functions instead of all $sin(nx),cos(nx)$ to model a regression function. We can give some interpretation of each term after fitting the regression model. Say for $k=4$, the term $a_0 + a_1\cdot sin(x) + b_1\cdot cos(x)$ is considered as "main effect" while the higer order term "$a_2\cdot sin(2x)+b_2\cdot cos(2x)$" as "interaction". However true, what we care is the regression function itself instead of the meaning of each term.

The very basics can be found in Rigollet's Notes[3] which talked a lot about some penalized methods and the bound imposed on them.

Now, what I am talking about is that we can penalize the form of regression function for its smoothness like [2]. A readable introduction can be found at [1].

What I meant in the comment to your question in the previous comment is [2], or we can even use the Bayesian model selection mentioned in [1]. If you have already determine that you must select from those 3 forms of interactions, that is find since you can run a model selection procedure directly and see how well each model fits; if you have not yet determine that, you would probably want to write a model in the form $Y=\alpha X_1 +\beta X_2 +f +\epsilon$ where the regression function is actually $\alpha X_1 +\beta X_2 +f$ in which $\alpha X_1,\beta X_2$ are the main-effect indicators. Then you can regard $f$ as some "interaction" in $\mathcal{C}^2$ or $L^{2}$ depending on how much freedom you want. The methods I referred to above will work neatly to give you the "best" $f\in\mathcal{C}^2$ or $L^{2}$ under MSE risk.

As you indicated, logistic model is also treated quite well in this framework. But I think such a selection with a risk other than MSE is still under research.

However, @Frank Harrell pointed out something worth noticing. If you decided to interpolate that, close attetion to the choice of cuttings is needed.

Reference

[1]http://www.stat.ufl.edu/archived/casella/Papers/BL-Final.pdf

[2]https://cran.r-project.org/web/packages/penalized/vignettes/penalized.pdf

[3]http://www-math.mit.edu/~rigollet/PDFs/RigNotes15.pdf

Answer 3

I would try approach 2 or 3, as those seem to be the most robust given your description of the interaction. It seems that including X1*X2 would not effectively the manner in which the interaction is occurring.

A fourth option would be to determine what the interaction effect of X1 and X2 is mathematically and use critical points along that trend (i.e. where changes in slope occur or max and min) to divide your data rather than arbitrarily splitting your model into 3 pieces.