proc genmod poisson using bonferroni corrections on 95% confidence interval

[Cruise]

Hi SAS enthusiasts,

I'm running poisson regression analysis for 40 different cancer types where association of each cancer type was tested against the 6 distinct predictors such as: PTRAF, PNPL, PRMP, PTSDF, NPL_CNT and TSDF_CNT.  At the end of the macro, I have 40*6=240 separate SAS annotated outputs. I'd like to do bonferroni correction to the 95% CI. 

Do you know how to reflect bonferron correction in the code?  Shall I manually manipulate Lower and Upper Wald statistics?

Any hints, insights, brainstorming or direct help will be apprceiated!

Thanks  for your time indeed!

%SYMDEL; 
%macro models (parest,ejvar,parej);
ODS OUTPUT   ParameterEstimates=&parest;
PROC GENMOD DATA=MYDATA;
BY CANCER_TYPE NOTSORTED;
CLASS AGECAT(REF='10')/PARAM=REF;
MODEL N_BLOCK = &ejvar AGECAT/ TYPE3 DIST=POISSON LINK=LOG OFFSET=LN;
RUN;

ODS RTF CLOSE;
DATA &parej; SET &parest;
RR= EXP(ESTIMATE);
LOW_RR=EXP(LowerWaldCL);
UP_RR=EXP(UpperWaldCL); 
RUN;

%mend models;
%MODELS(PAREST_PTRAF , X_PTRAF , PAREJ_PTRAF );
%MODELS(PAREST_PNPL  , X_PNPL  , PAREJ_PNPL  );
%MODELS(PAREST_PRMP  , X_PRMP  , PAREJ_PRMP  );
%MODELS(PAREST_PTSDF , X_PTSDF , PAREJ_PTSDF );
%MODELS(PAREST_NPL   , X_NPL_CNT , PAREJ_NPL);
%MODELS(PAREST_TSDF  , X_TSDF_CNT, PAREJ_TSDF);

[Rick_SAS]

Maybe I am misunderstanding, but if you want to compare the effects of the cancer types, wouldn't you want to put

CANCER_TYPE on the CLASS statement and include it in the model? That enables you to compare the different types, estimate difference, adjust for multiplicity, and so on.

As a reminder, here is a discussion of the difference between the BY statement and the CLASS statement for statistical anal...

[Cruise]

Hi @Rick_SAS, yes, this is the web resource that I consulted with before I develop my models. Thanks for making this resource available which helped me tremendously like any other publications of yours. 

In this case, my cancer types are independent from each other and I'm not comparing them each other. Each cancer types are subset data-sets on the given cancer and I'm using BY statement with notsorted statement to run multiple tests. 

With that said, do you have nay suggestions how to do Bonferroni corrections on my 95% CI intervals?

[Rick_SAS · Accepted solution]

Each BY group analysis is independent, so the cancer type is not part of the analysis.

So are you asking how to estimate differences between levels of the AGECAT variable? That's the only categorical variable in your model.

You can use the LSMEANS statement to compare the levels of AGECAT. You can use the ILINK option to get the estimates on the data scale and CL to request confidence intervals. You can use the ADJUST= option to adjust for multiple comparisons by using Bonferroni or other methods. Here's an example for a logistic model, but the syntax for the Poisson model will be similar.

[Cruise]

@Rick_SAS

Thanks Rick. That makes sense. 

I'm conducting multiple analyses for the number of cancer types as they're separate datasets. Don't I have to correct on the level of p-value since that I'm conducting that many number of tests for each cancer types? I guess, I have to understand what Bonferroni correction is. 

[Rick_SAS]

Yes, I think that is a good idea. You might want to discuss this with a colleague/mentor/advisor. Either CANCER_TYPE is part of the model or it isn't. If it is part of the model, then it belongs on the CLASS statement and you can test for differences between cancer types. If it isn't part of the model, you can't test for differences.

To use a Bonferroni adjustment: you use alpha/k as the significance level, where k is the number of hypotheses that you are testing. So if you are testing 10 hypotheses, you would only count a test as "significant" when the p-value is less than 0.005.  You can adjust CIs the same way: the significance level used to form the standard errors is adjusted. For example, if you are using a t statistic, you would compute a CI by using t_crit = quantile("T", 0.005, df) instead of t_crit = quantile("T", 0.05, df) .

[Cruise]

Got it. All agree. Thank you!