Rozdział 6 Rozkład ujemny dwumianowy
6.1 Funkcja gęstości
Rozkład ujemny dwumianowy (zwany też rozkładem Pascala) został zaimplementowany do funkcji scipy.stats.nbinom i można go przedstawić za pomocą wzoru:
\[\begin{equation}
f(x\;|\;r,p)={x+r-1\choose r-1}p^r(1-p)^x,\quad x=0,1,...,\quad r\in N
\tag{6.1}
\end{equation}\]
gdzie: \(r\) jest liczbą sukcesów, \(x\) jest liczbą niepowodzeń tj. liczba zdarzeń poprzedzających \(r\) sukcesów, a \(p\) jest prawdopodobieństwem niepowodzeń.
Do powyższego wzoru (6.1) można zastosować alternatywny zapis współczynnika dwumianu: \[\begin{equation} f(x\;|\;r,p)={x+r-1\choose x}p^r(1-p)^x,\quad x=0,1,...,\quad r>0 \tag{6.2} \end{equation}\] dzięki któremu w rozkładzie ujemnym dwumianowym (zwanym też rozkładem Polya) można przyjąć, że parametr \(r>0\). Dodatkowo współczynnik dwumianowy można zapisać w oparciu o funkcję gamma: \[\begin{equation} f(x\;|\;r,p)=\frac{\Gamma(x+r)}{\Gamma(r)\Gamma(x+1)}p^r(1-p)^x,\quad x=0,1,...,\quad r>0 \tag{6.3} \end{equation}\] gdzie: \(E(X)=r(1-p)/p\) oraz \(V(X)=r(1-p)/p^2\).
import scipy.stats as stats
import numpy as np
x = stats.nbinom.rvs(n=5.7, p=0.3, size=10000, random_state=2305)
r = np.mean(x)**2/(np.var(x,ddof=1)-np.mean(x))
p = np.mean(x)/np.var(x,ddof=1)
print("MOM: r= %.4f, p= %.4f" % (r,p))## MOM: r= 5.5556, p= 0.2948Jeżeli przyjmiemy, że parametr \(r=\phi\) oraz \(p=\frac{\phi}{\mu+\phi}\) to mieszanka rozkładu Poissona-Gamma będzie miała postać: \[\begin{equation} f(x\;|\;\mu,\phi) =\frac{\Gamma(x+\phi)}{\Gamma(\phi)\Gamma(x+1)}\left(\frac{\phi}{\mu+\phi}\right)^{\phi}\left(\frac{\mu}{\mu+\phi}\right)^{x},\quad x=0,1,...,\quad \phi>0 \tag{6.4} \end{equation}\] gdzie: \(E(X)=\mu\) oraz \(V(X)=\mu+\phi^{-1}\mu^2\).
import scipy.stats as stats
import numpy as np
from scipy.optimize import minimize
def rn(mu, phi, n, rand):
r = phi # phi_nb2
p = r/(mu+r)
return stats.nbinom.rvs(n=r, p=p, size=n, random_state=rand)
x = rn(mu = 1.5, phi = 5, n = 10000, rand = 2305)
mu = np.mean(x)
phi = np.mean(x)**2/(np.var(x,ddof=1)-np.mean(x))
print("MOM: mean= %.4f, phi_NB2= %.4f" % (mu,phi))
def L_nb2(par):
phi = par[0]
mu = par[1]
logLik = -np.sum( stats.nbinom.logpmf(x, n=phi, p=phi/(phi+mu)) )
return(logLik)
initParams = [1,1]
res = minimize(L_nb2, initParams, method= "Nelder-Mead")
print("MLE: mean= %.4f, phi_NB2= %.4f, logLik= %.2f" % (res.x[1],res.x[0],L_nb2(res.x)))## MOM: mean= 1.5083, phi_NB2= 5.1984
## MLE: mean= 1.5083, phi_NB2= 5.1768, logLik= 16106.60Gdy do wzoru (6.4) podstawimy \(\phi=\alpha^{-1}\) i dokonamy prostych przekształceń to otrzymamy: \[\begin{equation} f(x\;|\;\mu,\alpha)=\frac{\Gamma(x+\alpha^{-1})}{\Gamma(\alpha^{-1})\Gamma(x+1)}\left(\frac{1}{\alpha\mu+1}\right)^{\alpha^{-1}}\left(\frac{\alpha\mu}{\alpha\mu+1}\right)^x \tag{6.5} \end{equation}\] gdzie: \(E(X)=\mu\) oraz \(V(X)=\mu+\alpha\mu^2\).
import scipy.stats as stats
import numpy as np
def rn(mu, alpha, n, rand):
r = 1/alpha # phi_nb2
p = r/(mu+r)
return stats.nbinom.rvs(n=r, p=p, size=n, random_state=rand)
x = rn(mu = 1.5, alpha = 9, n = 10000, rand = 2305)
mu = np.mean(x)
alpha = (np.var(x,ddof=1)-np.mean(x))/np.mean(x)**2
print("MOM: mean= %.4f, alpha_NB2= %.4f" % (mu,alpha))## MOM: mean= 1.5355, alpha_NB2= 9.1284Po zlogarytmowaniu wyrażenia (6.5) otrzymamy funkcję logarytmu wiarygodności o postaci: \[\begin{equation} L(x\;|\;\mu,\alpha)=\ln\Gamma(x+\alpha^{-1})-\ln\Gamma(\alpha^{-1})-\ln\Gamma(x+1)-\alpha^{-1}\ln(1+\alpha\mu)-x\ln(1+\alpha\mu)+x\ln(\alpha\mu) \tag{6.6} \end{equation}\]
6.2 Liniowy model ujemnej dwumianowej regresji
Do modelowania zmiennych licznikowych można stosować regresję Poissona jeśli zostało spełnione założenie równości średniej i wariancji. W przypadku wystąpienia zjawiska zawyżonej dyspersji tj. \(E(X)\leq V(X)\) nie można prawidłowo wyznaczyć błędów standardowych ocen parametrów. Rozwiązaniem tego problemu może być zastosowanie takich rozkładów które są przeznaczone do modelowania nadmiernie rozproszonych danych. Jedną z wielu propozycji (Coly i in. 2016) jest rozkład ujemny dwumianowy z liniową (NB1) lub kwadratową (NB2) zależnością między średnią a wariancją (Cameron i Trivedi 1998): \[\begin{equation} V(X)=\mu+\alpha\mu^p \tag{6.7} \end{equation}\]
model NB1 dla \(p=1\): \[\begin{equation} E[(y_i-\mu_i)^2]=\phi\mu_i\quad\longrightarrow\quad\phi = E[(y_i-\mu_i)^2/\mu_i] \tag{6.8} \end{equation}\] Estymator parametru \(\phi\) po zastosowaniu korekty: \[\begin{equation} \hat{\phi}_{\mathrm{NB1}}=\frac{1}{n-k}\sum_{i=1}^{n}\frac{(y_i-\hat{\mu}_i)^2}{\hat{\mu}_i}\quad\mathrm{gdzie}\quad\hat{\alpha}_{\mathrm{NB1}}=\hat{\phi}_{\mathrm{NB1}}-1 \tag{6.9} \end{equation}\]
model NB2 dla \(p=2\): \[\begin{equation} E[(y_i-\mu_i)^2-\mu_i]=\alpha\mu_i^2\quad\longrightarrow\quad\alpha = E[\{(y_i-\mu_i)^2-\mu_i\}/\mu_i^2] \tag{6.10} \end{equation}\] Estymator parametru \(\alpha\) po zastosowaniu korekty: \[\begin{equation} \hat{\alpha}_{\mathrm{NB2}}=\frac{1}{n-k}\sum_{i=1}^{n}\frac{(y_i-\hat{\mu}_i)^2-\hat{\mu}_i}{\hat{\mu}^2_i}\quad\mathrm{gdzie}\quad\hat{\phi}_{\mathrm{NB2}}=1/\hat{\alpha}_{\mathrm{NB2}} \tag{6.11} \end{equation}\]
W modelu który uwzględnia nadmierną dyspersje (model quasi-Poissona) błędy standardowe są modyfikowane w oparciu o wzór: \[\begin{equation} SE_{Q}(\beta)=SE_{Pois}(\beta)\cdot \sqrt{\hat{\phi}_{\mathrm{NB1}}} \tag{6.12} \end{equation}\] Warto podkreślić, że w równaniu (6.12) jest wykorzystany estymator \(\hat{\phi}_{\mathrm{NB1}}\) który jest powszechnie stosowany do szacowania nadmiernej dyspersji w modelu Poissona (McCullagh i Nelder 1989). Dodajmy jeszcze, że estymator dyspersji (6.9) można zapisać w alternatywny sposób: \[\begin{equation} \hat{\phi}_{\mathrm{NB1}}=\frac{\chi^2_P}{n-k} \tag{6.13} \end{equation}\] gdzie: \(\chi^2_{P}\) to suma kwadratów reszt Pearsona która jest często stosowana do oceny dobroci dopasowania modelu. Reszty Pearsona wyznaczamy na bazie regresji Poissona za pomocą wzoru: \[\begin{equation} r_{i}^P=\frac{y_i-\hat{u}_i}{\sqrt{\hat{u}_i}} \tag{6.14} \end{equation}\]
Wykorzystanie funkcji statsmodels.discrete.discrete_model.NegativeBinomial umożliwia oszacowanie modelu NB1 lub NB2 (opcja domyślna) oraz parametru \(\alpha\). Zastosowanie liniowego modelu ujemnej dwumianowej regresji zostanie zaprezentowane na przykładzie zestawu danych nb_data.
import pandas as pd
import statsmodels.api as sm
import patsy
df = pd.read_stata('https://stats.idre.ucla.edu/stat/stata/dae/nb_data.dta')
df['prog'].replace([1.0,2.0,3.0],['General','Academic','Vocational'],inplace=True)
model = 'daysabs ~ math + C(prog, Treatment(reference="General"))'
y, x = patsy.dmatrices(model, df, return_type='dataframe')
x.columns = ['Intercept', 'Academic', 'Vocational','math']
nb = sm.NegativeBinomial(y,x,loglike_method='nb2').fit(disp=0)
a = nb.params.values[4]
print(nb.summary())
print("\nphi: ",round(1/a, 8), ", sqrt_phi: ", round((1/a)**0.5, 8))## NegativeBinomial Regression Results
## ==============================================================================
## Dep. Variable: daysabs No. Observations: 314
## Model: NegativeBinomial Df Residuals: 310
## Method: MLE Df Model: 3
## Date: Mon, 19 Aug 2019 Pseudo R-squ.: 0.03441
## Time: 19:17:04 Log-Likelihood: -865.63
## converged: True LL-Null: -896.47
## Covariance Type: nonrobust LLR p-value: 2.563e-13
## ==============================================================================
## coef std err z P>|z| [0.025 0.975]
## ------------------------------------------------------------------------------
## Intercept 2.6153 0.196 13.319 0.000 2.230 3.000
## Academic -0.4408 0.183 -2.414 0.016 -0.799 -0.083
## Vocational -1.2787 0.202 -6.331 0.000 -1.675 -0.883
## math -0.0060 0.003 -2.390 0.017 -0.011 -0.001
## alpha 0.9683 0.100 9.729 0.000 0.773 1.163
## ==============================================================================
##
## phi: 1.03271369 , sqrt_phi: 1.01622522Za pomocą wybranej pomocniczej regresji liniowej: \[\begin{equation} w_i =\hat{\alpha}_{\,\mathrm{NB1}}+\epsilon_i \tag{6.15} \end{equation}\] \[\begin{equation} w_i=\hat{\alpha}_{\,\mathrm{NB2}}\,\hat{u}_i+\epsilon_i \tag{6.16} \end{equation}\] można weryfikować hipotezy statystyczne: \[\begin{equation} H_0:\;\alpha_{\,\mathrm{NB1}}=0\quad\mathrm{vs}\quad H_1:\;\alpha_{\,\mathrm{NB1}}\neq0 \tag{6.17} \end{equation}\] \[\begin{equation} H_0:\;\alpha_{\,\mathrm{NB2}}=0\quad\mathrm{vs}\quad H_1:\;\alpha_{\,\mathrm{NB2}}\neq0 \tag{6.18} \end{equation}\] Warto podkreślić, że \(\phi_{\,\mathrm{NB1}}=\alpha_{\,\mathrm{NB1}}+1\) więc hipotezę (6.17) można przedstawić jako: \[\begin{equation} H_0:\;\phi_{\,\mathrm{NB1}}=1\quad\mathrm{vs}\quad H_1:\;\phi_{\,\mathrm{NB1}}\neq1 \tag{6.19} \end{equation}\] Dodatkowo wyniki uzyskane za pomocą regresji (6.15) są tożsame wynikami uzyskanymi na podstawie wzorów: \[\begin{equation} \hat{\alpha}_{\mathrm{NB1}}=E(w_i)\quad\mathrm{oraz}\quad SE_{\hat{\alpha}_{1NB}}=\sqrt{V(w_i)/n} \tag{6.20} \end{equation}\] gdzie: \[\begin{equation} w_i=\frac{(y_i-\hat{u}_i)^2-y_i}{\hat{u}_i} \tag{6.21} \end{equation}\]
import warnings
warnings.filterwarnings("ignore")
import pandas as pd
import statsmodels.api as sm
import patsy
df = pd.read_stata('https://stats.idre.ucla.edu/stat/stata/dae/nb_data.dta')
df['prog'].replace([1.0,2.0,3.0],['General','Academic','Vocational'],inplace=True)
model = 'daysabs ~ math + C(prog, Treatment(reference="General"))'
y, x = patsy.dmatrices(model, df, return_type='dataframe')
x.columns = ['Intercept', 'Academic', 'Vocational','math']
yp = sm.Poisson(y,x).fit(disp=0).predict()
w = ((y['daysabs']-yp)**2-y['daysabs'])/yp
n1 = sm.OLS(w,yp*0+1).fit(use_t=1) # regresja OLS dla alpha_NB1
n1.model.data.xnames = ['alpha_NB1']
n2 = sm.OLS(w,yp).fit(use_t=1) # regresja OLS dla alpha_NB2
n2.model.data.xnames = ['alpha_NB2']
print(n1.summary().tables[1])
print("phi_NB1: ",n1.params.values[0]+1,'\n')
print(n2.summary().tables[1])
print("phi_NB2: ",1/n2.params.values[0])## ==============================================================================
## coef std err t P>|t| [0.025 0.975]
## ------------------------------------------------------------------------------
## alpha_NB1 5.5105 0.768 7.178 0.000 4.000 7.021
## ==============================================================================
## phi_NB1: 6.51052636765949
##
## ==============================================================================
## coef std err t P>|t| [0.025 0.975]
## ------------------------------------------------------------------------------
## alpha_NB2 0.7986 0.117 6.835 0.000 0.569 1.028
## ==============================================================================
## phi_NB2: 1.2522191233195814Funkcja statsmodels.genmod.families.family.NegativeBinomial umożliwia estymację modelu NB2 dla ustalonej wartości \(\alpha\). Warto zwrócić uwagę, że za pomocą tej funkcji możemy w sposób symulacyjny dobrać odpowienik parametr \(\alpha\). Dodajmy jeszcze, że do modelowania danych licznikowych można wykorzystać złożony rozkład Poissona–gamma który jest szczególnym przypadkiem rozkładu Tweedie:
\[\begin{equation}
E(X) = \mu \quad\mathrm{oraz}\quad Var(X) = \phi\mu^p
\tag{6.22}
\end{equation}\]
W zależności od wartości parametru kształtu \(p\) można otrzymać
kilka znanych rozkładów jako szczególne przypadki dystrybucji Tweedie:
\(p = 0\) - rozkład normalny,
\(0 < p < 1\) - rozkład nie jest zdefiniowany,
\(p = 1\) - rozkład Poissona,
\(1 <p <2\) - rozkład Poissona–gamma,
\(p = 2\) - rozkład gamma,
\(2 <p <3\) - dodatnie rozkłady stabilne,
\(p = 3\) - odwrotny rozkład Gaussa / rozkład Walda,
\(p> 3\) - dodatnie rozkłady stabilne,
\(p =\infty\) - ekstremalne stabilne rozkłady.
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm
df = pd.read_stata('https://stats.idre.ucla.edu/stat/stata/dae/nb_data.dta')
df['prog'].replace([1.0,2.0,3.0],['General','Academic','Vocational'],inplace=True)
model = 'daysabs ~ math + C(prog, Treatment(reference="General"))'
B = 100
X = np.linspace(1.0001,1.9999,B)
n = [smf.glm(model, data=df,\
family = sm.families.Tweedie(link = sm.families.links.log, var_power=i)).fit() for i in X]
res = [n[i].deviance for i in range(B)]
sol = n[np.argmin(res)]
sol.model.data.xnames = ['Inercept','Academic','Vocational','math']
print(sol.summary())
print("\np: ",X[np.argmin(res)])## Generalized Linear Model Regression Results
## ==============================================================================
## Dep. Variable: daysabs No. Observations: 314
## Model: GLM Df Residuals: 310
## Model Family: Tweedie Df Model: 3
## Link Function: log Scale: 2.3160
## Method: IRLS Log-Likelihood: nan
## Date: Mon, 19 Aug 2019 Deviance: 902.82
## Time: 19:17:09 Pearson chi2: 718.
## No. Iterations: 10
## Covariance Type: nonrobust
## ==============================================================================
## coef std err z P>|z| [0.025 0.975]
## ------------------------------------------------------------------------------
## Inercept 2.6258 0.192 13.699 0.000 2.250 3.001
## Academic -0.4410 0.178 -2.484 0.013 -0.789 -0.093
## Vocational -1.2779 0.203 -6.305 0.000 -1.675 -0.881
## math -0.0062 0.003 -2.423 0.015 -0.011 -0.001
## ==============================================================================
##
## p: 1.6363363636363637Bibliografia
Cameron, A.C., i P.K. Trivedi. 1998. Regression Analysis of Count Data 2nd Edition. Cambridge University Press. https://assets.cambridge.org/97805216/32010/frontmatter/9780521632010_frontmatter.pdf.
Coly, Sylvain, Anne-Franoise Yao, David Abrial, i Myriam Charras-Garrido. 2016. „Distributions to model overdispersed count data”. Journal de la Societe Francaise de Statistique 157 (2). Societe Francaise de Statistique: 39–64. http://journal-sfds.fr/article/view/556.
McCullagh, P., i J.A. Nelder. 1989. Generalized Linear Models, Second Edition. Chapman and Hall/CRC Monographs on Statistics and Applied Probability Series. Chapman & Hall. https://www.amazon.com/Generalized-Chapman-Monographs-Statistics-Probability/dp/0412317605.