# 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
代表觀測值的相異程度不平均，在大樣本的卡方時會有問題
需要調整 standard error
但我們還是可以從 HL test 得知，他的 p 值 < 0.05 $\Rightarrow$ reject $H_0$
代表應該有過度離散的情況
另外，在解釋的部分，先從這裡看他的正負知道他們的關係是正向還是負向影響
可以先做出模型(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}$
然後也可以發現係數都是負的，代表都為負向影響
所以做解釋的時候都要以負向影響來解釋
所以當有改變(change = 1)時，勝算會下降 0.588 倍
有風(wind = 1)時，勝算會下降 0.509 倍
再來，我們要來看配對
從 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
可以看到 敏感度 跟 特異性 幾乎相等
*/
```
我們的切點是希望 `敏感度` 跟 `特異性` 幾乎相等
補充 :
:::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 值來看是否顯著*/
```
可以看到 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;
```
可以看出 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;
```
可以先得出 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;
```
解釋方法跟上面一模一樣
可以發現 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;
```
解釋方法跟上面一模一樣
可以發現 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 而已，顯示資料不平衡*/
```
但從上方的回應概況可以發現 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 倍
*/
```
我們可以先從上方的回應概況看到是以較低的排序值進行累積
並且從上方的設定可以知道我們設定 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;
```
可以看到 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;
```
可以看到 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 查卡方的表會發現超出臨界值，為過度離散，所以要調整尺度*/
```
可以看到目前縮放是 $1$
可以看到卡方的自由度為 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
並且也能知道經過縮放之後 **估計值** 不會有改變
---
設定資料檔
```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;
```
從圖上可以看到 $\beta$ 的係數是 $-0.0329$
有負向相關
也就是說跟 1975 年比較有越來越少
接著來看過度離散的情況
從卡方的自由度為 27 去找值，值為 40.11
所以可以看到 pearson 卡方有過度離散
調整過度離散
```sas=
proc genmod data = class14d3a;
model Train_road_collision = x / link = log dist = poi offset = offset scale = p; /*因為過度離散，所以調整scale*/
run;
```
可以看到已經有調整完成了
---
**零膨脹模型**
設定資料檔
```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 的部分
可以看出截距項的部分是顯著的，也就是說，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 的報表
所以可以得到以下解釋 :
**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
*/
```
可以先建立出零膨脹的模型 :
$\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$
然後在建立出可以釣多少魚的模型
$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;
```
會得到這個結果
建立邊際獨立模型
會分為 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 (模型的好壞')*/
```
| 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;
```
可以看到
其中，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 這項變數是重要的
剩下的就以此類推
接著，可以算他們是否有過度配適的問題
在課本 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;
```
最後可以得出結論，使用 homogeneous 的模型就好，因為 GLS 的模型已經過度配適了
Sign in with Wallet
Connect another wallet