libname in1 'C:\Projects\Web Tutorial\variance'; 

/* Multcicd - program to calculated several 95 % CI's
each set of code is a stand alone module
    */

proc format ;
value statf
1='interviewed'
2='interviewed and examined'
;
value racethf
1='nh-white'
2='nh-black'
3='mex amer'
4='all others'
;

value hinsxf
1='priv/sing'
2='gov'
3='oth/unk'
4='none'
;
run;
                 
 
 
/*Variables on dataset in1.multci(keep=seqn sdmvstra sdmvpsu
wtmec6yr wtmec4yr wtmec2yr sddsrvyr
varpos varposx vartst 
ridageyr sex raceth2 hinsx); */
Data multci;
set in1.multci ;

if sddsrvyr in (1,2,3) ;

 
wtmec6yr=. ;
if sddsrvyr in (1,2) then wtmec6yr=(2/3)*wtmec4yr ;
else if sddsrvyr=3 then wtmec6yr=(1/3)*wtmec2yr ;

run;

data temp1 ;
set multci ;
run;

proc sort data=temp1 out=temps ;
by sdmvstra sdmvpsu ;
run;

PROC DESCRIPT DATA=TEmpS DESIGN=WR MEANS
ATLEVEL1=1 ATLEVEL2=2  deft2 noprint ;
NEST sdmvstra sdmvpsu/ MISSUNIT ;
WEIGHT wtmec6yr ;
VAR varposx ;
SUBPOPN ridageyr > 19  and ridageyr < 50 and vartst=100 ;
  SUBGROUP  raceth2  hinsx    ;
 LEVELS 4   4  ;
 TABLES raceth2*hinsx ;
  RTITLE "tutci - run to compare   calc 95% CI  ";
  OUTPUT   NSUM MEAN SEMEAN  atlev2 atlev1   /
FILENAME=temp
 FILETYPE=SAS REPLACE;
 
rformat raceth2 racethf.;
rformat hinsx hinsxf. ;
 run;

proc sort data=temp out=temps;
by raceth2 hinsx  ;
run;

 DATA tempa; SET temps;

  *usual WALD confidence interval;

* first calculate the nominal degrees of freedom=# of strata-# of PSUs
 and t-statistic ;
df=atlev2-atlev1;
tlow=tinv(.025,df);
tup=tinv(.975,df); 

ll=round((mean+tlow*semean),.1);
ul=round((mean+tup*semean),.1);
percent=round(mean,.1);
sepercent=round(semean,.1);

format   nsum    5.0;
format   df atlev1 atlev2 2.0 ;
format tlow 4.1 ;
format ll ul percent sepercent 6.1 ;
 

**  COLLETT EXACT Binomial / Clopper-Pearson CONF INT  ** ;
 
* calculate the nominal degrees of freedom=# of strata-# of PSUs
 and t-statistic although not required for this method ;

df=atlev2-atlev1;
tlow=tinv(.025,df);
tup=tinv(.975,df); 

P=MEAN/100 ;
IF MEAN>99.999 THEN NTOT=NSUM ;
ELSE IF P<1 THEN NTOT=(P*(1-P))/(SEMEAN/100)**2 ;
IF NTOT<1 THEN NTOT=1 ;
Y=P*NTOT ;
NDFL=2*(NTOT-Y+1) ;
DDFL=2*Y;
ALPHA=.975;
NDFU=2*(Y+1) ;
DDFU=2*(NTOT-Y) ;
FL=FINV(ALPHA,NDFL,DDFL) ;
FU=FINV(ALPHA,NDFU,DDFU) ;
PL=Y*(1/(Y+((NTOT-Y+1)*FL))) ;
PU=(Y+1)*(1/(Y+1+((NTOT-Y)/FU))) ;

LOWCI=round(PL*100,.1) ;
UPCI=round(PU*100,.1) ;

percent=round(mean,.1);
sepercent=round(semean,.1);

format   nsum    5.0;
format   df atlev1 atlev2 2.0 ;
format tlow 4.1 ;
format lowci upci percent sepercent 6.1 ;

* arcsine transformation;

* first calculate the nominal degrees of freedom=# of strata-# of PSUs
 and t-statistic ;
df=atlev2-atlev1;
tlow=tinv(.025,df);
tup=tinv(.975,df); 


mean=mean/100 ;
semean=semean/100;
meanz=arsin(sqrt(mean));
semean1=semean;
dy=1/(2*sqrt(mean-mean*mean));semeanz=dy*semean;
llar=meanz+tlow*semeanz;
ular=meanz+tup*semeanz;
lltrar=round(((sin(llar)**2)*100),.1);
ultrar=round(((sin(ular)**2)*100),.1);
mean=mean*100 ;	* convert back to percent ;
semean=semean*100 ;  

