options linesize=150;
*******************************************************************************************************;
* Example SUDAAN code to replicate mean ratios, NCHS Data Brief No. 270, Figures 1 - 5.               *;
* Alternative method of obtaining the estimates using ratio of means also shown.                      *;                      *;  
*                                                                                                     *;
* Rosinger A, Herrick K, Gahche J, Park S. Sugar-sweetened Beverage Consumption Among U.S. Adults,    *;
* 2011-2014. Data Brief. No 270. Hyattsville, MD: National Center for Health Statistics. 2017.        *;                                                                           *;
*                                                                                                     *;
*  Available at: https://www.cdc.gov/nchs/products/databriefs/db270.htm                               *;
*******************************************************************************************************;

options nocenter nodate nonumber pagesize=100 linesize=150;
OPTIONS FORMCHAR="|----|+|---+=|-/\<>*";


%put Run in SAS &sysver (maintenance release and release year: &sysvlong4)
 and SUDAAN Release 11.0.1 (SAS-Callable, 32 bit version);

** 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 xptIn url "https://wwwn.cdc.gov/nchs/nhanes/2017-2018/&DS..xpt"; %end;
    %else %if (&Suffix = _i) %then %do; filename xptIn url "https://wwwn.cdc.gov/nchs/nhanes/2015-2016/&DS..xpt"; %end;
    %else %if (&Suffix = _h) %then %do; filename xptIn url "https://wwwn.cdc.gov/nchs/nhanes/2013-2014/&DS..xpt"; %end;
    %else %if (&Suffix = _g) %then %do; filename xptIn url "https://wwwn.cdc.gov/nchs/nhanes/2011-2012/&DS..xpt"; %end;
    %else %if (&Suffix = _f) %then %do; filename xptIn url "https://wwwn.cdc.gov/nchs/nhanes/2009-2010/&DS..xpt"; %end;
    %else %if (&Suffix = _e) %then %do; filename xptIn url "https://wwwn.cdc.gov/nchs/nhanes/2007-2008/&DS..xpt"; %end;
    %else %if (&Suffix = _d) %then %do; filename xptIn url "https://wwwn.cdc.gov/nchs/nhanes/2005-2006/&DS..xpt"; %end;
    %else %if (&Suffix = _c) %then %do; filename xptIn url "https://wwwn.cdc.gov/nchs/nhanes/2003-2004/&DS..xpt"; %end;
    %else %if (&Suffix = _b) %then %do; filename xptIn url "https://wwwn.cdc.gov/nchs/nhanes/2001-2002/&DS..xpt"; %end;
    %else %do; filename xptIn url "https://wwwn.cdc.gov/nchs/nhanes/1999-2000/&DS..xpt"; %end;
    libname xptIn xport;

    data &DS;
      set xptIn.&DS;
    run;

    libname xptIn clear;
    filename xptIn clear;

    %let i = %eval(&i+1);
    %let DS = %scan(&myDS, &i);
  %end;
%mend CreateDS;

***************;
** Data prep **;
***************;

** Call macro to import datafiles for multiple individual components **;
* Demographic (DEMO) *;
* Individual food files (dr1iff) *;
* Reproductive health (RHQ) *;
%CreateDS(demo_g demo_h dr1iff_g dr1iff_h dr1tot_g dr1tot_h rhq_g rhq_h);


** Append data files from multiple survey cycles **;

      
*DEMOGRAPHIC DATA FILES;
data demo (keep=seqn sddsrvyr riagendr RIDEXPRG ridageyr ridreth3 sdmvpsu sdmvstra);
  set 
      demo_g
      demo_h;
run;

proc sort data=demo;
by seqn;
run;

PROC CONTENTS DATA=DEMO; RUN;


*TOTAL DIETARY INTAKE FILES;
data diet;
  set dr1tot_g dr1tot_h;
; RUN;

proc sort data=diet;
by seqn;
run;


*Reproductive health Questionnaire for lactation status;
data rhq;
  set 
    rhq_g (keep=seqn rhq200) 
    rhq_h (keep=seqn rhq200); RUN;

proc sort data=rhq;
by seqn;
run;

*INDIVIDUAL FOODS FILE -- Identify sugar-sweetened beverages ;
data SSB;
  set 
    dr1iff_g 
    dr1iff_H ; 
  IF DR1IFDCD IN (42401010,42403010,42404010,64200100,64201010,64201500,64202010,64203020,64204010,64205010,64210010,64213010,64215010,64221010,92101600,92101800,92101920,92101930,92101950,92101960,92121000,92121010,92121020, 
92130000,92130001,92130020,92301060,92301130,92301160,92301190,92302200,92302400,92302600,92302800,92305000,92305040,92305050,92305800,92306020,92306040,92306090,92306100,92400000,92410110,92410310,92410315,92410330,92410340,92410360, 
92410390,92410410,92410510,92410550,92410610,92410710,92410810,92411510,92417010,92431000,92432000,92433000,92510610,92510650,92510720,92510730,92511010,92511250,92512050,92512090,92512110,92530410,92530510,92530610,92530950,92531030, 
92541010,92542000,92550030,92550110,92550350,92552030,92560000,92582100,92582110,92582120,92610010,92610110,92611010,92611100,92611510,92611600,92612010,92613010,92803000,92804000,94100300,94210200,11561010,92101650,92101670,92101925,
92101935,92161005,92162005,92204000,92306200,92306610,92307500,92411520,92513000,92550610,92613510,95310200,95310400,95310500,95310550,95310555,95310560,95310600,95310700,95310750,95310800,95311000,95320200,95320500,95321000,95342000, 
74302000,92101820,92101921,92101923,92101926,92101928,92101931,92101933,92101936,92101938,92101955,92101965,92101970,92101975,92102000,92102010,92102020,92102030,92102040,92102050,92102060,92102070,92102080,92102090,92102100,92102110, 
92102450,92102600,92102601,92102602,92102610,92102611,92102612,92121001,92121050,92130001,92305910,92307400,92307510,92308000,92308030,92308500,92308530,92309000,92309010,92309500,92510955,92510960,92511015,92513010,92550200,92550360,
92550370,92550380,92610020,92610030)

 THEN DO;
    SSB_sugrG=DR1Isugr;
    SSB_KCAL=DR1IKCAL;
    SSB_YN=1;
