*--- 3.1節使用するデータセット ;

*--- 3.1節使用するデータセット ;
data QOL ;
input GROUP $ QOL EVENT DAY PREDRUG $ DURATION ;
cards ;
A 15 1 50 NO 1
A 13 1 200 NO 3
A 11 1 250 NO 2
A 11 1 300 NO 4
A 10 1 350 NO 2
A 9 1 400 NO 2
A 8 1 450 NO 4
A 8 1 550 NO 2
A 6 1 600 NO 5
A 6 1 100 NO 7
A 4 2 250 NO 4
A 3 2 500 NO 6
A 3 2 750 NO 3
A 3 2 650 NO 7
A 1 2 1000 NO 8
A 6 1 150 YES 6
A 5 1 700 YES 5
A 4 2 800 YES 7
A 2 2 900 YES 12
A 2 2 950 YES 10
B 13 1 380 NO 9
B 12 1 880 NO 5
B 11 1 940 NO 2
B 4 2 20 NO 7
B 4 2 560 NO 2
B 5 1 320 YES 11
B 5 1 940 YES 3
B 4 2 80 YES 6
B 3 2 140 YES 7
B 3 2 160 YES 13
B 3 2 240 YES 15
B 2 2 280 YES 9
B 2 2 440 YES 8
B 2 2 520 YES 7
B 2 2 620 YES 9
B 2 2 740 YES 8
B 2 2 860 YES 2
B 1 2 880 YES 10
B 0 2 920 YES 8
B 0 2 960 YES 4
C 9 1 170 NO 1
C 7 1 290 NO 4
C 5 1 430 NO 2
C 3 2 610 NO 4
C 2 2 110 NO 5
C 2 2 410 NO 2
C 1 2 530 NO 7
C 1 2 580 NO 2
C 0 2 810 NO 3
C 0 2 990 NO 10
C 6 1 30 YES 1
C 5 1 830 YES 6
C 3 2 70 YES 16
C 2 2 310 YES 9
C 2 2 370 YES 18
C 1 2 490 YES 7
C 1 2 690 YES 10
C 0 2 730 YES 3
C 0 2 770 YES 12
C 0 2 910 YES 8
;
run ;

*************************************************;

*--- 3.2節 SASによる統計解析の流れ ;

data MYDATA ;
set QOL ;
where GROUP="A" ;
run ;

proc print data=MYDATA ;
var QOL ; *--- 出力する変数名 ;
run ;

title "ヒストグラムと密度推定" ;
proc sgplot data=MYDATA ;
histogram QOL ; *--- ヒストグラム ;
density QOL / type=kernel ; *--- 密度推定曲線 ;
run ;

title "棒グラフ";
proc sgplot data=MYDATA ;
vbar PREDRUG ;
run ;
title "" ; *--- タイトルを初期化 ;

proc univariate data=MYDATA ;
var QOL ; *--- 要約統計量を求める変数名 ;
run ;

proc univariate data=MYDATA ;
class PREDRUG ; *--- カテゴリ変数名 ;
var QOL ; *--- 要約統計量を求める変数名 ;
run ;

title "箱ひげ図" ;
proc sgplot data=MYDATA ;
vbox QOL ; *--- hbox にすると横向きになる ;
run ;
title "" ; *--- タイトルを初期化 ;

proc univariate data=MYDATA mu0=3 ; *--- 平均が 3 かどうか ;
var QOL ; *--- 検定の対象とする変数 ;
run ;

ods html file="C:\temp\ファイル名.html" ; *--- html ;
title "要約統計量" ;
proc univariate data=MYDATA ;
var QOL ;
run;
ods html close;

ods rtf file="C:\temp\ファイル名.rtf" ; *--- rtf(WORD) ;
title "要約統計量" ;
proc univariate data=MYDATA ;
var QOL ;
run;
ods rtf close ;

ods pdf file="C:\temp\ファイル名.pdf" ; *--- pdf ;
title "要約統計量" ;
proc univariate data=MYDATA ;
var QOL ;
run;
ods pdf close ;

ods rtf file="C:\temp\ファイル名.rtf" ; *--- pdfも同様に出力可 ;
title "密度推定" ;
proc sgplot data=MYDATA ;
density QOL / type=kernel ;
run ;
ods rtf close;

