FILENAME goptions URL "http://192.38.117.59/~linearpredictors/datafiles/goptions.sas";
%include goptions;

FILENAME url URL "http://192.38.117.59/~linearpredictors/datafiles/readVitaminD.sas";
%include url;

/*
Women from Ireland 
*/ 

DATA irishwomen;
	SET vitamind (WHERE = (country=4 and category=2));
    if bmi > 0 then bmigroup = 0;
    if bmi >= 25 then bmigroup = 1;
RUN;

PROC MEANS N MEAN STDERR DATA=irishwomen NWAY NOPRINT;
    CLASS bmigroup;
    VAR vitd;
    OUTPUT OUT=summary N=count MEAN=mvitd STDDEV=sdvitd STDERR=sdmvitd;
RUN;

DATA summary;
SET summary;

ci_lower = mvitd - probit(0.975)*sdmvitd;
ci_upper = mvitd + probit(0.975)*sdmvitd;
RUN;

title 'individual confidence limits: (ci_lower, ci_upper)';
proc print data=summary;
    var bmigroup count mvitd sdmvitd ci_lower ci_upper;
run;

DATA calculate;
MERGE summary (where=(bmigroup=0) rename=(count=n1 mvitd=mvitd1 sdvitd=sdvitd1))
      summary (where=(bmigroup=1) rename=(count=n2 mvitd=mvitd2 sdvitd=sdvitd2))
    ;

mean_diff=mvitd2-mvitd1;
sd_pool=sqrt(((n1-1)*sdvitd1**2+(n2-1)*sdvitd2**2)/(n1+n2-2));
sd_mdiff=sd_pool*sqrt(1/n1+1/n2);
ci_lower = mean_diff - probit(0.975)*sd_mdiff;
ci_upper = mean_diff + probit(0.975)*sd_mdiff;
RUN;

title 'confidence limit for mean difference: (ci_lower, ci_upper)';
proc print data=calculate;
    var _FREQ_ mean_diff sd_mdiff ci_lower ci_upper;
run;

/* genmod can make the normality approximation to the t */
/* but cannot do reml */ 
/*
PROC GENMOD DATA=irishwomen;
    CLASS bmigroup;
    MODEL vitd=bmigroup / dist=normal link=identity;
RUN;
*/
    
/* ttest will always make a CI based on the t-distribution */
/* or an approximation, if the SD's differ between groups */    
/*
PROC TTEST DATA=irishwomen;
    CLASS bmigroup;
    VAR vitd;
RUN;
*/