In het zelfstudie document en het boek is een aantal toetsen behandeld: de t-toets voor één groep, de gepaarde t-toets, en de t-toets voor twee onafhankelijke groepen. Ook is aandacht besteed aan de power van een toets en multiple testing. In dit werkcollege gaan we verder oefenen met deze onderwerpen aan de hand van R uitdraai. Aan het eind zijn er aanwijzingen voor het verkrijgen van de uitdraai in R als je dat zelf wilt oefenen.
Stappen van data analyse
In voorgaande cursussen hebben we de stappen van data analyse gebruikt. Hier pakken we die weer op en breiden we die uit met de specifieke handelingen bij een statistische toets. Ter herinnering, de vier stappen:
State. Wat is de onderzoeksvraag? Geef de praktische context erbij.
Plan. Welke specifieke statistische handelingen/berekeningen hebben we nodig? Wanneer we toetsen, noemen we hier expliciet de nul- en alternatieve hypotheses.
Solve. Gebruik geschikte visualisaties en/of berekeningen om de vraag te beantwoorden. Bij het toetsen (met betrouwbaarheidsinterval of p-waarde), checken we eerst (zo mogelijk) de aannames voor inferentie.
Conclude. Geef je conclusie in de context van de onderzoeksvraag.
Vraag 1
We hebben de gewichtsverandering in kg van 20 mannen na het gebruik van een bepaald medicijn gedurende twee maanden gemeten:
0, -6, -6, 2, -9, -4, 1, -2, -7, -6, 7, -9, 0, -4, -2, 4, -5, 0, 2, 6
De gegevens zijn in R geanalyseerd. Gebruik de uitdraai hieronder om de volgende vragen te beantwoorden.
# Data invoeren
gv <- c(0, -6, -6, 2, -9, -4, 1, -2, -7, -6, 7, -9, 0, -4, -2, 4, -5, 0, 2, 6)
# Beschrijvende statistiek
library(psych)##
## Attaching package: 'psych'## The following object is masked from 'package:SummarizedExperiment':
##
## distance## The following object is masked from 'package:qqman':
##
## manhattan## The following object is masked from 'package:randomForest':
##
## outlier## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha## The following object is masked from 'package:GenomicRanges':
##
## distance## The following objects are masked from 'package:IRanges':
##
## distance, reflectdescribe(gv, na.rm=TRUE, skew=FALSE)## vars n mean sd median min max range se
## X1 1 20 -1.9 4.71 -2 -9 7 16 1.05# Boxplot
boxplot(gv)
# t-toets
t.test(gv)##
## One Sample t-test
##
## data: gv
## t = -1.8034, df = 19, p-value = 0.08721
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -4.1051377 0.3051377
## sample estimates:
## mean of x
## -1.9a. Begin met de vier stappen. Wat is de onderzoeksvraag, en hoe wil je deze vraag beantwoorden?
Stap 1. onderzoeksvraag: verandert het gemiddelde gewicht over twee maanden bij mannen die dit medicijn gebruiken?
Stap 2. We hebben de gewichtsveranderingen van de mannen, dat is één steekproef. We zijn geïnteresseerd in de gemiddelde gewichtsverandering, dus we gebruiken een boxplot van de gewichtsveranderingen en een t-toets voor 1 steekproef.
b. Wat is de nulhypothese bij een toets om na te gaan of dit medicijn gewichtsverandering als bijwerking heeft? En wat is de alternatieve hypothese?
H0: gem. gewichtsverandering (μ)=0 en H1: μ ≠0
c. Aan welke aanname(s) moet de data voldoen om de toets uit te mogen voeren (of een betrouwbaarheidsinterval te mogen gebruiken)? Wat is je conclusie ten aanzien van die aanname(s), en hoe heb je dat gecheckt?
De gewichtsveranderingen moeten uit een normaal verdeelde populatie komen, en de 20 gewichtsveranderingen in onze steekproef moeten onafhankelijk van elkaar zijn (i.e. het moet een aselecte steekproef zijn). Als we kijken naar de boxplot van gewichtsveranderingen dan zien we dat het niet perfect, maar nagenoeg symmetrisch is. Dan is het toegestaan om de t-toets voor 1 steekproef te gebruiken.
d. Wat is de gemiddelde gewichtsverandering in de steekproef, en wat is de standaardfout van de gemiddelde gewichtsverandering?
De gemiddelde gewichtsverandering in de steekproef is -1,9. De standaardfout van de gemiddelde gewichtsverandering is 1,05.
e. Wat is het 95% BHI voor de gewichtsverandering in de populatie? Hoe werd dat berekend?
Het 95% BHI voor de gewichtsverandering in de populatie is [-4.1; 0.31]. De formule die daarvoor gebruikt werd, is \(\bar{x} \pm t^* \frac{s}{\sqrt{n}}\), waar \(\bar{x}\) de gemiddelde gewichtsverandering in de steekproef is, \(\frac{s}{\sqrt{n}}\) de standaardfout van de gemiddelde, en \(t^*\) een kritieke waarde is die opgezocht kan worden in de t-tabel bij 19 vrijheidsgraden.
f. Welke conclusie trek je op basis van het 95% BHI?
Het 95% BHI bevat (o.a.) de waarde 0, dat is het gemiddelde gewichtsverandering onder H0. We kunnen de nulhypothese niet verwerpen, we hebben onvoldoende bewijs voor een significante gewichtsverandering in mannen die dit medicijn twee maanden lang gebruiken.
g. Vind de toetsingsgrootheid (t-waarde) in de uitdraai.
De toetsingsgrootheid is -1.8034.
h. Wat betekent ‘df’ in de uitdraai, en waarom is het 19?
‘df’ staat voor ‘degrees of freedom’ (vrijheidsgraden), en met 1 steekproef is het aantal vrijheidsgraden gelijk aan n-1, in dit geval 20-1=19.
i. Zoek in de t-tabel in PSLS de kritieke t-waarde op bij α =0.05 en het juiste aantal vrijheidsgraden.
\(t^* = t_{0.05/2,19} = 2.093\).
j. Vergelijk de t- toetsingsgrootheid met de kritieke t-waarde en formuleer een conclusie.
|t| < t* dus H0 wordt niet verworpen.
k. Schat uit dezelfde t-tabel de p-waarde. Wat is hier je conclusie? Kun je ook de precieze p-waarde vinden in de R uitdraai?
Volgens de t-tabel is de kans dat een t-verdeling met 19 vrijheidsgraden een waarde heeft van -1.803 of minder, tussen de 0.05 en 0.10. Uit R krijgen we een preciezere p-waarde, namelijk 0.08721. In ieder geval is p > α, dus H0 wordt niet verworpen.
l. Trek een conclusie op basis van de p-waarde.
We hebben onvoldoende bewijs voor een significante gewichtsverandering in mannen die dit medicijn twee maanden lang gebruiken.
m. Vergelijk je conclusies uit f., i. en k.
De twee methoden (t vergelijken met t* of p vergelijken met α) moeten altijd dezelfde resultaat geven, hier in beide gevallen: H0 kan niet worden verworpen.
n. Welke toets is hier uitgevoerd? Hoe had je dit probleem anders kunnen aanpakken?
Dit is een t-toets voor 1 steekproef, uitgevoerd op verschilscores. Als we de gewichten van de mannen hadden voor en na twee maanden medicijngebruik, hadden we een t-toets voor 2 gepaarde steekproeven kunnen gebruiken.
Vraag 2
In de volgende twee opdrachten worden resultaten gebruikt van een aantal longfunctiemetingen, met name de zgn. piekstroomwaarde (PEF = ‘peak expiratory flow’ = maximale uitademingssnelheid in liter/min).
De vraag wordt gesteld of kleinere personen een andere PEF waarde hebben dan langere personen. Om dit te bepalen zijn de PEF waarden bij 24 vrouwen gemeten. Van deze groep vrouwen hadden 12 vrouwen een lengte vanaf 1,65 m en 12 vrouwen een lengte onder 1,65 m. Van elke vrouw wordt driemaal de piekstroomwaarde (PEF) (in L/min) gemeten; de gemiddelden van deze metingen zijn weergegeven in onderstaande tabel.
Piekstroomwaarden (L/min)
onder 1,65 m: 430, 415, 468, 463, 390, 401, 465, 393, 403, 440, 467, 458
vanaf 1,65 m: 497, 512, 415, 503, 499, 493, 478, 448, 465, 482, 479, 445
De gegevens zijn in R geanalyseerd.
# Data invoeren
piekstroom_onder <- c(430, 415, 468, 463, 390, 401, 465, 393, 403, 440, 467, 458)
piekstroom_vanaf <- c(497, 512, 415, 503, 499, 493, 478, 448, 465, 482, 479, 445)
# Dataframe maken
d1 <- data.frame(
piekstroom = c(piekstroom_onder, piekstroom_vanaf),
groep = c(rep(1, 12), rep(2, 12))
)
# Boxplot
boxplot(piekstroom~groep, data=d1, names=c("<1,65m",">=1,65m"))
# Beschrijvende statistiek
library(psych)
describeBy(piekstroom~groep, na.rm=TRUE, skew=FALSE, data=d1)##
## Descriptive statistics by group
## groep: 1
## vars n mean sd median min max range se
## piekstroom 1 12 432.75 31.17 435 390 468 78 9
## ------------------------------------------------------------
## groep: 2
## vars n mean sd median min max range se
## piekstroom 1 12 476.33 28.47 480.5 415 512 97 8.22t.test(piekstroom~groep, data=d1)##
## Welch Two Sample t-test
##
## data: piekstroom by groep
## t = -3.5764, df = 21.823, p-value = 0.001701
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
## -68.86815 -18.29851
## sample estimates:
## mean in group 1 mean in group 2
## 432.7500 476.3333t.test(piekstroom_onder, piekstroom_vanaf, paired=TRUE)##
## Paired t-test
##
## data: piekstroom_onder and piekstroom_vanaf
## t = -3.1686, df = 11, p-value = 0.00894
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -73.85698 -13.30968
## sample estimates:
## mean difference
## -43.58333Gebruik de juiste output uit de uitdraai hierboven en de vier stappen om de volgende vragen te beantwoorden.
a. State. Wat is de onderzoeksvraag?
Hebben vrouwen met een lengte onder 1,65 meter gemiddeld een andere piekstroomwaarde dan vrouwen met een lengte vanaf 1,65 meter?
b. Plan. Welke specifieke statistische handelingen/berekeningen hebben we nodig? Denk specifiek aan de volgende zaken:
- Zijn de steekproeven onafhankelijk of zijn de waarnemingen gepaard?
De steekproeven zijn onafhankelijk, dus we gebruiken een t-toets voor 2 onafhankelijke steekproeven.
- Hoe luiden de nulhypothese en de alternatieve hypothese die horen bij deze toets?
H0: geen verschil in gemiddelde PEF tussen beide groepen of H0: μ1-μ2=0 vs. H1: wel verschil in gemiddelde PEF tussen beide groepen of H1: μ1-μ2≠0
- Is hier sprake van een eenzijdige of tweezijdige toets?
In de formulering van het probleem staat ‘De vraag wordt gesteld of kleinere personen een andere PEF waarde hebben’, dat wordt vertaald in een tweezijdige onderzoeksvraag en dus in tweezijdige nul- en alternatieve hypothesen.
c. Solve.
- Begin bij de beschrijvende statistiek: wat zie je intuïtief? Is er op het eerste gezicht een (significant) verschil tussen beide groepen te zien?
Ja, misschien wel. In de boxplots zien we dat de uitkomsten voor vrouwen ≥1.65 m zijn in het algemeen groter zijn dan die voor vrouwen onder 1.65 m, en het gemiddelde is meer dan 40 L/min hoger bij de langere vrouwen.
- Aan welke aanname(s) moet de data voldoen om statistische inferentie te doen in deze studie? Hoe kun je de aanname(s) checken?
We nemen aan dat de twee steekproeven afkomstig zijn uit normaal verdeelde populaties. De side-by-side boxplots helpen om dit te checken. Hoewel ze niet perfect zijn, zouden allebei de boxplots afkomstig kunnen zijn uit normaal verdeelde populaties.
- Welke uitdraai, en welke delen daarvan, heb je nodig om tot een conclusie te komen?
We gebruiken de eerste toets, die de t-toets voor 2 onafhankelijke steekproeven. De tweezijdige p-waarde is 0,001701. Dit is veel kleiner dan α=0,05, dus we verwerpen de nulhypothese.
d. Conclude. Welke conclusie trek je?
We hebben sterk bewijs dat kleinere vrouwen gemiddeld een lagere piekstroomwaarde hebben dan langere vrouwen.
Vraag 3
Nadat de piekstroomwaarden voor de vrouwen in vraag 2 bepaald waren, trok een griepgolf over Nederland. Een tiental vrouwen werd ziek. Tijdens de ziekteperiode werd bij deze vrouwen nogmaals de piekstroomwaarde gemeten, men verwacht dat die lager ligt. De resultaten van deze metingen zijn weergegeven in de tabel. De waarden zijn een gemiddelde van drie metingen.
Piekstroomwaarden (in L/min)
| Gezond | Ziek |
|---|---|
| 430 | 512 |
| 468 | 463 |
| 445 | 503 |
| 499 | 478 |
| 448 | 465 |
| 415 | 489 |
| 465 | 432 |
| 444 | 505 |
| 489 | 478 |
| 450 | 434 |
De gegevens zijn in R geanalyseerd.
# Data invoeren
gezond <- c(430, 468, 445, 499, 448, 415, 465, 444, 489, 450)
ziek <- c(512, 463, 503, 478, 465, 489, 432, 505, 478, 434)
d2 <- data.frame(gezond = gezond, ziek = ziek)
d2$verschil <- d2$gezond - d2$ziek
library(psych)
describe(d2, skew=FALSE)## vars n mean sd median min max range se
## gezond 1 10 455.3 25.56 449 415 499 84 8.08
## ziek 2 10 475.9 27.95 478 432 512 80 8.84
## verschil 3 10 -20.6 43.80 -6 -82 33 115 13.85t.test(d2$gezond, d2$ziek)##
## Welch Two Sample t-test
##
## data: d2$gezond and d2$ziek
## t = -1.7198, df = 17.858, p-value = 0.1028
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -45.779693 4.579693
## sample estimates:
## mean of x mean of y
## 455.3 475.9t.test(d2$gezond, d2$ziek, paired = TRUE)##
## Paired t-test
##
## data: d2$gezond and d2$ziek
## t = -1.4874, df = 9, p-value = 0.1711
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -51.9294 10.7294
## sample estimates:
## mean difference
## -20.6t.test(d2$verschil)##
## One Sample t-test
##
## data: d2$verschil
## t = -1.4874, df = 9, p-value = 0.1711
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -51.9294 10.7294
## sample estimates:
## mean of x
## -20.6Gebruik de uitdraai hieronder en de vier stappen om de onderzoeksvraag te beantwoorden.
a. Zijn de piekstroomwaarden, zoals weergegeven in de tabel, onafhankelijk van elkaar?
Nee, de piekstroomwaarden zijn gepaard: voor elke vrouw hebben we zowel een meting wanneer ze gezond is als wanneer ze ziek is.
b. Welke toets is geschikt voor deze vraagstelling?
Een gepaarde t-toets is geschikt voor deze vraagstelling.
c. Hoe luidt de nulhypothese en de alternatieve hypothese voor deze vraagstelling?
H0: μd≤0 vs H1: μd>0 waar μd het gemiddeld verschil aangeeft tussen de gezonde piekstroomwaarde en piekstroomwaarde bij griep voor de 10 vrouwen.
d. Wat zijn de aanname(s) waarop de toets is gebaseerd? Zijn/Is daar aan voldaan?
Om de gepaarde t-toets te mogen gebruiken, moeten we eerst kijken of de verschillen uit een normale verdeling afkomstig kunnen zijn. We gebruiken dus de boxplot van de verschilscores. Hoewel die duidelijk rechts-scheef is, moeten we onthouden dat deze steekproef klein is (n=10 verschillen) en dat een steekproef van deze omvang makkelijk gegenereerd had kunnen worden uit een normaal verdeelde populatie. We gaan door met de gepaarde t-toets (tweede uitdraai). Het gemiddeld verschil in piekstroomwaarden is 11 L/min, en dit is ook precies de richting waarin je het verschil had verwacht. De toetsingsgrootheid is t=2.622, en de tweezijdige p-waarde is 0.02772, dus de éénzijdige p-waarde is 0.02772/2 = 0.01386
e. Zou je deze vraag ook met een 95% BHI kunnen beantwoorden?
Het is tricky om een éénzijdige vraag met een BHI te beantwoorden. Beter is om R te vragen een éénzijdige toets uit te voeren:
# Data invoeren
gezond <- c(430, 468, 445, 499, 448, 415, 465, 444, 489, 450)
ziek <- c(512, 463, 503, 478, 465, 489, 432, 505, 478, 434)
d2 <- data.frame(gezond = gezond, ziek = ziek)
d2$verschil <- d2$gezond - d2$ziek
library(psych)
describe(d2, skew=FALSE)## vars n mean sd median min max range se
## gezond 1 10 455.3 25.56 449 415 499 84 8.08
## ziek 2 10 475.9 27.95 478 432 512 80 8.84
## verschil 3 10 -20.6 43.80 -6 -82 33 115 13.85t.test(d2$gezond, d2$ziek, paired = TRUE, alternative="greater")##
## Paired t-test
##
## data: d2$gezond and d2$ziek
## t = -1.4874, df = 9, p-value = 0.9145
## alternative hypothesis: true mean difference is greater than 0
## 95 percent confidence interval:
## -45.98742 Inf
## sample estimates:
## mean difference
## -20.6We krijgen dan gelijk de éénzijdige p-waarde, en een éénzijdige BHI: [3.3; ∞], waar 0 (de waarde van het gemiddeld verschil onder de nulhypothese) niet in valt.
f. Er staan drie t-toetsen in de uitdraai. Verklaar de overeenkomsten en verschillen tussen de drie. Welke van de drie moet je hebben voor deze vraagstelling?
De eerste t-toets is de t-toets voor 2 onafhankelijke steekproeven, die gaat er van uit dat de metingen bij de gezonde vrouwen totaal onafhankelijk zijn van de metingen bij de zieke vrouwen, dat is hier niet het geval. De tweede toets is de t-toets voor gepaarde data, die gaat uit van afhankelijkheid: de eerste waarde onder ‘gezond’ is van dezelfde vrouw als de eerste waarde onder ‘ziek’. De derde toets is eigenlijk precies hetzelfde als de tweede: we maken een verschilscore per vrouw (PEFgezond - PEFziek) en gebruiken een t-toets voor één steekproef. Je ziet ook dat het gemiddeld verschil, de toetsingsgrootheid en de p-waarde identiek zijn voor deze twee toetsen.
g. Conclude. Welke conclusie trek je?
We hebben bewijs gevonden dat piekstroomwaarden van vrouwen lager zijn tijdens de griep dan wanneer ze gezond zijn.
Vraag 4
In de volgende vragen worden steeds twee studies uitgevoerd. Beargumenteer je antwoord op de casussen met woorden, plaatjes, of allebei.
a. In de eerste studie, A, worden twee behandelingen vergeleken in twee groepen van 30 patiënten; in de tweede studie, B, worden dezelfde twee behandelingen vergeleken in twee groepen van 20 en 30 patiënten, resp. Het verwachte effect is in beide studies gelijk. Welke studie heeft meer onderscheidingsvermogen (power)?
Studie A heeft meer onderscheidingsvermogen omdat de steekproef groter is, en dus de SEM (\(\frac{s}{\sqrt{n}}\)) kleiner zal zijn.
b. In studies C en D worden twee behandelingen vergeleken in twee groepen van 30 patiënten; in studie C was de standaarddeviatie 6 en in studie D 5. Het verwachte effect is in beide studies gelijk. Welke studie heeft meer onderscheidingsvermogen?
Studie D heeft meer onderscheidingsvermogen omdat de SD kleiner is, dus de SEM zal ook kleiner zijn.
c. In studies E en F worden twee behandelingen vergeleken in twee groepen van 30. De standaarddeviaties zijn gelijk. Het verwachte effect is groter in de eerste studie dan in de tweede. Welke studie heeft meer onderscheidingsvermogen?
In studie E; een groter effect geeft meer onderscheidingsvermogen.
d. In studies G en H worden twee behandelingen vergeleken in twee groepen. De standaarddeviaties zijn gelijk. Het verwachte effect is twee keer zo groot in studie G dan in studie H. Het onderscheidingsvermogen moet gelijk zijn. Welke studie heeft een grotere steekproef nodig?
Als studie G een groter effect verwacht dan heeft die ook meer onderscheidingsvermogen bij gelijke aantallen. Om hetzelfde onderscheidingsvermogen te bereiken, moet studie H een grotere omvang hebben.
Vraag 5
In veel biomedische artikelen is de eerste tabel een beschrijving van de steekproef. In het geval van een clinical trial worden de twee gerandomiseerde groepen beschreven. Een voorbeeld is de volgende tabel:
Bron: Guidon M & McGee H: One-year effect of a supervised exercise programme on functional capacity and quality of life in peripheral arterial disease. Disability & Rehabilitation 2013;35:39-404.
In deze tabel worden er resultaten van 21 toetsen gepresenteerd. Men kan hier meerdere kritische kanttekeningen bij zetten. Noem hier twee van.
Problemen:
Het niet verwerpen van de nulhypothese betekent niet dat de nulhypothese waar is! Dus deze toetsen geven geen bewijs dat de steekproeven gelijk zijn bij randomisatie.
Multiple testing: 21 toetsen zal de kans op een type I fout (α/alfa) flink doen oplopen! In principe verwacht je bij α=0,05 dat je gemiddelde 1 ware nulhypothese verwerpt wanneer je 20 toetsen uitvoert (1/20 = 0,05).
Het testen van baseline variabelen past niet bij een clinical trial, want we gaan er van uit dat de randomisatie ons gebalanceerde groepen geeft. De vraagstelling bij een trial gaat over het effect van de interventie, niet om verschillen in de randomisatie.