topic Re: Relative risks between multinomial categories within a population in proc gee in Statistical Procedures

Relative risks between multinomial categories within a population in proc gee

quentinread — Thu, 22 Sep 2022 13:29:22 GMT

Disclaimer: I am really an R user but helping a client with their SAS code. I understand the basics of SAS but have trouble with customizing output to give me exactly the tests I want, hence this question!

I have a model I fit with proc gee that estimates the probability of insects being attracted to various types of plant tissue as a function of time. At any given time, an insect may choose between one of four tissue types or not respond at all. There were 4 trials, 10 replicates within each trial, and each replicate consisted of 20 insects. Within each replicate, at 4 time points, the number of insects on each of the tissue types was counted. So we have repeated measures data with a non-ordered multinomial outcome. I chose to use proc gee to fit this model. Following is the code, with a small subset of the data, enough to provide a reproducible example. Each row in the data is an individual insect at a single time point (interval), with the character column "response" having five possible values: Leaftip, Male, Female, Berries, and No response.

I was able to get the contrasts between different levels of "interval" (the time point) using the lsmeans statement, but I do not know how to get contrasts and the appropriate statistical test, to compare (1) the relative preference (relative risk) between the response categories within each interval and (2) the relative preference between the response categories, averaged across all four intervals. The hypothesis of interest to test would be, for example, "do insects have a higher preference for leaf tips relative to berries at time interval 1?" I found this note on calculating relative risks in a multinomial response model which seems to be doing something similar but I am not sure how to adapt the macro to work for this particular case. The examples given are for different procedures, not proc gee. Any help would be greatly appreciated!

data example;
	input trial rep interval response $12.; 
	datalines;
