*================================================;
*--- 3.8節 例数設計 ;
proc power;
onesamplemeans
mean = 3
stddev = 4
alpha = 0.05
power = 0.8
ntotal = . ;
run ;
proc power;
onesamplemeans
mean = 3
stddev = 4
alpha = 0.05
power = .
ntotal = 20 ;
run ;
proc power;
pairedmeans
test = diff
meandiff = 3
stddev = 4
alpha = 0.05
corr = 0.6
power = 0.8
npairs = . ;
run ;
proc power ;
twosamplemeans
test = diff
meandiff = 3
stddev = 4
alpha = 0.05
power = 0.8
ntotal = . ;
run ;
proc power ;
twosamplemeans
test = diff
meandiff = 3
stddev = 4
alpha = 0.05
power = 0.8
npergroup = . ;
run ;
proc power ;
twosamplemeans
test = diff
alpha = 0.025
sides = u
groupmeans = (2,3)
groupweights = (1,2)
nulldiff = -1
stddev = 4
power = 0.8
ntotal = . ;
run ;
proc power;
onesamplefreq
test = z
method = normal
alpha = 0.05
nullproportion = 0.1
proportion = 0.2
power = 0.8
ntotal = . ;
run ;
proc power;
twosamplefreq
test = pchi
alpha = 0.05
groupproportions = (0.1 0.2)
nullproportiondiff = 0.00
power = 0.8
npergroup = . ;
run ;
proc power ;
twosamplefreq
test = pchi
alpha = 0.025
sides = u
groupproportions = (0.1 0.2)
nullproportiondiff = -0.1
groupweights = (1,2)
power = 0.8
ntotal = . ;
run ;
proc power;
twosamplesurvival
test = logrank
alpha = 0.05
curve("薬剤A") = (2):(0.9)
curve("薬剤B") = (2):(0.8)
groupsurvival = "薬剤A" | "薬剤B"
accrualtime = 0.001
followuptime = 5
power = 0.8
ntotal = . ;
run ;
proc power;
twosamplesurvival
test = logrank
alpha = 0.05
groupsurvexphazards = (0.05268 0.11157)
accrualtime = 0.001
followuptime = 5
power = 0.8
ntotal = . ;
run ;
proc power;
onewayanova
test = contrast
alpha = 0.05
groupmeans = 0 | 10 | 20 | 30
stddev = 25
power = 0.8
contrast = (-3 -1 1 3)
ntotal = . ;
run ;
proc power;
onewayanova
test = contrast
alpha = 0.05
groupmeans = 0 | 10 | 20 | 30
groupweights = (1 2 2 2)
stddev = 25
power = 0.8
contrast = (-3 -1 1 3)
ntotal = . ;
run ;
ods output OUTPUT=MYSAMPLESIZE ;
proc power ;
twosamplemeans
test = diff
meandiff = 1 2 3
stddev = 4 5 6
alpha = 0.05
power = 0.8 0.9
npergroup = . ;
run ;
ods output close ;
*================================================;
*--- 3.9節 乱数とシミュレーション ;
data RDATA01(keep=X) ;
X = ranuni(777) ; *--- 一様乱数を 1 個生成する ;
run ;
data RDATA02(keep=X) ;
do I=1 to 5 ;
X = ranuni(777) ; *--- 一様乱数を 5 個生成する ;
output ;
end ;
run ;
data RDATA03(keep=X1) ;
do I=1 to 10;
X1 = ranuni(100) ;
output ;
end ;
run ;
data RDATA04(keep=X2 X3) ;
do I=1 to 5;
X2 = ranuni(100) ;
X3 = ranuni(100) ;
output ;
end ;
run ;
data RDATA05(keep=X2 X3) ;
do I=1 to 5;
X2 = ranuni(100) ;
X3 = ranuni(200) ;
output ;
end ;
run ;
data RDATA06(keep=X2 X3) ;
SEED2 = 100 ;
SEED3 = 200 ;
do I=1 to 5 ;
call ranuni(SEED2, X2) ;
call ranuni(SEED3, X3) ;
output ;
end ;
run ;
data X ;
do I=1 to 10000 ;
X = ranuni(100) ; output ; end ;
data Y ;
do I=1 to 10000 ;
Y = ranuni(200) ; output ; end ;
data RDATA07 ;
merge X Y ;
run ;
proc gplot data=RDATA07 ;
plot Y*X ;
run ;
data RDATA08(keep=X1) ;
call streaminit(100) ;
do I=1 to 10;
X1 = rand('UNIFORM') ;
output ;
end ;
run ;
data RDATA09(keep=X2 X3) ;
call streaminit(100) ;
do I=1 to 5;
X2 = rand('UNIFORM') ;
X3 = rand('UNIFORM') ;
output ;
end ;
run ;
data RDATA10(keep=X2 X3) ;
do I=1 to 5;
call streaminit(100) ;
X2 = rand('UNIFORM') ;
call streaminit(200) ;
X3 = rand('UNIFORM') ;
output ;
end ;
run ;
data X ;
call streaminit(100) ;
do I=1 to 10000 ;
X = rand('UNIFORM') ; output ; end ;
data Y ;
call streaminit(200) ;
do I=1 to 10000 ;
Y = rand('UNIFORM') ; output ; end ;
data RDATA11 ;
merge X Y ;
run ;
proc gplot data=RDATA11 ;
plot Y*X ;
run ;
data RDATA12 ;
call streaminit(777) ;
do I=1 to 5;
X = rand("HYPER", 40, 3, 20) ;
Y = rand("TABLE",0.1,0.2,0.3,0.4) ;
output ;
end ;
run ;
/*
data RDATA12 ;
call streaminit(100) ;
do I=1 to 5;
A = rand("BERNOULLI", 0.1) ;
B = rand("BETA", 2, 3) ;
C = rand("BINOMIAL", 0.4, 50) ;
D = rand("CAUCHY") ;
E = rand("CHISQUARE", 6) ;
F = rand("ERLANG", 7) ;
G = rand("EXPONENTIAL") ;
H = rand("F", 8, 9) ;
I = rand("GAMMA", 10) ;
J = rand("GEOMETRIC", 0.1) ;
K = rand("HYPER", 40, 3, 20) ;
L = rand("LOGNORMAL") ;
M = rand("NEGBINOMIAL", 0.5, 6) ;
N = rand("NORMAL", 7, 8) ;
O = rand("POISSON", 9) ;
P = rand("T", 10) ;
Q = rand("TABLE",0.1,0.2,0.3,0.4) ;
R = rand("TRIANGLE", 5) ;
S = rand("UNIFORM") ;
T = rand("WEIBULL", 6, 7) ;
output ;
end ;
run ;
*/
data MV ;
input _TYPE_ $ _NAME_ $ X Y Z ;
cards ;
MEAN . 1.0 2.0 3.0
COV X 4.0 0.0 8.0
COV Y 0.0 9.0 0.0
COV Z 8.0 0.0 16.0
;
run ;
proc simnormal
data = MV(type=cov)
out = RDATA13
numreal = 1000
seed = 777 ;
var X Y Z ;
run ;
proc corr data=RDATA13 ;
var X Y Z ;
run ;
data RDATA14 ;
N = 10000 ;
call streaminit(777) ;
do I=1 to N ;
X = -log(rand('UNIFORM')) ;
output ;
end ;
run ;
proc univariate data=RDATA14 ;
histogram X / cfill=blue
midpoints=0 to 12 by 0.1 ;
run ;
data RDATA15 ;
N = 10000 ;
call streaminit(777) ;
do I=1 to N ;
U = rand('UNIFORM') ;
if (U < 1/2) then X = sqrt(2*U) ;
else X = 2-sqrt(2*(1-U)) ;
output ;
end ;
run ;
proc univariate data=RDATA15 ;
histogram X / cfill=blue
midpoints=-1 to 3 by 0.1 ;
run ;
data RDATA16 ;
I = 1 ;
N = 10000 ;
call streaminit(777) ;
do while (I <= N) ;
V = 2*rand('UNIFORM') ;
U = rand('UNIFORM') ;
if (0 <= V <= 1) then H = V ;
else if (1 < V <= 2) then H = 2-V ;
else H = 0 ;
if (U <= H) then do ;
X = V ;
I = I+1 ;
output ;
end ;
end ;
run ;
proc univariate data=RDATA16 ;
histogram X / cfill=blue
midpoints=-1 to 3 by 0.1 ;
run ;
data RDATA17 ;
N = 10000 ;
call streaminit(777) ;
do I=1 to N ;
X = rand('UNIFORM')+rand('UNIFORM')+rand('UNIFORM')+
rand('UNIFORM')+rand('UNIFORM')+rand('UNIFORM')+
rand('UNIFORM')+rand('UNIFORM')+rand('UNIFORM')+
rand('UNIFORM')+rand('UNIFORM')+rand('UNIFORM')-6 ;
output ;
end ;
run ;
proc univariate data=RDATA17 ;
histogram X / cfill=blue
midpoints=-4 to 4 by 0.5 ;
run ;
data SDATA01(keep=COIN) ;
call streaminit(777) ;
X = rand('UNIFORM') ;
if (X > 1/2) then COIN = 1 ; *--- 表 ;
else COIN = 0 ; *--- 裏 ;
run ;
data SDATA02(keep=DICE) ;
call streaminit(777) ;
X = rand('UNIFORM') ;
if (X < 1/6) then DICE = 1 ;
else if (X < 2/6) then DICE = 2 ;
else if (X < 3/6) then DICE = 3 ;
else if (X < 4/6) then DICE = 4 ;
else if (X < 5/6) then DICE = 5 ;
else DICE = 6 ;
run ;
data SDATA03(keep=COIN) ;
call streaminit(777) ;
COIN = rand("BERNOULLI", 0.5) ;
run ;
data SDATA04(keep=DICE) ;
call streaminit(777) ;
DICE = rand("TABLE",1/6,1/6,1/6,1/6,1/6,1/6) ;
run ;
data SDATA05(keep=I COIN) ;
do I=1 to 4 ;
COIN = rand("BERNOULLI", 0.5) ;
output ;
end ;
run ;
data SDATA06(keep=I DICE) ;
do I=1 to 4 ;
DICE = rand("TABLE",1/6,1/6,1/6,1/6,1/6,1/6) ;
output ;
end ;
run ;
ods output SUMMARY=SDATA07 ;
proc means data=SDATA06 nonobs sum ;
var DICE ;
run ;
ods output close ;
data SDATA08(keep=I DICE) ;
call streaminit(777) ;
do I=1 to 1000 ;
do j=1 to 4 ;
DICE = rand("TABLE",1/6,1/6,1/6,1/6,1/6,1/6) ;
output ;
end ;
end ;
run ;
ods listing close ;
ods output SUMMARY=SDATA09 ;
proc means data=SDATA08 nonobs sum ;
by I ;
var DICE ;
run ;
ods output close ;
ods listing ;
proc means data=SDATA09 nonobs n mean ;
var DICE_SUM ;
run ;
data NDATA01(keep=I N GROUP QOL);
call streaminit(777) ;
array MU(2) M1-M2 ;
M1=6.5; M2=4.0 ; SIGMA=3.0 ; N=20 ;
do I=1 to 1000 ;
do GROUP=1 to 2;
do J=1 to N ;
QOL = rand("NORMAL", MU(GROUP), SIGMA) ;
output ;
end ;
end ;
end ;
run ;
ods listing close ;
ods output ttests=RESULT01(where=(method="Pooled")) ;
proc ttest data=NDATA01 ;
by I ;
class GROUP ;
var QOL ;
run ;
ods output close ;
ods listing ;
data POWER01 ;
set RESULT01 ;
if (tValue > 0) and (Probt < 0.05)
then FLAG = 1 ; *--- 有意差あり ;
else FLAG = 0 ; *--- 有意差なし ;
run ;
proc freq data=POWER01 ;
table FLAG ; *--- FLAG=1 の割合を求める ;
run ;
proc power ;
twosamplemeans
test = diff
meandiff = 2.5
stddev = 3.0
alpha = 0.05
power = .
npergroup = 20 ;
run ;
data NDATA02(keep=I N GROUP EVENT K);
call streaminit(777) ;
array PROP(2) P1-P2 ;
P1=0.4; P2=0.2 ; N=100 ;
do I=1 to 1000 ;
do GROUP=1 to 2;
X = rand("BINOMIAL", PROP(GROUP), N) ;
EVENT = "Y"; K = X ; output ;
EVENT = "N"; K = N-X ; output ;
end ;
end ;
run ;
ods listing close ;
ods output ChiSq=ChiSq(where=(Statistic="カイ 2 乗値")) ;
proc freq data=NDATA02 ;
by I ;
weight K ;
tables GROUP*EVENT / chisq ;
run ;
ods output close ;
ods listing ;
ods listing close ;
ods output RiskDiffCol2=RiskDiff(where=(Row="差")) ; ;
proc freq data=NDATA02 ;
by I ;
weight K ;
tables GROUP*EVENT / riskdiff ;
run ;
ods output close ;
ods listing ;
data RESULT02 ;
merge ChiSq RiskDiff ;
by I ;
keep I Prob Risk ;
run ;
data POWER02 ;
set RESULT02 ;
if (Risk > 0) and (Prob < 0.05)
then FLAG = 1 ; *--- 有意差あり ;
else FLAG = 0 ; *--- 有意差なし ;
run ;
proc freq data=POWER02 ;
table FLAG ; *--- FLAG=1 の割合を求める ;
run ;
proc power;
twosamplefreq
test = pchi
alpha = 0.05
groupproportions = (0.2 0.4)
nullproportiondiff = 0.00
power = .
npergroup = 100 ;
run ;