How do I get positive predicted probabilities in GLIMMIX?

[plf515]

I am using SAS 9.4 on Windows 10

In PROC LOGISTIC, if you have a dependent variable with multiple levels and runs ordinal regression with cumulative logits, you can get the predicted probability of being at each level of the dependent variable for various combinations of the independent variables by using PREDPROBS = I.

How can I do the same for PROC GLIMMIX? I have this code:

title "Predicted Q 1";
proc glimmix data = jenny.allv2;
 class group resnumber h_c ressex pgy phase;
 model Q1_Make_at_ease = group pgy h_c phase group*phase ressex/s link = cumlogit dist = multinomial or;
 random int/sub = resnumber;
 output out = q1out pred(ilink) = p;
run;

Where Q1_Make_At_Ease has 3 levels. But the predicted value is a single probability, which is not particularly intuitive, even to me, much less to my clients. How can I get the predicted probability of being at each level of Q1..... ?

[sld · Accepted solution]

There may be more elegant ways, but this is how I do it:

/*  Example data from
    https://support.sas.com/documentation/cdl/en/statug/63962/HTML/default/viewer.htm#statug_glimmix_sect025.htm 
    modified by making zero counts be ones for tidiness */
data foot_mv;
  input yr b1 b2 b3 k1 k2 k3;
  sire = _n_;
  datalines;
  1  1  0  0  52 25  1
  1  1  0  0  49 17  1
  1  1  0  0  50 13  1
  1  1  0  0  42  9  1
  1  1  0  0  74 15  1
  1  1  0  0  54  8  1
  1  1  0  0  96 12  1
  1 -1  1  0  57 52  9
  1 -1  1  0  55 27  5
  1 -1  1  0  70 36  4
  1 -1  1  0  70 37  3
  1 -1  1  0  82 21  1
  1 -1  1  0  75 19  1
  1 -1 -1  0  17 12 10
  1 -1 -1  0  13 23  3
  1 -1 -1  0  21 17  3
 -1  0  0  1  37 41 23
 -1  0  0  1  47 24 12
 -1  0  0  1  46 25  9
 -1  0  0  1  79 32 11
 -1  0  0  1  50 23  5
 -1  0  0  1  63 18  8
 -1  0  0 -1  30 20  9
 -1  0  0 -1  31 33  3
 -1  0  0 -1  28 18  4
 -1  0  0 -1  42 27  4
 -1  0  0 -1  35 22  2
 -1  0  0 -1  33 18  3
 -1  0  0 -1  35 17  4
 -1  0  0 -1  26 13  2
 -1  0  0 -1  37 15  2
 -1  0  0 -1  36 14  1
 -1  0  0 -1  63 20  3
 -1  0  0 -1  41  8  1
;
data footshape; set foot_mv;
  array k{3};
  do Shape = 1 to 3;
     count = k{Shape};
     output;
  end;
  drop k:;
run;
/*  Fit model */
proc glimmix data=footshape method=quad;
   class sire;
   model Shape = yr b1 b2 b3 / s link=cumprobit dist=multinomial;
   random int / sub=sire;
   /* save BLUEs (predmupa) and BLUPs (predmu) */
   output out=glmm pred(noblup ilink)=predmupa pred(ilink)=predmu;
   freq count;
run;

/*  Process BLUEs */
proc sort data=glmm out=blue nodupkeys;
    by yr b1 b2 b3 _level_;
    run;
proc transpose data=blue out=blue_wide prefix=shape;
    by yr b1 b2 b3 ;
    id _level_;
    var predmupa;
run;
data blue_wide;
    set blue_wide;
    prob1 = shape1;
    prob2 = shape2 - shape1;
    prob3 = 1 - shape2;
    drop _name_ _label_;
    run;

/*  Process BLUPs */
proc sort data=glmm out=blup nodupkeys;
    by yr b1 b2 b3 sire _level_;
    run;
proc transpose data=blup out=blup_wide prefix=shape;
    by yr b1 b2 b3 sire;
    id _level_;
    var predmu;
run;
data blup_wide;
    set blup_wide;
    prob1 = shape1;
    prob2 = shape2 - shape1;
    prob3 = 1 - shape2;
    drop _name_ _label_;
    run;

[plf515]

Thanks.  That works. I wonder why SAS doesn't make this available directly?

[sld]

You could post the idea to https://communities.sas.com/t5/SASware-Ballot-Ideas/idb-p/sas_ideas 🙂