END;
ELSE DO;
    SSB_sugrG=0;
    SSB_YN=0;
    SSB_KCAL=0;

END;
run;

PROC SORT DATA=SSB;
  BY SEQN;
 RUN;

* Sum across individual foods to get total number of SSBs, grams of sugar from SSBs, and calories from SSBs *;
PROC MEANS DATA=SSB noprint;
by seqn;
VAR SSB_sugrG SSB_KCAL SSB_YN ;
output out=SSB_GH sum=SSB_TsugrG SSB_TKCAL SSB_TYN;
id DR1DRSTZ  ;
RUN;

*MERGE data files;

DATA NHANES;
   MERGE  DEMO
      diet
      rhq
      SSB_GH
      ;
   BY SEQN;
   RUN;


proc format;
VALUE AGEkidf
0='All'
1='2-5'
2='6-11'
3='12-19';
VALUE AGEG3F
0='All'
1='20-39'
2='40-59'
3='60+';
value sexf
0='Both'
1='men' 
2='women';
VALUE race2fmt   
0 = 'ALL'
1 = 'NHW'
2 = 'NHB'
3 = 'Hispanic'
4 = 'NHA'
5 = 'Other';

run;

Data nhanes; set nhanes;

* calculate Dietary Day 1 4-year weight *;
if sddsrvyr in (7,8)then 
   weight  = 1/2 * wtdrd1 ; /* for 2011-2014 */

AGEG3=.;
IF (20<=ridageyr<40) THEN AGEG3=1; *20-39 yr;
IF (40<=ridageyr<60) THEN AGEG3=2; *40-59;
IF (ridageyr>=60) THEN AGEG3=3; *60+;

agekid=.;
IF (2<=ridageyr<6) THEN AGEkid=1; *2-5 yr;
IF (6<=ridageyr<12) THEN AGEkid=2; *6-11;
IF (12<=ridageyr<20) THEN AGEkid=3; *12-19;


sex=RIAGENDR;

*RACE;
* recode race/ethnicity; 
if sddsrvyr in (7,8) then do;

  if RIDRETH3=3 then race2=1;         *Non-Hispanic White;
  if RIDRETH3=4 then race2=2;         *Non-Hispanic Black;
  if RIDRETH3 in (1,2) then race2=3;  *HISPANIC**********;
  if RIDRETH3=6 then race2=4;     *Asian*************;
  if RIDRETH3=7 then race2=5;     *Other*************;
end;

*Adult cohort, excluding pregnant and lactating women;
incoh=.;
if sex=1 then do;
  if DR1DRSTZ=1 AND (RIDAGEYR>=20) then incoh=1; else incoh=0;
end;
if sex=2 and 20<=ridageyr<45 then do;
if DR1DRSTZ=1 AND RIDEXPRG ne 1 and rhq200 ne 1 then incoh=1; else incoh=0; 
end;
if sex=2 and ridageyr>=45 then do;
if DR1DRSTZ=1 and rhq200 ne 1 then incoh=1; else incoh=0;
end;

*Adult cohort, no exclusions beyond age and diet;
incohadult=.;
  if DR1DRSTZ=1 AND (RIDAGEYR>=20) then incohadult=1; else incohadult=0;

*Kids cohort, exclusions by age and diet;
incohK=.;
if DR1DRSTZ=1 AND (.<RIDAGEYR<20) then incohk=1; else incohk=0;

SSBn=SSB_TYN;
if SSBn>=3 then ssbn=3;

pE_SSB=(SSB_TKCAL/dr1tkcal)*100;

SSB_any=.;
if SSB_TYN>0 then SSB_any=100;
if SSB_TYN=0 then SSB_any=0;

FORMAT sex sexf. race2 race2fmt. AGEG3 AGEG3F. agekid agekidf.;
run;



***********************************************************************************************;
                       **PROPORTIONS***;
***********************************************************************************************;

*Figure 1 Percent of adults who consume any SSBs by sex;
*****************************************************************************;

PROC SORT DATA=WORK.nhanes ;
BY sdmvstra sdmvpsu ;                      
RUN;

PROC DESCRIPT DATA=WORK.nhanes atlevel1=1 atlevel2=2 DESIGN=WR;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
subpopx incohadult=1;
* 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 weight;
* Subgroup statement: specify other categorical variable(s), and number of categories of the variable(s)*;
SUBGROUP SEX;
LEVELS 2 ;
* Var statement: specify the analysis variable(s),and the levels of the analysis variable(s)*;
VAR ssb_any; 
catlevel 100;
* 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 SEX;
* Print statement: Prints the output*; 
PRINT NSUM="SAMPLE SIZE" percent="percent" sepercent="SE"/ NOHEAD NOTIME STYLE=NCHS  NSUMFMT=F8.0 percentFMT=F8.1
SEpercentFMT=F8.1;
rtitle "Percent of adults who consume ANY SSBs by gender";
RTITLE "NHANES 2011-2014";
output nsum percent sepercent  atlev1 atlev2 /filename=SSBany_Adult_sex replace;
RUN;


*******************************************************************************************;
                        ** Korn and Graubard confidence intervals **;
        *NOTE: percent and SEpercent are specified, based on the output statement above***
            *** NOTE: this code is not appropriate for age-standardized proportions ***         
*******************************************************************************************;
proc format;
  value varlabel 
  1="Any SSB" ;
run;

*To construct confidence intervals using the Korn and Graubard method;
*See pp 364 and 365 of the SUDAAN 11 manual;

