Het cursuspakket

Tijdens onderstaande uitleg wordt gebruikt gemaakt van een aantal databestanden. Download nu het cursuspakket en plaats de data bestanden in de working directory.

Voorbeeld van een “data management script”:

## Data import ##
olddata <- read.csv("mortaliteit.csv", header=TRUE, sep=",") # Behoud altijd een kopie van het origineel 
mydata <- read.csv("mortaliteit.csv", header=TRUE, sep=",")  # maak een data.frame om vanuit te werken

## Data cleanup ##
# converteer datum en tijd naar class date
mydata$gebdat <- as.Date(mydata$gebdat)
mydata$opndat <- as.Date(mydata$opndat)
mydata$ontdat <- as.Date(mydata$ontdat)
mydata$datdood <- as.Date(mydata$datdood)

## DATA MANIPULATIE ##
# Maak een nieuwe variabele leeftijd op basis van de geboortedatum en de opnamedatum
mydata$age <- round(as.numeric(mydata$opndat - mydata$gebdat)/365.25,0)

# Verwijder rijen met technisch incorrecte leeftijd outliers (behoud alle leeftijden > 1 jaar)
mydata <- subset(mydata, age > 1)

Dit script is een gecondenseerde versie van data cleanup en data manipulatie. Zonodig worden verdere transformaties van de data ten behoeve van de statistische analyse besproken.


Inleiding

De statistische analyse van de data begint met een onderzoeksvraag. Hiermee wordt een hypothese opgesteld en een bijbehorende statistische toets gekozen.

De procedure van een hypothese toetsende statistische analyse:

1. Stel op basis van de onderzoeksvraag een statistische hypothese op: de NULL-Hypothese.
2. Stel een tweede hypothese op: de alternatieve hypothese.
3. Kies een statistische significantie 'alpha' (meestal: 0.05)
4. Verricht de analyse: Verwerp de null hypothese als de p-waarde kleiner is als het gekozen statistische significantie niveau.

Bijvoorbeeld:
Onderzoeksvraag:  Is er een verschil tussen de gemiddelde leeftijd van groep A en B?
H0:               Het verschil tussen de gemiddelde leeftijd van groep A en B = 0
H1:               Het verschil in gemiddelde leeftijd van groep A en B is meer dan 0

Parametrisch

One Sample T-test

De one sample T-test is een statistische toets die de gevonden gemiddelde ( \(mu = \mu\) ) van de steekproef vergelijkt met de gekozen nullhypothese ( \(\mu_{0}\) ). Om te bepalen of deze significant verschilt wordt gebruikt gemaakt van de standaard deviatie en standard error of the mean ( \(SEM = \frac{SD}{\sqrt{n}}\) ). Hoe groter de steekproefgrootte hoe minder de steekproef gemiddelde afwijkt van de populatie gemiddelde. Dus hoe groter de steekproef, hoe kleiner de SEM en hoe betrouwbaarder de meting. Voor een significantieniveau van 5% valt grofweg te stellen dat de als de t buiten -1.96 en 1.96 valt (ookwel de Z-waarde genoemd) de null-hypothese verworpen mag worden.

Wiskundige formule: \(t=\frac{\overline{x}-u_{0}}{\frac{SD}{\sqrt{n}}}\)

t.test(x = mydata$age, mu= 80,                # nullhypothese is dat de 'echte' gemiddelde leeftijd 80 jaar is
       alternative = "two.sided",             # alternatieven: "less" of "greater" (default = two.sided)
       paired = FALSE,                        # gepaard of ongepaarde t-test (default = FALSE)
       conf.level = 0.95                      # de Confidence Interval (default = 95%)
       )
## 
##  One Sample t-test
## 
## data:  mydata$age
## t = -6.4992, df = 592, p-value = 1.715e-10
## alternative hypothesis: true mean is not equal to 80
## 95 percent confidence interval:
##  77.93363 78.89268
## sample estimates:
## mean of x 
##  78.41315
# De gemiddelde leeftijd van de steekproef is 78.4 met een 95% CI van ~78 tot ~79. De Null-hypothese dat de gemiddelde leeftijd 80 jaar is kan met een 95% betrouwbaarheid verworpen worden.

Two sample T-test

De two-sample T-test wordt gebruikt om de null-hypothese \(\mu_{0}\) te toetsen dat beide steekproeven afkomstig zijn uit dezelfde populatie. Er zijn twee methodes voor het testen van gemiddelden: met en zonder aanname van homogeniteit van variantie. De equal versie neemt aan dat de spreiding (standaard deviatie) tussen de groepen niet teveel verschilt. Dit kan worden getest met de Levine test. De variant van equal variance heeft in principe de voorkeur, het geeft een nauwkeurigere p-waarde.

