Tema 6 Contrastos de bondat d’ajust

Els contrastos d’un paràmetre del tema anterior els empràvem bàsicament per cercar evidència que a un paràmetre poblacional (una mitjana, una proporció) “li passava qualque cosa”:

  • Tenim evidència que la mitjana poblacional és diferent de 2?

  • Tenim evidència que la proporció poblacional d’èxits és més gran que el 10%?

Però també podem entendre que hi cercàvem si podem acceptar que el paràmetre poblacional de la variable que havia produit la mostra era igual a qualque cosa:

  • Podem acceptar que la mitjana poblacional és igual de 2?

  • Podem acceptar que la proporció poblacional d’èxits és del 10%?

Als contrastos de bondat d’ajust anam una passa més enllà, i ens demanam si podem acceptar que la variable poblacional té una distribució concreta, o si per contra seria molt estrany que la nostra mostra hagués estat produïda per una variable amb aquella distribució. Per exemple:

  • Llençam un dau en l’aire moltes vegades, apuntam els resultats, i ens demanam si podem acceptar que el dau és honrat (és a dir, que tots els resultats són equiprobables) o si per contra hi ha evidència que està trucat.

  • Anotam els ingressos diaris per una determinada malaltia en un hospital i ens demanam si podem acceptar que segueixen una distribució de Poisson o si per contra hi ha evidència que no és el cas, la qual cosa seria senyal que hi ha un brot estrany de la malaltia.

  • Volem emprar una mostra petita en un t-test. Podem acceptar que prové d’una distribució normal?

De fet, ja hem estudiat un tipus de contrast de bondat d’ajust: els contrastos bilaterals d’una proporció, on ens demanàvem si podem acceptar que la mostra prové de tal variable Bernoulli o si per contra hi ha evidència que no.

6.1 Bondat d’ajust

Considerau l’exemple següent:

Exemple 6.1 Les freqüències de les darreres xifres dels primers premis de la Grossa de la loteria de Nadal als 211 sortejos de la seva història han estat les següents:

\[ \begin{array}{r|cccccccccc} \hline \text{xifra} & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 \\ \hline \text{freqüència} & 22 & 8 & 13 & 21 & 27 & 31 & 27 & 22 & 23 & 17\\ \hline \end{array} \]

Hi ha evidència que no totes les xifres tenen la mateixa probabilitat de sortir com a darrera xifra de la Grossa?

Aquesta és la típica pregunta d’un contrast de bondat d’ajust.

  • En principi, estam disposats a acceptar que la nostra mostra prové d’una variable amb una certa distribució:

    Que totes les darreres xifres de la Grossa tenen la mateixa probabilitat d’aparèixer.

  • Rebutjarem aquesta hipòtesi si seria molt estrany haver obtingut la nostra mostra si la variable seguís aquesta distribució:

    Si les darreres xifres tenguessin totes la mateixa probabilitat d’aparèixer, esperaríem una certa uniformitat en les seves freqüències a la nostra mostra. Si les freqüències reals se separen molt d’aquesta uniformitat, tendrem motius per dubtar que tots els dígits tenien la mateixa probabilitat d’aparèixer.

