--- title: "PC-practicum 3: gepaarde en ongepaarde testen" output: html_document --- # Voorbeeld 1 Een studie naar de kwaliteit van rivierwater geeft voor een reeks stromen een indicatie van de nitraatwaarde in het water. De stromen worden eerst ingedeeld in twee groepen naargelang ze initieel een lage (groep 1) of hoge (groep 2) nitraatwaarde (in mg N/l) hebben. Voor elke stroom is er tevens een meting in 1986 en in 1996 a anwezig. Sinds 1986 werden een aantal nieuwe maatregelen toegepast om te trachten de nitraatwaarde te drukken. Het doel is om te onderzoeken of er enige indicatie is of deze maatregelen effect hebben gehad. Lees de dataset nitrate2.dat in. ```{r} #nitrate <- read.table("nitrate2.dat", header = TRUE) nitrate <- read.table("nitrate2.dat",header = TRUE) attach(nitrate) ``` ``` #GROEP <- nitrate$GROEP #GVOOR <- nitrate$GVOOR #GNA <- nitrate$GNA ``` ```{r} nitrate ``` Bereken de gemiddelde wijziging in nitraatgehalte in de groep rivieren met initieel een lage nitraatwaarde (groep 1). ```{r} diff <- GVOOR-GNA mean(diff[GROEP == 1]) mean(GVOOR[GROEP == 1]) - mean(GNA[GROEP == 1]) ## gemiddelde van verschil = verschil van gemiddelde ``` De gemiddelde afname in nitraatwaarde in de rivieren met een initieel lage nitraatwaarde (groep 1) tussen 1986 en in 1996 is gelijk aan 0,701 mg N/l. Bereken de gemiddelde wijziging in nitraatgehalte in de groep rivieren met initieel een hoge nitraatwaarde (groep 2). ```{r} mean(diff[GROEP == 2]) mean(GVOOR[GROEP == 2]) - mean(GNA[GROEP == 2]) ``` De gemiddelde afname in nitraatwaarde in de rivieren met een initieel hoge nitraatwaarde (groep 2) tussen 1986 en in 1996 is gelijk aan 0,639 mg N/l. Bereken het gemiddelde verschil: het verschil van de gemiddelden! ```{r} (mean(GVOOR[GROEP == 1]) - mean(GNA[GROEP == 1])) - (mean(GVOOR[GROEP == 2]) - mean(GNA[GROEP == 2])) ``` De gemiddelde afname in nitraatwaarde in de rivieren met een initieel lage nitraatwaarde (groep 1) tussen 1986 en in 1996 is **0,062 mg N/l groter** dan de gemiddelde afname in nitraatwaarde in rivieren met een initieel hoge nitraatwaarde (groep 2). Maak een dot-plot van de gepaarde metingen voor beide groepen. ```{r} min(c(GVOOR,GNA)) ## vind de limieten om mooie figuur te maken max(c(GVOOR,GNA)) plot(GVOOR[GROEP == 1], ylim = c(14,31)) points(GNA[GROEP == 1], col = "Red") plot(GVOOR[GROEP == 2], ylim = c(14,31)) points(GNA[GROEP == 2], col = "Red") ``` ## Vraag A Test voor de rivieren met een initieel lage nitraatwaarde (groep 1) of de gemiddelde nitraatwaarde in 1996 anders is dan in 1986 op 5%-significantieniveau. **Hypotheses?** $H_0$: Voor rivieren met een initieel lage nitraatwaarde (groep 1) is de gemiddelde nitraatwaarde in 1986 gelijk aan de gemiddelde nitraatwaarde in 1996. $H_A$: Voor rivieren met een initieel lage nitraatwaarde (groep 1) is de gemiddelde nitraatwaarde in 1986 niet gelijk aan de gemiddelde nitraatwaarde in 1996. Test: We hebben gepaarde gegevens. We kiezen dus voor een gepaarde t-test. (de data is gepaard aangezien we de nitraatconcentraties hebben van exact dezelfde lokatie voor de jaren 1986 and 1996). Deze twee datapunten voor eenzelfde lokatie zijn allesbehalve onafhankelijk Voorwaarden: - gepaarde gegevens - onafhankelijke observaties - verschil normaal verdeeld ```{r} qqnorm(diff[GROEP == 1]) qqline(diff[GROEP == 1]) ``` ```{r} par(mfrow=c(3,3)) for(i in 1:9){ x <- rnorm(n=10, mean=mean(diff[GROEP == 1]), sd=sd(diff[GROEP == 1])) qqnorm(x) qqline(x) } ``` Voer de gepaarde test uit. ```{r} t.test(GVOOR[GROEP == 1], GNA[GROEP == 1], paired = TRUE) ``` Alternatieve methode; ongepaarde test op het verschil ```{r} t.test(diff[GROEP == 1]) ## is het verschil tussen 1986 en 1996 verschillend van nul ``` Interpreteer het betrouwbaarheidsinterval en formuleer een conclusie. - Interpretatie betrouwbaarheidsinterval: met 95% waarschijnlijkheid ligt de gemiddelde afname in nitraatwaarde voor rivieren met een initieel lage nitraatwaarde (groep 1) in 1986 vs. 1996 in het interval 0,39 mg N/l tot 1,01 mg N/l. - De gemiddelde afname in nitraatgehalte tussen 1986 en 1996 voor rivieren met een initieel lage nitraatwaarde is gelijk aan 0,701 mg N/l. De p-waarde is gelijk aan 0,0006. Deze afname in nitraatgehalte is dus extreem significant (p < 0,001) op het 5\% significantieniveau. ## Vraag B Test voor de rivieren met een initieel hoge nitraatwaarde (groep 2) of de gemiddelde nitraatwaarde na de maatregel **LAGER LIGT** dan ervoor op **1%-significantieniveau**. => eenzijdige test $H_0$: Voor rivieren met een initieel hoge nitraatwaarde (groep 2) is de gemiddelde nitraatwaarde in 1986 gelijk aan de gemiddelde nitraatwaarde in 1996. $H_A$: Voor rivieren met een initieel hoge nitraatwaarde (groep 2) is de gemiddelde nitraatwaarde in 1996 kleiner dan de gemiddelde nitraatwaarde in 1986. Test: Gepaarde t-test Voorwaarden: - gepaarde gegevens - onafhankelijke observaties (rivieren) - verschil normaal verdeeld ```{r} qqnorm(diff[GROEP == 2]) qqline(diff[GROEP == 2]) ``` De normaliteitsassumptie is moeilijk te beoordelen aangezien er weinig gegevens zijn! ```{r} ?t.test t.test(GNA[GROEP == 2], GVOOR[GROEP == 2], paired=TRUE, alternative = "less", conf.level = 0.99) ``` **Conclusie** Het 99\%-betrouwbaarheidsinterval bevat 0 niet en de p-waarde is kleiner dan 0.01, dus we zullen de nulhypothese verwerpen op het 1\% significantieniveau. Onze conclusie luidt dan als volgt: De gemiddelde afname in nitraatgehalte in rivieren met een initeel hoge nitraatwaarde (groep 2) in 1986 vs. 1996 is gelijk aan 0,64 mg N/l (95\% BI: min oneindig tot 0,23 mg N/l). Dit verschil is extreem significant (p = 0,0008) op het 1\% significantieniveau. **Merk op:** de keuze voor eenzijdig of tweezijdig testen moet altijd gemaakt worden voor je naar je data kijkt! Eenzijdig testen is enkel te verantwoorden indien je op voorhand bereid bent aan te nemen dat een zeer kleine p-waarde in één bepaalde richting altijd door willekeurig toeval veroorzaakt is en je in dat geval nog altijd de nulhypothese zult aanvaarden. In de praktijk komt dit zelden voor. Bijvoorbeeld in dit geval; We gaan maatregelen nemen om de nitraatconcentraties te verlagen. We achten het bijgevolg niet mogelijk dat ze zullen verhogen, of we verwachten dat een verhouding puur aan toeval te wijten zou zijn. Merk op dat in de praktijk we vaak voorzichting willen zijn, en er niet vanuit willen gaan dat een ingreep zal leiden tot het verlagen van de nitraatconcentraties. ## Vraag C Test of het gemiddeld verschil in nitraatgehalte in 1986 vs. 1996 verschillend is tussen groep 1 en 2. Bereken het gemiddelde verschil: het verschil van de gemiddelden! ```{r} (mean(GVOOR[GROEP == 1]) - mean(GNA[GROEP == 1])) - (mean(GVOOR[GROEP == 2]) - mean(GNA[GROEP == 2])) ``` Hypotheses? $H_0$: Het gemiddelde verschil in nitraatgehalte in 1986 vs. 1996 is gelijk voor rivieren in groep 1 en groep 2. $H_A$: Het gemiddelde verschil in nitraatgehalte in 1986 vs. 1996 is verschillend voor rivieren in groep 1 en groep 2 Test: Ongepaarde t-test (two-sample t-test) (Ongepaard aangezien de stromen die we nu bekijken onafhankelijk zijn van mekaar, m.a.w. riveiren van groep 1 en groep 2 zijn onafhankelijk.) Voorwaarden: - onafhankelijke groepen - onafhankelijke observaties in elke groep - verschil normaal verdeeld in elke groep ```{r} qqnorm(diff[GROEP == 1]) qqline(diff[GROEP == 1]) qqnorm(diff[GROEP == 2]) qqline(diff[GROEP == 2]) ``` Opnieuw is er zeer weinig data beschikbaar. Het is zeer moeilijk om met zo weinig data de normaliteitsassumptie na te gaan. Gelijke of ongelijke varianties? ```{r} var.test(diff~GROEP) ## = F-test, let op, heel gevoelig aan sample size, als deze groot is zegt de test vaak dat de variances verschillen, terwijl als de sample size klein is gaat de test de assumptie van gelijke varianties bijna nooit verwerpen ## Daarom soms beter om gewoon naar de plots te kijken boxplot(diff~GROEP) ``` Interkwartiel afstand van rechtse box is groter dan die van de linkse. Maar zolang dit niet extreem is, i.e. 3 a 4 keer groter, zullen we uitgaan van gelijke varianties. Merk op dat de var.test een objectievere maat is. Voer de test uit. ```{r} t.test(diff~GROEP, var.equal = TRUE, paired = FALSE) ## verschil tussen de verschillen (0.701 en 0.639) is niet significant ## de rivieren waar oorspronkelijk veel stikstof was hebben niet meer of minder afname van stikstof gehad dan rivieren die initieel een lage concentratie stikstof hadden ``` Interpretatie van het 95% betrouwbaarheidsinterval: - Het verschil in gemiddelde afname in nitraatwaarde in 1986 vs. 1996 tussen de twee groepen rivieren ligt met 95% waarschijnlijkheid tussen -0,36 mg N/l en 0,48 mg N/l. Conclusie: - Het verschil in gemiddelde afname in nitraatwaarde in 1986 vs. 1996 tussen de twee groepen rivieren is gelijk aan 0,062 mg N/l (95\% BI: -0,36 mg N/l tot 0,48 mg N/l). De p-waarde is gelijk aan 0,76. Er is dus geen significant verschil in de afname van de nitraatwaarde in 1986 vs. 1996 tussen beide groepen rivieren op het 5\% significantieniveau. # Voorbeeld 2 De data die we hier beschouwen zijn afkomstig van een studie waarbij men een experimentele groep van 16 ratten met de pertussis bacteria heeft geïnfecteerd. Men mat vervolgens de concentratie van de witte bloedcellen (WBC) in aantal per kubieke mm. In een controlegroep van 24 ratten werd eveneens de concentratie van witte bloedcellen gemeten. De data is beschikbaar in het bestand wbcon.dat. Dit bestand bevat twee variabelen, de variabele WBC die de `concentratie` witte bloedcellen weergeeft en de variabele `groep` die twee mogelijke waarden kan aannemen: - 0: de rat behoort tot de controlegroep - 1: de rat behoort tot de experimentele groep Lees de dataset wbcon.dat in. ```{r} wbcon <- read.table("wbcon.dat", header = TRUE) attach(wbcon) ``` ## Vraag A Is de gemiddelde concentratie witte bloedcellen groter in de controlegroep of in de experimentele groep? ```{r} mean(WBC[groep == 0]) mean(WBC[groep == 1]) ## wanneer je bacteriën toevoegd verdriedubbelt het aantal witte bloedcellen blijkbaar. Maar is dit verschil significant ``` ## Vraag B Figuur om concentratie witte bloedcellen bij de 2 groepen weer te geven. ```{r} boxplot(WBC~groep) ## variantie tussen groepen is erg verschillend (sommige ratten reageren sterk, sommigen amper) ``` ## Vraag C Welke test en welke voorwaarden om te bepalen of de gemiddelde concentratie witte bloedcellen verschillend is in controle vs experimentele groep? => **ongepaarde two-sample t-test** Veronderstellingen ongepaarde t-test: - onafhankelijke groepen - onafhankelijke observaties in elke groep - concentratie witte bloedcellen in elke groep normaal verdeeld ** testen voorwaarden ** ```{r} par(mfrow = c(1,2)) qqnorm(WBC[groep == 0]) qqline(WBC[groep == 0]) qqnorm(WBC[groep == 1]) qqline(WBC[groep == 1]) ``` - Zijn de varianties van de concentraties witte bloedcellen in de twee groepen gelijk? ```{r} var.test(WBC~groep) ## zeer significant verschil in varianties ``` Er is een extreem significant verschil in de variantie tussen de controlegroep en de experimentele groep (p < 2,2e-16) op het 5\% significantieniveau. ## Vraag D Wat zijn de nulhypotheses en alternatieve hypothese voor een tweezijdige, links eenzijdige en rechts eenzijdige test? Voer deze testen ook uit en vergelijk. Stel $\mu_1$ = gemiddelde concentratie witte bloedcellen in controlegroep en $\mu_2$ = gemiddelde concentratie witte bloedcellen in experimentele groep - tweezijdige ongepaarde t-test In symbolen: $H_0: \mu_1 = \mu_2$ In woorden: $H_0:$ de gemiddelde concentratie witte bloedcellen in controlegroep is gelijk aan de gemiddelde concentratie witte bloedcellen in experimentele groep. ## merk op dat de nulhypothese dezelfde is voor alledrie de testen, enkel de alternatieve hypothese verandert In symbolen: $H_A: \mu_1 \neq \mu_2$ In woorden: $H_A:$ de gemiddelde concentratie witte bloedcellen in controlegroep is niet gelijk aan de gemiddelde concentratie witte bloedcellen in experimentele groep. ```{r} t.test(WBC~groep) ``` - links eenzijdige ongepaarde t-test In symbolen: $H_0: \mu_1 = \mu_2$ In woorden: $H_0:$ de gemiddelde concentratie witte bloedcellen in de controlegroep is gelijk aan de gemiddelde concentratie witte bloedcellen in experimentele groep. In symbolen: $H_A: \mu_1 < \mu_2$ In woorden: $H_A:$ de gemiddelde concentratie witte bloedcellen in de controlegroep is kleiner dan de gemiddelde concentratie witte bloedcellen in experimentele groep. ```{r} t.test(WBC~groep, alternative = c("less")) ## significanter (dubbel zo significant), je kijkt maar aan 1 kant van de t-distributie ## t-test statistiek is dezelfde ``` - rechts eenzijdige ongepaarde t-test In symbolen: $H_0: \mu_1 = \mu_2$ In woorden: $H_0:$ de gemiddelde concentratie witte bloedcellen in de controlegroep is gelijk aan de gemiddelde concentratie witte bloedcellen in experimentele groep. In symbolen: $H_A: \mu_1 > \mu_2$ In woorden: $H_A:$ de gemiddelde concentratie witte bloedcellen in de controlegroep is groter dan de gemiddelde concentratie witte bloedcellen in experimentele groep. ```{r} t.test(WBC~groep, alternative = c("greater")) ``` ## Vraag E Welke van de drie testen is volgens u het meest geschikt? De tweezijdige test is het meest geschikt. Hoewel men biologisch gezien kan verwachten dat de concentratie aan witte bloedcellen zal stijgen na een infectie, is een onverwacht resultaat in de tegenovergestelde richting zeker ook interessant. Misschien onderdrukt de bacterie dan wel het immuunsysteem op een manier die nog niet gekend was, misschien zou er een andere reden zijn? In elk geval gaan we in de praktijk bij een onverwacht resultaat in de tegenovergestelde richting niet bereid zijn dit zomaar af te doen als louter toeval. Daarom testen we in de praktijk bijna altijd tweezijdig.