*------------------------------------------------------------------;
* 昭和大学　第９回実践臨床統計学セミナー・ハンズオン2 			                   ;
* 診断法・予測モデルの性能評価		 	                                   ;
* 2020年9月7日 　　                                                    ;
*------------------------------------------------------------------;


*------------------------------------------------------------------;
* 1. FREQプロシジャ： 感度・特異度・陽性的中率・陰性的中率・診断精度の計算			　  ;
*------------------------------------------------------------------;

data asthma;
	input Test Disease n;
	datalines;
	0 0  198
	0 1  52
	1 0 85
	1 1 173
	;

proc freq data=asthma;
    where Disease=1;
    weight n;
    tables Test / binomial(level="1");
    exact binomial;
    run;

 *-- 感度の計算。;
 *-- proc freqは、1列のデータを対象とした場合には、2項確率の信頼区間を計算します。;


proc freq data=asthma;
    where Disease=0;
    weight n;
    tables Test / binomial(level="0");
    exact binomial;
    run;

 *-- 特異度の計算。;
 *-- proc freqは、1列のデータを対象とした場合には、2項確率の信頼区間を計算します。;


proc freq data=asthma;
    where Test=1;
    weight n;
    tables Disease / binomial(level="1");
    exact binomial;
    run;

 *-- 陽性的中率の計算。;
 *-- proc freqは、1列のデータを対象とした場合には、2項確率の信頼区間を計算します。;


proc freq data=asthma;
    where Test=0;
    weight n;
    tables Disease / binomial(level="0");
    exact binomial;
    run;

 *-- 陰性的中率の計算。;
 *-- proc freqは、1列のデータを対象とした場合には、2項確率の信頼区間を計算します。;


data acc; 
	set asthma;
	ACC = (Test = Disease); 
	run;

proc freq;
	weight n;
	tables ACC / binomial(level="1");
	exact binomial;
	run;

 *-- 診断精度の計算。;
 *-- proc freqは、1列のデータを対象とした場合には、2項確率の信頼区間を計算します。;

 
 


*------------------------------------------------------------------;
* 2. CSVデータファイルの読み込み　　　　 　　　　　　　　　                       			;
*------------------------------------------------------------------;

proc import out=work.aSAH
datafile = "C:\SAS2\aSAH.csv"
dbms = csv replace;
getnames = yes;
datarow = 2;
run;


 *-- datafile=... のフォルダ・ファイルのパスは適当に変更してください。;
 *-- Cドライブ直下に "SAS2" というフォルダを作っていただくと、以後の実習でも、プログラムの改変なしに、配布したコードでデモを実行することができます。;





*------------------------------------------------------------------;
* 3. ROC曲線の描画とAUCの計算　　　　　　　　　 　　　　　　　　　                        ;
*------------------------------------------------------------------;

proc logistic data=aSAH;
	class outcome / ref=first param=glm;
	model outcome = s100b / nofit;
	roc 'S100b' s100b;
run;

 *-- S100b についてのROC曲線を作成します。;



proc logistic data=aSAH;
	class outcome / ref=first param=glm;
	model outcome = s100b ndka / nofit;
	roc 'S100b' s100b;
	roc 'NDKA' ndka;
	roccontrast reference("s100b") / estimate e;
run;

 *-- S100b と NDKA についてのROC曲線を作成します。;
 *-- roccontrast ステートメントをつけることによって、2つのROC曲線のAUCを比較した検定を行うことができます。;



proc logistic data=aSAH;
	class outcome gender / ref=first param=glm;
	model outcome = age gender ndka s100b wfns / nofit;
	roc 'S100b' s100b;
	roc 'Multivariate model' age gender ndka s100b wfns;
	roccontrast reference("s100b") / estimate e;
run;

 *-- 多変量ロジスティック回帰モデルによって、複数の予測因子を利用した多変量モデルに基づくリスクスコアを作成することができます。;
 *-- このリスクスコアによるROC曲線を作成することができます。;





*------------------------------------------------------------------;
* 4. Hosmer-Lemeshow検定　　　　　　　　　 　　　　　　　　　                        ;
*------------------------------------------------------------------;

proc logistic data=aSAH;
	class outcome gender / ref=first param=glm;
	model outcome = age gender ndka s100b wfns / lackfit;
run;






*------------------------------------------------------------------;
* 5. Calibrationプロットの作成　　　　　　　　　 　　　　　　　　　                        ;
*------------------------------------------------------------------;

proc logistic data=aSAH noprint;
	class outcome gender / ref=first param=glm;
	model outcome = age gender ndka s100b wfns;
	output out=LogiOut predicted=PredProb;
run;

/* To construct the decile calibration plot, identify deciles of the predicted prob. */
proc rank data=LogiOut out=LogiDecile groups=10;
   var PredProb;
   ranks Decile;
run;
 
/* Then compute the mean predicted prob and the empirical proportions (and CI) for each decile */
proc means data=LogiDecile noprint;
   class Decile;
   types Decile;
   var y PredProb;
   output out=LogiDecileOut mean=yMean PredProbMean
          lclm=yLower uclm=yUpper;
run;
 
title "Calibration Plot";
proc sgplot data=LogiDecileOut noautolegend aspect=1;
   lineparm x=0 y=0 slope=1 / lineattrs=(color=grey pattern=dash);
   *loess x=PredProbMean y=yMean;  /* if you want a smoother based on deciles */
   series x=PredProbMean y=yMean;  /* if you to connect the deciles */
   scatter x=PredProbMean y=yMean / yerrorlower=yLower yerrorupper=yUpper;
   yaxis label="Observed Probability of Outcome";
   xaxis label="Predicted Probability of Outcome";
   xaxis min=0;
   xaxis max=1;
run;