library(ggplot2)
library(xtable)
library(dplyr)
library(plotROC)
library(ROCR)
library(MASS)

Do戼㸳戼㸹czanie pakietu: 㤼㸱MASS㤼㸲

Nast攼㹡puj戼㸹cy obiekt zosta戼㸳 zakryty z 㤼㸱package:dplyr㤼㸲:

    select

Rozkład wielomianowy

Estymatory NW parametrów modelu

Statystyka \(X^2\) Pearsona

\[ \begin{array}{rcl} H_0 & = & \pi_j=\pi_j^0 \: (j=1,2,\dots,I)\\ H_1 & = & \neg H_0\\ \end{array} \]

Niech \(\mu_j=n\pi_j^0\). Statystyka \(X^2\) Pearsona: \[ X^2=\sum_{j=1}^{I}\frac{\left( n_j-\mu_j\right)^2}{\mu_j} \]

Gdy \(n\to\infty\) to \(X^2\to \chi^2_{I-1}\)

Uzasadnienie Fishera

Statystyka Youle’a

\[ \phi^2 = \sum_{j=1}^{I}\pi_j^0 \left( 1-\frac{\widehat{\pi_j}}{\pi_j^0} \right)^2 \]

\[ \phi=\sqrt{\frac{X^2}{n}} \] Dowód

Test \(\chi^2\) dla dużych \(n\) często odrzuca \(H_0\) !!

Przykład

Test ilorazu wiarygodności

\[ H_0:\pi_j=\pi_j^0 \]

\[ \log(l_0)=\sum_j n_j \log \left( \pi_j^0 \right) \]

\[ \log(l_1)=\sum_j n_j \log \left( \widehat{\pi_j} \right) \]

\[ G^2=2\log\left( \frac{l_1}{l_0}\right)=2\sum_j n_j \log \left(\frac{\widehat{\pi_j}}{ \pi_j^0} \right) \]

\(G^2\) ma asymptotycznie rozkład \(\chi^2\) z \(I-1\) stopniami swobody. Rozmiar przestrzeni \(H_0\) jest 0 (punkt) rozmiar \(H_1\) jest \(I-1\) bo \(\sum_{j=1}^{I}\pi_j=1\)
Statystyki \(G^2\) i \(X^2\) są asympttotycznie równe (różnica zbiega według prawdopodobieństwa do 0). \(X^2\) szybciej zbiega do \(\chi^2\) niż \(G^2\).
Przykład

W 60 rzutach kością było 15 jedynek i po 9 dwójek, trójek itd. Oblicz statystyki \(X^2\) i \(G^2\) dla hipotezy, że kość nie jest symetryczna .

\(\mu_j=10\)

\[ X^2=\frac{(15-10)^2}{10}+5\frac{(9-10)^2}{10}=2.5+0.5=3 \]

\(\widehat{\pi_1}=0.25\), \(\widehat{\pi_j}=0.15\) dla \(j=2,...,6\)

\[ G^2=2*\left( 15\log(0.25*6)+5*9*\log(0.15*6) \right) = 2.6815 \]

Wartość krytyczna testu \(\chi^2\) z 5 stopniami swobody wynosi

qchisq(0.05,5,lower.tail = FALSE)
[1] 11.0705

Testowanie gdy \(\mu_j\) są estymowane

Gdy \(\pi_j^0=h_j(\theta)\) to estymator największej wiarygodności \(\widehat{\pi_j^0}=h_j(\widehat{\theta})\) gdzie \(\widehat{\theta}\) jest estymatorem NW parametru \(\theta\).
Zamieniając \(\mu_j\) w statystyce \(X^2\) przez \(\widehat{\mu_j}=n \widehat{\pi_j^0}\) zmieniamy rozkład \(X^2\). Ma on asymptotycznie rozkład \(\chi^2\) z \(I-1-p\) stopniami swobody, gdzie \(p=\dim(\theta)\).

Przykład

Reszty Pearsona dla parametru jednowymiarowego

