Modelowanie liczby bramek w Serie-A - rozkład Poissona i rozkład Skellam

Problem

Wszystkie dotychczas przedstawione sposoby modelowania wyników wydarzeń sportowych oparte były na regresji logistycznej. Modelowaliśmy siłę drużyn, także rozróżniając ich jakość w zależności od tego czy dana drużyna była gospodarzem. Dzisiejszy wpis zakłada inny sposób mierzenia jakości drużyn, w którym bierzemy pod uwagę liczbę strzelonych i straconych bramek. Liczba bramek jaką strzela drużyna i-ta jest uzależniona od siły ofensywnej drużyny $\alpha_i$, słabości defensywny przeciwnika $\beta_j$, ewentualnego wsparcia kibiców $\gamma$ oraz innych czynników, których nie znamy albo za krótko zastanawialiśmy się nad problemem żeby zjawisko do końca zrozumieć. Przy modelowaniu liczby bramek musimy założyć właściwy rozkład, który byłby dyskretny i nieujemny - rozkład Poissona. Rozkładem Poissona opisujemy intensywność badanych zdarzeń w określonej jednostce czasu. W naszym przypadku intensywnością jest liczba bramek strzelonych przez gospodarze $X_{i,j}$ oraz liczba bramek strzelonych przez gościa $Y_{i,j}$.

$X_{i,j} \sim Poisson(\alpha_i \beta_j \gamma)$

$Y_{i,j} \sim Poisson(\alpha_j \beta_i)$

Powyższy jest modelem bazowym, zaproponowanym przez Maher (1982), w którym liczba strzelonych i liczba straconych bramek przez gospodarza są niezależne od siebie i podążają rozkładem Poissona. Warto mieć na uwadze dodatkowe założenia $\alpha_i, \beta_i, \gamma > 0, \forall_i$.

Tak to robi R.

Rozwiązanie

Pomijam tym razem pierwszą część kody, która wiele nie różni się od poprzednich. Ponownie zaciągnięto dane z www.football-data.co.uk i po kilku operacjach doprowadzamy zbiór do takiej formy:

id	goals	atak	obrona	home
1	0	Chievo	Juventus	1
1	1	Juventus	Chievo	0
2	2	Roma	Fiorentina	1
2	0	Fiorentina	Roma	0
3	0	Atalanta	Verona	1
3	0	Verona	Atalanta	0

Dane w swojej formie odpowiadają formule, którą zamieściłem we wstępie, gdzie liczba goli w jednym meczu ($X_{i,j}$ i $Y_{i,j}$) to dwa rzędy z zbiorze. W związku z tym, że $X_{i,j}$ i $Y_{i,j}$ są od siebie niezależne przy szacowaniu są oddzielnymi obserwacjami. Do analizy użyjemy uogólnionego modelu liniowego z rozkładem Poissona (GLM). Anova zastosowana na modelu glm jest testem LR, który już opisywałem wcześniej. Szacujemy model i porządkujemy parametry wrzucając je do jednego data.frame. (Tak nawiasem pisząc - I <3 maggritr.)

Do prognozy wyników 23. kolejki Serie-A użyjemy wbudowanej funkcji predict(). Dane, na których dokonamy prognozy powinny mieć taką samą formę jak te, na których model szacowano. Wynik funkcji predict to logarytm intensywności rozkładu Poissona dla każdego $X_{i,j}$ $Y_{j,i}$ prognozowanych meczów. Intensywność rozkładu Poissona to inaczej średnia liczba bramek strzelonych przez daną drużynę w danym meczu. Dla każdego z prognozowanych meczów intensywności (h i a) wyglądają następująco.

Na podstawie intensywności możemy określić np. prawodpodobieństwo strzelenia określonej liczby bramek przez daną drużynę w danym meczu (dpois). Ważniejsze być może będzie dla nas rozstrzygnięcie spotkania czyli różnica w strzelonych golach. Ponieważ rozstrzygnięcie jest różnicą pomiędzy dwoma rozkładami Poissona, rozkład wyniku będzie miał postać rozkładu Skellama. Skorzystamy z funkcji dskellam() w pakiecie VGAM do określenia prawdopodobieństwa danej różnicy bramek. Ponieważ zwycięstwo gospodarza to różnica $X_{i,j} - Y_{j,i} \geqslant 1$ musimy zsumować gęstości dla każdego z możliwych wyników, w którym wygrywa gospodarz (od 1 do 30 powinno wystarczyć). Podobnie robimy dla porażek i dla zwycięstw.

id	home	away	h	a	gosc_p	remis_p	gosp_p
1	Sassuolo	Fiorentina	1,09494072	1,341957729	42,36%	27,15%	30,49%
2	Palermo	Napoli	1,660903345	2,07955973	47,76%	21,05%	31,19%
3	Milan	Empoli	1,468506037	0,816874635	20,89%	26,56%	52,55%
4	Atalanta	Inter	0,952321309	1,290657643	44,27%	28,26%	27,47%
5	Genoa	Verona	1,990845919	0,850722982	15,23%	20,77%	64,00%
6	Roma	Parma	2,955615659	0,501645062	3,82%	9,87%	86,31%
7	Torino	Cagliari	1,814261385	0,968426876	19,88%	22,87%	57,24%
8	Udinese	Lazio	1,237176484	1,454853572	42,17%	25,68%	32,15%
9	Chievo	Sampdoria	0,669219856	0,994739866	42,22%	33,81%	23,97%
10	Cesena	Juventus	0,382502684	2,742712365	86,40%	10,36%	3,24%

Tak prezentują się parametry siły dla poszczególnych zespołów.

link do pełnego skryptu

link do funkcji

Literatura:

Maher, M.J. (1982) Modelling association football scores.Statistica Neerlandica36, 109-118