*------------------------------------------------------------------;
* 昭和大学　第15回実践臨床統計学セミナー・ハンズオン2 			           ;
* 傾向スコア重み付き解析 	 	                           ;
* 2021年3月8日 　　                                                  ;
*------------------------------------------------------------------;


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

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

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





*------------------------------------------------------------------;
* 2. FREQ プロシジャ： 分割表の解析（未調整） 　　　　　　　　　    			   ;
*------------------------------------------------------------------;

proc freq data=cattaneo2;
table mismoke*Y / nopercent nocol riskdiff relrisk chisq;
exact or;
run;

 *-- 単純な分割表の解析によって、交絡調整をしていないリスク差・リスク比を計算します。;
 
 
 
 
 
*------------------------------------------------------------------;
* 3. CAUSALTRT プロシジャ： 傾向スコア重み付け解析　　　　　　　    			   ;
*------------------------------------------------------------------;

proc causaltrt data=cattaneo2 desc covdiffps poutcomemod;
class Y mismoke;
psmodel mismoke(ref='0')=mage medu mmarried mrace mhisp foreign alcohol fbaby fage fedu / plots = all;
model Y(ref='0') / dist=bin;
run;

 *-- CAUSALTRTプロシジャを使うことによって、傾向スコア重み付け解析を行うことができます。;
 *-- ただし、アウトカム指標について、リスク差しか出力することができません。;
 *-- 相対リスクの指標として、リスク比の傾向スコア重み付け推定量を計算したいのであれば、3番のマクロをご使用ください。;






*------------------------------------------------------------------;
* 4. PSMATCH プロシジャ： 重み付け前後の共変量の分布の比較　    			   ;
*------------------------------------------------------------------;

proc psmatch data=cattaneo2;
class mismoke mmarried fbaby;
psmodel mismoke(treated='1') = mage medu mmarried mrace mhisp foreign alcohol fbaby fage fedu;
assess var = (mage medu mmarried mrace mhisp foreign alcohol fbaby fage fedu)
  / weight=atewgt varinfo;
run;

 *-- PSMATCHプロシジャを使うことによって、重み付けの前後での要約統計量を計算することができます。;







*------------------------------------------------------------------;
* 5. マクロ： 傾向スコア重み付けリスク比の計算			　    			   ;
*------------------------------------------------------------------;


*-- https://github.com/t-yui/causaltrtRiskRatioCI/blob/master/causaltrtRiskRatioCI.sas ;
*-- に、リスク比の推定のためのマクロがある。;

%macro causaltrtRiskRatioCI(outdata_causaltrt, column_no, outdata_result, alpha=0.05);
	/*
	parameters
	----------
	outdata_causaltrt : dataset
		Specify CausalEffects dataset from CAUSALTRT procedure.
	column_no : integer
		Specify column number of standard error to be used.
	outdata_result : dataset
		Specify name of dataset to be output.
	alpha : float
		Specify significance level.
	*/
	proc iml;
		/* get estimates and btstd */
		use &outdata_causaltrt.;
		read all into data;
		risk_pom1 = data[1,1];
		risk_pom0 = data[2,1];
		std_pom1 = data[1,&column_no.];
		std_pom0 = data[2,&column_no.];
		std_diff = data[3,&column_no.];
		
		/* calculate var and covar */
		var_pom1 = std_pom1 ** 2;
		var_pom0 = std_pom0 ** 2;
		var_diff = std_diff ** 2;
		covar = (var_pom1 + var_pom0 - var_diff) / 2;
		
		/* calculate var of risk ratio using delta method */
		var_ratio = 0;
		var_ratio = var_ratio + var_pom1 / (risk_pom1 ** 2);
		var_ratio = var_ratio + var_pom0 / (risk_pom0 ** 2);
		var_ratio = var_ratio - 2 * covar / (risk_pom1 * risk_pom0);
		var_ratio = var_ratio * ((risk_pom1 / risk_pom0) ** 2);
		
		/* calculate 95%CI */
		risk_ratio = risk_pom1 / risk_pom0;
		alpha = &alpha.;
		percent = 1 - alpha / 2;
		percentile = quantile("normal", percent);
		low_ci = risk_ratio - percentile * sqrt(var_ratio);
		up_ci = risk_ratio + percentile * sqrt(var_ratio);
		
		/* output */
		print
			risk_pom1
			risk_pom0
			risk_ratio
			var_pom1
			var_pom0
			var_ratio
			low_ci
			up_ci
		;
		
		/* output as dataset*/
		create &outdata_result. 
			var {
				risk_pom1
				risk_pom0
				risk_ratio
				var_pom1
				var_pom0
				var_ratio
				low_ci
				up_ci
			};
		append;
		close &outdata_result.;
	quit;
%mend;



ODS Output CausalEffects=out_ce1;
proc causaltrt data=cattaneo2 desc covdiffps poutcomemod;
class Y mismoke;
psmodel mismoke(ref='0') = mage medu mmarried mrace mhisp foreign alcohol fbaby fage fedu / plots = all;
model Y(ref='0') / dist=bin;
run;

%causaltrtRiskRatioCI(out_ce1, 4, out_rr1)