/*Test of proportions for PSU student drinking example*/
/*To test H0:p=0.5 for high-risk drinkers */

options ls=90 nocenter nodate;


/************************binomial likelihood function plot******************************/

data likelhd;
do p=0.0 to 0.99 by 0.001;
like=PDF('BINOMIAL',630,p,1315);
keep p like;
output;
end;
run;

proc print;run;

goption reset=all colors=(black);
symbol1 i=spline line=1;
proc gplot data=likelhd;
plot like*p   ;
run;

/************************binomial log-likelihood function plot****************************/

data loglike;
do x=0.01 to 0.99 by 0.001;
y=lfact(1315)-lfact(630)-lfact(685)+630*log(x)+685*log(1-x);   *the log-likelihood function;
output;
end;
run;

data cutoff;
cutoff=lfact(1315)-lfact(630)-lfact(685)+630*log(0.48)+685*log(1-0.48)-1.92;
run;

proc print data=loglike;run;    *cutoff=-5.738;

data aa;
set loglike;
if round(y,0.1) in (-5.7, -5.8);
run;
proc print;run;  *get whereabout the cutoff line intersects the loglikelihood;


goption reset=all colors=(black);
symbol1 i=spline line=1;
proc gplot data=loglike;
plot y*x / vref=-5.738 href=0.452 0.506 ; /*add vertical or horizontal lines to the plot*/
run;

quit;



/*******************************MLE*****************************************************/

proc sort data=likelhd out=srtlikelhd;
by  descending like;
run;

proc print data=srtlikelhd (obs=1);run;   *first record is the maximum;



/******************************frequency procedure***************************************/

DATA drinking;
        INPUT drinking $ count;
        DATALINES;
        high-risk 630
        low-risk 685
        ;
RUN;

PROC FREQ DATA=drinking ORDER=data;
        WEIGHT count;
        table drinking/binomial(p=0.5);
RUN;


data one;
phat=630/(630+685);
n=630+685;
ASE=sqrt(phat*(1-phat)/n);
Z=(phat-0.5)/ASE;
run;
proc print;run;