Wiskundige formules: Equal variance \(t=\frac{\overline{x}_{1}-\overline{x}_{2}}{\sqrt{\frac{(n_1-1)(SD_1)^2+(n_2-1)(SD_2)^2}{n_1+n_2-2}}}\) en Unequal variance \(t=\frac{\overline{x}_{1}-\overline{x}_{2}}{\sqrt{\frac{(SD_1)^2}{n_1}+\frac{(SD_2)^2}{n_2}}}\)

# Test de homogeniteit van de variantie
library(car)
leveneTest(mydata$age,mydata$delier)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value  Pr(>F)  
## group   1  5.9115 0.01534 *
##       591                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# De null-hypothese stelt dat het verschil van variantie in de groepen door 'random sampling' komt. de p-waarde is SIGNIFICANT dus de null-hypothese moet verworpen worden. Er mag GEEN gebruik worden gemaakt van de aanname van "equal variance". 

# nullhypothese is dat de gemiddelde leeftijd voor mensen met en zonder een delier hetzelfde zijn.
t.test(mydata$age~mydata$delier,              # age = uitkomstvariabele. delier = groepen.
       var.equal = FALSE                      # keuze van type formule (default = FALSE)
       )
## 
##  Welch Two Sample t-test
## 
## data:  mydata$age by mydata$delier
## t = 5.8628, df = 84.955, p-value = 8.455e-08
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  3.136942 6.356501
## sample estimates:
##      mean in group delier mean in group geen delier 
##                  82.59155                  77.84483
# De gemiddelde leeftijd en distributie van patiënten met en zonder een delier zijn significant verschillend van elkaar.

Paired T-test

De gepaarde t-test is een vorm van de de t-test waarbij twee metingen van dezelfde variabele met elkaar vergeleken worden.

## Voorbeeld:
meting_1 <- c(120,140,147,130,125,110,170,100,99)     # systolische bloeddruk VOOR het drinken van een energydrank
meting_2 <- c(130,145,168,120,140,120,150,110,120)    # systolische bloeddruk  NA  het drinken van een energydrank

# nullhypothese is dat het verschil in gemiddelde bloeddruk van de 1ste en 2de meting NIET significant verschillen
t.test(meting_1,meting_2, paired=TRUE)  # gebruik de parameter "paired" om aan te geven dat het gaat om een GEPAARDE t-toets. 
## 
##  Paired t-test
## 
## data:  meting_1 and meting_2
## t = -1.5068, df = 8, p-value = 0.1703
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -17.43145   3.65367
## sample estimates:
## mean of the differences 
##               -6.888889
# de null-hypothese is dat het verschil tussen de gemiddelden gelijk is aan 0.
# de null-hypothese mag NIET verworpen worden want: p > 0.05 EN/OF de "0" valt in de 95% Conference Interval

ANOVA

ANOVA staat voor Analysis of variance. Het is een statistische test door de analyse van de variantie \(\sigma\). De spreiding rondom het gemiddelde kan per groep (within) of tussen groepen (between) verschillen. ANOVA test of de spreiding tussen de groepen (\(\mu_{1}\), \(\mu_{2}\), \(\mu_{3}\), etc…) verschilt van de spreiding binnen een groep. Als de spreiding tussen de groepen significant groter is dan de spreiding binnen de groepen moet de null-hypothese worden verworpen. Voor een video uitleg klik hier.

In wiskundige termen: \(F = \frac{between\:group\:variance}{within\:group\:variance}\) wat neerkomt op de de ratio tussen de verklaarbare systemiche variantie en de onverklaarbare variantie. Het verwachte verschil in gemiddelden gedeeld door het gevonden verschil is de F-ratio (Fisher). De waarschijnlijkheid (probability) dat de gevonden verschil significant afwijkt van het verwachte verschil wordt met een significantie-niveau ‘alpha’ (de p-waarde) vastgesteld.

F distributie

Voorwaarden:

  1. ANOVA is een parametrische test: sensitief voor non-normaliteit. Indien niet voldaan wordt aan de aanname van normaliteit kan gekozen worden voor een non-parametrische alternatief: Kruskal-Wallis-test.
  2. De variantie van de uitkomstvariabele (DV) is gelijk in alle groepen (IV): homogeniteit van variantie. Doe een a priori test.

ANOVA formule

onderzoeksvorm ANOVA vorm Formule

DV (ongepaard) ~ IV

One-way ANOVA

DV ~ IV

DV (gepaard) ~ IV

Repeated measures ANOVA

DV ~ IV + Error(N/IV)

Overige ANOVA varianten (worden niet besproken):

onderzoeksvorm ANOVA vorm Formule

onderzoeken van covariaat

One-way ANCOVA met één covariaat (x)

DV ~ IV + x

onderzoeken van 2 IV

Two-way Factorial ANOVA

DV ~ IV1 * IV2

2 IV + 2 covariaten

Two-way ANCOVA met twee covariaten (x1,x2)