Més en general, un contrast de bondat d’ajust té la forma següent: \[ \left\{ \begin{array}{l} H_0: \text{la mostra prové de la distribució esperada}\\ H_1: \text{la mostra NO prové de la distribució esperada} \end{array} \right. \] I, com sempre,

  • Si obtenim evidència que ens permeti rebutjar la hipòtesi nul·la, conclourem que la mostra no prové de la distribució esperada

  • Si no obtenim evidència que ens permeti rebutjar la hipòtesi nul·la, donarem per bo que la mostra prové de la distribució esperada

Els tests de bondat d’ajust es basen bàsicament en:

  1. Comparar les freqüències observades amb les freqüències teòriques de la distribució que contrastam

  2. Determinar si les freqüències observades són prou diferents de les teòriques com per fer inversemblant la hipòtesi nul·la

Per exemple, continuant amb el nostre exemple de la Grossa, si totes les darreres xifres tinguessin la mateixa probabilitat d’aparèixer, la seva probabilitat seria 1/10, i per tant el valor esperat d’aparicions de cadascuna d’elles seria 211/10=21.1.

Les desviacions observades respecte d’aquest valor, representat per la línia vermella al gràfic anterior, són prou grans com per dubtar que aquestes xifres són equiprobables, o podem acceptar que es deuen a l’atzar?

6.2 Test \(\chi^2\) de Pearson

6.2.1 El test bàsic

El test més popular per a variables aleatòries qualitatives, ordinals i quantitatives discretes (o contínues agrupades) és el test \(\chi^2\) de Pearson.

Suposem que tenim una mostra aleatòria simple de mida \(n\). Agrupam tots els resultats possibles en \(k\) classes, \(C_1,\ldots,C_k\) (poden ser tots els resultats individuals possibles, si només n’hi ha un conjunt finit). Així, al nostre exemple de les terminacions de la loteria, hi tenim 10 classes, definides pels 10 valors possibles de la darrera xifra: \(C_1=\{0\}\), \(C_2=\{1\}\),…, \(C_{10}=\{9\}\).

Volem contrastar si les observacions segueixen una distribució totalment coneguda, per a la qual coneixem la probabilitat que una observació caigui dins cada una de les classes.

El contrast és, recordem, \[ \left\{\begin{array}{l} H_{0}: \text{La població té aquesta distribució }\\ H_{1}: \text{La població no té aquesta distribució} \end{array} \right. \]

Per a cada classe \(C_i\), diguem

  • \(obs_i\): la freqüència absoluta observada d’aquesta classe. Per exemple, en el nostre exemple de la loteria, \(obs_1\) és el nombre de vegades que la Grossa ha acabat en 0: 21; \(obs_2\) és el nombre de vegades que la Grossa ha acabat en 1: 8; etc.

  • \(p_i\): la probabilitat teòrica que una observació pertanyi a aquesta classe si \(H_0\) és vertadera. En el nostre exemple de la loteria, totes les probabilitats \(p_i\) valen el mateix, 1/10.

  • \(esp_i\): la freqüència absoluta esperada, o teòrica, d’aquesta classe si \(H_0\) és vertadera: \(esp_i=p_i\cdot n\). En el nostre exemple de la loteria, totes aquestes freqüències absolutes \(esp_i\) valen el mateix, 211/10=21.1.

Rebutjarem \(H_0\) si les \(obs_i\) són prou diferents de les \(esp_i\). Per mesurar-ho, empram el teorema següent:

Teorema 6.1 Si \(H_0\) és vertadera, la distribució contrastada està completament determinada i se satisfan les condicions següents:

  • La mida \(n\) de la mostra és gran (ho fixarem en \(n\geqslant 30\))

  • Les classes cobreixen tots els resultats possibles: \(\sum\limits_{i=1}^kesp_i= n\)

  • Totes les classes tenen prou probabilitat teòrica com per tenir-les en compte (nosaltres emprarem la Regla de Cochran: totes les \(esp_i\geqslant 5\))

aleshores l’estadístic de contrast \[ \chi^2=\sum_{i=1}^k \frac{(obs_{i}-esp_{i})^2}{esp_{i}} \] té (aproximadament) una distribució \(\chi_{k-1}^2\).

Si la distribució contrastada no estava completament determinada i se n’han estimat \(m\) paràmetres a partir de la mostra, s’ha de prendre com a distribució de l’estadístic de contrast una \(\chi_{k-1-m}^2\).

Fixau-vos que l’extadístic de contrast \(\chi^2\) s’obté de la manera següent:

  1. Per cada classe \(C_i\):

    1. Se resta \(obs_{i}-esp_{i}\)
    2. S’eleva el resultat al quadrat: \((obs_{i}-esp_{i})^2\)
    3. Es divideix el resultat per \(esp_i\): \((obs_{i}-esp_{i})^2/esp_{i}\)
  2. Se sumen aquests valors per a totes les classes \(C_i\).

Aleshores, el contrast se fa de la manera següent:

  1. Se calcula el valor que pren l’estadístic de contrast sobre la nostra mostra, diguem-li \(\chi_0^2\)

  2. El p-valor del contrast és \(P(\chi^2\geqslant\chi_0^2)\), calculada tenint en compte que \(\chi^2\) té distribució \(\chi^2\) amb nombre de graus de llibertat:

    • Si no s’ha estimat cap paràmetre de la distribució amb la mostra, el nombre de classes menys 1: \(\chi^2_{k-1}\)

    • Si s’ha estimat qualque paràmetre de la distribució amb la mostra, el nombre de classes menys 1 i menys el nombre de paràmetres estimats: \(\chi^2_{k-1-m}\) amb \(m\) el nombre de paràmetres estimats

Com més gran sigui el nombre de classes emprades, més potència té el contrast, per tant procurau emprar-ne el màxim nombre possible tot respectant la regla de Cochran.

Exemple 6.2 Tornem a l’Exemple 6.1. Recordem que la nostra mostra era

\[ \begin{array}{r|cccccccccc} \hline \text{xifra} & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 \\ \hline \text{freqüència} & 22 & 8 & 13 & 21 & 27 & 31 & 27 & 22 & 23 & 17\\ \hline \end{array} \] Com que totes les feqüències esperades valen 21.1, el valor de l’estadístic de contrast sobre aquesta mostra és \[ \begin{array}{rl} \chi_0^2& \displaystyle=\frac{(22-21.1)^2}{21.1}+\frac{(8-21.1)^2}{21.1}+\cdots\\[2ex] & \qquad\quad\displaystyle +\frac{(23-21.1)^2}{21.1} +\frac{(17-21.1)^2}{21.1}=20.23 \end{array} \]

Calculat amb R:

## [1] 20.23223

En resum, per resoldre la qüestió plantejada en l’Exemple 6.1 procedim de la manera següent:

Contrast: \[ \left\{\begin{array}{l} H_{0}: \text{Les darreres xifres de la Grossa són equiprobables}\\ H_{1}: \text{Les darreres xifres de la Grossa NO són equiprobables} \end{array} \right. \]

Estadístic de contrast: \[ \chi^2=\sum_{i=1}^k \frac{(obs_{i}-esp_{i})^2}{esp_{i}} \]

Estam en les condicions del Teorema 6.1:

  • \(n=211\), gran

  • Les classes, que corresponen als dígits 0,…,9, cobreixen tots els resultats possibles; \(k=10\)

  • Cada \(esp_i=21.1\geqslant 5\)

Per tant, com que no estimam cap paràmetre de la distribució, aquest estadístic de contrast segueix una llei \(\chi^2_9\).

Valor de l’estadístic de contrast: Ja l’hem calculat, 20.23.

p-valor: \(P(\chi_9^2\geqslant 20.23)=\texttt{1-pchisq(20.23,9)}= 0.0165\)

Conclusió: Hem obtingut evidència estadísticament significativa que les darreres xifres de la Grossa no apareixen de manera equiprobable (test \(\chi^2\) de Pearson, p-valor 0.0165).

Exemple 6.3 Anem a comprovar amb una simulació aquest p-valor. El que farem serà el següent: Repetirem 10000 vegades, amb un replicate, el procés d’extreure a l’atzar i de manera equiprobable 211 dígits, calcular-ne la taula de freqüències de cada una d’aquestes mostres i calcular l’estadístic de contrast \(\chi^2\) sobre cada una d’aquestes taules. Finalment, calcularem la fracció de vegades que \(\chi^2\) ha donat per sobre del valor de la nostra mostra, que era 20.23.

## [1] 0.0158

Només un 1.6% de les 10000 simulacions han donat un valor de \(\chi^2\) més gran o igual que el de la mostra de Grosses de Nadal. Per tant, el valor de la \(\chi^2\) de la mostra de Grosses de Nadal és anormalment gran si les darreres xifres tenen totes la mateixa probabilitat d’aparèixer. Naturalment, anormal no vol dir impossible, però…

La funció per realitzar un test \(\chi^2\) amb R és

on obs és el vector de les freqüències observades de les classes, \((obs_i)_i\) i p s’ha d’igualar al vector de les probabilitats teòriques de les classes, \((p_i)_i\).

Si no s’especifica aquest paràmetre p, R entén que totes les probabilitats són iguals.
La suma de les probabilitats entrades a p ha de ser igual a 1. Recordau que les classes han de cobrir tots els casos possibles.

Al nostre exemple de la Grossa de Nadal, entraríem simplement

## 
##  Chi-squared test for given probabilities
## 
## data:  Grossa
## X-squared = 20.232, df = 9, p-value = 0.01653

Obtenim el valor de \(\chi^2\), X-squared, els graus de llibertat, df, i el p-valor, p-value.

Exemple 6.4 Un tècnic de medi ambient vol estudiar l’augment de temperatura de l’aigua a dos quilòmetres dels abocaments d’aigua autoritzats d’una planta industrial.

El responsable de l’empresa afirma que

aquests augments de temperatura segueixen una llei normal amb \(\mu=3.5\) dècimes de grau C i \(\sigma=0.7\) dècimes de grau C.

El tècnic ho posa en dubte. Per decidir-ho, pren una mostra aleatòria d’observacions de l’augment de les temperatures (en dècimes de grau). Les dades són les de la taula següent, ja agrupades en classes:

\[ \begin{array}{c|c} \text{Rang d'augments} & \text{Freqüències}\\ \hline 1.45\text{ a }1.95 & 2 \\ 1.95\text{ a }2.45 & 1 \\ 2.45\text{ a }2.95 & 4 \\ 2.95\text{ a }3.45 & 15 \\ 3.45\text{ a }3.95 & 10 \\ 3.95\text{ a }4.45 & 5 \\ 4.45\text{ a }4.95 & 3 \\ \hline \text{Total} & 40\\ \end{array} \]

Hi ha evidència que la sospita del tècnic sigui vertadera, amb un nivell de significació del 5%?

Contrast: \[ \left\{ \begin{array}{l} H_{0}:\text{La distribució dels augments de temp. és $N(3.5,0.7)$}\\ H_{1}: \text{La distribució dels augments de temp. NO és $N(3.5,0.7)$} \end{array} \right. \]

Estadístic de contrast: Volem emprar el test \(\chi^2\), així que l’estadístic de contrast serà \[ \chi^2=\sum_{i=1}^k \frac{(obs_{i}-esp_{i})^2}{esp_{i}} \]

Cal comprovar que se satisfan les condicions del Teorema 6.1.

  • La mida de la mostra és \(n=40\), suficient.

  • Les classes, cobreixen tots els resultats possibles?

No, perquè els possible valors d’una variable normal són tots els nombres reals, i les nostres classes només cobreixen l’interval [1.45,4.95]. El que farem serà afegir les cues a les dues classes dels extrems, i considerar les classes \[ \begin{array}{l|c} \text{Rang d'augments} & \text{Freqüències}\\ \hline \text{menys de 1.95} & 2 \\ 1.95 \text{ a } 2.45 & 1 \\ 2.45\text{ a }2.95 & 4 \\ 2.95\text{ a }3.45 & 15 \\ 3.45\text{ a }3.95 & 10 \\ 3.95\text{ a }4.45 & 5 \\ \text{més de 4.45} & 3 \\ \hline \text{Total} & 40\\ \end{array} \] Ara ja cobreixen tots els resultats possibles.

  • Totes les freqüències esperades són com a mínim 5?

Doncs no ho sabem, les haurem de calcular. Recordem que són les freqüències esperades de cada classe suposant que la variable poblacional és \(N(3.5,0.7)\), i s’obtenen multiplicant la probabilitat de cada interval (amb aquesta distribució) per la mida de la mostra, 40.

La probabilitat de la 1a classe és \[ p_1 =P(X\leqslant 1.95)=\texttt{pnorm(1.95,3.5,0.7)}=0.0134 \] i la seva freqüència esperada és \[ esp_1=p_1\cdot n= 0.0134\cdot 40=0.536 \]

La probabilitat de la 2a classe és \[ \begin{array}{rl} p_2 & =P(1.95\leqslant X\leqslant 2.45)\\ & =\texttt{pnorm(2.45,3.5,0.7)-pnorm(1.95,3.5,0.7)}=0.0534 \end{array} \] i per tant la seva freqüència esperada és \[ esp_2=p_2\cdot n= 0.0534\cdot 40= 2.136 \]

Anam a calcular amb R totes les freqüències esperades d’aquesta manera i d’un sol cop:

## [1]  0.536  2.136  5.968 10.220 10.733  6.912  3.495

Per tant, tenim la taula \[ \begin{array}{r|c|c} \text{Rang de temperatures} & obs_i & esp_i\\ \hline \text{menys de 1.95} & 2 & 0.536 \\ 1.95\text{ a }2.45 & 1 & 2.136 \\ 2.45\text{ a }2.95 & 4 & 5.968 \\ 2.95\text{ a }3.45 & 15 & 10.220 \\ 3.45\text{ a }3.95 & 10 & 10.733 \\ 3.95\text{ a }4.45 & 5 & 6.912 \\ \text{més de 4.45} & 3 & 3.495 \end{array} \] Tenim freqüències esperades per davall de 5, per tant no podem aplicar el test \(\chi^2\) tal qual. El que farem serà agrupar classes contigües a fi d’obtenir freqüències esperades més grans o iguals que 5 amb el màxim de classes. En concret, haurem d’agrupar per un costat les tres primeres classes, i per un altre costat les dues darreres classes \[ \begin{array}{r|c|c} \text{Rang de temperatures} & obs_i & esp_i\\\hline \text{menys de 2.95} & 7 & 8.64\\ 2.95\text{ a }3.45 & 15 & 10.22 \\ 3.45\text{ a }3.95 & 10 & 10.733 \\ \text{més de 3.95} & 8& 10.407 \end{array} \]

Ara ja podem aplicar el text \(\chi^2\) amb aquestes \(k=4\) classes. Com que no hem estimat cap paràmetre, prendrem com a distribució de l’estadístic de contrast una \(\chi^2_3\).

Valor de l’estadístic de contrast: \[ \begin{array}{rl} \chi_0^2 &\displaystyle=\frac{(7-8.64)^2}{8.64}+\frac{(15-10.22)^2}{10.22}\\ &\displaystyle\qquad\quad +\frac{(10-10.733)^2}{10.733}+\frac{(8-10.407)^2}{10.407}=3.154 \end{array} \]

p-valor: \[ P(\chi_{3}^2\geqslant 3.154)=\texttt{1-pchisq(3.154,3)}=0.37 \]

Conclusió: Podem acceptar que els augments de temperatures observats segueixen una llei normal \(N(3.5,0.7)\) (test \(\chi^2\) de Pearson, p-valor=0.37).

Amb R, si aplicam directament la funció chisq-test als vectors de freqüències i probabilitats de les classes originals (bé, de les esteses per cobrir tot \(\mathbb{R}\)) R ens avisa que el resultat no és de fiar:

## Warning in chisq.test(freq.obs, p = prob.esp): Chi-squared approximation may be
## incorrect
## 
##  Chi-squared test for given probabilities
## 
## data:  freq.obs
## X-squared = 8.1337, df = 6, p-value = 0.2285

Haurem d’agrupar primer a mà les classes, abans d’aplicar chisq.test:

## 
##  Chi-squared test for given probabilities
## 
## data:  freq.obs.agrup
## X-squared = 3.1531, df = 3, p-value = 0.3686

6.2.2 Mètode de Montecarlo

Quan volem realitzar un test \(\chi^2\) i no es compleixen les condicions per que el p-valor tengui sentit, es pot emprar el mètode de Montecarlo (o basat en simulacions) per estimar el p-valor:

  • Es genera un conjunt molt gran de mostres aleatòries amb la distribució esperada i de la mateixa mida que la nostra mostra

  • Es calcula l’estadístic de contrast \(\chi^2\) sobre cada mostra

  • S’estima el p-valor mitjançant la fracció de mostres que han donat un \(\chi^2\) més gran que el de la nostra mostra

Recordau que ja hem emprat aquest mètode a l’Exemple 6.3 per “comprovar” el p-valor obtingut a l’Exemple 6.2 amb el test \(\chi^2\) “exacte”. Però no cal complicar-se tant la vida per realitzar-lo. La funció chisq.test també permet usar el mètode de Montecarlo, entrant-hi el paràmetre simulate.p.value=TRUE i, si es vol, igualant el paràmetre B al nombre de simulacions desitjat (per defecte, en farà 2000). Naturalment, com que es basa en una simulació aleatòria, pot ser que en cada execució doni un p-valor diferent.

Per exemple, podíem emprar el mètode de Montecarlo per realitzar el contrast de l’Exemple 6.4 amb les classes originals esteses, sense agrupar-les.

## 
##  Chi-squared test for given probabilities with simulated p-value (based
##  on 5000 replicates)
## 
## data:  freq.obs
## X-squared = 8.1337, df = NA, p-value = 0.2206

Conclusió: Podem acceptar que els augments de temperatures observats segueixen una llei normal \(N(3.5,0.7)\) (test \(\chi^2\) de Montercarlo, p-valor=0.22).

Arribam, doncs, a la mateixa conclusió que amb el test \(\chi^2\) aplicat a les dades agrupades.

6.2.3 Test \(\chi^2\) amb paràmetres poblacionals desconeguts

De vegades ens interessarà contrastar si les observacions segueixen algun tipus determinat de distribució (Poisson, normal…) amb algun paràmetre indeterminat. En aquest cas, primer estimam els paràmetres desconeguts a partir de les observacions i a continuació contrastam l’ajustament de la mostra a la distribució amb aquests paràmetres.

El test és exactament el mateix, excepte que ara cal restar al nombre \(k-1\) de graus de llibertat de la \(\chi^2\) el nombre \(m\) de paràmetres que hem estimat.

Exemple 6.5 Volem determinar si el nombre de vegades que apareix la seqüència GATACA en una cadena d’ADN humà de longitud 1000 segueix una llei de Poisson.

Per fer-ho, prenem 576 mostres de cadenes d’ADN humà de longitud 1000 i hi comptam els nombres de GATACA. Els resultats són els de la taula següent: \[ \begin{array}{r|rrrrrr} \hline \text{nombre $x_i$ de vegades} & 0 & 1 & 2 & 3 & 4 & 5 \\[-1ex] \text{que hi apareix GATACA} & & & & & & \\ \hline \text{freqüència $obs_{i}$} & 229 & 211 & 93 & 35 & 7 & 1 \\ \hline \end{array} \]

Contrast: \[ \left\{ \begin{array}{ll} H_0: \text{La mostra prové d'una distribució $Po(\lambda)$}\\ H_1: \text{La mostra no prové d'aquesta distribució} \end{array} \right. \]

Necessitam estimar el paràmetre \(\lambda\), que és el valor esperat d’una variable Poisson i per tant l’estimarem amb la mitjana mostral: \[ \lambda =\dfrac{229 + 211 + 93 + 35 +7 + 1 }{576}=\dfrac{576}{576}=1 \]

Aquesta no és la mitjana de la mostra!

Ara sí: \[ \begin{array}{rl} \lambda & =\dfrac{229\cdot 0+ 211\cdot 1+ 93\cdot 2+ 35\cdot 3+7\cdot 4+ 1\cdot 5}{229+ 211+ 93 +35 + 7 + 1}\\ & =\dfrac{535}{576}=0.929 \end{array} \]

Estadístic de contrast: Volem emprar el test \(\chi^2\), així que l’estadístic de contrast serà \[ \chi^2=\sum_{i=1}^k \frac{(obs_{i}-esp_{i})^2}{esp_{i}} \]

Cal comprovar que se satisfan les condicions del Teorema 6.1.

  • La mida de la mostra és \(n=576\), ben gran.

  • Per veure si les freqüències esperades són totes més grans o iguals que 5, les hem de calcular.

Recordem que, si \(X\) és una variable aleatòria de Poisson \(Po(\lambda)\), \[ P(X=k)=e^{-\lambda}\cdot \frac{\lambda^k}{k!} \]

Les probabilitats \(p_i\) d’obtenir \(i=0,1,2,3,4,5\) són

## [1] 0.395 0.367 0.170 0.053 0.012 0.002

i les freqüències esperades \(esp_i=p_i\cdot n\) d’aquests valors són

## [1] 227.49 211.34  98.17  30.40   7.06   1.31

Obtenim la taula \[ \begin{array}{l|rrrrrr} \hline x_i & 0 & 1 & 2 & 3 & 4 & 5 \\ \hline obs_{i} & 229 & 211 & 93 & 35 & 7 & 1 \\ \hline p_i &0.395 & 0.367 &0.170 & 0.053 & 0.012 & 0.002 \\ \hline esp_i & 227.49 & 211.34 & 98.17 & 30.40 & 7.06 & 1.31\\ \hline \end{array} \]

Aquestes classes estan malament! Recordau que les classes han de cobrir tots els resultats possibles! Els resultats possibles d’una variable Poisson són tots els nombres naturals.

Per cobrir tots els resultats possibles, cal canviar el 5 per “5 o més” i recalcular la probabilitat i la freqüència esperada d’aquesta classe: la seva probabilitat és

## [1] 0.003

i la seva freqüència esperada és

## [1] 1.55

Obtenim la taula \[ \begin{array}{l|rrrrrr} \hline x_i & 0 & 1 & 2 & 3 & 4 & \geqslant 5 \\ \hline obs_{i}& 229 & 211 & 93 & 35 & 7 & 1 \\ \hline p_i &0.395 & 0.367 &0.170 & 0.053 & 0.012 & 0.003 \\ \hline esp_i & 227.49 & 211.34 & 98.17 & 30.40 & 7.06 & 1.55\\ \hline \end{array} \]

Per que totes les classes tenguin freqüència esperada més gran o igual que 5, agruparem les dues darreres classes

\[ \begin{array}{l|rrrrr} \hline x_i & 0 & 1 & 2 & 3 & \geqslant 4 \\ \hline obs_{i}& 229 & 211 & 93 & 35 & 8 \\ \hline p_i &0.395 & 0.367 &0.170 & 0.053 & 0.015 \\ \hline esp_i & 227.49 & 211.34 & 98.17 & 30.40 & 8.61\\ \hline \end{array} \]

Ara totes les freqüències esperades són prou grans i les classes cobreixen tots els valors possibles de la variable. L’estadístic de contrast té distribució \(\chi_{k-1}^2\) amb \(k=5\), per tant \(\chi_{4}^2\).

Com que hem estimat la \(\lambda\), hem de considerar l’estadístic de contrast amb distribució \(\chi_{k-1-m}^2\) amb \(k\) el nombre de classes i \(m\) el nombre de paràmetres estimats.

Perdó. L’estadístic de contrast té distribució \(\chi_{3}^2\). Podem continuar.

Valor de l’estadístic de contrast:

## [1] 1.02

p-valor: \(P(\chi_3^2\geqslant 1.02)=\texttt{1-pchisq(1.02,3)}=0.796\)

Conclusió: Podem acceptar que les observacions trobades segueixen una llei de Poisson (test \(\chi^2\) de Pearson, p-valor 0.796).

Amb R simplement ens estalviam el càlcul explícit de \(\chi_0^2\) i del p-valor, perquè tota la saragata de calcular freqüències esperades i agrupar a fi que se satisfaci la llei de Cochran no ens la podem estalviar.

## 
##  Chi-squared test for given probabilities
## 
## data:  obs.gataca
## X-squared = 1.0215, df = 4, p-value = 0.9065

Conclusió: No podem rebutjar que les observacions trobades segueixin una llei de Poisson (test \(\chi^2\) de Pearson, p-valor 0.9065????).

No havíem quedat que el nombre de graus de llibertat era 3? R n’ha emprat 4, df=4, perquè no sap que hem estimat un paràmetre

R ha calculat el p-valor prenent \(\chi_4^2\), nosaltres l’hem de calcular amb \(\chi_3^2\)

## [1] 0.7960498

Ara ja tenim el mateix p-valor que abans.

6.3 Test de Kolmogorov-Smirnov

El test de Kolgomorov-Smirnov (test K-S, per abreviar) és el test més popular per contrastar si una mostra segueix o no una distribució contínua concreta, sense restriccions sobre la mida de la mostra. Es pot emprar amb tota distribució contínua sempre que estigui completament especificada.

Per a distribucions concretes, mides de mostres dins marges concrets o quan estimam els paràmetres, existeixen tests específics millors que el K-S, però no els veurem aquí. Estan a la lliçó corresponent de R i potser que les hàgiu d’emprar als tallers.

Tot seguit explicam com funciona aquest test.

Partim d’una mostra \(x_1,x_2,\ldots,x_n\) amb tots els valors diferents i volem contrastar si ha estat produïda per una variable \(X\) contínua amb funció de distribució \(F_X\).

Per exemple, suposem que volem constrastar si podem acceptar que els valors del vector següent provenen d’una distribució normal \(N(3,1.5)\).

Per tant, el contrast és \[ \left\{ \begin{array}{l} H_0: \text{ aquesta mostra prové d'una $X\sim N(3,1.5)$}\\ H_0: \text{ aquesta mostra no prové d'una $X\sim N(3,1.5)$} \end{array} \right. \]

  1. Ordenam la mostra: \(x_{(1)}< x_{(2)}<\cdots< x_{(n)}\).
## [1] 1.34 1.54 2.25 3.58 4.57 5.84

  1. Consideram la funció de distribució mostral \(F_{n}\) d’aquesta mostra: la funció de distribució d’una variable \(n\) amb domini tot \(\mathbb{R}\) per a la qual cada \(x_{(i)}\) tingués probabilitat \(1/n\) i la resta de reals tingués probabilitat 0: \[ F_n(x)=\left\{\begin{array}{ll} 0 &\text{ si } x< x_{(1)} \\ \dfrac{k}{n}&\text{ si } x_{(k)}\leqslant x < x_{(k+1)}\\ 1 & \text{ si } x_{(n)} \leqslant x \end{array} \right. \]

    A la nostra mostra ja ordenada, seria \[ F_6(x)=\left\{\begin{array}{ll} 0 &\text{ si } x< 1.34 \\ 1/6 &\text{ si } 1.34\leqslant x <1.54\\ 2/6 &\text{ si } 1.54\leqslant x <2.25\\ 3/6 &\text{ si } 2.25\leqslant x <3.58\\ 4/6 &\text{ si } 3.58\leqslant x <4.57\\ 5/6 &\text{ si } 4.57\leqslant x <5.84\\ 1 &\text{ si } 5.84\leqslant x \end{array} \right. \] El gràfic d’aquesta distribució és el següent:

  1. Comparam \(F_n(x)\) amb \(F_X(x)\). Si són molt diferents, rebutjarem que la mostra s’ha produït amb distribució \(F_X\).

Per comparar \(F_n(x)\) amb \(F_X(x)\), calcularem \[ \max\{|F_n(x)-F_X(x)|\mid x\in \mathbb{R}\} \] Com que \(F_X\) és creixent, aquest màxim s’assoleix a qualque escaló. Per tant, per calcular aquest màxim primer calculam, per a cada \(x_{(i)}\), la seva discrepància: la distància entre \(F_X(x_{(i)})\) i els extrems de l’escaló corresponent, \[ D_n(x_{(i)})=\max\Big\{\Big| F_X(x_{(i)})-\frac{i-1}{n}\Big|, \Big|F_{X}(x_{(i)})-\frac{i}{n}\Big|\Big\} \] i prenem la discrepància màxima, és a dir, el màxim d’aquestes discrepàncies: \[ D_n=\max\big\{D_n(x_{(i)})\mid i=1,\ldots, n\big\} \]

Al nostre exemple, els valors \(F_X(x_{(i)})\) són:

## [1] 0.134 0.165 0.309 0.650 0.852 0.971

i per tant les discrepàncies són \[ \begin{array}{rl} D_6(x_{(1)}) & =\max\{| F_X(1.34)-0|, |F_X(1.34)-1/6|\}\\ & =\max\{|0.134-0|, |0.134-1/6|\}\\ & =\max\{ 0.134, 0.033\}=0.134\\[2ex] D_6(x_{(2)}) & =\max\{| F_X(1.54)-1/6|, |F_X(1.54)-2/6|\}\\ & =\max\{| 0.165-1/6|, |0.165-2/6|\}\\ & =\max\{ 0.002, 0.168\}=0.168\\[2ex] D_6(x_{(3)}) & =\max\{| F_X(2.25)-2/6|, |F_X(2.25)-3/6|\}\\ & =\max\{|0.309-2/6|, |0.309-3/6|\}\\ & =\max\{ 0.024, 0.191\}=0.191\\[2ex] D_6(x_{(4)}) & =\max\{| F_X(3.58)-3/6|, |F_X(3.58)-4/6|\}\\ & =\max\{| 0.65-3/6|, |0.65-4/6|\}\\ & =\max\{ 0.15, 0.017\}=0.15\\[2ex] D_6(x_{(5)}) & =\max\{| F_X(4.57)-4/6|, |F_X(4.57)-5/6|\}\\ & =\max\{|0.852-4/6|, |0.852-5/6|\}\\ & =\max\{ 0.185, 0.019\}=0.185\\[2ex] D_6(x_{(6)}) & =\max\{| F_X(5.84)-5/6|, |F_X(5.84)-6/6|\}\\ & =\max\{| 0.971-5/6|, |0.971-1|\}\\ & =\max\{ 0.138, 0.029\}=0.138 \end{array} \] i la discrepància màxima és \[ D_6=\max\{0.134,0.168,0.191,0.15,0.185,0.138\}=0.191 \]

  1. Aquesta discrepància màxima segueix una distribució coneguda (que no depèn de la \(X\) mentre sigui contínua) que permet calcular el p-valor com la probabilitat que \(D_n\) prengui un valor més gran o igual que l’obtingut sobre la nostra mostra. La corresponent funció de distribució acumulada és donada per la funció pkolm del paquet kolmim. Per tant, el p-valor és
## [1] 0.951

Conclusió: Podem acceptar que la mostra prové d’una variable aleatòria normal amb \(\mu=3\) i \(\sigma=1.5\) (test de Kolmogorov-Smirnov, p-valor 0.95)

Per realitzar un test K-S amb R, disposam de la instrucció

on x és el vector que conté la mostra, la pdistribució és la funció de distribució que contrastam, i els paràmetres són els de la distribució. Al nostre exemple, simplement entraríem

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  valors
## D = 0.19146, p-value = 0.95
## alternative hypothesis: two-sided

Dóna el valor de l’estadístic (la D) i el p-valor

6.4 Test de Kolmogorov-Smirnov-Lilliefors

Quan es vol emprar el test K-S per contrastar si una mostra prové de qualque distribució normal, i no d’una concreta, l’opció més popular és:

  • Estimar els paràmetres de la normal a partir de la mostra

  • Calcular l’estadístic del test K-S amb aquests paràmetres

  • Emprar la distribució del test K-S-Lilliefors per calcular el p-valor

Amb R, aquest test està implementat en la funció lillie.test del paquet nortest.

Per exemple, si ens haguéssim demanat d’entrada si la nostra mostra prové de qualque variable normal, haguéssim entrat

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  valors
## D = 0.19912, p-value = 0.6425
L’estadístic de contrast no val el mateix que abans, perquè ara la variable normal que contrasta té paràmetres \(\mu=\texttt{mean(valors)}=3.187\) i \(\sigma=\texttt{sd(valors)}=1.795\).