1 1 1 Leaftip
1 1 1 Leaftip
1 1 1 Leaftip
1 1 1 Leaftip
1 1 1 Male
1 1 1 Male
1 1 1 No response
1 1 1 No response
1 1 1 No response
1 1 1 No response
1 1 1 No response
1 1 1 No response
1 1 1 No response
1 1 1 No response
1 1 1 No response
1 1 1 No response
1 1 1 No response
1 1 1 No response
1 1 1 No response
1 1 1 No response
1 1 2 Leaftip
1 1 2 Leaftip
1 1 2 Leaftip
1 1 2 Leaftip
1 1 2 Leaftip
1 1 2 Leaftip
1 1 2 Leaftip
1 1 2 Leaftip
1 1 2 Leaftip
1 1 2 Male
1 1 2 Male
1 1 2 No response
1 1 2 No response
1 1 2 No response
1 1 2 No response
1 1 2 No response
1 1 2 No response
1 1 2 No response
1 1 2 No response
1 1 2 No response
1 1 3 Leaftip
1 1 3 Leaftip
1 1 3 Leaftip
1 1 3 Leaftip
1 1 3 Leaftip
1 1 3 Leaftip
1 1 3 Leaftip
1 1 3 Leaftip
1 1 3 Leaftip
1 1 3 Leaftip
1 1 3 Male
1 1 3 Male
1 1 3 Berries
1 1 3 No response
1 1 3 No response
1 1 3 No response
1 1 3 No response
1 1 3 No response
1 1 3 No response
1 1 3 No response
1 1 4 Leaftip
1 1 4 Leaftip
1 1 4 Leaftip
1 1 4 Leaftip
1 1 4 Leaftip
1 1 4 Leaftip
1 1 4 Leaftip
1 1 4 Leaftip
1 1 4 Leaftip
1 1 4 Leaftip
1 1 4 Leaftip
1 1 4 Leaftip
1 1 4 Leaftip
1 1 4 Leaftip
1 1 4 Leaftip
1 1 4 Leaftip
1 1 4 No response
1 1 4 No response
1 1 4 No response
1 1 4 No response
1 2 1 Leaftip
1 2 1 Leaftip
1 2 1 Leaftip
1 2 1 Leaftip
1 2 1 Male
1 2 1 Female
1 2 1 No response
1 2 1 No response
1 2 1 No response
1 2 1 No response
1 2 1 No response
1 2 1 No response
1 2 1 No response
1 2 1 No response
1 2 1 No response
1 2 1 No response
1 2 1 No response
1 2 1 No response
1 2 1 No response
1 2 1 No response
1 2 2 Leaftip
1 2 2 Leaftip
1 2 2 Leaftip
1 2 2 Leaftip
1 2 2 Leaftip
1 2 2 Leaftip
1 2 2 Berries
1 2 2 No response
1 2 2 No response
1 2 2 No response
1 2 2 No response
1 2 2 No response
1 2 2 No response
1 2 2 No response
1 2 2 No response
1 2 2 No response
1 2 2 No response
1 2 2 No response
1 2 2 No response
1 2 2 No response
1 2 3 Leaftip
1 2 3 Leaftip
1 2 3 Leaftip
1 2 3 Leaftip
1 2 3 Leaftip
1 2 3 Male
1 2 3 Berries
1 2 3 No response
1 2 3 No response
1 2 3 No response
1 2 3 No response
1 2 3 No response
1 2 3 No response
1 2 3 No response
1 2 3 No response
1 2 3 No response
1 2 3 No response
1 2 3 No response
1 2 3 No response
1 2 3 No response
1 2 4 Leaftip
1 2 4 Leaftip
1 2 4 Leaftip
1 2 4 Leaftip
1 2 4 Leaftip
1 2 4 Leaftip
1 2 4 Leaftip
1 2 4 Male
1 2 4 Male
1 2 4 No response
1 2 4 No response
1 2 4 No response
1 2 4 No response
1 2 4 No response
1 2 4 No response
1 2 4 No response
1 2 4 No response
1 2 4 No response
1 2 4 No response
1 2 4 No response
2 1 1 Leaftip
2 1 1 Leaftip
2 1 1 Leaftip
2 1 1 Berries
2 1 1 Berries
2 1 1 Berries
2 1 1 Berries
2 1 1 No response
2 1 1 No response
2 1 1 No response
2 1 1 No response
2 1 1 No response
2 1 1 No response
2 1 1 No response
2 1 1 No response
2 1 1 No response
2 1 1 No response
2 1 1 No response
2 1 1 No response
2 1 1 No response
2 1 2 Leaftip
2 1 2 Leaftip
2 1 2 Leaftip
2 1 2 Leaftip
2 1 2 Leaftip
2 1 2 Leaftip
2 1 2 Leaftip
2 1 2 Berries
2 1 2 Berries
2 1 2 Berries
2 1 2 Berries
2 1 2 No response
2 1 2 No response
2 1 2 No response
2 1 2 No response
2 1 2 No response
2 1 2 No response
2 1 2 No response
2 1 2 No response
2 1 2 No response
2 1 3 Leaftip
2 1 3 Leaftip
2 1 3 Leaftip
2 1 3 Leaftip
2 1 3 Leaftip
2 1 3 Leaftip
2 1 3 Leaftip
2 1 3 Berries
2 1 3 Berries
2 1 3 No response
2 1 3 No response
2 1 3 No response
2 1 3 No response
2 1 3 No response
2 1 3 No response
2 1 3 No response
2 1 3 No response
2 1 3 No response
2 1 3 No response
2 1 3 No response
2 1 4 Leaftip
2 1 4 Leaftip
2 1 4 Leaftip
2 1 4 Leaftip
2 1 4 Leaftip
2 1 4 Leaftip
2 1 4 Leaftip
2 1 4 Male
2 1 4 Female
2 1 4 Female
2 1 4 Berries
2 1 4 Berries
2 1 4 Berries
2 1 4 No response
2 1 4 No response
2 1 4 No response
2 1 4 No response
2 1 4 No response
2 1 4 No response
2 1 4 No response
2 2 1 Leaftip
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 1 No response
2 2 2 Leaftip
2 2 2 Leaftip
2 2 2 Leaftip
2 2 2 Leaftip
2 2 2 Male
2 2 2 Female
2 2 2 Female
2 2 2 Female
2 2 2 Female
2 2 2 Female
2 2 2 Berries
2 2 2 No response
2 2 2 No response
2 2 2 No response
2 2 2 No response
2 2 2 No response
2 2 2 No response
2 2 2 No response
2 2 2 No response
2 2 2 No response
2 2 3 Leaftip
2 2 3 Leaftip
2 2 3 Male
2 2 3 Female
2 2 3 Female
2 2 3 Female
2 2 3 Female
2 2 3 Female
2 2 3 No response
2 2 3 No response
2 2 3 No response
2 2 3 No response
2 2 3 No response
2 2 3 No response
2 2 3 No response
2 2 3 No response
2 2 3 No response
2 2 3 No response
2 2 3 No response
2 2 3 No response
2 2 4 Leaftip
2 2 4 Leaftip
2 2 4 Leaftip
2 2 4 Leaftip
2 2 4 Leaftip
2 2 4 Leaftip
2 2 4 Leaftip
2 2 4 Male
2 2 4 Male
2 2 4 Female
2 2 4 Berries
2 2 4 Berries
2 2 4 Berries
2 2 4 No response
2 2 4 No response
2 2 4 No response
2 2 4 No response
2 2 4 No response
2 2 4 No response
2 2 4 No response
;
run;

