how to plot confidence interval

[yuwentaiwan]

克洛珀·皮爾森（Clopper Pearson），

這是n = 19，r = 3。

我該如何繪製這張照片。

謝謝！！

data propci;

if sided = 1 then do;
zalpha=probit(1-(alpha));
end;
else if sided = 2 then do;
zalpha=probit(1-(alpha/2));
end;

** Exact Clopper Pearson;
x = round (n*p,0.1);

* Calculate the lower limit.;

v1 = 2*(n-x+1);
v2 = 2*x;

if sided = 1 then do;
a = 1-(alpha);
end;
else if sided = 2 then do;
a = 1-(alpha/2); 
end;

coef = (n-x+1)/x;
fscore = finv(a,v1,v2);
exact_lcl = 1/(1+coef*fscore);

* Calculate the upper limit.;

v11 = 2*(x+1);
v22 = 2*(n-x);
fscore2 = finv(a,v11,v22);
coef2 = (x+1)/(n-x);
numer = coef2*fscore2;
exact_ucl = numer/(1+numer);
run;

[FreelanceReinh · Accepted solution]

Hello @yuwentaiwan,

Thanks for posting the nice graph.

This is my reproduction with SAS. I've added a label to the x-axis.

Here's the program which created it:

/* Set parameters */

%let n=19; /* number of Bernoulli trials */
%let ph=3; /* observed number of successes ("p hat") */
%let alpha=0.05;

/* Compute probabilities */

data prob;
do _n_=0 to 1000;
  p=_n_/1000;
  p1=cdf('binom',&ph,p,&n);
  p2=cdf('binom',&ph-1,p,&n);
  output;
end;
run;

/* Compute Clopper-Pearson confidence interval */

data fdat;
c=1; n=&ph;    output;
c=2; n=&n-&ph; output;
run;

ods select none;
ods output binomial=bcl;
proc freq data=fdat;
weight n;
tables c / binomial alpha=α
run;
ods select all;

data _null_;
set bcl;
if name1=:'XL' then call symputx('lcl',put(nvalue1,6.4));
if name1=:'XU' then call symputx('ucl',put(nvalue1,9.7));
run;

/* Create the graph */

ods html style=journal2;
title bold "h(p) = P[X <= &ph] or h(p) = P[X <= %eval(&ph-1)] for X ~ Bin(&n, p)";

proc sgplot data=prob;
series x=p y=p1 / legendlabel="P[X<=&ph]";
series x=p y=p2 / legendlabel="P[X<=%eval(&ph-1)]" lineattrs=(pattern=2);
xaxis offsetmin=0.04 offsetmax=0.04;
yaxis offsetmin=0.04 offsetmax=0.04 label='h(p)';
dropline x=&lcl y=%sysevalf(1-&alpha/2) / dropto=both label="&lcl" lineattrs=(color=black);
dropline x=&ucl y=%sysevalf(&alpha/2)   / dropto=both label="&ucl" lineattrs=(color=black);
keylegend / location=inside position=top across=1 linelength=20;
run;

As you can see, the parameters n, ph and alpha are macro variable values. For certain parameter combinations the placement, e.g., of the legend should be modified in order to avoid collisions. The above code implements the two-sided Clopper-Pearson confidence interval (I used PROC FREQ rather than a DATA step to compute it). Please let me know if you need help with the one-sided case.

[yuwentaiwan]

Thank you!!:)