******************************************************************************************************************; * Example SUDAAN code to produce population totals associated with Figures 1 & 2, NCHS Data Brief No. 405 *; *; * Osteoporosis of Low Bone Mass in Older Adults:United States, 2017-2018 *; * *; * Sarafrazi N, Wambogo EA, Shephered JA. Osteoporosis of Low Bone Mass in Older Adults:United States, 2017-2018 *; * NCHS Data Brief. No 405. Hyattsville, MD: National Center for Health Statistics. 2021. *; * *; * Available at: https://www.cdc.gov/nchs/products/databriefs/db405.htm *; ******************************************************************************************************************; ** Indicate location of your files **; %global home macros titles liblabel; %let home = .\;/*Insert path to Home folder where the macro is saved*/ /*** Include macros ***/ /*** This code assumes all macros are stored in a subfolder of home called 'macros' ***/ %include "&home\macros\KG_macro.sas"; options nocenter nodate nonumber pagesize=100 linesize=150; OPTIONS FORMCHAR="|----|+|---+=|-/\<>*"; %put Run in SAS &sysver (maintenance release and release year: &sysvlong4); ** Macro To Download Data from NHANES website **; %macro CreateDS(myDS); %let i = 1; %let DS = %scan(&myDS, &i); %do %until(&DS = %nrstr()); %let Suffix = %lowcase(%substr(&DS, %eval(%length(&DS)-1))); %if (&Suffix = _j) %then %do; filename &DS url "https://wwwn.cdc.gov/nchs/nhanes/2017-2018/&DS..xpt"; %end; %else %if (&Suffix = _i) %then %do; filename &DS url "https://wwwn.cdc.gov/nchs/nhanes/2015-2016/&DS..xpt"; %end; %else %if (&Suffix = _h) %then %do; filename &DS url "https://wwwn.cdc.gov/nchs/nhanes/2013-2014/&DS..xpt"; %end; %else %if (&Suffix = _g) %then %do; filename &DS url "https://wwwn.cdc.gov/nchs/nhanes/2011-2012/&DS..xpt"; %end; %else %if (&Suffix = _f) %then %do; filename &DS url "https://wwwn.cdc.gov/nchs/nhanes/2009-2010/&DS..xpt"; %end; %else %if (&Suffix = _e) %then %do; filename &DS url "https://wwwn.cdc.gov/nchs/nhanes/2007-2008/&DS..xpt"; %end; %else %if (&Suffix = _d) %then %do; filename &DS url "https://wwwn.cdc.gov/nchs/nhanes/2005-2006/&DS..xpt"; %end; %else %if (&Suffix = _c) %then %do; filename &DS url "https://wwwn.cdc.gov/nchs/nhanes/2003-2004/&DS..xpt"; %end; %else %if (&Suffix = _b) %then %do; filename &DS url "https://wwwn.cdc.gov/nchs/nhanes/2001-2002/&DS..xpt"; %end; %else %do; filename &DS url "https://wwwn.cdc.gov/nchs/nhanes/1999-2000/&DS..xpt"; %end; libname &DS xport; data &DS; set &DS..&DS; run; %let i = %eval(&i+1); %let DS = %scan(&myDS, &i); %end; %mend CreateDS; ** download datasets **; %CreateDS( demo_j DXX_j DXXFEM_J DXXSPN_J ); /*Import edited population totals data - Unedited data located at https://wwwn.cdc.gov/nchs/nhanes/ResponseRates.aspx*/ PROC IMPORT OUT=ACS DBMS=XLSX DATAFILE=".\ACS-Population-Totals-For-2017-2018" REPLACE; GETNAMES=YES; sheet="Male_Female"; RUN; PROC SORT DATA=dxxspn_j;by seqn;run; PROC SORT DATA=dxxfem_j;by seqn;run; PROC SORT DATA=demo_j;by seqn;run; data dexa; merge DEMO_J(keep=seqn riagendr ridageyr ridreth1 RIDRETH3 ridexprg sdmvstra sdmvpsu wtmec2yr sddsrvyr) DXXSPN_J(keep=seqn DXXL1BCC DXXL2BCC DXXL3BCC DXXL4BCC DXXOSBMD DXXL1BMD DXXL2BMD DXXL3BMD DXXL4BMD) DXXFEM_J(keep=seqn DXXNKBMD); by seqn; run; proc format; value sexfmt 1="Males" 2="Females" 0="All"; value agecatfmt 0='ALL' 1="50-59" 2="60-69" 3="70-79" 4="80+"; value agecat2fmt 0='ALL' 1="50-59" 2="60-69" 3="70+"; value agecat3fmt 0='ALL' 1="50-64" 2="65+"; VALUE racefmt 0 = 'ALL' 1 = 'NHW' 2 = 'NHB' 3 = 'Hispanic' 4 = 'Asian' 5 = 'Other'; value race2fmt 0='All race/ethnic groups' 1='NHW' 2='NHB' 3='HISPANIC' 4='Other, inc Asian'; VALUE race3fmt 0 = 'ALL' 1 = 'NHW' 2 = 'NHB' 3 = 'Hispanic' 4 = 'Asian' 5 = 'Other'; value kcyc6fmt 0='ALL' 1='NHANES 2005-2006' 2='NHANES 2007-2008' 3='NHANES 2009-2010' 4='NHANES 2013-2014' 5='NHANES 2017-2018'; value femurfmt 0 ='ALL' 1 ='Osteoporosis' 2 ='Low_bone_mass' 3 ='Normal'; value spnstatfmt 0 ='ALL' 1 ='Osteoporosis_spine' 2 ='LBM_spine' 3 ='Normal_spine'; value femstatfmt 0 ='ALL' 1 ='Osteoporosis_femur' 2 ='LBM_femur' 3 ='Normal_femur'; value bonstatfmt 0 ='ALL' 1 ='Osteoporosis_either_site' 2 ='LBM_both_or_either_site' 3 ='Normal_both_sites' ; run; data DEXA2; set DEXA; /*SAS code to calculate femoral neck T-Score, Lumbar spine BMD and Lumbar spine T-score Selected variables of interest to create: DXXNKBMD - Femoral neck BMD DXXL1BCC - L1 BMD invalidity code DXXL2BCC - L2 BMD invalidity code DXXL3BCC - L3 BMD invalidity code DXXL4BCC - L4 BMD invalidity code DXXOSBMD - Total spine BMD DXXL1BMD - L1 BMD DXXL2BMD - L2 BMD DXXL3BMD - L3 BMD DXXL4BMD - L4 BMD*/ /*1. FN_T:*/ fn_t= (dxxnkbmd - 0.86)/0.12; /*Total_femur_T_score= (dxxofbmd - 0.94)/0.122; */ /*2. LSBMD and LS_T*/ *create counting variable for invalid vertebra; if dxxl1bcc gt 0 then invalid_l1=1; if dxxl1bcc eq 0 then invalid_l1=0; if dxxl2bcc gt 0 then invalid_l2=1; if dxxl2bcc eq 0 then invalid_l2=0; if dxxl3bcc gt 0 then invalid_l3=1; if dxxl3bcc eq 0 then invalid_l3=0; if dxxl4bcc gt 0 then invalid_l4=1; if dxxl4bcc eq 0 then invalid_l4=0; invalid_sum = invalid_l1 + invalid_l2 + invalid_l3 + invalid_l4; *Calculate spine t scores using Hologic white female AP spine reference data; ****AFTER RECALCULATING TOTAL SPINE BMD TO ALLOW FOR 4 VALID VERT!!!!!!!!; *l1,2,3,4; if invalid_sum eq 0 then LS_T = (dxxosbmd - 1.047)/0.11; if invalid_sum eq 0 then LSBMD=DXXOSBMD; *L2,3,4; if invalid_sum eq 1 and invalid_l1 eq 1 THEN TOTSPINEBMD_234=(dxxl2bmd + dxxl3bmd + dxxl4bmd)/3; if invalid_sum eq 1 and invalid_l1 eq 1 then LS_T = (totspinebmd_234 - 1.079)/0.11; if invalid_sum eq 1 and invalid_l1 eq 1 THEN LSBMD= TOTSPINEBMD_234; *L1,3,4; if invalid_sum eq 1 and invalid_l2 eq 1 THEN TOTSPINEBMD_134=(dxxl1bmd + dxxl3bmd + dxxl4bmd)/3; if invalid_sum eq 1 and invalid_l2 eq 1 THEN LS_T = (totspinebmd_134 - 1.053)/0.11; if invalid_sum eq 1 and invalid_l2 eq 1 THEN LSBMD=TOTSPINEBMD_134; *L1,2,4; if invalid_sum eq 1 and invalid_l3 eq 1 THEN TOTSPINEBMD_124=(dxxl1bmd + dxxl2bmd + dxxl4bmd)/3; if invalid_sum eq 1 and invalid_l3 eq 1 then LS_T = (totspinebmd_124 - 1.034)/0.11; if invalid_sum eq 1 and invalid_l3 eq 1 THEN LSBMD=TOTSPINEBMD_124; *L1,2,3; if invalid_sum eq 1 and invalid_l4 eq 1 THEN TOTSPINEBMD_123=(dxxl1bmd + dxxl2bmd + dxxl3bmd)/3; if invalid_sum eq 1 and invalid_l4 eq 1 then LS_T = (totspinebmd_123 - 1.018)/0.11; if invalid_sum eq 1 and invalid_l4 eq 1 THEN LSBMD=TOTSPINEBMD_123; *L3,4; if invalid_sum eq 2 and invalid_l1 eq 1 and invalid_l2 eq 1 THEN TOTSPINEBMD_34=(dxxl3bmd + dxxl4bmd )/2; if invalid_sum eq 2 and invalid_l1 eq 1 and invalid_l2 eq 1 then LS_T = (totspinebmd_34 - 1.101)/0.11; if invalid_sum eq 2 and invalid_l1 eq 1 and invalid_l2 eq 1 THEN LSBMD=TOTSPINEBMD_34; *L2,4; if invalid_sum eq 2 and invalid_l1 eq 1 and invalid_l3 eq 1 THEN TOTSPINEBMD_24=(dxxl2bmd + dxxl4bmd)/2; if invalid_sum eq 2 and invalid_l1 eq 1 and invalid_l3 eq 1 then LS_T = (totspinebmd_24 - 1.077)/0.11; if invalid_sum eq 2 and invalid_l1 eq 1 and invalid_l3 eq 1 THEN LSBMD=TOTSPINEBMD_24; *L2,3; if invalid_sum eq 2 and invalid_l1 eq 1 and invalid_l4 eq 1 THEN TOTSPINEBMD_23=(dxxl2bmd + dxxl3bmd)/2; if invalid_sum eq 2 and invalid_l1 eq 1 and invalid_l4 eq 1 then LS_T = (totspinebmd_23 - 1.058)/0.11; if invalid_sum eq 2 and invalid_l1 eq 1 and invalid_l4 eq 1 THEN LSBMD=TOTSPINEBMD_23; *L1,4; if invalid_sum eq 2 and invalid_l2 eq 1 and invalid_l3 eq 1 THEN TOTSPINEBMD_14=(dxxl1bmd + dxxl4bmd)/2; if invalid_sum eq 2 and invalid_l2 eq 1 and invalid_l3 eq 1 then LS_T =(totspinebmd_14 - 1.037)/0.11; if invalid_sum eq 2 and invalid_l2 eq 1 and invalid_l3 eq 1 THEN LSBMD=TOTSPINEBMD_14; *L1,3; if invalid_sum eq 2 and invalid_l2 eq 1 and invalid_l4 eq 1 THEN TOTSPINEBMD_13=(dxxl1bmd + dxxl3bmd)/2; if invalid_sum eq 2 and invalid_l2 eq 1 and invalid_l4 eq 1 then LS_T = (totspinebmd_13 - 1.013)/0.11; if invalid_sum eq 2 and invalid_l2 eq 1 and invalid_l4 eq 1 THEN LSBMD=TOTSPINEBMD_13; *L1,2; if invalid_sum eq 2 and invalid_l3 eq 1 and invalid_l4 eq 1 THEN TOTSPINEBMD_12=(dxxl1bmd + dxxl2bmd)/2; if invalid_sum eq 2 and invalid_l3 eq 1 and invalid_l4 eq 1 then LS_T = (totspinebmd_12 - 0.979)/0.11; if invalid_sum eq 2 and invalid_l3 eq 1 and invalid_l4 eq 1 THEN LSBMD=TOTSPINEBMD_12; if ridageyr GE 50 and dxxnkbmd >. or (invalid_sum ne . and invalid_sum le 2) then selb=1; else selb=0; /*Create variables to define skeletal status at the femur neck or lumbar spine*/ /*Create a 3-category variable for femur neck status:*/ /*NOTE: If FN_T eq . then femur neck status variable = .;*/ /*femstat=1 (osteoporosis) femstat=2 (low bone mass) femstat=3 (normal)*/ if FN_T le -2.5 then femstat= 1; if FN_T gt -2.5 and FN_T lt -1.0 then femstat= 2; if FN_T ge -1.0 then femstat= 3; if FN_T eq . then femstat=.; /*Create a 3-category variable for lumbar spine status*/ /*NOTE: if LS_T = . then lumbar spine status variable should = .;*/ /*spnstat=1 (osteoporosis) spnstat=2 (low bone mass) spnstat=3 (normal)*/ if LS_T le -2.5 then spnstat = 1; if LS_T gt -2.5 and LS_T lt -1.0 then spnstat= 2; if LS_T ge -1.0 then spnstat= 3; if LS_T = . then spnstat=.; /*Create a 3-category variable to combine the femur neck and lumbar spine status variables */ /*NOTE: This analysis should be limited to respondents with non-missing data for femur neck status and lumbar spine status */ /*Category 1: Osteoporosis at EITHER femur neck or lumbar spine: */ /*Category 2: Low bone mass at both skeletal sites OR low bone mass at one site, normal at other site: */ /*NOTE: we do not want to count someone with osteoporosis at one site and low bone mass at the other site in this category*/ /*Category 3: Normal at both skeletal sites:*/ if femstat= 1 or spnstat= 1 then bonstat=1;/*osteo_either_sites*/ if (femstat= 2 and spnstat= 2) OR (femstat= 2 and spnstat=3) OR (femstat= 3 and spnstat= 2) then bonstat=2;/*LBM either site*/ if femstat= 3 and spnstat= 3 then bonstat=3;/*Normal both site*/ ***Construct sample MEC weight for the 5 cycles; if sddsrvyr in (4,5,6,8,10) then wtmec2yr2= wtmec2yr*1/5; /* for 2005-2018 */ if sddsrvyr=4 then kcyc6=1; if sddsrvyr=5 then kcyc6=2; if sddsrvyr=6 then kcyc6=3; if sddsrvyr=8 then kcyc6=4; if sddsrvyr=10 then kcyc6=5; agecat=.; if 50<=ridageyr le 59 then ageCat=1; if 60<=ridageyr le 69 then ageCat=2; if 70<=ridageyr le 79 then ageCat=3; if ridageyr ge 80 then ageCat=4; agecat2=.; if 50<=ridageyr le 59 then ageCat2=1; if 60<=ridageyr le 69 then ageCat2=2; if ridageyr ge 70 then ageCat2=3; agecat3=.; if 50<=ridageyr le 64 then ageCat3=1; if ridageyr ge 65 then ageCat3=2; sex=riagendr; *RACE; * recode race/ethnicity; race=.; if RIDRETH3=3 then race=1; *Non-Hispanic White; if RIDRETH3=4 then race=2; *Non-Hispanic Black; if RIDRETH3 in (1,2) then race=3; *HISPANIC**********; if RIDRETH3=6 then race=4; *Asian*************; if RIDRETH3=7 then race=5; *Other*************; race2=.; if RIDRETH1=3 then race2=1; *Non-Hispanic White; if RIDRETH1=4 then race2=2; *Non-Hispanic Black; if RIDRETH1 in (1,2) then race2=3; *HISPANIC**********; if RIDRETH1 in (5) then race2=4; *Other, including Asian*************; race3=.; if RIDRETH1=3 then race3=1; *Non-Hispanic White; if RIDRETH1=4 then race3=2; *Non-Hispanic Black; if RIDRETH1 in (1,2) then race3=3; *HISPANIC**********; if RIDRETH3 in (6) then race3=4; * Asian*************; if RIDRETH3=7 then race3=5; *Other*************; LABEL ageCat = 'AGE GROUP' race = 'Race ethnicity' riagendr = 'Gender' ; format agecat agecatfmt. agecat2 agecat2fmt. sex sexfmt. race racefmt. race2 race2fmt. race3 race3fmt. bonstat bonstatfmt. spnstat spnstatfmt. femstat femstatfmt. kcyc6 kcyc6fmt. agecat3 agecat3fmt. ; RUN; *************************************************************************************; *CRUDE PROPORTIONS AND SAMPLE CHI-SQUARE TESTS* *************************************************************************************; proc sort data = dexa2; by sdmvstra sdmvpsu; run; proc crosstab data=dexa2 design=wr; * Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*; subpopx selb=1 /NAME=" "; * Nest statement: PSUs nested within Strata accounts for the design effects*; nest SDMVSTRA SDMVPSU; * Weight statement: specify appropriate weight, accounts for the unequal probability of sampling and non-response.*; weight wtmec2yr ; * Class statement: specify the analysis categorical variable(s)*; class bonstat agecat3 sex; * Table(s) statement: specifies cross-tabulations for which estimates are requested. If a table statement is not present, a one-dimensional distribution is generated.*; tables (agecat3 sex agecat3*sex)*bonstat; *test statement: requests chi-square statistics*; test chisq llchisq; setenv rowwidth=12 colwidth=10 labwidth=25; print nsum="samsize" rowper serow colper secol lowrow uprow/stest=default; output nsum rowper serow colper secol lowrow uprow /replace nsumfmt=f8.0 wsumfmt=f15.10 rowperfmt=f15.10 serowfmt=f15.10 colperfmt=f15.10 secolfmt=f15.10 uprowfmt=f15.10 lowrowfmt=f15.10 filename="Crosstabs"; TITLE "NHANES 2017-2018 "; rformat sex sexfmt.; rformat ageCat3 agecat3fmt.; rformat spnstat spnstatfmt.; rformat femstat femstatfmt.; rformat bonstat bonstatfmt.; RTITLE "Estimates of prevalence of osteoporosis" "Test Null Hyp: Prevalence of osteoporosis" Footnote "NHANES"; run; proc print data=Crosstabs noobs label; run; **********************************************************************; *PROPORTIONS**; **********************************************************************; *****STEP 1*******; PROC DESCRIPT DATA=dexa2 atlevel1=1 atlevel2=2 DESIGN=WR; * Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*; subpopx selb=1 /NAME=" "; * Nest statement: PSUs nested within Strata accounts for the design effects*; nest SDMVSTRA SDMVPSU; * Weight statement: specify appropriate weight, accounts for the unequal probability of sampling and non-response.*; *weight wtmec2yr2; weight wtmec2yr; * Subgroup statement: specify categorical variable(s), and number of categories of the variable(s)*; SUBGROUP agecat3 sex; LEVELS 2 2; * Var statement: specify the analysis variable(s),and the levels of the analysis variable(s)*; var bonstat ; catlevel 1; * Table(s) statement: specifies cross-tabulations for which estimates are requested. If a table statement is not present, a one-dimensional distribution is generated.*; TABLE agecat3 sex agecat3*(sex) ; PRINT NSUM="n" percent="%" SEpercent="SE"/ NOHEAD NOTIME STYLE=NCHS NSUMFMT=F8.0 percentFMT=F8.1 SEpercentFMT=F8.2 uppctFMT=F8.1 lowpctFMT=F8.2; rtitle "Prevalence of Osteoporosis"; RTITLE "NHANES 2017-2018"; * output statement to output the number of observations (nsum), percent, standard error of the percent (sepercent), number of strata (atlev1), and number of PSUs (atlev2) to a SAS file named osteo_mean**; output nsum percent sepercent uppct lowpct atlev1 atlev2 /filename=Osteo replace; rformat agecat3 AGEcat3FMT.; rformat sex SEXFMT.; rformat race racefmt.; RUN; ******************************************************; * Calculating Wald Confidence Intervals*; ******************************************************; DATA Osteo2; SET Osteo; df=atlev2-atlev1; * Use tcritl and tcritu to calculate the t-statistics based on the calculated degrees of freedom (df).*; tcritl=tinv( .025 ,df); tcritu=tinv(.975 ,df); * Calculate the lower limit (ll), and upper limit (ul) for the Wald 95% confidence intervals.*; ll=round((percent+tcritl*sepercent),.01); ul=round((percent+tcritu*sepercent),.01); * Percent statement rounds the estimates and renames them to percent and sepercent (standard error of percent).*; percent=round(percent,.01); sepercent=round(sepercent,.01); run ; proc format; value varlabel 1="Osteoporosis" ; run; proc print data=Osteo2 noobs label; format variable varlabel.; run; quit; ods tagsets.ExcelXP close; ****************************************************************; *SAMPLE POPULATION TOTALS FOR OSTEOPOROSIS PREVALENCE***; ****************************************************************; ******STEP 2**************; PROC descript data= ACS design=srs noprint means; * Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*; SUBPOPX selb = 1; * Subgroup statement: specify categorical variable(s), and number of categories of the variable(s)*; SUBGROUP sex agecat3; LEVELS 2 2; * Var statement: specify the analysis variable(s)*; VAR Total; * Table(s) statement: specifies cross-tabulations for which estimates are requested. If a table statement is not present, a one-dimensional distribution is generated.*; TABLES sex*agecat3 ; OUTPUT total / FILENAME=pt1718 FILETYPE=SAS REPLACE; Run ; proc sort data =pt1718; by sex agecat3; run ; ****************************************************************; * *POPULATION COUNTS WITH KG Confidence Intervals *; ****************************************************************; **STEP 3****; *Data from previous step will be fed into a SAS program where degrees of freedom, t-statistics, and 95% confidence limits will be calculated for each prevalence estimate, as shown below. Values also will be labeled and formatted. NOTE: The population totals are supplied as an Excel dataset.***; proc sort data =Osteo2;by sex agecat3; run ; data comb; merge Osteo2( in =a drop=ul ll rename=(uppct=ul lowpct=ll)) pt1718 ; by sex agecat3 ; if a ; if sex=-2 or agecat3=-2 then delete; popmean=(percent/100 )*total ; popl=ll/100 *total ; popu=ul/100 *total ; poplr=round(popl,1000); popur=round(popu,1000); popmeanr=round(popmean,1000); totalr=round(total,1000) ;run; *Call the KG macro*; %calculate_KornGraubard_CIs(roundingIncrement=0.01); proc sort data =Osteo_kg;by sex agecat3; run ; data comb_kg; merge Osteo_kg( in =a rename=(lowerCL=ll upperCL=ul)) pt1718 ; by sex agecat3 ; if a ; if sex=-2 or agecat3=-2 then delete; popmean=(percent/100 )*total ; popl=ll/100 *total ; popu=ul/100 *total ; poplr=round(popl,1000); popur=round(popu,1000); popmeanr=round(popmean,1000); totalr=round(total,1000) ; run; %let CItype=KG; proc print noobs split = '/' double ; var sex agecat3 percent sepercent ll ul df nsum totalr popmeanr poplr popur ; format sex sexfmt. agecat3 agecat3fmt. nsum 5.0 percent 5.2 sepercent 5.2 ll ul 5.2 df 2.0 totalr popmeanr poplr popur comma12.0; label percent='%' / 'with' / 'Osteoporosis' nsum='Num' / 'Osteoporosis' / 'status' sepercent='Std' / 'error' ll='Lower' / '95%' / "&CItype." / 'CI' ul='Upper' / '95%' / "&CItype." / 'CI' df='degrees' / 'of' / 'freedom' popmeanr='Pop' / 'Est' / 'US' / 'with' / 'Osteoporosis' totalr='Pop' / 'total' / 'US' poplr='Pop Est' / 'Lower' / '95%' / "&CItype." / 'CI' popur='Pop Est' / 'Upper' / '95 %' / "&CItype." / 'CI' ; title1 'Prevalence of persons with Osteoporosis - US, 2017-18' ; title2 "Percent and population estimates of number with Osteoporosis-&CItype. CI" ; run ;