filename mygraph "c:\temp\ファイル名.emf" ; *--- emf ファイル;
goptions reset=all gsfname=mygraph device=emf ;
title "密度推定" ;
proc sgplot data=MYDATA ;
density QOL / type=kernel ;
run ;
filename mygraph clear ;

ods graphics on / imagefmt=emf ;

title "" ; *--- タイトルを初期化 ;

ods output Moments=MYOUTPUT ;
proc univariate data=MYDATA ;
var QOL ;
run;
ods output close ;

ods trace on /listing; *--- 統計量のラベルを出力する命令 ;
proc univariate data=MYDATA ;
var QOL ;
run ;

ods trace off ;

ods output TestsForLocation=TTEST(where=(Test="Student の t 検定")) ;
proc univariate data=MYDATA ;
var QOL ;
run;
ods output close ;

proc ttest data=MYDATA h0=3 sides=2 ; *--- sides=2（両側） ;
var QOL ;
run ;

proc ttest data=MYDATA h0=3 sides=l ; *--- sides=l（下側） ;
var QOL ;
run ;

proc ttest data=MYDATA h0=3 sides=u ; *--- sides=u（上側） ;
var QOL ;
run ;

quit ;

*************************************************;

*--- 3.3節散布図，回帰直線，相関係数 ;

data MYDATA ;
set QOL ;
where GROUP="A" ;
run ;

proc means data=MYDATA nonobs ;
class GROUP ;
var DURATION QOL ;
run ;

title "散布図と相関係数" ;
proc sgplot data=MYDATA ;
scatter x=DURATION y=QOL ; *--- 散布図 ;
reg x=DURATION y=QOL ; *--- 回帰直線 ;
run ;

proc reg data=MYDATA ;
model QOL = DURATION ;
run ;

proc corr data=MYDATA pearson spearman ;
var DURATION QOL ;
run ;

quit ;

*************************************************;

*--- 3.4節平均値の比較と回帰分析 ;

data MYDATA ;
set QOL ;
where GROUP ne "C" ;
run ;

title "棒グラフ（平均値）" ;
proc sgplot data=MYDATA ;
vbar GROUP / response=QOL stat=mean ;
run ;

title "棒グラフ（平均値±標準偏差）" ;
proc sgplot data=MYDATA ;
hbar GROUP / response=QOL stat=mean numstd=1
limitstat=stddev limits=upper ;
run ;

title "棒グラフ" ;
proc sgpanel data=MYDATA ;
panelby PREDRUG / rows=1 columns=2 ;
vbar GROUP / response=QOL stat=mean numstd=1
limitstat=stddev limits=upper ;
run ;

title "" ;

proc summary data=MYDATA print ;
class GROUP ;
var QOL ;
run ;

/*
proc means data=MYDATA NONOBS N MEAN STDDEV ;
class GROUP ;
var QOL ;
run ;
*/

ods graphics on / imagefmt=emf ;
proc ttest data=MYDATA ;
class GROUP ;
var QOL ;
run ;

proc glm data=MYDATA ;
class GROUP ;
model QOL = GROUP / solution ss3 ;
run ;

proc glm data=MYDATA ;
class GROUP ;
model QOL = GROUP / noint solution ss3 ;
run ;

ods graphics on / imagefmt=emf ;
proc freq data=MYDATA ;
table GROUP*PREDRUG / nocol nopercent ;
run ;

proc glm data=MYDATA ;
class GROUP PREDRUG ;
model QOL = GROUP PREDRUG GROUP*PREDRUG / solution ss3 ;
lsmeans GROUP ;
run ;

proc glm data=MYDATA ;
class GROUP ;
model QOL = GROUP DURATION GROUP*DURATION / solution ss3 ;
lsmeans GROUP ;
run ;

proc means data=MYDATA NONOBS N MEAN STDDEV ;
class GROUP ;
var DURATION ;
run ;

proc means data=MYDATA NONOBS N MEAN STDDEV ;
var DURATION;
run ;

proc mixed data=MYDATA ;
class GROUP ;
model QOL = GROUP DURATION GROUP*DURATION / solution ;
lsmeans GROUP / at DURATION=10 ; *--- 10 を代入 ;
lsmeans GROUP / at means ; *--- 全体の平均を代入 ;
run ;

ods graphics on / imagefmt=emf ;
proc glm data=MYDATA ;
class PREDRUG GROUP ;
model QOL = GROUP PREDRUG GROUP*PREDRUG ;
run ;