data SSBany_Adult_sex_CI;
  set SSBany_Adult_sex;
  *p and its standard error should be decimal numbers between 0 and 1 ;
  p=percent/100; 
  sep=sepercent/100; 
  * complementary proportion *;
  q=1-p; 

  df_flag=0;
  df=atlev2-atlev1;
  if df<8 then df_flag=1; 

  * Effective sample size *;
  if (0<p<1) then n_eff=(p*(1-p))/(sep**2); 
  else n_eff=nsum;

  *Ratio of ts: adjustment to sample size suggested by Korn and Graubard for complex survey data; 
  *A two-sided a (0.05/2 or 0.025) is used in the equation below: 1-0.025 = 0.975 ; 
  if df > 0 then rat_squ=(tinv(.975,nsum-1)/tinv(.975,df))**2; 
  else rat_squ=0; *limit case: set to zero;
  *df-adjusted effective sample size (can be no greater than the sample size);
  if p > 0 then n_eff_df=min(nsum,rat_squ*n_eff); 
  else n_eff_df=nsum; *limit case: set to sample size;

  *Parameters for beta confidence limits; 
  x=n_eff_df*p; 
  v1=x; 
  v2=n_eff_df-x+1; 
  v3=x+1; 
  v4=n_eff_df-x;
  *lower and upper confidence limits for Korn and Graubard interval ;
  * Note: Using inverse beta instead of ratio of Fs for numerical efficiency ;
  *if (0<p<1), otherwise set lower limit to 0 when p=0 and upper limit to 1 when p=1 
  * A two-sided a (0.05/2 or 0.025) is used in the 
  equations below: 0.025 and 0.975; 
  if (v1=0) then kg_l=0; 
  else kg_l=betainv(.025,v1,v2); 
  if (v4=0) then kg_u=1; 
  else kg_u=betainv(.975,v3,v4);

  *Korn and Graubard CI absolute width;
  kg_wdth=kg_u - kg_l;
  *Korn and Graubard CI relative width for p; 
  if (p>0) then kg_relw_p=100*(kg_wdth/p); 
  else kg_relw_p=.;
  *Korn and Graubard CI relative width for q; 
  if (q>0) then kg_relw_q=100*(kg_wdth/q); 
  else kg_relw_q=.;

  * Application of the NCHS Data Presentation Standards for Proportions *;
  *Proportions with CI width <= 0.05 are reliable, unless...; 
  p_reliable=1; 
  *...effective sample size or nominal sample size is less than 30; 
  if min(nsum, n_eff) < 30 then p_reliable=0;
  *...absolute CI width is greater than or equal 0.30; 
  else if kg_wdth ge 0.30 then p_reliable=0;
  *...relative CI width is greater than 130%; 
  else if (kg_relw_p > 130 and kg_wdth > 0.05) then p_reliable=0;

  *Determine if estimate should be flagged as having an unreliable complement; 
  if (p_reliable=1) then do;
    *Complementary proportions are reliable, unless...; 
    q_reliable=1;
    *...relative CI width is greater than 130% ; 
    if (kg_relw_q > 130 and kg_wdth > 0.05) then q_reliable=0; 
  end;

  *Determine if estimate should be flagged for statistical review ; 
  p_statistical=0; 
  if p_reliable=1 then do;
    *Estimates with df < 8 or percents = 0 or 100 or unreliable complement are flagged for clerical or ADS review; 
    if df_flag=1 or p=0 or p=1 or q_reliable=0 then p_statistical =1; 
  end;
  format variable varlabel.;
run;

 proc print data=SSBany_Adult_sex_CI noobs label;
var  sex variable nsum percent sepercent p sep q kg_l kg_u p_reliable q_reliable p_statistical;
  format percent sepercent p sep q kg_l kg_u 8.3;
    /*format variable varlabel.;*/
footnote "Confidence limits are calculated by the method of Korn and Graubard using alternative example code";
run; 
quit;
ods tagsets.ExcelXP close;

title;
footnote;


*****************************************************************************;
   *Figure 1 Percent of adults who consume 0, 1, 2, 3 or more SSBs by sex;
*****************************************************************************;

PROC DESCRIPT DATA=WORK.nhanes atlevel1=1 atlevel2=2 DESIGN=WR;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
subpopx incohadult=1;
* 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 weight;
* Subgroup statement: specify other categorical variable(s), and number of categories of the variable(s)*;
SUBGROUP SEX;
LEVELS 2 ;
* Var statement: specify the analysis variable(s),and the levels of the analysis variable(s)*;
VAR SSBN ssbn ssbn ssbn; 
catlevel 0 1 2 3;
* 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 SEX;
* Print statement: Prints the output*; 
PRINT NSUM="SAMPLE SIZE" percent="percent" sepercent="SE"/ NOHEAD NOTIME STYLE=NCHS  NSUMFMT=F8.0 percentFMT=F8.1 
SEpercentFMT=F8.1;
rtitle "Percent of adults who consume 0, 1, 2, 3 or more SSBs by gender";
RTITLE "NHANES 2011-2014";
output nsum percent sepercent  atlev1 atlev2 /filename=SSBAdult_sex replace;
RUN;


*******************************************************************************************;
                           ** Korn and Graubard confidence intervals **;
         *NOTE: Percent and SEpercent are specified, based on the output statement above***
            *** NOTE: this code is not appropriate for age-standardized proportions ***         
*******************************************************************************************;
 proc format;
  value varlabel 
  1="No SSB"
  2="One SSB"
  3="Two SSB"
  4="three or more SSB";
run;

*To construct confidence intervals using the Korn and Graubard method;
*See pp 364 and 365 of the SUDAAN 11 manual;

