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.
  • $\begingroup$ A possible way of doing this is by allowing more flexible forms of interaction terms. For example you might want the interaction be, say $\mathcal{C}^{2}$, and then add a penalized term in the risk function when you evaluate the regression model. This approach yields satisfying results and is still considered as a research area. $\endgroup$
    – Henry.L
    Jan 31, 2016 at 14:30
  • $\begingroup$ @Henry.L: would you mind helping me understand what you mean by $C^2$ and risk function in your comment? A reference or even a worked example would be very helpful! Thanks! $\endgroup$ Feb 1, 2016 at 6:22

3 Answers 3


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.

  • $\begingroup$ Thank you Dr Harrell for the reply! I found some description of the tensor spline in 2.7.2 of your 2015 book and 2.9.2 of your lecture note. Is there a worked example of the same in either of the sources above? Y = 1 is about 1.5% of the total sample (2 million observations), I suspect that frequency should be sufficiently large (is rarity of the event a problem?)? Lastly, is your last sentence meant to say that the percentile-cut is arbitrary because it is not individual characteristic? Thanks! $\endgroup$ Feb 1, 2016 at 5:09
  • $\begingroup$ See expanded answer $\endgroup$ Feb 1, 2016 at 12:31
  • $\begingroup$ Thanks! I read the narrative in the book (Chap 10) carefully. I am wondering why the decision was made to take the simpler model containing only $X1 \cdot X2$ as interaction term? Also, since the process is not blind to Y, will we run the risk of multiple comparison when we test out different non-linearity and interaction specification? $\endgroup$ Feb 1, 2016 at 22:44
  • $\begingroup$ That was a teaching example but you are right it involved some cheating. I wanted to show that you could get more evidence for some kind of interaction if you concentrated the interaction effect into fewer parameters. $\endgroup$ Feb 1, 2016 at 22:48

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.






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.

  • 2
    $\begingroup$ Dividing data by observed associations with $Y$ creates bias ("double dipping"). $\endgroup$ Jan 31, 2016 at 14:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.