130 " Lead Benefits ( Sub 12). sas" Using data from the lead specific data set, changes in blood lead levels are estimated and benefits calculated for reductions in cases of hypertension, CHD, brain infarctions, and mortality for recreational fishermen and their families. 131 /* Program: Lead Benefits ( Sub 12). sas Author: Stuart Smith, Abt Associates, Inc. Date: 07/ 05/ 2000 Updated: 11/ 08/ 2002 This program estimates benefits from reduced consumption of fish tissue containing lead. Reductions in MP& M facility lead loadings result in lower probabilities of lead related illnesses. */

title2 ' Lead Related Benefits Recreational ( Sub 12)';

data lead; set rundat. lead;

/* Combine the exposed populations recreational populations, fresh and saltwater */

exppop = fwsptpop + swsptpop;

/* Calculate the blood lead levels under the baseline and post compliance scenarios */

PbB0m = 4.55; /* Typical PbB conc. in absence of fish cons. for men */ PbB0f = 3.45; /* Typical PbB conc. in absence of fish cons. for women */ bcf = 49; /* Bioconcentration factor for lead L/ KG */ bksf = 0.4; /* Biokinetic slope factor for lead */ ingm = 17.5; /* Average daily fish consumption ( g/ day) males age 15 44 */ ingf = 17.5; /* Average daily fish consumption ( g/ day) females age 15 44 */ ingm2 = 17.5; /* Average daily fish consumption ( g/ day) males over age 45 */ ingf2 = 17.5; /* Average daily fish consumption ( g/ day) females over age 45 */ af = 0.03; /* Absorbtion factor */ ef = 365; /* Exposure Frequency */ at = 365; /* Average timing, or duration of exposure in days */ cf = 0.001; /* Conversion factor kg/ g */

/* Uses ingestion for males ages 15 44 */ PbB_ bm = PbB0m + (( cnc_ base* bcf* ingm* af* bksf* ef* cf)/ at); PbB_ 1m = PbB0m + (( cnc_ trt1* bcf* ingm* af* bksf* ef* cf)/ at); PbB_ 2m = PbB0m + (( cnc_ trt2* bcf* ingm* af* bksf* ef* cf)/ at);

/* Uses ingestion for males ages 45 and older */ PbB_ bme = PbB0m + (( cnc_ base* bcf* ingm2* af* bksf* ef* cf)/ at); PbB_ 1me = PbB0m + (( cnc_ trt1* bcf* ingm2* af* bksf* ef* cf)/ at); PbB_ 2me = PbB0m + (( cnc_ trt2* bcf* ingm2* af* bksf* ef* cf)/ at);

/* Uses ingestion for females ages 15 44 */ PbB_ bf = PbB0f + (( cnc_ base* bcf* ingf* af* bksf* ef* cf)/ at); PbB_ 1f = PbB0f + (( cnc_ trt1* bcf* ingf* af* bksf* ef* cf)/ at); PbB_ 2f = PbB0f + (( cnc_ trt2* bcf* ingf* af* bksf* ef* cf)/ at);

/* Uses ingestion for females ages 45 and older */ PbB_ bfe = PbB0f + (( cnc_ base* bcf* ingf2* af* bksf* ef* cf)/ at); 132 PbB_ 1fe = PbB0f + (( cnc_ trt1* bcf* ingf2* af* bksf* ef* cf)/ at); PbB_ 2fe = PbB0f + (( cnc_ trt2* bcf* ingf2* af* bksf* ef* cf)/ at);

drop fwsptpop fwsubpop swsptpop swsubpop; run;

data men; set lead;

/* Set up some constants */ /* All percentages based on the 1996 Census numbers */

pctmale = 0.319518; /* Percentage of U. S. Pop. who are male, age 20 74 */ pctm4054 = 0.098405; /* Percentage of U. S. Pop. who are male, age 40 54 */ pctm4059 = 0.118945; /* Percentage of U. S. Pop. who are male, age 40 59 */ pctm4574 = 0.12933; /* Percentage of U. S. Pop. who are male, age 45 74 */ pctm5564 = 0.03829; /* Percentage of U. S. Pop. who are male, age 55 64 */ pctm6064 = 0.01775; /* Percentage of U. S. Pop. who are male, age 60 64 */ pctm6574 = 0.031427; /* Percentage of U. S. Pop. who are male, age 65 74 */

bbpm4054 = 79.9; /* Mean diastolic bp from NHANES III, males, age 40 54 */ bbpm4059 = 80.0; /* Mean diastolic bp from NHANES III, males, age 40 59 */ bbpm4574 = 78.8; /* Mean diastolic bp from NHANES III, males, age 45 74 */ bbpm5564 = 79.0; /* Mean diastolic bp from NHANES III, males, age 55 64 */ bbpm6064 = 77.8; /* Mean diastolic bp from NHANES III, males age 60 64 */ bbpm6574 = 76.4; /* Mean diastolic bp from NHANES III, males age 65 74 */

/**************************************************************************/ /* Begin Benefits for Males */ /**************************************************************************/

/* Calculate the change in the probability of hypertension */

pr_ hyp1 = ( 1/( 1 + exp( 2.744  0.793* log( PbB_ bm))))  ( 1/( 1 + exp( 2.744  0.793* log( PbB_ 1m))));

pr_ hyp2 = ( 1/( 1 + exp( 2.744  0.793* log( PbB_ bm))))  ( 1/( 1 + exp( 2.744  0.793* log( PbB_ 2m))));

ben_ hyp1 = pr_ hyp1 * exppop * pctmale * 1141; ben_ hyp2 = pr_ hyp2 * exppop * pctmale * 1141; hyppop = exppop * pctmale; hyp_ red1 = pr_ hyp1 * exppop * pctmale; hyp_ red2 = pr_ hyp2 * exppop * pctmale;

/* Calculate Benefits from Changes in CHD, non fatal cases only, which are 2/ 3 of the cases */

dbpm1 = 1.4 * log( PbB_ bm/ PbB_ 1m); dbpm2 = 1.4 * log( PbB_ bm/ PbB_ 2m);

dbpm1e = 1.4 * log( PbB_ bme/ PbB_ 1me); dbpm2e = 1.4 * log( PbB_ bme/ PbB_ 2me); 133 /* Ages 40 59, change in the 10 year prob of occ. */ d1m4059 = bbpm4059  dbpm1; d2m4059 = bbpm4059  dbpm2;

pr1_ 4059 = ( 1/( 1 + exp( 4.996  0.030365* bbpm4059)))  ( 1/( 1 + exp( 4.996  0.030365* d1m4059)));

pr2_ 4059 = ( 1/( 1 + exp( 4.996  0.030365* bbpm4059)))  ( 1/( 1 + exp( 4.996  0.030365* d2m4059)));

/* Ages 60 64, change in the 2 year prob of occ. */ d1m6064 = bbpm6064  dbpm1e; d2m6064 = bbpm6064  dbpm2e;

pr1_ 6064 = ( 1/( 1 + exp( 5.19676  0.02351* bbpm6064)))  ( 1/( 1 + exp( 5.19676  0.02351* d1m6064)));

pr2_ 6064 = ( 1/( 1 + exp( 5.19676  0.02351* bbpm6064)))  ( 1/( 1 + exp( 5.19676  0.02351* d2m6064)));

/* Ages 65 74, change in the 2 year prob of occ. */ d1m6574 = bbpm6574  dbpm1e; d2m6574 = bbpm6574  dbpm2e;

pr1_ 6574 = ( 1/( 1 + exp( 4.90723  0.02031* bbpm6574)))  ( 1/( 1 + exp( 4.90723  0.02031* d1m6574)));

pr2_ 6574 = ( 1/( 1 + exp( 4.90723  0.02031* bbpm6574)))  ( 1/( 1 + exp( 4.90723  0.02031* d2m6574)));

b14059 = ( pr1_ 4059 * exppop * pctm4059 * 76347 * 0.66)/ 10; b24059 = ( pr2_ 4059 * exppop * pctm4059 * 76347 * 0.66)/ 10; b16064 = ( pr1_ 6064 * exppop * pctm6064 * 76347 * 0.66)/ 2; b26064 = ( pr2_ 6064 * exppop * pctm6064 * 76347 * 0.66)/ 2; b16574 = ( pr1_ 6574 * exppop * pctm6574 * 76347 * 0.66)/ 2; b26574 = ( pr2_ 6574 * exppop * pctm6574 * 76347 * 0.66)/ 2;

/* Number of reduced cases Multiply by the 2/ 3 (. 66) and divide by the years of occurance associated with the probability to get cases for one year. */

r14059 = ( pr1_ 4059 * exppop * pctm4059 * 0.66)/ 10; r24059 = ( pr2_ 4059 * exppop * pctm4059 * 0.66)/ 10; r16064 = ( pr1_ 6064 * exppop * pctm6064 * 0.66)/ 2; r26064 = ( pr2_ 6064 * exppop * pctm6064 * 0.66)/ 2; r16574 = ( pr1_ 6574 * exppop * pctm6574 * 0.66)/ 2; r26574 = ( pr2_ 6574 * exppop * pctm6574 * 0.66)/ 2;

/* Calculate the exposed populations at each age group */ popm4059 = exppop * pctm4059; popm6064 = exppop * pctm6064; 134 popm6574 = exppop * pctm6574;

/* Calculate Benefits from Reduced Cerebrovascular Accidents and Brain Infarctions, affects men aged 45 74 Non fatal only, which is 70% of all cases */

/* Ages 45 74, change in the 2 year prob of occ. */ d1m4574 = bbpm4574  dbpm1e; d2m4574 = bbpm4574  dbpm2e;

pr1_ cba = ( 1/( 1 + exp( 8.58889  0.04066* bbpm4574)))  ( 1/( 1 + exp( 8.58889  0.04066* d1m4574)));

pr2_ cba = ( 1/( 1 + exp( 8.58889  0.04066* bbpm4574)))  ( 1/( 1 + exp( 8.58889  0.04066* d2m4574)));

pr1_ bi = ( 1/( 1 + exp( 9.9516  0.04840* bbpm4574)))  ( 1/( 1 + exp( 9.9516  0.04840* d1m4574)));

pr2_ bi = ( 1/( 1 + exp( 9.9516  0.04840* bbpm4574)))  ( 1/( 1 + exp( 9.9516  0.04840* d2m4574)));

/* Number of reduced cases Multiply by the 70% (. 7) and divide by the years of occurance associated with the probability to get cases for one year. */

b1_ cba = ( pr1_ cba * exppop * pctm4574 * 335135 * 0.7)/ 2; b2_ cba = ( pr2_ cba * exppop * pctm4574 * 335135 * 0.7)/ 2; b1_ bi = ( pr1_ bi * exppop * pctm4574 * 335135 * 0.7)/ 2; b2_ bi = ( pr2_ bi * exppop * pctm4574 * 335135 * 0.7)/ 2;

r1_ cba = ( pr1_ cba * exppop * pctm4574 * 0.7)/ 2; r2_ cba = ( pr2_ cba * exppop * pctm4574 * 0.7)/ 2; r1_ bi = ( pr1_ bi * exppop * pctm4574 * 0.7)/ 2; r2_ bi = ( pr2_ bi * exppop * pctm4574 * 0.7)/ 2;

/* Calculate the exposed population*/ popm4574 = exppop * pctm4574;

/* Calculate Benefits from Premature Mortality */

/* Ages 40 54, change in the 12 year prob of death */ d1m4054 = bbpm4054  dbpm1; d2m4054 = bbpm4054  dbpm2;

m1_ 4054 = ( 1/( 1 + exp( 5.3158  0.03516* bbpm4054)))  ( 1/( 1 + exp( 5.3158  0.03516* d1m4054))); 135 m2_ 4054 = ( 1/( 1 + exp( 5.3158  0.03516* bbpm4054)))  ( 1/( 1 + exp( 5.3158  0.03516* d2m4054)));

/* Ages 55 64, change in the 2 year prob of death */ d1m5564 = bbpm5564  dbpm1e; d2m5564 = bbpm5564  dbpm2e;

m1_ 5564 = ( 1/( 1 + exp( 4.89528  0.01866* bbpm5564)))  ( 1/( 1 + exp( 4.89528  0.01866* d1m5564)));

m2_ 5564 = ( 1/( 1 + exp( 4.89528  0.01866* bbpm5564)))  ( 1/( 1 + exp( 4.89528  0.01866* d2m5564)));

/* Ages 65 74, change in the 2 year prob of death */ d1m6574 = bbpm6574  dbpm1e; d2m6574 = bbpm6574  dbpm2e;

m1_ 6574 = ( 1/( 1 + exp( 3.05723  0.00547* bbpm6574)))  ( 1/( 1 + exp( 3.05723  0.00547* d1m6574)));

m2_ 6574 = ( 1/( 1 + exp( 3.05723  0.00547* bbpm6574)))  ( 1/( 1 + exp( 3.05723  0.00547* d2m6574)));

/* Number of reduced cases and Benefits Divide by the years of occurance associated with the probability to get cases for one year. */

m14054l = ( m1_ 4054 * exppop * pctm4054 * 2500000)/ 12; m24054l = ( m2_ 4054 * exppop * pctm4054 * 2500000)/ 12; m14054h = ( m1_ 4054 * exppop * pctm4054 * 13600000)/ 12; m24054h = ( m2_ 4054 * exppop * pctm4054 * 13600000)/ 12; mid40541 = ( m1_ 4054 * exppop * pctm4054 * 6500000)/ 12; mid40542 = ( m2_ 4054 * exppop * pctm4054 * 6500000)/ 12;

m15564l = ( m1_ 5564 * exppop * pctm5564 * 2500000)/ 2; m25564l = ( m2_ 5564 * exppop * pctm5564 * 2500000)/ 2; m15564h = ( m1_ 5564 * exppop * pctm5564 * 13600000)/ 2; m25564h = ( m2_ 5564 * exppop * pctm5564 * 13600000)/ 2; mid55641 = ( m1_ 5564 * exppop * pctm5564 * 6500000)/ 2; mid55642 = ( m2_ 5564 * exppop * pctm5564 * 6500000)/ 2;

m16574l = ( m1_ 6574 * exppop * pctm6574 * 2500000)/ 2; m26574l = ( m2_ 6574 * exppop * pctm6574 * 2500000)/ 2; m16574h = ( m1_ 6574 * exppop * pctm6574 * 13600000)/ 2; m26574h = ( m2_ 6574 * exppop * pctm6574 * 13600000)/ 2; mid65741 = ( m1_ 6574 * exppop * pctm6574 * 6500000)/ 2; mid65742 = ( m2_ 6574 * exppop * pctm6574 * 6500000)/ 2;

r14054l = ( m1_ 4054 * exppop * pctm4054)/ 12; r24054l = ( m2_ 4054 * exppop * pctm4054)/ 12; r14054h = ( m1_ 4054 * exppop * pctm4054)/ 12; r24054h = ( m2_ 4054 * exppop * pctm4054)/ 12; 136 r15564l = ( m1_ 5564 * exppop * pctm5564)/ 2; r25564l = ( m2_ 5564 * exppop * pctm5564)/ 2; r15564h = ( m1_ 5564 * exppop * pctm5564)/ 2; r25564h = ( m2_ 5564 * exppop * pctm5564)/ 2; r16574l = ( m1_ 6574 * exppop * pctm6574)/ 2; r26574l = ( m2_ 6574 * exppop * pctm6574)/ 2; r16574h = ( m1_ 6574 * exppop * pctm6574)/ 2; r26574h = ( m2_ 6574 * exppop * pctm6574)/ 2;

/* Calculate the exposed populations */ popm4054 = exppop * pctm4054; popm5564 = exppop * pctm5564; run;

data women; set lead;

/* Set up some constants */ /* All percentages based on the 1996 Census numbers */ pctfem = 0.511735; /* Percentage of U. S. Pop. who are female */ pctf4574 = 0.143282; /* Percentage of U. S. Pop. who are female, age 45 74 */

bbpf4574 = 74.8; /* Mean diastolic bp from NHANES III, females age 45 74 */

/**************************************************************************/ /* Begin Benefits for Females */ /**************************************************************************/

/* No dose response function for female hypertension */

/* Calculate Benefits from Changes in CHD, non fatal cases only which are 2/ 3 of the cases */

dbpf1 = ( 0.6 * 1.4) * log( PbB_ bfe/ PbB_ 1fe); dbpf2 = ( 0.6 * 1.4) * log( PbB_ bfe/ PbB_ 2fe);

/* Ages 45 74, change in the 2 year prob of occ. */ d1f4574 = bbpf4574  dbpf1; d2f4574 = bbpf4574  dbpf2;

pr1_ 4574 = ( 1/( 1 + exp( 6.9401  0.03072* bbpf4574)))  ( 1/( 1 + exp( 6.9401  0.03072* d1f4574)));

pr2_ 4574 = ( 1/( 1 + exp( 6.9401  0.03072* bbpf4574)))  ( 1/( 1 + exp( 6.9401  0.03072* d2f4574)));

b1fchd = ( pr1_ 4574 * exppop * pctf4574 * 76347 * 0.66)/ 2; b2fchd = ( pr2_ 4574 * exppop * pctf4574 * 76347 * 0.66)/ 2;

red_ chd1 = ( pr1_ 4574 * exppop * pctf4574 * 0.66)/ 2; red_ chd2 = ( pr2_ 4574 * exppop * pctf4574 * 0.66)/ 2; 137 popf4574 = exppop * pctf4574;

/* Calculate Benefits from Reduced Cerebrovascular Accidents and Brain Infarctions, affects women aged 45 74 Non fatal only which is 70% of all cases */

/* Ages 45 74, change in the 2 year prob of occ. */

pr1f_ bi = ( 1/( 1 + exp( 10.6716  0.0544* bbpf4574)))  ( 1/( 1 + exp( 10.6716  0.0544* d1f4574)));

pr2f_ bi = ( 1/( 1 + exp( 10.6716  0.0544* bbpf4574)))  ( 1/( 1 + exp( 10.6716  0.0544* d2f4574)));

pr1f_ cba = ( 1/( 1 + exp( 9.07737  0.04287* bbpf4574)))  ( 1/( 1 + exp( 9.07737  0.04287* d1f4574)));

pr2f_ cba = ( 1/( 1 + exp( 9.07737  0.04287* bbpf4574)))  ( 1/( 1 + exp( 9.07737  0.04287* d2f4574)));

b1f_ cba = ( pr1f_ cba * exppop * pctf4574 * 251351 * 0.7)/ 2; b2f_ cba = ( pr2f_ cba * exppop * pctf4574 * 251351 * 0.7)/ 2; b1f_ bi = ( pr1f_ bi * exppop * pctf4574 * 251351 * 0.7)/ 2; b2f_ bi = ( pr2f_ bi * exppop * pctf4574 * 251351 * 0.7)/ 2;

r1f_ cba = ( pr1f_ cba * exppop * pctf4574 * 0.7)/ 2; r2f_ cba = ( pr2f_ cba * exppop * pctf4574 * 0.7)/ 2; r1f_ bi = ( pr1f_ bi * exppop * pctf4574 * 0.7)/ 2; r2f_ bi = ( pr2f_ bi * exppop * pctf4574 * 0.7)/ 2;

/* Calculate Benefits from Premature Mortality */

/* Ages 45 74, change in the 2 year prob of death */

m1_ 4574 = ( 1/( 1 + exp( 5.40374  0.01511* bbpf4574)))  ( 1/( 1 + exp( 5.40374  0.01511* d1f4574)));

m2_ 4574 = ( 1/( 1 + exp( 5.40374  0.01511* bbpf4574)))  ( 1/( 1 + exp( 5.40374  0.01511* d2f4574)));

m14574l = ( m1_ 4574 * exppop * pctf4574 * 2500000)/ 2; m24574l = ( m2_ 4574 * exppop * pctf4574 * 2500000)/ 2; m14574h = ( m1_ 4574 * exppop * pctf4574 * 13600000)/ 2; m24574h = ( m2_ 4574 * exppop * pctf4574 * 13600000)/ 2; mid45741 = ( m1_ 4574 * exppop * pctf4574 * 6500000)/ 2; mid45742 = ( m2_ 4574 * exppop * pctf4574 * 6500000)/ 2;

r14574l = ( m1_ 4574 * exppop * pctf4574)/ 2; r24574l = ( m2_ 4574 * exppop * pctf4574)/ 2; r14574h = ( m1_ 4574 * exppop * pctf4574)/ 2; r24574h = ( m2_ 4574 * exppop * pctf4574)/ 2; run; 138 /* Sum benefits and cases and print out results */

proc means sum data= men; var ben_ hyp1 ben_ hyp2 hyppop b14059 b24059 b16064 b26064 b16574 b26574 popm4059 popm6064 popm6574 b1_ cba b2_ cba b1_ bi b2_ bi popm4574 m14054l m14054h m24054l m24054h m15564l m15564h m25564l m25564h m16574l m16574h m26574l m26574h mid40541 mid40542 mid55641 mid55642 mid65741 mid65742 popm4054 popm5564 hyp_ red1 hyp_ red2 r14059 r24059 r16064 r26064 r16574 r26574 r1_ cba r2_ cba r1_ bi r2_ bi r14054l r24054l r15564l r25564l r16574l r26574l; output out = rundat. msptsamp( drop = _ type_ _ freq_) sum = ben_ hyp1 ben_ hyp2 hyppop b14059 b24059 b16064 b26064 b16574 b26574 popm4059 popm6064 popm6574 b1_ cba b2_ cba b1_ bi b2_ bi popm4574 m14054l m14054h m24054l m24054h m15564l m15564h m25564l m25564h m16574l m16574h m26574l m26574h mid40541 mid40542 mid55641 mid55642 mid65741 mid65742 popm4054 popm5564 hyp_ red1 hyp_ red2 r14059 r24059 r16064 r26064 r16574 r26574 r1_ cba r2_ cba r1_ bi r2_ bi r14054l r24054l r15564l r25564l r16574l r26574l; title3 ' Mens Lead Related Benefits ( Sample Basis)'; quit;