data SSBAdult_sex_CI;
  set SSBAdult_sex;
  *p and its standard error should be decimal numbers between 0 and 1 ;
  p=percent/100; 
  sep=sepercent/100; 
  * complementary proportion *;
  q=1-p; 

  df_flag=0;
  df=atlev2-atlev1;
  if df<8 then df_flag=1; 

  * Effective sample size *;
  if (0<p<1) then n_eff=(p*(1-p))/(sep**2); 
  else n_eff=nsum;

  *Ratio of ts: adjustment to sample size suggested by Korn and Graubard for complex survey data; 
  *A two-sided a (0.05/2 or 0.025) is used in the equation below: 1-0.025 = 0.975 ; 
  if df > 0 then rat_squ=(tinv(.975,nsum-1)/tinv(.975,df))**2; 
  else rat_squ=0; *limit case: set to zero;
  *df-adjusted effective sample size (can be no greater than the sample size);
  if p > 0 then n_eff_df=min(nsum,rat_squ*n_eff); 
  else n_eff_df=nsum; *limit case: set to sample size;

  *Parameters for beta confidence limits; 
  x=n_eff_df*p; 
  v1=x; 
  v2=n_eff_df-x+1; 
  v3=x+1; 
  v4=n_eff_df-x;
  *lower and upper confidence limits for Korn and Graubard interval ;
  * Note: Using inverse beta instead of ratio of Fs for numerical efficiency ;
  *if (0<p<1), otherwise set lower limit to 0 when p=0 and upper limit to 1 when p=1 
  * A two-sided a (0.05/2 or 0.025) is used in the 
  equations below: 0.025 and 0.975; 
  if (v1=0) then kg_l=0; 
  else kg_l=betainv(.025,v1,v2); 
  if (v4=0) then kg_u=1; 
  else kg_u=betainv(.975,v3,v4);

  *Korn and Graubard CI absolute width;
  kg_wdth=kg_u - kg_l;
  *Korn and Graubard CI relative width for p; 
  if (p>0) then kg_relw_p=100*(kg_wdth/p); 
  else kg_relw_p=.;
  *Korn and Graubard CI relative width for q; 
  if (q>0) then kg_relw_q=100*(kg_wdth/q); 
  else kg_relw_q=.;

  * Application of the NCHS Data Presentation Standards for Proportions *;
  *Proportions with CI width <= 0.05 are reliable, unless...; 
  p_reliable=1; 
  *...effective sample size or nominal sample size is less than 30; 
  if min(nsum, n_eff) < 30 then p_reliable=0;
  *...absolute CI width is greater than or equal 0.30; 
  else if kg_wdth ge 0.30 then p_reliable=0;
  *...relative CI width is greater than 130%; 
  else if (kg_relw_p > 130 and kg_wdth > 0.05) then p_reliable=0;

  *Determine if estimate should be flagged as having an unreliable complement; 
  if (p_reliable=1) then do;
    *Complementary proportions are reliable, unless...; 
    q_reliable=1;
    *...relative CI width is greater than 130% ; 
    if (kg_relw_q > 130 and kg_wdth > 0.05) then q_reliable=0; 
  end;

  *Determine if estimate should be flagged for statistical review ; 
  p_statistical=0; 
  if p_reliable=1 then do;
    *Estimates with df < 8 or percents = 0 or 100 or unreliable complement are flagged for clerical or ADS review; 
    if df_flag=1 or p=0 or p=1 or q_reliable=0 then p_statistical =1; 
  end;
  format variable varlabel.;
run;



proc print data=SSBAdult_sex_CI noobs label;
var  sex variable nsum percent sepercent p sep q kg_l kg_u p_reliable q_reliable p_statistical;
  format percent sepercent p sep q kg_l kg_u 8.3;
    /*format variable varlabel.;*/
footnote "Confidence limits are calculated by the method of Korn and Graubard using alternative example code";
run; 
quit;
ods tagsets.ExcelXP close;

title;
footnote;


*********************************************************************************;
*T-TESTS;
*tests between sexes: ANY SSB;
**********************************************************************************;

PROC DESCRIPT data=WORK.NHANES design=wr;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
SUBPOPX incohadult=1;
* 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  weight;
* Class statement: specify categorical variable(s)*;
class sex  SSB_any/nofreq; 
* Var statement: specify the analysis variable(s),and the levels of the analysis variable(s)*;
var  SSB_any ; 
catlevel 100;
* Pairwise statement: Requests all pairwise contrasts*; 
pairwise sex;
SETENV LABWIDTH=35 COLWIDTH=10 ROWWIDTH=10 DECWIDTH=2 COLSPCE=2;
* Print statement: Prints the output*; 
print  nsum= "Sample Size" percent="percent" sepercent="SE" t_pct  p_pct="P-value"/nohead notime style=NCHS  
nsumfmt=f8.0  percentfmt=f8.1  sepercentfmt=f8.2 p_pctfmt=f8.4;
run; 

*********************************************************************************;
*T-TESTS;
*tests between sexes: SSBs= 1, 2, 3+;
**********************************************************************************;

PROC DESCRIPT data=WORK.NHANES design=wr;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
SUBPOPX incohadult=1;
* 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  weight;
* Class statement: specify categorical variable(s)*;
class sex  SSBn/nofreq;  
* Var statement: specify the analysis variable(s),and the levels of the analysis variable(s)*;
var  SSBn SSBn SSBn ; 
catlevel 1 2 3;
* Pairwise statement: Requests all pairwise contrasts*; 
pairwise 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.*;
*table SSBn;
SETENV LABWIDTH=35 COLWIDTH=10 ROWWIDTH=10 DECWIDTH=2 COLSPCE=2;
* Print statement: Prints the output*; 
print  nsum= "Sample Size" percent="percent" sepercent="SE" t_pct  p_pct="P-value"/nohead notime style=NCHS  
nsumfmt=f8.0  percentfmt=f8.1  sepercentfmt=f8.2 p_pctfmt=f8.4;
run; 


********************************************************************************************;
* Figure 2: Mean calories from SSBs, by sex and age, adults 20 years and older, US, 2011-14; 
*********************************************************************************************;

PROC DESCRIPT DATA=WORK.nhanes atlevel1=1 atlevel2=2 DESIGN=WR;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
subpopx incohadult=1;
* 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 weight;
* Subgroup statement: specify categorical variable(s), and the number of categories of the categorical variable(s)*;
SUBGROUP ageg3 SEX;
LEVELS 3 2 ;
* Var statement: specify the analysis variable(s)*;
VAR ssb_tkcal; 
* 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 SEX*ageg3;
* Print statement: Prints the output*; 
PRINT NSUM="SAMPLE SIZE" MEAN="MEAN" SEMEAN="SE"/ NOHEAD NOTIME STYLE=NCHS  NSUMFMT=F8.0 MEANFMT=F8.0 SEMEANFMT=F8.1;
rtitle "Mean total kcal from SSB by race and gender";
RTITLE "NHANES 2011-2014";
output nsum mean semean  atlev1 atlev2 /filename=SSBkcalAdult_ageg3_sex replace;
RUN;

