有時候應變數是以「次數」為單位,例如一年內的族群數目、一個月的確診人數、一周內的案件次數、一天內的發生次數等等。泊松(Siméon Poisson)在1838年發表的一個離散機率分布,適合於描述一個給定時間內,隨機事件發生次數的機率分布,這項分布以他的名字命名為泊松分配。泊松分配在自然科學與社會科學都有廣泛的應用:
- 美國總統任命大法官的分布(King, 1987)
- 貨船被大浪襲擊的次數(McCullagh and Nelder, 1989)
- 精神病患治療後半年發生暴力事件的次數(Gardner, Mulvey and Shaw, 1995)
- 飛機遭遇鳥擊的次數(Vaira, 2012)
- 二次感染Covid-19的人數(Kremer, Torneri and Boesmenas et al., 2021)
泊松迴歸的應變數是可以計算的次數(count),如果自變數是連續變數,稱之為泊松回歸模型(Poisson regression model);如果自變數是類別變數,傳統上稱之為對數線性模型(Log-linear model)。所以也有人稱泊松回歸為泊松對數線性回歸,之所以稱為對數線性,其實只不過是把應變數取\(\log\)而已:
\[\log(\mu)=\alpha x + \beta \]
經由對數-指數轉換,可改寫為:
\[\mu=\exp(\alpha x + \beta) = e^{\alpha x} e^{\beta}\]
泊松迴歸與簡單回歸的差別僅有左邊的\(log\),所以模型解釋的方法和簡單回歸一樣,當\(\alpha\)係數為正數時,代表自變數對應變數有正面影響,反之亦然。\(x\)每增加一個單位,\(\mu\)會增加\(e^\alpha\)倍。至於在做預測的時候,則是要將原本的對數模型,轉換為指數型態。
泊松回歸:計數資料
NBA近年來的攻守步調加快、得分屢創新高,其中最明顯的改變就是三分球投射比例增加。多數球評都認為,造成這種情況與防守規則鬆綁有關,其中最為大家所討論的就是hand-checking。
Hand-checking也就是手部接觸防守規則,指的是防守方用手接觸對方身體,通常是放在對方的腰、背、手肘等處來幫助防守,身體接觸讓防守方更容易追蹤進攻球員的動作,對持球進攻的球員構成很大的防守威脅。然而這項利於防守的規則,在2004-2005賽季被禁止,只有在禁區內面對背框進攻的時候,才允許防守方手部接觸。這項規則被視為是造成現在大三分時代最關鍵的因素。
為了驗證這項說法,我們從Basketball Reference取得開始有三分球以來1979-2021年全聯盟的平均三分球出手數作為應變數。因為Hand-checking不好量化,但我們假設手部接觸越頻繁,球員犯規及失誤也會相對提高,所以姑且以個人犯規及失誤兩項指標代表聯盟的防守強度,作為我們的自變數來建立泊松回歸模型。整理好的資料3pa.csv可以直接下載。
> three_pointer<-read.csv("c:/Users/USER/downloads/3pa.csv", header=T, sep=",")
> str(three_pointer)
'data.frame': 42 obs. of 4 variables:
$ season : chr "2020-2021" "2019-2020" "2018-2019" "2017-2018" ...
$ year : int 2021 2020 2019 2018 2017 2016 2015 2014 2013 2012 ...
$ a3p : int 2494 2408 2625 2378 2214 1975 1838 1766 1636 1213 ...
$ m3p : int 914 862 932 860 792 698 643 635 587 423 ...
$ foul : int 1389 1467 1714 1628 1632 1662 1658 1697 1626 1291 ...
$ turnover: int 996 1027 1155 1170 1144 1179 1177 1201 1192 962 ...
> attach(three_pointer)
載入上述資料後,用glm()指令並設定機率分配為poisoon與連結函數為log,建立泊松回歸模型。glm()是廣義線性模型的縮寫,標準語法是:glm(formula=, data=, family=機率分配(link=連結函數))。除了高斯常態分配、泊松分配外,還有下列參數可以搭配glm()使用:
Family | Link Function |
---|---|
binomial | link=logit |
gaussian | link=identity |
Gamma | link=inverse |
inverse.gaussian | link=1/mu^2 |
poisoon | link=log |
quasi | link=identity, variance=constant |
quasibinominal | link=logit |
quasipoisson | link=log |
泊松回歸模型
我們將建立三個泊松回歸模型,分別是三分球出手數與個人犯規、三分球出手數與失誤、三分球出手數與犯規+失誤。首先是三分球出手數與個人犯規的泊松回歸:
> model_1<-glm(formula=a3p~foul, data=three_pointer, family=poisson(link=log)) #模型1:三分球出手與個人犯規
> summary(model_1)
Call:
glm(formula = a3p ~ foul, family = poisson(link = log), data = three_pointer)
Deviance Residuals:
Min 1Q Median 3Q Max
-53.148 -15.120 3.015 10.216 34.538
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 9.733e+00 3.158e-02 308.24 <2e-16 ***
foul -1.528e-03 1.821e-05 -83.93 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 18664 on 41 degrees of freedom
Residual deviance: 12356 on 40 degrees of freedom
AIC: 12723
Number of Fisher Scoring iterations: 5
模型1三分球出手數與個人犯規的關係式為:
\[\log(\mu)=9.73-0.00153x_{foul}\]
模型1顯示個人犯規數與三分球出手數達顯著,係數是負數表示聯盟的犯規越少,整體三分球出手就越多。關於泊松模型的解釋,可以參考Poisson Step by Step。根據模型1可以計算當個人犯規達到1500次時的三分球出手數:
> foul1500<-model_1$coefficients[1]+model_1$coefficients[2]*1500 #依據模型1計算犯規達1500次時的三分球出手數
> exp(foul1500) #指數轉換
(Intercept)
1703.941
接下來我們來看看模型2:
> model_2<-glm(formula=a3p~turnover, data=three_pointer, family=poisson(link=log)) #模型2:三分球出手與失誤
> summary(model_2)
Call:
glm(formula = a3p ~ turnover, family = poisson(link = log), data = three_pointer)
Deviance Residuals:
Min 1Q Median 3Q Max
-57.492 -14.450 1.247 10.516 30.662
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 9.940e+00 3.142e-02 316.40 <2e-16 ***
turnover -2.365e-03 2.611e-05 -90.59 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 18664 on 41 degrees of freedom
Residual deviance: 11019 on 40 degrees of freedom
AIC: 11385
Number of Fisher Scoring iterations: 5
模型2三分球出手數與失誤的關係式為:
\[\log(\mu)=9.94-2.00237x_{turnover}\]
模型2的失誤一樣達到統計顯著,失誤的係數是負數,也表示失誤越高,三分球出手就越少。我們同樣可以根據模型2計算當失誤達到1500次時的三分球出手數:
> turnover1500<-model_2$coefficients[1]+model_1$coefficients[2]*1500 #依據模型2計算失誤達1500次時的三分球出手數
> exp(turnover1500) #指數轉換
(Intercept)
2096.378
接來看模型3:
> model_3<-glm(formula=a3p~foul+turnover, data=three_pointer, family=poisson(link=log)) #模型3:三分球出手與犯規+失誤
> summary(model_3)
Call:
glm(formula = a3p ~ foul + turnover, family = poisson(link = log), data = three_pointer)
Deviance Residuals:
Min 1Q Median 3Q Max
-56.765 -12.812 -0.104 9.412 29.044
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 9.873e+00 3.236e-02 305.07 <2e-16 ***
foul 1.042e-03 6.776e-05 15.38 <2e-16 ***
turnover -3.809e-03 9.797e-05 -38.88 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 18664 on 41 degrees of freedom
Residual deviance: 10781 on 39 degrees of freedom
AIC: 11150
Number of Fisher Scoring iterations: 5
模型3的關係式為:
\[\log(\mu)=9.87+0.00104x_{foul}-0.00381x_{turnover}\]
模型3犯規與失誤都達到統計顯著,但個人犯規的係數變成正數,失誤的係數一樣是負數。依據模型,我們可以利用predict.glm()來估計球季犯規與失誤平均達到1500次防守強度時的三分球出手數:
> predict.glm(model_3, type="response", newdata=data.frame(foul=1500, turnover=1500)) #預測犯規與失誤各達1500次時的三分球出手數
1
305.6643
殘差與模型配適度
觀測值跟擬合值之間的差距,也就是殘差,是決定模型良莠的關鍵。在診斷回歸模型時,經常使用皮爾森殘差(Pearson residuals)和偏差殘差(deviance residuals),當模型充分擬合數據時,它們等效且近似為標準常態分佈。。偏差殘差指的是當前模型的偏差與理想模型的最大偏差之間的差異。如果偏差夠小,代表模型擬合度高,模型擬合的卡方檢定將會不顯著。residuals()可用來檢視皮爾森殘差殘差與偏差偏差。
> residuals(model_3, type="pearson")
> residuals(model_3, type="deviance")
> with(model_3, cbind(res.deviance=deviance, df=df.residual, p.value=pchisq(deviance, df.residual, lower.tail=FALSE)))
res.deviance df p.value
[1,] 10781.13 39 0
泊松回歸:比例資料
有時候應變數並非計數而是比例,例如三分球命中率。只要用\(\mu\)除以\(t\),泊松迴歸很容易可以轉換為比例模型:
\[\log(\frac{\mu}{t})=\alpha x + \beta \]
\[\log(\mu)-\log(t)=\alpha x + \beta \]
\[\log(\mu)=\log(t) + \alpha x + \beta \]
其中\(\log(t) \)被稱為調整項(offset),上述式子經指數轉換後可得:
\[(\mu)=t \exp (\alpha x + \beta) \]
\[(\mu)=te^{\alpha x}e^{\beta}\]
先前我們已經有了三分球出手數,現在以三分球進球數當作調整項,在glm()中設定offset(),擬合三分球命中率與球季的泊松回歸模型:
> model_4<-glm(formula=a3p~year+offset(log(m3p)), data=three_pointer, family=poisson(link=log))
> summary(model_4)
Call:
glm(formula = a3p ~ year + offset(log(m3p)), family = poisson(link = log), data = three_pointer)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.2404 -0.7672 0.0544 1.1871 3.8609
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 8.4402757 0.9037810 9.339 < 2e-16 ***
year -0.0036823 0.0004502 -8.179 2.86e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 166.77 on 41 degrees of freedom
Residual deviance: 100.30 on 40 degrees of freedom
AIC: 466.73
Number of Fisher Scoring iterations: 3
計算單尾卡方檢定的模型配適度:
> 1-pchisq(q=model_4$deviance, df=model_4$df.residual)
[1] 4.35533e-07
泊松回歸的比例模型為:
\[\log(\mu)=\log(3p)-0.0037x_{year} + 8.44 \]