Adding Relative Risk and Incidence Ratio columns in Modified Poisson Regression

[MGhab]

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;

[Ksharp]

You mean this  Rate Ratio ?

http://support.sas.com/kb/24/188.html

[Ksharp]

I just check the result . It is ODDS not RELRISK .

[MGhab]

Hi. thanks for your response. I don't understand what you mean by this line: "I just check the result . It is ODDS not RELRISK ."

I tried the code and link that you provided above. But adding an lsmeans line gives me the error:
"The model does not have a GLM parameterization. This parameterization is required for the
TEST, LSMEANS, LSMESTIMATE, and SLICE statement. These statements are ignored."

Here is my 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;
run;

[Ksharp]

I mean Rate Ratio Estimates in that URL is ODDS Ratio ,not RelRisk .

Your LOG said you need another design matrix GLM to get Rate Ratio:
/param=ref ref=first;
-->
/param=glm ;

[MGhab]

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!

[SAS_Rob]

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;

[MGhab]

Thanks, Rob.
I used your code, but it does not run to completion. First, I still receive the log warning: "The model does not have a GLM parameterization. This parameterization is required for the TEST, LSMEANS, LSESTIMATE, and SLICE statement. These statements are ignored."
The 20 imputations run through, but the data set work.est_ds has 0 observations and 12 variables so the proc sort command does not work (explanation: data set is empty). And mianalyze does not run because work.est_ds has no observations.

[SAS_Rob]

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?

[MGhab]

Thanks, Rob. This actually brings up another one of my questions. I'm unclear on the meaning of these coefficients in the estimates statements. I'm not sure why there is sometimes "1", or "1 -1", or "1 0". I merely guessed these numbers based on examples I found online. But I'm unable to find any text that explains them. Could you explain that to me as well? Log below.
Thank you!!

20597 Proc GenMod data=irnfcs ;
20598 title3 'MI. Mod Pois Regress. Healthcare provider w covs w/o moderators';
20599 Class id ethnorac (ref='0') b5gp (ref='1') b2 (ref='1') provreg (ref='1') a13status (ref='0')
20599! d3_reasons (ref='0')
20600 a8_7cat (ref='2') e1_3cat (ref='1') d8lgbtcom (ref='0') poverty_2cat (ref='0') a14_5cat
20600! (ref='1')
20601 /param=glm ref=first;
20602 Model f1_2cat = ethnorac b5gp b2 provreg a13status d3_reasons a8_7cat d1_length d8lgbtcom a1
20602! e1_3cat poverty_2cat
20603 a14_5cat / Dist=poisson Link=log;
20604 lsmeans ethnorac b5gp b2 provreg a13status d3_reasons a8_7cat d8lgbtcom /ilink cl;
20605 weight tpcweights;
20606 by _imputation_ ;
20607 Repeated subject=id /printmle type=Ind;
20608 ods output GEEEmpPest= regparms;
20609 Estimate 'ethno' ethnorac 1 /exp;
20610 Estimate 'gender' b5gp 1 -1/exp;
20611 Estimate 'sab' b2 1 /exp;
20612 Estimate 'status' a13status 1 -1/exp;
20613 Estimate 'reasons' d3_reasons 1 /exp;
20614 Estimate 'regions of origin' a8_7cat 1 -1/exp;
20615 Estimate 'length' d1_length 1 /exp;
20616 Estimate 'd8lgbtcom' d8lgbtcom 1 /exp;
20617 ods output estimates=est_ds(where=(index(label,'exp')=1)); *keeps the EXP values and their SE;
20618 run;