data SSBkcalAdult_ageg3_sex_CI;
        set SSBkcalAdult_ageg3_sex;
        df=atlev2-atlev1; /*Calculate (DF) from the number of PSU (ATLEV2) minus the number of strata (ATLEV1).*/
        rsese=round((1/sqrt(df))*100,.01);
        rse=round((semean/mean)*100,.01);
        if rse=>40 then rseast="**";
        else if rse=>30 then rseast="*";
        if rsese=>40 then rseseast="**";
        else if rsese=>30 then rseseast="*";
        lowpercent=round(mean+tinv(0.025, df)*SEmean,0.01);
        uppercent=round(mean+tinv(0.975, df)*SEmean,0.01);
        drop ATLEV1 ATLEV2;       
run;

proc print  data=SSBkcalAdult_ageg3_sex_CI noobs label;
run; 
quit;


*********************************************************************************;
*T-TESTS;
*tests between sexes ;
*********************************************************************************;

PROC DESCRIPT data=WORK.NHANES design=wr;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
SUBPOPX incoh=1 ;
* 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  weight ;
* Class statement: specify categorical variable(s)*;
class sex ageg3 /nofreq; 
* Var statement: specify the analysis variable(s)*;
var ssb_tkcal;
* 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 ageg3;
* Pairwise statement: Requests all pairwise contrasts*;
pairwise sex;
SETENV LABWIDTH=35 COLWIDTH=10 ROWWIDTH=10 DECWIDTH=2 COLSPCE=2;
* Print statement: Prints the output*; 
Print nsum= "Sample Size"  mean="mean" semean="SE" t_mean p_mean="P-value"/nohead notime style=NCHS  
nsumfmt=f8.0 meanfmt=f8.1  semeanfmt=f8.2 p_meanfmt=f8.4;
run; 

*********************************************************************************;
*ANALYSIS OF TRENDS: by age group;
*********************************************************************************;

PROC DESCRIPT data=WORK.NHANES design=wr;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
SUBPOPX incoh=1 ;
* 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  weight ;
* Class statement: specify categorical variable(s)*;
class sex ageg3/nofreq; 
* Var statement: specify the analysis variable(s)*;
var ssb_tkcal ; 
* Polynomial statement: Requests linear and quadratic trends*;
polynomial ageg3=2;
* 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 sex;
SETENV LABWIDTH=45 COLWIDTH=10 ROWWIDTH=10 DECWIDTH=2 COLSPCE=2;
* Print statement: Prints the output*; 
print nsum= "Sample Size" mean="mean" semean="SE" t_mean p_mean="P-value"/nohead notime style=NCHS 
nsumfmt=f8.0 meanfmt=f8.1 semeanfmt=f8.1 p_meanfmt=f8.4;
run; 

*************************************************************************************************************;
* Figure 3: Mean calories from SSBs, by sex and race/Hispanic origin, adults 20 years and older, US, 2011-14;
*************************************************************************************************************;

PROC DESCRIPT DATA=WORK.nhanes atlevel1=1 atlevel2=2 DESIGN=WR;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
subpopx incohadult=1;
* 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 weight;
* Subgroup statement: specify categorical variable(s), and the number of categories of the categorical variable(s)*;
SUBGROUP race2 SEX;
LEVELS 5 2 ;
* Var statement: specify the analysis variable(s)*;
VAR ssb_tkcal ;
* 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 SEX*race2;
* Print statement: Prints the output*; 
PRINT NSUM="SAMPLE SIZE" MEAN="MEAN" SEMEAN="SE"/ NOHEAD NOTIME STYLE=NCHS  NSUMFMT=F8.0 MEANFMT=F8.0 SEMEANFMT=F8.1;
rtitle "Mean kcal SSB by race and gender";
RTITLE "NHANES 2009-2012";
output nsum mean semean  atlev1 atlev2 /filename=SSBkcalAdult_race_sex replace;
RUN;

data SSBkcalAdult_race_sex_CI;
        set SSBkcalAdult_race_sex;
        df=atlev2-atlev1; /*Calculate (DF) from the number of PSU (ATLEV2) minus the number of strata (ATLEV1).*/
        rsese=round((1/sqrt(df))*100,.01);
        rse=round((semean/mean)*100,.01);
        if rse=>40 then rseast="**";
        else if rse=>30 then rseast="*";
        if rsese=>40 then rseseast="**";
        else if rsese=>30 then rseseast="*";
        lowpercent=round(mean+tinv(0.025, df)*SEmean,0.01);
        uppercent=round(mean+tinv(0.975, df)*SEmean,0.01);
        drop ATLEV1 ATLEV2;       
run;

proc print  data=SSBkcalAdult_race_sex_CI noobs label;
run; 
quit;

*********************************************************************************;
*T-TESTS;
* test for differences by race and sex;
*********************************************************************************;

PROC DESCRIPT data=WORK.NHANES design=wr;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
SUBPOPX incoh=1;
* 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  weight ;
* Class statement: specify categorical variable(s)*;
class sex ageg3 race2/nofreq; 
* Var statement: specify the analysis variable(s)*;
var ssb_tkcal ;
* 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 ageg3 sex;
* Pairwise statement: Requests all pairwise contrasts*;
pairwise race2;
SETENV LABWIDTH=35 COLWIDTH=10 ROWWIDTH=10 DECWIDTH=2 COLSPCE=2;
* Print statement: Prints the output*; 
Print nsum= "Sample Size" mean="mean" semean="SE" t_mean p_mean="P-value"/nohead notime style=NCHS 
nsumfmt=f8.0 meanfmt=f8.1  semeanfmt=f8.2 p_meanfmt=f8.4;
run; 