proc means sum data= men; weight event_ wt; var ben_ hyp1 ben_ hyp2 hyppop b14059 b24059 b16064 b26064 b16574 b26574 popm4059 popm6064 popm6574 b1_ cba b2_ cba b1_ bi b2_ bi popm4574 m14054l m14054h m24054l m24054h m15564l m15564h m25564l m25564h m16574l m16574h m26574l m26574h mid40541 mid40542 mid55641 mid55642 mid65741 mid65742 popm4054 popm5564 hyp_ red1 hyp_ red2 r14059 r24059 r16064 r26064 r16574 r26574 r1_ cba r2_ cba r1_ bi r2_ bi r14054l r24054l r15564l r25564l r16574l r26574l; output out = rundat. msptwght( drop = _ type_ _ freq_) sum = 139 ben_ hyp1 ben_ hyp2 hyppop b14059 b24059 b16064 b26064 b16574 b26574 popm4059 popm6064 popm6574 b1_ cba b2_ cba b1_ bi b2_ bi popm4574 m14054l m14054h m24054l m24054h m15564l m15564h m25564l m25564h m16574l m16574h m26574l m26574h mid40541 mid40542 mid55641 mid55642 mid65741 mid65742 popm4054 popm5564 hyp_ red1 hyp_ red2 r14059 r24059 r16064 r26064 r16574 r26574 r1_ cba r2_ cba r1_ bi r2_ bi r14054l r24054l r15564l r25564l r16574l r26574l; title3 ' Mens Lead Related Benefits ( National Basis)'; quit;

% macro maleleadtable( tbl_ in=, tbl_ out=); data rundat.& tbl_ out ( keep= gender categ agecat type case1 ben1l ben1m ben1h case2 ben2l ben2m ben2h); set rundat.& tbl_ in;

length gender $ 50 categ $ 50 agecat $ 50 type $ 50;

array _ hyper[ 4] ben_ hyp1 ben_ hyp2 hyp_ red1 hyp_ red2; array _ chd[ 12] b14059 b24059 b16064 b26064 b16574 b26574 r14059 r24059 r16064 r26064 r16574 r26574; array _ cba[ 4] b1_ cba b2_ cba r1_ cba r2_ cba; array _ bi[ 4] b1_ bi b2_ bi r1_ bi r2_ bi; array _ mort[ 24] m14054l mid40541 m14054h m24054l mid40542 m24054h m15564l mid55641 m15564h m25564l mid55642 m25564h m16574l mid65741 m16574h m26574l mid65742 m26574h r14054l r24054l r15564l r25564l r16574l r26574l; array ben[ 6] ben1l ben1m ben1h ben2l ben2m ben2h; array cases[ 2] case1 case2;

do i = 1 to 9; if i = 1 then do; ben[ 1] = _ hyper[ 1]; ben[ 2] = _ hyper[ 1]; ben[ 3] = _ hyper[ 1]; ben[ 4] = _ hyper[ 2]; ben[ 5] = _ hyper[ 2]; ben[ 6] = _ hyper[ 2]; cases[ 1] = _ hyper[ 3]; cases[ 2] = _ hyper[ 4]; categ = " Hypertension"; agecat = " All"; type = " Sport"; gender = " Male"; output; end; if i = 2 then do; ben[ 1] = _ chd[ 1]; ben[ 2] = _ chd[ 1]; ben[ 3] = _ chd[ 1]; ben[ 4] = _ chd[ 2]; ben[ 5] = _ chd[ 2]; ben[ 6] = _ chd[ 2]; cases[ 1] = _ chd[ 7]; cases[ 2] = _ chd[ 8]; 140 categ = " CHD"; agecat = " Age 40 59"; type = " Sport"; gender = " Male"; output; end; if i = 3 then do; ben[ 1] = _ chd[ 3]; ben[ 2] = _ chd[ 3]; ben[ 3] = _ chd[ 3]; ben[ 4] = _ chd[ 4]; ben[ 4] = _ chd[ 4]; ben[ 6] = _ chd[ 4]; cases[ 1] = _ chd[ 9]; cases[ 2] = _ chd[ 10]; categ = " CHD"; agecat = " Age 60 64"; type = " Sport"; gender = " Male"; output; end; if i = 4 then do; ben[ 1] = _ chd[ 5]; ben[ 2] = _ chd[ 5]; ben[ 3] = _ chd[ 5]; ben[ 4] = _ chd[ 6]; ben[ 4] = _ chd[ 6]; ben[ 6] = _ chd[ 6]; cases[ 1] = _ chd[ 11]; cases[ 2] = _ chd[ 12]; categ = " CHD"; agecat = " Age 65 74"; type = " Sport"; gender = " Male"; output; end; if i = 5 then do; ben[ 1] = _ cba[ 1]; ben[ 2] = _ cba[ 1]; ben[ 3] = _ cba[ 1]; ben[ 4] = _ cba[ 2]; ben[ 5] = _ cba[ 2]; ben[ 6] = _ cba[ 2]; cases[ 1] = _ cba[ 3]; cases[ 2] = _ cba[ 4]; categ = " CBA"; agecat = " All"; type = " Sport"; gender = " Male"; output; end; if i = 6 then do; ben[ 1] = _ bi[ 1]; ben[ 2] = _ bi[ 1]; ben[ 3] = _ bi[ 1]; ben[ 4] = _ bi[ 2]; ben[ 5] = _ bi[ 2]; ben[ 6] = _ bi[ 2]; cases[ 1] = _ bi[ 3]; cases[ 2] = _ bi[ 4]; categ = " BI"; agecat = " All"; type = " Sport"; gender = " Male"; output; end; if i = 7 then do; ben[ 1] = _ mort[ 1]; ben[ 2] = _ mort[ 2]; ben[ 3] = _ mort[ 3]; ben[ 4] = _ mort[ 4]; ben[ 5] = _ mort[ 5]; ben[ 6] = _ mort[ 6]; cases[ 1] = _ mort[ 19]; cases[ 2] = _ mort[ 20]; categ = " Mortality"; agecat = " Age 40 54"; type = " Sport"; 141 gender = " Male"; output; end; if i = 8 then do; ben[ 1] = _ mort[ 7]; ben[ 2] = _ mort[ 8]; ben[ 3] = _ mort[ 9]; ben[ 4] = _ mort[ 10]; ben[ 5] = _ mort[ 11]; ben[ 6] = _ mort[ 12]; cases[ 1] = _ mort[ 21]; cases[ 2] = _ mort[ 22]; categ = " Mortality"; agecat = " Age 55 64"; type = " Sport"; gender = " Male"; output; end; if i = 9 then do; ben[ 1] = _ mort[ 13]; ben[ 2] = _ mort[ 14]; ben[ 3] = _ mort[ 15]; ben[ 4] = _ mort[ 16]; ben[ 5] = _ mort[ 17]; ben[ 6] = _ mort[ 18]; cases[ 1] = _ mort[ 23]; cases[ 2] = _ mort[ 24]; categ = " Mortality"; agecat = " Age 65 74"; type = " Sport"; gender = " Male"; output; end; end; run; % mend;

% maleleadtable( tbl_ in= msptsamp, tbl_ out= tblmssmp); % maleleadtable( tbl_ in= msptwght, tbl_ out= tblmswgt);

proc means sum data= women; var b1fchd b2fchd popf4574 b1f_ cba b2f_ cba b1f_ bi b2f_ bi m14574l m14574h m24574l m24574h mid45741 mid45742 red_ chd1 red_ chd2 r1f_ cba r2f_ cba r1f_ bi r2f_ bi r14574l r14574h r24574l r24574h; output out = rundat. fsptsamp( drop = _ type_ _ freq_) sum = b1fchd b2fchd popf4574 b1f_ cba b2f_ cba b1f_ bi b2f_ bi m14574l m14574h m24574l m24574h mid45741 mid45742 red_ chd1 red_ chd2 r1f_ cba r2f_ cba r1f_ bi r2f_ bi r14574l r14574h r24574l r24574h; title3 ' Female Lead Related Benefits ( Sample Basis)'; quit;

proc means sum data= women; weight event_ wt; var b1fchd b2fchd popf4574 b1f_ cba b2f_ cba b1f_ bi b2f_ bi m14574l m14574h m24574l m24574h mid45741 mid45742 red_ chd1 red_ chd2 142 r1f_ cba r2f_ cba r1f_ bi r2f_ bi r14574l r14574h r24574l r24574h; output out = rundat. fsptwght( drop = _ type_ _ freq_) sum = b1fchd b2fchd popf4574 b1f_ cba b2f_ cba b1f_ bi b2f_ bi m14574l m14574h m24574l m24574h mid45741 mid45742 red_ chd1 red_ chd2 r1f_ cba r2f_ cba r1f_ bi r2f_ bi r14574l r14574h r24574l r24574h; title3 ' Female Lead Related Benefits ( National Basis)'; quit;

% macro femaleleadtable( tbl_ in=, tbl_ out=); data rundat.& tbl_ out ( keep= gender categ agecat type case1 ben1l ben1m ben1h case2 ben2l ben2m ben2h); set rundat.& tbl_ in;

length gender $ 50 categ $ 50 agecat $ 50 type $ 50;

array _ chd[ 4] b1fchd b2fchd red_ chd1 red_ chd2; array _ cba[ 4] b1f_ cba b2f_ cba r1f_ cba r2f_ cba; array _ bi[ 4] b1f_ bi b2f_ bi r1f_ bi r2f_ bi; array _ mort[ 8] m14574l mid45741 m14574h m24574l mid45742 m24574h r14574l r24574l; array ben[ 6] ben1l ben1m ben1h ben2l ben2m ben2h; array cases[ 2] case1 case2;

do i = 1 to 4; if i = 1 then do; ben[ 1] = _ chd[ 1]; ben[ 2] = _ chd[ 1]; ben[ 3] = _ chd[ 1]; ben[ 4] = _ chd[ 2]; ben[ 5] = _ chd[ 2]; ben[ 6] = _ chd[ 2]; cases[ 1] = _ chd[ 3]; cases[ 2] = _ chd[ 4]; categ = " CHD"; agecat = " All"; type = " Sport"; gender = " Female"; output; end; if i = 2 then do; ben[ 1] = _ cba[ 1]; ben[ 2] = _ cba[ 1]; ben[ 3] = _ cba[ 1]; ben[ 4] = _ cba[ 2]; ben[ 5] = _ cba[ 2]; ben[ 6] = _ cba[ 2]; cases[ 1] = _ cba[ 3]; cases[ 2] = _ cba[ 4]; categ = " CBA"; agecat = " All"; type = " Sport"; gender = " Female"; output; end; if i = 3 then do; ben[ 1] = _ bi[ 1]; ben[ 2] = _ bi[ 1]; ben[ 3] = _ bi[ 1]; ben[ 4] = _ bi[ 2]; ben[ 5] = _ bi[ 2]; ben[ 6] = _ bi[ 2]; cases[ 1] = _ bi[ 3]; cases[ 2] = _ bi[ 4]; 143 categ = " BI"; agecat = " All"; type = " Sport"; gender = " Female"; output; end; if i = 4 then do; ben[ 1] = _ mort[ 1]; ben[ 2] = _ mort[ 2]; ben[ 3] = _ mort[ 3]; ben[ 4] = _ mort[ 4]; ben[ 5] = _ mort[ 5]; ben[ 6] = _ mort[ 6]; cases[ 1] = _ mort[ 7]; cases[ 2] = _ mort[ 8]; categ = " Mortality"; agecat = " Age 45 74"; type = " Sport"; gender = " Female"; output; end; end; run; % mend;

% femaleleadtable( tbl_ in= fsptsamp, tbl_ out= tblfssmp); % femaleleadtable( tbl_ in= fsptwght, tbl_ out= tblfswgt);

data rundat. tblspt_ s; set rundat. tblmssmp rundat. tblfssmp; run;

data rundat. tblspt_ w; set rundat. tblmswgt rundat. tblfswgt; run;

proc means sum data= men; var event_ wt; title3 ' Number of Affected Reaches Nationally'; quit;

proc sql;

create table malepop as select sum( exppop) as tpops, sum( exppop* pctmale) as pm2074s, sum( exppop* pctm4574) as pm4574s, sum( exppop* event_ wt) as tpopn, sum( exppop* pctmale* event_ wt) as pm2074n, sum( exppop* pctm4574* event_ wt) as pm4574n from men;

create table fempop as select sum( exppop) as tpops, sum( exppop* pctf4574) as pf4574s, sum( exppop* event_ wt) as tpopn, sum( exppop* pctf4574* event_ wt) as pf4574n from women; 144 quit;

proc print data= malepop; var tpops pm2074s pm4574s tpopn pm2074n pm4574n; label tpops = ' Total Exposed Pop ( male + fem) ( Smpl)' pm2074s = ' Exposed Males Age 20 74 ( Smpl)' pm4574s = ' Exposed Males Age 45 74 ( Smpl)' tpopn = ' Total Exposed Pop ( male + fem) ( Natl)' pm2074n = ' Exposed Males Age 20 74 ( Natl)' pm4574n = ' Exposed Males Age 45 74 ( Natl)'; title3 ' Exposed Population  Males'; quit;

proc print data= fempop; var tpops pf4574s tpopn pf4574n; label tpops = ' Total Exposed Pop ( male + fem) ( Smpl)' pf4574s = ' Exposed Females Age 45 74 ( Smpl)' tpopn = ' Total Exposed Pop ( male + fem) ( Natl)' pf4574n = ' Exposed Females Age 45 74 ( Natl)'; title3 ' Exposed Population  Females'; quit;

/* End Lead Benefits ( Sub 12). sas */ 145 Program Output 146 MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 94 Lead Related Benefits Recreational ( Sub 12) Mens Lead Related Benefits ( National Basis) 00: 44 Sunday, February 9, 2003 The MEANS Procedure Variable Sum                          ben_ hyp1 460499.98 ben_ hyp2 553567.91 hyppop 513469.79 b14059 18274.49 b24059 22009.81 b16064 5135.78 b26064 6185.43 b16574 7945.64 b26574 9569.43 popm4059 191146.24 popm6064 28524.49 popm6574 50503.62 b1_ cba 42475.92 b2_ cba 51161.33 b1_ bi 23926.04 b2_ bi 28819.42 popm4574 207835.08 m14054l 761904.57 m14054h 4144760.84 m24054l 917654.37 m24054h 4992039.77 m15564l 415061.53 m15564h 2257934.71 m25564l 499881.04 m25564h 2719352.87 m16574l 202836.08 m16574h 1103428.29 m26574l 244272.35 m26574h 1328841.56 mid40541 1980951.87 mid40542 2385901.36 mid55641 1079159.97 mid55642 1299690.71 mid65741 527373.81 mid65742 635108.10 popm4054 158138.18 popm5564 61532.55 hyp_ red1 403.5933249 hyp_ red2 485.1603105 r14059 0.2393610 r24059 0.2882866 r16064 0.0672689 r26064 0.0810173 r16574 0.1040728 147 r26574 0.1253413 r1_ cba 0.1267427 r2_ cba 0.1526589 r1_ bi 0.0713922 r2_ bi 0.0859935 r14054l 0.3047618 r24054l 0.3670617 r15564l 0.1660246 r25564l 0.1999524 r16574l 0.0811344 r26574l 0.0977089                          MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 96 Lead Related Benefits Recreational ( Sub 12) Female Lead Related Benefits ( National Basis) 00: 44 Sunday, February 9, 2003 The MEANS Procedure Variable Sum                          b1fchd 12490.59 b2fchd 15058.89 popf4574 230256.13 b1f_ cba 18054.88 b2f_ cba 21768.25 b1f_ bi 11050.22 b2f_ bi 13323.46 m14574l 437169.27 m14574h 2378200.82 m24574l 527031.61 m24574h 2867051.98 mid45741 1136640.10 mid45742 1370282.19 red_ chd1 0.1636029 red_ chd2 0.1972428 r1f_ cba 0.0718314 r2f_ cba 0.0866050 r1f_ bi 0.0439633 r2f_ bi 0.0530074 r14574l 0.1748677 r14574h 0.1748677 r24574l 0.2108126 r24574h 0.2108126                          MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 97 Lead Related Benefits Recreational ( Sub 12) Number of Affected Reaches Nationally 00: 44 Sunday, February 9, 2003 The MEANS Procedure Analysis Variable : EVENT_ WT Event Weight 148 Sum              198.9830931              MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 98 Lead Related Benefits Recreational ( Sub 12) Exposed Population  Males 00: 44 Sunday, February 9, 2003 Obs tpops pm2074s pm4574s tpopn pm2074n pm4574n 1 269981 86263.79 34916.64 1607013.67 513469.79 207835.08 MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 99 Lead Related Benefits Recreational ( Sub 12) Exposed Population  Females 00: 44 Sunday, February 9, 2003 Obs tpops pf4574s tpopn pf4574n 1 269981 38683.42 1607013.67 230256.13 149 " Subsistence Lead Benefits ( Sub 17). sas" Using data from the lead specific data set, changes in blood lead levels are estimated and benefits calculated for reductions in cases of hypertension, CHD, brain infarctions, and mortality for subsistence fishermen and their families. 150 /* Program: Subsistence Lead Benefits ( Sub 17). sas Author: Stuart Smith, Abt Associates, Inc. Date: 07/ 11/ 2000 Updated: 11/ 08/ 2002 This program estimates benefits from reduced consumption of fish tissue containing lead. Reductions in MP& M facility lead loadings result in lower probabilties of lead related illnesses. */

% let tbl_ path = % sysfunc( pathname( rundat));

title2 ' Subsistence Lead Related Benefits ( Sub 17)';

data lead; set rundat. lead;

/* Combine the exposed populations for subsistence; fresh and saltwater. */

exppop = fwsubpop + swsubpop;

/* Calculate the blood lead levels under the baseline and post compliance scenarios */

PbB0m = 5.40; /* Typical PbB conc. in absence of fish cons. for men */ PbB0f = 4.12; /* Typical PbB conc. in absence of fish cons. for women */ bcf = 49; /* Bioconcentration factor for lead L/ KG */ bksf = 0.4; /* Biokinetic slope factor for lead */ ingm = 142.4; /* Average daily fish consumption ( g/ day) males age 15 44 */ ingf = 142.4; /* Average daily fish consumption ( g/ day) females age 15 44 */ ingm2 = 142.4; /* Average daily fish consumption ( g/ day) males over age 45 */ ingf2 = 142.4; /* Average daily fish consumption ( g/ day) females over age 45 */ af = 0.03; /* Absorbtion factor */ ef = 365; /* Exposure Frequency */ at = 365; /* Average timing, or duration of exposure in days */ cf = 0.001; /* Conversion factor kg/ g */

/* Uses ingestion for males ages 15 44 */ PbB_ bm = PbB0m + (( cnc_ base* bcf* ingm* af* bksf* ef* cf)/ at); PbB_ 1m = PbB0m + (( cnc_ trt1* bcf* ingm* af* bksf* ef* cf)/ at); PbB_ 2m = PbB0m + (( cnc_ trt2* bcf* ingm* af* bksf* ef* cf)/ at);

/* Uses ingestion for males over age 45 */ PbB_ bme = PbB0m + (( cnc_ base* bcf* ingm2* af* bksf* ef* cf)/ at); PbB_ 1me = PbB0m + (( cnc_ trt1* bcf* ingm2* af* bksf* ef* cf)/ at); PbB_ 2me = PbB0m + (( cnc_ trt2* bcf* ingm2* af* bksf* ef* cf)/ at);

/* Uses ingestion for females ages 15 44 */ PbB_ bf = PbB0f + (( cnc_ base* bcf* ingf* af* bksf* ef* cf)/ at); PbB_ 1f = PbB0f + (( cnc_ trt1* bcf* ingf* af* bksf* ef* cf)/ at); PbB_ 2f = PbB0f + (( cnc_ trt2* bcf* ingf* af* bksf* ef* cf)/ at); 151 /* Uses ingestion for females over age 45 */ PbB_ bfe = PbB0f + (( cnc_ base* bcf* ingf2* af* bksf* ef* cf)/ at); PbB_ 1fe = PbB0f + (( cnc_ trt1* bcf* ingf2* af* bksf* ef* cf)/ at); PbB_ 2fe = PbB0f + (( cnc_ trt2* bcf* ingf2* af* bksf* ef* cf)/ at);

drop fwsptpop fwsubpop swsptpop swsubpop; run;

/* Print reaches in the lead dataset */ proc print; var reachnmb PbB_ bm PbB_ 1m PbB_ 2m cnc_ base cnc_ trt1 cnc_ trt2 exppop event_ wt; quit;

data men; set lead;

/* Set up some constants */ /* All percentages based on the 1996 Census numbers */

pctmale = 0.319518; /* Percentage of U. S. Pop. who are male, age 20 74 */ pctm4054 = 0.098405; /* Percentage of U. S. Pop. who are male, age 40 54 */ pctm4059 = 0.118945; /* Percentage of U. S. Pop. who are male, age 40 59 */ pctm4574 = 0.12933; /* Percentage of U. S. Pop. who are male, age 45 74 */ pctm5564 = 0.03829; /* Percentage of U. S. Pop. who are male, age 55 64 */ pctm6064 = 0.01775; /* Percentage of U. S. Pop. who are male, age 60 64 */ pctm6574 = 0.031427; /* Percentage of U. S. Pop. who are male, age 65 74 */

bbpm4054 = 81.9; /* Mean diastolic bp from NHANES III, males, age 40 54 */ bbpm4059 = 81.8; /* Mean diastolic bp from NHANES III, males, age 40 59 */ bbpm4574 = 81.1; /* Mean diastolic bp from NHANES III, males, age 45 74 */ bbpm5564 = 80.6; /* Mean diastolic bp from NHANES III, males, age 55 64 */ bbpm6064 = 79.5; /* Mean diastolic bp from NHANES III, males age 60 64 */ bbpm6574 = 79.5; /* Mean diastolic bp from NHANES III, males age 65 74 */

