lecture
website icon
累積羅吉斯回歸
勝算比 2x2交叉表 2x2xK交叉表 泊松回歸 羅吉斯回歸 多項羅吉斯回歸 累積羅吉斯回歸 Bradley-Terry模型 配對模型
×
website icon 資料管理 統計分析 相關資源 巨人肩膀 語法索引 關於作者

簡介

本部分介紹累積羅吉斯回歸,使用到的指令包含:

Facebook Icon Twitter Icon LinkedIn Icon LINE Icon

累積羅吉斯迴歸模型(Cumulative Logit Model)又稱順序羅吉斯回歸模型(Ordered Logit Model),是一種累加機率模型,適用於應變數為順序尺度的情況。其模型如下:

\[log \frac{P(Y \leq j)}{P(Y > j )}=logit \left( P(Y \leq j) \right) =\alpha x + \beta_{j} \]

不過在R的polr()中,預設的累積羅吉斯回歸的參數是:

\[log \frac{P(Y \leq j)}{P(Y > j )}=logit \left( P(Y \leq j) \right) =-\alpha x + \beta_{j} \]

其中j代表應變數Y順序尺度的項目,由於是累加機率所以應變數最後一個項目累加機率恆為1。j-1決定累積模型的數目。如果應變數Y有3個順序尺度,可以計算出3-1=2個累積羅吉斯迴歸模型。

鳥類也有護花使者?

織布鳥(Ploceus cucullatus)是一種原產於非洲的鳥類,有鳥界建築師之稱。每當繁殖季節,公鳥會用樹枝、樹葉、雜草當作材料築巢,並將築好的巢送給母鳥示愛,就像是展示自己有錢有房一樣,母鳥如果接受公鳥則會負起鳥巢內部裝修的責任。

織布鳥是一夫多妻制(polygynous)。鳥巢通常呈現球狀且開口向下,用於避開包含蛇、猴子、老鷹等天敵,當公鳥透過同伴呼叫感知到可能的威脅的時候,會有拋家棄子逃離鳥巢的避難行為。現在問題來了,我們通常都有男生應該是護花使者的刻板印象,尤其是遭遇危險的時候,更是展現男子氣概與騎士精神的時候,遇到危難先烙跑對男生更是大扣分的行為。如果人類的在追求女性的行為是這樣,那麼鳥類求偶的時候是不是也是如此呢? 在Male risk-taking is related to number of mates in a polygynous bird的研究中,試圖回答的就是公鳥會不會像人類一樣,在危急時刻會為了伴侶願意承擔更多的風險?

為了解答這項問題,首先先從DRYAD下載原始資料order_of_return_data.csv。其中Number of Females代表的是公織布鳥巢中的伴侶數、Order of Return是公織布鳥受到驚嚇逃離後,再度飛回鳥巢的次序。越早飛回鳥巢,我們假設公織布鳥會為了保護伴侶,而願意冒更多的生命風險。

> weaverbird<-read.csv("c:/Users/USER/downloads/order_of_return_data.csv", header=T, sep=",")
> return<-c("")
> return[weaverbird$Order.of.Return<=3]<-1
> return[weaverbird$Order.of.Return>=4 & weaverbird$Order.of.Return<=6]<-2
> return[weaverbird$Order.of.Return>=7]<-3
> return<-factor(return, levels=c(1,2,3), labels=c("Fast", "Normal", "Slow"))
> weaverbird<-cbind(weaverbird, return)
> attach(weaverbird)
> table(return)
return
  Fast Normal   Slow
   102     80     36

資料準備完成後,以MASS套件中的polr()進行累積羅吉斯迴歸分析。應變數是飛回鳥巢的順序,自變數是伴侶數:

> library(MASS)
> polr_model<-polr(formula=as.ordered(return)~Number.of.Females, data=weaverbird, method=c("logistic"))
> summary(polr_model)

Re-fitting to get Hessian

Call:
polr(formula = as.ordered(return) ~ Number.of.Females, data = weaverbird, method = c("logistic"))

Coefficients:
                    Value Std. Error t value
Number.of.Females -0.5213     0.1765  -2.953

Intercepts:
            Value   Std. Error t value
Fast|Normal -0.9700  0.3148    -3.0811
Normal|Slow  0.8352  0.3160     2.6426

Residual Deviance: 436.0094
AIC: 442.0094

輸出報表分為模型、係數、截距以及偏差殘差與AIC等四大部分。其中最重要的就是係數與截距,可得出Fast與Normal兩個羅吉斯回歸模型:

\(\text{Fast: }logit (P(Y \leq 1) = -(-0.52_{Female})-0.97\)

\(\text{Normal: }logit (P(Y \leq 2) = -(-0.52_{Female})+0.84\)

polr()雖然計算出t值但沒有p值,無法得知顯著與否。一個簡單計算p值的方式是像z檢定一樣,把t值放到標準常態分配:

> coef<-coef(summary(polr_model)) #輸出polr_model的係數並存放在coef
> p_value<-pnorm(abs(coef[,"t value"]), lower.tail=FALSE)*2 #利用coef的t值,計算雙尾檢定p值
> p_value #輸出結果
Number.of.Females       Fast|Normal       Normal|Slow
      0.003144886       0.002062031       0.008226470
> p_value<-(1-pnorm(abs(coef[,"t value"])))*2 #直接用1減去累積機率,也能獲得相同結果
Number.of.Females       Fast|Normal       Normal|Slow
      0.003144886       0.002062031       0.008226470
> coef<-cbind(coef, "p value"=round(p_value, digits=3)) #合併報表
> coef
                       Value Std. Error   t value p value
Number.of.Females -0.5212964  0.1765186 -2.953209   0.003
Fast|Normal       -0.9700407  0.3148308 -3.081149   0.002
Normal|Slow        0.8351630  0.3160347  2.642631   0.008

模型檢定

羅吉斯回歸沒有\(R^2\),polr()報表也沒有pseudo \(R^2\),可以下載DescTools套件利用PseudoR2()來計算:

> PseudoR2(polr_model, which=c("CoxSnell", "Nagelkerke", "McFadden"))
  CoxSnell Nagelkerke   McFadden
0.04043662 0.04647125 0.02022070

將擬合出來的模型與常數項(only intercept model)模型相比較,可以獲得加入自變數之後,模型是否有顯著差異的資訊:

> intercept_model<-polr(formula=as.ordered(return)~1, data=weaverbird, method=c("logistic")) #常數項模型
> anova(intercept_model, polr_model) #兩種模型的變異數分析
Likelihood ratio tests of ordinal regression models

Response: as.ordered(return)
              Model Resid. df Resid. Dev   Test    Df LR stat.     Pr(Chi)
1                 1       216   445.0078
2 Number.of.Females       215   436.0094 1 vs 2     1 8.998367 0.002702209

模型解釋

要注意的是,polr()的累積羅吉斯模型x項係數是\(-\alpha\),所以模型中Number of Females的係數其實是\(-(-0.52)\)。依據模型,對於那些巢中有越多伴侶的公織布鳥來說,慢速飛回vs.以正常速度或快速飛回鳥巢的對數勝率(log odds)是0.52。對數勝率在解釋上畢竟比較難懂,所以大多會將再將係數用exp()還原成勝算比:

> exp(coef(polr_model))
Number.of.Females
        0.5937503

上述伴侶數勝算比是\(exp(-\alpha)\)的結果,也就是\(exp(-0.52)=\frac{1}{exp(\alpha)}=0.59\)。這其實是應變數第三個類別慢速飛回(Slow)的勝算比。通常我們會希望直接求得\(exp(\alpha)\),也就是第一個類別快速飛回(Fast)的勝算比,在模型解釋上會比較直觀,所以我們取倒數\(\frac{1}{0.59}=1.69\),或者直接計算\(exp(0.52)\):

> 1/exp(coef(polr_model))
Number.of.Females
          1.68421

上述結果顯示,對於一夫多妻的公織布鳥來說,鳥巢中每多出一隻母鳥,他快速(fast)飛回巢中的勝算,是正常(normal)或慢速(slow)飛回的1.68倍,也就是增加68%。我們可以用predict()來實際預測當伴侶數分別是1、2、3隻的時候,公織布鳥飛回的機率:

> forecast<-data.frame(Number.of.Females=c(1,2,3))
> predict(polr_model,forecast, type="probs")
       Fast    Normal      Slow
1 0.3896594 0.4055243 0.2048163
2 0.5181301 0.3492231 0.1326469
3 0.6442477 0.2725073 0.0832450

上表可以看到,當伴侶數增加的時候,公織布鳥快速飛回巢中的機率也增加,而正常與慢速飛回巢中的機率則是越來越低。快速飛回的勝算比1.68,則是這麼計算出來的:

\(\dfrac{^{0.518}/_{1-0.518}}{^{0.390}/_{1-0.939}}=1.68\)

同理,可以計算以正常速度及慢速飛回巢中的勝算比:

\(\dfrac{^{0.349}/_{1-0.349}}{^{0.406}/_{1-0.406}}=0.78\)

\(\dfrac{^{0.133}/_{1-0.133}}{^{0.205}/_{1-0.205}}=0.59\)

累積機率

既然是累積機率模型,我們可以根據polr()計算出的模型,分別計算快速、正常、慢速三種飛回順序的累積機率:

\(\text{Fast: }(P \leq 1 )=\dfrac{exp(0.52 \times 1 -0.97)}{1-exp(0.52 \times 1 -0.97)}=81.6\% \)

\(\text{Normal: }(P \leq 2 )=\dfrac{exp(0.52 \times 2 +0.84)}{1-exp(0.52 \times 2 +0.84)}=5.2\% \)

\(\text{Slow: }(P \leq 3 )=1-\dfrac{exp(0.52 \times 2 +0.84)}{1-exp(0.52 \times 2 +0.84)}=13.2\% \)

所以大部分的公織布鳥在逃離鳥巢後,其實返回鳥巢的速度是快的;累積羅吉斯回歸模型顯示,伴侶數越多公織布鳥飛回的順序越快,鳥類其實也有護花使者的行為。