PROC DESCRIPT data=WORK.NHANES design=wr;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
SUBPOPX incoh=1;
* 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  weight ;
* Class statement: specify categorical variable(s)*;
class sex ageg3 race2/nofreq; 
* Var statement: specify the analysis variable(s)*;
var  ssb_tkcal ;
* 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 ageg3 race2;
* Pairwise statement: Requests all pairwise contrasts*;
pairwise sex;
SETENV LABWIDTH=35 COLWIDTH=10 ROWWIDTH=10 DECWIDTH=2 COLSPCE=2;
* Print statement: Prints the output*; 
print  nsum= "Sample Size" mean="mean" semean="SE" t_mean p_mean="P-value"/nohead notime style=NCHS nsumfmt=f8.0 
meanfmt=f8.1  semeanfmt=f8.2 p_meanfmt=f8.4;
run; 


***********************************************************************************************;
**OBTAINING ESTIMATES FOR FIGURES 4 AND 5 USING THE MEAN RATIOS METHOD USED IN THE DATA BRIEF**;
***********************************************************************************************;

* Figure 4: SSB calories as a percent of total caloric intake, by sex and age, adults 20 years 
* and older, US, 2011-14; 
**************************************************************************************************************;

PROC DESCRIPT DATA=WORK.nhanes atlevel1=1 atlevel2=2 DESIGN=WR;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
subpopx incohadult=1;
* 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 weight;
* Subgroup statement: specify categorical variable(s), and the number of categories of the variables*;
SUBGROUP SEX ageg3;
LEVELS  2 3;
* Var statement: specify the analysis variable(s),uses precalculated percentage for each individual*;
VAR pE_SSB; 
* 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  sex*ageg3;
PRINT NSUM="SAMPLE SIZE" MEAN="MEAN" SEMEAN="SE"/ NOHEAD NOTIME STYLE=NCHS NSUMFMT=F8.0 MEANFMT=F8.1 SEMEANFMT=F8.2;
rtitle "Mean percent energy from SSB by gender and age";
RTITLE "NHANES 2011-2014";
output nsum mean semean  atlev1 atlev2 /filename=pE_SSBAdult_age_sex replace;
RUN;

data pE_SSBAdult_age_sex_CI;
        set pE_SSBAdult_age_sex;
        df=atlev2-atlev1; /*Calculate (DF) from the number of PSU (ATLEV2) minus the number of strata (ATLEV1).*/
        rsese=round((1/sqrt(df))*100,.01);
        rse=round((semean/mean)*100,.01);
        if rse=>40 then rseast="**";
        else if rse=>30 then rseast="*";
        if rsese=>40 then rseseast="**";
        else if rsese=>30 then rseseast="*";
        lowpercent=round(mean+tinv(0.025, df)*SEmean,0.01);
        uppercent=round(mean+tinv(0.975, df)*SEmean,0.01);
        drop ATLEV1 ATLEV2;       
run;

proc print  data=pE_SSBAdult_age_sex_CI noobs label;
run; 
quit;


*************************************;
*T-TESTS;
*************************************;

PROC DESCRIPT data=WORK.NHANES design=wr;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
SUBPOPX incohadult=1;
* 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  weight;
* Class statement: specify categorical variable(s)*;
class sex ageg3 race2/nofreq;   
* Var statement: specify the analysis variable(s),uses precalculated percentage for each individual*;
var pE_SSB;
* 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 ageg3;
* Pairwise statement: Requests all pairwise contrasts*; 
pairwise sex;
SETENV LABWIDTH=35 COLWIDTH=10 ROWWIDTH=10 DECWIDTH=2 COLSPCE=2;
print nsum= "Sample Size"  mean="mean" semean="SE" t_mean p_mean="P-value"/nohead notime style=NCHS  
nsumfmt=f8.0  meanfmt=f8.1  semeanfmt=f8.2  p_meanfmt=f8.4;
run; 

*************************************;
*ANALYSIS OF TRENDS;
*************************************;

PROC DESCRIPT data=WORK.NHANES design=wr;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
SUBPOPX incohadult=1;
* 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  weight;
* Class statement: specify categorical variable(s)*;
class sex ageg3/nofreq;   
* Var statement: specify the analysis variable(s),uses precalculated percentage for each individual*;
var  pE_SSB; 
* Polynomial statement: Requests linear and quadratic trends*; 
polynomial ageg3=2;
* 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 sex;
SETENV LABWIDTH=45 COLWIDTH=10 ROWWIDTH=10 DECWIDTH=2 COLSPCE=2;
print nsum= "Sample Size" mean="mean" semean="SE" t_mean p_mean="P-value"/nohead notime style=NCHS  nsumfmt=f8.0 
meanfmt=f8.1   semeanfmt=f8.1 p_meanfmt=f8.4;
run;


**************************************************************************************************************;
* Figure 5: SSB calories as a percent of total caloric intake, by sex and race/hispanic origin, adults 20 years 
* and older, US, 2011-14; 
**************************************************************************************************************;

PROC DESCRIPT DATA=WORK.nhanes atlevel1=1 atlevel2=2 DESIGN=WR;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
subpopx incohadult=1;
* 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 weight;
* Subgroup statement: specify categorical variable(s), and the number of categories of the variable(s)*;
SUBGROUP race2 SEX ageg3;
LEVELS 5 2 3;
* Var statement: specify the analysis variable(s),uses precalculated percentage for each individual*;
VAR pE_SSB; 
* 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 SEX*race2;
PRINT NSUM="SAMPLE SIZE" MEAN="MEAN" SEMEAN="SE"/ NOHEAD NOTIME STYLE=NCHS  NSUMFMT=F8.0 MEANFMT=F8.1 SEMEANFMT=F8.2;
rtitle "Mean percent energy from SSB by gender and race";
RTITLE "NHANES 2011-2014";
output nsum mean semean  atlev1 atlev2 /filename=pE_SSBAdult_race_sex replace;
RUN;

