Niet parametrische statistiek

Het is geweten dat PCBs (Polychlorinated biphenyls, o.a. gebruikt in koelelementen en smeermiddel) makkelijk opstapelen in het vetweefsel van garnalen. In dit experiment werden twee groepen van elk 18 samples van garnalen gekweekt in verschillende condities. De garnalen werden gerandomiseerd over de twee groepen, en de PCB concentratie werd gemeten door het vetweefsel te extraheren van 100g garnalen in een sample. De onderzoeksvraag is om na te gaan of de PCB concentratie beïnvloed wordt door de groeicondities. De PCB concentraties werden gemeten in pg/g vet.

shrimps <- read.table("/Users/koenvandenberge/Dropbox/PhD/Onderwijs/Statistiek Biochemie/201819/dropboxStats1819/class6_nonPar/full/shrimps.txt")
shrimps$group <- factor(shrimps$group)

Data-exploratie

table(shrimps$group)

 1  2 
18 18 

boxplot(PCB.conc ~ group, data=shrimps)

De ``outliers’’ in de data werden dubbel gecheckt door de onderzoekers en zijn correcte metingen.

par(mfrow = c(1,2))
qqnorm(shrimps$PCB.conc[shrimps$group == 1])
qqline(shrimps$PCB.conc[shrimps$group == 1])

qqnorm(shrimps$PCB.conc[shrimps$group == 2])
qqline(shrimps$PCB.conc[shrimps$group == 2])

par(mfrow = c(1,1))

Conclusie:

• De data lijken niet afkomstig te zijn uit een Normale verdeling;
• Er zijn onvoldoende observaties om te berusten op de centrale limiet stelling;
• De varianties in beide groepen zijn niet gelijk;
• Er zijn grote outliers.

Geen enkel van de parametrische toetsen die we reeds gezien hebben kunnen gebruikt worden voor deze data analyse. We zullen dus beroep doen op niet-parametrische methoden.

Permutatie-test

Bij een permutatietest testen we of er gelijkheid is in distributies tussen groepen. De nulhypothese kan men dus schrijven als:

\[ H_0: f_1(y) = f_2(y) \]

In de alternatieve hypothese wordt verondersteld dat beide groepen een gelijke distributie hebben, maar met een verschillend gemiddelde:

\[ H_1: f_1(y) = f_2(y + \Delta) \]

Hierbij is \(f_1(y)\) de distributie voor de PBC concentraties in groeiconditie 1, en \(f_2(y)\) voor groeiconditie 2. Onder de alternatieve hypothese, stelt \(\Delta\) het verschil in gemiddelde voor tussen groep 1 en groep 2. Merk op dat \(\Delta\) zowel positief als negatief kan zijn.

We zullen eerst de geobserveerde test-statistiek van de two-sample t-test berekenen.

t.test(PCB.conc~group, data=shrimps)

    Welch Two Sample t-test

data:  PCB.conc by group
t = 0.085098, df = 19.805, p-value = 0.933
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -14.90106  16.16773
sample estimates:
mean in group 1 mean in group 2 
       46.63889        46.00556 

t.obs <- t.test(PCB.conc~group, data=shrimps)$statistic

Aangezien niet aan de assumpties voor de two-sample t-test voldaan is, zal de theoretische nuldistributie niet opgaan voor de geobserveerde teststatistiek. Hierdoor zullen we zelf een nuldistributie berekenen op basis van de data.

N <- 10000 #1000 permutaties
t.star <- c() #initialiseer een vector met permutatie teststatistieken
group.labels <- shrimps$group #groep labels: om te permuteren
perm.data <- shrimps #kopieer originele dataset voor permutatie test
for(i in 1:N) { #doe 1 tot N keer
  perm.labels <- sample(group.labels,replace=F) # willekeurige herschikking van de groep labels
  perm.data$group <- perm.labels #voeg willekeurige groeplabels toe aan originele dataset
  t.star <- c(t.star,t.test(PCB.conc~group,data=perm.data)$statistic) #voeg permutatie teststatistiek toe aan vector.
}
hist(t.star) #de nuldistributie
abline(v=t.obs, col="red") #geobserveerde teststatistiek = rode lijn

We zien dat de geobserveerde \(t\)-teststatistiek niet extreem is in vergelijking met de nuldistributie. We kunnen de nuldistributie gebruiken om de (tweezijdige) p-waarde te berekenen:

pvalue <- 2*mean(t.star >= t.obs)
pvalue

[1] 0.926

We bekomen een heel hoge p-waarde. We vinden aldus onvoldoende bewijs voor een verschil in gemiddelde PBC concentratie tussen de twee groeicondities op het 5% significantieniveau (p = 0.02). Merk op dat iedereen een andere p-waarde zal uitkomen vanwege het willekeurig herordenen van de groep-labels in de simulatie. Om een reproduceerbare uitkomst te bekomen kan je een seed gebruiken a.d.h.v. set.seed (zie cursus)

Rank test

