LS means on glimmix multinomial dist

[apereira]

Hello everyone,

I'd like help with my code once the following error has happened: "ERROR: Least-squares means are not available for the multinomial distribution."

Data body_score;
infile datalines dlm="09"x;
input ID block room trtm$ score;
datalines;

. . . . . 

. . . . 

;
Run;

--

proc glimmix data=body_score method=quad;
class trtm room block ;
model score=trtm*room/dist=multinomial link=cumlogit;
random intercept/subject=block(room);
lsmeans trtm*room/adjust=tukey ilink;
Run;

Any thoughts on what could be happening?

Thanks in advance

[PaigeMiller]

The error message seems to be clear. You can't compute LSMEANS for multinomial models.

[apereira]

Hi Paige, thanks for your reply.

I'm a beginner on SAS, and I recently participated in a workshop where LS means was used in a multinomial distribution of GLIMMIX, and I saw some answers in other topics that were also used it as well.

Do you have any thoughts on how I get the means if I remove LS means?

[PaigeMiller]

Perhaps those workshops have a replacement method. The SAS documentation is very clear (as was the error message), you cannot compute LSMEANS on multinomial models.

LSMEANS Statement

• LSMEANS fixed-effects </ options>;

The LSMEANS statement computes least squares means (LS-means) of fixed effects. As in the GLM and the MIXED procedures, LS-means are predicted population margins—that is, they estimate the marginal means over a balanced population. In a sense, LS-means are to unbalanced designs as class and subclass arithmetic means are to balanced designs. The  matrix constructed to compute them is the same as the  matrix formed in PROC GLM; however, the standard errors are adjusted for the covariance parameters in the model. Least squares means computations are not supported for multinomial models.

[StatDave · Accepted solution]

While LS-means for the ordinal logistic model are not available using PROC GLIMMIX, you could instead fit an ordinal logistic model using the GEE (Generalized Estimating Equations) method. That can be done in PROC GEE using the REPEATED statement instead of a RANDOM statement. The same can be done in PROC GENMOD but PROC GEE is now the recommended procedure for fitting GEE models. For example, the following fits a cumulative logit model to repeated measurements on the subjects (SUBJID) and gets LS-means for one of the predictors.

proc gee;
class subjid a b;
model y=a b / dist=mult;
repeated subject=subjid;
lsmeans a / ilink; 
run;

[jiltao]

While PROC GEE proposed by @StatDave_sas is a straghtforward way to get LSMEANS for your model, you might use the ESTIMATE statement in PROC GLIMMIX to compute the LSMEANS. The exact syntax for the ESTIMATE statement depends on your data and model, but the following example code might be helpful to you --

data test;
   do id=1 to 1000;
      trt=ceil(ranuni(123)*2);
      do time=1 to 3;
      n=rantbl(123,.10,.20,.25,.30,.15);
      output;
      end;
   end;
run;

proc glimmix data=test;
   class trt time;
   model n=trt time trt*time / dist=multinomial link=clogit s;
   estimate 'p(n=1) for trt=1' int 1 trt 1 / ilink e;
   estimate 'p(n=1) for trt=2' int 1 trt 0 1 / ilink e;
   estimate 'p(n=1,2) for trt=1' int 0 1 trt 1 time .33333 .33333 .33333
                                 trt*time .33333 .33333 .33333/ ilink e;
   estimate 'p(n=1,2) for trt=2' int 0 1 trt 0 1 time .33333 .33333 .33333
                                 trt*time 0 0 0 .33333 .33333 .33333/ ilink e;
   estimate 'p(n=1,2) for trt=1,time=2' int 0 1 trt 1 time 0 1 0 trt*time 0 1 0 0 0 0 / ilink;
   estimate 'p(n=1,2) for trt=2,time=2' int 0 1 trt 0 1 time 0 1 0 trt*time 0 0 0 0 1 0 / ilink;
   estimate 'difference' trt -1 1 time 0 0 0 trt*time 0 -1 0 0 1 0;
   estimate 'p(n=1,2,3) for trt=1,time=2' int 0 0 1 trt 1 time 0 1 0 trt*time 0 1 0 0 0 0 / ilink;
   estimate 'p(n=1,2,3) for trt=2,time=2' int 0 0 1 trt 0 1 time 0 1 0 trt*time 0 0 0 0 1 0 / ilink;
   estimate 'difference' trt -1 1 time 0 0 0 trt*time 0 -1 0 0 1 0;
   estimate 'p(n=1,2) for trt=2' int 0 1 trt 0 1 time .33333 .33333 .33333
                                 trt*time 0 0 0 .33333 .33333 .33333/ ilink;
   estimate 'p(n=1,2,3) for trt=1' int 0 0 1 trt 1 time .33333 .33333 .33333
                                 trt*time .33333 .33333 .33333/ ilink;
   estimate 'p(n=1,2,3) for trt=2' int 0 0 1 trt 0 1 time .33333 .33333 .33333
                                 trt*time 0 0 0 .33333 .33333 .33333/ ilink;
run;

Jill