# Ch2 countingency
這就是列聯表，也是在後續會用到的範例 :
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Myocardial\ Infarction$
| | Fatal attack | Nonfatal attack | No attack |
| ------- | ------------ | --------------- | --------- |
| Placebo | 18 | 171 | 10845 |
| Aspirin | 5 | 99 | 10933 |
---
交叉列表，也稱做列聯表(countingency table / cross table)
通常交叉列表都是一個 2 維表格，行、列分別代表不同的變數，單元格中是 2 個變數的組合的觀察次數，目的是為了分析這些變數的相互關係
以範例來舉例 :
我們會希望找到 Aspirin 跟 MI 的關聯性
接著，我們要設定一些參數 :
* $Y$ : 反映變數，有 $J$ 類類別
* $X$ : 解釋變數，有 $I$ 類類別
$\Rightarrow$ 所以分類的總和是 $IJ$
$i$ 所表示的就是 **列**
$j$ 所表示的就是 **行**
$n$ 所表示的就是 **人數**
然後因為我們希望可以看到 2 個類別之間的關係，因此我們透過機率來看出他們的關聯
我們設定機率為 $\pi_{ij} = P[X = i,\ Y = j]$ (這是為了看出聯合的機率)
並且為了看出 $X$ 跟 $Y$ 的關聯，所以要把變數分開來看，因此設定邊際分配的符號
$\pi_{i+} = \Sigma_j \pi_{ij}$
$\pi_{+j} = \Sigma_i \pi_{ij}$
$\Rightarrow \Sigma_i \pi_{i+} = \Sigma_j \pi_{+j} = \Sigma_i\Sigma_j \pi_{ij} = 1$
:::danger
樣本估計
$\pi_{ij}$ 算出來的就是樣本的比例
$p_{ij} = \cfrac{n_{ij}}{n}$ ($p$ 是估計)
$\Rightarrow p_{j|i} = \cfrac{p_{ij}}{p_{i+}} = \cfrac{n_{ij}}{n_{i+}}$
:::
:::danger
條件推論
$\pi_{ij} = P[Y=j|X=i] = \cfrac{\pi_{ij}}{\pi_{i+}}$
目的是為了知道在給定 $j$ 的情況下 $i$ 的改變情況
:::
---
接著我們要先做假設 :
假設 $X$ 和 $Y$ 獨立，我們定義為 :
$\pi_{ij} = \pi_{i+} \pi_{+j}$
因此變成條件機率時會變成 :
$\pi_{j|i} = \cfrac{\pi_{ij}}{\pi_{i+}} = \cfrac{\pi_{i+} \cdot \pi_{+j}}{\pi_{i+}} = \pi_{j+}$
::: success
小補充(獨立&同質的區別)(在統計層面不太影響)
* 獨立性檢定
也就是 $X$ 和 $Y$ 相互獨立
另外，在取樣時 $X$ 跟 $Y$ 都是隨機變數
* 同質性檢定
條件分布不會因為 $X$ 而有不同
另外，在取樣時 $X$ 是控制變數(固定)，$Y$ 是隨機變數
:::
---
**study design**
接下來我們要來找出取樣方法(但其實後續計算並不會著重於這部分) :
會根據不同的分配而有不同的取樣方法(下方的 $m$ 表示固定數量)
* 固定期間的取樣 : **Poisson**
| $n_{11}$ | $n_{12}$ |
| -------- |:-------- |
| $n_{21}$ | $n_{22}$ |
4 個都是獨立的 poisson
* 分成 2 類都取樣固定數量 : **Binomial**
| $n_{11}$ | $n_{12}$ |
| -------- |:-------- |
| $n_{21}$ | $n_{22}$ |
| $m$ | $m$ |
會有 2 個獨立的 binomial
*在實務應用上常用到(尤其在醫學界，又被稱做 case-control study)*
* 總共取樣的數量是固定的 : **Multinomial**
| $n_{11}$ | $n_{12}$ | |
| -------- |:-------- |:--- |
| $n_{21}$ | $n_{22}$ | |
| | | $m$ |
把符合條件的分別放入格子內
::: success
小補充
Poisson、Binomial、Multinomail 都是取後可以放回的(母體數 N 通常較大)
而 Hypergeometric 是取後**不**放回(母體數 N 通常較小)
:::
---
* retrospective study--回朔性研究
**(主要控制的是直排的行)**
事件已經發生
* prospective study--前瞻性研究
**(主要控制的是橫排的列)**
事件尚未發生，會有研究設計
包含了 :
- clinical trails : 臨床試驗，對不同受試者採取 2 種不同的方法，持續追蹤觀察是否有差別
- cohort study : 因為對某一群體感興趣，持續追蹤觀察是否跟其他的不一樣
---
我們接下來要設定一些參數 :
如果只有 2 組的話 :
$\pi_{2|i} = 1 - \pi_{1|i}$ 可以被簡化為 $\pi_{1|i} = \pi_i$ (因為發生 2 的情況會被 1 所決定)
想要比較 2 者發生的機率，我們會有一些方法 :
**第一種方法--相減 :**
定義 : $\pi_2 - \pi_1$
得到的值 : 界在 $-1 \sim 1$ 之間
假設 : $\pi_1 = \pi_2$ 代表 $Y$ 獨立
:::info
但這個方法如果差距很小的時候就可能沒辦法完整呈現出他們之間的差異
:::
**第二種方法--相對風險(RR) :**
相對風險(Relative Risk)，簡稱為 RR
定義 : $RR = \cfrac{\pi_1}{\pi_2}$
得到的值 : $RR > 0$
假設 : $RR = 1$ 代表 $Y$ 獨立
直接觀察 $RR$ 就能知道 $\pi_1$ 跟 $\pi_2$ 的大小關係
$RR < 1 \Rightarrow \pi_1 < \pi_2$
$RR > 1 \Rightarrow \pi_1 > \pi_2$
:::info
但這個方法沒辦法直接透過現有資料做出估算，需要通過額外的資訊來得知估算值
(因此只能使用在前瞻性的研究)
而且會受到 study design 的影響
:::
**第三種方法--勝算率(OR) :**
勝算率(Odds Ratio)，簡稱為 OR
假設成功的機率是 $\pi$，勝算(Odds)是 $\Omega = \cfrac{\pi}{1 - \pi}$
那我們就可以得到勝算率
定義 : $OR = \theta = \cfrac{\Omega_1}{\Omega_2} = \cfrac{\frac{\pi_1}{1-\pi_1}}{\frac{\pi_2}{1-\pi_2}}$
得到的值 :
$1 < \theta < \infty \Rightarrow \pi_1 > \pi_2 \Rightarrow$ row_1 相較於 row_2 較容易**成功**
$0 < \theta < 1 \Rightarrow \pi_1 < \pi_2 \Rightarrow$ row_1 相較於 row_2 較容易**失敗**
假設 : $\theta = 1$ 代表 $X$ 跟 $Y$ 相互獨立
如果現在是 2 $\times$ 2 的列表的話則定義會變成 : $\theta = \cfrac{\Omega_1}{\Omega_2} = \cfrac{{\pi_{11}}/{\pi_{12}}}{{\pi_{21}}/{\pi_{22}}} = \cfrac{{\pi_{11}} \cdot {\pi_{22}}}{{\pi_{21}} \cdot {\pi_{12}}}$
證明 : $\because \pi_1 = \cfrac{\pi_{11}}{\pi_{1+}},\ \pi_2 = \cfrac{\pi_{21}}{\pi_{2+}}$
$\Rightarrow \Omega_1 = \cfrac{\pi_1}{1 - \pi_1} = \cfrac{{\pi_{11}}/{\pi_{1+}}}{1 - ({\pi_{11}}/{\pi_{1+})}} = \cfrac{\cfrac{\pi_{11}}{\pi_{1+}}}{\cfrac{\pi_{11} + \pi_{12} -{\pi_{11}}}{\pi_{1+}}} = \cfrac{\pi_{11}}{\pi_{12}},\ \Omega_2 \ 同理 = \cfrac{\pi_{21}}{\pi_{22}}$
同時也能直接觀察 $\theta$ 知道 $\Omega_1$ 跟 $\Omega_2$ 的大小關係
$\theta < 1 \Rightarrow \Omega_1 < \Omega_2$
$\theta > 1 \Rightarrow \Omega_1 > \Omega_2$
除此之外，也能知道當 $\theta > 1$ 時表示 $\pi_{11}$ 或 $\pi_{22}$ 是大的，也就表示不論怎樣 $\pi_{11}$ 這格的勝算都比較高
$\hat{\theta} = \cfrac{n_{11} \cdot n_{22}}{n_{12} \cdot n_{21}}$
:::info
不會受到 study design 的影響，不論回朔性或前瞻性的研究都可使用
:::
:::danger
小重點--RR 跟 OR 的關係
$\theta = \cfrac{\pi_1 / (1 - \pi_1)}{\pi_2 / (1 - \pi_2)} = \cfrac{1 - \pi_2}{1 - \pi_1} \times RR$
因此如果 $\pi_i \rightarrow 0$ 的話(prevalence 小)，$\theta$ 會跟 $RR$ 很接近
並且如果 $RR>1 \Leftrightarrow \theta >1, \ \because 方向一致，只有數值差異而已$
:::
:::warning
解釋技巧 :
* 假設 $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$)
:::
---
**信賴區間(CI)**
$\hat{\pi} \ 或 \ \hat{p} = \cfrac{n_1}{n} \sim Bin(n,p)$
by CLT(中央極限定理) $Bin(n,p) = \Sigma^n_{i=1} y_i = n_1,\ when \ y_i \sim Bernoulli(p)$
$\Rightarrow \cfrac{\hat{p} - p}{\sqrt{p(1-p)/n}} \stackrel{\text{d}}{\longrightarrow} N(0,1)$
因為我們不知道 $p$ 只知道 $\hat{p}$ 所以要使用到 Wald 的信賴區間法
**Wald CI**
$P[\cfrac{|{\hat{p} - p}|}{\sqrt{\hat{p}(1-\hat{p})/n}} < z_{\alpha/2}] = 1 - \alpha$
$\hat{p} - z_{\alpha/2} \cdot \sqrt{\cfrac{\hat{p}(1-\hat{p})}{n}} \le p \le \hat{p} + z_{\alpha/2} \cdot \sqrt{\cfrac{\hat{p}(1-\hat{p})}{n}}$
$p_1 - p_2$ (風險差) 二組獨立
$\cfrac{\hat{p_1} - \hat{p_2} - (p_1 - p_2)}{\sqrt{\cfrac{p_1 (1-p_1)}{n_1} + \cfrac{p_2 (1-p_2)}{n_2}}} \stackrel{\text{d}}{\longrightarrow} N(0,1)$
(如果分母的 $p$ 都改成 $\hat{p}$ 的話那就是 Wald，另外分母本身也等同於 SE(Standard Error))
$\therefore \Rightarrow p_1 - p_2 \pm z_{\alpha/2}(SE)$
:::danger
**RR 的信賴區間** :
因為 $(log(\pi_1) - log(\pi_2)) \pm z_{\alpha/2} \sigma (log(RR))$ 是以 $log(RR)$ 所計算的，所以最後只要取 $exponential$ 就可以了
$\Rightarrow e^{(log(\pi_1) - log(\pi_2)) \pm z_{\alpha/2} \sigma (log(RR))}$
**OR 的信賴區間** :
證明過程跟 $RR$ 類似
因為 $log(\theta) \pm z_{\alpha/2} (SE)$ 是以 $log(\theta)$ 所計算的，所以最後只要取 $exponential$ 就可以了
$\Rightarrow e^{log(\theta) \pm z_{\alpha/2} (SE)}$
:::
---
這是證明，參考用
:::spoiler
**相對風險 (RR) 的漸進分配**
by CLT : $\cfrac{\hat{\pi_1} - \pi_1}{\sqrt{\pi_1 (1-\pi_1)/n}} \stackrel{\text{a}}{\longrightarrow} N(0, 1)$
$\Rightarrow 因為 RR 不好算，所以我們可以先取對數$
$log(\cfrac{\hat{\pi_1}}{\hat{\pi_2}}) = log(\hat{\pi_1}) - log(\hat{\pi_2})$
現在已知 $\hat{\pi_i} \stackrel{\text{a}}{\longrightarrow} N(\pi_i, \cfrac{\pi_i(1-\pi_i)}{n_{i+}})$
接著透過下方所介紹的 Delta Method，把 $log(\pi_i)$ 當成 $f(\theta)$
會得到 $log(\hat{\pi_i}) \stackrel{\text{a}}{\longrightarrow} N \left ( log(\pi_i), \cfrac{\pi_i(1-\pi_i)}{n_{i+}} \times \cfrac{1}{\pi_i^2} \right) = N \left ( log(\pi_i), \cfrac{1-\pi_i}{n_{i+} \pi_i} \right)$
$\therefore \Rightarrow \sigma(log \ \hat{RR}) \approx \sqrt{\cfrac{1-\pi_1}{n_{1+} \pi_1} + \cfrac{2-\pi_2}{n_{2+} \pi_2}}$
**Delta Method**
令 $T(\#)$ 為一個統計量($\theta$)
我們可以知道這個統計量為大樣本時會趨近於常態 $T(\#) \stackrel{\text{d}}{\longrightarrow} N(\theta, \sigma^2)$($\sigma$ 是已知的)
現在還有另一個 **可微(differentiable)** 函數 $f(x)$
我們希望找到這個函數的分配 $f(T(\#))$ 的分配
$f(T(\#)) \stackrel{\text{d}}{\longrightarrow} N(f(\theta), f^{\prime}(\theta)^2 \cdot \sigma^2)$
而 Delta Method 是由 **one-term Taylor expansion** 來的
$f(x) = f(\theta) + f^{\prime}(\theta)(x - \theta) + R^{\prime}_m$
我們可以透過此 **泰勒展開式** 得出上方的 $f^{\prime}(\theta)^2$
$Var(f(T(\#))) = Var(f^{\prime}(\theta)(T(\#)- \theta)) = f^{\prime}(\theta)^2 \underbrace{Var(T(\#))}_{\sigma^2}$
:::
---
sampling zeros :
在取樣時因為樣本數抽得太少所以沒有抽到導致有一格為 0
為了避免計算時 0 被放在分母所以我們會做修正
像是把 SE 當中的所有格子都加上 $0.5$
:::info
小補充 :
* sampling zeros : 因為樣本數抽得太少所以沒抽到
* structure zeros : 該格一定不會有數據(像是自主免疫之類的調查就會出現)
:::
## 程式碼實作
**我們使用 SAS 來執行**
首先，設定我們會用到的 data :
```
data class2a;
input trt status count;
cards;
0 1 189
0 0 10845
1 1 104
1 0 10933
run;
```
接下來這是把我們的資料作成交叉列表(contingency table)
```
proc freq data = class2a;
weight count; /*weight : count 做加權(會把每列都當一個人)*/
table trt*status; /*table 後面接上的"*"是把列表做成交叉列表(列*行)*/
run;
```
然後為了把表做成符合我們要求，所以加上 `order` 讓資料順序是我們給定的
```
proc freq data = class2a order = data; /*按照給定的資料順序呈現*/
weight count; /*加權*/
table trt*status; /*列*行*/
run;
```
小補充
(如果位置是錯誤的，可以先做排序(位置通常會是 上方:1 0，左方:0 1))
```
proc sort data = class2a;
by decending status; /*因為位置是錯誤的，要修正(通常會是 上面:1 0 左邊 : 0 1)*/
```
接著，更進一步可以計算出風險差值
```
proc freq data = class2a order = data; /*按照給定的資料順序呈現*/
weight count; /*加權*/
table trt*status/ riskdiff; /*列*行*/ /*riskdiff : 計算風險差值*/
run;
```
這邊是把風險差值改成相對風險
```
proc freq data = class2a order = data; /*order : 按照我們給的資料順序呈現*/
weight count; /*加權*/
table trt*status/ relrisk; /*列*行*/ /*relrisk : 計算相對風險*/
run;
```
這邊是改成卡方檢定
```
proc freq data = class2a order = data; /*order : 按照我們給的資料順序呈現*/
weight count; /*加權*/
table trt*status/ chisq; /*列*行*/ /*chisq : 計算卡方檢定*/
run;
```
這邊是想要看 $\mu_{ij}$ test 的值
```
proc freq data = class2a order = data; /*order : 按照我們給的資料順序呈現*/
weight count; /*加權*/
table trt*status/ expected deviation; /*列*行*/
/*expected : expected cell count*/
/*deviation : n_{ij} - \mu_{ij}*/
run;
```
這邊是想要看 Fisher Exact test
```
proc freq data = class2a order = data; /*order : 按照我們給的資料順序呈現*/
weight count; /*加權*/
table trt*status/ exact; /*列*行 exact:Fisher Exact test*/
run;
```
Sign in with Wallet
Connect another wallet