NOTE: Class levels for some variables were not printed due to excessive size.
NOTE: Algorithm converged.
NOTE: The scale parameter was held fixed.
NOTE: Algorithm converged.
NOTE: The empirical covariance matrix estimate is used in the LSMEANS statement.
NOTE: The empirical covariance matrix estimate is used in the ESTIMATE statement.
NOTE: The structure of the LSMeans table has changed from an earlier release of SAS.
NOTE: The above message was for the following BY group:
Imputation Number=1
NOTE: Class levels for some variables were not printed due to excessive size.
NOTE: Algorithm converged.
NOTE: The scale parameter was held fixed.
NOTE: Algorithm converged.
NOTE: The empirical covariance matrix estimate is used in the LSMEANS statement.
NOTE: The empirical covariance matrix estimate is used in the ESTIMATE statement.
NOTE: The structure of the LSMeans table has changed from an earlier release of SAS.
NOTE: The above message was for the following BY group:
Imputation Number=2
NOTE: Class levels for some variables were not printed due to excessive size.
NOTE: Algorithm converged.
NOTE: The scale parameter was held fixed.
NOTE: Algorithm converged.
NOTE: The empirical covariance matrix estimate is used in the LSMEANS statement.
NOTE: The empirical covariance matrix estimate is used in the ESTIMATE statement.
NOTE: The structure of the LSMeans table has changed from an earlier release of SAS.
NOTE: The above message was for the following BY group:
Imputation Number=3
/*...etc. for all 20 imputations*/
NOTE: The data set WORK.EST_DS has 0 observations and 12 variables.
NOTE: The data set WORK.REGPARMS has 740 observations and 9 variables.
NOTE: PROCEDURE GENMOD used (Total process time):
real time 47.20 seconds
cpu time 27.35 seconds


20619 proc sort data=est_ds;
20620 by label _imputation_;
20621 run;

NOTE: Input data set is empty.
NOTE: The data set WORK.EST_DS has 0 observations and 12 variables.
NOTE: PROCEDURE SORT used (Total process time):
real time 0.01 seconds
cpu time 0.01 seconds


20622 proc mianalyze data=est_ds;
20623 by label;
20624 modeleffects LBetaEstimate;
20625 stderr stderr;
20626 run;

NOTE: No observations in data set WORK.EST_DS.
NOTE: PROCEDURE MIANALYZE used (Total process time):
real time 0.00 seconds
cpu time 0.00 seconds

[SAS_Rob]

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));

[SAS_Rob]

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.

[MGhab]

Okay. I made that change. The MIANALYZE now produces output for 4 (gender,length,regions,status) of 8 of my predictors - in other words, no output for "ethno", "sab", "reasons", or "d8lgbtcom". These are all of my dichotomous variables, so I wonder if it's an issue with the coefficients, I selected. I changed the code from param=ref to param=glm as per your instructions. I'm not sure of the implications of this.
Based on your explanation of coefficients above, does this mean that a variable with multiple levels will need multiple estimates lines? For example, the b5gp variable has 3 levels:
Estimate 'gender' b5gp 1 -1/exp;*
Estimate 'gender' b5gp 1 /exp;*

Additionally, the mianalyze output provides one single-row LBetaEstimate Variance table and one single-row LBetaEstimate Parameter Estimates table). I'm not sure what this tells me that the original mianalyze output wasn't telling me before.
Thank you!

here is the log for mianlayze.
NOTE: The data set WORK.EST_DS has 80 observations and 12 variables.
NOTE: The data set WORK.REGPARMS has 740 observations and 9 variables.
NOTE: PROCEDURE GENMOD used (Total process time):
real time 44.16 seconds
cpu time 27.37 seconds


21008 proc sort data=est_ds;
21009 by label _imputation_;
21010 run;

NOTE: There were 80 observations read from the data set WORK.EST_DS.
NOTE: The data set WORK.EST_DS has 80 observations and 12 variables.
NOTE: PROCEDURE SORT used (Total process time):
real time 0.01 seconds
cpu time 0.01 seconds


21011 proc mianalyze data=est_ds;
21012 by label;
21013 modeleffects LBetaEstimate;
21014 stderr stderr;
21015 run;

NOTE: PROCEDURE MIANALYZE used (Total process time):
real time 0.07 seconds
cpu time 0.07 seconds