Afsnit 3.3: GG-teststørrelsen

Jeg vil nu formulere en hypotese i multinomialfordelingen generelt. Udgangspunktet er modellen
M0:(A1,…,Ak)∼multinom(n,(π1,…,πk)),πj≥0,π1+⋯+πk=1.(3.3.1) M_0:\enspace(A_1,\ldots,A_k)\sim\text{multinom}(n,(\pi_1,\ldots,\pi_k)),\enspace \pi_j\geq 0,\enspace\pi_1+\cdots+\pi_k=1.\quad \tag{3.3.1}
Hypotesen lægger begrænsninger på variationsområdet for (π1,…,πk),(\pi_1,\ldots,\pi_k), idet
πj=pj(θ),j=1,…,k,θ⊆Θ.(3.3.2) \pi_j=p_j(\theta),\enspace j=1,\ldots,k,\quad\theta\subseteq \Theta. \tag{3.3.2}
Her er pj(⋅)p_j(\cdot) kendte funktioner, θ\theta er en ukendt parameter, der skal estimeres ud fra data, og θ\theta kan variere i området Θ,\Theta, som indeholder et åbent område af Rd.\mathbf{R}^d. Det sidste udtrykker vi sprogligt på den måde, at θ\theta har dd frie parametre. Under hypotesen betegnes den statistiske model med M1,M_1, og man kan enten sige, at vi ønsker at teste hypotesen, eller at vi ønsker at teste reduktion fra model M0M_0 til model M1.M_1.
Som skøn over θ\theta bruges den værdi, der giver maksimum af likelihoodfunktionen L1(θ)L_1(\theta):
L1(θ)=L(p1(θ),…,pk(θ)), L_1(\theta)=L\big(p_1(\theta),\ldots,p_k(\theta)\big),
hvor L(⋅)L(\cdot) er givet i ligning (3.1.1).
Vi kan nu beregne likelihood ratio teststørrelsen Q,Q, som er forholdet mellem den maksimale værdi af likelihoodfunktionen under model M1M_1 og den maksimale værdi af likelihoodfunktionen under model M0.M_0. Ved samme beregning som i foregående afsnit finder vi
Q=L(p1(θ^),…,pk(θ^))L(A1n,…,Akn)=1(A1np1(θ^))A1⋯(Aknpk(θ^))Ak, Q=\frac{L\big(p_1(\hat\theta),\ldots,p_k(\hat\theta)\big)} {L\big(\frac{A_1}{n},\ldots,\frac{A_k}{n}\big)} =\frac{1}{\big(\frac{A_1}{np_1(\hat\theta)}\big)^{A_1} \cdots \big(\frac{A_k}{np_k(\hat\theta)}\big)^{A_k}},
og dermed
G=−2log⁡(Q)=2∑j=1kAjlog⁡(Ajej),ej=npj(θ^).(3.3.3) G=-2\log(Q)=2\sum_{j=1}^k A_j\log\big(\frac{A_j}{e_j}\big),\quad e_j=np_j(\hat\theta). \tag{3.3.3}
Her kaldes eje_j det forventede antal i kasse jj under hypotesen (under model M1M_1).
En lille værdi af QQ betyder, at data beskrives meget dårligere under model M1M_1 end under model M0.M_0. Jo mindre værdi af QQ jo mere kritisk for hypotesen. Dette er det samme som, at jo større GG er, jo mere kritisk. P-P\text{-}værdien for et test baseret på GG er derfor sandsynligheden for ved gentagelse af eksperimentet at få en værdi af G,G, der er større end eller lig med den faktisk observerede værdi af G.G. PP-værdien kan ikke beregnes eksakt, og i stedet benyttes en approksimation baseret på ki-i-anden-fordelingen, som I kender fra calculuskurset. Hvis den stokastiske variabel VV er ki-i-anden fordelt med ff frihedsgrader, skriver jeg dette som V∼χ2(f).V\sim\chi^2(f). Sandsynligheden for at V≤vV\leq v skrives som χcdf2(v,f)\chi^2_{\text{cdf}}(v,f) og beregnes i python med kommandoen st.chi2.cdf(v,f).
Resultat 3.3.1. (G-test)
Betragt multinomialmodellen (A1,…,Ak)∼multinom(n,(π1,…,πk)),(A_1,\ldots,A_k)\sim\text{multinom}(n,(\pi_1,\ldots,\pi_k)), πj≥0,\pi_j\geq 0, π1+⋯+πk=1\pi_1+\cdots+\pi_k=1 (model M0M_0) og hypotesen πj=pj(θ),\pi_j=p_j(\theta), j=1,…,k,j=1,\ldots,k, hvor θ\theta har dd frie parametre (model M1M_1). Betragt teststørrelsen G=−2log⁡(Q)=2∑j=1kAjlog⁡(Ajej),G=-2\log(Q)=2\sum_{j=1}^k A_j\log\big(\frac{A_j}{e_j}\big), ej=npj(θ^),e_j=np_j(\hat\theta), og lad GobsG_{\text{obs}} være den observerede værdi af teststørrelsen. Hvis alle de forventede ej=npj(θ^)e_j=np_j(\hat\theta) er større end eller lig med 5, har vi approksimativt
p-værdi=P(G≥Gobs)=1−χcdf2(Gobs,k−1−d). p\text{-værdi}=P(G\geq G_{\text{obs}})= 1-\chi^2_{\text{cdf}}(G_{\text{obs}},k-1-d).
Beviset for dette resultat er ikke nemt. Beviset bygger på den centrale grænseværdisætning (se afsnit 2.4) og en andenordens taylorudvikling af likelihoodfunktionen. På denne måde bliver teststørrelsen en sum af kvadrerede led, hvor hvert led approksimativt er normalfordelt. Antallet af frihedsgrader k−1−dk-1-d i χ2\chi^2-fordelingen er generelt d(M0)−d(M1),d(M_0)-d(M_1), hvor d(M0)d(M_0) og d(M1)d(M_1) er antallet af frie parametre i henholdsvis model M0M_0 og model M1.M_1. I model M0M_0 har vi bindingen, at π1+⋯+πk=1,\pi_1+\cdots+\pi_k=1, hvorfor antallet af frie parametre er k−1.k-1.
Resultatet ovenfor er første gang, vi støder på χ2\chi^2-fordelingen (ki-i-anden fordelingen) i dette kursus. Fordelingen optræder i sandsynlighedsdelen af jeres calculuskursus. For at I kan have en fornemmelse for denne fordeling, nævner jeg lige, at hvis U1,…,UfU_1,\ldots,U_f er ff uafhængige standard normalfordelte variable, så har U12+⋯+Uf2U_1^2+\cdots+U_f^2 en χ2-\chi^2\text{-}fordeling med ff frihedsgrader.
Eksempel 3.3.2. (Konstant dødsrate)
Vi vender tilbage til underafsnit 3.2.1 med data omkring dødstidspunkt for zebrafisk i en opløsning med sølvnanopartikler, og laver GG-testet for hypotesen beskrevet der. Først skal vi finde et skøn over parameteren θ,\theta, hvor πj=(1−θ)j−1θ\pi_j=(1-\theta)^{j-1}\theta, j≤4,j\leq 4, og π5=(1−θ)4.\pi_5=(1-\theta)^4. Likelihoodfunktionen bliver
L1(θ)=(9033,15,9,12,21)θ33((1−θ)θ)15((1−θ)2θ)9((1−θ)3θ)12((1−θ)4)21=(9033,15,9,12,21)θ69(1−θ)153.\begin{aligned} L_1(\theta)&=\binom{90}{33,15,9,12,21} \theta^{33} \big((1-\theta)\theta\big)^{15} \big((1-\theta)^2\theta\big)^{9} \big((1-\theta)^3\theta\big)^{12} \big((1-\theta)^4\big)^{21} \\ & =\binom{90}{33,15,9,12,21} \theta^{69} (1-\theta)^{153}. \end{aligned}
Ved sammenligning med likelihoodfunktionen i binomialmodellen i afsnit 2.1 ses, at θ^=69/(69+153)=0.3108.\hat\theta=69/(69+153)=0.3108. Dernæst beregnes de forventede antal som ej=90⋅(1−θ^)j−1θ^e_j=90\cdot (1-\hat\theta)^{j-1}\hat\theta, j≤4,j\leq 4, og e5=90⋅(1−θ^)4.e_5=90\cdot (1-\hat\theta)^4. Dette giver følgende tabel (forventede er afrundet til eˊ\acute{\text{e}}n decimal).
Tidsinterval0−2424−4848−7272−9696−Antal Døde331591221Forventede28.019.313.39.220.3 \begin{array}{lccccc} \hline \text{Tidsinterval} & 0-24 & 24-48 & 48-72 & 72-96 & 96- \\ \text{Antal Døde} & 33 & 15 & 9 & 12 & 21 \\ \text{Forventede} & 28.0 & 19.3 & 13.3 & 9.2 & 20.3 \\ \hline \end{array}
Da alle de forventede er større end eller lig med 5, beregner vi GG-teststørrelsen og den approksimative pp-værdi fra en χ2\chi^2-fordeling med 5−1−15-1-1 frihedsgrader. Ved beregning af antal frihedsgrader benyttes, at multinomialmodellen her deler op i 5 kasser, og den hypotese, der testes, har 1 parameter (nemlig θ\theta). Beregningen i kodevinduet nedenfor giver G=4.27G=4.27 og en pp-værdi på 0.234. Da pp-værdien er noget over 0.05, siger vi, at data ikke strider mod hypotesen om samme dødsrate i perioden 0-96 timer.

3.3.3 Beregning i python af G-test

Den følgende kode er kun semi-generel. Man skal selv indskrive en datavektor (aa), antal parametre (dpar) og beregne de forventede antal (ex).

χ2-fordelingen i python

I det følgende kodevindue tegnes tætheden for en χ2(df)\chi^2(\mathit{df}) fordeling, og 95%-fraktilen markeres med en lodret streg. Fraktilen angiver punktet, hvor 95% af sandsynligheden i fordelingen ligger til venstre for punktet og 5% ligger til højre for punktet. Prøv at køre koden med forskellige valg af antallet af frihedgrader df.\mathit{df}. Prøv også i kodevinduet at beregne sandsynligheden for at ligge til højre for 5.99 i en χ2\chi^2-fordeling med 2 frihedsgrader.

ForegåendeNæste