卡方檢驗的多組比較

        類似於對連續型變量做方差分析(ANOVA)，當對分類變量的頻數做RxC列聯表(R>2)的卡方檢驗時，如果p<0.05則拒絕零假設，認為至少有兩組頻數分布是有統計學差異的，下一步是確定哪兩組間有統計學差異。SAS在ANOVA中可以直接指定post hoc分析，但對於卡方檢驗，SAS沒有自帶程序可以直接處理，這期講一下有哪些比較實用的方法與代碼~

方法1：挨個做每兩組的卡方檢驗

        存在的問題是多重比較，會增大Type I Error，需要校正p值，比如用簡單但保守的Bonferroni method（假設有4組，則會做4*3/2=6次的兩兩比較，把p值*6，再跟0.05進行比較）。

        在查資料時，發現了一種將Tukey HSD和permutation相結合的方法，筆者介紹其統計學效能較高(連結：https://support.sas.com/resources/papers/proceedings14/1544-2014.pdf)。

        Tukey HSD方法是ANOVA中post hoc分析的一種方法，全稱是Tukey Honestly Significantly Difference，思想是利用q分布(studentarized range distribution)計算出HSD值，一個「可信的顯著差異」，當兩組均值差的絕對值大於HSD，可認為significant。

      由於卡方檢驗的'q分布'是未知的，所以採用permutation——置換，針對每次隨機置換原數據集而產生的新數據集做兩兩比較卡方檢驗，保留最大卡方值。重複上述過程很多很多次，如5000次來模擬組間無差異時的最大卡方值的分布。計算原組間卡方值在這個分布中右側的面積，即為近似的p值。

        資料中給出了兩個macro，第一個嵌套在第二個內，使用上述方法直接調用第二個macro即可。最後可以得到原始p value，bonferroni和Tukey校正後的p value。

       這裡想說的是第一個macro，%Pairwise_Chisq，即使不用上述方法，也可以利用這個宏來實現每兩組卡方檢驗的自動化。在使用這個宏前，需要對數據集進行調整。

        比如現有一數據集A，我們關注BMI和sport兩個分類變量，頻數分布如下(模擬數據)：     

       將一人一行的原始數據集A轉換成頻數數據集_A，建議在_A中用具體分類代替數字分類，後面比較的結果可以直觀看出是哪兩組，而不是1 vs 2。然後運行macro：

proc freq data=A noprint;
  table BMI*sport/out=_A;
run;
proc format;
  value bmi 1='underweight' 2='normal' 3='overweight' 4='obesity';
  value sport 1='light' 2='medium' 3='heavy' 4='violent';
run;
data _A;
    set _A(rename=(BMI=_BMI sport=_sport));
    length BMI sport $40.;
    BMI=put(_BMI,bmi.);
    sport=put(_sport,sport.);
run;

%pairwise_chisq(data=_A, group=bmi, outcome=sport, count=count);
        結果如下，針對這個結果，可以手動進行bonferroni校正或者用proc multtest校正最後一列的P_value。

方法2：構建poisson log-linear model
        poisson分布是針對頻數資料的常見分析方法，卡方檢驗的結果是可以通過構建Poisson對數線性模型的方法得到的。 比如對前面生成的頻數數據集_A，運行Proc genmod程序，這裡只關注主效應：
proc genmod data=_A;
    class bmi sport;
    model count=bmi sport/dist=poisson;
run;
        在擬合優度的結果中，Pearson卡方值與卡方檢驗中的卡方值一致，偏差與卡方檢驗中似然比卡方檢驗結果一致：        Proc freq卡方檢驗結果：
        做組間比較需要引入交互項bmi*sport，用lsmestimate語句進行聯合卡方檢驗(選項joint)。可以先用lsmeans語句查看交互項的順序，方便確認lsmestimate語句的係數排列，0不納入，兩個比較組的係數為1和-1，各係數和=0。
        比如第一個lsmestimate語句，normal組與obesity組運動頻率無差異，需要這兩組間heavy vs light， heavy vs medium，heavy vs violate無差異。故做聯合卡方檢驗，其餘同理。
proc genmod data=_A;
  class bmi sport;
  model count=bmi sport bmi*sport/ dist=poisson type3 wald;
  lsmeans bmi*sport;
  lsmestimate bmi*sport 'normal vs obesity' 1 -1 0 0 -1 1 0 0,
		 	'normal vs obesity' 1 0 -1 0 -1 0 1 0,
			'normal vs obesity' 1 0 0 -1 -1 0 0 1/joint(only);
  lsmestimate bmi*sport 'normal vs overweight' 1 -1 0 0 0 0 0 0 -1 1 0 0,
	                'normal vs overweight' 1 0 -1 0 0 0 0 0 -1 0 1 0,
	                'normal vs overweight' 1 0 0 -1 0 0 0 0 -1 0 0 1/joint(only);
  lsmestimate bmi*sport 'normal vs underweight' 1 -1 0 0 0 0 0 0 0 0 0 0 -1 1 0 0,
	                'normal vs underweigh' 1 0 -1 0 0 0 0 0 0 0 0 0 -1 0 1 0,
                        'normal vs underweigh' 1 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 1/joint(only);
  lsmestimate bmi*sport 'obesity vs overweight' 0 0 0 0 1 -1 0 0 -1 1 0 0,
	                'obesityl vs overweight' 0 0 0 0 1 0 -1 0 -1 0 1 0,
	                'obesity vs overweight' 0 0 0 0 1 0 0 -1 -1 0 0 1/joint(only);
  lsmestimate bmi*sport 'obesity vs underweight' 0 0 0 0 1 -1 0 0 0 0 0 0 -1 1 0 0,
		        'obesity vs underweight' 0 0 0 0 1 0 -1 0 0 0 0 0 -1 0 1 0,
		        'obesity vs underweight' 0 0 0 0 1 0 0 -1 0 0 0 0 -1 0 0 1/joint(only);
  lsmestimate bmi*sport 'overweight vs underweight' 0 0 0 0 0 0 0 0 1 -1 0 0 -1 1 0 0,
	                'overweight vs underweight' 0 0 0 0 0 0 0 0 1 0 -1 0 -1 0 1 0,
	                'overweight vs underweight' 0 0 0 0 0 0 0 0 1 0 0 -1 -1 0 0 1/joint(only);
run;
        給出的比較結果如下，可以對比這裡的卡方值，近似等於上面用%pairwise_chisq計算的兩組卡方檢驗值。
        使用這種方法，最後仍需校正p值，校正方法同上所述。