data pE_SSBAdult_race_sex_CI;
        set pE_SSBAdult_race_sex;
        df=atlev2-atlev1; /*Calculate (DF) from the number of PSU (ATLEV2) minus the number of strata (ATLEV1).*/
        rsese=round((1/sqrt(df))*100,.01);
        rse=round((semean/mean)*100,.01);
        if rse=>40 then rseast="**";
        else if rse=>30 then rseast="*";
        if rsese=>40 then rseseast="**";
        else if rsese=>30 then rseseast="*";
        lowpercent=round(mean+tinv(0.025, df)*SEmean,0.01);
        uppercent=round(mean+tinv(0.975, df)*SEmean,0.01);
        drop ATLEV1 ATLEV2;       
run;
proc print  data=pE_SSBAdult_race_sex_CI noobs label;
run; 
quit;


*************************************;
*T-TESTS;
*************************************;

PROC DESCRIPT data=WORK.NHANES design=wr;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
SUBPOPX incohadult=1;
* 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  weight;
* Class statement: specify categorical variable(s)*;
class sex ageg3 race2/nofreq; 
* Var statement: specify the analysis variable(s),uses precalculated percentage for each individual*;
var  pE_SSB;
* 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 sex;
* Pairwise statement: Requests all pairwise contrasts*; 
pairwise race2;
SETENV LABWIDTH=35 COLWIDTH=10 ROWWIDTH=10 DECWIDTH=2 COLSPCE=2;
print  nsum= "Sample Size" mean="mean" semean="SE" t_mean p_mean="P-value"/nohead notime style=NCHS  nsumfmt=f8.0 
meanfmt=f8.1   semeanfmt=f8.2 p_meanfmt=f8.4;
run;


PROC DESCRIPT data=WORK.NHANES design=wr;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
SUBPOPX incohadult=1;
* 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  weight;
* Class statement: specify categorical variable(s)*;
class sex ageg3 race2/nofreq;  
* Var statement: specify the analysis variable(s),uses precalculated percentage for each individual*;
var  pE_SSB ;
* 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 race2;
* Pairwise statement: Requests all pairwise contrasts*; 
pairwise sex;
SETENV LABWIDTH=35 COLWIDTH=10 ROWWIDTH=10 DECWIDTH=2 COLSPCE=2;
print nsum= "Sample Size" mean="mean" semean="SE" t_mean p_mean="P-value"/nohead notime style=NCHS  nsumfmt=f8.0  
meanfmt=f8.1   semeanfmt=f8.2 p_meanfmt=f8.4;
run;


**************************************************************************************************************;
        **ALTERNATIVE METHOD FOR OBTAINING ESTIMATES FOR FIGURES 4 AND 5 - RATIO OF MEANS METHOD**;
**************************************************************************************************************;

* Figure 4: Percentage of total daily kilocalories consumed from sugar-sweetened beverages on a given day*; 
*by adults aged 20 and over, by sex and age: United States, 2011–2014;                                        *; 
***************************************************************************************************************;

proc ratio data =nhanes  filetype = sas atlevel1=1 atlevel2=2 design=WR;
* Weight statement: specify appropriate weight, accounts for the unequal probability of sampling and non-response.*;
weight weight;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
subpopx incohadult=1  /NAME=" ";
* Nest statement: PSUs nested within Strata accounts for the design effects*;
nest SDMVSTRA SDMVPSU;
* Numer statement: specify the numerator analysis variable(s)*;
numer SSB_TKCAL;
* Denom statement: specify the denominator analysis variable(s), must match the number of numerator variable(s)*;
denom dr1tkcal;
* Subgroup statement: specify categorical variable(s), and number of categories of the variable(s)*;
SUBGROUP SEX ageg3;
LEVELS  2 3;
* 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  sex*ageg3;
SETENV LABWIDTH=20 COLSPCE=1 colwidth=10 decwidth=3;
* Print statement: Prints the output*; 
PRINT NSUM="Sample" WSUM="PopSize" WYSUM="Food" WXSUM="Total"
RHAT="Ratio" / NSUMFMT=F6.0 WSUMFMT=F9.0 RHATFMT=F6.5 STYLE=NCHS;
PRINT RHAT="Ratio" SERHAT="SE" LOWRHAT="Lower 95% Limit" UPRHAT="Upper 95% Limit"
/ STYLE=NCHS;
output nsum RHAT SERHAT LOWRHAT UPRHAT atlev1 atlev2/filename=PIR  replace;
run;

*************************************;
*T-TESTS;
*************************************;

proc ratio data =NHANES  filetype = sas design=WR;
* Weight statement: specify appropriate weight, accounts for the unequal probability of sampling and non-response.*;
weight weight;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
subpopx incohadult=1  /NAME=" ";
* Nest statement: PSUs nested within Strata accounts for the design effects*;
nest SDMVSTRA SDMVPSU;
* Numer statement: specify the numerator analysis variable(s)*;
numer SSB_TKCAL;
* Denom statement: specify the denominator analysis variable(s), must match the number of numerator variable(s)*;
denom dr1tkcal;
* Subgroup statement: specify categorical variable(s), and number of categories of the variable(s)*;
SUBGROUP SEX ageg3;
LEVELS  2 3;
* 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 ageg3;
* Pairwise statement: Requests all pairwise contrasts*; 
pairwise sex;
setenv colwidth = 1;
setenv decwidth = 4;
* Print statement: Prints the output*; 
PRINT NSUM="SamSize" RHAT="Estimate" SERHAT="SE" T_RHAT="T_Stat" P_RHAT="P-Value" / NSUMFMT=F8.0 
RHATFMT=F8.3 SERHATFMT=F7.4 T_RHATFMT=F7.2 P_RHATFMT=F9.5 STYLE=NCHS;
output nsum RHAT SERHAT T_RHAT P_RHAT /filename=Sex_Contrasts replace;
run;

