Vorige hoofdstukken: parametrische methoden
Inferentie enkel correct als voldaan aan param. veronderstellingen:
De \(p\)-waarde: \(\text{P}_0\left[ \vert T\vert \geq \vert t \vert \right]\).
\(95\%\) BI steunt eveneens op veronderstellingen. Niet voldaan: geen garantie dat intervallen populatie parameterwaarde omvatten met \(95\%\) kans.
Asymptotische theorie moeilijker te plaatsen: je kan stellen dat \(t\)-test asymptotisch niet-parametrisch is omdat bij erg grote steekproefgroottes de distributionele veronderstelling van normaliteit niet meer belangrijk is.
Parametrische aanpak:
chol <- read_tsv("https://raw.githubusercontent.com/statOmics/sbc/master/data/chol.txt")
chol <- chol |>
mutate(group = as.factor(group))
nGroups <- chol |>
pull(group) |>
table()
n <- nGroups |>
sum()
chol# A tibble: 10 × 2
group cholest
<fct> <dbl>
1 1 244
2 1 206
3 1 242
4 1 278
5 1 236
6 2 188
7 2 212
8 2 186
9 2 198
10 2 160
Vraagstelling Cholestorol voorbeeld vertaling naar nulhypothese. Klassiek: \[H_0: \mu_1=\mu_2 \text{ versus } H_1: \mu_1\neq \mu_2.\]
We gaan de voorwaarden na:
\(Y_{1j}\) en \(Y_{2j}\) uitkomsten uit respectievelijk groep 1 en 2: \[Y_{1j} \text{ iid } N(\mu_1,\sigma^2) \;\;\;\text{ en }\;\;\; Y_{2j} \text{ iid } N(\mu_2,\sigma^2).\]
Onder \(H_0:\mu_1=\mu_2\), wordt dit (stel \(\mu=\mu_1=\mu_2\) onder \(H_0\)) \[ Y_{ij} \text{ iid } N(\mu,\sigma^2),\]
Drukt uit dat alle \(n=n_1+n_2\) uitkomsten uit zelfde normale distributie komen en onafhankelijk zijn
Laat toe om oorspronkelijke nulhypothese anders te schrijven: [H_0: f_1(y) = f_2(y) y ] met
Onder de alternatieve hypothese wordt een locatie-shift verondersteld: \[H_1: f_1(y)=f_2(y-\Delta) \;\;\;\text{ voor alle } y\] met
We illustreren dit in R voor \(f_1\sim N(0,1)\) en \(f_2\sim N(1,1)\) en \(\Delta=-1\).
mu1 <- 0; mu2 <- 1; sigma1 <- sigma2 <- 1
y <- -2:2
delta <- mu1-mu2; delta[1] -1rbind(dnorm(y,mu1,sigma1), dnorm(y-delta,mu2,sigma2)) [,1] [,2] [,3] [,4] [,5]
[1,] 0.05399097 0.2419707 0.3989423 0.2419707 0.05399097
[2,] 0.05399097 0.2419707 0.3989423 0.2419707 0.05399097De permutatietesten die we ontwikkelen kunnen gebruikt worden voor het testen van [ H_0: F_1=F_2 ;;;;;; H_0:f_1=f_2. ] maar zonder de Normaliteitsveronderstellingen.
We weten dat onder \(H_0\) geldt dat :
We kunnen alle m=252 permutaties in R genereren a.d.h.v. de functie combn(n,n_1). Dit wordt geïllustreerd in de onderstaande R code:
G <- combn(n,nGroups[1])
dim(G)[1] 5 252G[,1:10] [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 1 1 1 1 1 1 1 1 1
[2,] 2 2 2 2 2 2 2 2 2 2
[3,] 3 3 3 3 3 3 3 3 3 3
[4,] 4 4 4 4 4 4 5 5 5 5
[5,] 5 6 7 8 9 10 6 7 8 9We berekenen nu de teststatistiek voor elke permutatie
tOrig <- t.test(cholest ~ group, chol)$statistic
tOrig t
3.664425 tStar <- combn(
n,
nGroups[1],
function(g,y) {
t.test(y[g],y[-g])$statistic
},
y = chol |> pull(cholest)
)
head(tStar)[1] 3.6644253 1.6397911 2.3973923 1.5876250 1.9217173 0.9967105length(tStar)[1] 252
Permutatienuldistributie \(\longrightarrow\) hypothesetesten uitvoeren.
tweezijdige permutatie \(p\)-waarde \[p=\text{P}_0\left[\vert T\vert \geq \vert t\vert \mid \mathbf{y}\right].\]
p-waarde geconditioneerd op geobserveerde cholesterolwaarden \(\mathbf{y}=(y_{11},\ldots,y_{51},y_{12},\ldots,y_{52})^T\).
Gezien permutatienuldistrubutie van \(T\) bepaald wordt door \(t^*_g\), \(g \in{\cal{G}}\), berekenen we \[p = \frac{\#\{g\in {\cal{G}}: \vert t^*_g\vert \geq \vert t \vert \}}{m}\]
pval <- mean(abs(tStar)>=abs(tOrig))
pval[1] 0.01587302Op het \(5\%\) significantieniveau besluiten we dat de distributies van de cholesterol concentraties niet gelijk zijn bij hartpatiënten en bij gezonde personen. (\(p=\) 0.0159).
De \(p\)-waarde op basis van alle permutaties wordt een exacte \(p\)-waarde genoemd.
De permutatienuldistributie wordt een exacte nuldistributie genoemd.
De term exact betekent dat de resultaten correct zijn voor iedere steekproefgrootte \(n\).
De kans op een type I fout wordt dus gecontroleerd door een permutatietest, maar wel conditioneel op de geobserveerde uitkomsten data \(\mathbf{y}\).
We kunnen ons nu afvragen of we de conclusies kunnen veralgemenen naar de populatie toe? Het antwoord is ja, als de subjecten at random getrokken zijn uit de populatie.
Het bewijs hiervan valt buiten het bestek van de cursus.
Soms probleem omdat het aantal permutaties \(m=\#{\cal{G}}\) erg groot is
Daarom niet alle \(g \in {\cal{G}}\) te beschouwen, maar groot aantal random permutaties uitvoeren (b.v. 10000)
Ter illustratie staat de code hiervan ook in de cursus.
We permuteren immers t statistiek, effectgrootte locatie-shift
Als locatie-shift model niet opgaat is een significant resultaat moeilijk te interpreteren, tegenover welk alternatief testen we dan???
Locatie-shift kan niet worden nagegaan in kleine steekproeven…
Rank testen starten vanuit rank-getransformeerde uitkomsten.
chol |>
pull(cholest) [1] 244 206 242 278 236 188 212 186 198 160chol |>
pull(cholest) |>
rank() [1] 9 5 8 10 7 3 6 2 4 1Soms komen ties voor in de data, i.e. minstens twee observaties hebben dezelfde numerieke waarde. Een klein voorbeeld:
metTies <- c(403, 507, 507, 610, 651, 651, 651, 830, 900)
metTies |> rank()[1] 1.0 2.5 2.5 4.0 6.0 6.0 6.0 8.0 9.0Ties: 507 komt tweemaal voor, 651 komt driemaal voor. Dit zijn voorbeelden van ties.
Wanneer ties voorkomen, worden midranks gebruikt.
midrank van observatie \(Y_i\) wordt \[\begin{eqnarray*} R_i &=& \frac{ \#\{Y_j: Y_j\leq Y_i\} + ( \#\{Y_j: Y_j < Y_i\} +1)}{2}. \end{eqnarray*}\]
t(chol) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
group "1" "1" "1" "1" "1" "2" "2" "2" "2" "2"
cholest "244" "206" "242" "278" "236" "188" "212" "186" "198" "160"z <- chol |>
pull(cholest)
z [1] 244 206 242 278 236 188 212 186 198 160rank(z) [1] 9 5 8 10 7 3 6 2 4 1\(R_{ij}=R(Y_{ij})\) is de rank van uitkomst \(Y_{ij}\) in de gepoolde steekproef. [ T = {i=1}^{n_1} R(Y{i1}) - {i=1}^{n_2} R(Y{i2}) . ]
Onder \(H_0\) verwachten we dat gemiddelde rank in de eerste groep ongeveer gelijk is aan de gemiddelde rank in de tweede groep en \(T\) dicht bij nul.
Als \(H_1\) waar is dan verwachten we dat gemiddelde ranks verschillen en dus dat \(T\) niet dicht bij nul zal liggen.
Eigenlijk voldoende om enkel \[S_1=\sum_{i=1}^{n_1} R(Y_{i1})\].
\(S_1\) is som van de ranks van observaties uit groep 1 vandaar de naam rank sum test.
Want [ S_1+S_2 = = 1+2++ n=n(n+1). ]
Vaak wordt gestandaardiseerde teststatistiek gebruikt [ T = , ]
met \(\text{E}_{0}\left[S_1\right]\) en \(\text{Var}_{0}\left[S_1\right]\) de verwachtingswaarde en variantie van \(S_1\) onder \(H_0\).
Dit zijn dus het gemiddelde en variantie van de permutatienuldistributie van \(S_1\).
Onder \(H_0\) geldt [ {0}= n_1(n+1) ;;;;;;;; {0}=n_1n_2(n+1). ]
Verder kan men onder \(H_0\) en als \(\min(n_1,n_2)\rightarrow \infty\) opgaat aantonen dat, [ T = N(0,1).
]
Asymptotisch volgt gestandaardiseerde teststatistiek een standaardnormaal verdeling!
We illustreren de WMW test aan de hand van de R functie wilcox.test.
wilcox.test(cholest ~ group, data = chol)
Wilcoxon rank sum exact test
data: cholest by group
W = 24, p-value = 0.01587
alternative hypothesis: true location shift is not equal to 0We verwerpen \(H_0\) (\(p=\) 0.016 \(<0.05\))
De output geeft de teststatistiek \(W=\) 24?
In volgende lijnen berekenen we \(S_1\) en \(S_2\) manueel voor de dataset.
S1 <- chol |>
mutate(cholRank=rank(cholest)) |>
filter(group==1) |>
pull(cholRank) |>
sum()
S2 <- chol |>
mutate(cholRank = rank(cholest)) |>
filter(group == 2) |>
pull(cholRank) |>
sum()
S1[1] 39S2[1] 16Waar komt \(W=\) 24 vandaan?
Mann en Whitney test in afwezigheid van ties: [ U_1 = {i=1}{n_1}_{k=1}{n_2} {Y{i1}Y_{k2}}. ]
\(\text{I}\left\{.\right\}\): indicator is 1 is als uitdrukking waar is en 0 als dit niet het geval is.
U telt hoeveel keer observatie uit groep 1 \(\geq\) observatie uit groep 2.
y1 <- chol |>
filter(group==1) |>
pull("cholest")
y2 <- chol |>
filter(group==2) |>
pull("cholest")
u1Hlp <- sapply(
y1,
function(y1i,y2) {y1i>=y2},
y2=y2)
colnames(u1Hlp) <- y1
rownames(u1Hlp) <- y2u1Hlp 244 206 242 278 236
188 TRUE TRUE TRUE TRUE TRUE
212 TRUE FALSE TRUE TRUE TRUE
186 TRUE TRUE TRUE TRUE TRUE
198 TRUE TRUE TRUE TRUE TRUE
160 TRUE TRUE TRUE TRUE TRUEU1 <- sum(u1Hlp)
U1[1] 24Er kan worden aangetoond dat
\[U_1 = S_1 - \frac{1}{2}n_1(n_1+1).\]
S1-nGroups[1]*(nGroups[1]+1)/2 1
24 \(U_1\) heeft interpretatievoordeel
Stel \(Y_j\) een willekeurige uitkomst uit behandelingsgroep \(j\) (\(j=1,2\)). Dan geldt \[\begin{eqnarray*} \frac{1}{n_1n_2}\text{E}\left[U_1\right] &=& \text{P}\left[Y_1 \geq Y_2\right]. \end{eqnarray*}\]
Intuïtief voelen we dit aan:
Op basis van steekproef kunnen we kans: gemiddelde van alle indicator waarden \(\text{I}\left\{Y_{i1}\geq Y_{k2}\right\}\).
\(n_1 \times n_2\) vergelijkingen
mean(u1Hlp)[1] 0.96U1/(nGroups[1]*nGroups[2]) 1
0.96 [] wordt een probabilistische index (Engels: probabilistic index) genoemd.
Kans dat een uitkomst uit groep1 \(\geq\) dan uitkomst uit groep 2
Onder \(H_0\): [=]
wilcox.test geeft niet de Wilcoxon rank sum statistiek, maar wel de Mann-Whitney statistiek \(U_1\).wTest <- wilcox.test(cholest~group, data = chol)
wTest
Wilcoxon rank sum exact test
data: cholest by group
W = 24, p-value = 0.01587
alternative hypothesis: true location shift is not equal to 0U1[1] 24probInd <- wTest$statistic/prod(nGroups)
probInd W
0.96 Aangezien \(p=\) 0.0159 \(<0.05\) besluiten we op het \(5\%\) significantieniveau dat de gemiddelde cholesterolconcentratie groter is bij hartpatiënten kort na een hartaanval dan bij gezonde personen. (We nemen aan dat locatie-shift opgaat)
We weten ook dat een cholesterolwaarde van hartpatiënten met een kans van \(U1/(n_1\times n_2)=\) 96% groter is die van gezonde personen.
We zouden de veronderstelling van de locatie-shift moeten nagaan, maar met slechts 5 observaties in elke behandelingsgroep is dit zinloos.
Zonder locatie-shift veronderstelling blijft de conclusie in termen van de probabilistische index correct!
Dus wanneer we geen locatie-shift veronderstellen en een tweezijdige test uitvoeren testen we eigenlijk
[H_0: F_1=F_2 (Y_1 Y_2) .]
Conclusie: Er is een significant verschil in de distributie van de cholesterolconcentraties bij hartpatiënten 2 dagen na hun hartaanval en gezonde individuen (\(p=\) 0.0159). Het is meer waarschijnlijk om een hogere cholesterolconcentraties te observeren hartpatiënten dan bij gezonde individuen. De puntschatting voor deze kans bedraagt 96%.
1,2-dimethylhydrazine dihydrochloride (DMH) testen op genotoxiciteit (EU directive)
24 ratten
Vier groepen volgens dagelijkse DMH dosis
Genotoxiciteit in de lever a.d.h.v. een comet assay op 150 levercellen per rat.
De onderzoekers wensen na te gaan of verschillen zijn in de DNA schade tengevolge van de DMH dosis.
Comet Assay:
dna <- read_delim("https://raw.githubusercontent.com/statOmics/sbc/master/data/dna.txt",delim=" ")
dna <- dna |>
mutate(dose = as.factor(dna$dose))
dna# A tibble: 24 × 3
id length dose
<chr> <dbl> <fct>
1 Rat1 19.3 0
2 Rat2 18.9 0
3 Rat3 18.6 0
4 Rat4 19.0 0
5 Rat5 17.1 0
6 Rat6 18.8 0
7 Rat7 55.4 1.25
8 Rat8 59.2 1.25
9 Rat9 59.1 1.25
10 Rat10 52.1 1.25
# ℹ 14 more rows
De Kruskal-Wallis Rank Test (KW-test) is een niet-parametrisch alternatief voor de ANOVA F-test.
klassieke \(F\)-teststatistiek kan geschreven worden als [ F = = , ]
met \(g\) het aantal groepen.
SSTot hangt enkel af van uitkomsten \(\mathbf{y}\) en zal niet variëren bij permutaties.
Voldoende om SST als teststatistiek te gebruiken. \[\text{SST}=\sum_{j=1}^t n_j(\bar{Y}_j-\bar{Y})^2\]
Exacte KW-test via permutatienuldistributie (die enkel afhangt van \(n_1, \ldots, n_g\)) \[H_0: f_1=\ldots=f_g \text{ vs } H_1: \text{ minstens twee gemiddelden verschillend}.\]
Locatie-shift
Indien locatie-shift niet opgaat, moet \(H_1\) geformuleerd worden in termen van probabilistische indexen: \[H_0: f_1=\ldots=f_g \text{ vs } H_1: \exists\ j,k \in \{1,\ldots,g\} : \text{P}\left[Y_j\geq Y_k\right]\neq 0.5\]
kruskal.test(length ~ dose, data = dna)
Kruskal-Wallis rank sum test
data: length by dose
Kruskal-Wallis chi-squared = 14, df = 3, p-value = 0.002905Op \(5\%\) significantieniveau kan nulhypothese worden verworpen.
R-functie kruskal.test enkel asymptotische benadering voor berekening van \(p\)-waarden.
Met slechts 6 observaties per groep, is dit geen optimale benadering van de exacte \(p\)-waarde!
coin R package kan de exacte \(p\)-waarde wel berekenenlibrary(coin)
kwPerm <- kruskal_test(length ~ dose, data = dna,
distribution = approximate(nresample = 100000))
kwPerm
Approximative Kruskal-Wallis Test
data: length by dose (0, 1.25, 2.5, 5)
chi-squared = 14, p-value = 0.00047Welke groepen zijn verschillend???
Posthoc analyse a.d.h.v WMW testen
Eigenlijk niet optimaal na KW-test want de data worden opnieuw gerangschikt!
Betere alternatieven bestaan maar vallen buiten bestek van de cursus
pairwise.wilcox.test(dna$length,dna$dose)
Pairwise comparisons using Wilcoxon rank sum exact test
data: dna$length and dna$dose
0 1.25 2.5
1.25 0.013 - -
2.5 0.013 0.818 -
5 0.013 0.721 0.788
P value adjustment method: holm pairwise.wilcox.test output. Puntschatter op de kans op hogere DNA-schade?nGroep <- table(dna$dose)
probInd <- combn(
levels(dna$dose),
2,
function(x) {
test <- wilcox.test(
length~dose,
dna |> filter(dose%in%x)
)
return(test$statistic/prod(nGroep[x]))
})
names(probInd) <- combn(
dna |> pull(dose) |> levels(),
2,
paste,
collapse="vs")
probInd 0vs1.25 0vs2.5 0vs5 1.25vs2.5 1.25vs5 2.5vs5
0.0000000 0.0000000 0.0000000 0.4444444 0.2777778 0.3333333 Omdat er twijfels zijn of het locatie-shift model geldig is, doen we enkel uitspraken in termen van de probabilistische index.
we veronderstellen afwezigheid van ties↩︎