*--- 3.10節 ベイズ統計の基礎 ;
data MYDATA ;
input Y N ;
cards ;
2 5
;
run ;
ods graphics on ;
proc mcmc data=MYDATA seed=777 nmc=1000 ;
parm THETA 0.5 ;
prior THETA ~ beta(1, 1) ;
model Y ~ binomial(N, theta) ;
run ;
ods graphics off ;
data MYDATA ;
input Y @@ ;
cards ;
19 24 26 27 29
;
run ;
ods graphics on ;
proc mcmc data=MYDATA seed=777 nmc=1000
outpost=RESULT ;
parm mu 25 ;
prior mu ~ normal(mean=27, var= 9) ;
model Y ~ normal(mean=mu, var=16) ;
run ;
ods graphics off ;
data MYDATA ;
input n y x1 x2 ;
cards ;
5 2 1 0
15 10 1 1
15 2 0 0
5 3 0 1
;
run ;
ods graphics on ;
proc mcmc data=MYDATA seed=777 nmc=5000 nbi=5000 thin=5 ;
parms (beta0 beta1 beta2) 0 ;
prior beta0 beta1 beta2 ~ normal(0, var = 1000000) ;
p = logistic(beta0 + beta1*x1 + beta2*x2) ;
model y ~ binomial(n,p) ;
run ;
ods graphics off ;
data MYDATA ;
input X Y @@ ;
cards ;
1 1 2 2 3 3 4 4 5 5.1
;
run ;
ods graphics on;
proc mcmc data=MYDATA seed=777 nmc=2000
outpost=RESULT monitor=(_parms_ sigma) ;
array tau[2] ;
array beta[2] ;
array sigma[2] ;
parms tau: 1 ;
parms beta: 0 ;
beginprior ;
prior tau1 ~ gamma(0.001, iscale=0.001) ;
hyper tau2 ~ gamma(0.001, iscale=0.001) ;
prior beta:~ normal(mean=0, var=1/tau2) ;
sigma1 = sqrt(1/tau1) ;
sigma2 = sqrt(1/tau2) ;
endprior ;
mu = beta1 + beta2*X ;
model Y ~ normal(mu, var=1/tau1) ;
run;
ods graphics off ;
*************************************************;
*--- 3章 演習問題 ;
* (1) GroupごとにWeightの平均値に関するグラフを作成してください。;
data EXAMPLE ;
call streaminit(7777) ;
do GROUP=1 to 2 ;
do N=1 to 100 ;
WEIGHT = rand("NORMAL", 1-3*GROUP, 10) ;
ADR = 2-rand("BINOMIAL", 0.1*GROUP, 1) ;
TIME = ceil(100*rand("UNIFORM")/GROUP) ;
GENDER = rand("BINOMIAL", 0.4, 1) ;
output ;
end ;
end ;
run ;
* (2) GroupごとにWEIGHTのようやく統計量を算出した後、2標本t検定と2標本Wilcoxon検定を用いてGroup間のWeightの比較をしてください。;
;
proc sgplot data=EXAMPLE ;
hbar GROUP / response=WEIGHT stat=mean numstd=1
limitstat=stddev limits=upper ;
run ;
* (3) GroupごとにADRの頻度を集計した後、X2検定を用いてGroup間の比較をしてください。;
;
proc means data=EXAMPLE ;
class GROUP ;
var WEIGHT ;
run ;
proc ttest data=EXAMPLE ;
class GROUP ;
var WEIGHT ;
run ;
proc npar1way wilcoxon data=EXAMPLE correct=no ;
class GROUP ;
var WEIGHT ;
run ;
* (4)GroupごとにADRの頻度を集計した後、X2検定を用いてGroup間の比較をしてください。;
;
proc glm data=EXAMPLE ;
class GROUP GENDER ;
model WEIGHT = GROUP GENDER / solution ss3 ;
lsmeans GROUP ;
run ;
* (5)GroupごとにADRの頻度を集計した後、X2検定を用いてGroup間の比較をしてください。;
;
proc freq data=EXAMPLE ;
tables GROUP*ADR / riskdiff nocol nopercent chisq ;
exact or ;
run ;
* (6) GroupごとにGenderで調整した上でのADRに関する調整オッズ比を算出してください。;
proc logistic data=EXAMPLE order=internal ;
class GROUP GENDER / ref=last param=ref ;
model ADR = GROUP GENDER / expb ;
run ;
* (7)Groupごとに「副作用が発言するまでの日数」について、カプラン・マイヤー法を用いて累積発言割合の推定を計算した後、ログランク検定を用いてGroup間の比較をしてください。;
proc lifetest data=EXAMPLE plots=(s) ;
time TIME * ADR(2) ;
strata GROUP ;
run;
* (8)Groupごとに、Genderで調整した上でのADRに関する調整ハザード比を算出してください。 ;
proc phreg data=EXAMPLE ;
class GROUP GENDER ;
model TIME * ADR(2) = GROUP GENDER / ties=exact;
run ;
* (9) 検定を4回行いました。（T1:p=0.001, T2:p=0.002, T3:p=0.003, T4:p=0.004)、この時ボンフェローニの方法で調整p値を算出してください。
;
data PVALUE ;
input TEST $ RAW_P ;
cards ;
T1 0.001
T2 0.002
T3 0.030
T4 0.400
;
run ;
proc multtest pdata=PVALUE bon out=OUTDATA ;
run ;
* (10) 「平均値の差=4, 標準偏差=10, α=5%, 検出力90%」と設定したときに「帰無仮説H0:平均値の差が0である」という帰無仮説に関する2標本t検定を行ったときに必要となる例数を算出してください。
;
proc power ;
twosamplemeans
test = diff
meandiff = 4
stddev = 10
alpha = 0.05
power = 0.9
ntotal = . ;
run ;
* (11)以下の文章をＳＡＳプログラムのコードで記述してください。
サイコロを3個振った時の目の合計が10～13になる確率を、趣味レーションにて計算してください。;
data SIMULATION;*(keep=I SUCCESS) ;
call streaminit(777) ;
do I=1 to 100000 ;
D1 = rand("TABLE",1/6,1/6,1/6,1/6,1/6,1/6) ;
D2 = rand("TABLE",1/6,1/6,1/6,1/6,1/6,1/6) ;
D3 = rand("TABLE",1/6,1/6,1/6,1/6,1/6,1/6) ;
if (10 <= sum(D1,D2,D3) <=13) then SUCCESS=1; else SUCCESS=0 ;
output ;
end ;
run ;
proc means data=SIMULATION nonobs mean ;
var SUCCESS ;
run ;