proc gee data=example;
   class response trial rep interval;
   model response=interval / dist=multinomial link=glogit;
   repeated subject=rep(trial);
   lsmeans interval / diff ilink oddsratio;
run;

Re: Relative risks between multinomial categories within a population in proc gee

sbxkoenk — Thu, 22 Sep 2022 16:15:57 GMT

Hello,

I think you need to :

Use the NLMeans macro OR
Use the NLEST/NLEstimate macro

@StatDave knows them well.

Do not try to change the source code of the macro. That is not needed.

Thanks,

Koen

Re: Relative risks between multinomial categories within a population in proc gee

quentinread — Thu, 22 Sep 2022 17:38:11 GMT

Hi Koen, yes, thanks very much for your response. I am hacking away at the NLMeans macro right now (don't worry I did not try to edit the source code, haha), and I think I am close to the correct solution. I believe I have all the correct contrasts. I have one other question: is it possible within the NLMeans macro to automatically adjust the p-values for multiple comparisons, or must I do this manually by calculating my desired adjustment to alpha and inputting that value to the alpha argument?

This is what I have so far:

proc gee data=example;
   class response(ref="No response") trial rep interval;
   model response=interval / dist=multinomial link=glogit;
   repeated subject=rep(trial);
   lsmeans interval / e diff ilink oddsratio;
   ods output lsmeans=gee1_lsmeans coef=gee1_coeffs;
   store out=gee1_mod;
run;

data tissue_contrasts;
        infile datalines missover;
        input label $40. k1-k20;
        set=1;
        datalines;
      Interval 1 P(leaftip)/P(male)      0  0  0  0    0  0  0  0    1  0  0  0   -1  0  0  0    0  0  0  0 
      Interval 1 P(leaftip)/P(female)    0  0  0  0   -1  0  0  0    1  0  0  0    0  0  0  0    0  0  0  0 
      Interval 1 P(leaftip)/P(berries)  -1  0  0  0    0  0  0  0    1  0  0  0    0  0  0  0    0  0  0  0 
      Interval 1 P(male)/P(female)       0  0  0  0   -1  0  0  0    0  0  0  0    1  0  0  0    0  0  0  0 
      Interval 1 P(male)/P(berries)     -1  0  0  0    0  0  0  0    0  0  0  0    1  0  0  0    0  0  0  0 
      Interval 1 P(female)/P(berries)   -1  0  0  0    1  0  0  0    0  0  0  0    0  0  0  0    0  0  0  0 
	  Interval 2 P(leaftip)/P(male)      0  0  0  0    0  0  0  0    0  1  0  0    0 -1  0  0    0  0  0  0 
      Interval 2 P(leaftip)/P(female)    0  0  0  0    0 -1  0  0    0  1  0  0    0  0  0  0    0  0  0  0 
      Interval 2 P(leaftip)/P(berries)   0 -1  0  0    0  0  0  0    0  1  0  0    0  0  0  0    0  0  0  0 
      Interval 2 P(male)/P(female)       0  0  0  0    0 -1  0  0    0  0  0  0    0  1  0  0    0  0  0  0 
      Interval 2 P(male)/P(berries)      0 -1  0  0    0  0  0  0    0  0  0  0    0  1  0  0    0  0  0  0 
      Interval 2 P(female)/P(berries)    0 -1  0  0    0  1  0  0    0  0  0  0    0  0  0  0    0  0  0  0
      Interval 3 P(leaftip)/P(male)      0  0  0  0    0  0  0  0    0  0  1  0    0  0 -1  0    0  0  0  0 
      Interval 3 P(leaftip)/P(female)    0  0  0  0    0  0 -1  0    0  0  1  0    0  0  0  0    0  0  0  0 
      Interval 3 P(leaftip)/P(berries)   0  0 -1  0    0  0  0  0    0  0  1  0    0  0  0  0    0  0  0  0 
      Interval 3 P(male)/P(female)       0  0  0  0    0  0 -1  0    0  0  0  0    0  0  1  0    0  0  0  0 
      Interval 3 P(male)/P(berries)      0  0 -1  0    0  0  0  0    0  0  0  0    0  0  1  0    0  0  0  0 
      Interval 4 P(female)/P(berries)    0  0 -1  0    0  0  1  0    0  0  0  0    0  0  0  0    0  0  0  0
      Interval 4 P(leaftip)/P(male)      0  0  0  0    0  0  0  0    0  0  0  1    0  0  0 -1    0  0  0  0 
      Interval 4 P(leaftip)/P(female)    0  0  0  0    0  0  0 -1    0  0  0  1    0  0  0  0    0  0  0  0 
      Interval 4 P(leaftip)/P(berries)   0  0  0 -1    0  0  0  0    0  0  0  1    0  0  0  0    0  0  0  0 
      Interval 4 P(male)/P(female)       0  0  0  0    0  0  0 -1    0  0  0  0    0  0  0  1    0  0  0  0 
      Interval 4 P(male)/P(berries)      0  0  0 -1    0  0  0  0    0  0  0  0    0  0  0  1    0  0  0  0 
      Interval 4 P(female)/P(berries)    0  0  0 -1    0  0  0  1    0  0  0  0    0  0  0  0    0  0  0  0
;
run;

%NLMeans(instore=gee1_mod, coef=gee1_coeffs, link=glogit, options=ratio, 
          contrasts=tissue_contrasts, null=1, title=Contrasts between tissue types within each interval);

Re: Relative risks between multinomial categories within a population in proc gee

sbxkoenk — Thu, 22 Sep 2022 19:04:40 GMT

@quentinread wrote:

I have one other question: is it possible within the NLMeans macro to automatically adjust the p-values for multiple comparisons, or must I do this manually by calculating my desired adjustment to alpha and inputting that value to the alpha argument?

No idea. 🤔☹️

Can you specify the alpha to use by the NLMeans macro? Probably yes (?).

And how would you calculate the alpha value to use? By using the simple Bonferroni correction? That might be too conservative.

There's a procedure for multiple comparisons and correction for inflation of the type I error rate.

It's called PROC MULTTEST.

SAS/STAT 15.2 User's Guide
The MULTTEST Procedure
https://go.documentation.sas.com/doc/en/pgmsascdc/9.4_3.5/statug/statug_multtest_toc.htm

Cheers,

Koen

Re: Relative risks between multinomial categories within a population in proc gee

StatDave — Thu, 22 Sep 2022 20:15:45 GMT

Good job putting together the contrasts data set for use in NLMeans. That looks correct. As far as adjusting the p-values for multiple testing... no, the macro does not incorporate that capability but it doesn't need to because it is easily done after the fact with PROC MULTTEST. NLMeans will create data set EST which contains the relative risk estimates. The p-values are in variable PR. So, you could use PROC MULTTEST as follows to adjust them using Holm's stepdown method. You could, of course, choose a different method, but Holm's method is more powerful than the ordinary Bonferroni method and is generally applicable.

proc multtest inpvalues(pr)=est holm; 
run;

Re: Relative risks between multinomial categories within a population in proc gee

quentinread — Mon, 26 Sep 2022 13:23:09 GMT

Thanks @StatDave , that is definitely what I was looking for. I hate to keep asking questions on this same thread, but do you have any idea how to use Hans-Peter Piepho's mult macro to get the multiple comparison letters for those adjusted p-values? I have tried to use it but it gives me nonsensical output.

Here is what I have tried: first, label the comparisons by interval so that the proc multtest does the Holm adjustment within each interval.

data my_nlmeans; set est;
	interval=ceil(_n_/6);
run;
proc multtest inpvalues(pr)=my_nlmeans holm; 
	by interval; id label; 
	ods output pvalues=mult_pvals;
run;

Next, manually label the p-value data object and the nlmeans data object with columns for the treatments being compared, as I believe is required by the mult macro. Then run the mult macro for an individual interval. This gives me "a" for every group even though some of the adjusted p-values are below 0.05. Maybe I am using it incorrectly?

data tissue_contrast_labels;
	infile datalines missover;
	input interval _interval response $10. _response $10.;
	datalines;
1 1 berries    female
1 1 berries    leaftip
1 1 berries    male
1 1 female     leaftip
1 1 female     male
1 1 leaftip    male
2 2 berries    female
2 2 berries    leaftip
2 2 berries    male
2 2 female     leaftip
2 2 female     male
2 2 leaftip    male
3 3 berries    female
3 3 berries    leaftip
3 3 berries    male
3 3 female     leaftip
3 3 female     male
3 3 leaftip    male
4 4 berries    female
4 4 berries    leaftip
4 4 berries    male
4 4 female     leaftip
4 4 female     male
4 4 leaftip    male
;
run;

data mult_pvals; merge mult_pvals tissue_contrast_labels; run;
data my_nlmeans; merge my_nlmeans tissue_contrast_labels; run;

%mult(diffs=mult_pvals, lsmeans=my_nlmeans, trt=response, p=stepdownbonferroni, by=interval, level=1);