Outline

Logistic regression
Multicatory logit model
Zero model
Loglinear model for contingency table
Models for matched pair

之前的重點小提醒 :

RR 跟 OR 的解釋技巧 :

假設
$R R = \frac{π_{1} \Rightarrow P}{π_{2} \Rightarrow A} = 1.82$ : P 會 …<題目>… 的機率是 A 的 1.82 倍
(
$∵ π_{1} = R R \times π_{2}$ )
假設
$O R = θ = \frac{Ω_{1} \Rightarrow P}{Ω_{2} \Rightarrow A} = 1.83$ : 使用 P 得到 …<題目>… 的勝算是使用 A 的 1.83 倍
(
$∵ Ω_{1} = θ \times Ω_{2}$ )

羅吉斯迴歸

主要是 2 元會使用到
首先會先選模，有以下這些方式(利用 SAS 自動選模) :

backward
stepwise
forward
這些方法就是把一些無關的特徵或雜訊值刪除掉

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

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

程式碼實作–test4

設定資料檔









%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) 取出可解釋的變數









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 來決定是否保留該筆資料，
因為有些資料的殘差太大(跟其他的差距過大)會讓模型配飾不準確







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
這些都是殘差過大的，之後需要被刪除





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 的模型





















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 Not Showing Possible Reasons

The image was uploaded to a note which you don't have access to
The note which the image was originally uploaded to has been deleted

Learn More →

代表觀測值的相異程度不平均，在大樣本的卡方時會有問題
需要調整 standard error
但我們還是可以從 HL test 得知，他的 p 值 < 0.05

\Rightarrow

reject

H_{0}

Image Not Showing Possible Reasons

The image was uploaded to a note which you don't have access to
The note which the image was originally uploaded to has been deleted

Learn More →

代表應該有過度離散的情況

另外，在解釋的部分，先從這裡看他的正負知道他們的關係是正向還是負向影響

Image Not Showing Possible Reasons

The image was uploaded to a note which you don't have access to
The note which the image was originally uploaded to has been deleted

Learn More →

可以先做出模型(logistic regression model) :

\begin{aligned} l o g i t (\hat{π (x)}) & = \hat{β_{0}} + \hat{β_{1}} x_{1} + \hat{β_{2}} x_{2} + . . . \\ = 7.0858 - 0.0297 * w e e k - 0.5309 * (c h a n g e = 1) - 0.6753 * (w i n d = 1) - 0.1312 * (d i s t a n c e = 1) \end{aligned}

然後也可以發現係數都是負的，代表都為負向影響
所以做解釋的時候都要以負向影響來解釋

Image Not Showing Possible Reasons

The image was uploaded to a note which you don't have access to
The note which the image was originally uploaded to has been deleted

Learn More →

所以當有改變(change = 1)時，勝算會下降 0.588 倍
有風(wind = 1)時，勝算會下降 0.509 倍

再來，我們要來看配對

Image Not Showing Possible Reasons

The image was uploaded to a note which you don't have access to
The note which the image was originally uploaded to has been deleted

Learn More →

從 c 值可以看到配適效果還不錯(接近 1，遠離 0.5)

然後我們想要知道未來我們要做預測該如何設定切點
因此使用 ctable 來看










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 Not Showing Possible Reasons

The image was uploaded to a note which you don't have access to
The note which the image was originally uploaded to has been deleted

Learn More →

Image Not Showing Possible Reasons

The image was uploaded to a note which you don't have access to
The note which the image was originally uploaded to has been deleted

Learn More →

我們的切點是希望 敏感度 跟 特異性 幾乎相等

補充 :





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

程式碼實作–class

假設現在資料檔有其中一個變數沒辦法估計的話，就要改成使用 exact 來看它的

p

值

設定資料檔










/*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;

接著做出羅吉斯迴歸
















/*從結果可以發現我們沒辦法估計 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 個變數是顯著的

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

設定資料檔

















/*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;

接著做出羅吉斯迴歸







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 不適用於加權過的資料，所以先做轉換

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 來計算，結果應該會是一樣的






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

多元模型

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

程式碼實作–test5

設定資料檔




















%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} :

l o g i t (\frac{π_{s p} (x)}{π_{h} (x)}) = \hat{β_{01}} + \hat{β_{11}} x_{1} + \hat{β_{21}} x_{2} + \hat{β_{31}} x_{3} + \hat{β_{41}} x_{4}

{type = scab} :

l o g i t (\frac{π_{s c} (x)}{π_{h} (x)}) = \hat{β_{02}} + \hat{β_{12}} x_{1} + \hat{β_{22}} x_{2} + \hat{β_{32}} x_{3} + \hat{β_{42}} x_{4}

接著執行多元羅基斯迴歸的模型來看該填入甚麼係數 :








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} :

l o g i t (\frac{π_{s p} (x)}{π_{h} (x)}) = 19.2378 - 0.4286 * (c l a s s = s r w) - 14.8195 * d e n s i t y - 0.0184 * h a r d n e s s - 0.00899 * w e i g h t

{type = scab} :

l o g i t (\frac{π_{s c} (x)}{π_{h} (x)}) = 31.9655 - 0.3808 * (c l a s s = s r w) - 21.3453 * d e n s i t y - 0.0107 * h a r d n e s s - 0.2304 * w e i g h t

並且可以發現 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) 跟其他的比較(它會影響的因素是甚麼) :

l o g i t (\frac{π_{s c} (x)}{π_{s p} (x) + π_{h} (x)}) = \hat{β_{01}} + \hat{β_{11}} x_{1} + \hat{β_{21}} x_{2} + \hat{β_{31}} x_{3} + \hat{β_{41}} x_{4}

Sprout(sp) 跟 Healthy(h) 的比較

l o g i t (\frac{π_{s p} (x)}{π_{h} (x)}) = \hat{β_{01}} + \hat{β_{11}} x_{1} + \hat{β_{21}} x_{2} + \hat{β_{31}} x_{3} + \hat{β_{41}} x_{4}

因此，分別執行多元羅吉斯迴歸
Scab(sc) 跟其他的比較 :






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) 的比較





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

設定資料檔














/*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;

接著一樣是跑羅吉斯迴歸





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 個而已，代表資料不平衡
這樣子會讓最後分類出來的結果不理想

設定資料檔







/*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;

接著一樣是跑羅吉斯迴歸
















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 倍
然後還可以得到模型是 :

l o g i t (P [Y \leq j | x]) = α_{j} + β (x)

x = 0 (d e m o) ； x = 1 (r e p u b)

零膨脹模型

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

程式碼實作–class

設定資料檔


























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





/*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





/*R = 2 vs R = 1*/
proc logistic data = class14d1;    
weight count;    
model r2 =c;
run;

