create forest plot using output from multiple models

[jdr1]

Dear SAS Community, 

I've been using SAS for several years and hate to admit that I typically cut and paste model results and graph or create tables using a different program. But I've decided to take my SAS graphing skills farther. 

Here is the problem. I'm looking at temperature sensitive (exposure) in a severe maternal morbidity (outcome) across 5 critical windows of pregnancy: (1) pre-conception, (2) Trimester 1, (3) Trimester 2, (4) Trimester 3, and (5) last month of pregnancy. For each critical period, I'm modeling relative risk for 4 separate exposure outcomes : very cold vs moderate, cold vs moderate, very hot vs moderate, and hot vs moderate across three temperature variables Average temperature, maximum temperature, and minimum temperature. 

For example, results examining average temperature exposure and SMM relative risk across 5 critical periods will produce results from 20 models. How do I create a forest plot for maximum temperature and illustrate Relative Risk (RR) for 4 separate exposure outcomes: very cold vs moderate, cold vs moderate, very hot vs moderate, and hot vs moderate across 5 critical periods of pregnancy? 

Here is example code that produces RR for SMM in response to 4 exposure categories of maximum temperature:

**LMvcold_mod LMvhot_mod LMcold_mod LMhot_mod;

proc genmod data= SMMheat.Tmaxdup;

*by preterm race_ethn;

class RFA_ID  ;

model SMM20 (Event='1')= LMvcold_mod  / dist=poisson type3 link = log dscale;

repeated subject= RFA_ID   /  type=exch modelse;

estimate 'Beta LMvcold_mod' LMvcold_mod 1 / exp;

run;

proc genmod data= SMMheat.Tmaxdup;

*by preterm race_ethn;

class RFA_ID  ;

model SMM20 (Event='1')= LMvhot_mod  / dist=poisson type3 link = log dscale;

repeated subject= RFA_ID   /  type=exch modelse;

estimate 'Beta LMvhot_mod' LMvhot_mod 1 / exp;

run;

proc genmod data= SMMheat.Tmaxdup;

*by preterm race_ethn;

class RFA_ID  ;

model SMM20 (Event='1')= LMcold_mod  / dist=poisson type3 link = log dscale;

repeated subject= RFA_ID   /  type=exch modelse;

estimate 'Beta LMcold_mod' LMcold_mod 1 / exp;

run;

proc genmod data=  SMMheat.Tmaxdup;

*by preterm race_ethn;

class RFA_ID  ;

model SMM20 (Event='1')= LMhot_mod  / dist=poisson type3 link = log dscale;

repeated subject= RFA_ID   /  type=exch modelse;

estimate 'Beta LMhot_mod' LMhot_mod 1 / exp;

run;

[Reeza · Accepted solution]

1. Capture your Estimate output into a data table
   1. Here's some instructions and explanations on how to capture output that is shown.
      https://blogs.sas.com/content/sastraining/2017/03/31/capturing-output-from-any-procedure-with-an-ods...
2. Combine those data sets into one using a data step:

   1. data combined;
      set estimate1 - estimate4;
      run;

3. Create your graph using the data from step 2
   1. Find the plot closest to what you want here and then follow the examples https://blogs.sas.com/content/?s=forest+plot

---

@jdr1 wrote:

Dear SAS Community, 

I've been using SAS for several years and hate to admit that I typically cut and paste model results and graph or create tables using a different program. But I've decided to take my SAS graphing skills farther. 

Here is the problem. I'm looking at temperature sensitive (exposure) in a severe maternal morbidity (outcome) across 5 critical windows of pregnancy: (1) pre-conception, (2) Trimester 1, (3) Trimester 2, (4) Trimester 3, and (5) last month of pregnancy. For each critical period, I'm modeling relative risk for 4 separate exposure outcomes : very cold vs moderate, cold vs moderate, very hot vs moderate, and hot vs moderate across three temperature variables Average temperature, maximum temperature, and minimum temperature. 

For example, results examining average temperature exposure and SMM relative risk across 5 critical periods will produce results from 20 models. How do I create a forest plot for maximum temperature and illustrate Relative Risk (RR) for 4 separate exposure outcomes: very cold vs moderate, cold vs moderate, very hot vs moderate, and hot vs moderate across 5 critical periods of pregnancy? 

Here is example code that produces RR for SMM in response to 4 exposure categories of maximum temperature:

**LMvcold_mod LMvhot_mod LMcold_mod LMhot_mod;

proc genmod data= SMMheat.Tmaxdup;

*by preterm race_ethn;

class RFA_ID  ;

model SMM20 (Event='1')= LMvcold_mod  / dist=poisson type3 link = log dscale;

repeated subject= RFA_ID   /  type=exch modelse;

estimate 'Beta LMvcold_mod' LMvcold_mod 1 / exp;