/**************************************************************************/ /* Begin Benefits for Males */ /**************************************************************************/

/* Calculate the change in the probability of hypertension */

pr_ hyp1 = ( 1/( 1 + exp( 2.744  0.793* log( PbB_ bm))))  ( 1/( 1 + exp( 2.744  0.793* log( PbB_ 1m))));

pr_ hyp2 = ( 1/( 1 + exp( 2.744  0.793* log( PbB_ bm))))  ( 1/( 1 + exp( 2.744  0.793* log( PbB_ 2m))));

ben_ hyp1 = pr_ hyp1 * exppop * pctmale * 1141; ben_ hyp2 = pr_ hyp2 * exppop * pctmale * 1141; hyppop = exppop * pctmale; hyp_ red1 = pr_ hyp1 * exppop * pctmale; hyp_ red2 = pr_ hyp2 * exppop * pctmale; 152 /* Calculate Benefits from Changes in CHD, non fatal cases only which are 2/ 3 of the cases */

dbpm1 = 1.4 * log( PbB_ bm/ PbB_ 1m); dbpm2 = 1.4 * log( PbB_ bm/ PbB_ 2m);

dbpm1e = 1.4 * log( PbB_ bme/ PbB_ 1me); dbpm2e = 1.4 * log( PbB_ bme/ PbB_ 2me);

/* Ages 40 59, change in the 10 year prob of occ. */ d1m4059 = bbpm4059  dbpm1; d2m4059 = bbpm4059  dbpm2;

pr1_ 4059 = ( 1/( 1 + exp( 4.996  0.030365* bbpm4059)))  ( 1/( 1 + exp( 4.996  0.030365* d1m4059)));

pr2_ 4059 = ( 1/( 1 + exp( 4.996  0.030365* bbpm4059)))  ( 1/( 1 + exp( 4.996  0.030365* d2m4059)));

/* Ages 60 64, change in the 2 year prob of occ. */ d1m6064 = bbpm6064  dbpm1e; d2m6064 = bbpm6064  dbpm2e;

pr1_ 6064 = ( 1/( 1 + exp( 5.19676  0.02351* bbpm6064)))  ( 1/( 1 + exp( 5.19676  0.02351* d1m6064)));

pr2_ 6064 = ( 1/( 1 + exp( 5.19676  0.02351* bbpm6064)))  ( 1/( 1 + exp( 5.19676  0.02351* d2m6064)));

/* Ages 65 74, change in the 2 year prob of occ. */ d1m6574 = bbpm6574  dbpm1e; d2m6574 = bbpm6574  dbpm2e;

pr1_ 6574 = ( 1/( 1 + exp( 4.90723  0.02031* bbpm6574)))  ( 1/( 1 + exp( 4.90723  0.02031* d1m6574)));

pr2_ 6574 = ( 1/( 1 + exp( 4.90723  0.02031* bbpm6574)))  ( 1/( 1 + exp( 4.90723  0.02031* d2m6574)));

b14059 = ( pr1_ 4059 * exppop * pctm4059 * 76347 * 0.66)/ 10; b24059 = ( pr2_ 4059 * exppop * pctm4059 * 76347 * 0.66)/ 10; b16064 = ( pr1_ 6064 * exppop * pctm6064 * 76347 * 0.66)/ 2; b26064 = ( pr2_ 6064 * exppop * pctm6064 * 76347 * 0.66)/ 2; b16574 = ( pr1_ 6574 * exppop * pctm6574 * 76347 * 0.66)/ 2; b26574 = ( pr2_ 6574 * exppop * pctm6574 * 76347 * 0.66)/ 2;

/* Number of reduced cases */ r14059 = ( pr1_ 4059 * exppop * pctm4059 * 0.66)/ 10; r24059 = ( pr2_ 4059 * exppop * pctm4059 * 0.66)/ 10; r16064 = ( pr1_ 6064 * exppop * pctm6064 * 0.66)/ 2; r26064 = ( pr2_ 6064 * exppop * pctm6064 * 0.66)/ 2; r16574 = ( pr1_ 6574 * exppop * pctm6574 * 0.66)/ 2; r26574 = ( pr2_ 6574 * exppop * pctm6574 * 0.66)/ 2; 153 popm4059 = exppop * pctm4059; popm6064 = exppop * pctm6064; popm6574 = exppop * pctm6574;

/* Calculate Benefits from Reduced Cerebrovascular Accidents and Brain Infarctions, affects men aged 45 74 Non fatal only which is 70% of all cases */

/* Ages 45 74, change in the 2 year prob of occ. */ d1m4574 = bbpm4574  dbpm1e; d2m4574 = bbpm4574  dbpm2e;

pr1_ cba = ( 1/( 1 + exp( 8.58889  0.04066* bbpm4574)))  ( 1/( 1 + exp( 8.58889  0.04066* d1m4574)));

pr2_ cba = ( 1/( 1 + exp( 8.58889  0.04066* bbpm4574)))  ( 1/( 1 + exp( 8.58889  0.04066* d2m4574)));

pr1_ bi = ( 1/( 1 + exp( 9.9516  0.04840* bbpm4574)))  ( 1/( 1 + exp( 9.9516  0.04840* d1m4574)));

pr2_ bi = ( 1/( 1 + exp( 9.9516  0.04840* bbpm4574)))  ( 1/( 1 + exp( 9.9516  0.04840* d2m4574)));

b1_ cba = ( pr1_ cba * exppop * pctm4574 * 335135 * 0.7)/ 2; b2_ cba = ( pr2_ cba * exppop * pctm4574 * 335135 * 0.7)/ 2; b1_ bi = ( pr1_ bi * exppop * pctm4574 * 335135 * 0.7)/ 2; b2_ bi = ( pr2_ bi * exppop * pctm4574 * 335135 * 0.7)/ 2;

r1_ cba = ( pr1_ cba * exppop * pctm4574 * 0.7)/ 2; r2_ cba = ( pr2_ cba * exppop * pctm4574 * 0.7)/ 2; r1_ bi = ( pr1_ bi * exppop * pctm4574 * 0.7)/ 2; r2_ bi = ( pr2_ bi * exppop * pctm4574 * 0.7)/ 2;

popm4574 = exppop * pctm4574;

/* Calculate Benefits from Premature Mortality */

/* Ages 40 54, change in the 12 year prob of death */ d1m4054 = bbpm4054  dbpm1; d2m4054 = bbpm4054  dbpm2;

m1_ 4054 = ( 1/( 1 + exp( 5.3158  0.03516* bbpm4054)))  ( 1/( 1 + exp( 5.3158  0.03516* d1m4054)));

m2_ 4054 = ( 1/( 1 + exp( 5.3158  0.03516* bbpm4054)))  ( 1/( 1 + exp( 5.3158  0.03516* d2m4054)));

/* Ages 55 64, change in the 2 year prob of death */ d1m5564 = bbpm5564  dbpm1e; d2m5564 = bbpm5564  dbpm2e; 154 m1_ 5564 = ( 1/( 1 + exp( 4.89528  0.01866* bbpm5564)))  ( 1/( 1 + exp( 4.89528  0.01866* d1m5564)));

m2_ 5564 = ( 1/( 1 + exp( 4.89528  0.01866* bbpm5564)))  ( 1/( 1 + exp( 4.89528  0.01866* d2m5564)));

/* Ages 65 74, change in the 2 year prob of death */ d1m6574 = bbpm6574  dbpm1e; d2m6574 = bbpm6574  dbpm2e;

m1_ 6574 = ( 1/( 1 + exp( 3.05723  0.00547* bbpm6574)))  ( 1/( 1 + exp( 3.05723  0.00547* d1m6574)));

m2_ 6574 = ( 1/( 1 + exp( 3.05723  0.00547* bbpm6574)))  ( 1/( 1 + exp( 3.05723  0.00547* d2m6574)));

m14054l = ( m1_ 4054 * exppop * pctm4054 * 2500000)/ 12; m24054l = ( m2_ 4054 * exppop * pctm4054 * 2500000)/ 12; m14054h = ( m1_ 4054 * exppop * pctm4054 * 13600000)/ 12; m24054h = ( m2_ 4054 * exppop * pctm4054 * 13600000)/ 12; mid40541 = ( m1_ 4054 * exppop * pctm4054 * 6500000)/ 12; mid40542 = ( m2_ 4054 * exppop * pctm4054 * 6500000)/ 12;

m15564l = ( m1_ 5564 * exppop * pctm5564 * 2500000)/ 2; m25564l = ( m2_ 5564 * exppop * pctm5564 * 2500000)/ 2; m15564h = ( m1_ 5564 * exppop * pctm5564 * 13600000)/ 2; m25564h = ( m2_ 5564 * exppop * pctm5564 * 13600000)/ 2; mid55641 = ( m1_ 5564 * exppop * pctm5564 * 6500000)/ 2; mid55642 = ( m2_ 5564 * exppop * pctm5564 * 6500000)/ 2;

m16574l = ( m1_ 6574 * exppop * pctm6574 * 2500000)/ 2; m26574l = ( m2_ 6574 * exppop * pctm6574 * 2500000)/ 2; m16574h = ( m1_ 6574 * exppop * pctm6574 * 13600000)/ 2; m26574h = ( m2_ 6574 * exppop * pctm6574 * 13600000)/ 2; mid65741 = ( m1_ 6574 * exppop * pctm6574 * 6500000)/ 2; mid65742 = ( m2_ 6574 * exppop * pctm6574 * 6500000)/ 2;

r14054l = ( m1_ 4054 * exppop * pctm4054)/ 12; r24054l = ( m2_ 4054 * exppop * pctm4054)/ 12; r14054h = ( m1_ 4054 * exppop * pctm4054)/ 12; r24054h = ( m2_ 4054 * exppop * pctm4054)/ 12; r15564l = ( m1_ 5564 * exppop * pctm5564)/ 2; r25564l = ( m2_ 5564 * exppop * pctm5564)/ 2; r15564h = ( m1_ 5564 * exppop * pctm5564)/ 2; r25564h = ( m2_ 5564 * exppop * pctm5564)/ 2; r16574l = ( m1_ 6574 * exppop * pctm6574)/ 2; r26574l = ( m2_ 6574 * exppop * pctm6574)/ 2; r16574h = ( m1_ 6574 * exppop * pctm6574)/ 2; r26574h = ( m2_ 6574 * exppop * pctm6574)/ 2;

popm4054 = exppop * pctm4054; popm5564 = exppop * pctm5564; 155 run;

data women; set lead;

/* Set up some constants */ /* All percentages based on the 1996 Census numbers */ pctfem = 0.511735; /* Percentage of U. S. Pop. who are female */ pctf4574 = 0.143282; /* Percentage of U. S. Pop. who are female, age 45 74 */ bbpf4574 = 76.5; /* Mean diastolic bp from NHANES III, females age 45 74 */

/**************************************************************************/ /* Begin Benefits for Females */ /**************************************************************************/

/* No dose response function for female hypertension */

/* Calculate Benefits from Changes in CHD, non fatal cases only which are 2/ 3 of the cases */

dbpf1 = ( 0.6 * 1.4) * log( PbB_ bfe/ PbB_ 1fe); dbpf2 = ( 0.6 * 1.4) * log( PbB_ bfe/ PbB_ 2fe);

/* Ages 45 74, change in the 2 year prob of occ. */ d1f4574 = bbpf4574  dbpf1; d2f4574 = bbpf4574  dbpf2;

pr1_ 4574 = ( 1/( 1 + exp( 6.9401  0.03072* bbpf4574)))  ( 1/( 1 + exp( 6.9401  0.03072* d1f4574)));

pr2_ 4574 = ( 1/( 1 + exp( 6.9401  0.03072* bbpf4574)))  ( 1/( 1 + exp( 6.9401  0.03072* d2f4574)));

b1fchd = ( pr1_ 4574 * exppop * pctf4574 * 76347 * 0.66)/ 2; b2fchd = ( pr2_ 4574 * exppop * pctf4574 * 76347 * 0.66)/ 2;

red_ chd1 = ( pr1_ 4574 * exppop * pctf4574 * 0.66)/ 2; red_ chd2 = ( pr2_ 4574 * exppop * pctf4574 * 0.66)/ 2;

popf4574 = exppop * pctf4574;

/* Calculate Benefits from Reduced Cerebrovascular Accidents and Brain Infarctions, affects women aged 45 74 Non fatal only which is 70% of all cases */

/* Ages 45 74, change in the 2 year prob of occ. */

pr1f_ bi = ( 1/( 1 + exp( 10.6716  0.0544* bbpf4574)))  ( 1/( 1 + exp( 10.6716  0.0544* d1f4574)));

pr2f_ bi = ( 1/( 1 + exp( 10.6716  0.0544* bbpf4574)))  156 ( 1/( 1 + exp( 10.6716  0.0544* d2f4574)));

pr1f_ cba = ( 1/( 1 + exp( 9.07737  0.04287* bbpf4574)))  ( 1/( 1 + exp( 9.07737  0.04287* d1f4574)));

pr2f_ cba = ( 1/( 1 + exp( 9.07737  0.04287* bbpf4574)))  ( 1/( 1 + exp( 9.07737  0.04287* d2f4574)));

b1f_ cba = ( pr1f_ cba * exppop * pctf4574 * 251351 * 0.7)/ 2; b2f_ cba = ( pr2f_ cba * exppop * pctf4574 * 251351 * 0.7)/ 2; b1f_ bi = ( pr1f_ bi * exppop * pctf4574 * 251351 * 0.7)/ 2; b2f_ bi = ( pr2f_ bi * exppop * pctf4574 * 251351 * 0.7)/ 2;

r1f_ cba = ( pr1f_ cba * exppop * pctf4574 * 0.7)/ 2; r2f_ cba = ( pr2f_ cba * exppop * pctf4574 * 0.7)/ 2; r1f_ bi = ( pr1f_ bi * exppop * pctf4574 * 0.7)/ 2; r2f_ bi = ( pr2f_ bi * exppop * pctf4574 * 0.7)/ 2;

/* Calculate Benefits from Premature Mortality */

/* Ages 45 74, change in the 2 year prob of death */

m1_ 4574 = ( 1/( 1 + exp( 5.40374  0.01511* bbpf4574)))  ( 1/( 1 + exp( 5.40374  0.01511* d1f4574)));

m2_ 4574 = ( 1/( 1 + exp( 5.40374  0.01511* bbpf4574)))  ( 1/( 1 + exp( 5.40374  0.01511* d2f4574)));

m14574l = ( m1_ 4574 * exppop * pctf4574 * 2500000)/ 2; m24574l = ( m2_ 4574 * exppop * pctf4574 * 2500000)/ 2; m14574h = ( m1_ 4574 * exppop * pctf4574 * 13600000)/ 2; m24574h = ( m2_ 4574 * exppop * pctf4574 * 13600000)/ 2; mid45741 = ( m1_ 4574 * exppop * pctf4574 * 6500000)/ 2; mid45742 = ( m2_ 4574 * exppop * pctf4574 * 6500000)/ 2;

r14574l = ( m1_ 4574 * exppop * pctf4574)/ 2; r24574l = ( m2_ 4574 * exppop * pctf4574)/ 2; r14574h = ( m1_ 4574 * exppop * pctf4574)/ 2; r24574h = ( m2_ 4574 * exppop * pctf4574)/ 2; run;

/* Sum up benefits and reduced cases and print results */

proc means sum data= men; var ben_ hyp1 ben_ hyp2 hyppop b14059 b24059 b16064 b26064 b16574 b26574 popm4059 popm6064 popm6574 b1_ cba b2_ cba b1_ bi b2_ bi popm4574 m14054l m14054h m24054l m24054h m15564l m15564h m25564l m25564h m16574l m16574h m26574l m26574h mid40541 mid40542 mid55641 mid55642 mid65741 mid65742 popm4054 popm5564 hyp_ red1 hyp_ red2 157 r14059 r24059 r16064 r26064 r16574 r26574 r1_ cba r2_ cba r1_ bi r2_ bi r14054l r24054l r15564l r25564l r16574l r26574l; output out = rundat. msubsamp( drop = _ type_ _ freq_) sum = ben_ hyp1 ben_ hyp2 hyppop b14059 b24059 b16064 b26064 b16574 b26574 popm4059 popm6064 popm6574 b1_ cba b2_ cba b1_ bi b2_ bi popm4574 m14054l m14054h m24054l m24054h m15564l m15564h m25564l m25564h m16574l m16574h m26574l m26574h mid40541 mid40542 mid55641 mid55642 mid65741 mid65742 popm4054 popm5564 hyp_ red1 hyp_ red2 r14059 r24059 r16064 r26064 r16574 r26574 r1_ cba r2_ cba r1_ bi r2_ bi r14054l r24054l r15564l r25564l r16574l r26574l; title3 ' Mens Lead Related Benefits ( Sample Basis)'; quit;

proc means sum data= men; weight event_ wt; var ben_ hyp1 ben_ hyp2 hyppop b14059 b24059 b16064 b26064 b16574 b26574 popm4059 popm6064 popm6574 b1_ cba b2_ cba b1_ bi b2_ bi popm4574 m14054l m14054h m24054l m24054h m15564l m15564h m25564l m25564h m16574l m16574h m26574l m26574h mid40541 mid40542 mid55641 mid55642 mid65741 mid65742 popm4054 popm5564 hyp_ red1 hyp_ red2 r14059 r24059 r16064 r26064 r16574 r26574 r1_ cba r2_ cba r1_ bi r2_ bi r14054l r24054l r15564l r25564l r16574l r26574l; output out = rundat. msubwght( drop = _ type_ _ freq_) sum = ben_ hyp1 ben_ hyp2 hyppop b14059 b24059 b16064 b26064 b16574 b26574 popm4059 popm6064 popm6574 b1_ cba b2_ cba b1_ bi b2_ bi popm4574 m14054l m14054h m24054l m24054h m15564l m15564h m25564l m25564h m16574l m16574h m26574l m26574h mid40541 mid40542 mid55641 mid55642 mid65741 mid65742 popm4054 popm5564 hyp_ red1 hyp_ red2 r14059 r24059 r16064 r26064 r16574 r26574 r1_ cba r2_ cba r1_ bi r2_ bi r14054l r24054l 158 r15564l r25564l r16574l r26574l; title3 ' Mens Lead Related Benefits ( National Basis)'; quit;

% macro maleleadtable( tbl_ in=, tbl_ out=); data rundat.& tbl_ out ( keep= gender categ agecat type case1 ben1l ben1m ben1h case2 ben2l ben2m ben2h); set rundat.& tbl_ in;

length gender $ 50 categ $ 50 agecat $ 50 type $ 50;

array _ hyper[ 4] ben_ hyp1 ben_ hyp2 hyp_ red1 hyp_ red2; array _ chd[ 12] b14059 b24059 b16064 b26064 b16574 b26574 r14059 r24059 r16064 r26064 r16574 r26574; array _ cba[ 4] b1_ cba b2_ cba r1_ cba r2_ cba; array _ bi[ 4] b1_ bi b2_ bi r1_ bi r2_ bi; array _ mort[ 24] m14054l mid40541 m14054h m24054l mid40542 m24054h m15564l mid55641 m15564h m25564l mid55642 m25564h m16574l mid65741 m16574h m26574l mid65742 m26574h r14054l r24054l r15564l r25564l r16574l r26574l; array ben[ 6] ben1l ben1m ben1h ben2l ben2m ben2h; array cases[ 2] case1 case2;

do i = 1 to 9; if i = 1 then do; ben[ 1] = _ hyper[ 1]; ben[ 2] = _ hyper[ 1]; ben[ 3] = _ hyper[ 1]; ben[ 4] = _ hyper[ 2]; ben[ 5] = _ hyper[ 2]; ben[ 6] = _ hyper[ 2]; cases[ 1] = _ hyper[ 3]; cases[ 2] = _ hyper[ 4]; categ = " Hypertension"; agecat = " All"; type = " Sub"; gender = " Male"; output; end; if i = 2 then do; ben[ 1] = _ chd[ 1]; ben[ 2] = _ chd[ 1]; ben[ 3] = _ chd[ 1]; ben[ 4] = _ chd[ 2]; ben[ 5] = _ chd[ 2]; ben[ 6] = _ chd[ 2]; cases[ 1] = _ chd[ 7]; cases[ 2] = _ chd[ 8]; categ = " CHD"; agecat = " Age 40 59"; type = " Sub"; gender = " Male"; output; end; if i = 3 then do; ben[ 1] = _ chd[ 3]; ben[ 2] = _ chd[ 3]; ben[ 3] = _ chd[ 3]; ben[ 4] = _ chd[ 4]; ben[ 4] = _ chd[ 4]; ben[ 6] = _ chd[ 4]; cases[ 1] = _ chd[ 9]; cases[ 2] = _ chd[ 10]; categ = " CHD"; agecat = " Age 60 64"; 159 type = " Sub"; gender = " Male"; output; end; if i = 4 then do; ben[ 1] = _ chd[ 5]; ben[ 2] = _ chd[ 5]; ben[ 3] = _ chd[ 5]; ben[ 4] = _ chd[ 6]; ben[ 4] = _ chd[ 6]; ben[ 6] = _ chd[ 6]; cases[ 1] = _ chd[ 11]; cases[ 2] = _ chd[ 12]; categ = " CHD"; agecat = " Age 65 74"; type = " Sub"; gender = " Male"; output; end; if i = 5 then do; ben[ 1] = _ cba[ 1]; ben[ 2] = _ cba[ 1]; ben[ 3] = _ cba[ 1]; ben[ 4] = _ cba[ 2]; ben[ 5] = _ cba[ 2]; ben[ 6] = _ cba[ 2]; cases[ 1] = _ cba[ 3]; cases[ 2] = _ cba[ 4]; categ = " CBA"; agecat = " All"; type = " Sub"; gender = " Male"; output; end; if i = 6 then do; ben[ 1] = _ bi[ 1]; ben[ 2] = _ bi[ 1]; ben[ 3] = _ bi[ 1]; ben[ 4] = _ bi[ 2]; ben[ 5] = _ bi[ 2]; ben[ 6] = _ bi[ 2]; cases[ 1] = _ bi[ 3]; cases[ 2] = _ bi[ 4]; categ = " BI"; agecat = " All"; type = " Sub"; gender = " Male"; output; end; if i = 7 then do; ben[ 1] = _ mort[ 1]; ben[ 2] = _ mort[ 2]; ben[ 3] = _ mort[ 3]; ben[ 4] = _ mort[ 4]; ben[ 5] = _ mort[ 5]; ben[ 6] = _ mort[ 6]; cases[ 1] = _ mort[ 19]; cases[ 2] = _ mort[ 20]; categ = " Mortality"; agecat = " Age 40 54"; type = " Sub"; gender = " Male"; output; end; if i = 8 then do; ben[ 1] = _ mort[ 7]; ben[ 2] = _ mort[ 8]; ben[ 3] = _ mort[ 9]; ben[ 4] = _ mort[ 10]; ben[ 5] = _ mort[ 11]; ben[ 6] = _ mort[ 12]; cases[ 1] = _ mort[ 21]; cases[ 2] = _ mort[ 22]; categ = " Mortality"; agecat = " Age 55 64"; type = " Sub"; gender = " Male"; output; 160 end; if i = 9 then do; ben[ 1] = _ mort[ 13]; ben[ 2] = _ mort[ 14]; ben[ 3] = _ mort[ 15]; ben[ 4] = _ mort[ 16]; ben[ 5] = _ mort[ 17]; ben[ 6] = _ mort[ 18]; cases[ 1] = _ mort[ 23]; cases[ 2] = _ mort[ 24]; categ = " Mortality"; agecat = " Age 65 74"; type = " Sub"; gender = " Male"; output; end; end; run; % mend;

