************************************************************************* * CALCULATION OF THE NUMBER NEEDED TO BE EXPOSED (NNE) * * AND 95% CONFIDENCE INTERVALS * * WITH ADJUSTMENT FOR CONFOUNDING VARIABLES * * * * REFERENCE: BENDER & BLETTNER (2002), J. CLIN. EPIDEMIOL. 55: 525-530 * *************************************************************************; *------------------------------------------------------------------* | QUESTIONS ? | | | | ---> Dr. Ralf BENDER Phone: +49 221 35685-451 | | Fax: +49 221 35685-891 | | E-Mail: Ralf@rbsd.de | *------------------------------------------------------------------*; *------------------------------------------------------------------* !! NOTE: !! !! !! !! This code is appropriate for the SAS Version 8 or higher. !! !! You have to replace the simulated data set below with your !! !! own data. The SAS name of your data set should be "cohort". !! !! !! !! This data set must contain the following variables: !! !! 1) RESPONSE = adverse event of interest (1=yes, 0=no) !! !! 2) EXPOSURE = binary exposure variable (1=yes, 0=no) !! !! 3) COVAR_1 = the first confounding variable !! !! 4) If there are more than 1 confounding variables these !! !! variables must be named COVAR_2, COVAR_3, and so on. !! !! !! !! In line 39 you have to enter the number of confounding !! !! variables (e.g. %LET K=3) ! !! !! !! !! It is necessary to investigate the goodness-of-fit of the !! !! logistic regression model to get useful NNE estimates ! !! *------------------------------------------------------------------*; options linesize=105; *----->; %LET K=1 ; *** K = Number of covariables (excluding exposure); *--------------------------------------------------------------------* | Simulation of Cohort Data | *--------------------------------------------------------------------*; title1 'Calculation of Adjusted NNEs with 95% Confidence Intervals'; title4 'Simulated Cohort Data (n=5000, OR=3, 1 Covariate)'; title5 ' '; data sim; do i=1 to 5000; No=i; alpha=-10; gamma=1.0986123; beta1=0.1823216; Z=RANBIN(0,1,0.5); if Z=0 then X1=7*RANNOR(0)+45; if Z=1 then X1=7*RANNOR(0)+40; prob=exp(alpha+gamma*Z+beta1*X1)/(1+exp(alpha+gamma*Z+beta1*X1)); Y=RANBIN(0,1,prob); output; end; drop i alpha gamma beta1; run; *--------------------------------------------------------------------* | True Parameters of Simulated Cohort Data | *--------------------------------------------------------------------*; title4 'True Mean UER, Hypothetical EER, True OR, True ARI, and True NNEH'; title5 ' '; data UE; set sim; if Z=0; alpha=-10; gamma=1.0986123; beta1=0.1823216; UER=exp(alpha+beta1*X1)/(1+exp(alpha+beta1*X1)); EER=exp(alpha+beta1*X1+gamma)/(1+exp(alpha+beta1*X1+gamma)); ARI=EER-UER; keep UER EER ARI; run; proc means data=UE noprint; var UER EER ARI; output out=UEmean mean=UER_True EER_Hyp ARI_True; run; data TRUE; set UEmean; OR_True = exp(1.098612); NNEH_True = 1/ARI_True; format NNEH_True 16.8; run; proc print data=TRUE; var UER_True EER_Hyp OR_True ARI_True NNEH_True; run; *--------------------------------------------------------------------* | Descriptive Anlysis of Generated Cohort Data | *--------------------------------------------------------------------*; title4 'Descriptive Anlysis of Generated Cohort Data'; title5 ' '; data cohort; set sim; RESPONSE = Y; EXPOSURE = Z; COVAR_1 = X1; drop Y Z X1 PROB; run; proc means data=cohort n min max mean std sum maxdec=3 fw=8; proc means data=cohort n min max mean std sum maxdec=3 fw=8; class EXPOSURE; run; proc corr data=cohort spearman; var EXPOSURE COVAR_1; run; *--------------------------------------------------------------------* | Simple 2x2 Table | *--------------------------------------------------------------------*; title4 'Simple 2x2 Table with Crude Estimation of OR'; title5 ' '; proc freq data=cohort; tables EXPOSURE*RESPONSE / nopercent nocol cmh exact; run; *--------------------------------------------------------------------* | Multiple Logistic Regression Model | *--------------------------------------------------------------------*; title4 'Estimation of Adjusted Odds Ratio (Logistic Regression)'; title5 ' '; proc logistic data=cohort descending outest=est covout; model RESPONSE = EXPOSURE COVAR_1-COVAR_&K / rl lackfit rsq; output out=risks p=risk l=low u=up xbeta=linpred; run; *--------------------------------------------------------------------* | Data Sets Used for SAS/IML | *--------------------------------------------------------------------*; data COVARIAT; set cohort; keep COVAR_1-COVAR_&K; if EXPOSURE=0; run; data beta; set est; if _type_='PARMS'; drop _LINK_ _TYPE_ _NAME_ _LNLIKE_; run; data beta2_K; set beta; keep COVAR_1-COVAR_&K; run; data COVB; set est; if _type_='COV'; drop _LINK_ _TYPE_ _NAME_ _LNLIKE_; run; data covb2_k; set covb; keep COVAR_1-COVAR_&K; run; *--------------------------------------------------------------------* | PROC IML | *--------------------------------------------------------------------*; title4 'Crude and Adjusted Results Calculated with SAS/IML'; title5 ' '; proc IML; use cohort; read all var{RESPONSE} into RESPONSE; use cohort; read all var{RESPONSE} where (EXPOSURE=0) into RESP_UNE; use cohort; read all var{RESPONSE} where (EXPOSURE=1) into RESP_EXP; use cohort; read all var{EXPOSURE} into EXPOSURE; use covariat; read all var _all_ into COVARIAT; use beta; read all var{Intercept} into BETA0; use beta; read all var{EXPOSURE} into BETA1; use beta2_k; read all var _all_ into BETA2_k; use covb; read all var{Intercept EXPOSURE} into COV1; use covb2_k; read all var _all_ into COV2; COVB = (COV1 || COV2); N = NROW(RESPONSE); EVENT = SUM(RESPONSE); Ne = SUM(EXPOSURE); Nu = N-NE; Ee = SUM(RESP_EXP); Eu = SUM(RESP_UNE); K = NCOL(BETA2_k); UER = SUM(RESP_UNE)/NU; EER = SUM(RESP_EXP)/NE; ARI = EER-UER; NNEH = 1/ARI; print "GENERAL INFORMATION"; print "Sample size and number of events: " N " " EVENT; print "Number of unexposed and exposed subjects with events: " NU EU " " NE EE ,,,,; print "CRUDE (UNADJUSTED) RESULTS"; if ARI>0 then do; print "Crude event rates, ARI and NNEH: " UER[format=6.4] " " EER[format=6.4] " " ARI[format=6.4] " " NNEH[format=8.2] ,,,,; end; if ARI<0 then do; ARR = -ARI; NNEB = 1/ARR; print "Crude event rates, ARR and NNEB: " UER[format=6.4] " " EER[format=6.4] " " ARR[format=6.4] " " NNEB[format=8.2] ,,,,; end; if ARI=0 then do; NNEH="infinity"; print "Crude event rates, ARI and NNEH: " UER[format=6.4] " " EER[format=6.4] " " ARI[format=6.4] " " NNEH[format=8.2] ,,,,; end; D10 = 0; D11 = exp(BETA1); D12_k = 0#BETA2_k; D1 = (D10 || D11 || D12_k); LINPRED = BETA0 + COVARIAT*BETA2_k`; LOGIST2 = EXP(LINPRED)/((1+EXP(LINPRED))##2); DO J=1 to K; LOGISTJ = COVARIAT[,J] # EXP(LINPRED)/((1+EXP(LINPRED))##2); LOGIST = (LOGIST || LOGISTJ); END; D20 = SUM(LOGIST2)/NU; D21 = 0; D22_k = LOGIST[+,]/NU; D2 = (D20 || D21 || D22_k); D = (D1 // D2); COVF = D*COVB*D`; OR = exp(BETA1); RISK = EXP(LINPRED)/(1+EXP(LINPRED)); UER = SUM(RISK)/Nu; dv1 = UER#(1-UER) / ((1+UER#(OR-1))##2); dv2 = (OR-1)#(1-2#UER-(OR-1)#UER##2) / ((1+UER#(OR-1))##2); dv = (dv1 || dv2); ARI = (OR-1)#UER#(1-UER) / (1+UER#(OR-1)); EER = ARI+UER; print "ADJUSTED RESULTS (SIMPLE METHOD VIA ODDS RATIO)"; print "Attention: The simple method is only valid if the confounder distribution is not too wide!"; if ARI=0 then do; NNEH="infinity"; print "Absolute risk increase (ARI):" ARI[format=6.4]; print "Number needed to be exposed to be harmed (NNEH):" NNEH[format=8.2] ,,,; print "Sorry, but because ARI=0 the delta method is not applicable to calculate the standard error of ARI!"; abort; end; else NNEH = 1/ARI; VAR = dv*COVF*dv`; SE = SQRT(VAR); sn = 0.05; Zsn = PROBIT(1-SN/2); ARI_L = ARI - Zsn#SE; ARI_U = ARI + Zsn#SE; NNEH_L = 1/ARI_U; NNEH_U = 1/ARI_L; print "OR, observed unexposed (UER) and hypothetical exposed event rate (EER): " OR[format=7.3] " " UER[format=6.4] " " EER[format=6.4]; if ARI>0 then do; print "Absolute risk increase (ARI) with SE and 95% CI: " ARI[format=6.4] " " SE[format=8.6] " " ARI_L[format=6.4] ARI_U[format=6.4]; print "Number needed to be exposed to be harmed (NNEH) with 95% CI: " NNEH[format=8.2] " " NNEH_L[format=8.2] NNEH_U[format=8.2] ,,,,; end; if ARI<0 then do; ARR = -ARI; NNEB = 1/ARR; ARR_L = ARR - Zsn#SE; ARR_U = ARR + Zsn#SE; NNEB_L = 1/ARR_U; NNEB_U = 1/ARR_L; print "Absolute risk reduction (ARR) with SE and 95% CI: " ARR[format=6.4] " " SE[format=8.6] " " ARR_L[format=6.4] ARR_U[format=6.4]; print "Number needed to be exposed to benefit (NNEB) with 95% CI: " NNEB[format=8.2] " " NNEB_L[format=8.2] NNEB_U[format=8.2] ,,,,; end; print "ADJUSTED RESULTS (GENERAL METHOD VIA AVERAGING)"; print "Note: The general method is valid for narrow and wide confounder distributions"; LINPRED1 = BETA0 + COVARIAT*BETA2_k`; LINPRED2 = BETA0 + BETA1 + COVARIAT*BETA2_k`; LOGIST1 = EXP(LINPRED1) / ((1+EXP(LINPRED1))##2); LOGIST2 = EXP(LINPRED2) / ((1+EXP(LINPRED2))##2); d1 = SUM(LOGIST2-LOGIST1)/Nu; d2 = SUM(LOGIST2)/Nu; DO J=1 to K; LOGCOVJ = COVARIAT[,J] # (LOGIST2-LOGIST1); LOGCOV = (LOGCOV || LOGCOVJ); END; d3 = LOGCOV[+,]/Nu; d = (d1 || d2 || d3); RISK1 = EXP(LINPRED1)/(1+EXP(LINPRED1)); RISK2 = EXP(LINPRED2)/(1+EXP(LINPRED2)); UER = SUM(RISK1)/Nu; EER = SUM(RISK2)/Nu; ARI = EER-UER; if ARI=0 then do; NNEH="infinity"; print "Absolute risk increase (ARI):" ARI[format=6.4]; print "Number needed to be exposed to be harmed (NNEH):" NNEH[format=8.2] ,,,; print "Sorry, but because ARI=0 the delta method is not applicable to calculate the standard error of ARI!"; abort; end; else NNEH = 1/ARI; VAR = d*COVB*d`; SE = SQRT(VAR); sn = 0.05; Zsn = PROBIT(1-SN/2); ARI_L = ARI - Zsn#SE; ARI_U = ARI + Zsn#SE; NNEH_L = 1/ARI_U; NNEH_U = 1/ARI_L; print "Observed unexposed (UER) and hypothetical exposed event rate (EER): " UER[format=6.4] " " EER[format=6.4]; if ARI>0 then do; print "Absolute risk increase (ARI) with SE and 95% CI: " ARI[format=6.4] " " SE[format=8.6] " " ARI_L[format=6.4] ARI_U[format=6.4]; print "Number needed to be exposed to be harmed (NNEH) with 95% CI: " NNEH[format=8.2] " " NNEH_L[format=8.2] NNEH_U[format=8.2]; end; if ARI<0 then do; ARR = -ARI; NNEB = 1/ARR; ARR_L = ARR - Zsn#SE; ARR_U = ARR + Zsn#SE; NNEB_L = 1/ARR_U; NNEB_U = 1/ARR_L; print "Absolute risk reduction (ARR) with SE and 95% CI: " ARR[format=6.4] " " SE[format=8.6] " " ARR_L[format=6.4] ARR_U[format=6.4]; print "Number needed to be exposed to benefit (NNEB) with 95% CI: " NNEB[format=8.2] " " NNEB_L[format=8.2] NNEB_U[format=8.2]; end; quit; run;