Reszty Pearsona niestandaryzowane
\[ e_j=\frac{n_j-\widehat{\mu_j}}{\sqrt{\widehat{\mu_j}}} \]

e <- function(n,m){
  ln <- length(n)
  ew <- vector(length = ln)
  for (i in 1:ln) ew[i] <- (n[i]-m[i])/sqrt(m[i])
  chi <- sum(ew^2)
  fi <- sqrt(chi/sum(n))
  list(r=ew,chi=chi,fi=fi)
}
nn <- c(30,63,63)
mm <- c(38.1,39,78.9)
pp <- e(nn,mm)
pp$r
[1] -1.312268  3.843076 -1.790023
pp$chi
[1] 19.69546
pchisq(pp$chi,df = 1,lower.tail = FALSE)
[1] 9.081684e-06
pp$fi
[1] 0.3553209

Reszty Pearsona standaryzowane
\[ \widehat{\pi_j} = \pi_j(\widehat{\theta})\\ g_j=\frac{\partial}{\partial \theta}\left( \pi_j(\theta) \right) |_{\theta=\widehat{\theta}}\\ aa = \sum_j \frac{g_j^2}{\widehat{\pi_j}} \\ h_j=\widehat{\pi_j}+\frac{g_j^2}{aa*\widehat{\pi_j}}\\ r_j = \frac{e_j}{\sqrt{1-h_j}} \]

th <- 0.494
(pi <- c(th^2,th-th^2,1-th))
[1] 0.244036 0.249964 0.506000
(gpi <- c(2*th,1-2*th,-1))
[1]  0.988  0.012 -1.000
(aa <- 4+(1-2*th)^2/th/(1-th)+1/(1-th))
[1] 5.976861
h <- vector(length = 3)
for (i in 1:3) h[i] <- pi[i]+gpi[i]^2/pi[i]/aa
h
[1] 0.9132837 0.2500604 0.8366560
rst <- vector(length = 3)
for(i in 1:3) rst[i] <- pp$r[i]/sqrt(1-h[i])
rst
[1] -4.456276  4.437780 -4.429013

Tablice kontyngencji

Przykład 1 Physicians Health Study (choroba wieńcowa serca)

silny zwykły brak razem
placebo 18 171 10845 11034
aspiryna 5 99 10933 11037

Przykład 2 Palenie i płeć

Wylosowano 100 osób. Każdej z nich zadano dwa pytania: płeć i palenie

pali nie pali razem
kobieta 5 43 48
mężczyzna 8 44 52

Przykład 3 Rak i palenie

rak <- matrix(c(688,650,21,59),nrow=2,ncol=2,byrow=TRUE,
                dimnames=list(c("pali T","pali N"),(c("rak +","rak-"))))
#knitr::kable(rak, caption = "Rak i palenia")
# print(xtable(rak, digits=0),
#  type = "html", html.table.attributes = "border=2")
rak+ rak-
pali T 688 650
pali N 21 59

Modele dla tablic kontyngencji

Oznaczenia \[ \pi_{ij}=P\left( X=x_i,Y=y_j \right) \\ \pi_{i+} = P\left( X=x_i \right) = \sum_{j=1}^J \pi_{ij} \\ \pi_{+j} = P\left( Y=y_j \right) = \sum_{i=1}^I \pi_{ij} \\ \pi_{j|i} = P\left( Y=y_j | X=x_i\right) = \frac{\pi_{ij}}{\pi_{i+}} \] Są dwa możliwe sposoby zbierania informacji

  • Wartości zmiennej \(X\) są wybierane losowo (Przykład 2)
  • Wartości zmiennej \(X\) są wybierane deterministycznie (Przykład 1 i 3)

Niezależność i jednorodność

Gdy zmienna \(X\) jest wybierana deterministycznie mamy do czynienia z jednorodnością, która oznacza, że w każdym wierszu rozkład (warunkowy) zmiennej \(Y\) jest taki sam. Oznacza to, że dla kazdego \(j\): \[ \pi_{j|1} = \pi_{j|2} = \dots = \pi_{j|I} \]

Z tego wynika, że dla każdego \(i,j\)