percent=round(mean,.1);
sepercent=round(semean,.1);

format   nsum    5.0;
format   df atlev1 atlev2 2.0 ;
format tlow 4.1 ;
format   lltrar ultrar percent sepercent 6.1 ;

 
* log transformation;

* first calculate the nominal degrees of freedom=# of strata-# of PSUs
 and t-statistic ;
df=atlev2-atlev1;
tlow=tinv(.025,df);
tup=tinv(.975,df); 

mean=mean/100 ; 
semean=semean/100 ;
logl=log(mean/(1-mean));selog=semean/(mean*(1-mean));
lllog=logl+tlow*selog;ullog=logl+tup*selog;
lltr=round(exp(lllog)/(1+exp(lllog))*100,.1);
ultr=round(exp(ullog)/(1+exp(ullog))*100,.1);


mean=mean*100 ;
semean=semean*100 ;  
percent=round(mean,.1);
sepercent=round(semean,.1);

format   nsum    5.0;
format   df atlev1 atlev2 2.0 ;
format tlow 4.1 ;
format   lltr ultr percent sepercent 6.1 ;
   

**********ASYMETRIC 95% CI USING Y=LOG X TRANSFORMATION*****;
  **********REF: KLEINBAUM, KUPPER, MORGANSTERN (1982, PGS 296-299)
** Mantel/Haenszel (sp?) method **;

* first calculate the nominal degrees of freedom=# of strata-# of PSUs
 and t-statistic ;
df=atlev2-atlev1;
tlow=tinv(.025,df);
tup=tinv(.975,df); 
 
 LSEMEAN=SEMEAN/MEAN;
 mhLCL95=round(MEAN*(EXP(tlow*LSEMEAN)),.1);
 mhUCL95=round(MEAN*(EXP(tup*LSEMEAN)),.1);
 if mean=0 then mhucl95=(3/nsum)*100 ;

percent=round(mean,.1);
sepercent=round(semean,.1);

format   nsum    5.0;
format   df atlev1 atlev2 2.0 ;
format tlow 4.1 ;
format   mhlcl95 mhucl95 percent sepercent 6.1 ;


proc sort ;
by   raceth2  hinsx ;

proc print noobs split='/' data=tempa; 
where  raceth2=2 ;
by raceth2 ;
var   hinsx percent sepercent nsum atlev2 atlev1 df tlow  

 
;
format raceth2 racethf.  hinsx hinsxf. 
;
label 
atlev1='#'/'PSU'
atlev2='#'/'Strat'
raceth2='Race'/'ethn'  
nsum='Num'/'varicella'/'status'
semean='Std'/'err'/'%'
df='Deg'/'of'/'Freed' 
tlow='t-stat'/'using'/'exact DF'
 
ll='Lower'/'95 %'/'WALD'/'exact DF'/'CI'
ul='Upper'/'95 %'/'WALD'/'exact DF'/'CI'
lowci='Lower'/'95 %'/'Clopper'/'Pearson'/'Exact'/'CI'
upci='Upper'/'95 %'/'Clopper'/'Pearson'/'Exact'/'CI'
lltrar='Lower'/'95 %'/'Arcsine'/'transformed'/'CI'
ultrar='Upper'/'95 %'/'Arcsine'/'transformed'/'CI'
lltr='Lower'/'95 %'/'Log'/'transformed'/'CI'
ultr='Upper'/'95 %'/'Log'/'transformed'/'CI'
mhlcl95='lower'/'95 %'/'MH'/'asymmetric'/'CI'
mhucl95='upper'/'95 %'/'MH'/'asymmetric'/'CI'
;
title1 'Multicicd-Percent varicella positive from NHANES 99-04';
title2 'by race/health insurance to all  95% CIs';
run;

proc print noobs split='/' data=tempa; 
where  raceth2=2 ;
by raceth2 ;
var   hinsx
ll ul lowci upci  lltrar ultrar lltr ultr  mhlcl95 mhucl95
 
;
format raceth2 racethf.  hinsx hinsxf. 
;
label 
atlev1='#'/'PSU'
atlev2='#'/'Strat'
raceth2='Race'/'ethn'  
nsum='Num'/'varicella'/'status'
semean='Std'/'err'/'%'
df='Deg'/'of'/'Freed' 
tlow='t-stat'/'using'/'exact DF'
 
