/* edit this step to point to the location of the data */ proc import out=Adoptions datafile="C:\Users\\Documents\My SAS Files\Adoptions.csv" dbms=csv replace; getnames=yes; datarow=2; run; proc iml; use Adoptions; read all var {Events Trials}; close Adoptions; /* 1. Compute observed proportions */ Proportion = Events/ Trials; /* 2. compute overall proportion */ theta = sum(Events) / sum(Trials); Expected = theta * Trials; /* write theta to macro to use as reference line */ call symputx("AvgProp", theta); /* write results to data set */ results = Proportion || Expected; labels = {"Proportion" "Expected"}; create Stats from results[colname=labels]; append from results; close; /* 3. compute confidence limits for range of sample size */ /* Notice that these limits do not depend on the data, although we use theta = average proportion as a "target" */ /* plot limits at equally spaced points between min & max */ minN = min(Trials); maxN = max(Trials); n = T( do(minN, maxN, (maxN-minN)/20) ); n = round(n); /* binomial parameter must be integer */ p = {0.001 0.025 0.975 0.999}; /* lower/upper limits */ /* compute matrix with four columns, one for each CL */ limits = j(nrow(n), ncol(p)); do i = 1 to ncol(p); r = quantile("binom", p[i], theta, n); /* adjust quantiles according to method in Appendix A.1.1 Spiegelhalter (2005) p. 1197 */ numer = cdf("binom", r , theta, n) - p[i]; denom = cdf("binom", r , theta, n) - cdf("binom", r-1, theta, n); alpha = numer / denom; limits[,i] = (r-alpha)/n; end; /* write these control limits to data set */ results = n || limits; labels = {"N" "L3sd" "L2sd" "U2sd" "U3sd"}; create Limits from results[colname=labels]; append from results; close; quit; /* merge data with proportions */ data Adoptions; merge Adoptions Stats; run; /* append control limits */ data Funnel; set Adoptions Limits; run; /* 4. Plot proportions versus sample size. Overlay control limits */ proc sgplot data=Funnel; band x=N lower=L3sd upper=U3sd / nofill lineattrs=(color=lipk) /* fill fillattrs=GraphConfidence2 */ /* outline lineattrs=GraphPrediction */ legendlabel="99.8% limits" name="band99"; band x=N lower=L2sd upper=U2sd / nofill lineattrs=(color=gray) /* fill fillattrs=GraphConfidence*/ /* outline lineattrs=GraphPrediction */ legendlabel="95% limits" name="band95"; refline &AvgProp / axis=y; scatter x=Trials y=Proportion; keylegend "band95" "band99" / location=inside position=bottomright; yaxis grid values=(0.3 to 1 by 0.1); xaxis grid; run;