% maleleadtable( tbl_ in= msubsamp, tbl_ out= tblmbsmp); % maleleadtable( tbl_ in= msubwght, tbl_ out= tblmbwgt);

proc means sum data= women; var b1fchd b2fchd popf4574 b1f_ cba b2f_ cba b1f_ bi b2f_ bi m14574l m14574h m24574l m24574h mid45741 mid45742 red_ chd1 red_ chd2 r1f_ cba r2f_ cba r1f_ bi r2f_ bi r14574l r14574h r24574l r24574h; output out = rundat. fsubsamp( drop = _ type_ _ freq_) sum = b1fchd b2fchd popf4574 b1f_ cba b2f_ cba b1f_ bi b2f_ bi m14574l m14574h m24574l m24574h mid45741 mid45742 red_ chd1 red_ chd2 r1f_ cba r2f_ cba r1f_ bi r2f_ bi r14574l r14574h r24574l r24574h; title3 ' Female Lead Related Benefits ( Sample Basis)'; quit;

proc means sum data= women; weight event_ wt; var b1fchd b2fchd popf4574 b1f_ cba b2f_ cba b1f_ bi b2f_ bi m14574l m14574h m24574l m24574h mid45741 mid45742 red_ chd1 red_ chd2 r1f_ cba r2f_ cba r1f_ bi r2f_ bi r14574l r14574h r24574l r24574h; output out = rundat. fsubwght( drop = _ type_ _ freq_) sum = b1fchd b2fchd popf4574 b1f_ cba b2f_ cba b1f_ bi b2f_ bi m14574l m14574h m24574l m24574h mid45741 mid45742 red_ chd1 red_ chd2 r1f_ cba r2f_ cba r1f_ bi r2f_ bi r14574l r14574h r24574l r24574h; title3 ' Female Lead Related Benefits ( National Basis)'; quit; 161 % macro femaleleadtable( tbl_ in=, tbl_ out=); data rundat.& tbl_ out ( keep= gender categ agecat type case1 ben1l ben1m ben1h case2 ben2l ben2m ben2h); set rundat.& tbl_ in;

length gender $ 50 categ $ 50 agecat $ 50 type $ 50;

array _ chd[ 4] b1fchd b2fchd red_ chd1 red_ chd2; array _ cba[ 4] b1f_ cba b2f_ cba r1f_ cba r2f_ cba; array _ bi[ 4] b1f_ bi b2f_ bi r1f_ bi r2f_ bi; array _ mort[ 8] m14574l mid45741 m14574h m24574l mid45742 m24574h r14574l r24574l;

array ben[ 6] ben1l ben1m ben1h ben2l ben2m ben2h; array cases[ 2] case1 case2;

do i = 1 to 4; if i = 1 then do; ben[ 1] = _ chd[ 1]; ben[ 2] = _ chd[ 1]; ben[ 3] = _ chd[ 1]; ben[ 4] = _ chd[ 2]; ben[ 5] = _ chd[ 2]; ben[ 6] = _ chd[ 2]; cases[ 1] = _ chd[ 3]; cases[ 2] = _ chd[ 4]; categ = " CHD"; agecat = " All"; type = " Sub"; gender = " Female"; output; end; if i = 2 then do; ben[ 1] = _ cba[ 1]; ben[ 2] = _ cba[ 1]; ben[ 3] = _ cba[ 1]; ben[ 4] = _ cba[ 2]; ben[ 5] = _ cba[ 2]; ben[ 6] = _ cba[ 2]; cases[ 1] = _ cba[ 3]; cases[ 2] = _ cba[ 4]; categ = " CBA"; agecat = " All"; type = " Sub"; gender = " Female"; output; end; if i = 3 then do; ben[ 1] = _ bi[ 1]; ben[ 2] = _ bi[ 1]; ben[ 3] = _ bi[ 1]; ben[ 4] = _ bi[ 2]; ben[ 5] = _ bi[ 2]; ben[ 6] = _ bi[ 2]; cases[ 1] = _ bi[ 3]; cases[ 2] = _ bi[ 4]; categ = " BI"; agecat = " All"; type = " Sub"; gender = " Female"; output; end; if i = 4 then do; ben[ 1] = _ mort[ 1]; ben[ 2] = _ mort[ 2]; ben[ 3] = _ mort[ 3]; ben[ 4] = _ mort[ 4]; ben[ 5] = _ mort[ 5]; ben[ 6] = _ mort[ 6]; cases[ 1] = _ mort[ 7]; cases[ 2] = _ mort[ 8]; categ = " Mortality"; agecat = " Age 45 74"; type = " Sub"; 162 gender = " Female"; output; end; end; run; % mend;

% femaleleadtable( tbl_ in= fsubsamp, tbl_ out= tblfbsmp); % femaleleadtable( tbl_ in= fsubwght, tbl_ out= tblfbwgt);

data rundat. tblsub_ s; set rundat. tblmbsmp rundat. tblfbsmp; run;

data rundat. tblsub_ w; set rundat. tblmbwgt rundat. tblfbwgt; run;

proc means sum data= men; var event_ wt; title3 ' Number of Affected Reaches Nationally'; quit;

proc sql;

create table malepop as select sum( exppop) as tpops, sum( exppop* pctmale) as pm2074s, sum( exppop* pctm4574) as pm4574s, sum( exppop* event_ wt) as tpopn, sum( exppop* pctmale* event_ wt) as pm2074n, sum( exppop* pctm4574* event_ wt) as pm4574n from men;

create table fempop as select sum( exppop) as tpops, sum( exppop* pctf4574) as pf4574s, sum( exppop* event_ wt) as tpopn, sum( exppop* pctf4574* event_ wt) as pf4574n from women; quit;

proc print data= malepop; var tpops pm2074s pm4574s tpopn pm2074n pm4574n; label tpops = ' Total Exposed Pop ( male + fem) ( Smpl)' pm2074s = ' Exposed Males Age 20 74 ( Smpl)' pm4574s = ' Exposed Males Age 45 74 ( Smpl)' tpopn = ' Total Exposed Pop ( male + fem) ( Natl)' pm2074n = ' Exposed Males Age 20 74 ( Natl)' pm4574n = ' Exposed Males Age 45 74 ( Natl)'; title3 ' Exposed Population  Males'; quit; 163 proc print data= fempop; var tpops pf4574s tpopn pf4574n; label tpops = ' Total Exposed Pop ( male + fem) ( Smpl)' pf4574s = ' Exposed Females Age 45 74 ( Smpl)' tpopn = ' Total Exposed Pop ( male + fem) ( Natl)' pf4574n = ' Exposed Females Age 45 74 ( Natl)'; title3 ' Exposed Population  Females'; quit;

data rundat. leadtbls; if _ n_ = 1 then do; do i = 1 to 13; set rundat. tblspt_ s point = i; output; set rundat. tblsub_ s point = i; output; end; end; run;

data rundat. leadtblw; if _ n_ = 1 then do; do i = 1 to 13; set rundat. tblspt_ w point = i; output; set rundat. tblsub_ w point = i; output; end; end; run;

proc summary data= rundat. leadtbls nway; class gender categ; var case1 ben1l ben1m ben1h case2 ben2l ben2m ben2h; output out = rundat. leadsums ( drop= _ type_ _ freq_) sum = case1 ben1l ben1m ben1h case2 ben2l ben2m ben2h; run;

proc summary data= rundat. leadtblw nway; class gender categ; var case1 ben1l ben1m ben1h case2 ben2l ben2m ben2h; output out = rundat. leadsumw ( drop= _ type_ _ freq_) sum = case1 ben1l ben1m ben1h case2 ben2l ben2m ben2h; run;

proc summary data= rundat. leadtbls; var ben1l ben1m ben1h ben2l ben2m ben2h; output out = rundat. ldalls ( drop= _ type_ _ freq_) sum = ben1l ben1m ben1h ben2l ben2m ben2h; run;

proc summary data= rundat. leadtblw; 164 var ben1l ben1m ben1h ben2l ben2m ben2h; output out = rundat. ldallw ( drop= _ type_ _ freq_) sum = ben1l ben1m ben1h ben2l ben2m ben2h; run;

data exlleadw ( keep= titlecol col1 col10 label=" Adult Lead Benefits  Weighted"); set rundat. leadsumw rundat. ldallw; length titlecol $ 50 col1  col10 $ 128;

if _ n_ = 1 then do; titlecol = " Table: Adult Lead Benefits, Weighted Estimates"; output; titlecol = " Run on & sysdate"; output; titlecol = " "; output; col3 = " Option 1"; col7 = " Option 2"; output; col3 = " Reduced Cases"; col4 = " Low"; col5 = " Mid"; col6 = " High"; col7 = " Reduced Cases"; col8 = " Low"; col9 = " Mid"; col10 = " High"; output; end;

if _ n_ = 1 or _ n_ = 6 then do; col1 = gender; end;

if _ n_ = 10 then do; col1 = " Total"; end;

col2 = trim( categ); col3 = trim( case1); col4 = trim( ben1l); col5 = trim( ben1m); col6 = trim( ben1h); col7 = trim( case2); col8 = trim( ben2l); col9 = trim( ben2m); col10 = trim( ben2h); output; run;

proc export data = exlleadw outfile = " C:\ Output Files\ Adult Lead. xls" dbms= excel2000 replace; run;

/* End Subsistence Lead Benefits ( Sub 17). sas */ 165 Program Output 166 MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 102 Subsistence Lead Related Benefits ( Sub 17) Mens Lead Related Benefits ( National Basis) 00: 44 Sunday, February 9, 2003 The MEANS Procedure Variable Sum                          ben_ hyp1 161626.44 ben_ hyp2 196026.78 hyppop 25671.82 b14059 5938.97 b24059 7254.27 b16064 1654.13 b26064 2020.32 b16574 2615.16 b26574 3193.94 popm4059 9556.69 popm6064 1426.13 popm6574 2525.02 b1_ cba 14430.95 b2_ cba 17631.45 b1_ bi 8270.95 b2_ bi 10106.69 popm4574 10391.08 m14054l 250675.02 m14054h 1363672.09 m24054l 306212.15 m24054h 1665794.12 m15564l 132506.87 m15564h 720837.40 m25564l 161829.01 m25564h 880349.79 m16574l 63989.95 m16574h 348105.35 m26574l 78132.38 m26574h 425040.17 mid40541 651755.04 mid40542 796151.60 mid55641 344517.87 mid55642 420755.41 mid65741 166373.88 mid65742 203144.20 popm4054 7906.40 popm5564 3076.43 hyp_ red1 141.6533205 hyp_ red2 171.8026148 r14059 0.0777891 r24059 0.0950171 r16064 0.0216659 r26064 0.0264624 r16574 0.0342536 167 r26574 0.0418346 r1_ cba 0.0430601 r2_ cba 0.0526100 r1_ bi 0.0246794 r2_ bi 0.0301571 r14054l 0.1002700 r24054l 0.1224849 r15564l 0.0530027 r25564l 0.0647316 r16574l 0.0255960 r26574l 0.0312530                          MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 104 Subsistence Lead Related Benefits ( Sub 17) Female Lead Related Benefits ( National Basis) 00: 44 Sunday, February 9, 2003 The MEANS Procedure Variable Sum                          b1fchd 3981.51 b2fchd 4879.32 popf4574 11512.06 b1f_ cba 5873.34 b2f_ cba 7198.84 b1f_ bi 3663.28 b2f_ bi 4490.65 m14574l 135901.62 m14574h 739304.80 m24574l 166515.13 m24574h 905842.32 mid45741 353344.21 mid45742 432939.34 red_ chd1 0.0521502 red_ chd2 0.0639098 r1f_ cba 0.0233671 r2f_ cba 0.0286406 r1f_ bi 0.0145744 r2f_ bi 0.0178660 r14574l 0.0543606 r14574h 0.0543606 r24574l 0.0666061 r24574h 0.0666061                          MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 105 Subsistence Lead Related Benefits ( Sub 17) Number of Affected Reaches Nationally 00: 44 Sunday, February 9, 2003 The MEANS Procedure Analysis Variable : EVENT_ WT Event Weight 168 Sum              198.9830931              MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 106 Subsistence Lead Related Benefits ( Sub 17) Exposed Population  Males 00: 44 Sunday, February 9, 2003 Obs tpops pm2074s pm4574s tpopn pm2074n pm4574n 1 13500 4313.49 1745.96 80345.47 25671.82 10391.08 MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 107 Subsistence Lead Related Benefits ( Sub 17) Exposed Population  Females 00: 44 Sunday, February 9, 2003 Obs tpops pf4574s tpopn pf4574n 1 13500 1934.31 80345.47 11512.06 169 Lead Benefits for Children 170 " Cancer fish ( Sub 4). sas" From the data set of pollutant concentrations for each RF1 reach, observations for in stream lead concentrations are outputed into a separate data set. See Cancer from Fish Consumption for program code and output. 171 " NeoNatal Mortality Benefits ( Sub 15). sas" Using data from the lead specific data set, changes in blood lead levels for women of childbearing age are estimated. Reduced risk of neonatal mortality is calculated and monetized. 172 /* Program: NeoNatal Mortality Benefits ( Sub 15). sas Author: Stuart Smith, Abt Associates, Inc. Date: 07/ 11/ 2000 Updated: 11/ 08/ 2002 This program estimates benefits from reduced consumption of fish tissue containing lead. Reductions in MP& M facility lead loadings result in lower probabilties of lead related illnesses. */

title2 ' NeoNatal Mortality Lead Benefits ( Sub 15)';

data lead; set rundat. lead;

/* Combine the exposed recreational and subsistence populations; fresh and salt water */

sptpop = fwsptpop + swsptpop; subpop = fwsubpop + swsubpop;

/* Calculate the blood lead levels under the baseline and post compliance scenarios */

PbB0f = 2.05; /* Typical PbB conc. in absence of fish cons for women Aged 15 45 */ bcf = 49; /* Bioconcentration factor for lead L/ KG */ bksf = 0.4; /* Biokinetic slope factor for lead */ ing = 17.5; /* Average daily fish consumption ( g/ day) recreational */ ing2 = 142.4; /* Average daily fish consumption ( g/ day) subsistence */ af = 0.03; /* Absorbtion factor */ ef = 365; /* Exposure Frequency */ at = 365; /* Average timing, or duration of exposure in days */ cf = 0.001; /* Conversion factor kg/ g */

/* Recreational */ PbB_ bfr = ( PbB0f + ( cnc_ base* bcf* ing* af* bksf* ef* cf))/ at; PbB_ 1fr = ( PbB0f + ( cnc_ trt1* bcf* ing* af* bksf* ef* cf))/ at; PbB_ 2fr = ( PbB0f + ( cnc_ trt2* bcf* ing* af* bksf* ef* cf))/ at;

/* Subsistence */ PbB_ bfs = ( PbB0f + ( cnc_ base* bcf* ing2* af* bksf* ef* cf))/ at; PbB_ 1fs = ( PbB0f + ( cnc_ trt1* bcf* ing2* af* bksf* ef* cf))/ at; PbB_ 2fs = ( PbB0f + ( cnc_ trt2* bcf* ing2* af* bksf* ef* cf))/ at;

/* Changes in blood lead levels */ diff1r = PbB_ bfr  PbB_ 1fr; diff2r = PbB_ bfr  PbB_ 2fr;

diff1s = PbB_ bfs  PbB_ 1fs; diff2s = PbB_ bfs  PbB_ 2fs;

drop fwsptpop fwsubpop swsptpop swsubpop; 173 run;

proc print; var reachnmb PbB_ bfr PbB_ 1fr PbB_ 2fr PbB_ bfs PbB_ 1fs PbB_ 2fs diff1r diff2r diff1s diff2s cnc_ base cnc_ trt1 cnc_ trt2 sptpop subpop event_ wt; quit;

data neo; set lead;

pctfem = 0.22531; /* Percentage of Women Aged 15 45 ( Childbearing Age) */

expspt = sptpop* pctfem; expsub = subpop* pctfem;

red1r = diff1r * sptpop * pctfem * 0.0001; red2r = diff2r * sptpop * pctfem * 0.0001;

red1s = diff1s * subpop * pctfem * 0.0001; red2s = diff2s * subpop * pctfem * 0.0001;

ben1_ lr = red1r * 2500000; ben1_ hr = red1r * 13600000; ben2_ lr = red2r * 2500000; ben2_ hr = red2r * 13600000; mid1r = red1r * 6500000; mid2r = red2r * 6500000;

ben1_ ls = red1s * 2500000; ben1_ hs = red1s * 13600000; ben2_ ls = red2s * 2500000; ben2_ hs = red2s * 13600000; mid1s = red1s * 6500000; mid2s = red2s * 6500000;

tot1_ l = ben1_ lr + ben1_ ls; tot1_ h = ben1_ hr + ben1_ hs; tot2_ l = ben2_ lr + ben2_ ls; tot2_ h = ben2_ hr + ben2_ hs; mid1 = mid1r + mid1s; mid2 = mid2r + mid2s; run;

proc means sum data= neo; var red1r red2r red1s red2s ben1_ lr ben1_ ls ben1_ hr ben1_ hs ben2_ lr ben2_ ls ben2_ hr ben2_ hs tot1_ l tot1_ h tot2_ l tot2_ h mid1r mid1s mid2r mid2s mid1 mid2 expspt expsub; label red1r = ' Number of Reductions rec ( Opt 1)' red2r = ' Number of Reductions rec ( Opt 2)' red1s = ' Number of Reductions sub ( Opt 1)' red2s = ' Number of Reductions sub ( Opt 2)' 174 ben1_ lr = ' Sport Low ( Opt1)' ben1_ ls = ' Sub Low ( Opt1)' ben1_ hr = ' Sport High ( Opt1)' ben1_ hs = ' Sub High ( Opt1)' ben2_ lr = ' Sport Low ( Opt2)' ben2_ ls = ' Sub Low ( Opt2)' ben2_ hr = ' Sport High ( Opt2)' ben2_ hs = ' Sub High ( Opt2)' tot1_ l = ' Total Benefits  Low ( Opt 1)' tot1_ h = ' Total Benefits  High ( Opt 1)' tot2_ l = ' Total Benefits  Low ( Opt 2)' tot2_ h = ' Total Benefits  High ( Opt 2)' mid1 = ' Mid Benefits ( Opt1)' mid2 = ' Mid Benefits ( Opt2)' expspt = ' Exposed Sport Pop ( 15 45yrs)' expsub = ' Exposed Sub Pop ( 15 45yrs)'; title3 ' Sample Basis'; quit;

proc means sum data= neo; weight event_ wt; var red1r red2r red1s red2s ben1_ lr ben1_ ls ben1_ hr ben1_ hs ben2_ lr ben2_ ls ben2_ hr ben2_ hs tot1_ l tot1_ h tot2_ l tot2_ h mid1r mid1s mid2r mid2s mid1 mid2 expspt expsub; label red1r = ' Number of Reductions rec ( Opt 1)' red2r = ' Number of Reductions rec ( Opt 2)' red1s = ' Number of Reductions sub ( Opt 1)' red2s = ' Number of Reductions sub ( Opt 2)' ben1_ lr = ' Sport Low ( Opt1)' ben1_ ls = ' Sub Low ( Opt1)' ben1_ hr = ' Sport High ( Opt1)' ben1_ hs = ' Sub High ( Opt1)' ben2_ lr = ' Sport Low ( Opt2)' ben2_ ls = ' Sub Low ( Opt2)' ben2_ hr = ' Sport High ( Opt2)' ben2_ hs = ' Sub High ( Opt2)' tot1_ l = ' Total Benefits  Low ( Opt 1)' tot1_ h = ' Total Benefits  High ( Opt 1)' tot2_ l = ' Total Benefits  Low ( Opt 2)' tot2_ h = ' Total Benefits  High ( Opt 2)' mid1 = ' Mid Benefits ( Opt1)' mid2 = ' Mid Benefits ( Opt2)' expspt = ' Exposed Sport Pop ( 15 45yrs)' expsub = ' Exposed Sub Pop ( 15 45yrs)'; output out = rundat. chld_ neo ( drop= _ type_ _ freq_) sum = red1r red2r red1s red2s ben1_ lr ben1_ ls ben1_ hr ben1_ hs ben2_ lr ben2_ ls ben2_ hr ben2_ hs tot1_ l tot1_ h tot2_ l tot2_ h mid1r mid1s mid2r mid2s mid1 mid2 expspt expsub; title3 ' National Basis'; quit;

