---
# System prepended metadata

title: catedata_finalreview
tags: [categorical data analysis]

---

# Outline  
* [Logistic regression](#羅吉斯迴歸)  
* [Multicatory logit model](#多元模型)  
* [Zero model](#零膨脹模型)  
* [Loglinear model for contingency table](#列聯表的對數線性模型)  
* [Models for matched pair](#成對樣本模型)  

之前的重點小提醒 :  
:::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$)  
:::  


# 羅吉斯迴歸  
主要是 2 元會使用到  
首先會先選模，有以下這些方式(利用 SAS 自動選模) :  
* backward  
* stepwise  
* forward  
這些方法就是把一些無關的特徵或雜訊值刪除掉  

:::success  
無關的特徵 :  
也就是不會影響到模型的建置  
我們可以從報表看出來  
當變量的 p 值大於我們設定的(sls)時就是無關的特徵，可以被刪除  
:::  

刪除之後我們再用剩下的變數做模型  
通常題目有說到 **establish logistic regression model** 就是要建立羅吉斯迴歸模型  


## 程式碼實作--test4  
設定資料檔  
```sas=
%let dirclass = C:\Users\ACER\Desktop\類別 sas 檔;
proc import datafile = "&dirclass\catedata112fq4.csv"
	out = q4 DBMS = CSV replace;
run;

data q4a;
set q4;
id = _n_;
run;
```

先用 backward selection(SLS = 0.2) 取出可解釋的變數  

```sas=
proc logistic data = q4a desc;
class change type field wind / param = ref ref = first; 
/*
class : 設置類別型的變數(分為 0 或 1)
ref = first : 依照規定設置 ref
*/
model good = week change elap30 type field wind distance / 
	selection = backward sls = 0.2;
run;
```
可以發現篩選出 `week` `change` `wind` `distance` 這 4 個變數  
其他的都被刪掉了  

接著我們再用這 4 個變數做一次模型(為了讓整個模型是乾淨的)  
並且我們希望可以看到每一筆資料的 `influence` 來決定是否保留該筆資料，  
因為有些資料的殘差太大(跟其他的差距過大)會讓模型配飾不準確  

```sas=
proc logistic data = q4a desc;
class change wind / param = ref ref = first;
model good = week change wind distance / influence;
id id;
ods output influence = inf; 
/*把 influence 的影響輸出到新的檔案中方便我們使用*/
run;
```

然後我們要找出 Pearson residual > 6 以及 deleted Pearson residual > 30  
這些都是殘差過大的，之後需要被刪除  

```sas=
proc print data = inf;
where difchisq > 30 or abs(reschi) > 6;
/*這幾筆資料殘差太大(卡方值很大)*/ /*所以要刪除*/
run;
/*會出現這幾個 ID : 224 875 1084 219 251 383 515 567 604 1227 1324 1321*/
```

接著做出沒有包含這幾個 ID 的模型  

```sas=
proc logistic data = q4a desc;
class change wind / param = ref ref = first;
model good = week change wind distance / aggregate
	scale = none lackfit;
/*
aggregate : 合併組別，將分類變量的效果組合，以獲得每個分類變量的總體影響
scale = none : 不進行縮放
lackfit : 在唯一設定檔中，有可能有些組別樣本數太多，有些組別樣本數太少，避免造成檢定出錯
另外，最重要的，lackfit 會產生 HL test--因為這題有 2 個連續變數，不確定分布
*/
id id;
ods output influence = inf;
where id ^in (224 875 1084 219 251 383 515 567 604 1227 1324 1321);
run;
/*
可以看到卡方值一個大於 0.05 一個小於0.001()
可能觀測值相異程度不平均有過度離散
要調整 standard error
用 HL看
可以看到 HL = 0.0036 -->reject H_0，樣本有過度離散
*/
```

可以發現做出來的報表 p 值一個大於 0.05、另一個 < 0.0001  
![image](https://hackmd.io/_uploads/HJLMa9G_a.png)  
代表觀測值的相異程度不平均，在大樣本的卡方時會有問題  
需要調整 standard error  
但我們還是可以從 HL test 得知，他的 p 值 < 0.05 $\Rightarrow$ reject $H_0$  
![image](https://hackmd.io/_uploads/By66gjzdp.png)  
代表應該有過度離散的情況  

另外，在解釋的部分，先從這裡看他的正負知道他們的關係是正向還是負向影響  
![image](https://hackmd.io/_uploads/S1Hy7oM_6.png)  
可以先做出模型(logistic regression model) :  
$\begin{split} 
logit(\hat{\pi(x)}) &= \hat{\beta_0} + \hat{\beta_1}x_1 + \hat{\beta_2}x_2 + ... \\ 
&= 7.0858 -0.0297 * week - 0.5309 * (change=1) - 0.6753 * (wind=1) - 0.1312 * (distance=1)
\end{split}$  
然後也可以發現係數都是負的，代表都為負向影響  
所以做解釋的時候都要以負向影響來解釋  
![image](https://hackmd.io/_uploads/Hkdf7sGdp.png)  
所以當有改變(change = 1)時，勝算會下降 0.588 倍  
有風(wind = 1)時，勝算會下降 0.509 倍  

再來，我們要來看配對  
![image](https://hackmd.io/_uploads/BkTMGN7_a.png)  
從 c 值可以看到配適效果還不錯(接近 1，遠離 0.5)  

然後我們想要知道未來我們要做預測該如何設定切點  
因此使用 `ctable` 來看  
```sas=
proc logistic data = q4a desc;
class change wind / param = ref ref = first;
model good = week change wind distance / ctable;
id id;
where id ^in (224 875 1084 219 251 383 515 567 604 1227 1324 1321);
run;
/*
可以看出切點大概要切在 sample propotion 大約在 0.88  附近 --> c
可以看到 敏感度 跟 特異性 幾乎相等
*/
```
![image](https://hackmd.io/_uploads/B10HXNX_T.png)  
![image](https://hackmd.io/_uploads/ByoPXEXO6.png)  
我們的切點是希望 `敏感度` 跟 `特異性` 幾乎相等  


補充 :  
:::success
```sas=
/*做出 ROC 曲線*/
proc logistic data = class11 desc plots = (ROC);
class light spine / param = ref;
model y = light width / ctable; /*classification table*/
run;
```
:::


## 程式碼實作--class  
假設現在資料檔有其中一個變數沒辦法估計的話，就要改成使用 exact 來看它的 $p$ 值  

設定資料檔  
```sas=
/*ch5-p36*/
data class12a;
input trt age length case sample;
cards;
0 0 0 0 385
0 0 1 5 233
0 1 0 3 789
0 1 1 47 1081
1 1 1 5 5
run;
```

接著做出羅吉斯迴歸  
```sas=
/*從結果可以發現我們沒辦法估計 trt
也可以從警告標語看到 : 
Warning: The information matrix is singular and thus the convergence is questionable. 
*/
proc logistic data = class12a;
class trt age length/param = ref ref = first;
model case / sample = trt age length;
exact trt age length; 
/*
因為沒辦法進行估計
所以多加一個精準條件式子的檢定
exact 後面要放入我們想檢測的變數
*/
run;
/*這個條件式的檢定就是condition on sufficient statistics，然後做統計推論*/
/*最重要的就是要看精準當中的 p 值來看是否顯著*/
```
![image](https://hackmd.io/_uploads/HyzZeJVuT.png)  
可以看到 trt 跟 length 的 $p < 0.05$，因此這 2 個變數是顯著的  


---


現在的資料檔是列聯表的形式  

設定資料檔  
```sas=
/*ch5-p42*/
data class12b;
input race month prom count;
cards;
1 7 1 0
1 7 0 7
0 7 1 4 
0 7 0 16
1 8 1 0
1 8 0 7
0 8 1 4
0 8 0 13
1 9 1 0
1 9 0 8
0 9 1 2
0 9 0 13
run;
```

接著做出羅吉斯迴歸  
```sas=
proc logistic data = class12b desc;
weight count;
class race month / param = ref ref = first;
model prom = race month;
run;
/*會看到>999.999 以及 標準誤差是53.XXX的情況*/
/*也能看到警告訊息 : Warning: The maximum likelihood estimate may not exist. */
```

所以我們一樣要用 exact 來檢定  
又因為 exact 不適用於加權過的資料，所以先做轉換  
```sas
data class12b1;
set class12b;
do i = 1 to count ; /*依照 count 的這個變數把資料展開*/
	output;
end;
run;

proc logistic data = class12b1 desc;
class race month / param = ref ref = first;
model prom = race month;
exact race month; /*exact 不能用加權的*/
run;
```
![image](https://hackmd.io/_uploads/BkVSlJ4u6.png)  
可以看出 rece 的 $p < 0.05$，因此這個變數是顯著的  

另外，也可以使用 cmh 來計算，結果應該會是一樣的  
```sas=
proc freq data = class12b;
weight count;
tables month * race * prom / cmh;
run;
/*用logistic跟cmh所得到的結果應該是一樣的--但此資料可能因為細格具有太多0，而exact絕對是比較準確的*/
/*另外，cmh 比較適用於大樣本的資料*/
```




# 多元模型  
主要處理把很多不同的變數做整合，做成一個回歸模型  
也就是一般常看到的 logistic regression 的延伸版  
大部分都跟原本的一樣  
唯一不同的地方就是要把其中一個類別設為 "參照組"  
通常在題目看到 **estblish baseline category model** 就是要建立多元模型  


## 程式碼實作--test5  
設定資料檔  
```sas=
%let dirclass =  C:\Users\ACER\Desktop\類別 sas 檔;
proc import datafile = "&dirclass\catadata112fq5.csv"
	DBMS = CSV out = q5;
run;
data q5a;
set q5;
if type = 'Healthy' then r1 = 0;
else r1 = 1;
if type = 'Scab' then r2 = 0;
else r2 = 1;
if type = 'Sprout' then r3 = 0;
else r3 = 1;
if type = 'Healthy' then r4 = 0;
else if type = 'Scab' then r4 = 1;
if type = 'Healthy' then r5 = 0;
else if type = 'Sprout' then r5 = 1;
if type = 'Scab' then r6 = 0;
else if type = 'Sprout' then r6 = 1;
run;
proc print data = q5a (obs = 10); run;
```

首先我們會遇到這種問題 :  
Establish a **baseline category model** to model type and use
”Healthy” based on class, density, hardness, weight. Write down the
fitted model.  

根據題目我們可以知道題目希望我們根據 type 來做模型，然後 type 裡面有分成 *Healthy(h), Sprout(sp), Scab(sc)*，我們要以 *Health(h)* 做為基底  
{type = sprout} : $logit(\cfrac{\pi_{sp}(x)}{\pi_h(x)}) = \hat{\beta_{01}} + \hat{\beta_{11}}x_{1} + \hat{\beta_{21}}x_{2} + \hat{\beta_{31}}x_{3} + \hat{\beta_{41}}x_{4}$  
{type = scab} : $logit(\cfrac{\pi_{sc}(x)}{\pi_h(x)}) = \hat{\beta_{02}} + \hat{\beta_{12}}x_{1} + \hat{\beta_{22}}x_{2} + \hat{\beta_{32}}x_{3} + \hat{\beta_{42}}x_{4}$  
接著執行多元羅基斯迴歸的模型來看該填入甚麼係數 :  
```sas=
proc logistic data = q5a desc;
class class/ param = ref ref = first;
/*想要參考的組別是 hrw ，所以設 ref = first*/
/*從報表上方可以看到是否 hrw 為 0*/
model type = class density  hardness  weight/ link = glogit ;
/*link = glogit : 多分類 Logistic 迴歸，也稱為廣義Logistic迴歸*/
/*可從報表上方看到是否是以 type = 'Healthy' 做為參考累別*/
run;
```
![image](https://hackmd.io/_uploads/ryvS8p7dp.png)  
可以先得出 model(從最大概度估計值分析觀察) :  
{type = sprout} : $logit(\cfrac{\pi_{sp}(x)}{\pi_h(x)}) = 19.2378 -0.4286* (class=srw) -14.8195*density -0.0184*hardness -0.00899*weight$  
{type = scab} : $logit(\cfrac{\pi_{sc}(x)}{\pi_h(x)}) = 31.9655 -0.3808* (class=srw) -21.3453*density -0.0107*hardness -0.2304*weight$  
並且可以發現 class 並不顯著($p > 0.05$)，其餘的都算是顯著  
接著再觀察係數，可以發現都是負的，代表會有負向影響  
然後可以發現 density 是非常重要的變數，因為它的點估計值非常小($<0.001$)  
結合上面的觀察可以知道  
當(density)上升時，勝算會下降 0.001 倍  
有(hardness)上升時，勝算會下降 0.98 倍  


---


接著改成 **continuation ratio logit model**  
也就是需要分成不同組別來看  

我們遇到的問題 :  
Establish a continuation ratio logit model to model type and use ”Healthy” based on class, density, hardness, weight. Write down the fitted model.

根據題目我們可以知道題目希望我們根據 type 來做模型，然後 type 裡面有分成 *Healthy(h), Sprout(sp), Scab(sc)*，我們要以 *Health(h)* 做為基底，而因為這是 contiunous 的，是希望我們做出更 *Healthy(h)*  
而從題目可以知道 H > Sc > Sp  
所以模型改成 :  
*Scab(sc)* 跟其他的比較(它會影響的因素是甚麼) :  
$logit(\cfrac{\pi_{sc}(x)}{\pi_{sp}(x)+\pi_h(x)}) = \hat{\beta_{01}} + \hat{\beta_{11}}x_{1} + \hat{\beta_{21}}x_{2} + \hat{\beta_{31}}x_{3} + \hat{\beta_{41}}x_{4}$  
*Sprout(sp)* 跟 *Healthy(h)* 的比較  
$logit(\cfrac{\pi_{sp}(x)}{\pi_h(x)}) = \hat{\beta_{01}} + \hat{\beta_{11}}x_{1} + \hat{\beta_{21}}x_{2} + \hat{\beta_{31}}x_{3} + \hat{\beta_{41}}x_{4}$  

因此，分別執行多元羅吉斯迴歸  
*Scab(sc)* 跟其他的比較 :  
```sas=
proc logistic data = q5a; 
/*之所以沒有用 desc 是因為我們想看的就是 sc*/
/*上述設定檔當中 r2 = 0 就是 sc*/
class class/ param = ref ref = first;
model r2 = class density  hardness weight;
run;
```
![image](https://hackmd.io/_uploads/HkbzlAmda.png)  
解釋方法跟上面一模一樣  
可以發現 density 跟 weight 的 $p < 0.0001$，代表都是非常顯著的影響變數  
係數都是負的所以都是負向影響  

*Sprout(sp)* 跟 *Healthy(h)* 的比較  
```sas=
proc logistic data = q5a desc; 
/*加上了 desc 是因為 sp 在 r5 的設定中 = 1*/
class class / param = ref ref = first;
model r5 = class density hardness weight;
run;
```
![image](https://hackmd.io/_uploads/HJb2lRm_6.png)  
解釋方法跟上面一模一樣  
可以發現 class 、 density 、 weight 都是顯著的影響變數， $p < 0.05$  


## 程式碼實作--class  
設定資料檔  
```sas=
/*ch6-p7*/
data gator ;
input length  choice $ @@;
datalines;
1.24 I 1.30 I 1.30 I 1.32 F 1.32 F 1.40 F 1.42 I 1.42 F
1.45 I 1.45 O 1.47 I 1.47 F 1.50 I 1.52 I 1.55 I 1.60 I
1.63 I 1.65 O 1.65 I 1.65 F 1.65 F 1.68 F 1.70 I 1.73 O
1.78 I 1.78 I 1.78 O 1.80 I 1.80 F 1.85 F 1.88 I 1.93 I
1.98 I 2.03 F 2.03 F 2.16 F 2.26 F 2.31 F 2.31 F 2.36 F
2.36 F 2.39 F 2.41 F 2.44 F 2.46 F 2.56 O 2.67 F 2.72 I
2.79 F 2.84 F 3.25 O 3.28 O 3.33 F 3.56 F 3.58 F 3.66 F
3.68 O 3.71 F 3.89 F
;
run;
```

接著一樣是跑羅吉斯迴歸  
```sas=
proc logistic data = gator desc; /*因為 F 是比較小的，所以設定 desc*/
model choice = length / link = glogit;
output out = out_gator p = pred predprobs= i; /*看output*/
run;
/*從回應概況可以看到 O 只有 8 而已，顯示資料不平衡*/
```
![image](https://hackmd.io/_uploads/S1UNQkVdT.png)  
但從上方的回應概況可以發現 O 只有 8 個而已，代表資料不平衡  
這樣子會讓最後分類出來的結果不理想  


---

設定資料檔  
```sas=
/*ch6-p20*/
data class12c;
input party ideo count @@;
cards;
0 1 80 0 2 81 0 3 171 0 4 41 0 5 55
1 1 30 1 2 46 1 3 148 1 4 84 1 5 99
run;
```

接著一樣是跑羅吉斯迴歸  
```sas=
proc logistic data = class12c;
weight count; /*因為有整合過所以要加權*/
class party / param = ref ref = first; 
/*party 有 2 類所以設 class，而且因為想看民主黨所以ref = first*/
model ideo = party / link = clogit; 
/*cumulative logit (default)*/ /*累積勝算*/
run;
/*看上方回應概況可以知道是以較低的排序值進行累積*/
/*看資料，我們的設定是 0 : 民主黨 ， 1 : 共和黨*/
/*而且從下面的MLE分析可以看到 party 是負的*/
/*
logit(P[Y \le j|x]) = alpha_j + beta x
x = 0 (demo)； x = 1 (repub)
比較民主的勝算
共和黨比較民主的勝算是民主黨的 0.377 倍
*/
```
![image](https://hackmd.io/_uploads/HyJnBk4uT.png)  
我們可以先從上方的回應概況看到是以較低的排序值進行累積  
並且從上方的設定可以知道我們設定 0 : 民主黨 ； 1 : 共和黨  
接著從估計值可以看到 party 的係數是負的所以是有負向相關  
所以對勝算比做的推論就會是 :  
共和黨比較民主的勝算是民主黨的 0.377 倍  
然後還可以得到模型是 :  
$logit(P[Y \le j|x]) = \alpha_j + \beta (x)$  
$x = 0 (demo)； x = 1 (repub)$  




# 零膨脹模型  
抽出來的資料當中，分為 "有"或"無"，然後在 "有" 的當中又有可能包含 "0"  
這樣會導致取出來的資料當中包含了非常多的 "0"  

## 程式碼實作--class  
設定資料檔  
```sas=
data class14d1;
input c r count;
/* P[R=3]/P[R=1 or R=2] */
if r = 3 then r1 =0;
else r1 =1;
/*r1 : nonlive 比上其他*/
/* Give R<3, P[R=2]/P[R=1] */
if r = 2 then r2 =0;
else if r=1 then r2 =1;
cards;
0     3  15
0     2   1
0     1 281
62.5  3  17
62.5  2   0
62.5  1 225
125   3  22
125   2   7
125   1 283
250   3  38
250   2  59
250   1 202
500   3 144
500   2 132
500   1   9
run;
```

接著是要做羅吉斯迴歸，先做 $R = 3$ 比上 $R = 1, 2$  
```sas=
/*R = 3 vs R = 1, 2*/
proc logistic data = class14d1;    /*因為有興趣的是 0 所以不用設定 desc*/
weight count;    /*因為是列聯表所以需要做整合*/
model r1 =c;
run;
```
![image](https://hackmd.io/_uploads/SJB73zNuT.png)  
可以看到 c 大於 0 且 p < 0.001 是顯著的  
代表劑量增加的話，死亡的機率也會增加  

接著繼續做另一個羅吉斯迴歸比較 $R=2$ 跟 $R=1$  
```sas=
/*R = 2 vs R = 1*/
proc logistic data = class14d1;    
weight count;    
model r2 =c;
run;
```
![image](https://hackmd.io/_uploads/H1476GEO6.png)  
可以看到 c 大於 0 且 p < 0.001 是顯著的  
代表劑量增加的話，死亡的機率也會增加  

---

設定資料檔  
```sas=
data class14d2;
infile "C:\Users\ACER\Desktop\類別 sas 檔\crab.dat";
input color spine width sat weight;
y = (sat>0);
if width<23.25 then gp=1;
else if width<24.25 then gp=2;
else if width<25.25 then gp=3;
else if width<26.25 then gp=4;
else if width<27.25 then gp=5;
else if width<28.25 then gp=6;
else if width<29.25 then gp=7;
else gp=8;
run;
```

接著可以做 `proc genmod`  
```sas=
proc genmod data = class14d2;
model sat = width / link = log dist = poi;
run;
/*用 df 查卡方的表會發現超出臨界值，為過度離散，所以要調整尺度*/
```
![image](https://hackmd.io/_uploads/SJTw9X4Oa.png)  
可以看到目前縮放是 $1$  

![image](https://hackmd.io/_uploads/rkgv9m4dT.png)  
可以看到卡方的自由度為 171，對照查表找出來 0.95 的值大概為 202  
但我們算出來的值是 567，大於臨界值，reject $H_0$，所以是過度離散的  

因此我們下方就做調整，在 model 的後面加上 scale  
```sas=
proc genmod data = class14d2;
model sat = width / link = log dist = pois scale = p;
/*scale 是調整尺度*/
run;
```
接著就能夠從圖中看到縮放為 1.78  
![image](https://hackmd.io/_uploads/ByZsoQEd6.png)  
並且也能知道經過縮放之後 **估計值** 不會有改變  

---

設定資料檔  
```sas=
proc import datafile = "C:\Users\ACER\Desktop\類別 sas 檔\table3_4.xlsx" 
	out=class14d3 replace;
run;
```

因為現在要做的是跟 1975 年做比較，所以設了一個 x 變數  
等於是把 1975 年當 ref  
```sas=
data class14d3a;
set class14d3;
x = year-1975;
offset = log(train_km);  /*Index*/ /*設定補償項*/
run;
```



```sas=
proc genmod data = class14d3a;
model Train_road_collision = x / link = log dist = pois offset = offset;
/*
可以從評估配適度的準則那裏看過度離散的情況
自由度為 27 卡方值為 40 左右，看到pearson 卡方是過度離散
offset = offset 就是把變數設成 offset
*/
run;
```
![image](https://hackmd.io/_uploads/Hy0v6mVua.png)  
從圖上可以看到 $\beta$ 的係數是 $-0.0329$  
有負向相關  
也就是說跟 1975 年比較有越來越少  

接著來看過度離散的情況  
![image](https://hackmd.io/_uploads/B1936m4da.png)  
從卡方的自由度為 27 去找值，值為 40.11  
所以可以看到 pearson 卡方有過度離散  

調整過度離散  
```sas=
proc genmod data = class14d3a;
model Train_road_collision = x / link = log dist = poi offset = offset scale = p; /*因為過度離散，所以調整scale*/
run;
```
![image](https://hackmd.io/_uploads/SyOU0mEua.png)  
![image](https://hackmd.io/_uploads/SyaHR7E_p.png)  
可以看到已經有調整完成了  

---

**零膨脹模型**  

設定資料檔  
```sas=
libname class14  "C:\Users\ACER\Desktop\類別 sas 檔";

data class14d4;
set class14.fish;
run;
```

```sas=
proc genmod data = class14d4;
class livebait camper;
/*persons child當成是連續地所以不用設定class*/
model count = livebait camper persons child / dist = zip link = log;
zeromodel / link = logit;
run;
```
最後的輸出會分為 2 個部分  
最大概度參數的部分以及 zero model 的部分  
![image](https://hackmd.io/_uploads/B1RKvV4Oa.png)  
可以看出截距項的部分是顯著的，也就是說，0 的部分是需要被納入估計的  
而因為上述 zeromodel 沒有放入參數的關係，不確定哪些會影響 0 的機率  

因此下方我們加入參數  
```sas=
proc genmod data = class14d4;
class livebait camper;
model count = livebait camper persons child / dist = zip link = log;
zeromodel livebait camper persons child / link = logit;
/*假設出現 0 的機率跟這些變數也有關係*/
run;
```
可以看到下方的 zeromodel 的報表  
![image](https://hackmd.io/_uploads/B1-q_N4da.png)  
所以可以得到以下解釋 :  
**livebait** 的 p > 0.05 因此不會影響  
**camper** 的部分因為 p < 0.05，然後估計值的係數是正的，所以可以知道有 camper 的話釣魚勝算會比較高  
**person** 的部分因為 p < 0.05，然後估計值的係數是負的，所以可以知道 person 越多釣魚的勝算會越低  
**child** 的部分因為 p < 0.05，然後估計值的係數是正的，所以可以知道有 camper 的話釣魚勝算會比較高  

然後因為 livebait 不會影響所以可以刪除掉，再重新建立一個模型  
```sas=
proc genmod data = class14d4;
class livebait camper;
model count = livebait camper persons child / dist = zip link = log;
zeromodel  camper persons child / link = logit;
run;
/*
P[ Z = 0] = pi
mu = number of fish
x1 = bait, x2 = camper, x3 = pearson, x4 = child
logit(\hat{pi}(x2, x3)) = 0.7213 + 0.8389 x2 - 0.9456 x3 + 1.9726 x4
log(\hat{\mu}(x1, x2, x3, x4)) = -0.0301 - 1.7313 x1 - 0.5975 x2 + 0.8344 x3 - 1.1770 x4
*/
```
![image](https://hackmd.io/_uploads/B1c3YVVdT.png)  
可以先建立出零膨脹的模型 :  
$\begin{cases} P[Z = 0] = \pi \\ \mu = \text{number of fish} \\ x_1 = \text{bait},\ x_2 = \text{camper},\ x_3 = \text{pearson},\ x_4 = \text{child} \end{cases}$  
$logit(\hat{pi}(x_2, x_3)) = 0.7213 + 0.8389 x_2 - 0.9456 x_3 + 1.9726 x_4$  

![image](https://hackmd.io/_uploads/r1IpK4EOa.png)  
然後在建立出可以釣多少魚的模型  
$log(\hat{\mu}(x_1, x_2, x_3, x_4)) = -0.0301 - 1.7313 x_1 - 0.5975 x_2 + 0.8344 x_3 - 1.1770 x_4$  



# 列聯表的對數線性模型  


## 程式碼實作--class  
設定資料檔  
```sas=
data class15a;
do a = 1, 0;
	do c = 1, 0;
		do m = 1, 0;
			input count @@;
			output;
		end;
	end;
end;
cards;
911 538 44 456 3 43 2 279
run;
```

建立獨立模型  
```sas=
/*independent model*/
proc genmod data = class15a;
class a c m / ref = first;
/*如果沒有 ref = first : 會用最後當參考*/
model count = a c m / dist = poi link = log;
run;
```
會得到這個結果  
![image](https://hackmd.io/_uploads/SkgA7y8V_6.png)  

建立邊際獨立模型  
會分為 3 種 :  
AC M、AM C、A CM
```sas=
/*marginal independent*/
/* AC M*/
proc genmod data = class15a;
class a c m / ref = first;
model count = a c m a*c / dist = poi link = log;
run;
/*AM C*/
proc genmod data = class15a;
class a c m / ref = first;
model count = a c m a*m / dist = poi link = log;
run;
/*從偏差來看， AC 的模型會比 AM 好(偏差的值越小越好--因為比較接近飽和模型)*/
/*A CM*/
proc genmod data = class15a;
class a c m / ref = first;
model count = a c m c*m / dist = poi link = log;
run;
/* CM > AC > AM (模型的好壞')*/
```
![image](https://hackmd.io/_uploads/Sy5B18NOT.png)
![image](https://hackmd.io/_uploads/H1LOkI4ua.png)
![image](https://hackmd.io/_uploads/BJUZ7IVu6.png)  
| AC     | AM  |
|:------ | --- |
| **CM** |     |

評估模型好壞的方法就是看它偏差的值，越小越好  
所以 : $(CM = 534) < (AC = 843) < (AM = 939)$  
得到 CM 的模型是最好的  

接著我們繼續做條件獨立的模型  
```sas=
/*conditional model*/
/*因為從上面可以知道 CM 最重要--所以從 CM 開始*/
/*CM AC*/
proc genmod data = class15a;
class a c m / ref = first;
model count = a c m c*m a*c/ dist = poi link = log;
run;
/*CM AM*/
/*預期這個的偏差會比上面還要大一點點*/
proc genmod data = class15a;
class a c m / ref = first;
model count = a c m c*m a*m/ dist = poi link = log;
run;
/*AC AM*/
/*預期這個的偏差會比上面還要更大一點點*/
proc genmod data = class15a;
class a c m / ref = first;
model count = a c m a*c a*m/ dist = poi link = log;
run;
```
然後用一樣的方法評估模型好壞  

最後，接著做同質(齊一)模型  
```sas=
/*Homogeneous model*/
/*AC AM CM*/
proc genmod data = class15a;
class a c m / ref = first;
model count = a c m a*c a*m c*m/ dist = poi link = log;
run;
```

最後的最後，做出飽和模型  
```sas=
/*ACM*/
proc genmod data = class15a;
class a c m / ref = first;
model count = a c m a*c a*m c*m a*c*m/ dist = poi link = log;
run;
```

---

設定資料檔  
```sas=
data ch7d2;
do g=0, 1;
	do l=0, 1;
		do s=0, 1;
			do i=0, 1;
				input count @@;
				output;
			end;
		end;
	end;
end;
cards;
7287 996 11587 759
3246 973 6134 757
10381 812 10969 380 
6123 1084 6693 513
run;
```

接著可以做 homogeneous model  
```sas=
/*homogeneous model*/
proc genmod data = ch7d2;
class g l s i / ref = first;
model count = g l s i g*l g*s g*i l*s l*i s*i / dist = poi link = log;
output out = out_hom p = p;
run;
```
可以看到  
![image](https://hackmd.io/_uploads/r1xfn8VOT.png)  
其中，df 之所以是 5 是因為(組合數 = 2^4 = 16) - (主要效果 = 4) - (交互作用 = 6) - (截距項 = 1) = 5
另外，也可以看到 `完整對數概度` = -88，而 AIC = 196  
從定義可以知道 $AIC = -2log(L) = -88 * (-2) = 196$  

接著做 3 維度的交互作用  
```sas=
/*GLS GI LI SI*/
/*
因為 I 沒有出現在 3 維度的交互作用中
所以 2 維度的要標上去
*/
proc genmod data = ch7d2;
class g l s i / ref = first;
model count = g l s i g*l g*s g*i l*s l*i s*i g*l*s / dist = poi link = log;
output out = out_gls p = p;
run;
/*
-2 log L = 80.4062*2
AIC = 184.9239
LR test G^2 = 88*2 - 80*2 = 16 > chi_1^2 = 3.84
代表 GLS 是顯著的(reject H0)
*/

/*GLI GS LS SI*/
proc genmod data = ch7d2;
class g l s i / ref = first;
model count = g l s i g*l g*s g*i l*s l*i s*i g*l*i / dist = poi link = log;
output out = out_gli p = p;
run;
/*
-2 log L = 86*2
AIC = 196.0288
LR test G^2 = 88*2 - 86*2 = 4
*/

/*GSI GL SL LI*/
proc genmod data = ch7d2;
class g l s i / ref = first;
model count = g l s i g*l g*s g*i l*s l*i s*i g*s*i / dist = poi link = log;
output out = out_gsi p = p;
run;
/*
-2 log L = 88*2
AIC = 200.3062
LR test G^2 = 88*2 - 88*2 = 0
*/

/*SLI GS GL GI*/
proc genmod data = ch7d2;
class g l s i / ref = first;
model count = g l s i g*l g*s g*i l*s l*i s*i s*l*i / dist = poi link = log;
output out = out_sli p = p;
run;
/*
-2 log L = 87*2
AIC = 198.0928
LR test G^2 = 88*2 - 87*2 = 2
*/
```

可以看他們各自的 AIC 判斷是否是顯著的  
以 GLS 舉例  
-2 log L = 80.4062*2  
AIC = 184.9239  
LR test G^2 = 88*2 - 80*2 = 16 > chi_1^2 = 3.84  
代表 GLS 是顯著的(reject H0)  
除此之外也可以從他們的報表中看到  
只有在最下面 g*l*s 全部都是 1 的時候有數值  
就看那一列當中的 p 值可以看到 < 0.05  
代表 GLS 這項變數是重要的  
![image](https://hackmd.io/_uploads/B1C7JvV_a.png)  
剩下的就以此類推  

接著，可以算他們是否有過度配適的問題  
在課本 ch8-p35 那裏  
計算 dissimilarity, $\hat{\Delta}$  
使用課本中的算法  
在這裡使用程式幫我們做運算  
```sas=
/*計算總樣本數 n = 68694*/
data out_hom1;
set out_hom;
diff = abs(count - p)/(2*68694);
run;
proc means data = out_hom1 sum;
var diff;
run;

data out_gls1;
set out_gls;
diff = abs(count - p)/(2*68694);
run;
proc means data = out_gls1 sum;
var diff;
run;
```

![image](https://hackmd.io/_uploads/H1yGxvVd6.png)
![image](https://hackmd.io/_uploads/H1KMeDVd6.png)  
最後可以得出結論，使用 homogeneous 的模型就好，因為 GLS 的模型已經過度配適了  






















