Visualizing an interaction between 2 continuous predictors w/confidence bands (Proc Mixed model)

[Poppi21]

Hi, 

I have data from an experience sampling study (daily observations nested in individuals) and have been using proc mixed to analyze this multilevel data. I have an interaction that I want to plot: between a continuous (day level) predictor (X) and a continuous (day level) moderator (M). 

I would like to graph the interaction w/confidence bands in sas (I previously graphed the interaction in excel +/- 1 standard deviation for both X and M (so plotting at :low-low, low-high, high-low, high-high).

I used the "store" statement in my model and then Proc plm w/effectplot but haven't been able to get my intended result. 

Can anyone help with this?

Thank you!

[Ksharp]

Where is Original Post ?

@Rick_SAS  has written a blog about this topic .

https://blogs.sas.com/content/iml/2019/12/05/longitudinal-data-mixed-model.html

https://blogs.sas.com/content/iml/2018/12/19/visualize-mixed-model.html

https://blogs.sas.com/content/iml/2019/05/30/visualize-interaction-effects-regression.html

[Poppi21]

Hi Thank you,
I read these pages and actually posted the question there as well Rick suggested I post it on communities:-)

[JacobSimonsen]

I think you are right by using "store" and proc plm. Can you provide your program?

[Poppi21]

Hi Thank you!
After storing the Mixed model used the plm program below:
* XT1 = continuous predictor and MT1 = continuous moderator.

proc plm source=MixedModel1;                 /* or use restore statement */
effectplot slicefit(x=XT1 sliceby=MT1=4.9273927 5.398148)/ clm
run; 

However, I'm actually trying to plot it at the following 4 points +/- 1 SD for XT1 and MT1:

	Low XT1	High XT1
Low MT1	4.927392	5.95138
High MT1	5.398148	6.886211

And to ultimately produce the below graph with confidence bands (this one was created in excel) 

Any help would be greatly appreciated!

Thank you!

[Rick_SAS]

To clarify, you want the predicted value and prediction intervals for the four points in the table? And XT1 and MT1 are both fixed effects, right?

Please provide the SAS code that specifies the model that you fit.

[Poppi21]

Yes exactly, 

This is the code for the model: 

PROC MIXED COVTEST DATA=Work.EsmDataNew METHOD=reml cl;
 CLASS ID Day;
 MODEL Y=Day_midc
Work_AP_c  
AM_motivation_pc
XT1_pm 
XT1_pc 
MT1_pm
MT1_pc
XT1_pm *MT1_pm
XT1_pc *MT1_pm
XT1_pm *MT1_pc
XT1_pc *MT1_pc/S; 
 RANDOM intercept XT1_pc /SUBJECT=id TYPE=UN GCORR cl;
 REPEATED Day/subject=id type=ar(1);
 store out=MixedModel1;                       /* create item store */
run;

[*pc represents the within person component (person centered) and _pm represents between persons components (grand mean centered)].

Thanks again!

[Rick_SAS]

I could be wrong, but I don't think the EFFECTPLOT statement supports creating an INTERACTION plot for arbitrary values of a continuous variable (at mean +/- stdDev). You can get a SLICEFIT or CONTOUR plot, which both display more information for continuous variables. Maybe plot BY day?

Anyway, If you really want the graph that you display, I think you will have to create your own scoring set in which you specify the mean value (or some other reference value)  for most of the covariates. I've shown how to do this in two articles:

Use the missing value trick to score and visualize a regression model

How to create a sliced fit plot in SAS

Both of those articles assume that you want to see a sliced fit plot. To get the interaction plot that you describe, try something like this:

/* sample model with interaction terms */
proc glm data=sashelp.cars;
model mpg_city = weight | horsepower ;
store glmModel;
run;

/* form mean +/- std dev and create the scoring data set */
/* find mean and std dev for continuous covariates */
proc means data=sashelp.cars;
var weight horsepower;
output out=MeanOut mean  =weight_mean horsepower_mean 
                   stddev=weight_std  horsepower_std;
run;

/* form the scoring data set */
data ScoreData;
set MeanOut;
wmin = weight_mean - weight_std;
wmax = weight_mean + weight_std;
hmin = horsepower_mean - horsepower_std;
hmax = horsepower_mean + horsepower_std;
xLabel = "Low XT1 "; cLabel="Low MT1 "; weight=wmin; horsepower = hmin; output;
xLabel = "Low XT1 "; cLabel="High MT1"; weight=wmin; horsepower = hmax; output;
xLabel = "High XT1"; cLabel="Low MT1 "; weight=wmax; horsepower = hmin; output;
xLabel = "High XT1"; cLabel="High MT1"; weight=wmax; horsepower = hmax; output;
run;

/* evaluate model on scoring data set */
proc plm restore=glmModel;
score data=ScoreData out=ScoreOut PREDICTED LCLM UCLM;
run;

proc print data=ScoreOut;run;

proc sgplot data=ScoreOut;
series x=xLabel y=Predicted / group=cLabel;
scatter x=xLabel y=Predicted / group=cLabel yerrorlower=LCLM yerrorupper=UCLM;
xaxis offsetmin=0.2 offsetmax=0.2;
run;

[JacobSimonsen]

Hi again,

With some simpler simulated I was able to plot the effect with the proc plm using the effectplot with slicefit option. 

*simulate a dataset with predictors c and x, and outcome y;
data simulation;
  do i=1 to 1000;
    c=int(rand('uniform',0,5));
    x=int(rand('uniform',0,5));
	y=rand('normal',(1/2-c)*x+x*x+2,1);
    output; 
  end;
run; 

*estimate the model and save the model to "mymodel";
proc mixed data=simulation;
  model y= c*x c c*x*x;
  store out=mymodel;
run;

*plot the prediction with confidence limits;
ods graphics;
ods hmtl;
proc plm restore=mymodel;
  effectplot slicefit(x=x sliceby= c=2,3) /clm;
run;

 Since you have many more covariates in your model, all the other covariates will be set at some level. I hope you find it usefull.

[Ksharp]

Jacob, there is interaction plot in EFFECTPLOT like:
proc plm restore=mymodel;
effectplot interaction(x=...........) /clm;
run;