Hello,
I am trying to obtain relative risk output for a modified poisson regression after conducting multiple imputation. The resource I found (www(dot)lexjansen(dot)com/wuss/2013/81_Paper(dot)pdf ) provides code that should add "Relative Risk" (RR) and "Incidence Rate Ratio" (IRR) columns to my output, but it isn't working.
Without the RR and IRR columns, I'm not sure what informative output to report. Based on my observations, researchers do not appear to report the estimates column. Here is my code, which was adapted from the code I found on the lex jansen page.
Thank you!
/*Step 1: Imputation Phase*/
proc mi data=tpc2019irn seed=573382369 nimpute=20 out=irnfcs ;
title3 'Imputing Phase irn ppct';
class a14_5cat ethnorac a13status d3_reasons a14 pp_immigrant provreg a8_7cat
b5gp b2 poverty_2cat e1_3cat f1_2cat i2recode d8lgbtcom ;
var a14_5cat a1 ethnorac a13status d3_reasons a14 pp_immigrant provreg a8_7cat
b5gp b2 poverty_2cat e1_3cat f1_2cat i2recode tpcweights Indi_a d1_length
d8lgbtcom;
fcs logistic (provreg = Indi_a d3_reasons /details);
fcs logistic (a8_7cat = Indi_a poverty_2cat d3_reasons ethnorac/ details );
fcs logistic (b5gp = b2 i2recode Indi_a poverty_2cat / details );
fcs discrim (b2 = b5gp Indi_a /classeffects=include);
fcs discrim (poverty_2cat = a1 a13status b5gp b2 Indi_a /classeffects=include);
fcs logistic (e1_3cat = a1 a14 ethnorac a13status d3_reasons a8_7cat f1_2cat/details);
fcs discrim (i2recode = a1 b2 b5gp ethnorac d3_reasons a8_7cat f1_2cat/classeffects=include);
fcs regression (indi_a = a1 b2 b5gp ethnorac a8_7cat poverty_2cat provreg d1_length);
fcs regression (d1_length = a1 a8_7cat poverty_2cat a13status a14);
fcs logistic (d8lgbtcom = a1 a14 b5gp b2 ethnorac a13status d3_reasons /details);
run;
/*Step 2. Mod Poiss Reg phase*/
Proc GenMod data=irnfcs ;
title3 'MI. Mod Pois Regress';
Class id ethnorac (ref='0') b5gp (ref='1') b2 (ref='1') provreg (ref='1') a13status (ref='0') d3_reasons (ref='0')
a8_7cat (ref='2') e1_3cat (ref='1') d8lgbtcom (ref='0') poverty_2cat (ref='0') a14_5cat (ref='1') /param=ref ref=first;
Model f1_2cat = a1 ethnorac b5gp b2 provreg a13status d3_reasons a8_7cat e1_3cat d8lgbtcom poverty_2cat d1_length
a14_5cat / Dist=poisson Link=log;
weight tpcweights;
by _imputation_ ;
Repeated subject=id /printmle type=Ind;
ods output GEEEmpPest= regparms;
Estimate 'ethno' ethnorac 1 /exp;
Estimate 'gender' b5gp 1 -1/exp;
Estimate 'sab' b2 1 /exp;
Estimate 'status' a13status 1 -1/exp;
Estimate 'reasons' d3_reasons 1 0/exp;
Estimate 'regions of origin' a8_7cat 1 -1/exp;
Estimate 'length' d1_length 1 /exp;
Estimate 'd8lgbtcom' d8lgbtcom 1 /exp;
run;
/* Exponentiate parameter estimates */
data regparms_exp;
retain Parm Level1 Estimate expest;
set regparms;
expest = exp(Estimate);
run;
/* Keep the exponentiated Intercept estimate as the baseline incidence
rate */
data _null_;
set regparms_exp;
if Parm = "Intercept";
call symputx('baserate', expest);
run;
/* Calculate group-specific incidence rate from the baseline rate and the
relative risk */
data regparms_IR (drop=expest);
retain Parm Level1 Estimate expest IncidenceRatePer1000 RelativeRisk;
set regparms_exp;
if Parm = "Intercept" then RelativeRisk = .;
else RelativeRisk = expest;
if Parm = "Intercept" then IncidenceRatePer1000 = &baserate * 1000;
else IncidenceRatePer1000 = &baserate * expest * 1000;
run;
/*Proc Mianalyze*/
proc mianalyze parms(classvar=level)=regparms;
class ethnorac e1_3cat b5gp b2 provreg a13status d3_reasons a8_7cat d8lgbtcom poverty_2cat a14_5cat;
modeleffects intercept a1 ethnorac e1_3cat b5gp b2 provreg d1_length a13status d3_reasons a8_7cat
d8lgbtcom poverty_2cat a14_5cat;
run;
A note on my post:
My understanding is that, through proc genmod, I can obtain relative risk measures from the Contrast Estimates Results tables (the Mean Estimate column).
The issue is that I can not seem to obtain overall relative risk estimates, for example, through the amalgamating proc mianalyze stage. Which data from the mianalyze output will give me overall relative risk or incidence rate ratio estimates?
Thank you!
The issue with combining them in MIANALYZE is that there are two sets of estimates when you use the EXP option and you need to capture only the exponentiated estimates and standard errors. What you would need to do is add the following code:
Proc GenMod data=irnfcs ;
title3 'MI. Mod Pois Regress. Healthcare provider w covs w/o moderators';
Class id ethnorac (ref='0') b5gp (ref='1') b2 (ref='1') provreg (ref='1') a13status (ref='0') d3_reasons (ref='0')
a8_7cat (ref='2') e1_3cat (ref='1') d8lgbtcom (ref='0') poverty_2cat (ref='0') a14_5cat (ref='1')
/param=ref ref=first;
Model f1_2cat = ethnorac b5gp b2 provreg a13status d3_reasons a8_7cat d1_length d8lgbtcom a1 e1_3cat poverty_2cat
a14_5cat / Dist=poisson Link=log;
lsmeans ethnorac b5gp b2 provreg a13status d3_reasons a8_7cat d8lgbtcom /ilink cl;
weight tpcweights;
by _imputation_ ;
Repeated subject=id /printmle type=Ind;
ods output GEEEmpPest= regparms;
Estimate 'ethno' ethnorac 1 /exp;
Estimate 'gender' b5gp 1 -1/exp;
Estimate 'sab' b2 1 /exp;
Estimate 'status' a13status 1 -1/exp;
Estimate 'reasons' d3_reasons 1 /exp;
Estimate 'regions of origin' a8_7cat 1 -1/exp;
Estimate 'length' d1_length 1 /exp;
Estimate 'd8lgbtcom' d8lgbtcom 1 /exp;
ods output Estimates=est_ds(where=(index(label,'Exp')=1));*Keeps the EXP values and their standard errors;
run;
proc sort data=est_ds;
by label _imputation_;
run;
proc mianalyze data=est_ds;
by label;
modeleffects LBetaEstimate;
stderr stderr;
run;
Sorry for the confusion. I was only focused on the ESTIMATE statements. You will need to change the CLASS statement to PARAM=GLM if you want to use the LSMEANS statement. Of course if you do that, the coefficients on the ESTIMATE statements will also need to change.
Could you post the full LOG so I can see what you ran?
The WHERE statement is case sensitive--notice in my example I used 'Exp' while in your code you have 'exp' . If you correct that then the MIANALYZE part should run.
ods output estimates=est_ds(where=(index(label,'Exp')=1));
The coefficients reflect specific tests that you want to perform. It is dependent upon the PARAM= option, but assuming you kept it at PARAM=REF then
Estimate 'ethno' ethnorac 1 /exp;*-->would compare the first level of ETHNORAC vs the reference;
Estimate 'gender' b5gp 1 -1/exp;*-->would compare the first level of b5gp vs the second;
You could interpret the others likewise, but it is important that you understand what is going on and how to customize them to what you are actually wanting to report. For this I would suggest starting with the tutorial in this usage note.
Are you ready for the spotlight? We're accepting content ideas for SAS Innovate 2025 to be held May 6-9 in Orlando, FL. The call is open until September 25. Read more here about why you should contribute and what is in it for you!
ANOVA, or Analysis Of Variance, is used to compare the averages or means of two or more populations to better understand how they differ. Watch this tutorial for more.
Find more tutorials on the SAS Users YouTube channel.