# Nonlinear effect in an interaction term

If you have B, which is a 0/1 outcome variable, S, which is a continuous variable, and T, which is a treatment dummy variable, how can you show a hypothesized non-linear effect using regression results and a graph?

For example, I hypothesize that the treatment matters most for those in the middle of S's distribution.

The regression I have been running is B = S + T + ST.

Any cites on the topic would also be appreciated.

• What makes a source credible* in this instance? (* certainly no sources can be official) Jan 27, 2015 at 0:41

This is exactly why I switched from Stata to R and Frank's rms package (called Design back then) a few years ago.

Anyway, this somewhat hack-ish code will at least get you started. The syntax is a little outdated and there may be better ways to code this (haven't used Stata in a while), but here it goes

EDIT: Re-written after my morning coffee...

*** use automobile data
sysuse auto

*** create restricted cubic spline basis functions for mpg, with four knots
mkspline mpgsp = mpg, cubic nknots(4)

*** create the interactions
gen formpg1=foreign*mpgsp1
gen formpg2=foreign*mpgsp2
gen formpg3=foreign*mpgsp3

*** Regressing price on foreign and mpg allowing for non-linear interactions
xi: reg price i.foreign mpgsp* formpg*


To test the total interaction

test formpg1 formpg2 formpg3


Omit the first term for the test of any non-linear interaction terms, e.g.

test formpg2 formpg3


To get the global 4 d.f. test for T, which in this example is foreign, that Frank mentioned in his example above

test _Iforeign_1 formpg1 formpg2 formpg3


Just change reg to logit for logistic regression. To graph the result, you need to form the linear predictor, e.g. using predictnl, which I never managed to get right.

See a recent presentation by Patrick Royston at http://www.stata.com/meeting/germany12/abstracts/desug12_royston.pdf for some ideas.

Hope this helps.

• Great explanation on how to do the same in Stata. I would also throw a plug for Maarten Buis' mkspline2 (a user written Stata addon) which adds several post-regression conveniences. Jan 27, 2015 at 16:28
• I've never seen mkspline2 before, thanks for mentioning it. Jan 27, 2015 at 18:59

The following uses the R rms package using ordinary least squares modeling, and models the nonlinear effect smoothly using a restricted cubic spline with 4 knots at default knot locations. This generates one linear component and 2 nonlinear components for a total of 3 parameters per treatment group.

require(rms)
f <- ols(B ~ rcs(S, 4) * T, data=mydata)
anova(f)    # tests for interaction (shape differences across T, 3 d.f.)
# anova includes a test for nonlinear interaction
# also provides a global test for T, 4 d.f.
plot(Predict(f, S, T))   # shows 2 estimated curves for 2 values of T
ggplot(Predict(f, S, T))  # will be in next release; uses ggplot2


The plots include 0.95 pointwise confidence bands. There is an option to use simultaneous confidence bands instead.

Because I saw "ols" mentioned elsewhere I neglected to notice that the response variable is categorical. To fit the logistic regression model instead of an ols model, substitute lrm( ) for ols( ). No other code changes are needed. You can use summary(f, ...) to get odds ratios for T or S. By default the odds ratio for S will be the inter-quartile-range effect of S at the reference (most frequent) level of T.

• Whoops. Fixing now. Jan 27, 2015 at 4:42
• @FrankHarrell Can this be implemented in Stata?
– LF12
Jan 27, 2015 at 4:44
• Stata nicely handles restricted cubic splines but you may have to construct the interaction terms yourself, and the anova is not as comprehensive. Jan 27, 2015 at 4:46
• could you update the above in Stata, I can then accept
– LF12
Jan 27, 2015 at 6:50
• Writing this as an answer instead... Jan 27, 2015 at 8:21

Basically the model would be $$g(y) = X'\beta+\displaystyle\sum_j f_j(Z_j)+\epsilon$$

or in your specific case $$B = logit\left(f(S,T)\right)$$ In R, you could use the mgcv package, and run something like

library(mgcv)
m = gam(B~te(S,T),family=binomial)


which would give you a nonparametric interaction of the two variables. If you wanted to separate out main effects from the interaction effect, you could equivalently fit

m = gam(B~ti(S)+ti(T)+ti(S,T),family=binomial)


you can then look at contour plots of your estimated interaction via plot(m,pages=1, scheme=2) (I prefer the contour plots, myself), or you could use the vis.gam function to look at predicted values.

Or, if your treatment T is binary, you might fit

m = gam(B~s(S,by=as.factor(T)),family=binomial)


The textbook on all of this is made to go with the R package, and is here this, by Simon Wood.

Also you'll want to check ?te, ?ti, etc.

• I have, but really need to keep it to an OLS framework. Maybe a hazard model
– LF12
Jan 26, 2015 at 22:55
• Then why not just a polynomial expansion? Jan 26, 2015 at 23:11
• Multiple attempts at transformation so as to get good residual plots induces a distortion of inferential quantities such as $P$-values, standard errors, and confidence limits. Jan 27, 2015 at 17:05