run;

proc genmod data= SMMheat.Tmaxdup;

*by preterm race_ethn;

class RFA_ID  ;

model SMM20 (Event='1')= LMvhot_mod  / dist=poisson type3 link = log dscale;

repeated subject= RFA_ID   /  type=exch modelse;

estimate 'Beta LMvhot_mod' LMvhot_mod 1 / exp;

run;

proc genmod data= SMMheat.Tmaxdup;

*by preterm race_ethn;

class RFA_ID  ;

model SMM20 (Event='1')= LMcold_mod  / dist=poisson type3 link = log dscale;

repeated subject= RFA_ID   /  type=exch modelse;

estimate 'Beta LMcold_mod' LMcold_mod 1 / exp;

run;

proc genmod data=  SMMheat.Tmaxdup;

*by preterm race_ethn;

class RFA_ID  ;

model SMM20 (Event='1')= LMhot_mod  / dist=poisson type3 link = log dscale;

repeated subject= RFA_ID   /  type=exch modelse;

estimate 'Beta LMhot_mod' LMhot_mod 1 / exp;

run;

---

[ballardw]

Generic approach:

Create data sets for output from each of your Proc Genmod runs. I don't use this proc regularly so I'm not sure which options for an output data set to use. Either the OUTPUT statement or perhaps ODS OUTPUT would be preferred.

Combine the data sets using the INDSNAME option on the SET statement, create variable to hold the group/source information

Then, hopefully that data set has all the information you need to plot the data.

An example here : https://support.sas.com/kb/42/867.html

provides example data and code for one form of forest plot using SGPLOT.

[jdr1]

Many thanks to you both for the quick reply! I think I'm just stuck on how to output the data I need from the models. Here is sample model output and I'd like to be able to export the exponentiated beta for LMvcold_mod using ods output. Does this make sense?

Label	Mean Estimate	Mean
		Confidence Limits
Beta LMvcold_mod	0.8054	0.7236	0.8965

GEE Fit Criteria
QIC	737797.8182
QICu	737750.7264

Analysis Of GEE Parameter Estimates
Empirical Standard Error Estimates
Parameter	Estimate	Standard Error	95% Confidence Limits		Z	Pr > |Z|
Intercept	-4.8008	0.0153	-4.8308	-4.7708	-313.72	<.0001
LMvcold_mod	-0.2164	0.0547	-0.3235	-0.1093	-3.96	<.0001

Analysis Of GEE Parameter Estimates
Model-Based Standard Error Estimates
Parameter	Estimate	Standard Error	95% Confidence Limits		Z	Pr > |Z|
Intercept	-4.8008	0.0043	-4.8092	-4.7924	-1121.0	<.0001
LMvcold_mod	-0.2164	0.0157	-0.2471	-0.1857	-13.82	<.0001
Scale	0.2786	.	.	.	.	.

Note:

The scale parameter was held fixed.

Score Statistics For Type 3 GEE Analysis
Source	DF	Chi-Square	Pr > ChiSq
LMvcold_mod	1	21.69	<.0001

Contrast Estimate Results
Label	Mean Estimate	Mean		L'Beta Estimate	Standard Error	Alpha	L'Beta		Chi-Square	Pr > ChiSq
		Confidence Limits					Confidence Limits
Beta LMvcold_mod	0.8054	0.7236	0.8965	-0.2164	0.0547	0.05	-0.3235	-0.1093	15.67	<.0001
Exp(Beta LMvcold_mod)				0.8054	0.0440	0.05	0.7236	0.8965

[JeffMeyers]

When running any procedure run ODS TRACE ON; prior to it and your log will show the names of the potential datasets you can output with ODS OUTPUT.  For example in the PHREG procedure I often output the ParameterEstimates table that I can see in ODS TRACE.

When you run with ODS TRACE you should see dataset names in the order of the output and there should be labels to help you see which dataset name belongs to which output.  From there you take the dataset name and run ODS OUTPUT datasetname=TheNameYouWant; before your procedure and you will save the output as a dataset.

[Reeza]

Did you read through the first link in my response?

---

@jdr1 wrote:

Many thanks to you both for the quick reply! I think I'm just stuck on how to output the data I need from the models. Here is sample model output and I'd like to be able to export the exponentiated beta for LMvcold_mod using ods output. Does this make sense?

---

[jdr1]

Yes, thanks so much for sharing the two links. I sincerely apologize for the delayed response. It got hectic on my end. 

Here is what I tried:

1) Ran two models with ODS trace on:

ods trace on/listing;

proc genmod data= SMMheat.Tmaxfull_ICE ;

class RFA_ID  COUNTY_OF_RESIDENCE;

model SMM20 (Event='1')= LMvcold_mod  / dist=poisson type3 link = log dscale;