ods graphics on / imagefmt=emf ;
proc glm data=MYDATA ;
class GROUP ;
model QOL = GROUP DURATION GROUP*DURATION ;
run ;

proc glm data=MYDATA ; *--- モデル1 ;
class GROUP ;
model QOL = GROUP DURATION / solution ss1 ss2 ss3 ;
run ;

proc glm data=MYDATA ; *--- モデル2 ;
class GROUP ;
model QOL = DURATION GROUP / solution ss1 ss2 ss3 ;
run ;

proc glm data=MYDATA ; *--- モデル3 ;
class GROUP ;
model QOL = GROUP DURATION GROUP*DURATION / solution ss1 ss2 ss3 ;
run ;

ods graphics on / imagefmt=emf ;
proc glm data=MYDATA ;
class PREDRUG GROUP ;
model QOL = GROUP PREDRUG / solution ss3 ;
lsmeans GROUP ;
run ;

ods graphics on / imagefmt=emf ;
proc glm data=MYDATA ;
class GROUP ;
model QOL = GROUP DURATION / solution ss3 ;
lsmeans GROUP ;
run ;

proc npar1way hl wilcoxon data=MYDATA correct=no ;
class GROUP ;
var QOL ;
exact wilcoxon ;
run ;

data QOL1 ;
input GROUP $ X R ;
cards ;
A 3 1
B 5 2
A 6 3
A 8 4
B 11 5
;
run ;

proc npar1way data=QOL1 correct=no ;
class GROUP ;
var X ;
exact wilcoxon ;
run ;

proc npar1way data=QOL1 correct=yes ;
class GROUP ;
var X ;
exact wilcoxon ;
run ;

proc ttest data=MYDATA h0=-1 side=u ;
class GROUP ;
var QOL ;
run ;

proc mixed data=MYDATA ;
class GROUP PREDRUG ;
model QOL = GROUP / solution ;
random PREDRUG ;
run ;

data QOL0 ;
input PATIENT GROUP TREAT $ PERIOD QOL ;
cards ;
1 1 A 1 1
1 1 B 2 3
2 1 A 1 2
2 1 B 2 4
3 2 B 1 5
3 2 A 2 4
4 2 B 1 7
4 2 A 2 4
;
run ;

proc mixed data=QOL0 ;
class PATIENT GROUP TREAT PERIOD ;
model QOL = GROUP TREAT PERIOD ;
random PATIENT(GROUP) ;
lsmeans TREAT / pdiff;
run ;

proc glm data=QOL0 ;
class PATIENT GROUP TREAT PERIOD ;
model QOL = GROUP PATIENT(GROUP) TREAT PERIOD / SS3 ;
test H=GROUP E=PATIENT(GROUP) ;
lsmeans TREAT / pdiff ;
run ;

ods graphics on / imagefmt=emf ;
proc genmod data=MYDATA ;
class GROUP ;
model QOL = GROUP / dist=normal ;
bayes seed=777 ;
run;
ods graphics off ;

quit ;

*************************************************;

*--- 3.5節 2値データの比較とロジスティック回帰 ;

data MYDATA ;
set QOL ;
where GROUP ne "C" ;
run ;

proc format ;
value EVENTF 1="Yes" 2="No" ;
run ;

title "棒グラフ" ;
proc sgplot data=MYDATA ;
hbar GROUP / GROUP=EVENT ;
format EVENT EVENTF. ;
run ;

title "棒グラフ" ;
proc sgpanel data=MYDATA ;
panelby PREDRUG / rows=2 columns=1 ;
hbar GROUP / GROUP=EVENT ;
format EVENT EVENTF. ;
run ;

ods output CrossTabFreqs=OUT1(where=(GROUP ne "" and EVENT ne .)) ;
proc freq data=MYDATA ;
table GROUP*EVENT / nopercent nocol ;
run ;
ods output close ;

title "棒グラフ（改善の有無）" ;
proc sgplot data=OUT1 ;
hbar GROUP / GROUP=EVENT freq=ROWPERCENT ;
format EVENT EVENTF. ;
run ;

title "" ;

proc freq data=MYDATA ;
tables GROUP*EVENT / riskdiff relrisk nocol nopercent ;
exact or ;
format EVENT EVENTF. ;
run ;

data MYFREQ ;
input GROUP $ YN ;
cards ;
A 2
A 1
A 1
B 1
B 2
;
run ;

proc format ;
value YNF 1="Male"
2="Female" ;
run ;

