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.

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

4 Answers 4


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.

  • 2
    $\begingroup$ 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. $\endgroup$ Jan 27, 2015 at 16:28
  • $\begingroup$ I've never seen mkspline2 before, thanks for mentioning it. $\endgroup$
    – ErikL
    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.

dd <- datadist(mydata); options(datadist='dd')  # facilitates plotting
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.

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

Have you considered using a generalized additive model? Wikipedia link here

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

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.

  • $\begingroup$ I have, but really need to keep it to an OLS framework. Maybe a hazard model $\endgroup$
    – LF12
    Jan 26, 2015 at 22:55
  • $\begingroup$ Then why not just a polynomial expansion? $\endgroup$ Jan 26, 2015 at 23:11
  • $\begingroup$ say more please $\endgroup$
    – LF12
    Jan 26, 2015 at 23:13
  • $\begingroup$ y = b0 + b1T + b2S + b3S*T+ b4S^2 + b5S^2*T ..... etc $\endgroup$ Jan 26, 2015 at 23:15
  • $\begingroup$ Depending on why you need to "keep it in an OLS framework", you might also fit a gam to binary data as if your response were actually gaussian. that is equivalent to a penalized ols regression $\endgroup$ Jan 26, 2015 at 23:16

You may argue a non-linear effect by showing that a non-linear model fits better. For example you could implement a piecewise linear model to take into account changes in the influence of S. Dependent on your hypothesis, you could also linearise your factors. For example, a log transform of factors may reduce your residuals. This could be used to argue that the relationship between factors and the response is not linear, since the transformed variables fit better. I hope that helps.

  • $\begingroup$ 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. $\endgroup$ Jan 27, 2015 at 17:05
  • $\begingroup$ The point is only to linearise the variable on the assumption that if it has an exponential effect, it would fit better. It would just be part of an argument that there is actually a non-linear component. $\endgroup$
    – Emir
    Jan 27, 2015 at 21:17
  • $\begingroup$ That is fine if there is only one transformation and that transformation is pre-specified. It's not so fine if the transformation is determined by analyzing the data. $\endgroup$ Jan 27, 2015 at 21:35

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.