読者です 読者をやめる 読者になる 読者になる

Rを通じて統計学を学ぶ備忘録ブログ

SPSSからRに移行したい私のような人向けのR解説ブログ兼学習用備忘録。

第74回 ポアソン分布とサッカーのゴール数

R GLM サッカー

今回は、ポアソン分布とサッカーのゴール数の関係をもとに、基本的なポアソン分布の説明から、オフセットと過分散も検討したポアソン回帰までみていきます。

 

ポアソン分布

サッカーのゴール数がポアソン分布に従うことは有名な話ですね。そもそもポアソン分布とは、イベントが生じた回数を説明する確率分布として、よく使われます。いわゆるカウントデータですね。λを決めると形が決まり、λ=平均=分散という特徴を持ちます。詳しくは後で書きますが、平均≠分散となる場合が過分散という状態です。平均が大きくなると、分散も大きくなるので、分布の形が右に広がっていきます。

f:id:teruaki-sugiura:20160131132909p:plain   E[X]=\lambda \,   V[X]=\lambda \,

> par(mfrow = c(2, 2))
> x <- 0:10
> plot(dpois(x,1),type="o")
> plot(dpois(x,2),type="o")
> plot(dpois(x,5),type="o")
> plot(dpois(x,7),type="o")

f:id:teruaki-sugiura:20160131132727p:plain

ポアソン分布の導出については、p=λ/nの二項分布において、λ を一定にしてnを∞に飛ばすと、その分布がポアソン分布になります。これをポアソンの極限定理です。

\lim_{\lambda=np, n\to\infty}{n \choose k} p^k (1-p)^{n-k}= \frac{\lambda^k e^{-\lambda}}{k!}

\lim_{n\to\infty}\left(1-\frac{\lambda }{n}\right)^n=e^{-\lambda}

\lim_{n\to\infty} P(X=k)=\lim_{n\to\infty}{n \choose k} p^k (1-p)^{n-k}
=\lim_{n\to\infty}{n! \over (n-k)!k!} \left({\lambda \over n}\right)^k \left(1-{\lambda\over n}\right)^{n-k}

=\lim_{n\to\infty} \underbrace{\left({n \over n}\right)\left({n-1 \over n}\right)\left({n-2 \over n}\right) \cdots \left({n-k+1 \over n}\right)} \underbrace{\left({\lambda^k \over k!}\right)}\underbrace{\left(1-{\lambda \over n}\right)^n}\underbrace{\left(1-{\lambda \over n}\right)^{-k}}

nを∞に飛ばすと、最初の括弧は1に、2番目の括弧は n がないのでそのまま、3番目の括弧はe^-λ、最後の部分は1になります。では、実際にはどのような場面で使うのでしょうか。大味に仮定をまとめると、「時間を区切って、各時間帯で発生するイベントは1回だけで、他の時間帯に発生したイベントの影響は受けない」というものなのですが、いまいちわかりにくいので、具体例をみていきましょう。

  • 1時間に特定の交差点を通過する車両の台数。
  • 単位面積あたりの雨粒の数。
  • 1ページの文章を入力するとき、綴りを間違える回数。
  • 1日に受け取るメールの件数。
  • ある一定の時間内の店への来客数。
  • 1分間のwebサーバへのアクセス数。
  • 1キロメートルあたりのある通り沿いのコンビニとかレストランの軒数。
  • 1ヘクタールあたりのエゾマツの本数。
  • 1立方光年あたりの恒星の数。    etc

サッカーを例にするならば、「90分中に平均λ回生まれるゴールする事象」に相当します。1時間あたりのwikipediaの最近更新したページの編集数もポアソン分布に従うらしい、、、です。

ポアソン分布とサッカーのゴール数

では、Rを使いながらポアソン分布とサッカーのゴール数の関係をガンバ大阪の2015年のゴール数を例に見ていきましょう。このシーズンの分布は、以下のようになります。

ゴール数0の試合数:5

ゴール数1の試合数:10

ゴール数2の試合数:12