ll='Lower'/'95 %'/'WALD'/'exact DF'/'CI'
ul='Upper'/'95 %'/'WALD'/'exact DF'/'CI'
lowci='Lower'/'95 %'/'Clopper'/'Pearson'/'Exact'/'CI'
upci='Upper'/'95 %'/'Clopper'/'Pearson'/'Exact'/'CI'
lltrar='Lower'/'95 %'/'Arcsine'/'trans'/'CI'
ultrar='Upper'/'95 %'/'Arcsine'/'trans'/'CI'
lltr='Lower'/'95 %'/'Log'/'trans'/'CI'
ultr='Upper'/'95 %'/'Log'/'trans'/'CI'
mhlcl95='lower'/'95 %'/'MH'/'asym'/'CI'
mhucl95='upper'/'95 %'/'MH'/'asym'/'CI'
;
title1 'Multicicd-Percent varicella positive from NHANES 99-04';
title2 'by race/health insurance to all  95% CIs';
run;


 DATA tempb; SET temps;

meanneg=100-mean ;

  *usual WALD confidence interval;

* first calculate the nominal degrees of freedom=# of strata-# of PSUs
 and t-statistic ;
df=atlev2-atlev1;
tlow=tinv(.025,df);
tup=tinv(.975,df); 

ll=round((meanneg+tlow*semean),.1);
ul=round((meanneg+tup*semean),.1);
percneg=round(meanneg,.1);
sepercent=round(semean,.1);

format   nsum    5.0;
format   df atlev1 atlev2 2.0 ;
format tlow 4.1 ;
format ll ul percneg sepercent 6.1 ;
 

**  COLLETT EXACT Binomial / Clopper-Pearson CONF INT  ** ;
 
* calculate the nominal degrees of freedom=# of strata-# of PSUs
 and t-statistic although not required for this method ;

df=atlev2-atlev1;
tlow=tinv(.025,df);
tup=tinv(.975,df); 

P=meanneg/100 ;
IF meanneg>99.999 THEN NTOT=NSUM ;
ELSE IF P<1 THEN NTOT=(P*(1-P))/(semean/100)**2 ;
IF NTOT<1 THEN NTOT=1 ;
Y=P*NTOT ;
NDFL=2*(NTOT-Y+1) ;
DDFL=2*Y;
ALPHA=.975;
NDFU=2*(Y+1) ;
DDFU=2*(NTOT-Y) ;
FL=FINV(ALPHA,NDFL,DDFL) ;
FU=FINV(ALPHA,NDFU,DDFU) ;
PL=Y*(1/(Y+((NTOT-Y+1)*FL))) ;
PU=(Y+1)*(1/(Y+1+((NTOT-Y)/FU))) ;

LOWCI=round(PL*100,.1) ;
UPCI=round(PU*100,.1) ;

percneg=round(meanneg,.1);
sepercent=round(semean,.1);

format   nsum    5.0;
format   df atlev1 atlev2 2.0 ;
format tlow 4.1 ;
format lowci upci percneg sepercent 6.1 ;

* arcsine transformation;

* first calculate the nominal degrees of freedom=# of strata-# of PSUs
 and t-statistic ;
df=atlev2-atlev1;
tlow=tinv(.025,df);
tup=tinv(.975,df); 


meanneg=meanneg/100 ;
semean=semean/100;
meannegz=arsin(sqrt(meanneg));
semean1=semean;
dy=1/(2*sqrt(meanneg-meanneg*meanneg));semeanz=dy*semean;
llar=meannegz+tlow*semeanz;
ular=meannegz+tup*semeanz;
lltrar=round(((sin(llar)**2)*100),.1);
ultrar=round(((sin(ular)**2)*100),.1);
meanneg=meanneg*100 ;	* convert back to percneg ;
semean=semean*100 ;  

percneg=round(meanneg,.1);
sepercent=round(semean,.1);

format   nsum    5.0;
format   df atlev1 atlev2 2.0 ;
format tlow 4.1 ;
format   lltrar ultrar percneg sepercent 6.1 ;

 
* log transformation;

* first calculate the nominal degrees of freedom=# of strata-# of PSUs
 and t-statistic ;
df=atlev2-atlev1;
tlow=tinv(.025,df);
tup=tinv(.975,df); 

meanneg=meanneg/100 ; 
semean=semean/100 ;
logl=log(meanneg/(1-meanneg));selog=semean/(meanneg*(1-meanneg));
lllog=logl+tlow*selog;ullog=logl+tup*selog;
lltr=round(exp(lllog)/(1+exp(lllog))*100,.1);
ultr=round(exp(ullog)/(1+exp(ullog))*100,.1);


meanneg=meanneg*100 ;
semean=semean*100 ;  
percneg=round(meanneg,.1);
sepercent=round(semean,.1);