repeated subject= RFA_ID (COUNTY_OF_RESIDENCE)  / type=EXCH modelse;

estimate 'Beta LMvcold_mod' LMvcold_mod 1 / exp;

*repeated subject=zip  / type=EXCH;

*repeated subject=zip(COUNTY_OF_RESIDENCE)  / type=EXCH;

*repeated subject= COUNTY_OF_RESIDENCE  / type=EXCH;

*repeated subject=rfa_id(zip) / type=EXCH;

run;

proc genmod data= SMMheat.Tmaxfull_ICE ;

class RFA_ID  COUNTY_OF_RESIDENCE;

model SMM20 (Event='1')= LMvhot_mod  / dist=poisson type3 link = log dscale;

repeated subject= RFA_ID (COUNTY_OF_RESIDENCE)  / type=EXCH modelse;

estimate 'Beta LMvhot_mod' LMvhot_mod 1 / exp;

*repeated subject=zip  / type=EXCH;

*repeated subject=zip(COUNTY_OF_RESIDENCE)  / type=EXCH;

*repeated subject= COUNTY_OF_RESIDENCE  / type=EXCH;

*repeated subject=rfa_id(zip) / type=EXCH;

run;

ods trace off;

**I'm interesting in saving the Estimates output from several models into one file:

Output Added:

-------------

Name:       Estimates

Label:      ESTIMATE Statement Results

Template:   stat.genmod.estimates

Path:       Genmod.Estimates

2) Next I tried to rename the Estimates dataset and print, but got stuck here: 

ods output Estimates = test;

ods output close;

proc print data= test noobs;

run;

**And I receive this error message: 

WARNING: Output 'Estimates' was not created.  Make sure that the output object name, label, or path

         is spelled correctly.  Also, verify that the appropriate procedure options are used to

         produce the requested output object.  For example, verify that the NOPRINT option is not

         used.

386  ods listing;

387  ods output Estimates = test;

388  ods output close;

389  proc print data= test noobs;

ERROR: File WORK.TEST.DATA does not exist.

390  run;

NOTE: The SAS System stopped processing this step because of errors.

NOTE: PROCEDURE PRINT used (Total process time):

      real time           0.00 seconds

      cpu time

2a) I then tried this code:

ods listing close;

ods output  Estimates = two;

proc print data=two noobs ;

run;

ods listing;

And receive the same error message

[JeffMeyers]

Change your order up a little bit:

ods output Estimates = test;

proc genmod data= SMMheat.Tmaxfull_ICE ;

class RFA_ID  COUNTY_OF_RESIDENCE;

model SMM20 (Event='1')= LMvcold_mod  / dist=poisson type3 link = log dscale;

repeated subject= RFA_ID (COUNTY_OF_RESIDENCE)  / type=EXCH modelse;

estimate 'Beta LMvcold_mod' LMvcold_mod 1 / exp;

*repeated subject=zip  / type=EXCH;

*repeated subject=zip(COUNTY_OF_RESIDENCE)  / type=EXCH;

*repeated subject= COUNTY_OF_RESIDENCE  / type=EXCH;

*repeated subject=rfa_id(zip) / type=EXCH;

run;

proc print data= test noobs;

run;

[jdr1]

Changing the order up worked - thank you!

I'm working through Steps 1 to 3 below. If I'm capturing ODS output for multiple models at once, how do I perform Step 2? Given that each model doesn't have a unique Estimates1 - Estimates 20? Should I be running each model individually like you just showed me and rename each as Test1, Test2...Test20?

Really appreciative of your help with these novice questions? 

1. Capture your Estimate output into a data table
   1. Here's some instructions and explanations on how to capture output that is shown.
      https://blogs.sas.com/content/sastraining/2017/03/31/capturing-output-from-any-procedure-with-an-ods...
2. Combine those data sets into one using a data step:

   1. data combined;
      set estimate1 - estimate4;
      run;

3. Create your graph using the data from step 2
   1. Find the plot closest to what you want here and then follow the examples https://blogs.sas.com/content/?s=forest+plot

[JeffMeyers]

Yes, you'll have to give each model's ESTIMATE dataset a different name:

ODS OUTPUT estimates=test1;

**Model code;

ODS OUTPUT estimates=test2;

**Model code;

**combine them:

data all_models;

    set test1 (in=m1) test2 (in=m2);

    if m1 then model=1;

    else if m2 then model=2;

run;

The IN operators in the set statement creates a temporary logic variable so that observations in each dataset can be identified separately.  Doing this you can assign an identifier (I made the variable MODEL above) to distinguish between the different models.  Then if i have the same covariates in different models I can still sort and keep them separated in their respective model output.

[Reeza]

