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
Estymatory NW parametrów modelu
\[ \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
\[ 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
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 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
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 |
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
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} \]
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+}} \]
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
Badania retrospektywne (badanie kliniczno-kontrolne, case-control study): porównanie grupy osób u których punkt końcowy już wystąpił (np. nowotwór), z grupą kontrolną (osoby zdrowe), a następnie sprawdzamy częstość występowania jakiegoś czynnika w obu grupach (karta choroby, badanie cech genetycznych, etc.), np.: Używanie telefonu komórkowego a ryzyko nowotworów mózgu
Przykład 3
Badania przekrojowe (cross-sectional study) : W badaniach tych ocenia się chorobowość np. liczba chorych na gruźlicę w Polsce w danym momencie. Badanie to daje dużą niepewność co do związków przyczyno-skutkowych między ekspozycją-ocenianym czynnikiem, np.: Ryzyko choroby zakrzepowo-zatorowej oraz jej profilaktyka w warunkach opieki szpitalnej.
Przykład 2
ggplot(data.frame(x=c(0, 1)), aes(x)) + stat_function(fun=function(x) 688*x/(688*x+650*(1-x))) +
geom_abline(slope = 1,intercept = 0,col="blue",linetype=2) +
labs(title="Prawdopodobieństwo raka, gdy pacjent pali",x="Prawdopodobieństwo raka w populacji",y="Prawdopodobieństwo raka dla palących")
\[ \pi_1=\pi_{1|1}=\frac{\pi_{11}}{\pi_{1+}} \\ \pi_2=\pi_{1|2}=\frac{\pi_{12}}{\pi_{2+}} \\ \]
Różnica i estymator różnicy w próbach o liczności \(n_1\) i \(n_2\) i liczbami sukcesóW \(n_{11}\) i \(n_{12}\) \[ \pi_1 - \pi_2 \\ \widehat{\pi_1}=\frac{n_{11}}{n_1}, \quad \widehat{\pi_2}=\frac{n_{12}}{n_2} \] Odchylenie standardowe różnicy \[ \sigma\left( \widehat{\pi_1} - \widehat{\pi_2} \right) = \sqrt { \frac{\pi_1(1-\pi_1)}{n_1} + \frac{\pi_2(1-\pi_2)}{n_2}} \]
\[ r = \frac{\widehat{\pi_1}}{\widehat{\pi_2}} \\ \sigma(\log(r)) = \sqrt{\frac{1-\pi_1}{\pi_1n_1}+\frac{1-\pi_2}{\pi_2n_2}} \]
Iloraz krzyżowy=1 \(\iff\) jednorodność.
Iloraz krzyżowy=1 \(\iff\) niezależność.
Estymator ilorazu krzyżowego
Iloraz krzyżowy dla doświadczeń typy case-control
Przykład. Rak i palenie
rak+ | rak- | |
---|---|---|
pali T | 688 | 650 |
pali N | 21 | 59 |
pi1 <- round(688/(688+650),digits=4)
pi2 <- round(21/(21+59),digits=4)
cat("\n",
"ryzyko raka wśród palących ",pi1,"\n",
"ryzyko raka wśród niepalących ",pi2,"\n",
"różnica ryzyk dla palących ",pi1-pi2,"\n",
"ryzyko względne ",round(pi1/pi2,digits=4),"\n",
"stosunek szans raka dla palących ",round(pi1/(1-pi1),digits=4),"\n",
"stosunek szans raka dla niepalących",round(pi2/(1-pi2),digits=4),"\n",
"iloraz krzyżowy ", round(pi1*(1-pi2)/pi2/(1-pi1),digits=4), "\n"
)
ryzyko raka wśród palących 0.5142
ryzyko raka wśród niepalących 0.2625
różnica ryzyk dla palących 0.2517
ryzyko względne 1.9589
stosunek szans raka dla palących 1.0585
stosunek szans raka dla niepalących 0.3559
iloraz krzyżowy 2.9738
Przykład
a1 <- 0.410; a2 <- 0.401
b1 <- 0.010; b2 <- 0.001
c(a1-a2,b1-b2)
[1] 0.009 0.009
c(a1/a2,b1/b2)
[1] 1.022444 10.000000
Stosunek szans dla palenia występowanie raka \[ \frac {P(X_1|Y_1)}{P(X_2|Y_1)}=\frac{\frac{\pi_{11}}{\pi_{+1}}}{\frac{\pi_{21}}{\pi_{+1}}}=\frac{\pi_{11}}{\pi_{21}} \]
Stosunek szans dla palenia niewystępowanie raka \[ \frac {P(X_1|Y_2)}{P(X_2|Y_2)}=\frac{\frac{\pi_{12}}{\pi_{+2}}}{\frac{\pi_{22}}{\pi_{+2}}}=\frac{\pi_{12}}{\pi_{22}} \]
Niepotrzebna jest znajomość \(\pi_{+1}\) - prawdopodobieństwa zachorowania na raka!
Estymatory ilorazów szans
c(688/21,650/59)
[1] 32.76190 11.01695
Gdybyśmy mieli w próbie pięciokrotnie mniejszą grupę z osób nie mających raka (z zachowanymi proporcjami)
rak+ | rak- | |
---|---|---|
pali T | 688 | 130 |
pali N | 21 | 12 |
Estymatory ilorazów szans
c(688/21,130/12)
[1] 32.76190 10.83333
Ilorazy krzyżowe w obu przypadkach.
688*59/(650*21)
[1] 2.973773
688*12/(130*21)
[1] 3.024176
Można dowolnie sterować liczbą osób w próbie.
Jeżeli jest reprezentatywna, a więc zachowały się proporcje, to stosunki szans i iloczyny krzyżowe nie zmienią się.
Przedziały ufności dla ilorazu krzyżowego
\[ OR=\frac{n_{11}*n_{22}}{n_{12}*n_{21}} \\ OR= \frac{(n_{11}+0.5)*(n_{22}+0.5)}{(n_{12}+0.5)*(n_{21}+0.5)} \\ \sigma(\log(OR)) = \sqrt{\frac{1}{n_{11}}+\frac{1}{n_{12}}+\frac{1}{n_{21}}+\frac{1}{n_{22}}} \]
oddsratioWald.proc <- function(n00, n01, n10, n11, alpha = 0.05){
#
# Compute the odds ratio between two binary variables, x and y,
# as defined by the four numbers nij:
#
# n00 = number of cases where x = 0 and y = 0
# n01 = number of cases where x = 0 and y = 1
# n10 = number of cases where x = 1 and y = 0
# n11 = number of cases where x = 1 and y = 1
#
OR <- (n00 * n11)/(n01 * n10)
#
# Compute the Wald confidence intervals:
#
siglog <- sqrt((1/n00) + (1/n01) + (1/n10) + (1/n11))
zalph <- qnorm(1 - alpha/2)
logOR <- log(OR)
loglo <- logOR - zalph * siglog
loghi <- logOR + zalph * siglog
#
ORlo <- exp(loglo)
ORhi <- exp(loghi)
#
oframe <- data.frame(LowerCI = ORlo, OR = OR, UpperCI = ORhi, alpha = alpha)
oframe
}
Przykład rak i palenie
oddsratioWald.proc(688,650,21,59)
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 |
Po połączeniu dwóch kolumn otrzymamy tabelę binarną
choroba | brak | razem | |
---|---|---|---|
placebo | 189 | 10845 | 11034 |
aspiryna | 104 | 10933 | 11037 |
pi1 <- round(189/(189+10845),digits=4)
pi2 <- round(104/(104+10933),digits=4)
cat("\n",
"ryzyko choroby wśród placebo ",pi1,"\n",
"ryzyko choroby wśród aspiryna ",pi2,"\n",
"różnica ryzyk dla placebo ",pi1-pi2,"\n",
"ryzyko względne ",round(pi1/pi2,digits=4),"\n",
"stosunek szans choroby dla placebo ",round(pi1/(1-pi1),digits=4),"\n",
"stosunek szans choroby dla aspiryny",round(pi2/(1-pi2),digits=4),"\n",
"iloraz krzyżowy ", round(pi1*(1-pi2)/pi2/(1-pi1),digits=4), "\n"
)
ryzyko choroby wśród placebo 0.0171
ryzyko choroby wśród aspiryna 0.0094
różnica ryzyk dla placebo 0.0077
ryzyko względne 1.8191
stosunek szans choroby dla placebo 0.0174
stosunek szans choroby dla aspiryny 0.0095
iloraz krzyżowy 1.8334
Gdy ryzyko zdarzenia małe (w obu przypadkach) to ryzyko względne = iloraz krzyżowy
Przykład 1 Physicians Health Study
n1 <- 189+10845
n2 <- 104+10933
roznica <- pi1-pi2 + c(-1,1)*1.96*sqrt(pi1*(1-pi1)/n1 + pi2*(1-pi2)/n2)
stosunek <- log(pi1/pi2) + c(-1,1)*1.96*sqrt((1-pi1)/pi1/n1+(1-pi2)/pi2/n2)
stosunek <- exp(stosunek)
OR <- oddsratioWald.proc(189,10845,104,10933)
ORpu <- c(OR$LowerCI,OR$UpperCI)
cat("\n",
"różnica ryzyk dla placebo ",round(roznica,digits=4),"\n",
"ryzyko względne ",round(stosunek,digits=4),"\n",
"iloraz krzyżowy ",round(ORpu,digits=4), "\n"
)
różnica ryzyk dla placebo 0.0047 0.0107
ryzyko względne 1.4337 2.3082
iloraz krzyżowy 1.44 2.3308
\[ \chi^2_{(I-1)(J-1)} \\ V = \sqrt{\frac{\chi^2}{n(\min(I,J)-1)}} \\ p_{i+} = \frac{n_{i+}}{n} \quad p_{+j} = \frac{n_{+j}}{n} \\ \widehat{\mu_{ij}}=\frac{n_{i+} n_{+j}}{n} \\ e_{ij}=\frac{n_{ij}-\widehat{\mu_{ij}}}{\sqrt{\widehat{\mu_{ij}}}} \\ r_{ij} = \frac{e_{ij}}{\sqrt{(1-p_{i+})(1-p_{+j})}} \]
Przykład 1 Physicians Health Study
silny | zwykły | brak | razem | |
---|---|---|---|---|
placebo | 18 | 171 | 10845 | 11034 |
aspiryna | 5 | 99 | 10933 | 11037 |
razem | 23 | 270 | 21778 | 22071 |
phs <- matrix(c(18,171,10845,5,99,10933),nrow = 2,ncol = 3,byrow = TRUE,
dimnames=list(c("placebo","aspiryna"),c("silny","zwykły","brak")))
ct <- chisq.test(phs)
ct$statistic
X-squared
26.90301
cat("\n","V = ",sqrt(ct$statistic/22071),"\n",
"p-value = ",ct$p.value,"\n")
V = 0.03491318
p-value = 1.439084e-06
# mij#
ct$expected
silny zwykły brak
placebo 11.49844 134.9817 10887.52
aspiryna 11.50156 135.0183 10890.48
# eij
ct$residuals
silny zwykły brak
placebo 1.917337 3.100177 -0.4075003
aspiryna -1.917076 -3.099755 0.4074449
# rij
ct$stdres
silny zwykły brak
placebo 2.712753 4.411078 -5.001388
aspiryna -2.712753 -4.411078 5.001388
FEV1 % normalnej pojemnosci płuc skorygowanej dla wieku i wzrostu pylica - “+” występuje, “-” brak
[Dane: Pylica u górników, Campbell, Machin, Medical Statistics, str.44]
FEV1 <- c(40,43,47,49,50,50,53,57,58,58,58,62,65,69,71,73,74,75,75,77,78,79,80,87,90,100,105,60,67,73,75,79,80,83,87,89,100,105,109,115)
pylica <- c(rep("+",27),rep("-",13))
pylica <- data.frame(FEV1=FEV1,pylica=pylica)
pylica
qplot(pylica, FEV1, data = pylica, geom = "violin", fill = pylica)
qplot(pylica, FEV1, data = pylica,
geom= "boxplot", fill = pylica)
Diagnoza na podstawie FEV1: od kiedy uznawać wystąpienie pylicy? Wprowadźmy dwa typy diagnozy:
D80: FEV1<=80 i D90: FEV1<=90
pylica$D80 <- case_when(
FEV1<=80 ~ "+",
TRUE ~ "-"
)
pylica$D80 <- as.factor(pylica$D80)
pylica$D90 <- case_when(
FEV1<=90 ~ "+",
TRUE ~ "-"
)
pylica$D90 <- as.factor(pylica$D90)
pylica$pylica <- factor(pylica$pylica, levels=rev(levels(pylica$pylica)))
pylica$D80 <- factor(pylica$D80, levels=rev(levels(pylica$D80)))
pylica$D90 <- factor(pylica$D90, levels=rev(levels(pylica$D90)))
Tabele kontyngencji
t80 <- table(pylica$D80,pylica$pylica,dnn=list("diagnoza","choroba"))
t80
choroba
diagnoza + -
+ 23 6
- 4 7
t90 <- table(pylica$D90,pylica$pylica,dnn=list("diagnoza","choroba"))
t90
choroba
diagnoza + -
+ 25 9
- 2 4
choroba + | choroba - | ||
---|---|---|---|
diagnoza + | TP | FP | DP |
diagnoza - | FN | TN | DN |
CP | CN | S |
Wskaźniki oceny tablicy diagnostycznej \[ \text{czułość:}\quad TPR=\frac{TP}{CP}\\ \text{specyficzność:}\quad TNR=\frac{TN}{CN}\\ FPR=\frac{FP}{CN}=1-TNR \]
Dla diagnozy D80 i D90
parD <- data.frame(TPR=c(23/27,25/27),TNR=c(7/13,4/13),FPR=c(6/13,9/13))
rownames(parD) <- c("D80","D90")
round(parD,digits=3)
Globalne oceny tablicy diagnostycznej \[ \text{dokładność:}\quad ACC=\frac{TP+TN}{S} \\ \text{informatywność:} \quad INF=TPR-FPR \\ \text{korelacja Matthewsa} \quad MCC = \phi=\frac{TP*TN-FP*FN}{\sqrt{DP*DN*CP*CN}} \]
Niech \(X=1\) gdy diagnoza jest “+” i \(X=0\) gdy diagnoza jest “-”. Podobnie, niech \(Y=1\) gdy choroba jest “+” i \(Y=0\) gdy choroba jest “-”. Wtedy współczynnik \(\phi\) jest współczynnikiem korelacji między \(X\) i \(Y\). Ponadto \(|\phi|\) jest równe statystyce \(V\) Cramera dla testowania zależności między \(X\) i \(Y\).
Dla diagnozy D80 i D90
gloD <- data.frame(ACC=c((23+7)/40,(25+4)/40),INF=parD$TPR-parD$FPR,
MCC=c((23*7-6*4)/sqrt(11*29*13*27),(4*25-4*9)/sqrt(6*34*13*27))
)
rownames(gloD) <- c("D80","D90")
round(gloD,digits=3)
Jak wybrać najlepszy próg diagnostyczny?
pylica$pylica1 <- 2-as.numeric(pylica$pylica)
ggroc <- ggplot(pylica, aes(m =-FEV1, d = pylica1)) + geom_roc()+style_roc()
ggroc
round(calc_auc(ggroc),digits=2)
pred <- prediction(-pylica$FEV1,2-as.numeric(pylica$pylica))
perf <- performance(pred,measure="tpr", x.measure="fpr")
prfdf <- data.frame(cutoff=unlist(perf@alpha.values),fpr=unlist(perf@x.values),tpr=unlist(perf@y.values))
prfdf <- cbind(prfdf,inf=prfdf$tpr-prfdf$fpr)
prfdf
Najlepszy punkt podziału: 78
pylica$D78 <- case_when(
pylica$FEV1<=78 ~ "+",
TRUE ~ "-"
)
pylica$D78 <- as.factor(pylica$D78)
pylica$D78 <- factor(pylica$D78, levels=rev(levels(pylica$D78)))
t78 <- table(pylica$D78,pylica$pylica,dnn=list("diagnoza","choroba"))
t78
choroba
diagnoza + -
+ 21 4
- 6 9
parD <- data.frame(TPR=c(23/27,25/27,21/27),TNR=c(7/13,4/13,9/13),FPR=c(6/13,9/13,4/13))
rownames(parD) <- c("D80","D90","D78")
round(parD,digits=3)
gloD <- data.frame(ACC=c((23+7)/40,(25+4)/40,(21+9)/40),INF=parD$TPR-parD$FPR,
MCC=c((23*7-6*4)/sqrt(11*29*13*27),(4*25-4*9)/sqrt(6*34*13*27),
(9*21-4*6)/sqrt(25*15*13*27))
)
rownames(gloD) <- c("D80","D90","D78")
round(gloD,digits=3)
W ankiecie wzięło udział 5 kobiet i 10 mężczyzn. Na pytanie o poparcie partii XX 3 kobiety i 5 mężczyzn odpowiedziały TAK. Czy poparcie partii XX jest większe u kobiet?
Tablica kontyngencji
tfis <- matrix(c(3,2,5,5,5,10,8,7,15),nrow=3,ncol=3,byrow = TRUE,dimnames=list(c("K","M","suma"),c("T","N","suma")))
tfis
T N suma
K 3 2 5
M 5 5 10
suma 8 7 15
Pytanie to można wyrazić w języku testowania: \[ H_0: \theta=1 \\ H_1: \theta>1 \]
Liczności w tablicy są tak małe, że poprzednie testy, jako asymptotyczne, są nieadekwatne.
Pomysł FIshera
Gdy \(H_0\) jest prawdziwa to prawdopodobieństwo wyprodukowania tablicy kontyngencji
\(y_1\) | \(y_2\) | |
---|---|---|
\(x_1\) | \(n_{11}\) | \(n_{12}\) |
\(x_2\) | \(n_{21}\) | \(n_{22}\) |
o ustalonych wartościach brzegowych \(n_{1+}\), \(n_{2+}\), \(n_{+1}\), \(n_{+2}\) wynosi \[ \frac{\binom{n_{+1}}{n_{11}}\binom{n_{+2}}{n_{12}}}{\binom{n_{++}}{n_{+1}}} = \frac{n_{+1}!\,n_{+2}!\,n_{1+}!\,n_{2+}!}{n_{11}!\,n_{12}!\,n_{21}!\,n_{22}!\,n!} \]
ffisher <- function(n11,n12,n21,n22,cyfr=4) {
n <- n11+n12+n21+n22
fff <- factorial(n11+n12)*factorial(n21+n22)* factorial(n11+n21)*factorial(n22+n12)/(
factorial(n11)*factorial(n12)* factorial(n21)*factorial(n22)*factorial(n)
)
round(fff,digits=cyfr)
}
ffisher(3,2,5,5)
[1] 0.3916
Warunek \(\theta>1\) jest równoważny warunkowi: \[ n_{11}>\widehat{n_{11}}=\frac{n_{1+}n_{+1}}{n} \] Warunek ten w naszym przykładzie ma postać
\[ n_{11}>\frac{5*8}{15}=\frac{8}{3}\\ n_{11} \geq 3 \]
Prawdopodobieństwo, że \(n_{11} \geq 3\) jest równe sumie prawdopodobieństw dla \(n_{11} = 3,4,5\) i wynosi
ffisher(3,2,5,5)
[1] 0.3916
ffisher(4,1,4,6)
[1] 0.1632
ffisher(5,0,3,7)
[1] 0.0186
ffisher(3,2,5,5) + ffisher(4,1,4,6) + ffisher(5,0,3,7)
[1] 0.5734
fisher.test(tfis[1:2,1:2],alternative = "greater")
Fisher's Exact Test for Count Data
data: tfis[1:2, 1:2]
p-value = 0.5734
alternative hypothesis: true odds ratio is greater than 1
95 percent confidence interval:
0.1543702 Inf
sample estimates:
odds ratio
1.460011
https://rstudio-pubs-static.s3.amazonaws.com/65435_a8b26773b5d64138b6411a5aa306a5b9.html https://css.ncifcrf.gov/services/alvord/R_Tutorial_Partition_LR_Chi-Square_IxJ_Contingency_Table.pdf
Dane [1996 (USA) General Social Survey, National Opinion Research Center.]
wiersze: wykształcenie“podstawowe”,“średnie”,“wyższe”
kolumny: stosunek do religii “fundamentalista”,“umiarkowany”,“liberalny”
rel.dat <- as.table(matrix(c(178,138,108,570,648,442,138,252,252),nrow=3,ncol=3,byrow=TRUE,dimnames = list(edu=c("P","S","W"),rel=c("F","U","L"))))
addmargins(rel.dat)
rel
edu F U L Sum
P 178 138 108 424
S 570 648 442 1660
W 138 252 252 642
Sum 886 1038 802 2726
Zadanie polega na wybraniu jednorodnych podmacierzy.
Twierdzenie o sumowaniu
Warunki dobrego podziału
t1 <- loglm( ~ edu + rel,data = rel.dat)
t1
Call:
loglm(formula = ~edu + rel, data = rel.dat)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 69.81162 4 2.486900e-14
Pearson 69.15676 4 3.419487e-14
chisq.test(rel.dat)$stdres
rel
edu F U L
P 4.5349181 -2.5521511 -1.9417050
S 2.5532682 1.2859224 -3.9946969
W -6.8097502 0.7009713 6.2525467
t2 <- as.table(rel.dat[c(1,2),c(2,3)])
t2
rel
edu U L
P 138 108
S 648 442
t2m <- loglm( ~ edu + rel,data = t2)
t2m
Call:
loglm(formula = ~edu + rel, data = t2)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 0.9267211 1 0.3357164
Pearson 0.9310773 1 0.3345831
chisq.test(t2)$stdres
rel
edu U L
P -0.9649235 0.9649235
S 0.9649235 -0.9649235
apply(t2,1,sum)
P S
246 1090
t3 <- as.table(matrix(c(246,1090,c(178,570)),nrow=2,ncol=2,dimnames = list(edu=c("P","S"),rel=c("U+L","F"))))
t3
rel
edu U+L F
P 246 178
S 1090 570
t3m <- loglm( ~ edu + rel,data = t3)
t3m
Call:
loglm(formula = ~edu + rel, data = t3)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 8.440411 1 0.003669733
Pearson 8.575911 1 0.003406395
chisq.test(t3)$stdres
rel
edu U+L F
P -2.928466 2.928466
S 2.928466 -2.928466
apply(t2,2,sum)
U L
786 550
t4 <- as.table(matrix(c(786,252,c(550,252)),nrow=2,ncol=2,dimnames = list(edu=c("P+S","W"),rel=c("U","L"))))
t4
rel
edu U L
P+S 786 550
W 252 252
t4m <- loglm( ~ edu + rel,data = t4)
t4m
Call:
loglm(formula = ~edu + rel, data = t4)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 11.55507 1 0.0006756460
Pearson 11.61005 1 0.0006559629
chisq.test(t4)$stdres
rel
edu U L
P+S 3.407353 -3.407353
W -3.407353 3.407353
sum(t2)
[1] 1336
t5 <- as.table(matrix(c(1336,252+252,c(178+570,138)),nrow=2,ncol=2,dimnames = list(edu=c("P+S","W"),rel=c("U+L","F"))))
t5
rel
edu U+L F
P+S 1336 748
W 504 138
t5m <- loglm( ~ edu + rel,data = t5)
t5m
Call:
loglm(formula = ~edu + rel, data = t5)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 48.88941 1 2.708056e-12
Pearson 46.37270 1 9.776846e-12
chisq.test(t5)$stdres
rel
edu U+L F
P+S -6.80975 6.80975
W 6.80975 -6.80975