format   nsum    5.0;
format   df atlev1 atlev2 2.0 ;
format tlow 4.1 ;
format   lltr ultr  percneg sepercent 6.1 ;
   

**********ASYMETRIC 95% CI USING Y=LOG X TRANSFORMATION*****;
  **********REF: KLEINBAUM, KUPPER, MORGANSTERN (1982, PGS 296-299)
** Mantel/Haenszel (sp?) method **;

* first calculate the nominal degrees of freedom=# of strata-# of PSUs
 and t-statistic ;
df=atlev2-atlev1;
tlow=tinv(.025,df);
tup=tinv(.975,df); 
 
 Lsemean=semean/meanneg;
 mhLCL95=round(meanneg*(EXP(tlow*Lsemean)),.1);
 mhUCL95=round(meanneg*(EXP(tup*Lsemean)),.1);
 if meanneg=0 then mhucl95=(3/nsum)*100 ;

percneg=round(meanneg,.1);
sepercent=round(semean,.1);

format   nsum    5.0;
format   df atlev1 atlev2 2.0 ;
format tlow 4.1 ;
format   mhlcl95 mhucl95 percneg sepercent 6.1 ;


proc sort ;
by   raceth2  hinsx ;

proc print noobs split='/' data=tempb ; 
where  raceth2=2 ;
by raceth2 ;
var   hinsx percneg sepercent nsum atlev2 atlev1 df tlow  
;
format raceth2 racethf.  hinsx hinsxf. 
;
label 
atlev1='# of'/'PSUs  '
atlev2='# of'/'Strata  '
raceth2='Race'/'Ethnicity'  
nsum='Number with'/'varicella'/'status'
semean='Standard'/'error'/'of the percent'
percneg='Percent'/'Varicella'/'Negative'
sepercent='Standard'/'error'/'of the'/'Percent'
df='Degrees'/'of'/'Freedom' 
tlow='t-stat'/'using'/'exact degrees of'/'freedom'
ll='Lower'/'95 %'/'WALD'/'exact DF'/'CI'
ul='Upper'/'95 %'/'WALD'/'exact DF'/'CI'
lowci='Lower'/'95 %'/'Clopper'/'Pearson'/'Exact'/'CI'
upci='Upper'/'95 %'/'Clopper'/'Pearson'/'Exact'/'CI'
lltrar='Lower'/'95 %'/'Arcsine'/'transformed'/'CI'
ultrar='Upper'/'95 %'/'Arcsine'/'transformed'/'CI'
lltr='Lower'/'95 %'/'Log'/'transformed'/'CI'
ultr='Upper'/'95 %'/'Log'/'transformed'/'CI'
mhlcl95='lower'/'95 %'/'MH'/'asymmetric'/'CI'
mhucl95='upper'/'95 %'/'MH'/'asymmetric'/'CI'
;
title1 'Multicicd-Percent varicella negative from NHANES 99-04';
title2 'by race/health insurance to all  95% CIs';
run;

proc print noobs split='/' data=tempb ; 
where  raceth2=2 ;
by raceth2 ;
var   hinsx  
ll ul lowci upci  lltrar ultrar lltr ultr  mhlcl95 mhucl95
 
;
format raceth2 racethf.  hinsx hinsxf. 
;
label 
atlev1='# of'/'PSUs  '
atlev2='# of'/'Strata  '
raceth2='Race'/'Ethnicity'  
nsum='Number with'/'varicella'/'status'
semean='Standard'/'error'/'of the percent'
percneg='Percent'/'Varicella'/'Negative'
sepercent='Standard'/'error'/'of the'/'Percent'
df='Degrees'/'of'/'Freedom' 
tlow='t-stat'/'using'/'exact degrees of'/'freedom'
ll='Lower'/'95 %'/'WALD'/'exact DF'/'CI'
ul='Upper'/'95 %'/'WALD'/'exact DF'/'CI'
lowci='Lower'/'95 %'/'Clopper'/'Pearson'/'Exact'/'CI'
upci='Upper'/'95 %'/'Clopper'/'Pearson'/'Exact'/'CI'
lltrar='Lower'/'95 %'/'Arcsine'/'trans'/'CI'
ultrar='Upper'/'95 %'/'Arcsine'/'trans'/'CI'
lltr='Lower'/'95 %'/'Log'/'trans'/'CI'
ultr='Upper'/'95 %'/'Log'/'trans'/'CI'
mhlcl95='lower'/'95 %'/'MH'/'asym'/'CI'
mhucl95='upper'/'95 %'/'MH'/'asym'/'CI'
;
title1 'Multicicd-Percent varicella negative from NHANES 99-04';
title2 'by race/health insurance to all  95% CIs';
run;