data rundat. neo_ tbl ( keep= gender categ type exppop case1 case2 benlow1 benlow2 benmid1 benmid2 benhigh1 benhigh2); 175 set rundat. chld_ neo;

length gender $ 50 categ $ 50 type $ 50;

gender = " Children"; categ = " Neo Natal Mortality"; type = " Sport"; exppop = expspt; case1 = red1r; case2 = red2r; benlow1 = ben1_ lr; benlow2 = ben2_ lr; benmid1 = mid1r; benmid2 = mid2r; benhigh1 = ben1_ hr; benhigh2 = ben2_ hr; output;

gender = " Children"; categ = " Neo Natal Mortality"; type = " Sub"; exppop = expsub; case1 = red1s; case2 = red2s; benlow1 = ben1_ ls; benlow2 = ben2_ ls; benmid1 = mid1s; benmid2 = mid2s; benhigh1 = ben1_ hs; benhigh2 = ben2_ hs; output; run;

/* End NeoNatal Mortality Benefits ( Sub 15). sas */ 176 Program Output 177 MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 110 NeoNatal Mortality Lead Benefits ( Sub 15) National Basis 00: 44 Sunday, February 9, 2003 The MEANS Procedure Variable Label Sum                                                              red1r Number of Reductions rec ( Opt 1) 1.1394593 red2r Number of Reductions rec ( Opt 2) 1.3672214 red1s Number of Reductions sub ( Opt 1) 0.4632762 red2s Number of Reductions sub ( Opt 2) 0.5559741 ben1_ lr Sport Low ( Opt1) 2848648.33 ben1_ ls Sub Low ( Opt1) 1158190.54 ben1_ hr Sport High ( Opt1) 15496646.90 ben1_ hs Sub High ( Opt1) 6300556.52 ben2_ lr Sport Low ( Opt2) 3418053.52 ben2_ ls Sub Low ( Opt2) 1389935.18 ben2_ hr Sport High ( Opt2) 18594211.14 ben2_ hs Sub High ( Opt2) 7561247.36 tot1_ l Total Benefits  Low ( Opt 1) 4006838.87 tot1_ h Total Benefits  High ( Opt 1) 21797203.43 tot2_ l Total Benefits  Low ( Opt 2) 4807988.70 tot2_ h Total Benefits  High ( Opt 2) 26155458.50 mid1r 7406485.65 mid1s 3011295.40 mid2r 8886939.15 mid2s 3613831.46 mid1 Mid Benefits ( Opt1) 10417781.05 mid2 Mid Benefits ( Opt2) 12500770.61 expspt Exposed Sport Pop ( 15 45yrs) 362076.25 expsub Exposed Sub Pop ( 15 45yrs) 18102.64                                                              178 Integrated Exposure Uptake Biokinetic Model ( IEUBK) The model used a range of lead concentrations in fish tissue to estimate the corresponding blood lead level associated with a given lead concentration in the fish tissue. 179 " Gen4. txt" Output data set from the IEUBK model run. See DCN 37746 for data set. 180 " Create Lookup Table. sas" Takes the output file from the IEUBK model run and generates another SAS program which when run creates a format used to transform lead concentrations in fish tissue into blood lead levels. 181 /* Program: Create Lookup Table. sas Author: Stuart Smith, Abt Associates Inc. Date: 05/ 18/ 2000 Revised: 03/ 28/ 2002 This program reads in the output file from the IEUBK model and generates a SAS data set with concentrations and predicted blood lead levels. A second file is produced which contains row numbers mapped to concentrations values for lookup by the benefits program. */

libname mpm2 ' H:\ ENVR\ Triton05 MP& M\ Benefits\ MP& M Phase 2 Analysis\ Data\ SAS Sets for Analysis'; /* Static Data Sets */ filename ubk ' H:\ ENVR\ IEUBK\ Gen4. txt'; /* IEUBK output file */

data readin; infile ubk firstobs= 11; input id 1 8 conc 56 61 bloodPb 72 76; run;

data mpm2. ubklook4; set readin; run;

/* Generate the format sas file */ data frmt; set readin; file ' H:\ ENVR\ IEUBK\ cnc_ fmt3. sas'; put conc ' = ' id; run;

/* End Create Lookup Table. sas */ 182 " Ublook4. sd2" Lookup table data set. See DCN 37669 for data set. 183 " Cnc_ fmt3. sas" Creates a library format for transforming lead concentrations in fish tissue into blood lead levels. See Appendix B: Programs and Estimation Results for the Proposed Option National Benefits Analysis, DCN 20479 for program code. 184 " Child Lead Benefits ( Sub 16). sas" Using data from the lead specific data set and results of the IEUBK lead model, child blood lead levels are estimated based on fish consumption in recreational and subsistence fishing families. 185 /* Program: Child Lead Benefits ( Sub 16). sas Author: Stuart Smith, Abt Associates, Inc. Date: 06/ 12/ 2000 Updated: 11/ 27/ 2002 This program estimates benefits from reduced consumption of fish tissue containing lead. Reductions in MP& M facility lead loadings result in lower probabilties of lead related illnesses. */

libname library ' H:\ ENVR\ Triton05 MP& M\ Benefits\ MP& M Phase 2 Analysis\ Data\ SAS Sets for Analysis\';

% let tbl_ path = % sysfunc( pathname( rundat));

title2 ' Child Lead Related Benefits ( Sub 16)';

data lead; set rundat. lead;

/* Combine the exposed populations, recreational and subsistence; fresh and salt water */

expspt = fwsptpop + swsptpop; expsub = fwsubpop + swsubpop;

/* Calculate the lead concentration in fish tissue available for intake */ /* = instream conc ( ug/ L)* biocon factor ( kg/ L)* ingestion rate( kg/ day) */

bcncspt = cnc_ base* 49* 0.00603; t1cncspt = cnc_ trt1* 49* 0.00603; t2cncspt = cnc_ trt2* 49* 0.00603;

bcncsub = cnc_ base* 49* 0.03033; t1cncsub = cnc_ trt1* 49* 0.03033; t2cncsub = cnc_ trt2* 49* 0.03033;

keep reachnmb expspt expsub bcncspt t1cncspt t2cncspt bcncsub t1cncsub t2cncsub event_ wt; run;

proc print; var reachnmb expspt expsub bcncspt t1cncspt t2cncspt bcncsub t1cncsub t2cncsub event_ wt; title3 ' Exposed Pops and Pb conc in fish tissue'; quit;

data lead2; set lead;

/* Determine sport blood levels */ /* cnc_ fmt. is a special format which translates the concentration in the fish tissue into an index value for a lookup table of corresponding blood lead levels */ 186 /* Get base pbp level */ lookup = input( put( round( bcncspt,. 1), cnc_ fmt.), 4.0); if lookup = 9999 then lookup = 1; set mpm2. ubklook4 point = lookup; pb_ bspt = bloodPb;

/* Get t1 pbp level */ lookup = input( put( round( t1cncspt,. 1), cnc_ fmt.), 4.0); if lookup = 9999 then lookup = 1; set mpm2. ubklook4 point = lookup; pb_ 1spt = bloodPb;

/* Get t2 pbp level */ lookup = input( put( round( t2cncspt,. 1), cnc_ fmt.), 4.0); if lookup = 9999 then lookup = 1; set mpm2. ubklook4 point = lookup; pb_ 2spt = bloodPb;

/* Determine subsistence blood levels */ /* Get base pbp level */ lookup = input( put( round( bcncsub,. 1), cnc_ fmt.), 4.0); if lookup = 9999 then lookup = 1; set mpm2. ubklook4 point = lookup; pb_ bsub = bloodPb;

/* Get t1 pbp level */ lookup = input( put( round( t1cncsub,. 1), cnc_ fmt.), 4.0); if lookup = 9999 then lookup = 1; set mpm2. ubklook4 point = lookup; pb_ 1sub = bloodPb;

/* Get t2 pbp level */ lookup = input( put( round( t2cncsub,. 1), cnc_ fmt.), 4.0); if lookup = 9999 then lookup = 1; set mpm2. ubklook4 point = lookup; pb_ 2sub = bloodPb;

run;

proc print; var reachnmb pb_ bspt pb_ 1spt pb_ 2spt pb_ bsub pb_ 1sub pb_ 2sub; title3 ' Blood lead levels'; quit;

data lead3; set lead2;

kidspct = 0.1177; /* Percent of kids aged 0 7, U. S. Statistical Abstract 1997 */

popspt = expspt* kidspct; popsub = expsub* kidspct;

GM1_ spt = pb_ bspt  pb_ 1spt; 187 GM2_ spt = pb_ bspt  pb_ 2spt;

GM1_ sub = pb_ bsub  pb_ 1sub; GM2_ sub = pb_ bsub  pb_ 2sub;

/* Avoided lost IQ points */ avIQ1spt = ( GM1_ spt)* 0.25*( popspt/ 7); avIQ2spt = ( GM2_ spt)* 0.25*( popspt/ 7);

avIQ1sub = ( GM1_ sub)* 0.25*( popsub/ 7); avIQ2sub = ( GM2_ sub)* 0.25*( popsub/ 7);

IQ1_ spt = avIQ1spt * 9419; IQ2_ spt = avIQ2spt * 9419;

IQ1_ sub = avIQ1sub * 9419; IQ2_ sub = avIQ2sub * 9419;

/* Fewer children with IQ below 70 */

if GM1_ spt > 0 then do; num1spt = ( 70  ( 100 + 0.25* GM1_ spt))/ 15; sv701spt = 0.02275  cdf(' NORMAL', num1spt); b70_ 1spt = sv701spt * ( popspt/ 7) * 58012; end; else do; sv701spt = 0; b70_ 1spt = 0; end;

if GM2_ spt > 0 then do; num2spt = ( 70  ( 100 + 0.25* GM2_ spt))/ 15; sv702spt = 0.02275  cdf(' NORMAL', num2spt); b70_ 2spt = sv702spt * ( popspt/ 7) * 58012; end; else do; sv702spt = 0; b70_ 2spt = 0; end;

if GM1_ sub > 0 then do; num1sub = ( 70  ( 100 + 0.25* GM1_ sub))/ 15; sv701sub = 0.02275  cdf(' NORMAL', num1sub); b70_ 1sub = sv701sub * ( popsub/ 7) * 58012; end; else do; sv701sub = 0; b70_ 1sub = 0; end;

if GM2_ sub > 0 then do; num2sub = ( 70  ( 100 + 0.25* GM2_ sub))/ 15; sv702sub = 0.02275  cdf(' NORMAL', num2sub); 188 b70_ 2sub = sv702sub * ( popsub/ 7) * 58012; end; else do; sv702sub = 0; b70_ 2sub = 0; end;

/* Number of cases of IQ< 70 saved */

c70_ 1spt = sv701spt * ( popspt/ 7); c70_ 2spt = sv702spt * ( popspt/ 7); c70_ 1sub = sv701sub * ( popsub/ 7); c70_ 2sub = sv702sub * ( popsub/ 7);

check1 = pdf(' LOGNORMAL', 20, pb_ bspt, 2.09); check2 = pdf(' LOGNORMAL', 20, pb_ bsub, 2.09);

spt1_ 20=( check1  pdf(' LOGNORMAL', 20, pb_ 1spt, 2.09))* 0.2*( popspt/ 7)* 16485; spt2_ 20=( check1  pdf(' LOGNORMAL', 20, pb_ 2spt, 2.09))* 0.2*( popspt/ 7)* 16485; sub1_ 20=( check2  pdf(' LOGNORMAL', 20, pb_ 1sub, 2.09))* 0.2*( popsub/ 7)* 16485; sub2_ 20=( check2  pdf(' LOGNORMAL', 20, pb_ 2sub, 2.09))* 0.2*( popsub/ 7)* 16485;

red201sp=( check1  pdf(' LOGNORMAL', 20, pb_ 1spt, 2.09))* 0.2*( popspt/ 7); red202sp=( check1  pdf(' LOGNORMAL', 20, pb_ 2spt, 2.09))* 0.2*( popspt/ 7); red201sb=( check2  pdf(' LOGNORMAL', 20, pb_ 1sub, 2.09))* 0.2*( popsub/ 7); red202sb=( check2  pdf(' LOGNORMAL', 20, pb_ 2sub, 2.09))* 0.2*( popsub/ 7); run;

proc means sum data= lead3; var avIQ1spt avIQ1sub avIQ2spt avIQ2sub IQ1_ spt IQ1_ sub IQ2_ spt IQ2_ sub b70_ 1spt b70_ 1sub b70_ 2spt b70_ 2sub sv701spt sv702spt sv701sub sv702sub c70_ 1spt c70_ 2spt c70_ 1sub c70_ 2sub spt1_ 20 spt2_ 20 sub1_ 20 sub2_ 20 red201sp red202sp red201sb red202sb popspt popsub; title3 ' Benefits ( sample)'; quit;

proc means sum data= lead3; weight event_ wt; var avIQ1spt avIQ1sub avIQ2spt avIQ2sub IQ1_ spt IQ1_ sub IQ2_ spt IQ2_ sub b70_ 1spt b70_ 1sub b70_ 2spt b70_ 2sub sv701spt sv702spt sv701sub sv702sub c70_ 1spt c70_ 2spt c70_ 1sub c70_ 2sub spt1_ 20 spt2_ 20 sub1_ 20 sub2_ 20 red201sp red202sp red201sb red202sb popspt popsub; output out = rundat. chld_ ld sum = avIQ1spt avIQ1sub avIQ2spt avIQ2sub IQ1_ spt IQ1_ sub IQ2_ spt IQ2_ sub 189 b70_ 1spt b70_ 1sub b70_ 2spt b70_ 2sub sv701spt sv702spt sv701sub sv702sub c70_ 1spt c70_ 2spt c70_ 1sub c70_ 2sub spt1_ 20 spt2_ 20 sub1_ 20 sub2_ 20 red201sp red202sp red201sb red202sb popspt popsub; title3 ' Benefits ( national)'; quit;

data rundat. chld_ tbl ( keep= gender categ type exppop case1 case2 benlow1 benlow2 benmid1 benmid2 benhigh1 benhigh2); set rundat. chld_ ld;

length gender $ 50 categ $ 50 type $ 50;

array IQ[ 8] IQ1_ spt avIQ1spt IQ2_ spt avIQ2spt IQ1_ sub avIQ1sub IQ2_ sub avIQ2sub; array RD[ 8] b70_ 1spt c70_ 1spt b70_ 2spt c70_ 2spt b70_ 1sub c70_ 1sub b70_ 2sub c70_ 2sub; array PB[ 8] spt1_ 20 red201sp spt2_ 20 red202sp sub1_ 20 red201sb sub2_ 20 red202sb;

gender = " Children";

if _ n_ = 1 then do; do i = 1 to 6; if i = 1 then do; categ = " Avoided IQ Loss"; type = " Sport"; case1 = IQ[ 2]; case2 = IQ[ 4]; benlow1 = IQ[ 1]; benlow2 = IQ[ 3]; benmid1 = benlow1; benmid2 = benlow2; benhigh1 = benlow1; benhigh2 = benlow2; output; end; if i = 2 then do; categ = " Avoided IQ Loss"; type = " Sub"; case1 = IQ[ 6]; case2 = IQ[ 8]; benlow1 = IQ[ 5]; benlow2 = IQ[ 7]; benmid1 = benlow1; benmid2 = benlow2; benhigh1 = benlow1; benhigh2 = benlow2; output; end; if i = 3 then do; categ = " Reduced IQ < 70"; type = " Sport"; case1 = RD[ 2]; case2 = RD[ 4]; benlow1 = RD[ 1]; benlow2 = RD[ 3]; benmid1 = benlow1; benmid2 = benlow2; benhigh1 = benlow1; benhigh2 = benlow2; output; end; if i = 4 then do; categ = " Reduced IQ < 70"; type = " Sub"; 190 case1 = RD[ 6]; case2 = RD[ 8]; benlow1 = RD[ 5]; benlow2 = RD[ 7]; benmid1 = benlow1; benmid2 = benlow2; benhigh1 = benlow1; benhigh2 = benlow2; output; end; if i = 5 then do; categ = " Reduced Pb > 20"; type = " Sport"; case1 = PB[ 2]; case2 = PB[ 4]; benlow1 = PB[ 1]; benlow2 = PB[ 3]; benmid1 = benlow1; benmid2 = benlow2; benhigh1 = benlow1; benhigh2 = benlow2; output; end; if i = 6 then do; categ = " Reduced Pb > 20"; type = " Sub"; case1 = PB[ 6]; case2 = PB[ 8]; benlow1 = PB[ 5]; benlow2 = PB[ 7]; benmid1 = benlow1; benmid2 = benlow2; benhigh1 = benlow1; benhigh2 = benlow2; output; end; end; end; run;

data rundat. chld_ all; set rundat. neo_ tbl rundat. chld_ tbl; run;

proc summary data= rundat. chld_ all; var benlow1 benmid1 benhigh1 benlow2 benmid2 benhigh2; output out = rundat. chld_ tot ( drop= _ type_ _ freq_) sum = benlow1 benmid1 benhigh1 benlow2 benmid2 benhigh2; run;

data exlchld ( keep= titlecol col1 col11 label=" Child Lead Benefits  Weighted"); set rundat. chld_ all rundat. chld_ tot; length titlecol $ 50 col1  col11 $ 128;

if _ n_ = 1 then do; titlecol = " Table: Child Lead Benefits, Weighted Estimates"; output; titlecol = " Run on & sysdate"; output; titlecol = " "; output; col4 = " Option 1"; col8 = " Option 2"; output; col2 = " Category"; col3 = " Type"; col4 = " Reduced Cases"; col5 = " Low"; col6 = " Mid"; 191 col7 = " High"; col8 = " Reduced Cases"; col9 = " Low"; col10 = " Mid"; col11 = " High"; output; end;

if _ n_ = 9 then col1 = " Total"; else col1 = gender;

col2 = trim( categ); col3 = trim( type); col4 = trim( case1); col5 = trim( benlow1); col6 = trim( benmid1); col7 = trim( benhigh1); col8 = trim( case2); col9 = trim( benlow2); col10 = trim( benmid2); col11 = trim( benhigh2); output; run;

proc export data = exlchld outfile = " C:\ Output Files\ Child Lead. xls" dbms= excel2000 replace; run;

/* End Child Lead Benefits ( Sub 16). sas */ 192 Program Output 193 MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 114 Child Lead Related Benefits ( Sub 16) Benefits ( national) 00: 44 Sunday, February 9, 2003 The MEANS Procedure Variable Sum                          avIQ1spt 796.2802134 avIQ1sub 282.1023758 avIQ2spt 905.6883457 avIQ2sub 333.8142973 IQ1_ spt 7500163.33 IQ1_ sub 2657122.28 IQ2_ spt 8530678.53 IQ2_ sub 3144196.87 b70_ 1spt 161062.44 b70_ 1sub 54944.28 b70_ 2spt 183822.10 b70_ 2sub 65535.40 sv701spt 0.0484820 sv702spt 0.0499591 sv701sub 0.2426528 sv702sub 0.2584448 c70_ 1spt 2.7763642 c70_ 2spt 3.1686909 c70_ 1sub 0.9471192 c70_ 2sub 1.1296869 spt1_ 20  14036.97 spt2_ 20  15134.93 sub1_ 20  5361.08 sub2_ 20  6368.32 red201sp  0.8514995 red202sp  0.9181031 red201sb  0.3252095 red202sb  0.3863099 popspt 189145.51 popsub 9456.66                          194 Recreational Benefits 195 Recreational Fishing 196 " Human Health AWQC Exceedences ( Sub 2). sas" The number of human health based AWQC exceedences for each RF1 reach are estimated and output to a data set. See Human Health Benefits Cancer from Fish Consumption for program code and output. 197 " Aquatic Life AWQC Exceedences ( Sub 3). sas" The number of aquatic life based AWQC exceedences for each RF1 reach are estimated and output to a data set. 198 /* Program: Aquatic Life AWQC ( Sub 3). sas Author: Stuart Smith ( based on Elena Besedin, Rob Sartain) , Abt. Associates, Inc. Date Modified: 10/ 31/ 2002 ( 1) IDENTIFIES MP& M DISCHARGE REACHES IN WHICH POLLUTANT CONCENTRATIONS EXCEED THE AWQC VALUES FOR AQUATIC LIFE; ( 2) EXTRAPOLATES THE RESULTS TO THE POPULATION LEVEL BY USING DIFFERENTIAL WEIGHTING.

THE TOTAL NUMBER OF EXCEEDENCE EVENTS IS THEN ESTIMATED BY SUMMING OVER THE RESULTS FROM THE FIRST AND SECOND RUN. */

title2 ' Output from Aquatic Life ( Sub 3)';

/* DISCHARGE EVENTS, ASSOCIATED LOADINGS AND WEIGHTS HAVE BEEN CALCULATED BY ANOTHER PROGRAM. ALL UNIT CONVERSIONS HAVE OCCURRED. SEE " Generate Discharge Events and Loadings ( Sub 1). sas" FOR SPECIFIC PROGRAM INFORMATION. */

/* PULL TOGETHER THE EVENTS AND CHEMICAL SPECIFIC PARAMETERS */

Proc sort data = mpm2. chemdata; by casno; quit;

Proc sql; create table stream as select reachnmb, event_ ld. casno, evnt_ flo, base, treat1, treat2, event_ wt, segl, mnflo as mean_ flo, svten as flo7q10, acute as adf, chronic as cdf, dcp, wbtype, biodegre as decay, surface, drainage, volume from rundat. event_ ld, mpm2. chemdata where event_ ld. casno = chemdata. casno order by reachnmb, casno; quit;

/******************************************************************************/ /* */ /* FRESHWATER REACHES  CONCENTRATIONS AND AWQC EXCEEDENCES */ /* */ /******************************************************************************/

/* Calculate chemical concentrations in lakes */

data Lake; set stream; if wbtype = ' lake'; if decay;