DV ~ IV1 * IV2 + x1 + x2


A priori testen

Om te testen voor de homogeniteit van de variantie kan gebruik worden gemaakt van de Levine, Brown-Forsythe en Bartlett testen. Lees meer.

Testen voor homogeniteit van de variantie: H0 = homogeniteit van de variantie

  • Bartlett’s test: testen van de homogeniteit van de variantie
  • Levene’s test: testen van de variantie middels de gemiddelden
  • Brown-Forsythe test: testen van de variantie middels de mediaan
# A priori analyse met de Bartlett test
bartlett.test(visus ~ medicijn, data=mydata)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  visus by medicijn
## Bartlett's K-squared = 3.4401, df = 2, p-value = 0.1791
# A priori analyse met de levene's Test (gemiddelde)
library(car)
leveneTest(visus ~ medicijn, data=mydata)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   2    1.28 0.2788
##       590
# A priori analyse met de Brown-€“Forsythe test (mediaan)
library(HH)
hov(visus ~ medicijn, data=mydata)
## 
##  hov: Brown-Forsyth
## 
## data:  visus
## F = 1.28, df:medicijn = 2, df:Residuals = 590, p-value = 0.2788
## alternative hypothesis: variances are not identical

Wat te doen als er niet wordt voldaan aan de H0: homogeniteit van de variantie:

Post-hoc testen

Met ANOVA wordt alleen getest of er een verschil is tussen de groepen, niet waar de de verschillen vandaan komen. Om dit te testen kan een post-hoc analyse worden verricht. Er zijn verschillende type post-hoc analyses beschikbaar afhankelijk van de onderzoekssituatie. Lees meer over post-hoc analyses.

Een aspect bij het uitvoeren van post-hoc testen is het probleem van Kanskapitalisme, waarbij de kans op een significant verschil dat op toeval berust (significantieniveau ‘alpha’) stijgt. Om hiervoor te corrigeren zijn een aantal correctiemethoden beschikbaar. Hieronder volgt een uitleg en bij de ‘one-way anova’ en de ‘repeated measures anova’ wordt ook een voorbeeld van de code in R gegeven.

Pairwise t-test met een Bonferroni correctie: Deze methode is gebaseerd op het idee dat de onderzoeker k hypothesen post-hoc (doorgaans na ANOVA) toetst. Om het oorspronkelijke significantieniveau ‘alpha’ te behouden moet volgens de Bonferroni-correctie ‘alpha’ verkleind worden tot ‘alpha’/k. Het nieuwe significantieniveau voor hypothesetoetsing is dan 1/k van de oorspronkelijke ‘alpha’. Lees meer

Tukey’s test: Vergelijken van alle factoren binnen een groep (a vs b | b vs c | a vs c) zonder te kijken naar combinatie’s van factoren (NIET; a vs b+c)

Scheffe’s Test: Vergelijken van alle mogelijke combinatie’s van groepen

Kiezen van de beste test:

  • Pairwise t-test met een Bonferroni correctie is meestal te strikt wanneer de hoeveelheid groepen om met elkaar te vergelijken (~ > 6) te groot wordt. Dan is deze test niet meer aanbevolen.
  • Tukey’s test: is stricter dan de bonferroni correctie bij kleine groepen maar minder strict als de groepen te groot worden voor een bonferroni correctie.
  • Scheffe’s test: Vergelijkt alle mogelijke combinaties en is daardoor de meest stricte test.

One way ANOVA

De one way anova is bedoel voor 1 uitkomstvariabele (DV) en 1 onafhankelijke variabele (IV) met >2 factor levels. Bijvoorbeeld: DV = visus-score (continue/interval) en IV = medicijn (factor met 3 levels)

Dit kan op verschillende manieren in R:

  • oneway.test(DV ~ IV, data=dataframe, var.equal=TRUE/FALSE)
  • aov(DV ~ IV, data=dataframe)

Het kan ook met de anova(lm(DV~IV)) combinatie functie maar dit is niet aanbevolen voor een one-way anova, wel voor andere vormen.

De oneway.test() functie uit de core-package stat:

# De *oneway.test()* verricht automatisch een Welsh correctie (var.equal=FALSE). Gebruik var.equal=TRUE om dit niet te doen. Daarnaast negeert het alle missing values in de analyse (zonder waarschuwing)!

# One Way ANOVA met de oneway.test() functie 
owANOVA <- oneway.test(visus ~ medicijn, data=mydata, # ANOVA van de visus (DV) per medicijn groep (IV)
                       var.equal=FALSE)               # var.equaL=FALSE want er wordt niet voldaan aan de H0: homogeniteit van de variantie

owANOVA
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  visus and medicijn
## F = 11.924, num df = 2.00, denom df = 390.09, p-value = 9.41e-06
# Conclusie: er zit een significante verschil in de variantie van de gemiddelde leeftijd (DV) per medicijn groep (IV)

