Plotting interaction effects on a nonlinear relationship

[Tadhg Ryan Charleton]

Hi – I will really appreciate any thoughts here. I expect (hope) to produce a plot that looks something like the attached.

I have been attempting variations of two approaches (where x1 is my continuous IV and x2 is my nonnegative integer moderator).

The first approach is to graph three qfit lines onto the same plot:
nbreg y c.x1##c.x1 c.x1##c.x2 …
predict yhat, xb … {which plots fitted values at the means}
*in theory I would create three fitted values for yhat at the three values of x2 I want to graph. And from there generate the qfit plot. Yet I'm not sure about this because predict does not allow the ‘at’ option and I haven’t found a workaround.

The second approach is marginsplot:
nbreg y c.x1##c.x1 c.x1##c.x2 …
margins, dydx(x1) at(x2 =0 0.33 1.31)
marginsplot
*I am getting the following error for this (even though I have specified x2 as a factor variable)
variable '0.33' not found in list of covariates

(The issue is complicated slightly because x2 is positively skewed. Therefore, instead of using + and – SD points, I am using min, mean, and mean +SD.)

Huge thanks in advance for any feedback!

[Phil Bromiley]

You'll increase your chances of a useful answer by following the FAQ on asking questions – provide Stata code in code delimiters, readable Stata output, and sample data using dataex.

You can probably do this with predictnl in margins. Alternatively, you can calculate it yourself using the parameter values.

[Tadhg Ryan Charleton]

Got it with margins. Thank you Phil.

[Clyde Schechter]

With respect to #1, I see a couple of additional problems.

margins, dydx(x1) at(x2 =0 0.33 1.31)

is producing the error message you got because the numlist 0 0.33 1.31 needs to be wrapped in its own parentheses, thus:

margins, dydx(x1) at(x2 =(0 0.33 1.31))

Also, you do not show your full -nbreg- command, but what you do show so far is an incomplete and mis-specified model.
It is wrong to have a quadratic term in x1 and then have an interaction between x1 and x2 without also having the interaction between x1 and the square of x2. So unless c.x1#c.x2#c.x2 appears somewhere in the ... portion of your model, you are just putting garbage in and you will just get garbage out.

I do not understand why so many people list out "main" effects and partial interactions separately, when the ## operator in factor-variable notation makes it shorter and simpler, and guarantees that you will never leave out something that needs to be there. So:

nbreg y c.x1##c.x1##c.x2 ...

will get it right and is simpler to write.