05-27-2016 10:43 AM - edited 05-27-2016 12:05 PM
I have a question about PROC POWER. For two groups, I can use the code below. However, for more than two groups (ex.
five groups), I am not sure what the option I need to add to it. Does anyone know?
proc power; twosamplefreq test=lrchi oddsratio = 2 refproportion = 0.4 npergroup = . power = 0.9; run;
05-27-2016 11:56 AM
Generally the options support multiple values and then proc power will calculated for combinations of the parameters provided.
proc power; twosamplefreq test=lrchi oddsratio = 2 2.5 refproportion = 0.4 0.5 npergroup = . power = 0.9 0.95; run;
provides 8 N per group results, 2 (oddsratio parameter values)*2(refproportion values)*2(power parameters).
You can provide different numbers of parameters for each
05-27-2016 04:15 PM
It seems to me that there is no such option available in PROC POWER (as of SAS/STAT version 14.1).
So, I see two possibilities to calculate the sample size for more than two groups in SAS:
With option 1 you have maximum flexibility, but option 2 requires less programming effort. The sample size formula I found in Sample Size Calculations in Clinical Research, Second Edition, p. 149, looks promising: I applied it to your two-group example (using macro Bisect found in this paper to solve for the non-centrality parameter involved) and obtained a result of 178.6..., i.e. n=179, while PROC POWER gave n=178.
Obviously, with five groups you have more degrees of freedom, so you would have to specify more input parameters (e.g. odds ratios) than in the two-group case.
05-27-2016 06:57 PM
05-30-2016 06:35 AM
Contrary to my initial perception, the features of the LOGISTIC statement of PROC POWER do not seem to exactly match your requirements.
But I think you can obtain the desired sample size using options 1 and 2 of my earlier post. Please find SAS code below for both options:
First, I determine the required sample size per group for a likelihood-ratio chi-square test to compare k=5 independent proportions with given reference proportion (0.4), odds ratios (3, 2.5, 1.8, 1.5), alpha (0.05) and power (0.90), using the formula mentioned in my earlier post. Then I perform a simulation (with 100,000 iterations) to determine the actual power of the test, given that sample size (run time: about 30 s on my workstation).
In the example, the result is n=87 per group and the estimated actual power is 0.9096 with an exact 95% CI of (0.9078, 0.9113). You can easily compute the power for other values of n by setting, say, %let n=85; and then re-running the simulation (i.e. fine-tune the sample size). As it turns out in this example, even with only n=85 per group the lower 95% confidence limit of the power estimate is 0.9000.
%let k=5; /* number of groups */ %let alpha=0.05; %let power=0.90; %let p0=0.4; /* reference proportion (in group 0) */ %let oddsr=3 2.5 1.8 1.5; /* list of odds ratios (group 1 vs. group 0, group 2 vs. group 0, ...) */ %let nsim=100000; /* number of iterations in simulation */ %let q=%sysfunc(quantile(chisq,%sysevalf(1-&alpha),&k-1)); %Bisect(cdf('chisq',&q,&k-1,x)=%sysevalf(1-&power),0,1E5,1E-12); /* macro source code from http://analytics.ncsu.edu/sesug/2007/PO21.pdf */ /* Sample size estimation based on formula from Chow, S.C., Shao, J., and Wang, H. (2008): Sample Size Calculations in Clinical Research, 2nd ed., Chapman & Hall/CRC Press/Taylor & Francis Group, New York, p. 149 */ data _null_; array r[0:%eval(&k-1)] _temporary_ (1 &oddsr); array p[0:1,0:%eval(&k-1)] _temporary_; do j=0 to &k-1; p[1,j]=1/(1+(1/&p0-1)/r[j]); P[0,j]=1-p[1,j]; end; do j=0 to &k-1; p0_+p[0,j]; p1_+p[1,j]; end; do j=0 to &k-1; s+(p[0,j]-p0_/&k)**2/p0_+(p[1,j]-p1_/&k)**2/p1_; end; n=ceil(&root/s/&k); call symputx('n',n); put 'n = ' n 'per group'; run; /* Simulation to determine the actual power, given the sample size determined above */ data sim; call streaminit(314159); array r[0:%eval(&k-1)] _temporary_ (1 &oddsr); array p[0:%eval(&k-1)] _temporary_; do j=0 to &k-1; p[j]=1/(1+(1/&p0-1)/r[j]); end; do m=1 to ≁ do g=0 to &k-1; y=1; f=rand('binom', p[g], &n); output; y=0; f=&n-f; output; end; end; run; ods listing close; ods output ChiSq=stats(keep=statistic prob where=(statistic=:'L')); proc freq data=sim; by m; weight f; tables y*g / chisq; run; ods listing; data lrchi; set stats(keep=prob); sig=(prob<&alpha); run; ods exclude BinomialTest; proc freq data=lrchi; tables sig / binomial(level='1'); run;