Can proc plm lsmeans give me all k categories in glimmix multinomial, not k-1 categories?

[quentinread]

The following code gives me the lsmeans by treatment for two of the three ordered categories in this multinomial model, slight and moderate. I understand that the probabilities conditioned on treatment for the three categories must sum to 1. However I am looking for the most concise syntax to get me the estimates with their standard errors for all k categories, not just k-1 . I have scoured the documentation and even tried the nlmeans macro to no avail (I am not a great SAS programmer). Any help would be greatly appreciated!

data ratings;
input obs blk trt rating $ y ; 
datalines;
1 1 0 slight 1 
2 1 0 moderate 4 
3 1 0 severe 23 
4 1 1 slight 2 
5 1 1 moderate 7 
6 1 1 severe 23 
7 1 2 slight 4 
8 1 2 moderate 7 
9 1 2 severe 18 
10 2 3 slight 8 
11 2 3 moderate 8 
12 2 3 severe 5 
13 2 4 slight 14 
14 2 4 moderate 14 
15 2 4 severe 7 
16 2 5 slight 14 
17 2 5 moderate 11 
18 2 5 severe 1 
19 3 0 slight 3 
20 3 0 moderate 6 
21 3 0 severe 11 
22 3 1 slight 5 
23 3 1 moderate 9 
24 3 1 severe 15 
25 3 2 slight 9 
26 3 2 moderate 5 
27 3 2 severe 11 
28 4 3 slight 11 
29 4 3 moderate 9 
30 4 3 severe 6 
31 4 4 slight 14 
32 4 4 moderate 9 
33 4 4 severe 3 
34 4 5 slight 14 
35 4 5 moderate 6 
36 4 5 severe 2 
37 5 0 slight 0 
38 5 0 moderate 6 
39 5 0 severe 13 
40 5 1 slight 4 
41 5 1 moderate 7 
42 5 1 severe 16 
43 5 2 slight 12 
44 5 2 moderate 12 
45 5 2 severe 4 
46 6 3 slight 6 
47 6 3 moderate 14 
48 6 3 severe 14 
49 6 4 slight 10 
50 6 4 moderate 13 
51 6 4 severe 8 
52 6 5 slight 16 
53 6 5 moderate 9 
54 6 5 severe 4 
55 7 0 slight 3 
56 7 0 moderate 9 
57 7 0 severe 17 
58 7 1 slight 2 
59 7 1 moderate 6 
60 7 1 severe 22 
61 7 2 slight 3 
62 7 2 moderate 10 
63 7 2 severe 9 
64 8 3 slight 8 
65 8 3 moderate 13 
66 8 3 severe 6 
67 8 4 slight 13 
68 8 4 moderate 10 
69 8 4 severe 1 
70 8 5 slight 17 
71 8 5 moderate 7 
72 8 5 severe 1 
73 9 0 slight 2 
74 9 0 moderate 3 
75 9 0 severe 23 
76 9 1 slight 1 
77 9 1 moderate 5 
78 9 1 severe 19 
79 9 2 slight 2 
80 9 2 moderate 13 
81 9 2 severe 4 
82 10 3 slight 4 
83 10 3 moderate 2 
84 10 3 severe 8 
85 10 4 slight 14 
86 10 4 moderate 5 
87 10 4 severe 11 
88 10 5 slight 12 
89 10 5 moderate 4 
90 10 5 severe 7 
;


proc glimmix data=ratings method=laplace;
  class blk trt (ref = "0");
  model rating(order=data)=trt / dist=multinomial;
  random intercept / subject=blk; 	
  freq y;
  store fit;
 run;

 proc plm restore=fit;
 	lsmeans trt / ilink cl;
run;

[quentinread]

By the way I understand that the lsmeans output for proc plm is giving me P(slight) and P(slight)+P(moderate) for each treatment. I understand that P(slight)+P(moderate)+P(severe) = 1. But what I really want is the estimate and standard error for P(slight), P(moderate), and P(severe) for each treatment and I am curious if this is easily achievable in proc plm. 

[quentinread]

I have managed to calculate the standard errors by exporting the lsmeans and their covariance matrix, importing them into R, and using the delta method with the function msm::deltamethod(). But it would still be great if I could find a way to do this in SAS!

[StatDave · Accepted solution]

This is easily done with the current release (v3.0) of the NLMeans macro. See the multinomial example in the Results tab of the macro documentation which helps to explain the syntax though it is applied to a nominal, not ordinal, multinomial model. For your example, the probability of SLIGHT in treatment 1 is the first mean from the LSMEANS statement, so it is referred to as mu1  in the macro. Similarly for the remaining probabilities. So, the estimate, standard error (using the delta method), and confidence interval are provided by the following macro call after adding the E option in your LSMEANS statement and also adding this ODS statements in your PLM step: ods output coef=coeffs;

   %nlmeans(instore=fit, coef=coeffs, link=clogit, f=mu1)

You can get estimates of the treatment 1 probabilities in all three levels by subtraction using the fdata= option and appropriate data set similar to what you see in the example in the macro documentation.

   data fd; 
      set=1;
      length label f $32767; 
      infile datalines delimiter=',';
      input label f; 
      datalines;
   P(slight trt1),   mu1 
   P(moderate trt1),  mu7-mu1
   P(severe trt1),  1-mu7
   ;
   %nlmeans(instore=fit, coef=coeffs, link=clogit, fdata=fd)

[quentinread]

Thank you this is extremely helpful! Just one thing I noticed. The parameter df must be explicitly provided, to get the same confidence intervals as the original lsmeans. According to the documentation, if df is omitted, NLMEANS will use large-sample Wald statistics. In this example it hardly makes any difference because df=768, but I wanted to point that out for future reference.