Annual Percentage Change with 95%CIs
inspired by
https://goo.gl/VPn5mo
https://communities.sas.com/t5/SAS-Procedures/Annual-Percentage-Change-with-95-CIs/m-p/329864
* Three techniques
1. proc freq exact (hypergeomethic?) 6.462% 6.478%
2. Normal approximation 6.461 6.479
3. Exact using an interative datastep solution 0.064622 0.064784
HAVE ( I restructure your data )
Up to 40 obs WORK.HAVE total obs=14
Obs YER EVENT TOT
1 2004 0 2313308 Have HIV
2 2004 1 33439457 Do not have AIDS
3 2005 0 2468735
4 2005 1 33516847
5 2006 0 2043291
6 2006 1 34203531
7 2007 0 1848984
8 2007 1 34703545
9 2008 0 1684509
10 2008 1 35171713
11 2009 0 1704744
12 2009 1 35372460
13 2010 0 1424794
14 2010 1 35893687
WANT
lower upper
Obs YER XL_BIN XU_BIN
Normal Approx see below Exact hand calculated
1 2004 6.462% 6.478% ( 6.461, 6.479 ) 0.064622 0.064784
2 2005 6.852% 6.869%
3 2006 5.630% 5.645%
4 2007 5.051% 5.066%
5 2008 4.564% 4.577%
6 2009 4.591% 4.605%
7 2010 3.812% 3.824%
SOLUTION
========
* I had to do some reformating;
data have;
input Yer pop pct lwr upr;
inc=int(pop*(lwr+upr)/200);
lwrupr=(lwr+upr)/2;
put lwrupr;
event=0;tot=inc;output;
event=1;tot=pop-inc;output;
keep yer event tot;
cards4;
2004 35752765 246 6.8805867 6.06 7.78
2005 35985582 262 7.2806937 6.44 8.2
2006 36246822 218 6.0143204 5.26 6.85
2007 36552529 198 5.4168619 4.7 6.21
2008 36856222 181 4.9109754 4.23 5.67
2009 37077204 183 4.9356473 4.26 5.69
2010 37318481 148 3.9658635 3.67 4.64
;;;;
run;quit;
proc freq data=have;
by yer;
tables event / list bin exact alpha=.05 out=cnt;
weight tot;
output out=rsk bin ;
run;quit;
proc print data=rsk;
format xl_bin xu_bin percent8.5;
var yer xl_bin xu_bin;
run;quit;
%macro utl_exbnci(n=,x=,p1=,p2=,alpha=);
%global &p1 &p2;
data _null_;
p1low=0.0;
p1mid=0.5;
p1high=1.0;
p2low=0.0;
p2mid=0.5;
p2high=1.0;
n=&n;
x=&x;
pval=&alpha/2;
if x=0 then p1mid=0;
else do until(p1high-p1low < 0.0001);
if 1-probbnml(p1mid,n,x-1) < pval then p1low=p1mid;
else p1high=p1mid;
p1mid=(p1low+p1high)/2;
end;
do until(p2high-p2low < 0.0001);
if probbnml(p2mid,n,x) < pval then p2high=p2mid;
else p2low=p2mid;
p2mid=(p2low+p2high)/2;
end;
call symput("&p1",put(100*p1mid,5.3));
call symput("&p2",put(100*p2mid,5.3));
run;
%mend utl_exbnci;
*----------*
| examples |
*----------*;
%utl_exbnci(n=35752765,x=2313308,p1=plow,p2=phigh,alpha=0.05)
data _null_;
put @5 "ci - ( &plow, &phigh )";
run;
/* T000101 DATASTEP INTERACTIVE METHOD FOR BINOMIOAL CONFIDENCE INTERVAL WORKS WITH 0 RESPONDERS */
%macro utl_bnmlcia(indata=, /* input data set name */
outdata=bnmlci /* output data set name */
);
data &outdata;
format n m best12. alpha p 8.6 l_limit u_limit 8.6 err_msg $55.;
set &indata;
if n=int(n) and m=int(m) and n>0 and m>=0 and m<=n and 0<alpha<1 then do;
p=m/n;
alpha2=alpha/2;
if m=0 then do;
l_limit=0;
link FIND_UL;
end;
else if m=n then do;
u_limit=1;
link FIND_LL;
end;
else do;
link FIND_LL;
link FIND_UL;
end;
end;
else do;
put n= m= alpha=;
if n ne int(n) or n<=0 then do;
err_msg='Variable n needs to be a positiv integer.';
put 'ERROR: ' err_msg;
end;
if m ne int(m) or m<0 then do;
err_msg='Variable m needs to be a nonnegative integer.';
put 'ERROR: ' err_msg;
end;
if m>n then do;
err_msg='Variable n needs to be equal or larger than variable m.';
put 'ERROR: ' err_msg;
end;
if not (0<alpha<1) then do;
err_msg='Variable alpha needs to be between 0 and 1.';
put 'ERROR: ' err_msg;
end;
put;
end;
RETURN;
FIND_LL: /* Bisection search for lower limit */
xu=p;
xl=0;
do until (crit < 0.000000001);
xm=(xu+xl)/2;
ptemp=probbnml(xm, n, m-1);
crit=abs(alpha2-(1-ptemp));
if ptemp<(1-alpha2) then do;
xl=xl;
xu=xm;
end;
else do;
xl=xm;
xu=xu;
end;
end;
l_limit=xm;
RETURN;
FIND_UL: /* Bisection search for upper limit */
xu=1;
xl=p;
do until (crit < 0.000000001);
xm=(xu+xl)/2;
ptemp=probbnml(xm, n, m);
crit=abs(alpha2-ptemp);
if ptemp<alpha2 then xu=xm;
else xl=xm;
end;
u_limit=xm;
RETURN;
label n = 'n'
m = 'm'
p = 'm/n'
alpha = 'Alpha Value'
l_limit = 'CI Lower Limit'
u_limit = 'CI Upper Limit';
drop alpha2 xu xl crit xm ptemp err_msg;
run;
%mend utl_bnmlcia;
****************** An Example **********************;
data t;
input n m alpha;
cards;
35752765 2313308 0.05
;
run;
%utl_bnmlcia(indata=t, outdata=t1);
proc print data=t1;
title1 "Binomial Distribution Confidence Intervals";
run;
Binomial Distribution Confidence Intervals
Obs N M ALPHA P L_LIMIT U_LIMIT
1 35752765 2313308 0.050000 0.064703 0.064622 0.064784