ゴール数3の試合数:6

ゴール数4の試合数:1

ゴール数5の試合数:0

> game <- c(0:5) #ゴール数iの試合数
> goal <- c(5, 10, 12, 6, 1, 0) #各ゴール数の試合数
> goal.sum <- sum(goal)
> s <- game %*% goal 
> r <- s / goal.sum #最尤推定法でλを求める

λの推定値r=1.64となりました。これをもとに近似的に95%信頼区間を求めてみます。

> conf <- c(r - 2 * sqrt(r / goal.sum), r + 2 * sqrt(r / goal.sum))

conf=1.20〜2.08となりましたので、λの推定値は1.20以上2.08以下であることがわかります。 そして、このλを使って、ガンバ大阪のゴールの期待値を求めます。観測度数を棒グラフとして、期待度数を線グラフで書きます。度数を求めるので、確率にゴール数の合計をかけることになります。ある程度、近似しています。

> names(goal) <- c(0:5)
> y <- dpois(c(0:5), lambda = r) * sum(goal)
> lines(barplot(goal, ylim = c(0, 15), col = "darkgreen"), y, type = "o" )  

f:id:teruaki-sugiura:20160131140454p:plain

この結果をもとにガンバ大阪のゴール数はポアソン分布に従うと言えるのでしょうか。適合度の検定を行って、チェックしていきましょう。

以下のように仮説を設定します。

・帰無仮説:この観測されたデータがポアソン分布に従う

・対立仮説:この観測されたデータがポアソン分布に従わない

goodfit()関数のxには行列(度数、ゴール数iの試合数)を入れます。ポアソン分布のチェックなのでtype = "poisson"、方法はカイ2乗値を最小にする方法を用いるので、method = "MinChisq"を指定。

> library(vcd)

> x <- matrix(c(5, 10, 12, 6, 1, 0, c(0:5)), nr =6)

> result <- goodfit(x, type = "poisson", method = "MinChisq")
> summary(result)
> plot(result)

 

f:id:teruaki-sugiura:20160131141229p:plain

 

 

カイ2乗値は3.06でp値は0.54であるため、帰無仮説(この観測されたデータがポアソン分布に従う)を棄却できず、対立仮説を採択することになります。しかし、実際にRを動かす、 Chi-squared approximation may be incorrectというエラーが出ます。つまり、期待度数が小さいカテゴリーがあるために、カイ2乗分布の近似がよくない可能性があります。では最後に、観測されたデータが、期待度数とどの程度異なるのかを視覚化してみましょう。つまり、期待度数に合うように観測されたデータの度数を上下させます。したがって、ポアソン分布に従わないデータであれば、この棒グラフの変動が大きくなることで、ズレを確認できますし、もはやそのデータはポアソン分布でなく、違う確率分布に従うということがわかります。

f:id:teruaki-sugiura:20160131142546p:plain

 

ポアソン回帰分析

次は、ガンバ大阪に限らず、ゴール数の期待値をポアソン回帰分析を用いて、計算していきましょう。ポアソン回帰分析とは目的変数がイベントの回数や個体などのカウントデータの時に利用されます。つまり、目的変数の誤差構造がポアソン分布に従う場合に利用します。そして、リンク関数には「対数」が用いられます。この理由は簡単で、カウントデータはマイナスの値を取りませんし、期待値がマイナスにならないようにするためです。

では、ゴール数をポアソン回帰分析で計算していきましょう。ここでのデータはJ1の18チームのゴール数をカウントしたものを用います。まずは「チーム」による違いに着目した分析を行います。

> j.one <- read.csv("j.one.csv", header=TRUE)

> attach(j.one)
> poisson <- glm(goal ~ team, family = "poisson", data = j.one)
> summary(poisson) 

