How do I use Effectplot for Mixed Model with Glimmix with Random Effect

[Sascoder]

For plotting regression results, I am loving the use of proc plm, effectplot and the ability to store results for use in sgplot, etc.

However, I am now stumped.

Is there an option to use effectplot for a mixed model where I want to plot the effects by the random effect?  It gives me the error that sliceby must be a main effect. 

Thanks!

Example code

proc glimmix data=tst IC=Q method=laplace;  
   class re_var;
   model resp (event="1") = var1 var2 /dist=binary link=logit solution ;   
    random  int var1 var1 /subject=re_var cl type = un;  
    store md1;

ods graphics / width=3in height=3in border = off imagename="tst2";
proc plm restore=md1 ;
  effectplot slicefit(x=var1 sliceby=re_var) / clm;
run;
ods graphics off;

[Rick_SAS]

Yes, the EFFECTPLOT statement is for main effects.

It sounds like you want to create an overlay (or panel) of predicted values for each subject (which might be a patient, a doctor, or a school). Here is the main idea:

1. Create a set of data on which to score the model for each level of re_var. This might contain evenly spaced points of var1 and mean (or median) values of var2 for each level or re_var.

2. Use the missing value trick or the SCORE statement in PROC PLM to evaluate the model on the scoring data.

3. Use PROC SGPLOT to overlay the predicted curves.

The following articles show the basic ideas:

Visualize multivariate regression models by slicing continuous variables

How to create a sliced fit plot in SAS

I have written about this topic several times, so also search for keywords and use

   site:blogs.sas.com/content/iml

to restrict the internet search to my blog.

[Sascoder]

Thanks very much!

The score option is new to me, so great tip.

Still, it seems like this is still generating only fixed effect.  I don't see any option in the score command to ask for predictions accounting for the random effect.

I created the x-axis dataset that has values for all variables in the model and for each level of re_all (random effect).   (I won't show this code for simplicity here).

Here is my plm code:

proc plm restore=an.md1 ; 
  score data=xvals out=tst predicted lclm uclm / ilink;
run;

and here is a sample of the output dataset tst for a subset of the data to compare whether I am getting the random slope. You will see that the predicted values change for the fixed categorical variable var1 but not for the random factor re_all.  And the plots I did show this also.

  Obs    re_all  var1          var2        Predicted      LCLM        UCLM

   31       A        0        -0.73061       0.012849    0.012849    0.012849
  231       A        1        -0.73061       0.079731    0.079731    0.079731
  431       B        0        -0.73061       0.012849    0.012849    0.012849
  631       B        1        -0.73061       0.079731    0.079731    0.079731

Any suggestions?  What am I missing? 

I think I could figure out how to hard code the slope (adding random and fixed effect betas) but I don't know how to get the confidence intervals.

And the plotting part I know how to do.  It is how to get the correct predicted values and confidence intervals to plot.

Thanks again!

[Rick_SAS]

Yes, the SCORE statement evaluates the fixed-effects model, which is also called the "marginal model." Sorry to mislead you, 

Mixed models have two different kinds of predictors. The marginal predictor does not incorporate the random effects, the other predictor (called "BLUP") does incorporate the random effects. See the last section of the article "Visualize a mixed model that has repeated measures or random coefficients" for a discussion.

Generalized mixed models enable you to make predictions on the scale of the data on the scale of the linear predictor. See the article "Predicted values in generalized linear models: The ILINK option in SAS."

Assuming you want the BLUP predictions on the scale of the data, I think you need to use the OUTPUT statement in PROC GLMMIX to get the predictions you want. I think the code will look like this (untested):

output out=out1 pred(blup ilink)=predicted lcl(blup ilink)=lower ucl(blup ilink)=upper;

[Sascoder]

So I think I figured this out.  Please let me know if this makes sense.

1. I created a dataset for x axis values for all levels of my random effect (essentially what I want to plot) and with appropriate values of all variables in the model (e.g., means for other continuous variables) AND with Response = . (missing)

2. I appended this to my real dataset

3. Run the model with the line you suggested

        output out=out1 pred(blup ilink)=predicted lcl(blup ilink)=lower ucl(blup ilink)=upper;

   This provides the predicted values for range of values I want to predict (and makes a nicer graph) and includes random component.

4. I then use sgplot to plot the predictions vs. x-axis values for the appended data.

Honestly, I was really hoping there was another way with proc plm because some of my models take a really long time to run, so if I forget a step, it is really handy to have a backup with plm so I don't have to rerun the whole model again. 

Do you happen to know if there is a way to hard code the confidence intervals if I have a stored covariance matrix?  I have seen examples like this in R, but wouldn't quite be able to copy that into sas with confidence.  I tend to store a lot of pieces of the output just because the models take so long to run.

(Apologies, I understand this stuff conceptually (and how to program) more than the specific math or always have the correct terminology.  I am learning tho.....)

Thanks much!

[Rick_SAS]

Yes, that is an accurate summary of the steps that I suggest.

Regarding obtaining the covariance matrix of the betas, I have written an article about this topic for fixed-effect models:

"3 ways to obtain the Hessian at the MLE solution for a regression model"

[Sascoder · Accepted solution]

Thanks very much!!

I'll check out your article on hard-coding the SE calcs and post separately if I have questions.

Thanks again!