Analiza danych jakościowych A

Literatura

• Alan Agresti, Categorical Data Analysis, https://mathdept.iut.ac.ir/sites/mathdept.iut.ac.ir/files/AGRESTI.PDF
  
• Laura A. Thompson R (and S-PLUS) Manual to Accompany Agresti’s Categorical Data Analysis (2002) https://home.comcast.net/~lthompson221/Splusdiscrete2.pdf

Typy zmiennych

• jakościowe
  • nominalne (czynnikowe (factor),categorical)
  • porządkowe (czynnikowe porządkowe (ordered factor))
• ilościowe

Rozkłady prawdopodobieństwa zmiennych dyskretnych

Bernouillego

\[\begin{array}{rll} P(Y=1)&=&\pi\\ P(Y=0)&=&1-\pi \end{array}\]

\[\begin{array}{rll} E(Y)&=&1*\pi+0*(1-\pi)=\pi\\ V(Y)&=&E(Y^2)-(E(Y))^2=\pi-\pi^2=\pi(1-\pi) \end{array}\]

Rozkład dwumianowy

\(Y_i\) są niezależne i mają rozkład Bernoulliego
\[\begin{array}{rll} Y&=&\sum_{i=1}^{n}Y_i\\ P(Y=k)&=&\left( \begin{array}{c} n \\ k \end{array} \right)\pi^{k}(1-\pi)^{n-k}\\ \mu&=&E(Y)=\sum_{i=1}^{n}E(Y_i)=n\pi\\ \sigma^2&=&V(Y)=\sum_{i=1}^{n}V(Y_i)=n\pi(1-\pi) \end{array}\]

Rozkład wielomianowy

