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.
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.
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)
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.
anova is not as comprehensive.
$\endgroup$
Commented
Jan 27, 2015 at 4:46
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
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.
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.
