Adjustment for multiple comparisons with LSMeans in multinomial generalized logit model

[Kastchei]

Hello,

I'm using PROC GEE to analyze a multinomial generalized logit model (4-category outcome = 3 logit models).  I am running LSMeans to get all pairwise differences with an adjustment for multiple comparisons.  I notice that the adjustments are made viewing the 3 logit models separately rather than jointly as a system; thus the adjustment is lot smaller than I expected.  For example, if I have 6 pairwise comparisons per logit model and 3 logit models, I have a total of 18 pairwise comparisons.  I figured the adjustment would be for 18 comparisons, but the adjustment only seems to be for 6 comparisons.  For example, if using Bonferroni*, adjusted p-value = p-value * 6, rather than p-value * 18.

Is this ok to do?  Can the 3 logit models really be viewed independently so I don't have to adjust for comparisons made in the other models?

This difference can be seen using LSMEstimate and coding all pairwise comparisons.  If I use category = joint, I get the expected adjustment for 18 comparisons, but if I use category = separate, then I get the same results as LSMeans, adjusted for only 6 comparisons.

(* Note, I am not actually using Bonferroni, but it was the easiest to use as an example.)

Warm regards,

Michael

proc gee data = cogFunc._07_Model descending;
    class ID outcomeClassN exposureC (ref = '1 = Low');

    model outcomeClassN = exposureC / dist = multinomial link = gLogit type3 wald;
    repeated subject = ID;

    lsMeans exposureC / diff adjust = bon;
    lsmEstimate exposureC 'Very high vs. Low'           0  0  1 -1,
                          'High vs. Low'                0  1  0 -1,
                          'Intermediate vs. Low'        1  0  0 -1,
                          'Very high vs. Intermediate' -1  0  1  0,
                          'High vs. Intermediate'      -1  1  0  0,
                          'Very high vs. High'          0 -1  1  0
        / category = joint adjust = bon;
run;

[PGStats]

You can combine as many test p-values for correction as you want with proc multtest. Just output the p-values to a dataset using ODS OUTPUT and refer to that dataset with the INPVALUES= option of proc multtest.