有時候類別變數不只分為有/無、是/否兩種類別，更多的時候是包含好幾種類別。例如地理可以分為北、中、南；宗教可分為回教、基督教、佛教；交通工具可分為汽車、機車、大眾運輸等等。這時候一般的羅吉斯回歸就無法派上用場，取而代之的是多項羅吉斯回歸((Multicategory Logit Model)。
多項羅吉斯回歸的概念其實與虛變數一樣，都是透過設定一個參考基準(baseline)來和其他類別做比較。當應變數的類別只分為兩種的時候，也就是最單純的羅吉斯回歸模型；應變數的類別在2個以上時，就可以利用設定參考類別的方式來建構模型。這裡要介紹的就是Baseline-Category多項羅吉斯迴歸模型，其數學式如下：
\[\log(\frac{\pi_{i}}{\pi_{j}})=\alpha x_{i} + \beta_{i}\]
其中\(\pi_{j}\)是參考類別，當\(\pi_{i}\)分別等於第1、第2類別時，與\(\pi_{j}\)相互比較可以得到下列logit關係式：
\[\log(\frac{\pi_{1}}{\pi_{2}})=\log\left(\dfrac{\frac{\pi_{1}}{\pi_{j}}}{\frac{\pi_{2}}{\pi_{j}}} \right)=\log(\frac{\pi_{1}}{\pi_{j}})-\log(\frac{\pi_{2}}{\pi_{j}})\]
哪一種鳥最不怕人? mlogit()
不同物種在不同情境下對人類活動有不同的敏感度，有些物種很容易親近人類，另一些物種則高度敏感。這類研究有助於環境影響評估，辨別出容易受人類活動干擾的物種。為了說明多項羅吉斯回歸的應用，我們引用A Bayesian optimal escape model reveals bird species differ in their capacity to habituate to humans。這篇論文比較七種水鳥近距離接觸人類的適應能力與容忍度。
所謂「驚飛距離」(Flight Initiation Distance, FID)是指鳥類對於接近中的威脅開始振翅起飛的距離，也就是人類可接近鳥類的最短距離。驚飛距離通常代表的鳥類對人類活動所建立的緩衝區，超過這個距離已經會讓鳥類感到危險而想逃跑。
這篇論文用複雜的模型量化黑天鵝(Black Swan)、家鴨(Domestic Duck)、黑水雞(Dusky Moorhen)、白冠雞(Eurasian Coot)、太平洋黑鴨(Pacfic Black Duck)、紫秧雞(Purple Swamphen)、澳洲紅嘴鷗(Silver Gull)對於人類活動的敏感度。我們從DRYAD取得原始資料Sutton et al Waterbird Data，用同一筆資料建立這七種水鳥對驚飛距離的羅吉斯回歸模型。首先先載入資料：
> bird<-read.csv("c:/Users/USER/Downloads/Sutton et al Waterbird Data.csv", header=T, sep=",")
> str(bird)
'data.frame': 848 obs. of 12 variables:
$ Site.Name : chr "ALBERT Park" "ALBERT Park" "ALBERT Park" "ALBERT Park" ...
$ Total.Number.of.Walkers: num 79 79 79 79 79 79 79 79 79 35 ...
$ Walker.Rate : num 30.6 30.6 30.6 30.6 30.6 ...
$ Species : chr "Black Swan" "Black Swan" "Black Swan" "Black Swan" ...
$ Flock.Size : int 1 3 1 1 1 6 2 2 2 1 ...
$ SD : num 42 29 57 57 34 63 40 40 42 65 ...
$ AD : num 17 23 40 48 34 37 40 40 42 51 ...
$ Response.Distance : num 17 21 22 32 34 37 40 40 42 37 ...
$ Towards..Away : chr "Towards" "Towards" "Towards" "Towards" ...
$ Seperation.Distance : num 0.5 0 3 0.5 6 1 1 6 1.5 12 ...
$ On.Shore.In.Water : chr "On Margin" "On Shore" "In Water" "In Water" ...
$ Distance.from.Edge : num 0.2 32 30 12 4 0.2 14 3 2 35 ...
> attach(bird)
> Species<-factor(Species, levels=c("Black Swan", "Domestic Duck",
+ "Dusky Moorhen", "Eurasian Coot", "Pacfic Black Duck", "Purple Swamphen",
+ "Silver Gull"), labels=c("黑天鵝", "家鴨", "黑水雞", "白冠雞", "太平洋黑鴨",
+ "紫秧雞", "澳洲紅嘴鷗")) #將字串轉換為類別變數
> table(Species)
Species
黑天鵝 家鴨 黑水雞 白冠雞 太平洋黑鴨 紫秧雞 澳洲紅嘴鷗
116 43 94 207 160 155 68
多項羅吉斯回歸
R有mlogit與nnet兩個延伸套件可以用來建構多項式羅吉斯回歸模型，這裡先嘗試用mlogit套件中的mlogit()指令：
> library(mlogit)
> waterbird<-mlogit.data(bird, choice="Species", shape="wide") #以mlogit.data()設定資料
> mlogit_model<-mlogit(formula=Species~0|Seperation.Distance, data=waterbird, reflevel="Black Swan") #設定羅吉斯迴歸以Black Swan為參考類別
> summary(mlogit_model)
Call:
mlogit(formula = Species ~ 0 | Seperation.Distance, data = waterbird, reflevel = "Black Swan", method = "nr")
Frequencies of alternatives:choice
Black Swan Domestic Duck Dusky Moorhen Eurasian Coot
0.137604 0.051008 0.111507 0.245552
Pacfic Black Duck Purple Swamphen Silver Gull
0.189798 0.183867 0.080664
nr method
7 iterations, 0h:0m:0s
g'(-H)^-1g = 2.94E-06
successive function values within tolerance limits
Coefficients :
Estimate Std. Error z-value
(Intercept):Domestic Duck -0.0565140 0.2442017 -0.2314
(Intercept):Dusky Moorhen 0.2334145 0.1783234 1.3089
(Intercept):Eurasian Coot 1.0169717 0.1503037 6.7661
(Intercept):Pacfic Black Duck 0.7971258 0.1584582 5.0305
(Intercept):Purple Swamphen 0.5994935 0.1551460 3.8641
(Intercept):Silver Gull 0.2305797 0.2044608 1.1277
Seperation.Distance:Domestic Duck -0.0741949 0.0209898 -3.5348
Seperation.Distance:Dusky Moorhen -0.0215533 0.0063332 -3.4032
Seperation.Distance:Eurasian Coot -0.0211345 0.0047293 -4.4689
Seperation.Distance:Pacfic Black Duck -0.0239006 0.0054981 -4.3471
Seperation.Distance:Purple Swamphen -0.0129986 0.0041143 -3.1593
Seperation.Distance:Silver Gull -0.0516382 0.0128779 -4.0098
Pr(>|z|)
(Intercept):Domestic Duck 0.8169859
(Intercept):Dusky Moorhen 0.1905548
(Intercept):Eurasian Coot 1.323e-11 ***
(Intercept):Pacfic Black Duck 4.892e-07 ***
(Intercept):Purple Swamphen 0.0001115 ***
(Intercept):Silver Gull 0.2594276
Seperation.Distance:Domestic Duck 0.0004081 ***
Seperation.Distance:Dusky Moorhen 0.0006660 ***
Seperation.Distance:Eurasian Coot 7.862e-06 ***
Seperation.Distance:Pacfic Black Duck 1.380e-05 ***
Seperation.Distance:Purple Swamphen 0.0015813 **
Seperation.Distance:Silver Gull 6.076e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Log-Likelihood: -1520.9
McFadden R^2: 0.021597
Likelihood ratio test : chisq = 67.144 (p.value = 1.573e-12)
mlogit()輸出了以黑天鵝為參考點的多項羅吉斯回歸報表。從分析結果中可以得知其他六種水鳥的羅吉斯回歸模型：
\(\text{家鴨： } \log(\frac{家鴨}{黑天鵝})=-0.074_{FID}-0.057\)
\(\text{黑水雞：}\log(\frac{黑水雞}{黑天鵝})=-0.022_{FID}+0.233\)
\(\text{白冠雞：}\log(\frac{白冠雞}{黑天鵝})=-0.021_{FID}+1.017\)
\(\text{太平洋黑鴨：}\log(\frac{太平洋黑鴨}{黑天鵝})=-0.024_{FID}+0.797\)
\(\text{紫秧雞：}\log(\frac{紫秧雞}{黑天鵝})=-0.013_{FID}+0.599\)
\(\text{澳洲紅嘴鷗：}\log(\frac{澳洲紅嘴鷗}{黑天鵝})=-0.052_{FID}+0.231\)
羅吉斯回歸並沒有像線性回歸一樣，有\(R^2\)判定係數，不過還是有學者專家試著以類\(R^2\)(pseudo R square)來替代。例如mlogit()的報表中輸出了McFadden \(R^2=0.021597\)模型解釋力不高。概似比檢定(likelihood ratio test)則可以看出模型中的自變數是否可以預測應變數。當模型中的自變數很多的時候，概似比檢定可以讓我們分辨哪一個自變數有解釋力。
模型解釋
與一般羅吉斯回歸一樣，多項羅吉斯回歸也是一種對數機率模型，模型係數是一種logit型態，經由指數轉算重新計算其勝算：
> coefficient<-exp(coef(mlogit_model))
> coefficient
(Intercept):Domestic Duck (Intercept):Dusky Moorhen
0.9450533 1.2629049
(Intercept):Eurasian Coot (Intercept):Pacfic Black Duck
2.7648094 2.2191535
(Intercept):Purple Swamphen (Intercept):Silver Gull
1.8211961 1.2593298
Seperation.Distance:Domestic Duck Seperation.Distance:Dusky Moorhen
0.9284907 0.9786773
Seperation.Distance:Eurasian Coot Seperation.Distance:Pacfic Black Duck
0.9790872 0.9763827
Seperation.Distance:Purple Swamphen Seperation.Distance:Silver Gull
0.9870855 0.9496724
在家鴨的模型中，上表顯示當驚飛距離增加一公尺，這隻水鳥是家鴨而不是黑天鵝的勝算是0.928，也就是每增加一公尺的驚飛距離，水鳥是家鴨的可能性就下降\((1-0.928) \times 100=7.2\%\)。驚飛距離越遠代表鳥類對人類越敏感、警覺性愈高，模型也顯示了和黑天鵝相比，家鴨更親近人類。
同樣地，我們可以依據模型計算各種水鳥的情況。在每增加一公尺驚飛距離的情況下，黑水雞的可能性下降2.2%、白冠雞下降2.1%、太平洋黑鴨下降2.4%、紫秧雞下降1.3%、澳洲紅嘴鷗下降5.0%。根據這份資料，我們大概可以推論這七種水鳥，家鴨對人類的警覺性最低，其次是澳洲紅嘴鷗，其他水鳥則是大同小異。
如果要比較家鴨與黑水雞，只需要把兩式相除，就能消掉參考類別黑天鵝，獲得家鴨與黑水雞的羅吉斯回歸方程式：
\(\log(\frac{\text{家鴨}}{\text{黑水雞}})=\log\left(\dfrac{\frac{\text{家鴨}}{\text{黑天鵝}}}{\frac{\text{黑水雞}}{\text{黑天鵝}}} \right)=\log(\frac{\text{家鴨}}{\text{黑天鵝}})-\log(\frac{\text{黑水雞}}{\text{黑天鵝}})\)
\(=(-0.074_{FID}-0.057)-(-0.022_{FID}+0.233)\)
\(=-0.052_{FID}-0.29\)
機率估計
Baseline-Category模型也可以用來估計機率，其估計方程式為：
\[\pi_{i}=\frac{e^{\alpha x_{i} + \beta_{i}}}{1 + \sum e^{\alpha x + \beta}}\]
因此我們可以列出七種水鳥的機率方程式：
\(\text{家鴨 }=\dfrac{e^{-0.074_{FID}-0.057}}{1+\sum e^{\alpha_{FID} + \beta}}\)
\(\text{黑水雞 }=\dfrac{e^{-0.022_{FID}+0.233}}{1+\sum e^{\alpha_{FID} + \beta}}\)
\(\text{白冠雞 }=\dfrac{e^{-0.021_{FID}+1.017}}{1+\sum e^{\alpha_{FID} + \beta}}\)
\(\text{太平洋黑鴨 }=\dfrac{e^{-0.024_{FID}+0.797}}{1+\sum e^{\alpha_{FID} + \beta}}\)
\(\text{紫秧雞 }=\dfrac{e^{-0.013_{FID}+0.599}}{1+\sum e^{\alpha_{FID} + \beta}}\)
\(\text{澳洲紅嘴鷗 }=\dfrac{e^{-0.052_{FID}+0.231}}{1+\sum e^{\alpha_{FID} + \beta}}\)
\(\text{黑天鵝 }=\dfrac{1}{1+\sum e^{\alpha_{FID} + \beta}}\)
根據上述機率方程式，當測得的驚飛距離為1的時候，這七種鳥類的機率依序分別為8%、11%、25%、20%、16%、11%、9%。所以，當人類靠近到1公尺的距離就飛走的鳥類，這種鳥最有可能是白冠雞，最不可能的是家鴨。
哪一種鳥最不怕人? multinom()
除了mlogit()之外，nnet套件中的multinom()不用轉換資料就能建立多項羅吉斯回歸模型，在應用上比mlogit()方便，但multinom()不會自動計算p值、\(R^2\)，要獲得這些統計量得需要額外計算。首先得用relevel()來設定參考類別：
> Species<-as.factor(bird$Species)
> Species_ref<-relevel(Species, ref="Black Swan")
多項羅吉斯回歸
設定黑天鵝為參考點候，用multinom()建立多項羅吉斯回歸模型：
> library(nnet)
> multinom_model<-multinom(Species_ref~Seperation.Distance, data=bird)
# weights: 21 (12 variable)
initial value 1640.402256
iter 10 value 1531.018174
iter 20 value 1520.916015
final value 1520.915560
converged
> summary(multinom_model)
Call:
multinom(formula = Species_ref ~ Seperation.Distance, data = bird)
Coefficients:
(Intercept) Seperation.Distance
Domestic Duck -0.0565549 -0.07418909
Dusky Moorhen 0.2334028 -0.02155325
Eurasian Coot 1.0169570 -0.02113435
Pacfic Black Duck 0.7971447 -0.02390176
Purple Swamphen 0.5994961 -0.01299853
Silver Gull 0.2305556 -0.05163492
Std. Errors:
Domestic Duck 0.2442000 0.020988477
Dusky Moorhen 0.1783236 0.006333205
Eurasian Coot 0.1503037 0.004729241
Pacfic Black Duck 0.1584582 0.005498200
Purple Swamphen 0.1551457 0.004114317
Silver Gull 0.2044594 0.012877283
Residual Deviance: 3041.831
AIC: 3065.831
multinom()的輸出報表沒有顯著性檢定，需要額外計算p值：
> z_score<-summary(multinom_model)$coefficients/summary(multinom_model)$standard.errors #計算z分數
> p_value<-(1- pnorm(abs(z_score),0,1))*2 #計算雙尾檢定p值
> p_value
(Intercept) Seperation.Distance
Domestic Duck 8.168545e-01 4.081560e-04
Dusky Moorhen 1.905777e-01 6.659817e-04
Eurasian Coot 1.323808e-11 7.863531e-06
Pacfic Black Duck 4.888652e-07 1.378886e-05
Purple Swamphen 1.115065e-04 1.581269e-03
Silver Gull 2.594741e-01 6.077839e-05
模型檢定
multinom()的輸出報表也沒有pseudo \(R^2\)，可以下載DescTools套件利用PseudoR2()計算：
> library(DescTools)
> PseudoR2(multinom_model, which=c("CoxSnell", "Nagelkerke", "McFadden"))
CoxSnell Nagelkerke McFadden
0.07655901 0.07852385 0.02159670
pseudo \(R^2\)顯示模型解釋力不高，所以我們將鳥群大小(Flock.Size)、離水岸邊的距離(Distance.from.Edge)兩個變數加入模型，藉此提高模型解釋力。合理推測當群體越大、離岸邊越遠的時候，鳥類會感覺相對安全，驚飛距離會縮小。
> multinom_model2<-multinom(Species_ref~Seperation.Distance+Flock.Size+Distance.from.Edge, data=bird)
> PseudoR2(multinom_model2, which=c("CoxSnell", "Nagelkerke", "McFadden"))
CoxSnell Nagelkerke McFadden
0.22291315 0.22863409 0.06838501
將擬合模型與只有常數項(only intercept model)的模型相比較，可以知道加入自變數之後模型是否有顯著差異：
> intercept_model<-multinom(Species_ref~1, data=bird) #常數項模型
> anova(intercept_model, multinom_model) #兩種模型的變異數分析
Likelihood ratio tests of Multinomial Models
Response: Species_ref
Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
1 1 5052 3108.975
2 Seperation.Distance 5046 3041.831 1 vs 2 6 67.14358 1.572964e-12
lmtest套件的lrtest()指令可以進行概似比檢定：
> library(lmtest)
> lrtest(multinom_model2, "Distance.from.Edge") #檢定岸邊距離
Likelihood ratio test
Model 1: Species_ref ~ Seperation.Distance + Flock.Size + Distance.from.Edge
Model 2: Species_ref ~ Seperation.Distance + Flock.Size
#Df LogLik Df Chisq Pr(>Chisq)
1 24 -1448.2
2 18 -1464.0 -6 31.704 1.859e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’
> lrtest(multinom_model2, "Flock.Size") #檢定鳥群大小
Likelihood ratio test
Model 1: Species_ref ~ Seperation.Distance + Flock.Size + Distance.from.Edge
Model 2: Species_ref ~ Seperation.Distance + Distance.from.Edge
#Df LogLik Df Chisq Pr(>Chisq)
1 24 -1448.2
2 18 -1504.7 -6 113.12 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’
機率估計
以multinom_model2進行機率估計，當驚飛距離=1，鳥群大小以及岸邊距離都等於平均數的時候：
> predict(multinom_model2, type="probs", newdata=data.frame(Seperation.Distance=1, Flock.Size=mean(Flock.Size), Distance.from.Edge=mean(Distance.from.Edge)))
Black Swan Domestic Duck Dusky Moorhen Eurasian Coot
0.08015436 0.08637194 0.09338987 0.26894800
Pacfic Black Duck Purple Swamphen Silver Gull
0.22768426 0.15637771 0.08707386
當可以更靠近鳥類，驚飛距離越來越近=0.5的時候：
> predict(multinom_model2, type="probs", newdata=data.frame(Seperation.Distance=0.5, Flock.Size=mean(Flock.Size), Distance.from.Edge=mean(Distance.from.Edge)))
Black Swan Domestic Duck Dusky Moorhen Eurasian Coot
0.07916200 0.08866341 0.09307404 0.26839112
Pacfic Black Duck Purple Swamphen Silver Gull
0.22754575 0.15518294 0.08798074
從機率估計可以看出來，家鴨與澳洲紅嘴鷗是比較不害怕人類的鳥類。multinom_model2模型預測的前10種水鳥，可以看見跟原始資料有明顯差異：
> head(predict(multinom_model2),10)
[1] Purple Swamphen Eurasian Coot Pacfic Black Duck Pacfic Black Duck
[5] Purple Swamphen Eurasian Coot Pacfic Black Duck Purple Swamphen
[9] Purple Swamphen Pacfic Black Duck
7 Levels: Black Swan Domestic Duck Dusky Moorhen ... Silver Gull