Uogólnienie rozkładu dwumianowego \[\begin{array}{rll} P(Y_i=1)&=&\pi_1\\ P(Y_i=2)&=&\pi_2\\ &\dots&\\ P(Y_i=c)&=&\pi_c \end{array}\].. \(Y_i\) \((i=1,2,\dots,n)\) są niezależne \[Y=(n_1,n_2,\dots,n_c)\iff \#(i:Y_i=1)=n_1,\#(i:Y_i=2)=n_2,\dots,\#(i:Y_i=c)=n_c\]
\[P(Y=(n_1,n_2,\dots,n_c))=\frac{n!}{n_1!n_2!\dots n_c!}\pi_1^{n_1}\pi_2^{n_2}\dots \pi_c^{n_c}\]

Rozkład Poissona

Prawo małych liczb

Wyprowadzenie
Związek z rozkładem wielomianowym

Metoda największej wiarygodności

Funkcja wiarygodności

Rzucono 10 razy monetą. Wyoadło 7 orłóW.Jakie może być prawdopodobieństwo wypadnięcia orła?

Gdyby prawdopodobieństwo wypadnięcia orła było równe \(\pi\), to zaobserwowany wynik pojawiłby się z prawdopodobieństwem:
\[\left( \begin{array}{c} 10 \\ 7 \end{array} \right)\pi^{7}(1-\pi)^{3}\]..

pi=0.4 0.0425 
pi=0.5 0.1172 
pi=0.6 0.215 
pi=0.7 0.2668 
pi=0.8 0.2013 

Metoda największej wiarygodności

Niech \(L(\beta)\) - logarytm funkcji wiarygodności. Estymator największej wiarygodności: \[\widehat{\beta}: \frac{\partial}{\partial\beta}L(\beta)=0\] Błąd standardowy estymatora NW parametru \(\beta\) jest pierwiastkiem kwadratowym przekątnej macierzy \(\Im^{-1}\), obliczonym w punkcie \(\widehat{\beta}\), gdzie \[\Im=E\left[-\frac{\partial^2L(\beta)}{\partial\beta_j\partial\beta_k}\right]\]
SE dla rozkładu dwumianowego na dwa sposoby

Jednowymiarowe, asymptotyczne testy istotności

\[\begin{array}{rcl} H_0&:&\beta=\beta_0\\ H_1&:&\beta\neq\beta_0 \end{array}\]

Test Walda

\[z=\frac{\widehat{\beta}-\beta_0}{SE(\widehat{\beta})}\asymp\mathcal{N}(0,1)\]

Test ilorazu wiarygodności

\[\begin{array}{rcl} l_0&=&max\{\beta\in H_0:l(\beta)\}\\ l_1&=&max\{\beta\in H_0\cup H_1:l(\beta)\} \end{array}\]

Tw (Wilks 1935)
\[-2 \log\left(\frac{l_0}{l_1}\right)\asymp\chi^2_k\] \[k=\dim(H_0\cup H_1)-\dim(H_0)\]

Test punktowy (Fisher,Rao)

\[\begin{array}{rcl} u(\beta_0) &=& \left. \frac{\partial L(\beta)}{\partial \beta}\right|_{\beta=\beta_0}\\ i(\beta_0) &=& -\left.E\left[\frac{\partial^2L(\beta)}{\partial\beta^2}\right] \right|_{\beta=\beta_0} \end{array}\]

Statystyka punktowa

\[ \begin{array}{rcl} u(\beta_0)\sqrt{i(\beta_0)} & \asymp & \mathcal{N}(0,1)\\ u(\beta_0)^2i(\beta_0) & \asymp & \chi^2_1 \end{array} \]

Przykłady dla rozkładu dwumianowego

Przedziały ufności

Przedziały ufności (dwustronne) dla \(\beta_0\) na poziomie ufności \(1-\alpha\) można otrzymać, gdy znana jest statystyka testowa \(T\) do testowania na poziomie istotności \(\alpha\) hipotez \[\begin{array}{rcl} H_0&:&\beta=\beta_0\\ H_1&:&\beta\neq\beta_0 \end{array}\]

Wystarczy rozwiązać wzgledem \(\beta_0\) równanie \[ P_{H_0} \left( \left| T \right| < u(\alpha) \right) = 1-\alpha \]

Uwaga. Statystyka \(T\) jest funkcją \(\beta_0\)

W 50 rzutach kostką 20 razy wypadło więcej niż 4 oczka. Skonstruuj przedział ufności na poziomie 0.95 dla prawdopodobieństwa wypadnięcia 5 lub 6 oczek. Czy można uważać, że kostka jest rzetelna?

Test Walda

Testujemy hipotezę, że \(\pi\ne \frac{1}{3}\)

pihat <- 20/50
wald <- function(pi.0) {
  (se <- sqrt(pihat*(1-pihat)/50))
  (T <- (pihat-pi.0)/se)
  pv <- 2*(1-pnorm(T))
  list(se=se,T=T,pv=pv)
}

wald(1/3)

$se
[1] 0.06928203

$T
[1] 0.9622504

$pv
[1] 0.3359238

wald(0.2)

$se
[1] 0.06928203

$T
[1] 2.886751

$pv
[1] 0.003892417

wald(0.2642)

$se
[1] 0.06928203

$T
[1] 1.960104

$pv
[1] 0.04998362

wald(0.2643)

$se
[1] 0.06928203

$T
[1] 1.958661

$pv
[1] 0.05015253

Przedział Walda

n <- 50
pihat + c(-1, 1) * qnorm(p = 0.975) * sqrt((pihat * (1 - pihat))/n) 

[1] 0.2642097 0.5357903

Przedział z testu ilorazu wiarygodności \[ LR=2 \left[ y\log \left( \frac{\widehat{\pi}}{\pi_0} \right) + (n-y)\log \left( \frac{1-\widehat{\pi}}{1-\pi_0} \right) \right] \]

Trzeba rozwiązać równanie \[ LR=\text{qchisq}(1-\alpha,df=1) \]

LR <- function(pi.0, y, n, alpha) {
  phat <- y/n
   2*(y*log(phat/pi.0) + (n-y)*log((1-phat)/(1-pi.0)))
}
LR0 <- function(pi.0, y, n, alpha){
  LR(pi.0, y, n, alpha)- qchisq(1-alpha,df=1)
} 

LR(0.25,20,50,0.05)

[1] 5.411532

(q <- qchisq(1-0.05,df=1))

[1] 3.841459

LR0(0.25,20,50,0.05)

[1] 1.570073

uniroot(f=LR0, interval=c(0.2,0.3), n=50, y=20, alpha=.05)

$root
[1] 0.271764

$f.root
[1] -0.001969518

$iter
[1] 4

$init.it
[1] NA

$estim.prec
[1] 6.103516e-05

uniroot(f=LR0, interval=c(0.5,0.6), n=50, y=20, alpha=.05)

$root
[1] 0.5382601

$f.root
[1] -0.0008796044

$iter
[1] 4

$init.it
[1] NA

$estim.prec
[1] 6.103516e-05

Przedział; \((0.2248,0.5948)\)
Przedział z testu punktowego
Należy rozwiązać, względem \(\pi_0\) równanie
\[ \left| \frac{\widehat{\pi}-\pi_0}{\sqrt{\frac{\pi_0(1-\pi_0)}{n}}} \right | \]

pt <- function(pi.0,y,n,alpha) {
  pihat <- y/n
  (pihat-pi.0)^2-qchisq(1-alpha,df=1)*(pi.0*(1-pi.0))/n
}

pt(0.4,20,50,0.05)

[1] -0.018439

uniroot(f=pt, interval=c(0.2,0.3), n=50, y=20, alpha=.05)

$root
[1] 0.2760895

$f.root
[1] -1.577236e-06

$iter
[1] 4

$init.it
[1] NA

$estim.prec
[1] 6.103516e-05

uniroot(f=pt, interval=c(0.5,0.6), n=50, y=20, alpha=.05)

$root
[1] 0.5381681

$f.root
[1] -4.957657e-06

$iter
[1] 4

$init.it
[1] NA

$estim.prec
[1] 6.103516e-05

Przedział; \((0.2761,0.5382)\)

Przedział Wilsona

\[ \begin{array}{rcl} l & = &\frac{z(\alpha)^2}{n+z(\alpha)^2}\\ p & = & (1-l)\widehat{\pi}+l\frac{1}{2}\\ w & = & (1-l)\widehat{\pi}(1-\widehat{\pi})+l\frac{1}{4}\\ p & \pm & z(\alpha)\sqrt{\frac{w}{n+z(\alpha)^2}} \end{array} \] Przedział` Wilsona 95%

l <- 1.96^2/(50+1.96^2)
p <- (1-l)*0.4+l*0.5
w <- (1-l)*0.4*0.6+l*0.25
round(p+c(-1,1)*sqrt(w/(50+1.96^2)),digits=4)

[1] 0.3403 0.4740

Gdy \(y=0\) przedział Walda nie ma sensu.
Przedział punktowy dla \(y=0\) ma postać \[\left(0,\frac{z(\alpha)}{n+z(\alpha)}\right)\] Przedział największej wiarygodności
Logarytm funkcji wiarygodności wynosi \(l(\pi)=n\log(1-\pi)\), Stąd \(l_1=l(0)=0\), \(l_0=n\log(1-\pi_0)\) \[ LR=-2n\log \left( 1-\pi_0 \right) \] Przedział LR będzie miał postać \[ \left( 0, 1-\exp(-\frac{\text{qchisq}(1-\alpha,df=1)}{2n}) \right) \] Przedział Wilsona zaś będzie równy \[ \begin{array}{rcl} l & = &\frac{z(\alpha)^2}{n+z(\alpha)^2}\\ p & = & l\frac{1}{2}\\ w & = & l\frac{1}{4}\\ p & \pm & z(\alpha)\sqrt{\frac{w}{n+z(\alpha)^2}} \end{array} \]

Dla dla \(n=50\) 95% przedział ufności punktowy ma postać

round(c(0,1.96/(50+1.96)),digits = 4)

[1] 0.0000 0.0377

Przedział LR

round(c(0,1-exp(-qchisq(0.95,1)/100)),digits = 4)

[1] 0.0000 0.0377

Przedział Wilsona

l <- 1.96^2/(50+1.96^2)
p <- l*0.5
w <- l*0.25
round(p+c(-1,1)*sqrt(w/(50+1.96^2)),digits=4)

[1] 0.0175 0.0539