基準カテゴリは「ガンバ大阪」です(基準カテゴリはアルファベット順もしくは値の小さいものがなるようです)。そのため、ガンバ大阪との違いをパラメタとして推定しています。結果は以下のような図(左:ポアソン回帰、右:負の二項分布のオフセット付)になりました。広島、鹿島、川崎、浦和の4チームがガンバ大阪よりもゴール数の期待値が大きくなり、山形、甲府あたりはかなり期待値が低いようです。しかし、この解析には少々問題があります。こんな単純なモデルが役に立つわけがない!!というような仮説の荒さではなく、統計モデルとしての理論的な問題があります。それが、オフセットと過分散です。この2つについてみていきましょう。

f:id:teruaki-sugiura:20160131163524p:plain

オフセット

ポアソン回帰分析の目的変数は、あるイベントがどれだけの時間や面積の大きさにおいて発生しているのかによって、意味が異なります。

例えば、サッカーのFK(フリーキック)を例に考えてみると、「10本蹴って1回成功の場合」と「30本蹴って2回成功の場合」があるとすると、「30本蹴って2回成功の場合」のほうが回数が多いが、「30本蹴って2回成功の場合」が「10本蹴って1回成功の場合」よりも頻繁に成功しているとは考えにくいですよね。なので、その回数などが異なっている場合、その回数をちゃんと評価してあげるために、オフセット項を説明変数(回帰係数は「1」)として追加する必要があります。ゴール数もシュートをたくさん打てば、たくさん入るということを無視して解析するのではなく、シュートの本数も加味した上で、ゴールの質を評価しましょうということです。

言い換えると、交通事故件数を評価する際に、人口も加味するというやつです。人口が多い県は、交通事故件数が増えるのは当たり前の話なので。

オフセットの考え方は以下のようなもので、わりと単純です。

f:id:teruaki-sugiura:20160131170357p:plain

対数の目的変数(分数)として、分析するのではなく、対数の性質を利用して、割算を引き算にして、移項するだけです。次は過分散についてみていきましょう。

過分散

ポアソン分布の基本的な性質として、期待値=分散というものがありますが、今回の例では、期待値と分散が等しくなく、分散がかなり大きい状態です。つまり過分散なのです。逆に分散が小さいことを過小分散といいますが、多くの場合、過分散ばかりです。それでは、なぜ過分散が問題なのでしょうか。これはポアソン回帰分析であれば分散を実際よりも小さいものとして計算していくことになるので、目的変数への影響を過大評価することなってしまいます。そのため、このような状態の場合、ポアソン回帰分析ではなく、擬似ポアソン回帰や負の二項回帰などを利用することになります。

 

負の二項回帰分析

負の二項分布はポアソン分布から導出されます。ポアソン分布のパラメタ λ がガンマ分布Ga(α^-1, μα) に従うとすると、以下のように書き換えることができます。

f:id:teruaki-sugiura:20160131192751p:plain

以上により、負の二項分布が導かれました。そして、この確率分布を使って、回帰分析を行うのが負の二項回帰分析です。

f:id:teruaki-sugiura:20160131183142p:plain

 μ(μ > 0)は Y の平均値で、α(α > 0)は不均一パラメータです。そして、負の二項回帰は以下のようなモデルです。

f:id:teruaki-sugiura:20160131183451p:plain

i行のXをx_iとして、指数を取ると、以下のように書き換えることができます。あとは対数尤度関数のα、βを最尤法を用いて最大化すると、パラメタの推定値が得られます。

f:id:teruaki-sugiura:20160131183649p:plain

> ngbinomial <- glm.nb(goal ~ team + offset(log(shoot)), data = j.one)
> summary(ngbinomial)

結果の図を再掲します。右側が、負の二項分布でシュート数をオフセット項として追加した場合の結果です。 

f:id:teruaki-sugiura:20160131163524p:plain

以上の結果よりも、ガンバ大阪よりも広島、柏、川崎、東京、浦和はゴールの数(質)が良く、山形はあまりよくないようですね。加えて、AICもかなり小さくなっています。補強するなら、広島、柏、川崎、東京、浦和あたりからストライカーを獲得した法が良いかもしれませんね。もっとも、より詳細に検討する必要はありますが。

 

以上で今回はおしまい。

広告を非表示にする