Bovenstaande permutatietest testte voor een verschil in gemiddelde. De test was echter computationeel intensief: we moesten duizenden keer opnieuw de teststatistiek berekenen. Voor deze studie was dat nog haalbaar, en verloopt het uitvoeren van de test over een tijdspanne van enkele seconden. Echter, in moderne hoog-dimensionale biologische datasets zijn permutatietesten vaak niet haalbaar. Denk maar aan genexpressie data, waar men een test zou kunnen uitvoeren voor meer dan \(50 000\) menselijke transcripten. Bovendien dient men in die context een erg accurate p-waarde te bekomen (dus gebaseerd op een groot aantal permutaties, zeker in de orde van miljoenen) aangezien de correctie voor meervoudig toetsen erg groot zal zijn. In deze sectie zullen we dus een alternatieve niet-parametrische toets uitvoeren, namelijk een rank test. Merk op dat permutatietesten ook enkel geldig zijn indien we het locatie-shift model veronderstellen, in andere woorden, indien we veronderstellen dat beide groepen een zelfde distributie hebben. Indien we deze assumptie ook nemen bij rank testen, kunnen we het resultaat van de rank test interpreteren in functie van een verschil in de mediaan. Het is logisch dat we hier geen interpretatie in functie van de gemiddelden zullen kunnen bekomen, aangezien we de gemiddelden niet gebruiken in het berekenen van de teststatistiek. Echter, rank-testen zijn ook geldig indien het locatie-shift model niet opgaat. In dat geval, kunnen we het resultaat interpreteren in functie van de probabilistische index \(P(Y_{1} \ge Y_2)\).

Bij rank tests wordt de data \(Y_i\) getransformeerd naar ranks \[ R_i=R(Y_i) = \#\{Y_j: Y_j\leq Y_i; j=1,\ldots, n\}. \] Ranks zijn erg robuust voor outliers: of de hoogste observatie nu een waarde van 100 of 10 heeft, het blijft dezelfde rank behouden. In het geval van ties, dit zijn gelijke waarden, worden midranks gebruikt: \[R(Y_i) = \frac{\sum\limits_{\forall j : Y_j=Y_i}R(Y_j)}{\#{\forall j:Y_j=Y_i}} \] dit is, de midrank is het gemiddelde van de ranks van de gelijke observaties.

Vervolgens kan men de gemiddelde rank vergelijken tussen de twee groepen aan de hand van de Wilcoxon-Mann-Whitney (WNW) test.

wilcox.test(PCB.conc~group, data = shrimps)

    Wilcoxon rank sum test

data:  PCB.conc by group
W = 88, p-value = 0.01871
alternative hypothesis: true location shift is not equal to 0

We vinden dat de test significant is (op het 5%-significantieniveau), maar wachten nog even met de interpretatie van het resultaat. De Wilcoxon test hierboven is standaard een exacte test (zie vierde paragraaf in sectie Details in ?wilcox.test). Dit wil zeggen dat de theoretische distributies voor onze steekproefgroottes gebruikt worden als nuldistributie, in plaats van permutaties. Hieronder gaan we na of het volgend resultaat ook bekomen wordt indien we zelf permuteren.

N <- 10000 #1000 permutaties
W.star <- c() #initialiseer een vector met permutatie teststatistieken
group.labels <- shrimps$group #groep labels: om te permuteren
perm.data <- shrimps #kopieer originele dataset voor permutatie test
for(i in 1:N) { #doe 1 tot N keer
  perm.labels<-sample(group.labels,replace=F) # willekeurige herschikking van de groep labels
  perm.data$group<-perm.labels #voeg willekeurige groeplabels toe aan originele dataset
  W.star[i] = wilcox.test(PCB.conc~group, data=perm.data)$statistic #voeg permutatie teststatistiek toe aan vector.
}
hist(W.star)
abline(v=88, col=2)

pvalue <- 2*mean(W.star<=88)
pvalue #via manuele permutaties vinden we inderdaad een gelijkaardige p-waarde

[1] 0.0172

De p-waarde is uiteraard niet exact gelijk omdat dit slechts een benadering is, maar komt wel dicht in de buurt! Merk op dat de teststatistiek kan bekomen worden door de som te nemen van de ranks van één van de groepen die men test. Aangezien we geen gemiddelden van de data gebruiken in de test, kunnen we de test uiteraard ook niet interpreteren in termen van deze gemiddelden. Via de ranks kunnen we het resultaat echter wel interpreteren aan de hand van de probabilistische index. De nulhypothese stelt dat de distributies van de probabiliteitsdensiteitsfuncties in beide groepen gelijk zijn.

\[ H_0: f_1 = f_2 \]

In woorden: de verdeling van PCB concentraties in garnalen zijn in beide condities gelijk. In termen van probabilistische index, kunnen we de nulhypothese ook formuleren als:

\[ H_0: P(Y_{i1} \geq Y_{i2}) = 1/2 \]

De alternatieve hypothese kan in termen van de probabilistische index kan men scrhijven als volgt:

\[ H_1: P(Y_{i1} \geq Y_{i2}) \ne 1/2 \]

In woorden: de kans dat een willekeurige observatie van de PCB concentratie van garnalen die opgegroeid werden in conditie 1 groter dan of gelijk is aan een willekeurige observatie van de PCB concentratie van garnalen die opgegroeid werden in conditie 2 is niet gelijk aan 50%.

Deze kans, die wordt aangenomen 50% te zijn onder de nulhypothese, kan men bekomen aan de hand van de geobserveerde teststatistiek:

n1 <- n2 <- 18 #18 observaties in elke groep
WObs <- wilcox.test(PCB.conc~group, data=shrimps)$statistic #geobserveerde teststatistiek
WObs/(n1*n2)

        W 
0.2716049 

De puntschatting voor deze kans is dus 27.2%.

We interpreteren dit als volgt: De kans dat de PCB concentratie in een willekeurige garnaal die werd opgegroeid in conditie 1 groter of gelijk is aan de PCB concentratie in een willekeurige garnaal die werd opgegroeid in conditie 2 is gelijk aan 27.2%.

Of nog:

De kans dat de PCB concentratie in een willekeurige garnaal die werd opgegroeid in conditie 1 kleiner is dan de PCB concentratie in een willekeurige garnaal die werd opgegroeid in conditie 2 is gelijk aan 72.8%.

Deze kans verschilt significant van 50% op het 5% significantieniveau (p=0.019).