keep reachnmb casno base treat1 treat2 adf cdf restime decay event_ wt;

/* drainage area and surface area available */ if surface ^= . and drainage ^= . then do; num = drainage; /* drainage area in sq. miles */ 199 denom = surface* 0.0015625; /* convert surface area from acres to sq. miles */ logdays = 4.077  ( 1.177* log10( num/ denom)); restime = 10** logdays; time1 = restime* 24; /* CONVERT HYDRAULIC TIME FROM DAYS INTO HOURS */ cnc_ aqcb = ( base/ evnt_ flo)/( 1+ time1* decay); cnc_ aqc1 = ( treat1/ evnt_ flo)/( 1+ time1* decay); cnc_ aqc2 = ( treat2/ evnt_ flo)/( 1+ time1* decay); good = 1; end;

/* missing area measure( s), do end of pipe value = loadings/ flow */ else do; cnc_ aqcb = ( base/ evnt_ flo); cnc_ aqc1 = ( treat1/ evnt_ flo); cnc_ aqc2 = ( treat2/ evnt_ flo); end;

cnc_ aqab = cnc_ aqcb; cnc_ aqa1 = cnc_ aqc1; cnc_ aqa2 = cnc_ aqc2;

keep reachnmb casno cnc_ aqcb cnc_ aqc1 cnc_ aqc2 cnc_ aqab cnc_ aqa1 cnc_ aqa2 event_ wt; run;

/* Great Lakes reaches */

data BigLake; set stream; if wbtype = ' glake';

if cdf ^=. then do; cnc_ aqcb = base/( evnt_ flo* cdf); cnc_ aqc1 = treat1/( evnt_ flo* cdf); cnc_ aqc2 = treat2/( evnt_ flo* cdf); end; else if dcp^=. then do; benchld= 10000* 907.2; /* benchmark load is 10000 per year; 907.2 is coversion factor from kg to ton*/ cnc_ aqcb=(( base/ 10** 9)* 365* dcp)/ benchld; cnc_ aqc1=(( treat1/ 10** 9)* 365* dcp)/ benchld; cnc_ aqc2=(( treat2/ 10** 9)* 365* dcp)/ benchld; end; else do; cnc_ aqcb = .; cnc_ aqc1 = .; cnc_ aqc2 = .; end; cnc_ aqab = cnc_ aqcb; cnc_ aqa1 = cnc_ aqc1; cnc_ aqa2 = cnc_ aqc2;

if cnc_ aqcb=. then delete; 200 keep reachnmb casno cnc_ aqcb cnc_ aqc1 cnc_ aqc2 cnc_ aqab cnc_ aqa1 cnc_ aqa2 event_ wt; run;

/* Stream reaches */ /* Calculate 1Q10. Calculate in stream concentrations based on 7Q10 and 1Q10*/

data stream2; set stream; if wbtype = ' stream';

/* Account for missing 7q10 data, averge */ if flo7q10 = . or flo7q10 = 0 then flo7q10 = 909.7836; if mean_ flo = . or mean_ flo = 0 then mean_ flo = 6711.2411;

/* Convert flo7q10 from ft3/ sec to l/ day using conversion factor 2446576 */

svnq10 = flo7q10* 2446576;

/* Calculate 1q10 flow */

flo1q10= 0.843*( svnq10** 0.993);

/* all flows in l/ day */

cnc_ aqcb = base/( evnt_ flo+ svnq10); cnc_ aqc1 = treat1/( evnt_ flo+ svnq10); cnc_ aqc2 = treat2/( evnt_ flo+ svnq10); cnc_ aqab = base/( evnt_ flo+ flo1q10); cnc_ aqa1 = treat1/( evnt_ flo+ flo1q10); cnc_ aqa2 = treat2/( evnt_ flo+ flo1q10);

keep reachnmb casno cnc_ aqcb cnc_ aqc1 cnc_ aqc2 cnc_ aqab cnc_ aqa1 cnc_ aqa2 event_ wt; run;

proc sort data= stream2; by reachnmb casno; quit;

/* Combine small and large lakes with streams */

data freshwat; set stream2 BigLake Lake; run;

proc sort data= freshwat; by reachnmb casno; quit;

/* Compare in stream pollutant concentrations with freshwater AWQCs for aquatic life */ proc sql;

/* Merge in freshwater aquatic life limits. */ 201 create table awqc1 as select reachnmb, freshwat. casno, chemical, cnc_ aqcb, cnc_ aqc1, cnc_ aqc2, cnc_ aqab, cnc_ aqa1, cnc_ aqa2, faval__ i as faval, fcval__ i as fcval, event_ wt from freshwat, mpm2. chemdata where freshwat. casno = chemdata. casno order by reachnmb, casno; quit;

/* Compare concentrations to awqc limits and determine exceedences */

data stream4; set awqc1;

if cnc_ aqcb > fcval and fcval > 0 then chronicb = 1; if cnc_ aqc1 > fcval and fcval > 0 then chronic1 = 1; if cnc_ aqc2 > fcval and fcval > 0 then chronic2 = 1; if cnc_ aqcb > faval and faval > 0 then acuteb = 1; if cnc_ aqc1 > faval and faval > 0 then acute1 = 1; if cnc_ aqc2 > faval and faval > 0 then acute2 = 1;

keep reachnmb casno chemical chronicb chronic1 chronic2 acuteb acute1 acute2 event_ wt; run;

/***********************************************************************/ /* END FRESHWATER REACHES */ /***********************************************************************/

/******************************************************************************/ /* */ /* SALTWATER REACHES  CONCENTRATIONS AND AWQC EXCEEDENCES */ /* */ /******************************************************************************/

/* Calculate pollutant concentrations in Oceans */

data marine; set stream; if wbtype = ' bay';

if cdf ^=. then do; cnc_ aqcb = base/( evnt_ flo* cdf); cnc_ aqc1 = treat1/( evnt_ flo* cdf); cnc_ aqc2 = treat2/( evnt_ flo* cdf); end; else if dcp ^=. then do; benchld= 10000* 907.2; /* benchmark load is 10000 per year; 907.2 is coversion factor from kg to ton */

/* facility loadings were given in ug/ day */

cnc_ aqcb=(( base/ 10** 9)* 365* dcp)/ benchld; 202 cnc_ aqc1=(( treat1/ 10** 9)* 365* dcp)/ benchld; cnc_ aqc2=(( treat2/ 10** 9)* 365* dcp)/ benchld; end; else do; cnc_ aqcb = .; cnc_ aqc1 = .; cnc_ aqc2 = .; end;

cnc_ aqab = cnc_ aqcb; cnc_ aqa1 = cnc_ aqc1; cnc_ aqa2 = cnc_ aqc2;

if cnc_ aqcb=. then delete; keep reachnmb casno cnc_ aqcb cnc_ aqc1 cnc_ aqc2 cnc_ aqab cnc_ aqa1 cnc_ aqa2 event_ wt; run;

/* Get AWQC limits for marine environment */

proc sql;

/* Merge in Aquatic life saltwater awqc limits */

create table marine1 as select reachnmb, marine. casno, chemical, cnc_ aqcb, cnc_ aqc1, cnc_ aqc2, cnc_ aqab, cnc_ aqa1, cnc_ aqa2, event_ wt, maval__ i as maval, mcval__ i as mcval, count( distinct reachnmb) as reaches from marine, mpm2. chemdata where marine. casno = chemdata. casno order by reachnmb, casno; quit;

/* Compare concentrations to awqc limits and determine exceedences */

data marine2; set marine1; if cnc_ aqcb > mcval and mcval > 0 then chronicb = 1; if cnc_ aqc1 > mcval and mcval > 0 then chronic1 = 1; if cnc_ aqc2 > mcval and mcval > 0 then chronic2 = 1; if cnc_ aqcb > maval and maval > 0 then acuteb = 1; if cnc_ aqc1 > maval and maval > 0 then acute1 = 1; if cnc_ aqc2 > maval and maval > 0 then acute2 = 1;

keep reachnmb casno chemical chronicb chronic1 chronic2 acuteb acute1 acute2 reaches event_ wt; run;

/* Combine exceedences from fresh and saltwater */ data allreach; set marine2 stream4;

keep reachnmb casno chemical chronicb chronic1 chronic2 acuteb acute1 acute2 reaches event_ wt; run; 203 proc sql;

/* Estimate AWQCs exceedance by chemicals on a sample and national basis*/ /* awt replaced with 1 */ create table chm_ awqc as select chemical, sum( chronicb) as crb, sum( chronic1) as cr1, sum( chronic2) as cr2, sum( acuteb) as acb, sum( acute1) as ac1, sum( acute2) as ac2, sum( chronicb* event_ wt* 1) as ntlbw, sum( chronic1* event_ wt* 1) as ntl1w, sum( chronic2* event_ wt* 1) as ntl2w, sum( acuteb* event_ wt* 1) as ntlbo, sum( acute1* event_ wt* 1) as ntl1o, sum( acute2* event_ wt* 1) as ntl2o, case when sum( chronicb) >= 1 then 1 else 0 end as chemcrb, case when sum( chronic1) >= 1 then 1 else 0 end as chemcr1, case when sum( chronic2) >= 1 then 1 else 0 end as chemcr2, case when sum( acuteb) >= 1 then 1 else 0 end as chemacb, case when sum( acute1) >= 1 then 1 else 0 end as chemac1, case when sum( acute2) >= 1 then 1 else 0 end as chemac2 from allreach where chronicb> 0 or acuteb> 0 group by chemical order by chemical;

/* Estimate aquatic life AWQC Exceedences by option/ reach */

create table aq_ life2 as select allreach. reachnmb, sum( chronicb) as test_ cb, sum( chronic1) as test_ c1, sum( chronic2) as test_ c2, sum ( acuteb) as test_ ab, sum( acute1) as test_ a1, sum ( acute2) as test_ a2, sum( chronicb* event_ wt* 1) as ntlb1, sum( chronic1* event_ wt* 1) as ntnl11, sum( chronic2* event_ wt* 1) as ntnl21, sum( acuteb * event_ wt* 1) as ntlb2, sum( acute1* event_ wt* 1) as ntnl12, sum( acute2* event_ wt* 1) as ntnl22, case when sum( chronicb) >= 1 then 1 else 0 end as rchcrb, case when sum( chronic1) >= 1 then 1 else 0 end as rchcr1, case when sum( chronic2) >= 1 then 1 else 0 end as rchcr2, case when sum( acuteb) >= 1 then 1 else 0 end as rchacb, case when sum( acute1) >= 1 then 1 else 0 end as rchac1, case when sum( acute2) >= 1 then 1 else 0 end as rchac2, sum( event_ wt) as reachwt from allreach where chronicb> 0 or acuteb> 0 group by reachnmb order by reachnmb;

/* CALCULATE THE TOTAL NUMBER OF AFFECTED REACHES ON A NATIONAL BASIS awt replaced with a 1 */ create table rch1 as select 204 sum( event_ wt) as wght1, reachnmb, chemical from allreach where chronicb > 0 group by reachnmb, chemical order by reachnmb, chemical;

create table chrnb as select max( wght1) as chrnc_ b, reachnmb from rch1 group by reachnmb order by reachnmb;

create table rch2 as select sum( event_ wt) as wght2, reachnmb, chemical from allreach where acuteb > 0 group by reachnmb, chemical order by reachnmb, chemical;

create table actb as select max( wght2) as acutr_ b, reachnmb from rch2 group by reachnmb order by reachnmb;

create table rch3 as select sum( event_ wt) as wght3, reachnmb, chemical from allreach where chronic1 > 0 group by reachnmb, chemical order by reachnmb, chemical;

create table chrn1 as select max( wght3) as chrnc_ 1, reachnmb from rch3 group by reachnmb order by reachnmb;

create table rch4 as select sum( event_ wt) as wght4, reachnmb, chemical from allreach where acute1 > 0 group by reachnmb, chemical order by reachnmb, chemical;

create table act1 as select max( wght4) as acutr_ 1, reachnmb from rch4 group by reachnmb order by reachnmb; 205 create table rch5 as select sum( event_ wt) as wght5, reachnmb, chemical from allreach where chronic2 > 0 group by reachnmb, chemical order by reachnmb, chemical;

create table chrn2 as select max( wght5) as chrnc_ 2, reachnmb from rch5 group by reachnmb order by reachnmb;

create table rch6 as select sum( event_ wt) as wght6, reachnmb, chemical from allreach where acute2 > 0 group by reachnmb, chemical order by reachnmb, chemical;

create table act2 as select max( wght6) as acutr_ 2, reachnmb from rch6 group by reachnmb order by reachnmb;

create table rch7 as select sum( event_ wt) as wght7, reachnmb, chemical from allreach where chronicb > 0 or acuteb > 0 group by reachnmb, chemical order by reachnmb, chemical;

create table multi_ b as select max( wght7) as multi_ b, reachnmb from rch7 group by reachnmb order by reachnmb;

create table rch8 as select sum( event_ wt) as wght8, reachnmb, chemical from allreach where chronic1 > 0 or acute1 > 0 group by reachnmb, chemical order by reachnmb, chemical;

create table multi_ 1 as select max( wght8) as multi_ 1, reachnmb from rch8 206 group by reachnmb order by reachnmb;

create table rch9 as select sum( event_ wt) as wght9, reachnmb, chemical from allreach where chronic2 > 0 or acute2 > 0 group by reachnmb, chemical order by reachnmb, chemical;

create table multi_ 2 as select max( wght9) as multi_ 2, reachnmb from rch9 group by reachnmb order by reachnmb; quit;

data rchwts; merge chrnb actb chrn1 act1 chrn2 act2 multi_ b multi_ 1 multi_ 2; by reachnmb;

if chrnc_ b=. then chrnc_ b= 0; if chrnc_ 1=. then chrnc_ 1= 0; if chrnc_ 2=. then chrnc_ 2= 0; if acutr_ b=. then acutr_ b= 0; if acutr_ 1=. then acutr_ 1= 0; if acutr_ 2=. then acutr_ 2= 0; if multi_ b=. then multi_ b= 0; if multi_ 1=. then multi_ 1= 0; if multi_ 2=. then multi_ 2= 0; run;

quit;

/* Output the affected reaches for the recreational benefits programs */

data rundat. rec_ al; set allreach; if chronicb > 0 or acuteb > 0; drop reaches; run;

/* Print out results */

proc means sum data= aq_ life2; var test_ cb test_ c1 test_ c2 test_ ab test_ a1 test_ a2 ntlb1 ntnl11 ntnl21 ntlb2 ntnl12 ntnl22 rchcrb rchcr1 rchcr2 rchacb rchac1 rchac2; label reachnmb =' Reach ID' test_ cb = ' Chronic AWQC Exeedences ( Smpl/ BL)' test_ c1 = ' Chronic AWQC Exeedences ( Smpl/ Opt 1)' test_ c2 = ' Chronic AWQC Exeedences ( Smpl/ Opt 2)' test_ ab = ' Acute AWQC Exeedences ( Smpl/ BL)' 207 test_ a1 = ' Acute AWQC Exeedences ( Smpl/ Opt 1)' test_ a2 = ' Acute AWQC Exeedences ( Smpl/ Opt 2)' ntlb1 = ' Chronic AWQC Exeedences ( Ntl/ BL)' ntnl11 = ' Chronic AWQC Exeedences ( Ntl/ OPt 1)' ntnl21 = ' Chronic AWQC Exeedences ( Ntl/ OPt 2)' ntlb2 = ' Acute AWQC Exeedences ( Ntl/ Bl)' ntnl12 = ' Acute AWQC Exeedences ( Ntl/ Opt 1)' ntnl22 = ' Acute AWQC Exeedences ( Ntl/ Opt 2)' rchcrb = ' Number of Chronic reaches ( BL)' rchcr1 = ' Number of Chronic reaches ( Opt 1)' rchcr2 = ' Number of Chronic reaches ( Opt 2)' rchacb = ' Number of Acute reaches ( BL)' rchac1 = ' Number of Acute reaches ( Opt 1)' rchac2 = ' Number of Acute reaches ( Opt 2)';

title3 ' Aquatic Life AWQC Exceedences by Option ( Sample/ National)';

proc means sum data= chm_ awqc; var crb cr1 cr2 acb ac1 ac2 ntlbw ntl1w ntl2w ntlbo ntl1o ntl2o chemcrb chemcr1 chemcr2 chemacb chemac1 chemac2; label crb = ' Chronic AWQC Exeedences ( Smpl/ BL)' cr1 = ' Chronic AWQC Exeedences ( Smpl/ Opt 1)' cr2 = ' Chronic AWQC Exeedences ( Smpl/ Opt 2)' acb = ' Acute AWQC Exeedences ( Smpl/ BL)' ac1 = ' Acute AWQC Exeedences ( Smpl/ Opt 1)' ac2 = ' Acute AWQC Exeedences ( Smpl/ Opt 2)' ntlbw = ' Chronic AWQC Exeedences ( Ntl/ BL)' ntl1w = ' Chronic AWQC Exeedences ( Ntl/ OPt 1)' ntl2w = ' Chronic AWQC Exeedences ( Ntl/ OPt 2)' ntlbo = ' Acute AWQC Exeedences ( Ntl/ Bl)' ntl1o = ' Acute AWQC Exeedences ( Ntl/ Opt 1)' ntl2o = ' Acute AWQC Exeedences ( Ntl/ Opt 2)' chemcrb = ' Number of Chronic pol. ( BL)' chemcr1 = ' Number of Chronic pol. ( Opt 1)' chemcr2 = ' Number of Chronic pol. ( Opt 2)' chemacb = ' Number of Acute pol. ( BL)' chemac1 = ' Number of Acute pol. ( Opt 1)' chemac2 = ' Number of Acute pol. ( Opt 2)';

title3 ' AWQC Exceedances by pollutant/ option ( Sample / National)';

proc means sum data= rchwts; var chrnc_ b chrnc_ 1 chrnc_ 2 acutr_ b acutr_ 1 acutr_ 2 multi_ b multi_ 1 multi_ 2; label chrnc_ b = '# of Reaches Exceeding Chronic AWQC ( Bl) ' chrnc_ 1 = '# of Reaches Exceeding Chronic AWQC ( Opt 1)' chrnc_ 2 = '# of Reaches Exceeding Chronic AWQC ( Opt 2)' acutr_ b = '# of Reaches Exceeding Acute AWQC ( Bl)' acutr_ 1 = '# of Reaches Exceeding Acute AWQC ( Opt 1)' acutr_ 2 = '# of Reaches Exceeding Acute AWQC ( Opt 2)' multi_ b = ' Either ( BL)' multi_ 1 = ' Either ( Opt 1)' multi_ 2 = ' Either ( Opt 2)'; title3 ' IL Reaches Exceeding AWQC Limits for Aquatic Life ( National)'; 208 run;

/* Sum the number of exceedence events by reach and event weight */ proc sql; create table preme as select reachnmb, event_ wt, sum( chronicb) as testb1, sum( chronic1) as test1_ 1, sum( chronic2) as test2_ 1, sum( acuteb) as testb2, sum( acute1) as test1_ 2, sum( acute2) as test2_ 2 from allreach group by reachnmb, event_ wt order by reachnmb, event_ wt; quit;

/* Determine the number of reaches which eliminate all exceedences and the number of reaches which only partially reduce */

data reduc; set preme; if testb1=. then testb1= 0; if test1_ 1=. then test1_ 1= 0; if test2_ 1=. then test2_ 1= 0; if testb2=. then testb2= 0; if test1_ 2=. then test1_ 2= 0; if test2_ 2=. then test2_ 2= 0;

if testb1 >= 1 and test1_ 1 = 0 then elim1_ 1 = 1; else elim1_ 1= 0; if testb1 >= 1 and test2_ 1 = 0 then elim2_ 1 = 1; else elim2_ 1= 0; if testb2 >= 1 and test1_ 2 = 0 then elim1_ 2 = 1; else elim1_ 2= 0; if testb2 >= 1 and test2_ 2 = 0 then elim2_ 2 = 1; else elim2_ 2= 0;

if testb1 >= 1 and test1_ 1 < testb1 and elim1_ 1 = 0 then red1_ 1 = 1; else red1_ 1= 0; if testb1 >= 1 and test2_ 1 < testb1 and elim2_ 1 = 0 then red2_ 1 = 1; else red2_ 1= 0; if testb2 >= 1 and test1_ 2 < testb2 and elim1_ 2 = 0 then red1_ 2 = 1; else red1_ 2= 0; if testb2 >= 1 and test2_ 2 < testb2 and elim2_ 2 = 0 then red2_ 2 = 1; else red2_ 2= 0;

if elim1_ 1= 1 or elim1_ 2= 1 then elim1 = 1; else elim1 = 0; if elim2_ 1= 1 or elim2_ 2= 1 then elim2 = 1; else elim2 = 0; if red1_ 1 = 1 or red1_ 2 = 1 then reduc1 = 1; else reduc1 = 0; if red2_ 1 = 1 or red2_ 2 = 1 then reduc2 = 1; else reduc2 = 0; run;

/* Print results */

proc means sum data= reduc; var elim1_ 1 elim2_ 1 elim1_ 2 elim2_ 2 red1_ 1 red2_ 1 red1_ 2 red2_ 2 reduc1 reduc2; label elim1_ 1 = ' Elim. Chronic ( Opt1)' elim2_ 1 = ' Elim. Chronic ( Opt2)' 209 elim1_ 2 = ' Elim. Acute ( Opt1)' elim2_ 2 = ' Elim. Acute ( Opt2)' red1_ 1 = ' Redu. Chronic ( Opt1)' red2_ 1 = ' Redu. Chronic ( Opt2)' red1_ 2 = ' Redu. Acute ( Opt1)' red2_ 2 = ' Redu. Acute ( Opt2)' reduc1 = ' Total Reaches Reducing ( Opt1)' reduc2 = ' Total Reaches Reducing ( Opt2)'; title3 ' AWQC Eliminations and Reductions ( Sample)'; title4 ' Number of Sample Reaches'; quit;