\[ \pi_{j|i} = \pi_{+j} \\ \pi_{ij} = \pi_{i+}\pi_{+j} \]

Dowód: \[ q_j=\pi_{j|1} = \pi_{j|2} = \dots = \pi_{j|I}\\ q_j = \frac{\pi_{ij}}{\pi_{i+}}\\ \pi_{ij} = q_j\pi_{i+} \\ \pi_{+j}= \sum_{i=1}^I \pi_{ij} = \sum_{i=1}^I q_j\pi_{i+}= q_j=\pi_{j|i} \]

Gdy zmienna \(X\) jest wybierana losowo mamy do czynienia z niezależnością, która oznacza, że Zmienne \(X\) i \(Y\) są niezależne (w sensie rachunku prawdopodobieństwa)

Z tego wynika, że dla każdego \(i,j\)

\[ \pi_{ij} = \pi_{i+}\pi_{+j} \]

Rozkłady prawdopodobieństwa w tablicach kontyngencji

Niech \(N_{ij}\) oznacza niezależne zmienne losowe, określające \(\#\{X=x_i,Y=y_j\}\) w próbie o \(n\) elementach (\(\sum_{ij}N_{ij}=n\)). Niech \(EN_{ij}=\mu_{ij}\). Rozważane są dwa rozkłady:
Poissona: \[ P\left(N_{ij} = n_{ij} \right)= e^{-\mu_{ij}}\frac{\mu_{ij}^{n_{ij}}}{n_{ij}!} \] Wielomianowy, gdy zmienna \(X\) ma wartości wybierane losowo: \[ P\left(N_{11} = n_{11},N_{12} = n_{12},\dots,N_{IJ} = n_{IJ} \right)= n!\prod_{ij}\frac{\pi_{ij}^{n_{ij}}}{n_{ij}!},\\ \pi_{ij} = \frac{\mu_{ij}}{n} \]

Wielomianowy, gdy zmienna \(X\) ma wartości wybierane deterministycznie: \[ P\left(N_{i1} = n_{i1},N_{i2} = n_{i2},\dots,N_{iJ} = n_{iJ} \right)= n_{i+}!\prod_{j}\frac{\pi_{j|i}^{n_{ij}}}{n_{ij}!},\\ \pi_{j|i} = \frac{\mu_{ij}}{\mu_{i+}} \]

Badania biomedyczne

http://www.authorstream.com/Presentation/panstudio-1176382-rodzaje-i-metodyka-bada-klinicznych/

Jak to się zaczęło?: JPo II Wojnie Światowej w czasopiśmie British Medical Journal ukazuje się publikacja opisująca wyniki z pierwszego badania klinicznego z losowym doborem pacjentów do grup. Jest to zarazem zwiastun nowoczesnej antybiotykoterapii, bowiem badanie dotyczy skuteczności streptomycyny w leczeniu gruźlicy. Oxman AD, Sackett DL, Guyatt GH . Users’ guides to the medical literature. I. How to get started. The Evidence-Based Medicine Working Group . JAMA ( 1993 ); 270: 2093-5 .

Badania obserwacyjne: tego rodzaju badania przeprowadza się głównie w celu oceny szkodliwości i analizy czynników rokowniczych; badania takie robi się z braku możliwości wykonania badań eksperymentalnych z powodów etycznych; dla wiarygodności tych wyników badań, podobnie jak w przypadku badań eksperymentalnych, zasadnicze znaczenie ma właściwe dobranie grupy kontrolnej, tak aby grupy różniły się tylko ocenianą ekspozycją; w badaniach obserwacyjnych łatwiej popełnić błąd systematyczny.

Badania prospektywne (kohortowe, cohort study): ocenia się prospektywnie wystąpienie określonego punktu końcowego w grupach (kohortach) osób narażonych i nienarażonych na dany czynnik lub interwencję, np.: Badanie umieralności związanej z paleniem tytoniu – 50-letnie obserwacje lekarzy brytyjskich BMJ. 2004 ; 328: 1519
Przykład 1