可以看到 c 大於 0 且 p < 0.001 是顯著的
代表劑量增加的話，死亡的機率也會增加

設定資料檔













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




proc genmod data = class14d2;
model sat = width / link = log dist = poi;
run;
/*用 df 查卡方的表會發現超出臨界值，為過度離散，所以要調整尺度*/

可以看到目前縮放是

1

可以看到卡方的自由度為 171，對照查表找出來 0.95 的值大概為 202
但我們算出來的值是 567，大於臨界值，reject

H_{0}

，所以是過度離散的

因此我們下方就做調整，在 model 的後面加上 scale




proc genmod data = class14d2;
model sat = width / link = log dist = pois scale = p;
/*scale 是調整尺度*/
run;

接著就能夠從圖中看到縮放為 1.78

並且也能知道經過縮放之後 估計值 不會有改變

設定資料檔



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

因為現在要做的是跟 1975 年做比較，所以設了一個 x 變數
等於是把 1975 年當 ref





data class14d3a;
set class14d3;
x = year-1975;
offset = log(train_km);  /*Index*/ /*設定補償項*/
run;








proc genmod data = class14d3a;
model Train_road_collision = x / link = log dist = pois offset = offset;
/*
可以從評估配適度的準則那裏看過度離散的情況
自由度為 27 卡方值為 40 左右，看到pearson 卡方是過度離散
offset = offset 就是把變數設成 offset
*/
run;

從圖上可以看到

β

的係數是

- 0.0329

有負向相關
也就是說跟 1975 年比較有越來越少

接著來看過度離散的情況

從卡方的自由度為 27 去找值，值為 40.11
所以可以看到 pearson 卡方有過度離散

調整過度離散



proc genmod data = class14d3a;
model Train_road_collision = x / link = log dist = poi offset = offset scale = p; /*因為過度離散，所以調整scale*/
run;

可以看到已經有調整完成了

零膨脹模型

設定資料檔





libname class14  "C:\Users\ACER\Desktop\類別 sas 檔";

data class14d4;
set class14.fish;
run;






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 的機率

因此下方我們加入參數






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 不會影響所以可以刪除掉，再重新建立一個模型












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] = π \\ μ = number of fish \\ x_{1} = bait, x_{2} = camper, x_{3} = pearson, x_{4} = child \end{cases}

l o g i t (\hat{p i} (x_{2}, x_{3})) = 0.7213 + 0.8389 x_{2} - 0.9456 x_{3} + 1.9726 x_{4}

然後在建立出可以釣多少魚的模型

l o g (\hat{μ} (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

設定資料檔












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;

建立獨立模型






/*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


















/*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

評估模型好壞的方法就是看它偏差的值，越小越好
所以 :

(C M = 534) < (A C = 843) < (A M = 939)

得到 CM 的模型是最好的

接著我們繼續做條件獨立的模型



















/*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;

然後用一樣的方法評估模型好壞

最後，接著做同質(齊一)模型






/*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;

最後的最後，做出飽和模型





/*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;

設定資料檔

















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






/*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
從定義可以知道

A I C = - 2 l o g (L) = - 88 * (- 2) = 196

接著做 3 維度的交互作用




















































/*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.40622
AIC = 184.9239
LR test G^2 = 882 - 802 = 16 > chi_1^2 = 3.84
代表 GLS 是顯著的(reject H0)
除此之外也可以從他們的報表中看到
只有在最下面 gl*s 全部都是 1 的時候有數值
就看那一列當中的 p 值可以看到 < 0.05
代表 GLS 這項變數是重要的

剩下的就以此類推

接著，可以算他們是否有過度配適的問題
在課本 ch8-p35 那裏
計算 dissimilarity,

\hat{Δ}

使用課本中的算法
在這裡使用程式幫我們做運算
















/*計算總樣本數 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 的模型已經過度配適了

Outline

羅吉斯迴歸

程式碼實作–test4

程式碼實作–class

多元模型

程式碼實作–test5

程式碼實作–class

零膨脹模型

程式碼實作–class

列聯表的對數線性模型

程式碼實作–class

Read more

高公局比賽問題

cryptography 圖文解析

常用的 latex

mathstat_midterm_review2