proc means sum data= reduc; weight event_ wt; var elim1_ 1 elim2_ 1 elim1_ 2 elim2_ 2 red1_ 1 red2_ 1 red1_ 2 red2_ 2 reduc1 reduc2; label elim1_ 1 = ' Elim. Chronic ( Opt1)' elim2_ 1 = ' Elim. Chronic ( Opt2)' elim1_ 2 = ' Elim. Acute ( Opt1)' elim2_ 2 = ' Elim. Acute ( Opt2)' red1_ 1 = ' Redu. Chronic ( Opt1)' red2_ 1 = ' Redu. Chronic ( Opt2)' red1_ 2 = ' Redu. Acute ( Opt1)' red2_ 2 = ' Redu. Acute ( Opt2)' reduc1 = ' Total Reaches Reducing ( Opt1)' reduc2 = ' Total Reaches Reducing ( Opt2)'; title3 ' AWQC Eliminations and Reductions ( National)'; title4 ' Number of National Reaches'; quit;

proc print data= aq_ life2; quit; proc print data= chm_ awqc; title3 ' Exceedences by Chemical'; quit;

/* End Aquatic Life AWQC ( Sub 3). sas */ 210 Program Output 211 MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 12 Output from Aquatic Life ( Sub 3) 00: 44 Sunday, February 9, 2003 Aquatic Life AWQC Exceedences by Option ( Sample/ National)

The MEANS Procedure Variable Label Sum                                                                  test_ cb Chronic AWQC Exeedences ( Smpl/ BL) 307.0000000 test_ c1 Chronic AWQC Exeedences ( Smpl/ Opt 1) 152.0000000 test_ c2 Chronic AWQC Exeedences ( Smpl/ Opt 2) . test_ ab Acute AWQC Exeedences ( Smpl/ BL) 90.0000000 test_ a1 Acute AWQC Exeedences ( Smpl/ Opt 1) 34.0000000 test_ a2 Acute AWQC Exeedences ( Smpl/ Opt 2) . ntlb1 Chronic AWQC Exeedences ( Ntl/ BL) 2582.49 ntnl11 Chronic AWQC Exeedences ( Ntl/ OPt 1) 1368.52 ntnl21 Chronic AWQC Exeedences ( Ntl/ OPt 2) . ntlb2 Acute AWQC Exeedences ( Ntl/ Bl) 630.5351378 ntnl12 Acute AWQC Exeedences ( Ntl/ Opt 1) 254.3570580 ntnl22 Acute AWQC Exeedences ( Ntl/ Opt 2) . rchcrb Number of Chronic reaches ( BL) 61.0000000 rchcr1 Number of Chronic reaches ( Opt 1) 29.0000000 rchcr2 Number of Chronic reaches ( Opt 2) 0 rchacb Number of Acute reaches ( BL) 33.0000000 rchac1 Number of Acute reaches ( Opt 1) 8.0000000 rchac2 Number of Acute reaches ( Opt 2) 0                                                                  MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 13 Output from Aquatic Life ( Sub 3) 00: 44 Sunday, February 9, 2003 AWQC Exceedances by pollutant/ option ( Sample / National)

The MEANS Procedure Variable Label Sum                                                                  crb Chronic AWQC Exeedences ( Smpl/ BL) 307.0000000 cr1 Chronic AWQC Exeedences ( Smpl/ Opt 1) 152.0000000 cr2 Chronic AWQC Exeedences ( Smpl/ Opt 2) . acb Acute AWQC Exeedences ( Smpl/ BL) 90.0000000 ac1 Acute AWQC Exeedences ( Smpl/ Opt 1) 34.0000000 ac2 Acute AWQC Exeedences ( Smpl/ Opt 2) . ntlbw Chronic AWQC Exeedences ( Ntl/ BL) 2582.49 ntl1w Chronic AWQC Exeedences ( Ntl/ OPt 1) 1368.52 ntl2w Chronic AWQC Exeedences ( Ntl/ OPt 2) . ntlbo Acute AWQC Exeedences ( Ntl/ Bl) 630.5351378 ntl1o Acute AWQC Exeedences ( Ntl/ Opt 1) 254.3570580 ntl2o Acute AWQC Exeedences ( Ntl/ Opt 2) . chemcrb Number of Chronic pol. ( BL) 47.0000000 chemcr1 Number of Chronic pol. ( Opt 1) 39.0000000 chemcr2 Number of Chronic pol. ( Opt 2) 0 chemacb Number of Acute pol. ( BL) 17.0000000 chemac1 Number of Acute pol. ( Opt 1) 12.0000000 chemac2 Number of Acute pol. ( Opt 2) 0 212                                                                  MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 14 Output from Aquatic Life ( Sub 3) 00: 44 Sunday, February 9, 2003 IL Reaches Exceeding AWQC Limits for Aquatic Life ( National)

The MEANS Procedure Variable Label Sum                                                                         chrnc_ b # of Reaches Exceeding Chronic AWQC ( Bl) 927.6219964 chrnc_ 1 # of Reaches Exceeding Chronic AWQC ( Opt 1) 539.0794437 chrnc_ 2 # of Reaches Exceeding Chronic AWQC ( Opt 2) 0 acutr_ b # of Reaches Exceeding Acute AWQC ( Bl) 329.9659873 acutr_ 1 # of Reaches Exceeding Acute AWQC ( Opt 1) 86.2192946 acutr_ 2 # of Reaches Exceeding Acute AWQC ( Opt 2) 0 multi_ b Either ( BL) 927.6219964 multi_ 1 Either ( Opt 1) 539.0794437 multi_ 2 Either ( Opt 2) 0                                                                         MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 16 Output from Aquatic Life ( Sub 3) 00: 44 Sunday, February 9, 2003 AWQC Eliminations and Reductions ( National) Number of National Reaches The MEANS Procedure Variable Label Sum                                                           elim1_ 1 Elim. Chronic ( Opt1) 388.5425526 elim2_ 1 Elim. Chronic ( Opt2) 948.6368971 elim1_ 2 Elim. Acute ( Opt1) 243.7466927 elim2_ 2 Elim. Acute ( Opt2) 342.1990236 red1_ 1 Redu. Chronic ( Opt1) 181.9521911 red2_ 1 Redu. Chronic ( Opt2) 0 red1_ 2 Redu. Acute ( Opt1) 43.6336308 red2_ 2 Redu. Acute ( Opt2) 0 reduc1 Total Reaches Reducing ( Opt1) 181.9521911 reduc2 Total Reaches Reducing ( Opt2) 0                                                           MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 18 Output from Aquatic Life ( Sub 3) 00: 44 Sunday, February 9, 2003 Exceedences by Chemical C H c c c c c c E h h h h h h M n n n n n n e e e e e e I t t t t t t m m m m m m O C c c c a a a l l l l l l c c c a a a b A r r r c c c b 1 2 b 1 2 r r r c c c s L b 1 2 b 1 2 w w w o o o b 1 2 b 1 2 1 Acenaphthene 1 1 . . . . 8.782 8.782 . . . . 1 1 0 0 0 0 213 2 Acrolein 9 6 . 6 4 . 43.572 33.283 . 33.283 26.1361 . 1 1 0 1 1 0 3 Aluminum 5 2 . 1 . . 31.794 12.233 . 10.177 . . 1 1 0 1 0 0 4 Ammonia as N 6 . . . . . 51.244 . . . . . 1 0 0 0 0 0 5 Aniline 9 8 . 1 . . 44.993 42.065 . 8.782 . . 1 1 0 1 0 0 6 Anthracene 10 5 . 10 5 . 63.916 29.218 . 63.916 29.2184 . 1 1 0 1 1 0 7 Biphenyl 1 1 . 1 . . 8.782 8.782 . 8.782 . . 1 1 0 1 0 0 8 Butyl benzyl phthalate 1 . . . . . 8.782 . . . . . 1 0 0 0 0 0 9 Cadmium 7 4 . 2 1 . 70.177 21.015 . 9.059 6.1309 . 1 1 0 1 1 0 10 Carbon disulfide 7 6 . . . . 37.771 34.428 . . . . 1 1 0 0 0 0 11 Chromium 7 2 . 2 . . 46.318 12.233 . 6.871 . . 1 1 0 1 0 0 12 Cobalt 2 2 . . . . 12.233 12.233 . . . . 1 1 0 0 0 0 13 Copper 40 8 . 34 8 . 343.790 69.234 . 241.245 69.2340 . 1 1 0 1 1 0 14 Cyanide 1 . . . . . 3.344 . . . . . 1 0 0 0 0 0 15 Di n butyl phthalate 1 1 . . . . 8.782 8.782 . . . . 1 1 0 0 0 0 16 Di n octyl phthalate 2 2 . . . . 12.233 12.233 . . . . 1 1 0 0 0 0 17 Dibenzofuran 3 2 . . . . 21.015 12.233 . . . . 1 1 0 0 0 0 18 Dibenzothiophene 3 2 . . . . 15.162 12.233 . . . . 1 1 0 0 0 0 19 Dimethylphenanthrene, 3,6 4 3 . . . . 23.943 21.015 . . . . 1 1 0 0 0 0 20 Dinitrophenol, 2,4 1 1 . 1 1 . 8.782 8.782 . 8.782 8.7819 . 1 1 0 1 1 0 21 Dinitrotoluene, 2,6 3 3 . . . . 21.015 21.015 . . . . 1 1 0 0 0 0 22 Diphenyl Ether 3 3 . . . . 21.015 21.015 . . . . 1 1 0 0 0 0 23 Fluoranthene 6 4 . 3 3 . 30.369 24.097 . 21.015 21.0149 . 1 1 0 1 1 0 24 Fluorene 5 3 . 3 3 . 27.026 21.015 . 21.015 21.0149 . 1 1 0 1 1 0 25 Fluoride 4 2 . . . . 53.870 12.846 . . . . 1 1 0 0 0 0 26 Iron 2 . . . . . 11.710 . . . . . 1 0 0 0 0 0 27 Isopropylnaphthalene, 2 3 2 . . . . 15.162 12.233 . . . . 1 1 0 0 0 0 28 Lead 28 11 . 2 1 . 244.307 82.780 . 9.059 6.1309 . 1 1 0 1 1 0 29 Magnesium 2 2 . . . . 12.233 12.233 . . . . 1 1 0 0 0 0 30 Manganese 3 . . . . . 31.707 . . . . . 1 0 0 0 0 0 31 Methylfluorene, 1 3 2 . . . . 15.162 12.233 . . . . 1 1 0 0 0 0 32 Methylnaphthalene, 2 2 2 . . . . 12.233 12.233 . . . . 1 1 0 0 0 0 33 Methylphenanthrene, 1 3 2 . . . . 15.162 12.233 . . . . 1 1 0 0 0 0 34 Molybdenum 12 3 . . . . 103.358 38.506 . . . . 1 1 0 0 0 0 35 Naphthalene 1 1 . . . . 8.782 8.782 . . . . 1 1 0 0 0 0 36 Nickel 16 3 . 5 . . 163.285 15.577 . 22.558 . . 1 1 0 1 0 0 37 Phenanthrene 4 3 . 1 1 . 23.943 21.015 . 8.782 8.7819 . 1 1 0 1 1 0 38 Phenol 1 . . . . . 8.782 . . . . . 1 0 0 0 0 0 39 Pyrene 3 3 . . . . 21.015 21.015 . . . . 1 1 0 0 0 0 40 Selenium 6 3 . 2 2 . 77.618 50.329 . 12.233 12.2330 . 1 1 0 1 1 0 41 Silver 20 15 . 6 2 . 166.324 131.148 . 60.122 12.2330 . 1 1 0 1 1 0 42 Styrene 1 . . . . . 8.782 . . . . . 1 0 0 0 0 0 43 Sulfide 12 9 . . . . 293.160 283.238 . . . . 1 1 0 0 0 0 44 Tin 17 4 . . . . 83.058 21.015 . . . . 1 1 0 0 0 0 45 Titanium 1 . . . . . 6.456 . . . . . 1 0 0 0 0 0 46 Vanadium 16 13 . . . . 156.682 141.727 . . . . 1 1 0 0 0 0 47 Zinc 10 3 . 10 3 . 84.855 33.447 . 84.855 33.4473 . 1 1 0 1 1 0 214 " Fishing Recreation ( Sub 7). sas" Data on human health based and aquatic life based AWQC exceedences for each RF1 reach are merged, and the total number of AWQC exceedences is generated. Based on the reduction in AWQC exceedences, the baseline value of a fishing day, and a range of percentage increases in value from AWQC reductions, benefits to recreational fishermen are estimated. 215 /* Program: Fishing Recreation ( Sub 7). sas Authors: Abt Associates Inc; Rob Sartain, Elena Besedin, Stuart Smith Date: 04/ 21/ 2000 Updated: 09/ 26/ 2002 This program estimates recreation benefits to sport anglers stemming from improvments in or the elimination of awqc exceedences for human health or aquatic life on reaches. */

title2 ' Fishing Recreation ( Sub 7)';

/* COUNTY LEVEL FISHING LICENSE DATA */

/* Load fishing license data here. Sum over the types by county. Sum over the counties containing a given reach.

Final data set contains reaches with the total possible number of affected anglers. This number will be multiplied by the ratio of reach length to total reach miles in affected counties to obtain an estimate of the number of affected anglers. */

data licens; set mpm2. fish_ lic;

fw_ lic = sum( fw_ res, fw_ reso, fw_ sres, fw_ smres, fw_ dres, trout, fw_ life, fw_ hunt1, fw_ 2hunt, fw_ hunts, fw_ hunt2, fw_ salt1, fw_ trout); sw_ lic = sum( fw_ salt1, sw_ res, sw_ dres);

if fw_ lic=. then fw_ lic= 0; if sw_ lic=. then sw_ lic= 0; cntyfip = fips;

keep cntyfip fw_ lic sw_ lic; run;

proc sort nodupkey data= licens; by cntyfip; quit;

/* DETERMINE COUNTIES AFFECTED BY REACH */ data facs; set mpm2. fac_ info; run;

proc sort data= facs; by siteid; quit;

data cnty; set facs;

cntyfip1= stco1; cntyfip2= stco2; cntyfip3= stco3; cntyfip4= stco4; cntyfip5= stco5; cntyfip6= stco6; cntyfip7= stco7; cntyfip8= stco8; 216 keep reachnmb cntyfip1 cntyfip8; run;

proc sort nodupkey data= cnty; by reachnmb; quit;

proc transpose data= cnty out= cnty2; by reachnmb; var cntyfip1 cntyfip8; quit;

data cnty3; set cnty2; cntyfip = col1/ 1; if cntyfip=. then delete; keep reachnmb cntyfip; run;

proc sort data= cnty3; by cntyfip; quit;

/* OBTAIN REACH MILES BY COUNTY */

data rchmile; set mpm2. us_ rchmi; cntyfip= fips; keep countymi cntyfip; run;

proc sort data= rchmile; by cntyfip;

/* PULL TOGETHER OUR REACHES WITH ASSOCIATED COUNTIES */

data tot_ mile; merge cnty3 rchmile; by cntyfip;

if reachnmb=. then delete; run;

proc sort data= tot_ mile; by reachnmb cntyfip; quit; proc sort nodupkey data= facs; by reachnmb; quit;

/* SUM REACH MILES OVER COUNTIES CONTAINING THE SAME REACH FINISHED DATA SET CONTAINS TOTAL NUMBER OF REACH MILES IN COUNTIES CONTAINING A GIVEN REACH. BRING IN THE FISHING LICNESE DATA BY COUNTY. */

proc sql;

create table miles as select tot_ mile. reachnmb, sum( countymi) as denom, min( segl) as leng, sum( fw_ lic) as fwlictot, sum( sw_ lic) as swlictot from tot_ mile, facs, licens where tot_ mile. reachnmb = facs. reachnmb and tot_ mile. cntyfip = licens. cntyfip 217 group by tot_ mile. reachnmb order by tot_ mile. reachnmb; quit;

proc sort data= miles; by reachnmb; quit;

/* ESTIMATE SPORT ANGLERS */

data pops; set miles; by reachnmb; if first. reachnmb;

/* no family adjustment, just anglers */ /* no subsistence anglers */

if leng=. then fwspt = fwlictot* 0.06; else fwspt = fwlictot*( leng/ denom);

if leng=. then swspt = swlictot* 0.06; else swspt = swlictot*( leng/ denom);

fwsptpop = round( fwspt); swsptpop = round( swspt);

sportpop = fwsptpop + swsptpop;

keep reachnmb sportpop fwsptpop swsptpop fwlictot swlictot; run;

proc sort data= pops; by reachnmb; quit;

/* Combine Population with average fishing days */

proc sql;

create table pops2 as select pops.*, fip from pops, mpm2. rchstate where pops. reachnmb = rchstate. reachnmb order by reachnmb;

create table pops3 as select pops2.*, avgfw, avgsw from pops2, mpm2. avgdfish where pops2. fip = avgdfish. fip order by reachnmb; quit;

/* Pull in the reach data sets with exceedence events and combine by reach, casno, event weight */ 218 proc sort data= rundat. rec_ hh; by reachnmb casno event_ wt; quit; proc sort data= rundat. rec_ al; by reachnmb casno event_ wt; quit;

data rchexc; merge rundat. rec_ hh rundat. rec_ al; by reachnmb casno event_ wt;

/* Calculate total exceedences for a reach, both sample and national Sum all exceedences for a reach, that's acute and chronic for AL and h2o/ org and org for HH */

if orgh2o_ b=. then orgh2o_ b= 0; if org_ b=. then org_ b= 0; if chronicb=. then chronicb= 0; if acuteb=. then acuteb= 0; if orgh2ot1=. then orgh2ot1= 0; if org_ t1=. then org_ t1= 0; if chronic1=. then chronic1= 0; if acute1=. then acute1= 0; if orgh2ot2=. then orgh2ot2= 0; if org_ t2=. then org_ t2= 0; if chronic2=. then chronic2= 0; if acute2=. then acute2= 0;

base = orgh2o_ b + org_ b + chronicb + acuteb; trt1 = orgh2ot1 + org_ t1 + chronic1 + acute1; trt2 = orgh2ot2 + org_ t2 + chronic2 + acute2; run;

proc sql; create table evntsum as select reachnmb, event_ wt, sum( base) as exc_ b, sum( trt1) as exc_ 1, sum( trt2) as exc_ 2 from rchexc group by reachnmb, event_ wt order by reachnmb, event_ wt; quit;

proc print data= evntsum; title3 ' Summary of Exceedences by Reach'; quit;

/* Calculate benefit ratios Ratio = # exceedences eliminated / # baseline exceedences If all exceedences eliminated then Ratio = 1 and full benefit is realized. Otherwise the fraction of the benefit realized is equal to the fraction of exceedences eliminated. */

data benrat; set evntsum; 219 bfrac_ 1 = ( exc_ b  exc_ 1)/ exc_ b; bfrac_ 2 = ( exc_ b  exc_ 2)/ exc_ b;

weight = event_ wt;

if weight=. then weight = 1;

keep reachnmb weight bfrac_ 1 bfrac_ 2 exc_ b exc_ 1 exc_ 2; run;

/* Pull together the benefit fractions and fishing days and populations */

proc sql;

create table benefit as select benrat.*, sportpop, fwsptpop, swsptpop, avgfw, avgsw from benrat, pops3 where benrat. reachnmb = pops3. reachnmb order by benrat. reachnmb; quit;

data benefit2; set benefit;

/* Base values for Fishing Days Values in $ 2001 Low ($ 99) 28.11 Upper ($ 99) 60.43 */

bvalfwl = fwsptpop * avgfw * 28.11; bvalfwh = fwsptpop * avgfw * 60.43; bvalswl = swsptpop * avgsw * 28.11; bvalswh = swsptpop * avgsw * 60.43; midfw = fwsptpop * avgfw * 42.12; midsw = swsptpop * avgsw * 42.12;

tripsfw = fwsptpop* avgfw; tripssw = swsptpop* avgsw;

/* Benefits under Average Values Increase range of 9.83% to 14.71% Midpoint 12.27% */

avg1_ lf = bvalfwl * bfrac_ 1 * 0.0983; avg1_ hf = bvalfwh * bfrac_ 1 * 0.1471; avg1_ ls = bvalswl * bfrac_ 1 * 0.0983; avg1_ hs = bvalswh * bfrac_ 1 * 0.1471; mid1_ f = midfw * bfrac_ 1 * 0.1227; mid1_ s = midsw * bfrac_ 1 * 0.1227; 220 avg2_ lf = bvalfwl * bfrac_ 2 * 0.0983; avg2_ hf = bvalfwh * bfrac_ 2 * 0.1471; avg2_ ls = bvalswl * bfrac_ 2 * 0.0983; avg2_ hs = bvalswh * bfrac_ 2 * 0.1471; mid2_ f = midfw * bfrac_ 2 * 0.1227; mid2_ s = midsw * bfrac_ 2 * 0.1227; run;

proc means sum; var avg1_ lf avg1_ hf avg1_ ls avg1_ hs avg2_ lf avg2_ hf avg2_ ls avg2_ hs mid1_ f mid1_ s mid2_ f mid2_ s fwsptpop swsptpop tripsfw tripssw bvalfwl bvalfwh bvalswl bvalswh; title3 ' Sample Benefits ( includes partial)'; quit;

proc means sum; weight weight; var avg1_ lf avg1_ hf avg1_ ls avg1_ hs avg2_ lf avg2_ hf avg2_ ls avg2_ hs mid1_ f mid1_ s mid2_ f mid2_ s fwsptpop swsptpop tripsfw tripssw bvalfwl bvalfwh bvalswl bvalswh; output out = rundat. rec_ fsh( drop = _ type_ _ freq_) sum = avg1_ lf avg1_ hf avg1_ ls avg1_ hs avg2_ lf avg2_ hf avg2_ ls avg2_ hs mid1_ f mid1_ s mid2_ f mid2_ s fwsptpop swsptpop tripsfw tripssw bvalfwl bvalfwh bvalswl bvalswh; title3 ' National Benefits ( includes partial)'; quit;