De aov() functie uit de package car:

# laad de package 'car'
library(car)

# One Way ANOVA met de aov() functie 
ANOVA = aov(visus ~ medicijn, data=mydata) # ANOVA van de visus (DV) per medicijn groep (IV)

# geef de summary van ANOVA
summary(ANOVA)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## medicijn      2  0.518 0.25891   10.85 2.36e-05 ***
## Residuals   590 14.082 0.02387                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Er is statistisch significant een verschil in de gemiddelde visus en variantie per medicijngroep. Met deze functie is niet duidelijk waar dit verschil vandaan komt. Hiervoor is een post-hoc analyse nodig.
## Post Hoc Analyse
# Vergelijk de gemiddelde visus-scores van de verschillende medicijn groepen
model.tables(ANOVA, type="means")
## Tables of means
## Grand mean
##           
## 0.4174283 
## 
##  medicijn 
##     geen med.   haldol placebo
##        0.4635   0.3961   0.401
## rep  172.0000 208.0000 213.000
# Doe een kruistabel t-test met correctie voor multiple testing volgens 'Holm'
pairwise.t.test(mydata$visus, mydata$medicijn )
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  mydata$visus and mydata$medicijn 
## 
##         geen med. haldol 
## haldol  7.9e-05   -      
## placebo 0.00018   0.74484
## 
## P value adjustment method: holm
# Doe een kruistabel t-test met correctie voor multiple testing volgens 'Bonferonni'
pairwise.t.test(mydata$visus, mydata$medicijn, p.adj="bonferroni" )
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  mydata$visus and mydata$medicijn 
## 
##         geen med. haldol 
## haldol  7.9e-05   -      
## placebo 0.00026   1.00000
## 
## P value adjustment method: bonferroni
# Verricht een POST-HOC analyse van de ANOVA volgens Scheffe
library(DescTools)                  
ScheffeTest(ANOVA, which="medicijn") # medicijn (IV)
## 
##   Posthoc multiple comparisons of means : Scheffe Test 
##     95% family-wise confidence level
## 
## $medicijn
##                          diff      lwr.ci      upr.ci    pval    
## haldol-geen med.  -0.06744074 -0.10651358 -0.02836790 0.00015 ***
## placebo-geen med. -0.06253712 -0.10140183 -0.02367241 0.00046 ***
## placebo-haldol     0.00490362 -0.03205355  0.04186079 0.94838    
## 
## ---
## Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
# Verricht een POST-HOC analyse van de ANOVA volgens Tukey
posthoc <- TukeyHSD(ANOVA)
posthoc
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = visus ~ medicijn, data = mydata)
## 
## $medicijn
##                          diff         lwr         upr     p adj
## haldol-geen med.  -0.06744074 -0.10485235 -0.03002913 0.0000784
## placebo-geen med. -0.06253712 -0.09974945 -0.02532479 0.0002594
## placebo-haldol     0.00490362 -0.03048227  0.04028951 0.9432429
plot(posthoc)

Lees meer over Oneway Anova


Repeated measures ANOVA

Een repeated measures ANOVA (within-subject design) wordt gebruikt voor het testen van een verschil in herhaalde metingen bij dezelfde patiënten. In onderstaand voorbeeld ligt de focus op één groep en drie metingen. Bijvoorbeeld: het meten van de bloeddruk voor en na het drinken van de een energydrank (within group) voor 3 risicogroepen (between group met 3 factor levels). LET OP: het is wel noodzakelijk eerst de data mogelijk te transformeren van BREED naar LANG (van een kolom per meting naar alle metingen in 1 kolom).

Bekijk eerst deze Video instructie

## Voorbeeld:
# Een dataframe van 9 patiënten met 3 metingen van de bloeddruk
ptid <- c(1,2,3,4,5,6,7,8,9)
meting_1 <- c(120,140,147,130,125,110,170,100,99)     # systolische bloeddruk VOOR het drinken van een energydrank
meting_2 <- c(130,145,168,120,140,120,150,110,120)    # systolische bloeddruk  NA  het drinken van een energydrank
meting_3 <- c(160,135,180,131,148,115,165,120,130)    # systolische bloeddruk  NA  het drinken van een 2de energydrank
data <- data.frame(ptid,meting_1,meting_2,meting_3)
head(data)
##   ptid meting_1 meting_2 meting_3
## 1    1      120      130      160
## 2    2      140      145      135
## 3    3      147      168      180
## 4    4      130      120      131
## 5    5      125      140      148
## 6    6      110      120      115
## Data Manipulatie
# Om een Repeated Measures Anova uit te voeren moet de data.frame eerst worden hergestructureerd van een BREED naar LANG
library(reshape)
newdata <- melt(data, id.vars = c("ptid"))
head(newdata,3); tail(newdata,3)
##   ptid variable value
## 1    1 meting_1   120
## 2    2 meting_1   140
## 3    3 meting_1   147
##    ptid variable value
## 25    7 meting_3   165
## 26    8 meting_3   120
## 27    9 meting_3   130
## a priori testen voor homogeniteit van de variantie
# testen voor homogeniteit in de variantie
bartlett.test(value ~ variable, data=newdata)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  value by variable
## Bartlett's K-squared = 0.3915, df = 2, p-value = 0.8222
## ANOVA
# Verricht in repeated measures anova en sla het resultaat op in een nieuw object
rmANOVA = aov(value ~ variable + Error(ptid/variable), data=newdata)    # value = de metingen (DV),  variable = de groepen (IV)

