How to graph (fixed effects) interaction for mixed model

[Fred Friesen]

Hi there, I am relatively new to STATA, and trying to learn to use it to analyze some of my longitudinal count data (Poisson distributed outcome, autoregressive covariance structure, some overdispersion).

Part 1 of my question has to do with how to write STATA code for the random effects and interaction term. For my model, based on previous analysis with another program (SAS), I know I have a significant random effect for the intercept, time (time0), and var2. This is as far as I have gotten creating the code, which does not yet include the random effect of time0 or var2 (I used the STATAMULTILEVEL MIXED EFFECTS REFERENC EMANUAL RELEASE 13):

 gllamm DV time0 var1 var2 var3, family(poisson) link(log) i(id) eform adapt

-where DV = my dependent variable (a count), time0 = day (day 1 to 13), var1 = level-1 predictor (continuous), var2 = level-1 predictor (continuous), var3 = level-2 predictor (continuous)

In addition, I believe there is a var1 x var2 interaction (level-1 interaction), but if I input var1##var2, I receive the error "interactions not allowed". I am guessing this may mean that gllamm doesn't allow interactions? How can I run this or a similar appropriate model that allows for random effects and interaction terms?

Part 2: I would also like to graph the interaction (predicted values), e.g., using a line graph or similar, preferably back in the units of the original scale (not in logarithmic form or odds ratios/incident rates), with the DV on the y-axis, var1 on the x-axis, and separate lines for var2 (is on a 7-point scale, so could be 3 lines, one at the mean and ones at +/1 SD from the mean). Any suggestions for how I could produce this graph would be greatly appreciated!

[Clyde Schechter]

If you are running a current version of Stata, you can do a mixed effects Poisson regression using the -mepoisson- command instead of -gllamm-. -mepoisson- will work with interactions on the bottom level. Note that you have to specify it as c.var1##c.var2. Without the c. prefixes, the default interpretation of the interaction term is that the variables are discrete.

To get the graph, you have to first pack some values of var1 and var2 that you are interested in. For the sake of demonstrating the code, I'll assume that an interesting set of values for var1 is 1, 3, 5, 7, and 9 and that an interesting set of values for var2 is 2.5 4 5.5

mepoisson DV time0 c.var1##c.var2 var3 || group:
margins, at(var1 = (1 3 5 7 9) var2 = (2.5 4 5.5))
marginsplot

Notes: In the above code, group is a stand-in for whatever is the variable that you are taking random intercepts for. The -marginsplot- command accepts nearly all -graph twoway- options, so you can modify the appearance of the graph to your taste. I think the plot you get from the above command will give you var1 on the horizontal axis, and separate curves for each value of var2. If I have that wrong and it's the other way around, just specify -xdimension(x1)- as an option in the -marginsplot- command.

[Fred Friesen]

Clyde, amazing, thank you so much for this clear and detailed answer. I tried it and it worked great! Apologies also for my "STATA" faux pas, I understand now that it should be Stata. If there is anything I can do to return the favor for the Stata community let me know.