You can skip that last part of the IN using the INDSNAME option.
If you use a systematic naming convention you'll be able to trace your model results back to your output.

data combined;
set estimate1 - estimate4 indsname= inputSrc;
model = inputSrc;
run;

FYI - I believe a PROC can have many ESTIMATE statements. If you use only one per PROC that's your choice - I assumed that's because they're from different models. But IF the model is the same, you can have more than one estimate statement and pipe them all to one file.

[jdr1]

Both ways work perfectly! Really appreciate your patience with me. So grateful to Reeza and JeffMeyers.

[jdr1]

Sorry, one last question. I ended up pulling two variables for each estimate statement:

Is there a way to remove all Exp(Beta x) variables (e.g., Exp(Beta LMvcold_mod)) at the end from the final table? 

The SAS System

Label	MeanEstimate	MeanLowerCL	MeanUpperCL	model
Beta LMvcold_mod	0.8054	0.7236	0.8965	WORK.ESTIMATE1
Exp(Beta LMvcold_mod)	_	_	_	WORK.ESTIMATE1
Beta LMvhot_mod	1.0221	0.9877	1.0576	WORK.ESTIMATE2
Exp(Beta LMvhot_mod)	_	_	_	WORK.ESTIMATE2
Beta LMcold_mod	0.7898	0.6862	0.9090	WORK.ESTIMATE3
Exp(Beta LMcold_mod)	_	_	_	WORK.ESTIMATE3
Beta LMhot_mod	1.0878	0.9922	1.1925	WORK.ESTIMATE4
Exp(Beta LMhot_mod)	_	_	_	WORK.ESTIMATE4
Beta PCvcold_mod	1.0730	0.8796	1.3089	WORK.ESTIMATE5
Exp(Beta PCvcold_mod)	_	_	_	WORK.ESTIMATE5
Beta PCvhot_mod	1.2473	1.0660	1.4596	WORK.ESTIMATE6
Exp(Beta PCvhot_mod)	_	_	_	WORK.ESTIMATE6
Beta LMcold_mod	1.1581	0.9091	1.4754	WORK.ESTIMATE7
Exp(Beta LMcold_mod)	_	_	_	WORK.ESTIMATE7
Beta LMhot_mod	1.1425	0.9042	1.4437	WORK.ESTIMATE8
Exp(Beta LMhot_mod)	_	_	_	WORK.ESTIMATE8
Beta T1vcold_mod	0.7853	0.6298	0.9792	WORK.ESTIMATE9
Exp(Beta T1vcold_mod)	_	_	_	WORK.ESTIMATE9
Beta T1vhot_mod	1.1497	0.9660	1.3684	WORK.ESTIMATE10
Exp(Beta T1vhot_mod)	_	_	_	WORK.ESTIMATE10
Beta LMcold_mod	0.7740	0.5924	1.0111	WORK.ESTIMATE11
Exp(Beta LMcold_mod)	_	_	_	WORK.ESTIMATE11
Beta T1hot_mod	0.9478	0.7345	1.2230	WORK.ESTIMATE12
Exp(Beta T1hot_mod)	_	_	_	WORK.ESTIMATE12
Beta T2vcold_mod	0.8470	0.6866	1.0449	WORK.ESTIMATE13
Exp(Beta T2vcold_mod)	_	_	_	WORK.ESTIMATE13
Beta T2vhot_mod	1.2549	1.0724	1.4685	WORK.ESTIMATE14
Exp(Beta T2vhot_mod)	_	_	_	WORK.ESTIMATE14
Beta T2cold_mod	0.8007	0.6195	1.0350	WORK.ESTIMATE15
Exp(Beta T2cold_mod)	_	_	_	WORK.ESTIMATE15
Beta T2hot_mod	1.2938	1.0297	1.6256	WORK.ESTIMATE16
Exp(Beta T2hot_mod)	_	_	_	WORK.ESTIMATE16
Beta T3vcold_mod	0.8379	0.6872	1.0217	WORK.ESTIMATE17
Exp(Beta T3vcold_mod)	_	_	_	WORK.ESTIMATE17
Beta T3vhot_mod	1.1442	1.0183	1.2855	WORK.ESTIMATE18
Exp(Beta T3vhot_mod)	_	_	_	WORK.ESTIMATE18
Beta T3cold_mod	0.9266	0.7202	1.1921	WORK.ESTIMATE19
Exp(Beta T3cold_mod)	_	_	_	WORK.ESTIMATE19
Beta T3hot_mod	1.3492	1.0897	1.6705	WORK.ESTIMATE20
Exp(Beta T3hot_mod)	_	_	_	WORK.ESTIMATE20

[Reeza]

Add the following as a data set option on your ODS statement.

where label not like 'Exp%';

or

where not missing(meanEstimate);

[jdr1]

THANK YOU!!!!