# gebruik de summary() functie voor de analyse van de rmANOVA
summary(rmANOVA)
## 
## Error: ptid
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  1   1805    1805               
## 
## Error: ptid:variable
##          Df Sum Sq Mean Sq
## variable  2  818.7   409.4
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## variable   2    335   167.6   0.388  0.683
## Residuals 21   9069   431.9
## Post-Hoc
# een post-hoc analyse is niet mogelijk met een aov() functie. Gebruik hiervoor de tegenhanger lme() uit de package 'nlme'
# lme() methode:
library(nlme)
Lme.mod <- lme(value ~  variable, random = ~1 | ptid/value, data = newdata)
summary(Lme.mod)
## Linear mixed-effects model fit by REML
##  Data: newdata 
##        AIC      BIC    logLik
##   217.2508 224.3191 -102.6254
## 
## Random effects:
##  Formula: ~1 | ptid
##         (Intercept)
## StdDev:    18.72684
## 
##  Formula: ~1 | value %in% ptid
##         (Intercept) Residual
## StdDev:     8.97015  4.73358
## 
## Fixed effects: value ~ variable 
##                      Value Std.Error DF   t-value p-value
## (Intercept)      126.77778  7.099020 16 17.858489  0.0000
## variablemeting_2   6.88889  4.781222 16  1.440822  0.1689
## variablemeting_3  15.88889  4.781222 16  3.323186  0.0043
##  Correlation: 
##                  (Intr) vrbl_2
## variablemeting_2 -0.337       
## variablemeting_3 -0.337  0.500
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -0.58875604 -0.29245297 -0.02148609  0.23429771  0.84470543 
## 
## Number of Observations: 27
## Number of Groups: 
##            ptid value %in% ptid 
##               9              27
# gebruik de 'multcomp' package voor de posthoc analyse
library(multcomp)
summary(glht(Lme.mod, linfct=mcp(variable = "Tukey")), test = adjusted(type = "bonferroni"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lme.formula(fixed = value ~ variable, data = newdata, random = ~1 | 
##     ptid/value)
## 
## Linear Hypotheses:
##                          Estimate Std. Error z value Pr(>|z|)   
## meting_2 - meting_1 == 0    6.889      4.781   1.441  0.44891   
## meting_3 - meting_1 == 0   15.889      4.781   3.323  0.00267 **
## meting_3 - meting_2 == 0    9.000      4.781   1.882  0.17936   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- bonferroni method)

Bekijk video instructie + bijbehorende tutorial.

Lees meer over repeated measures ANOVA met R. Zie ook Laerd Statistics voor een introductie in ‘Repeated Measures Anova’. Lees meer over multcomp package.


Non-parametrisch

One-Sample Signed Rank Test

Dit is de non-parametrische variant van de ‘One sample T-test’.

wilcox.test(mydata$age, mu=78.5)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  mydata$age
## V = 80730, p-value = 0.07867
## alternative hypothesis: true location is not equal to 78.5

Paired-Sample Wilcoxon Signed Rank Test

Dit is de non-parametrische variant van de gepaarde ‘Two sample T-test’.

## Voorbeeld:
meting_1 <- c(120,140,147,130,125,110,170,100,99)     # systolische bloeddruk VOOR het drinken van een energydrank
meting_2 <- c(130,145,168,120,140,120,150,110,120)    # systolische bloeddruk  NA  het drinken van een energydrank

# Nullhypothese is dat het verschil in mediane bloeddruk van de 1ste en 2de meting NIET significant verschillen
wilcox.test(meting_1,meting_2, paired=TRUE)  # gebruik de parameter "paired" om aan te geven dat het gaat om een GEPAARDE t-toets. 
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  meting_1 and meting_2
## V = 10.5, p-value = 0.1689
## alternative hypothesis: true location shift is not equal to 0
# De null-hypothese is dat het verschil tussen de distributie rondom de mediaan gelijk is aan 0.
# De null-hypothese mag NIET verworpen worden want: p > 0.05

Mann-Whitney-U

Dit is de non-parametrische variant van de ongepaarde ‘Two sample T-test’ en wordt ookwel de Mann-Whitney-Wilcoxon Test genoemd.

wilcox.test(mydata$age~mydata$delier)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  mydata$age by mydata$delier
## W = 26233.5, p-value = 1.229e-08
## alternative hypothesis: true location shift is not equal to 0

Kruskal-Wallis toets

Kruskal Wallis is de non-parametrische variant van de ‘One-Way Anova’.

# Kruskal-Wallis toets
kruskal.test(age ~ medicijn, data = mydata) 
## 
##  Kruskal-Wallis rank sum test
## 
## data:  age by medicijn
## Kruskal-Wallis chi-squared = 61.228, df = 2, p-value = 5.064e-14
library(PMCMR)
posthoc.kruskal.nemenyi.test(x=mydata$age, g=mydata$medicijn, method="Chisq")
## 
##  Pairwise comparisons using Tukey and Kramer (Nemenyi) test  
##                    with Tukey-Dist approximation for independent samples 
## 
## data:  mydata$age and mydata$medicijn 
## 
##         geen med. haldol
## haldol  4.8e-09   -     
## placebo 2.7e-13   0.31  
## 
## P value adjustment method: none
# Post hoc 
library(pgirmess)
kruskalmc(mydata$age, mydata$medicijn)
## Multiple comparison test after Kruskal-Wallis 
## p.value: 0.05 
## Comparisons
##                    obs.dif critical.dif difference
## geen med.-haldol  106.5675     42.27137       TRUE
## geen med.-placebo 131.0720     42.04620       TRUE
## haldol-placebo     24.5045     39.98251      FALSE
# Gepaarde wilcoxon test met correctie voor multiple testing volgens 'bonferroni'
pairwise.wilcox.test(mydata$age, mydata$medicijn, p.adj="bonferroni", exact=F)
## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  mydata$age and mydata$medicijn 
## 
##         geen med. haldol
## haldol  6.7e-09   -     
## placebo 1.4e-13   0.47  
## 
## P value adjustment method: bonferroni

Friedman toets

Friedman is de non-parametrische variant van de ‘Repeated Measures Anova’.

## Voorbeeld:
# Een dataframe van 9 patiënten met 3 metingen van de bloeddruk
ptid <- c(1,2,3,4,5,6,7,8,9)
meting_1 <- c(120,140,147,130,125,110,170,100,99)     # systolische bloeddruk VOOR het drinken van een energydrank
meting_2 <- c(130,145,168,120,140,120,150,110,120)    # systolische bloeddruk  NA  het drinken van een energydrank
meting_3 <- c(160,135,180,131,148,115,165,120,130)    # systolische bloeddruk  NA  het drinken van een 2de energydrank
data <- data.frame(ptid,meting_1,meting_2,meting_3)
head(data)
##   ptid meting_1 meting_2 meting_3
## 1    1      120      130      160
## 2    2      140      145      135
## 3    3      147      168      180
## 4    4      130      120      131
## 5    5      125      140      148
## 6    6      110      120      115
## Data Manipulatie
# Om een Repeated Measures Anova uit te voeren moet de data.frame eerst worden hergestructureerd van een BREED naar LANG
# Zie Data Management -> Data Manipulatie -> Reshape voor meer informatie
library(reshape)
newdata <- melt(data, id.vars = c("ptid"))
head(newdata,3); tail(newdata,3)
##   ptid variable value
## 1    1 meting_1   120
## 2    2 meting_1   140
## 3    3 meting_1   147
##    ptid variable value
## 25    7 meting_3   165
## 26    8 meting_3   120
## 27    9 meting_3   130
# Friedman toets kan op twee manieren:
friedman.test(newdata$value, newdata$variable, newdata$ptid) 
## 
##  Friedman rank sum test
## 
## data:  newdata$value, newdata$variable and newdata$ptid
## Friedman chi-squared = 5.5556, df = 2, p-value = 0.06218
friedman.test(value ~ variable | ptid, data= newdata)
## 
##  Friedman rank sum test
## 
## data:  value and variable and ptid
## Friedman chi-squared = 5.5556, df = 2, p-value = 0.06218
## Post hoc
# De post-hoc analyse van de Friedman-toets is het herhaaldelijk uitvoeren van een Mann-Whitney-U toets. 
# Dit is de non-parametrische equivalent van een 'pairwise t-test' (gebruik paired=TRUE.

# Correctie volgens 'bonferroni'
pairwise.wilcox.test(newdata$value, newdata$variable, p.adj="bonferroni", exact=FALSE, paired=TRUE)
## 
##  Pairwise comparisons using Wilcoxon signed rank test 
## 
## data:  newdata$value and newdata$variable 
## 
##          meting_1 meting_2
## meting_2 0.51     -       
## meting_3 0.17     0.13    
## 
## P value adjustment method: bonferroni
# Correctie volgens 'Holm'
pairwise.wilcox.test(newdata$value, newdata$variable, p.adj="holm", exact=FALSE, paired=TRUE)
## 
##  Pairwise comparisons using Wilcoxon signed rank test 
## 
## data:  newdata$value and newdata$variable 
## 
##          meting_1 meting_2
## meting_2 0.17     -       
## meting_3 0.13     0.13    
## 
## P value adjustment method: holm
## Post hoc met de package PMCMR
library(PMCMR)
posthoc.friedman.nemenyi.test(newdata$value, newdata$variable, newdata$ptid) 
## 
##  Pairwise comparisons using Nemenyi post-hoc test with q approximation for unreplicated blocked data 
## 
## data:  newdata$value , newdata$variable and newdata$ptid 
## 
##          meting_1 meting_2
## meting_2 0.466    -       
## meting_3 0.048    0.466   
## 
## P value adjustment method: none

Lees meer over Friedman toets.

Categoriale uitkomstvariabele

Fishers Exact (2 x 2 kruistabel)

De Fisher’s exact toets is een parametrische test die geschikt is voor kleine patiëntpopulaties. Grofweg is een Fisher’s exact toets te verkiezen boven een Chi-square bij een subgroeppopulatie (N per cel) onder de tien.

# Tussenstap: maak een kruistabel
kruistabel <- table(mydata$dood, mydata$delier)

## Fishers Exact test
# Gebruik conf.int om de confidence interval te berekenen. 
# Gebruik conf.level om de confidence level (default = 0.95) in te stellen.
fisher.test(kruistabel, conf.int = TRUE, conf.level = 0.99) 
## 
##  Fisher's Exact Test for Count Data
## 
## data:  kruistabel
## p-value = 1.103e-08
## alternative hypothesis: true odds ratio is not equal to 1
## 99 percent confidence interval:
##  0.1088797 0.4586787
## sample estimates:
## odds ratio 
##   0.227152
# Post Hoc
library(vcd) 
assocstats(kruistabel)
##                     X^2 df   P(> X^2)
## Likelihood Ratio 32.902  1 9.6918e-09
## Pearson          35.408  1 2.6739e-09
## 
## Phi-Coefficient   : 0.244 
## Contingency Coeff.: 0.237 
## Cramer's V        : 0.244
## Bereken ODDS-RATIO en bijbehorende confidence interval
# Gebruik package epiR
library(epiR)
epi.2by2(kruistabel, method="cohort.count", conf.level=0.95)
##              Outcome +    Outcome -      Total        Inc risk *
## Exposed +           26          375        401              6.48
## Exposed -           45          147        192             23.44
## Total               71          522        593             11.97
##                  Odds
## Exposed +      0.0693
## Exposed -      0.3061
## Total          0.1360
## 
## Point estimates and 95 % CIs:
## ---------------------------------------------------------
## Inc risk ratio                               0.28 (0.18, 0.43)
## Odds ratio                                   0.23 (0.13, 0.39)
## Attrib risk *                                -16.95 (-23.41, -10.50)
## Attrib risk in population *                  -11.46 (-18.00, -4.93)
## Attrib fraction in exposed (%)               -261.48 (-467.56, -130.23)
## Attrib fraction in population (%)            -95.75 (-132.86, -64.56)
## ---------------------------------------------------------
##  * Cases per 100 population units
# Geef de Confidence interval van de oddsratio per groep
epi.conf(kruistabel, ctype = "odds", method = "fleiss", conf.level = 0.95)
##                       est      lower      upper
## niet overleden 0.06933333 0.04427083 0.09863014
## overleden      0.30612245 0.21518987 0.42222222

LET OP: Bij sterk geskewde data is het beter om een Barnard’s Test te gebruiken.

Chi-squared (2 x 2+ kruistabel)

# Tussenstap: maak een kruistabel
kruistabel <- table(mydata$medicijn, mydata$delier)
kruistabel
##            
##             delier geen delier
##   geen med.      6         166
##   haldol        31         177
##   placebo       34         179
# Visuele inspectie van een relatieve kruistabel
barplot(prop.table(kruistabel), beside=TRUE, col = c(1,2,3))
legend("topleft", inset=.05, title="Legenda", levels(mydata$medicijn), fill = c(1,2,3))

# Chi-square test
CHI <- chisq.test(kruistabel, correct=T)  # verricht een Yates correctie
CHI
## 
##  Pearson's Chi-squared test
## 
## data:  kruistabel
## X-squared = 16.66, df = 2, p-value = 0.0002412
# Chi-square prediction tabel
CHI$expected
##            
##               delier geen delier
##   geen med. 20.59359    151.4064
##   haldol    24.90388    183.0961
##   placebo   25.50253    187.4975
# Confidence Interval van de 'prediction tabel'
epi.conf(CHI$expected, ctype = "odds", method = "fleiss", conf.level = 0.95)
##                 est      lower     upper
## geen med. 0.1360153 0.08176101 0.2027972
## haldol    0.1360153 0.08333333 0.1954023
## placebo   0.1360153 0.08673469 0.1966292
# Post Hoc analyse voor de coefficient
library(vcd) 
assocstats(kruistabel)
##                     X^2 df   P(> X^2)
## Likelihood Ratio 20.293  2 3.9222e-05
## Pearson          16.660  2 2.4117e-04
## 
## Phi-Coefficient   : 0.168 
## Contingency Coeff.: 0.165 
## Cramer's V        : 0.168

McNemar-Bowker

McNemar’s test is een gepaarde versie van de chi-square toets, vergelijkbaar met de gepaarde t-test.

Bijvoorbeeld: Is er een verschil in de het klinische aspect van de patiënt (0 = niet ziek, 1 = ziek) per meetmoment (1 = bij binnenkomst, 2 = na 24 uur).

# Een voorbeeld dataframe
meting_1 <- c(1,1,1,1,1,1,1,0,1,1)                     # 0 = niet ziek, 1 = ziek
meting_2 <- c(0,0,1,0,1,0,1,1,1,0)                     # 0 = niet ziek, 1 = ziek
df <- data.frame(meting_1,meting_2)

# Tussenstap: maak een kruistabel
table <- table(df$meting_1,df$meting_2)

# Verricht de mcnemar test
mcnemar.test(table)
## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  table
## McNemar's chi-squared = 1.5, df = 1, p-value = 0.2207

Cochran’s Q toets

Deze test onderzoekt of de ratio van een binomiale uitkomstvariabele (DV) verschilt per (>2) metingen (IV). Het is de variant van de repeated measures anova en friedman test voor binomiale gepaarde data Of anders verwoord de gepaarde variant van de chi-square. Post Hoc kan een mcnemar-test worden uitgevoerd.

Bijvoorbeeld: Patiënten ogen ‘ziek’ of ‘niet ziek’. Is er een verschil in de het klinische aspect van de patiënt (0 = niet ziek, 1 = ziek) per meetmoment (1 = bij binnenkomst, 2 = na 24 uur, 3 = na 48 uur).

# Een voorbeeld dataframe
ptid <- seq(1:10)                                  # rij getallen van 1 t/m 10
meting_1 <- c(1,1,0,1,0,1,1,0,1,1)                 # 0 = niet ziek, 1 = ziek
meting_2 <- c(1,1,0,0,1,0,1,0,1,1)                 # 0 = niet ziek, 1 = ziek
meting_3 <- c(0,0,0,1,0,0,0,0,0,0)                 # 0 = niet ziek, 1 = ziek
df <- data.frame(ptid,meting_1,meting_2,meting_3)

# tussenstap: transposeer met functie melt()
newdata <- melt(df, id.vars = c("ptid")  )

# Cochran's Q test
library(coin) 
symmetry_test(value ~ factor(variable) | factor(ptid), data = newdata, teststat = "quad")
## 
##  Asymptotic General Independence Test
## 
## data:  value by
##   factor(variable) (meting_1, meting_2, meting_3) 
##   stratified by factor(ptid)
## chi-squared = 7.75, df = 2, p-value = 0.02075
## Post-hoc analyse middels multiple mcnemar test
# meting 1 vs meting 2
p1 <- mcnemar.test(newdata[newdata$variable=="meting_1",]$value, newdata[newdata$variable=="meting_2",]$value)
p1
## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  newdata[newdata$variable == "meting_1", ]$value and newdata[newdata$variable == "meting_2", ]$value
## McNemar's chi-squared = 0, df = 1, p-value = 1
# meting 2 vs meting 3
p2 <- mcnemar.test(newdata[newdata$variable=="meting_2",]$value, newdata[newdata$variable=="meting_3",]$value)
p2
## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  newdata[newdata$variable == "meting_2", ]$value and newdata[newdata$variable == "meting_3", ]$value
## McNemar's chi-squared = 2.2857, df = 1, p-value = 0.1306
# meting 1 vs meting 3
p3 <- mcnemar.test(newdata[newdata$variable=="meting_1",]$value, newdata[newdata$variable=="meting_3",]$value)
p3
## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  newdata[newdata$variable == "meting_1", ]$value and newdata[newdata$variable == "meting_3", ]$value
## McNemar's chi-squared = 4.1667, df = 1, p-value = 0.04123
## correctie volgens bonferonni
# extraheer alle p waardes
p <- c(p1$p.value,p3$p.value,p2$p.value)
# correctie toepassen
p.adjust(p, method="bonferroni")
## [1] 1.0000000 0.1236805 0.3917101
# de 2de analyse: "meting 2 vs meting 3" is niet significant meer na een bonferroni correctie

Creative Commons-Licentie Dit werk valt onder een Creative Commons Naamsvermelding - Gelijk Delen 3.0 Internationaal-licentie .