*************************************;
*ANALYSIS OF TRENDS;
*************************************;
proc ratio data =NHANES  filetype = sas design=WR;
* Weight statement: specify appropriate weight, accounts for the unequal probability of sampling and non-response.*;
weight weight;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
subpopx incohadult=1  /NAME=" ";
* Nest statement: PSUs nested within Strata accounts for the design effects*;
nest SDMVSTRA SDMVPSU;
* Numer statement: specify the numerator analysis variable(s)*;
numer SSB_TKCAL;
* Denom statement: specify the denominator analysis variable(s), must match the number of numerator variable(s)*;
denom dr1tkcal;
* Subgroup statement: specify categorical variable(s), and number of categories of the variable(s)*;
SUBGROUP SEX ageg3;
LEVELS  2 3;
* 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 sex;
* Poly statement: Requests linear and quadratic trends*; 
POLY ageg3=2 / NAME="age TREND";
setenv colwidth = 1;
setenv decwidth = 4;
* Print statement: Prints the output*; 
PRINT NSUM="SamSize" RHAT="Estimate" SERHAT="SE" T_RHAT="T_Stat" P_RHAT="P-Value" / NSUMFMT=F8.0 
RHATFMT=F8.3 SERHATFMT=F7.4 T_RHATFMT=F7.2 P_RHATFMT=F9.5 STYLE=NCHS;
output nsum RHAT SERHAT T_RHAT P_RHAT /filename=Sex_Contrasts replace;
run;

**************************************************************************************************************;
* Figure 5: SSB calories as a percent of total caloric intake, by sex and race/hispanic origin, adults 20 years 
* and older, US, 2011-14; 
**************************************************************************************************************;

proc ratio data =nhanes  filetype = sas atlevel1=1 atlevel2=2 design=WR;
* Weight statement: specify appropriate weight, accounts for the unequal probability of sampling and non-response.*;
weight weight;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
subpopx incohadult=1  /NAME=" ";
* Nest statement: PSUs nested within Strata accounts for the design effects*;
nest SDMVSTRA SDMVPSU;
* Numer statement: specify the numerator analysis variable(s)*;
numer SSB_TKCAL;
* Denom statement: specify the denominator analysis variable(s), must match the number of numerator variable(s)*;
denom dr1tkcal;
* Subgroup statement: specify categorical variable(s), and number of categories of the variable(s)*;
SUBGROUP race2 SEX ageg3;
LEVELS 5 2 3;
* 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 SEX*race2 ;
SETENV LABWIDTH=20 COLSPCE=1 colwidth=10 decwidth=3;
* Print statement: Prints the output*; 
PRINT NSUM="Sample" WSUM="PopSize" WYSUM="Food" WXSUM="Total"
RHAT="Ratio" / NSUMFMT=F6.0 WSUMFMT=F9.0 RHATFMT=F6.5 STYLE=NCHS;
PRINT RHAT="Ratio" SERHAT="SE" LOWRHAT="Lower 95% Limit" UPRHAT="Upper 95% Limit"
/ STYLE=NCHS;
output nsum RHAT SERHAT LOWRHAT UPRHAT atlev1 atlev2/filename=PIR  replace;
run;

*************************************;
*T-TESTS;
*************************************;

proc ratio data =NHANES  filetype = sas design=WR;
* Weight statement: specify appropriate weight, accounts for the unequal probability of sampling and non-response.*;
weight weight;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
subpopx incohadult=1  /NAME=" ";
* Nest statement: PSUs nested within Strata accounts for the design effects*;
nest SDMVSTRA SDMVPSU;
* Numer statement: specify the numerator analysis variable(s)*;
numer SSB_TKCAL;
* Denom statement: specify the denominator analysis variable(s), must match the number of numerator variable(s)*;
denom dr1tkcal;
* Subgroup statement: specify categorical variable(s), and number of categories of the variable(s)*;
SUBGROUP SEX race2;
LEVELS  2 5;
* 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 sex;
* Pairwise statement: Requests all pairwise contrasts*; 
pairwise race2;
setenv colwidth = 1;
setenv decwidth = 4;
* Print statement: Prints the output*; 
PRINT NSUM="SamSize" RHAT="Estimate" SERHAT="SE" T_RHAT="T_Stat" P_RHAT="P-Value" / NSUMFMT=F8.0 
RHATFMT=F8.3 SERHATFMT=F7.4 T_RHATFMT=F7.2 P_RHATFMT=F9.5 STYLE=NCHS;
output nsum RHAT SERHAT T_RHAT P_RHAT /filename=Sex_Contrasts replace;
run;


proc ratio data =NHANES  filetype = sas design=WR;
* Weight statement: specify appropriate weight, accounts for the unequal probability of sampling and non-response.*;
weight weight;
* Subpopx statement: specify the subpopulation of interest (the inclusion criteria)*;
subpopx incohadult=1  /NAME=" ";
* Nest statement: PSUs nested within Strata accounts for the design effects*;
nest SDMVSTRA SDMVPSU;
* Numer statement: specify the numerator analysis variable(s)*;
numer SSB_TKCAL;
* Denom statement: specify the denominator analysis variable(s), must match the number of numerator variable(s)*;
denom dr1tkcal;
* Subgroup statement: specify categorical variable(s), and number of categories of the variable(s)*;
SUBGROUP SEX race2;
LEVELS  2 5;
* Table(s) statement: specifies cross-tabulations for which estimates are requested. If a table statement is not present, 
a one-dimensional distribution is generated for each variable in the subgroup statement.*;
table race2;
* Pairwise statement: Requests all pairwise contrasts*; 
pairwise sex;
setenv colwidth = 1;
setenv decwidth = 4;
* Print statement: Prints the output*; 
PRINT NSUM="SamSize" RHAT="Estimate" SERHAT="SE" T_RHAT="T_Stat" P_RHAT="P-Value" / NSUMFMT=F8.0 
RHATFMT=F8.3 SERHATFMT=F7.4 T_RHATFMT=F7.2 P_RHATFMT=F9.5 STYLE=NCHS;
output nsum RHAT SERHAT T_RHAT P_RHAT /filename=Sex_Contrasts replace;
run;