%let n1 = 10; /* group sizes*/
%let n2 = 10;
%let NumSamples = 5000; /* number of simulated samples */
data PowerSim(drop=i);
call streaminit(321);
do Delta = 0 to 2 by 0.1;
do SampleID = 1 to &NumSamples;
c = 1; /* Group 1 */
do i = 1 to &n1;
x = rand("Normal", 0, 1); /* x ~ N(0,1) */
output;
end;
c = 2; /* Group 2 */
do i = 1 to &n2;
x = rand("Normal", Delta, 1); /* x ~ N(Delta, 1) */
output;
end;
end;
end;
run;
/* 2. Compute (pooled) t test for each sample */
%ODSOff
proc ttest data=PowerSim;
by Delta SampleID;
class c;
var x;
ods output ttests=TTests(where=(method="Pooled"));
run;
%ODSOn
/* Construct indicator var for obs that reject H0 at 0.05 significance */
data Results;
set TTests;
RejectH0 = (Probt <= 0.05);
run;
/* 3. Compute proportion: (# that reject H0)/NumSamples and CI */
%ODSoff
proc freq data=Results;
by Delta;
tables RejectH0 / nocum binomial(level='1');
output out=Est binomial;
run;
%ODSOn
/* for comparison, compute exact power */
proc power;
twosamplemeans power = . /* missing ==> "compute this" */
meandiff= 0 to 2 by 0.1 /* delta = 0, 0.1, ..., 2 */
stddev=1 /* N(delta, 1) */
ntotal=20; /* 20 obs in the two samples */
ods output Output=Power; /* output results to data set */
run;
data Combine;
set Est Power;
run;
/* overlay the simulated results with the cuve of exact power */
title "Power of the t Test";
title2 "Samples are N(0,1) and N(delta,1), n1=n2=10";
proc sgplot data=Combine noautolegend;
series x=MeanDiff y=Power;
scatter x=Delta y=_BIN_ / yerrorlower=L_Bin yerrorupper=U_Bin;
inset ("Number of Samples"="&NumSamples") / border;
yaxis min=0 max=1 label="Power (1 - P[Type II Error])" grid;
xaxis label="Difference in Population Means" grid;
run;
扫码加好友，拉您进群
全部版块
我的主页