data rundat. rec_ fish( keep = activity benl1 benm1 benh1 benl2 benm2 benh2 pop); set rundat. rec_ fsh; activity = " Fishing"; pop = fwsptpop + swsptpop; benl1 = avg1_ lf + avg1_ ls; benm1 = mid1_ f + mid1_ s; benh1 = avg1_ hf + avg1_ hs; benl2 = avg2_ lf + avg2_ ls; benm2 = mid2_ f + mid2_ s; benh2 = avg2_ hf + avg2_ hs; run;

data reaches; set benrat;

if exc_ b > 0 then baseexc = 1; else baseexc = 0; if exc_ 1 > 0 then post1 = 1; else post1= 0; if exc_ 2 > 0 then post2 = 1; else post2= 0; if bfrac_ 1 = 1 then elm_ 1 = 1; else elm_ 1 = 0; if bfrac_ 1 < 1 and bfrac_ 1 > 0 then red_ 1 = 1; else red_ 1 = 0; if bfrac_ 2 = 1 then elm_ 2 = 1; else elm_ 2 = 0; if bfrac_ 2 < 1 and bfrac_ 2 > 0 then red_ 2 = 1; else red_ 2 = 0; run;

proc means sum data= reaches; var baseexc post1 post2 elm_ 1 elm_ 2 red_ 1 red_ 2; label 221 baseexc = ' Baseline Reaches Exceeding AWQCs' post1 = ' Post Reaches Exceeding AWQCs ( Opt1)' post2 = ' Post Reaches Exceeding AWQCs ( Opt2)' elm_ 1 = ' Reaches Eliminating ( Opt1)' elm_ 2 = ' Reaches Eliminating ( Opt2)' red_ 1 = ' Reaches Reducing ( Opt1)' red_ 2 = ' Reaches Reducing ( Opt2)'; title3 ' Reach Information'; title4 ' Sample Level'; quit;

proc means sum data= reaches; weight weight; var baseexc post1 post2 elm_ 1 elm_ 2 red_ 1 red_ 2; label baseexc = ' Baseline Reaches Exceeding AWQCs' post1 = ' Post Reaches Exceeding AWQCs ( Opt1)' post2 = ' Post Reaches Exceeding AWQCs ( Opt2)' elm_ 1 = ' Reaches Eliminating ( Opt1)' elm_ 2 = ' Reaches Eliminating ( Opt2)' red_ 1 = ' Reaches Reducing ( Opt1)' red_ 2 = ' Reaches Reducing ( Opt2)'; title3 ' Reach Information'; title4 ' National Level'; quit;

data check; set benefit2;

if bfrac_ 1 <= 0 then delete; run;

data check2; merge check pops; by reachnmb; if weight=. then delete;

keep reachnmb fwlictot swlictot fwsptpop swsptpop weight; run;

/* End Fishing Recreation ( Sub 7). sas */ 222 Program Output 223 MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 62 Fishing Recreation ( Sub 7) 00: 44 Sunday, February 9, 2003 National Benefits ( includes partial)

The MEANS Procedure Variable Sum                          avg1_ lf 53258143.71 avg1_ hf 171331386 avg1_ ls 638806.54 avg1_ hs 2055039.90 avg2_ lf 583300985 avg2_ hf 1876478588 avg2_ ls 14341905.01 avg2_ hs 46137891.70 mid1_ f 99610377.74 mid1_ s 1194779.93 mid2_ f 1090966139 mid2_ s 26824115.08 fwsptpop 12863228.23 swsptpop 543090.21 tripsfw 211095194 tripssw 5190300.21 bvalfwl 5933885906 bvalfwh 12756482580 bvalswl 145899339 bvalswh 313649842                          MP& M National Benefits Preferred( Opt1= Pref, Opt2= NULL) 64 Fishing Recreation ( Sub 7) 00: 44 Sunday, February 9, 2003 Reach Information National Level The MEANS Procedure Variable Label Sum                                                                 baseexc Baseline Reaches Exceeding AWQCs 5998.58 post1 Post Reaches Exceeding AWQCs ( Opt1) 5844.99 post2 Post Reaches Exceeding AWQCs ( Opt2) 0 elm_ 1 Reaches Eliminating ( Opt1) 153.5907860 elm_ 2 Reaches Eliminating ( Opt2) 5998.58 red_ 1 Reaches Reducing ( Opt1) 581.2015684 red_ 2 Reaches Reducing ( Opt2) 0                                                                 224 Near water Activity and Boating 225 " Human Health AWQC Exceedences ( Sub 2). sas" The number of human health based AWQC exceedences for each RF1 reach are estimated and output to a data set. See Human Health Benefits Cancer from Fish Consumption for program code and output. 226 " Aquatic Life AWQC Exceedences ( Sub 3). sas" The number of aquatic life based AWQC exceedences for each RF1 reach are estimated and output to a data set. See Recreational Fishing for program code and output. 227 " Other Recreation ( Sub 8). sas" Data on human health based and aquatic life based AWQC exceedences for each RF1 reach are merged, and the total number of AWQC exceedences is generated. Based on the reduction in AWQC exceedences, the baseline value of a boating or near water activity day, and a range of percentage increases in value from AWQC reductions, benefits to boaters and near water activity participants are estimated. 228 /* Program: Other Recreation ( Sub 8). sas Authors: Abt Associates Inc; Rob Sartain, Elena Besedin, Stuart Smith Date: 04/ 21/ 2000 Updated: 09/ 26/ 2002 This program estimates recreation benefits to activity participants stemming from improvements in or the elimination of awqc exceedences for human health or aquatic life on reaches. Valuations for viewing and boating are estimated. */

% let tbl_ path = % sysfunc( pathname( rundat));

title2 ' Other Recreation ( Sub 8)';

/* DETERMINE COUNTIES AFFECTED BY REACH */ data facs; set mpm2. fac_ info; run;

proc sort data= facs; by siteid; quit;

data cnty; set facs;

cntyfip1= stco1; cntyfip2= stco2; cntyfip3= stco3; cntyfip4= stco4; cntyfip5= stco5; cntyfip6= stco6; cntyfip7= stco7; cntyfip8= stco8;

keep reachnmb cntyfip1 cntyfip8; run;

proc sort nodupkey data= cnty; by reachnmb; quit;

proc transpose data= cnty out= cnty2; by reachnmb; var cntyfip1 cntyfip8; quit;

data cnty3; set cnty2; cntyfip = col1/ 1; if cntyfip=. then delete; keep reachnmb cntyfip; run;

proc sort data= cnty3; by cntyfip; quit;

/* OBTAIN REACH MILES BY COUNTY */

data rchmile; set mpm2. us_ rchmi; cntyfip= fips; keep countymi cntyfip; 229 run;

proc sort data= rchmile; by cntyfip;

/* PULL TOGETHER OUR REACHES WITH ASSOCIATED COUNTIES */

data tot_ mile; merge cnty3 rchmile; by cntyfip;

if reachnmb=. then delete; run;

proc sort data= tot_ mile; by reachnmb cntyfip; quit; proc sort nodupkey data= facs; by reachnmb; quit;

/* SUM REACH MILES OVER COUNTIES CONTAINING THE SAME REACH FINISHED DATA SET CONTAINS TOTAL NUMBER OF REACH MILES IN COUNTIES CONTAINING A GIVEN REACH. BRING IN THE FISHING LICNESE DATA BY COUNTY. */

proc sql; create table miles as select tot_ mile. reachnmb, sum( countymi) as denom, min( segl) as leng from tot_ mile, facs where tot_ mile. reachnmb = facs. reachnmb group by tot_ mile. reachnmb order by tot_ mile. reachnmb; quit;

/* Merge in the county populations and state abbreviations */ proc sql; create table rchpop as select cnty3. reachnmb, cnty3. cntyfip, pop2000 as totpop, stateabv as state, denom, leng from cnty3, mpm2. cntypop, mpm2. rchstate, miles where cnty3. cntyfip = cntypop. fips and cnty3. reachnmb = rchstate. reachnmb and cnty3. reachnmb = miles. reachnmb order by reachnmb;

/* ESTIMATE Number of Participants by Activity */

create table pops as select reachnmb, sum( totpop) as pop, min( state) as stateabv, avg( denom) as den, avg( leng) as lengt from rchpop group by reachnmb order by reachnmb;

create table pops2 as select reachnmb, stateabv, viewtrp, boattrp, swimtrp, pop, pop* pctview as totview, pop* pctboat as totboat, pop* pctswim as totswim, case 230 when lengt is missing then pop* pctview* 0.06 else pop* pctview*( lengt/ den) end as viewers, case when lengt is missing then pop* pctboat* 0.06 else pop* pctboat*( lengt/ den) end as boaters, case when lengt is missing then pop* pctswim* 0.06 else pop* pctswim*( lengt/ den) end as swimers from pops, mpm2. recpcts where pops. stateabv = recpcts. state order by reachnmb; quit;

/* Pull in the reach data sets with exceedence events and combine by reach, casno, event weight */

proc sort data= rundat. rec_ hh; by reachnmb casno event_ wt; quit; proc sort data= rundat. rec_ al; by reachnmb casno event_ wt; quit;

data rchexc; merge rundat. rec_ hh rundat. rec_ al; by reachnmb casno event_ wt;

/* Calculate total exceedences for a reach, both sample and national Sum all exceedences for a reach, that's acute and chronic for AL and h2o/ org and org for HH */

if orgh2o_ b=. then orgh2o_ b= 0; if org_ b=. then org_ b= 0; if chronicb=. then chronicb= 0; if acuteb=. then acuteb= 0; if orgh2ot1=. then orgh2ot1= 0; if org_ t1=. then org_ t1= 0; if chronic1=. then chronic1= 0; if acute1=. then acute1= 0; if orgh2ot2=. then orgh2ot2= 0; if org_ t2=. then org_ t2= 0; if chronic2=. then chronic2= 0; if acute2=. then acute2= 0;

base = orgh2o_ b + org_ b + chronicb + acuteb; trt1 = orgh2ot1 + org_ t1 + chronic1 + acute1; trt2 = orgh2ot2 + org_ t2 + chronic2 + acute2; run;

proc sql; create table evntsum as select reachnmb, event_ wt, sum( base) as exc_ b, sum( trt1) as exc_ 1, sum( trt2) as exc_ 2 from rchexc 231 group by reachnmb, event_ wt order by reachnmb, event_ wt; quit;

proc print data= evntsum; quit;

/* Calculate benefit ratios Ratio = # exceedences eliminated / # baseline exceedences If all exceedences eliminated then Ratio = 1 and full benefit is realized. Otherwise the fraction of the benefit realized is equal to the fraction of exceedences eliminated. */

data benrat; set evntsum;

bfrac_ 1 = ( exc_ b  exc_ 1)/ exc_ b; bfrac_ 2 = ( exc_ b  exc_ 2)/ exc_ b;

weight = event_ wt;

if weight=. then weight = 1;

keep reachnmb weight bfrac_ 1 bfrac_ 2 exc_ b exc_ 1 exc_ 2; run;

/* Pull together the benefit fractions and fishing days and populations */

proc sql;

create table benefit as select benrat.*, viewers, boaters, viewtrp, boattrp, swimers, swimtrp from benrat, pops2 where benrat. reachnmb = pops2. reachnmb order by benrat. reachnmb; quit;

data benefit2; set benefit;

/* Base values for Viewing and Boating and Swimming Days Viewing Day ($ 2001) 22.73  28.73 Boating Day ($ 2001) 37.30  59.26 midpoint view 26.28 midpoint boat 48.30 */

bvlow = viewers * 22.73 * viewtrp; bvhigh = viewers * 28.73 * viewtrp; bblow = boaters * 37.30 * boattrp; bbhigh = boaters * 59.26 * boattrp; 232 midview = viewers * 26.28 * viewtrp; midboat = boaters * 48.30 * boattrp;

tripsv = viewers* viewtrp; tripsb = boaters* boattrp;

/* Benefits from Average values Viewing Increase of 14.94% to 21.34% Boating Increase of 6.51% to 11.35%

mid viewing 18.14% mid boating 8.93% */

av1_ low = bvlow * bfrac_ 1 * 0.1494; av1_ high = bvhigh * bfrac_ 1 * 0.2134; av2_ low = bvlow * bfrac_ 2 * 0.1494; av2_ high = bvhigh * bfrac_ 2 * 0.2134;

ab1_ low = bblow * bfrac_ 1 * 0.0651; ab1_ high = bbhigh * bfrac_ 1 * 0.1135; ab2_ low = bblow * bfrac_ 2 * 0.0651; ab2_ high = bbhigh * bfrac_ 2 * 0.1135;

midvw1 = midview * bfrac_ 1 * 0.1814; midbt1 = midboat * bfrac_ 1 * 0.0893; midvw2 = midview * bfrac_ 2 * 0.1814; midbt2 = midboat * bfrac_ 2 * 0.0893; run;

proc means sum; var av1_ low av1_ high av2_ low av2_ high ab1_ low ab1_ high ab2_ low ab2_ high midvw1 midvw2 midbt1 midbt2 viewers boaters tripsv tripsb; label av1_ low = ' Avg View Benefits  Low ( Opt 1)' av1_ high = ' Avg View Benefits  High( Opt 1)' av2_ low = ' Avg View Benefits  Low ( Opt 2)' av2_ high = ' Avg View Benefits  High( Opt 2)' ab1_ low = ' Avg Boat Benefits  Low ( Opt 1)' ab1_ high = ' Avg Boat Benefits  High( Opt 1)' ab2_ low = ' Avg Boat Benefits  Low ( Opt 2)' ab2_ high = ' Avg Boat Benefits  High( Opt 2)' midvw1 = ' Mid View Ben ( Opt1)' midvw2 = ' Mid View Ben ( Opt2)' midbt1 = ' Mid Boat Ben ( Opt1)' midbt2 = ' Mid Boat Ben ( Opt2)' viewers = ' Number of Viewing Participants' boaters = ' Number of Boating Participants' tripsv = ' Number of Viewing Days' tripsb = ' Number of Boating Days' title3 ' Sample Benefits ( includes partial)'; title4 ' Viewing, Boating Activities '; 233 quit;

proc means sum; weight weight; var av1_ low av1_ high av2_ low av2_ high ab1_ low ab1_ high ab2_ low ab2_ high midvw1 midvw2 midbt1 midbt2 viewers boaters tripsv tripsb; output out = rundat. rec_ oth( drop = _ type_ _ freq_) sum = av1_ low av1_ high av2_ low av2_ high ab1_ low ab1_ high ab2_ low ab2_ high midvw1 midvw2 midbt1 midbt2 viewers boaters tripsv tripsb; label av1_ low = ' Avg View Benefits  Low ( Opt 1)' av1_ high = ' Avg View Benefits  High( Opt 1)' av2_ low = ' Avg View Benefits  Low ( Opt 2)' av2_ high = ' Avg View Benefits  High( Opt 2)' ab1_ low = ' Avg Boat Benefits  Low ( Opt 1)' ab1_ high = ' Avg Boat Benefits  High( Opt 1)' ab2_ low = ' Avg Boat Benefits  Low ( Opt 2)' ab2_ high = ' Avg Boat Benefits  High( Opt 2)' midvw1 = ' Mid View Ben ( Opt1)' midvw2 = ' Mid View Ben ( Opt2)' midbt1 = ' Mid Boat Ben ( Opt1)' midbt2 = ' Mid Boat Ben ( Opt2)' viewers = ' Number of Viewing Participants' boaters = ' Number of Boating Participants' tripsv = ' Number of Viewing Days' tripsb = ' Number of Boating Days' title3 ' National Benefits ( includes partial)'; title4 ' Viewing, Boating Activities '; quit;

data rundat. rec_ othr( keep = activity benl1 benm1 benh1 benl2 benm2 benh2 pop); set rundat. rec_ oth;

if _ n_ = 1 then do; activity = " Boating"; benl1 = ab1_ low; benm1 = midbt1; benh1 = ab1_ high; benl2 = ab2_ low; benm2 = midbt2; benh2 = ab2_ high; pop = boaters; output; activity = " Viewing"; benl1 = av1_ low; benm1 = midvw1; benh1 = av1_ high; benl2 = av2_ low; benm2 = midvw2; benh2 = av2_ high; pop = viewers; output; 234 end; run;

data allrec; set rundat. rec_ fish rundat. rec_ othr; run;

proc means sum data= allrec; var benl1 benm1 benh1 benl2 benm2 benh2; output out = rectot( drop = _ type_ _ freq_) sum = benl1 benm1 benh1 benl2 benm2 benh2; quit;

data totls( keep = activity benl1 benm1 benh1 benl2 benm2 benh2); set rectot; if _ n_ = 1 then do; activity = " Recreational Total"; output; activity = " Non Use Benefits"; benl1 = benl1 * 0.50; benm1 = benm1 * 0.50; benh1 = benh1 * 0.50; benl2 = benl2 * 0.50; benm2 = benm2 * 0.50; benh2 = benh2 * 0.50; output; activity = " Total Rec and Nonuse"; benl1 = benl1 + ( benl1 * 0.50); benm1 = benm1 + ( benm1 * 0.50); benh1 = benh1 + ( benh1 * 0.50); benl2 = benl2 + ( benl2 * 0.50); benm2 = benm2 + ( benm2 * 0.50); benh2 = benh2 + ( benh2 * 0.50); end; run;

data exlrec ( keep= titlecol col1 col8 label=" Recreational and Nonuse Benefits  Weighted"); set allrec totls; length titlecol $ 50 col1  col8 $ 128;

if _ n_ = 1 then do; titlecol = " Table: Recreational and Non Use Benefits, Weighted Estimates"; output; titlecol = " Run on & sysdate"; output; titlecol = " "; output; col2 = " Option 1"; col5 = " Option 2"; output; col1 = " Activity"; col2 = " Low"; col3 = " Mid"; col4 = " High"; col5 = " Low"; col6 = " Mid"; col7 = " High"; 235 col8 = " Population"; output; end;

col1 = trim( activity); col2 = trim( benl1); col3 = trim( benm1); col4 = trim( benh1); col5 = trim( benl2); col6 = trim( benm2); col7 = trim( benh2); col8 = pop; output; run;

proc export data = exlrec outfile = " C:\ Output Files\ Recreational and Nonuse Benefits. xls" dbms= excel2000 replace; run;

data reaches; set benrat;

if exc_ b > 0 then baseexc = 1; else baseexc = 0; if exc_ 1 > 0 then post1 = 1; else post1= 0; if exc_ 2 > 0 then post2 = 1; else post2= 0; if bfrac_ 1 = 1 then elm_ 1 = 1; else elm_ 1 = 0; if bfrac_ 1 < 1 and bfrac_ 1 > 0 then red_ 1 = 1; else red_ 1 = 0; if bfrac_ 2 = 1 then elm_ 2 = 1; else elm_ 2 = 0; if bfrac_ 2 < 1 and bfrac_ 2 > 0 then red_ 2 = 1; else red_ 2 = 0; run;

proc means sum data= reaches; var baseexc post1 post2 elm_ 1 elm_ 2 red_ 1 red_ 2; label baseexc = ' Baseline Reaches Exceeding AWQCs' post1 = ' Post Reaches Exceeding AWQCs ( Opt1)' post2 = ' Post Reaches Exceeding AWQCs ( Opt2)' elm_ 1 = ' Reaches Eliminating ( Opt1)' elm_ 2 = ' Reaches Eliminating ( Opt2)' red_ 1 = ' Reaches Reducing ( Opt1)' red_ 2 = ' Reaches Reducing ( Opt2)'; title3 ' Reach Information'; title4 ' Sample Level'; quit;

proc means sum data= reaches; weight weight; var baseexc post1 post2 elm_ 1 elm_ 2 red_ 1 red_ 2; label baseexc = ' Baseline Reaches Exceeding AWQCs' post1 = ' Post Reaches Exceeding AWQCs ( Opt1)' post2 = ' Post Reaches Exceeding AWQCs ( Opt2)' elm_ 1 = ' Reaches Eliminating ( Opt1)' elm_ 2 = ' Reaches Eliminating ( Opt2)' 236 red_ 1 = ' Reaches Reducing ( Opt1)' red_ 2 = ' Reaches Reducing ( Opt2)'; title3 ' Reach Information'; title4 ' National Level'; quit;

/* Determine the number of counties affected by the benefiting reaches */

data goodrch1; set benefit2;

if bfrac_ 1 <= 0 then delete; keep reachnmb weight; run;

data goodrch2; set benefit2;

if bfrac_ 2 <= 0 then delete;

keep reachnmb weight; run;

proc sql;

create table cntyb1 as select cnty3. reachnmb, cntyfip, weight from cnty3, goodrch1 where cnty3. reachnmb = goodrch1. reachnmb;

create table cntyb2 as select cnty3. reachnmb, cntyfip, weight from cnty3, goodrch2 where cnty3. reachnmb = goodrch2. reachnmb;

create table cntyb1a as select count( distinct cntyfip) as numcnty1 from cntyb1;

create table cntyb2a as select count( distinct cntyfip) as numcnty2 from cntyb2;

create table cntyb1b as select distinct cntyfip, avg( weight) as wght from cntyb1 group by cntyfip;

create table cntyb1c as select reachnmb, count( cntyfip) as cnties from cntyb1 group by reachnmb; 237 quit;

proc print data= cntyb1a; title3 ' Number of affected counties Option 1'; title4 ' Sample'; quit;

proc print data= cntyb2a; title3 ' Number of affected counties Option 2'; title4 ' Sample'; quit;

proc means sum data= cntyb1b; var wght; title3 ' Sum of Counties times Average weight'; quit;

proc means data= cntyb1b; var wght; title3 ' Mean County weight'; quit;

proc means data= cntyb1c; var cnties; title3 ' Average counties traversed by a reach'; quit;

data check; set benefit2;

if bfrac_ 1 <= 0 then delete; run;

data check2; merge check pops2; by reachnmb; if weight=. then delete;

keep reachnmb pop totview totboat viewers swimers weight; run;