proc freq data=MYFREQ order=data ;
tables GROUP*YN;
format YN YNF. ;
run ;

proc freq data=MYDATA order=formatted ;
tables GROUP*EVENT / riskdiff relrisk nocol nopercent ;
exact or ;
format EVENT EVENTF. ;
run ;

proc freq data=MYDATA ;
tables GROUP*EVENT / chisq ;
format EVENT EVENTF. ;
run ;

proc logistic data=MYDATA order=internal ;
class GROUP / ref=last param=ref ;
model EVENT = GROUP / expb ;
format EVENT EVENTF. ;
run ;

proc logistic data=MYDATA order=internal ;
class GROUP / param=glm ;
model EVENT = GROUP / expb noint ;
format EVENT EVENTF. ;
run ;

proc logistic data=MYDATA order=internal ;
model EVENT = DURATION / expb ;
units DURATION=5 ;
format EVENT EVENTF. ;
run ;

proc logistic data=MYDATA order=internal ;
class GROUP PREDRUG / ref=last param=ref ;
model EVENT = GROUP PREDRUG / expb ;
format EVENT EVENTF. ;
run ;

proc logistic data=MYDATA order=internal ;
class GROUP / ref=last param=ref ;
model EVENT = GROUP DURATION / expb ;
format EVENT EVENTF. ;
run ;

proc sort data=MYDATA ;
by PREDRUG ;
run ;

proc freq data=MYDATA order=internal ;
by PREDRUG ;
tables GROUP*EVENT / riskdiff relrisk nocol nopercent ;
exact or ;
format EVENT EVENTF. ;
run ;

proc logistic data=MYDATA order=internal ;
class GROUP PREDRUG / ref=last param=ref ;
model EVENT = GROUP PREDRUG GROUP*PREDRUG / expb ;
format EVENT EVENTF. ;
run ;

proc freq data=MYDATA order=data ;
tables PREDRUG*GROUP*EVENT / nocol nopercent chisq cmh ;
format EVENT EVENTF. ;
run ;

proc logistic data=MYDATA order=internal ;
class GROUP / ref=last param=ref ;
model EVENT = GROUP DURATION GROUP*DURATION / expb ;
format EVENT EVENTF. ;
run ;

data MYFREQ2 ;
input GROUP $ EVENT $ ;
cards ;
A Y
A N
A N
A N
A N
B Y
B Y
B N
B N
B N
;
run ;

data MYFREQ3 ;
input GROUP $ EVENT $ N ;
cards ;
A Y 1
A N 4
B Y 2
B N 3
;
run ;

proc freq data=MYFREQ3 ;
weight N ;
tables GROUP*EVENT ;
run ;

data MYFREQ4 ;
input PREDRUG $ GROUP $
EVENT NTOTAL ;
cards ;
Y A 2 5
Y B 2 13
N A 10 15
N B 3 5
;
run ;

proc logistic data=MYFREQ4 order=data ;
class GROUP / ref=last param=ref ;
model EVENT/NTOTAL = GROUP / expb ;
run ;

data MYFREQ5 ;
input PREDRUG $ GROUP $
EVENTYN $ N ;
cards ;
Y A Y 2
Y A N 3
Y B Y 2
Y B N 13
N A Y 10
N A N 5
N B Y 3
N B N 2
;
run ;

proc logistic data=MYFREQ5 order=data ;
freq N ;
class GROUP / ref=last param=ref ;
model EVENTYN = GROUP / expb ;
run ;

proc freq data=MYDATA ;
tables GROUP*EVENT / riskdiff(margin=0.1 method=FM noninf column=1)
alpha=0.025 ;
format EVENT EVENTF. ;
run ;

proc logistic data=MYDATA descending ;
model GROUP = QOL / outroc=ROC ;
run ;

data ROC (keep=X Y) ;
set ROC ;
X = 100*_1MSPEC_ ; *--- 1-特異度 ;
Y = 100*_SENSIT_ ; *--- 感度 ;
run ;

proc sgplot data=ROC ;
series x=X y=Y ;
run ;

ods graphics on ;
proc genmod data=MYDATA order=internal ;
class GROUP ;
model EVENT = GROUP / dist=bin link=logit ;
bayes seed=777 ;
run;
ods graphics off ;

quit ;

*************************************************;

*--- 3.1節 使用するデータセット ;

*--- 3.1節使用するデータセット ;