/* Simulation to determine coverage probabilities of the following three confidence regions for
the parameter vector (p1, p2, p3) of a Mult(100, p1, p2, p3) (multinomial) distribution with
p1=0.1, p2=0.18, p3=0.72 (example parameter values taken from SAS Usage Note 32609)
(1) 95% Wald confidence intervals for the individual probabilities (as computed by PROC CATMOD
in SAS Usage Note 32609, but the same as available from PROC FREQ except for the [trivial]
limitation of confidence bounds to [0, 1] in PROC FREQ)
(2) 98.333...% Wald confidence intervals as above (Bonferroni adjustment: (1-0.05/3)*100%), as
also suggested in R.G. Miller, Simultaneous Statistical Inference, p. 216, formula 10
(3) Bonferroni-adjusted CIs suggested in formula 14, p. 217, ibid., based on chi-square statistics
*/
%let n=100; /* parameters of multinomial distribution */
%let p1=0.1;
%let p2=0.18;
%let p3=0.72;
%let c=3; /* number of categories */
%let a=0.05; /* alpha */
data sim;
call streaminit(271828);
z=probit(1-&a/2); /* normal quantile without Bonferroni adjustment */
b=probit(1-&a/(2*&c)); /* normal quantile with Bonferroni adjustment */
x=cinv(1-&a/&c,1); /* chi-square quantile, df=1, with Bonferroni adjustment */
array n[3]; /* samples of the multinomial distribution */
array lcl[3,3]; /* lower confidence limits (indices e.g. [2,3] for interval 2, param. p3) */
array ucl[3,3]; /* upper confidence limits (indices as above) */
array y[3]; /* 0-1 indicators for "(p1,p2,p3) contained in conf. region (1, 2, 3)" */
do i=1 to 1e7; /* sample size for the simulation */
n1=0; n2=0; n3=0; /* initialization */
do j=1 to &n; /* generate sample from multinomial distribution */
v=rand('table', &p1, &p2, &p3);
n[v]+1;
end;
do k=1 to 3; /* calculate CIs for p1, p2, p3 */
p=n[k]/&n; /* point estimate of probability "pk" (k=1, 2, 3) */
s=sqrt(p*(1-p)/&n); /* approx. standard error of the above point estimate */
lcl[1,k]=p-z*s; /* CI 1 */
ucl[1,k]=p+z*s;
lcl[2,k]=p-b*s; /* CI 2 */
ucl[2,k]=p+b*s;
lcl[3,k]=(x+2*n[k]-sqrt(x*(x+4*n[k]*(&n-n[k])/&n)))/(2*(&n+x)); /* CI 3 */
ucl[3,k]=(x+2*n[k]+sqrt(x*(x+4*n[k]*(&n-n[k])/&n)))/(2*(&n+x));
do m=1 to 3; /* check if confidence regions cover the true parameter vector */
y[m]=(lcl[m,1] <= &p1 <= ucl[m,1] &
lcl[m,2] <= &p2 <= ucl[m,2] &
lcl[m,3] <= &p3 <= ucl[m,3]);
end;
end;
output;
end;
keep y:;
run;
/* Estimate coverage probabilities of the three confidence regions */
ods select BinomialCLs;
proc freq data=sim;
tables y: / binomial(level='1' cl=exact);
run;
/* Results:
Point est. Clopper-Pearson 95% CI
(1) 0.8551 0.8549 0.8554
(2) 0.9402 0.9401 0.9404
(3) 0.9633 0.9632 0.9634
*/
(Actually, this could be coded more elegantly by restricting the calculations to the finitely many [<<1E7] possible outcomes of the particular multinomial distribution and weighing the results by the corresponding probabilities. But the run time of the above code was only less than 1 minute.)
So, for the particular multinomial distribution Mult(100, 0.1, 0.18, 0.72), only the third confidence region achieves the nominal 95% confidence level. The second comes close, but the first fails, as expected.
Conclusion: The PROC CATMOD approach presented in SAS Usage Note 32609 does not seem to provide substantial benefits over using PROC FREQ (here with Bonferroni adjustment):
ods exclude BinomialTest;
proc freq data=a;
weight count / zeros;
tables y / binomial (level='1') alpha=0.01666666666667;
tables y / binomial (level='2') alpha=0.01666666666667;
tables y / binomial (level='3') alpha=0.01666666666667;
run;
There are alternative confidence regions in the literature which can be implemented easily in a data step.
