/
fpower.sas
201 lines (169 loc) · 6.73 KB
/
fpower.sas
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
/*-------------------------------------------------------------------*
* Name: fpower.sas *
* Title: Power computations for ANOVA designs *
Doc: http://www.datavis.ca/sasmac/fpower.html
* *
*-------------------------------------------------------------------*
* Author: Michael Friendly <FRIENDLY@YorkU.CA> *
* After : AJL Cary, Syntex Research; from: SAS SUGI, 1983 *
* Created: 14 Mar 1990 11:16:28 *
* Revised: 28 Aug 2013 09:05:51 *
* Version: 1.2-1 *
* - Default N extended to 90
* *
*-------------------------------------------------------------------*/
/*=
=Description:
The FPOWER macro computes power of an F-test for main effects or
interactions for between-subject effects in factorial designs with
or without repeated measures, assuming main effects are fixed.
The alternative used is a minimum power alternative. The program
can be used for ANY fixed effect in ANY crossed factorial design,
by designating the levels of the effect of interest as A, and the
levels of all other crossed factors as B.
It does *not* (correctly) handle power for repeated-measure factors,
split-plot designs, etc, for which power analysis is much more
complex.
The "effect size" is specified by DELTA. Assuming that only
the extreme factor levels have NON-ZERO effects, and
T1 = GM - DELTA/2
Tk = GM + DELTA/2
where DELTA is specified in units of SIGMA = SQRT(MSE)
The intraclass correlation (RHO) is assumed to be positive
and constant across all repeated measures.
The macro variable N provides:
DO-list of sample sizes to calculate power for.
An output dataset named PWRTABLE is created for plotting or
saving, and it contains an observation for each entry.
=*/
options nonotes;
*-------------- Define output format for small & -----------
--------------- and large probabilities -------------------;
Proc format;
picture PROBS /* .99901 - 1.0000 =' >.999'
0.0000 - 0.00089 =' <.001' */
other =' 9.999';
%macro fpower(
A=, /* Levels of factor for power calculation */
B=1, /* Levels of crossed factor B (default=1) */
C=1, /* Levels of crossed factor C (default=1) */
/* For >3 factors, make C=product of # of */
/* levels of factors D, E, etc. */
R=1, /* Levels of repeated measure factor */
ALPHA=.05, /* Significance level of test of factor A */
N =%str( 2 to 10 by 1, 12 to 18 by 2, 20 to 40 by 5, 50 to 90 by 10),
DELTA=.50 to 2.5 by 0.25, /* DO list of DELTA values */
RHO=0, /* Intraclass correlation for repeated */
/* measures (DO list of values) */
PTABLE=YES, /* Print a power table? */
PLOT=NO, /* Plot power*delta=N ? */
NTABLE=NO, /* Print a sample-size table ? */
OUT=PWRTABLE
);
%if &A=%str() %then %do;
%put ERROR : You must specify a value for A= in the macro call;
%goto DONE;
%end;
%let pfile = print;
%if %upcase(&PTABLE) ^=YES %then %do;
%if &sysscp = WIN
%then %str(filename ptable DUMMY 'nul:';);
%else %str(filename ptable DUMMY '/dev/null';);
%let pfile = ptable;
%end;
data &out;;
length A B C R N TRTDF 3;
drop fcrit ncol nalpha;
A = &A;
B = &B;
C = &C;
R = &R;
do ALPHA = &ALPHA;
nalpha+1;
file &pfile;
TRTDF = (A-1); *-- Treatment degrees of freedom--;
*-- Iterate thru RHO values supplied in RHOs --;
* If not repeated measures design, stop after ;
* one iteration. ;
do RHO = &RHO; *-- Get macro values;
put _PAGE_;
put @5 A 2. 'x' B 2. @;
if C>1 then put 'x' C 2. @;
put +2 ' layout Ha: T1=GM-Delta/2, T2=T3=...=T(k-1)'
'=GM, Tk=GM+Delta/2' +3 'tested at Alpha=' ALPHA 6.3;
if R>1 then put @10 'with' R 2.0
' Repeated measures and intraclass RHO = ' RHO 4.2 ;
else put ;
put @20 'DELTA (in units of sigma=Std. Dev.)'// @7 'N' @;
do DELTA = &DELTA; *-- Get macro values;
put DELTA 7.3 @;
end;
put;
ncol = 0;
do DELTA = &DELTA;
ncol = ncol+1;
put '-------' @;
end;
put '-------';
*-- Iterate through sample sizes --;
do N = &N; *-- Get macro values;
put / N 7.0 @;
*-- Compute error degrees of freedom & F crit --;
errdf = (a * b * c * (n-1));
fcrit = finv (1-alpha, trtdf, errdf, 0);
*-- Iterate through values of DELTA --;
do DELTA = &DELTA;
*-- Compute non-centrality parameter, & power --;
NC = N * B * C * R * DELTA**2/ (2*(1+(R-1)*RHO));
if NC > 140 then POWER = 1;
else POWER = 1 - probf(FCRIT, TRTDF, ERRDF, NC);
put POWER PROBS. @;
output;
end; *-- Do DELTA;
end; *-- Do N;
*-- Exit loop if not repeated measures;
if R=1 then STOP;
end; *-- do RHO;
end; *-- do ALPHA;
call symput(nalpha,put(nalpha,2.));
run;
%if %upcase(&PLOT)=YES %then %do;
proc plot data=&out;
plot POWER * DELTA = N /vaxis = 0 to 1 by .1 vref=1;
title2 'Power as a function of Effect Size (DELTA) and N';
run;
%end;
%if %upcase(&NTABLE)=YES %then %do;
data _ntable_;
set &out;
drop a b c r rho errdf nc trtdf;
where (power>=.5);
power = min(power,.9999);
* pwr = round(power,.05);
pwr = .1 * floor(10*power);
proc sort;
by alpha pwr delta n;
*-- Find the smallest n for given power, delta;
data _ntable_;
set _ntable_;
by alpha pwr delta n;
if first.delta then output;
/*
proc print data=_ntable_;
id pwr;
by pwr;
var delta alpha n power;
*/
proc tabulate;
class pwr delta;
var n;
table pwr='Power',
delta*n=' '*f=5.*sum=' ' / rts=9;
format n 3.;
title2 "Sample size to achive a given power, alpha= &alpha";
run;
%end;
%DONE:
options notes;
%mend FPOWER;
*------------------ End of FPOWER -----------------------;