# 基礎概念
* 做假設的技巧 :
幾乎做出的假設都是 $\Rightarrow$
$H_0$ : ...<表格-列>... 與 ...<表格-欄>... 相互獨立
$H_1$ : ...<表格-列>... 與 ...<表格-欄>... 相關
* compute risk different :
在 `sas` 中使用 `proc freq` 裡面的 `riskdiff`
* compute RR :
在 `sas` 中使用 `proc freq` 裡面的 `relrisk`
看表格最下面跑出來的 **相對風險(欄n)**
(n) 是代表我們想觀察的那一項在哪一欄
* compute OR :
在 `sas` 中使用 `proc freq` 裡面的 `relrisk`
看表格最下面跑出來的 **勝算比**
* 使用甚麼 test 去看是否符合假設 :
通常都是使用 $\chi^2$ 或是 $G^2$
這 2 個都是去測資料之間是否**獨立**
但如果資料的表格是比 2 $\times$ 2 還要大的就要改成使用 Pearson $\chi^2$
如果是使用 Pearson $\chi^2$ 時，程式碼還是會使用 `chisq` 來檢測
* 不同 test 的意義 :
$\begin{cases} \text{chisq : 是否相關} \\ \text{trend : 是否線性相關} \\ \text{monotone : 是否具有單調性(在 sas 中要打 measures)} \\ \text{BD(Breslow-Day) : 3 個變數時，其中 1 個是否會影響另外 2 個} \end{cases}$
:::danger
需要注意跑出來的報表我們想解釋的是哪一欄
如果不是我們要的需要使用 `proc sort`
把我們想解釋的放在左上角
**不然跑出來的 risk diff / RR / OR 會出問題**
:::
:::warning
解釋跑出來的值以及信賴區間 :
首先需要區分是 (risk diff) 還是 (RR 或 OR)
如果是 risk diff : 包含 $0$ 的話代表沒有顯著差異 / 不包含的話就是有顯著差異
如果是 RR 或 OR : 包含 $1$ 的話代表 2 者之間沒有顯著差異 / 不包含的話就是有顯著差異
:::
:::danger
RR 跟 OR 的解釋技巧 :
* 假設 $RR = \cfrac{\pi_1 \Rightarrow P}{\pi_2 \Rightarrow A} = 1.82$ : P 會 ...<題目>... 的機率是 A 的 1.82 倍
($\because \pi_1 = RR \times \pi_2$)
* 假設 $OR = \theta = \cfrac{\Omega_1 \Rightarrow P}{\Omega_2 \Rightarrow A} = 1.83$ : 使用 P 得到 ...<題目>... 的勝算是使用 A 的 1.83 倍
* ($\because \Omega_1 = \theta \times \Omega_2$)
:::
::: info
如果取樣不平均的話 :
total 得出來的結論跟分開討論會有差別
:::
# 程式碼實作
**我們使用 SAS 來執行**
## 風險
設定資料檔
```sas=
data q1_1;
input alc status count;
cards;
0 0 17066
0 1 48
1 0 15415
1 1 38
run;
```
把資料做成交叉列表(contingency table)(不常用到)
```sas=
proc freq data = q1_1;
weight count; /*weight : count 做加權(會把每列都當一個人)*/
table alc*status; /*table 後面接上的"*"是把列表做成交叉列表(列*行)*/
run;
```
把資料排序成我們好解釋的
```sas=
proc sort data = q1_1;
by decending status;
/*因為位置是錯誤的，要修正(通常會是 上面:1 0 左邊 : 0 1)*/
/*小補充 (如果位置是錯誤的，可以先做排序)*/
/*(位置通常會是 上方:1 0，左方:0 1)--這是我們想看的解釋變數為 1 的時候*/
```
計算風險差值
```sas=
proc freq data = q1_1 order = data; /*按照給定的資料順序呈現*/
weight count; /*加權*/
table alc*status/ riskdiff; /*列*行*/ /*riskdiff : 計算風險差值*/
run;
```
計算相對風險跟勝算比
```sas=
proc freq data = q1_1 order = data;
weight count;
table alc*status/ relrisk; /*列*行*/ /*relrisk : 計算相對風險*/
run;
```
## 獨立性檢定
設定資料檔
```sas=
data q2_1;
input snoring hd count;
cards;
0 1 24
0 0 1355
2 1 35
2 0 603
4 1 21
4 0 192
5 1 30
5 0 224
run;
```
測 2 者是否獨立(做卡方獨立檢定)
```sas=
proc freq data = q2_1 order = data;
weight count; /*加權*/
tables snoring*hd/chisq; /*做卡方檢定(測是否相關)*/
run;
```
看 Fisher Exact test (也是做獨立檢定，但樣本數比較少的時候)
```sas=
proc freq data = q2_1 order = data;
weight count;
table snoring*hd/ exact; /*列*行 exact:Fisher Exact test*/
run;
```
測 2 者是否**線性**獨立
```sas=
proc freq data = q2_1 order = data;
weight count;
tables snoring*hd/trend; /*做 trend test (測是否線性相關)*/
run;
```
測 2 者之間的關係是否為**單調**
```sas=
proc freq data = q2_1 order = data;
weight count;
tables snoring*hd/measures;
/*做 monotone test (測是否單調)*/
/*(要看 gamma)*/
run;
```
看 $\mu_{ij}$ test 的值
```sas=
proc freq data = q2_1 order = data;
weight count;
table snoring*hd/ expected deviation; /*列*行*/
/*expected : expected cell count*/
/*deviation : n_{ij} - \mu_{ij}*/
run;
```
## 合併與否
設定資料檔(**important**)
```sas=
/*哪些項目可以合併*/
data q3_1;
input color satellite count;
/* 設定合併類別的範例*/
/*先比 1 2*/
if color = 1 then color12 = 1;
else if color = 2 then color12 = 2;
/*再比 3 4*/
if color = 3 then color34 = 1;
else if color = 4 then color34 = 2;
/*最後綜合 1 2 跟 3 4*/
if color in ( 1 2 ) then color1234 = 1;
else if color in ( 3 4 ) then color1234 = 2;
Cards;
1 1 9
1 0 3
2 1 69
2 0 26
3 1 26
3 0 18
4 1 7
4 0 15
run;
/*之所以這樣分是因為要做有意義的分類*/
/*(題目的前 2 個都有 medium，後 2 個都有 dark)*/
```
```sas=
/*分完之後要做檢定*/
proc freq data = q3_1;
weight count;
tables (color12 color34 color1234)*satellite/chisq; /*括號可以同時做很多個(分配律)*/
run;
/*檢定結果如果不顯著才可以合併 (也就是 p > 0.05)*/
```
## 3 個變數時，是否互相影響
設定資料檔
```sas=
data q4_1;
input race age sat count;
cards;
0 1 0 109
0 1 1 57
0 2 0 108
0 2 1 47
0 3 0 141
0 3 1 41
1 1 0 31
1 1 1 13
1 2 0 7
1 2 1 5
1 3 0 10
1 3 1 3
run;
```
測是否會有交互作用
```sas=
proc freq data = q4_1 order = data;
weight count;
table age * race*sat/ cmh;
/*z 變數 * x 變數 * y 變數*/ /*順序很重要*/
/*測 z 變數是否會影響到 x 跟 y 之間的關係*/
run;
/*看到最下面的均齊性是 0.5979 在 0.05 底下不拒絕 H0 --> 會有共同勝算*/
/*往上會看到 CMH 的統計量(上方是給定條件--對照 age)*/
/*非零 : 檢測線性相關*/
/*列 : 假設 X 是類別，Y 是連續(跟 T / Wilcoxon 一樣)*/
/*一般關聯 : X Y 的關聯，沒有特別設限，把 X Y 假設是 nominal 的情況*/
/*都小於0.05 --> 拒絕H0 --> 條件不獨立*/
/*接著看下方的勝算比 : 就可以給出說明 -->*/
/*給定被害人的種族是白人的情況下，被告的種族是白人，被判死刑的勝算是黑人的 0.9730 倍*/
```
## 建立模型(model)
設定資料檔
```sas=
data q5_1;
input snoring hd count;
cards;
0 1 24
0 0 1355
2 1 35
2 0 603
4 1 21
4 0 192
5 1 30
5 0 224
run;
```
做 linear probability model
最後跑出的結果須寫成 : $\hat{\pi}(x) = \text{Intercept} + \text{第二項 }* x$
```sas=
/*linear probability model*/
proc genmod data = q5_1 desc; /*decending 的縮寫*/
/*因為有興趣的是有心臟病的機率*/
/*差別的部分會出現在回應概況--hd的排序值第一個會是 hd = 1*/
/*原始寫法是要先 sort 過才能執行*/
weight count; /*需分清楚是否做加權*/
model hd = snoring/ link = identity dist = bin;
/*link 選的是機率線性模型 identity*/
/*第二個選的是分配 binary*/
/*左邊是反應變數 右邊是解釋變數*/
run;
/*\hat{pi}(x) = 0.0172 + 0.0198 x*/
```
做 Logistic regression
最後跑出的結果須寫成 : $logit(\hat{\pi}(x)) = log \left( \cfrac{\hat{\pi}(x)}{1- \hat{\pi}(x)} \right) = \text{Intercept} + \text{第二項 }* x$
```sas=
/*Logistic regression*/
proc genmod data = q5_1 desc;
weight count;
model hd = snoring/ link = logit dist = bin;
run;
/*log(\hat{pi}(x) ) / (1- \hat{pi}(x) )= -3.8662 + 0.3973 x*/
/*x 前面係數的正負在 Logistic 跟 linear 應該要是相同的*/
```
做 Probit model
最後跑出的結果須寫成 : $\Phi(\hat{\pi}(x) ) = \text{Intercept} + \text{第二項 }* x$
```sas=
/*Probit model*/
proc genmod data = q5_1 desc;
weight count;
model hd = snoring/ link = probit dist = bin;
run;
/*\Phi(\hat{pi}(x) )= -2.06 + 0.19 x*/
```
做 Log linear model 或 Poisson model
最後跑出的結果須寫成 : $log(\hat{\mu}(x)) = -3.3 + 0.164 x$ 或 $e^ \hat{\beta} = exp(0.164)$
```sas=
/*Log linear model or Poisson model*/
proc genmod data = q5_1 desc;
/*不設權重，因為一筆就是一隻*/
model hd = snoring/ link = log dist = pois;
run;
/*log(\hat{\mu}(x)) = -3.3 + 0.164 x*/
/*e^ \hat{\beta} = exp(0.164)*/
```
## 檢測過度離散
設定資料檔
```sas=
data q6_1;
infile "C:\Users\ACER\Desktop\類別 sas 檔\crab.dat";
input color spine width sat weight;
y = (sat>0);
run;
```
從簡單的方法檢測是否過度離散(看平均值跟變異數)
```sas=
/*看平均值跟變異數是否相等 --> 從這個表可以看出有過度離散*/
proc means data = q6_1 n sum mean var;
class width;
var sat;
run;
```
用 `genmod` 看是否過度離散
從配適度的 `值/DF`，我們希望接近 1，如果不是的話就是過度離散 $\Rightarrow$ 需要做調整
```sas=
proc genmod data = q6_1;
model sat = width / link = log dist = pois
aggregate;
run;
/*aggregate 整合 likelihood (因為這個資料是連續的所以不太需要用到)*/
/*檢視是否過度離散(Examine overdispersion)*/
/*看配適度的 值/DF ，預期要接近 1，超過的話代表就是過度離散*/
/*這樣就要做調整*/
```
做調整
使用 `scale` 的 pearson $\chi^2$ test 做調整
可以從最下方的縮放看出確實有調整過
```sas=
/*做調整*/
proc genmod data = q6_1;
model sat = width / link = log dist = pois
aggregate scale = pearson;
/*scale 是做卡方獨立檢定的時候會用到的規模*/
/*因為樣本比較小所以要選擇 pearson test*/
run;
```
## 同質性分析、過度離散調整
設定資料檔
```sas=
data q7_1;
infile "C:\Users\ACER\Desktop\類別 sas 檔\crab.dat";
input color spine width sat weight;
color = color -1;
y = (sat>0);
run;
```
最簡單看是不是 monotone(單調) 的方法
```sas=
proc sgplot data = q7_1;
vline width/ response = y stat = mean markers;
run;
```
跟 `genmod` 差不多的程序
```sas=
proc logistic data = q7_1 desc; /*因為有興趣的事件是 1 所以使用 desc*/
model y = width;
/*因為是連續的所以不用調其他東西*/
run;
```
把資料檔從連續型的改成二元的
```sas=
data q7_2;
set q7_1;
if color in ( 1 2 ) then light = 1;
else if color in ( 3 4 ) then light = 0;
run;
```
然後看他的分析(預設 `param = effect`)
參考係數是(1, -1)
```sas=
/*Default : param = effect*/
proc logistic data = q7_2 desc; /*因為有興趣的事件是 1*/
class light;
/*因為想看的是 light 所以要設 class*/
model y = light;
run;
```
把預設的 `param = effect` 改成 `ref`(通常都會使用這個)
預設的參考組別是 (last)
參考係數是(1, 0)
```sas=
/*Reference : param = ref default : last*/
proc logistic data = q7_2 desc;
class light/ param = ref;
model y = light;
run;
```
把預設的參考係數 `param = effect` 改成 `ref`
把預設的參考組別 `ref = last` 改成 `first`
算出的 `勝算比點估計值` 會跟參考組別為 `last` 的互為倒數
```sas=
/*Reference : param = ref default : first*/
proc logistic data = q7_2 desc;
class light/ param = ref ref = first;
model y = light;
run;
/*淺色的 horsehoe crab 有 satellite 的勝算是深色的 2.689 倍*/
```
做同質性分析(看主成分)
```sas=
/*Homogeneous association -- main effect*/
proc logistic data = q7_2 desc;
class light spine/ param = ref ref = first;
model y = light spine;
run;
```
做異質性分析(看主成分)(然後看交互作用)
```sas=
/*Hetergeneous association -- add interaction*/
proc logistic data = q7_2 desc;
class light spine/ param = ref ref = first;
model y = light spine light * spine;
run;
```
因為從上面的連結性檢定跟全域檢定看出 2 個得出的結果不一致
因此跑出表格來看是不是係格的問題
我們期望係格都很多或都很少，而不是有多有少
接著我們看過度離散
如果偏差跟 pearson 都大於 0.05 代表不拒絕 $H_0$，沒有過度離散
```sas=
/*Overdispersion*/
proc logistic data = q7_2 desc;
class color spine/ param = ref ref = first;
model y = color spine/ aggregate scale = none;
run;
```
假設有過度離散需要做調整
把 `scale = none` 改成 `p` 就是做 pearson 調整
```sas=
proc logistic data = class9a desc;
class color spine/ param = ref ref = first;
model y = color spine/ aggregate scale = p;
/*Pearson chi-square adjustment*/
run;
```
# 程式碼的小筆記專區
* `weight` 的使用時機:
如果是原始資料的話就不需要加權，反之，其他情況都會需要
* 要記得加上 `order = data` 因為大部分的情況都會希望用我們給定的排序
* 記得判別我們想要解釋的變數在哪裏
如果位置是錯誤的記得要 `proc sort...by decending` 一下
至於怎麼判別位置 $\Rightarrow$ 通常會是 `上面:1 0 左邊 : 0 1` ，這樣的情況是我們想看的解釋變數為 1 的時候
<font color="#f00">左上方的那個就是我們想看的</font>
* 計算風險
$\begin{cases} \text{riskdiff} \Rightarrow 風險差值 \\ \text{relrisk} \Rightarrow \text{RR--相對風險 或 OR--勝算比} \end{cases}$
RR 就看報表中 **相對風險** 那一格，但因為還會分(欄1)(欄2)
通常已經整理好的話就會看(欄1)，也就是我們想知道的那一項
OR 就看報表中 **勝算比** 那一格
* 風險的解釋
首先看他們的信賴區間的值
風險差值 : 包含 $0$ 的話代表沒有顯著差異 / 不包含的話就是有顯著差異
RR 或 OR : 包含 $1$ 的話代表 2 者之間沒有顯著差異 / 不包含的話就是有顯著差異
然後下結論的時候 $\Rightarrow$
RR : ... 會是 ... 的 $p$ 倍
OR : ... 會 ... 的勝算會是 ... 的 $p$ 倍
* 測相關性
首先會先假設 :
$H_0$ : ...<表格-列>... 與 ...<表格-欄>... 相互獨立 / 無線性關聯
$H_1$ : ...<表格-列>... 與 ...<表格-欄>... 相關 / 有線性關聯
接著就能跑報表
$\begin{cases}
\text{chisq} \Rightarrow \text{是否獨立, p < 0.05 代表 reject }H_0 \text{, 2 變數 有關聯 / 不獨立} \\
\text{trend} \Rightarrow \text{是否線性相關, p < 0.05 代表 reject }H_0 \text{, 2 變數 有線性相關, 列百分比會隨 x 上升/下降 而 上升/下降} \\
\text{measures} \Rightarrow \text{是否單調, gamma/ASE = z 值，如果 |z| > 1.96 代表 reject }H_0 \text{, 2 變數 是單調關係}
\end{cases}$
* 看哪些類別可以做合併
假設現在有 1, 2, 3, 4 這四個類別
先分成 1 2 跟 3 4
接著分成 1+2 3+4
然後把這些都拿去做卡方檢定
當 $p< 0.05$ 就代表他們顯著，類別之間有差異，不能做合併
因此只能合併 $p>0.05$ 的分類
* 3 個變數時，測是否有同質關係
使用到 `cmh`
並且看到最下方均齊性的卡方值 & p 值
如果 $p > 0.05$ 代表 $\text{not reject } H_0$
也就代表第 3 項不會干擾到另外 2 者的關係
這樣才能說明均齊性符合，可以看表說明內容
反之，如果 $\text{reject } H_0$ 就不能看表
因為代表第 3 項與其他 2 者有關係
* 做 model 時
注意題目是想做哪種 model 以及注意最後寫出來的等式左方
並且 $x$ 前面係數的正負在 Logistic 跟 linear 應該要是**相同**的
除此之外，注意題目是否有需要加入 `desc` 把我們想看得放在第一項
::: danger
齊一性 : 會有共同勝率 $\theta$，因此可以共同討論
異質性 : 需要分開解釋
**會使用 cmh 去測**
:::
* 更改參數設定不會影響勝算比
差別只是把二元的參數設成(1, -1)還是(0, 1)
* 為甚麼要看是否單調(monotone)
如果資料的呈現形式是單調的代表我們可以把資料視為連續型
* `cmh` 跑出來的假設定義
**非零** :
*檢測* : 檢測線性相關
*假設* : 假設 2 者皆為連續變數
**列** :
*檢測* : 檢測 X 是類別，Y 是連續(跟 T / Wilcoxon 一樣)
*假設* : 假設 X 為類別變數、Y 是連續變數
**一般關聯** :
*檢測* : 檢測 X Y 的關聯，沒有特別設限，把 X Y 假設是 nominal 的情況
*假設* : 沒有額外的假設
* cmh 比較模型好壞
通常會比較有交互作用跟沒有交互作用的
看 2 者模型配適統計值中的 AIC，值越小代表該模型越好
* cmh 解釋
像這種的解釋就是 :
light = 1 的勝算會是 light = 0 勝算的 $e^{-0.2657}$ 倍
light = 1 且 spine = 2 的勝算會是 light = 0 且 spine = 1 勝算的 $e^{0.6712}$ 倍
Sign in with Wallet
Connect another wallet