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 het statistisch beschrijven van de data besproken.


Kiezen statistische test

Voor het kiezen van de correcte statistische test moeten er 3 vragen beantwoord woorden:

  • Is de uitkomstvariabele continue of categoriaal?
  • Is de uitkomstvariabele normaal verdeeld? (gaat alleen op voor continue data)
  • Zijn de metingen gepaard?

Parametrisch vs Non-parametrisch

De keuze voor een parametrische test heeft in principe altijd de voorkeur boven een non-parametrische alternatief. Echter moet er dan wel voldaan worden aan de aanname dat de data normaal verdeeld is. Indien niet aan deze voorwaarde wordt voldaan kan gepoogd worden een transformatie uit te voeren, bijvoorbeeld een logistische transformatie en kan indien de data dan wel normaal verdeeld is alsnog voor een parametrische test worden gekozen.

De analyse van continue uitkomstvariabelen

VERGELIJKINGEN

Onderzoekssituatie

Parametrisch

Non parametrisch

vergelijken met gemiddelde

one sample t-toets

One-Sample Signed Rank Test

vergelijken binnen 1 groep

gepaarde t-toets

Wilcoxon Signed Rank

vergelijken 2 groepen

ongepaarde t-toets

Mann-Whitney-U

2+ ongepaarde groepen

one way ANOVA

Kruskal-Wallis-toets

2+ gepaarde groepen

Repeated measures ANOVA

Friedman toets

ASSOCIATIES

Onderzoekssituatie

Parametrisch

Non parametrisch

crosssectionele relatie

Pearson (regressie co?fficient)
en lineare regressie analyse

Spearman-correlatie
en lineare regressie analyse

longitudinale relatie

Regressie analyse

Regressie analyse

De analyse van categoriale uitkomstvariabelen

Categoriaal

Onderzoekssituatie

Statistische Techniek

1 variabele vs referentie

Z-test voor proporties

2 variabelen (gepaard)

McNemar-Bowker of Cochran’??s Q toets

2 variabelen (ongepaard)

Chikwadraatttoets of Fisher exact toets

De analyse van survivaldata

Onderzoekssituatie Statistische Techniek

Time event studie

Kaplan Meier met Log-ranktoets

Time event met covariaten

Cox-regressieanalyse


Normaliteit en data transformatie

De analyse van de normaliteit van de data is samen te vatten in 6 stapen.

Descriptieve statistiek:

  1. Visuele inspectie met een scatterplot/histogram
  2. Statistische inspectie met de skewness en kurtosis

Inferenti?le statistiek:

  1. Visuele inspectie met een QQ-plot
  2. Statistische analyse met Kolgomorov-Smirnof en Shapiro-Wilks
  3. Overweeg data-transformatie:
    1. Data transformatie
    2. Herhaal statistische analyse
  4. Kies voor parametrisch of non-parametrische analyse

Visuele inspectie

# Visuele inspectie van de leeftijd
par(mfrow=c(2,2))  # maak een 2 bij 2 mutiplot
x <- mydata$age    # defineer de variabele om te onderzoeken

# 1. Scatterplot (linksboven)
plot(x)

# 2. Histogram (rechtsboven)
hist(x)

# 3. qq-plot (linksonder)
qqnorm(x)
qqline(x)

# 4. p-plot (rechtsboven)
library(e1071)
probplot(x, qdist=qnorm)

Lees meer

Statistische analyse

Goede packages uit R voor het testen van normaliteit zijn de packages: fBasics en nortest. Deze bieden een groot scale aan testen voor normaliteit. De twee meest gebruikelijke functies worden besproken.

Keuze voor statistische test Over het algemeen is de Kolmogorov-Smirnov test minder sensitief bij een kleine steekproef (N). Dan kan men beter voor de Shapiro-Wilks test kiezen. Is de steekproef erg groet (1000+) kan weer beter gekozen worden voor de kolmogorov. Lees meer.

Kolmogorov-Smirnov

# Kolmogorov-Smirnov test
library(nortest)
lillie.test(mydata$age)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  mydata$age
## D = 0.0948, p-value = 7.709e-14

Shapiro-Wilks

# Shapiro-Wilks test
shapiro.test(mydata$age)
## 
##  Shapiro-Wilk normality test
## 
## data:  mydata$age
## W = 0.9535, p-value = 1.012e-12

Data transformatie

Let op: Indien gekozen wordt voor een data transformatie is het nodig om uit te leggen waarom ervoor gekozen is en waarom de transformatie geschikt is.

Log transformatie:

# gebruik de functie transform()
# voorbeeld van een log transformatie
log <- transform(mydata, age = log(age))   
# Dit maakt een nieuw data.frame met als verschil dat de leeftijd nu een log transformatie is.

# Verricht opnieuw een visuele inspectie
par(mfrow=c(2,2)); x <- log$age; 
plot(x); hist(x); qqnorm(x); qqline(x); library(e1071); probplot(x, qdist=qnorm)

# Test opnieuw voor normaliteit
# Shapiro-Wilks test
shapiro.test(log$age)
## 
##  Shapiro-Wilk normality test
## 
## data:  log$age
## W = 0.9646, p-value = 9.193e-11
# Kolmogorov-Smirnov test
library(nortest)
lillie.test(log$age)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  log$age
## D = 0.0946, p-value = 9.187e-14

Rank transformatie:

# gebruik de functie transform()
# voorbeeld van een log transformatie
rank <- transform(mydata, age = rank(age))   
# Dit maakt een nieuw data.frame met als verschil dat de leeftijd nu een log transformatie is.

# Verricht opnieuw een visuele inspectie
par(mfrow=c(2,2)); x <- rank$age; 
plot(x); hist(x); qqnorm(x); qqline(x); library(e1071); probplot(x, qdist=qnorm)

# Test opnieuw voor normaliteit
# Shapiro-Wilks test
shapiro.test(rank$age)
## 
##  Shapiro-Wilk normality test
## 
## data:  rank$age
## W = 0.9506, p-value = 3.492e-13
# Kolmogorov-Smirnov test
library(nortest)
lillie.test(rank$age)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  rank$age
## D = 0.1002, p-value = 1.544e-15

Interpretatie:

De data is niet normaal verdeeld en hoewel het na een log transformatie meer op een binomiale verdeling lijkt blijven de statistische testen zeer significant en moet de null hypothese van normaal verdeelde data worden verworpen. Er kan nu beter gekozen worden voor een non parametrische test. Het is ook mogelijk om te testen op extreme outliers en overwegen deze te verwijderen. Vervolgens kan opnieuw getest worden of de data normaal verdeeld is.


Outliers

Outliers kunnen technisch incorrect zijn (denk aan een typefout) of statistisch onwaarschijnlijk. Indien je statistisch onwaarschijnlijke outliers wilt detecteren en recoderen of verwijderen alvorens een statistische analyse uit te voeren kan dat middels verschillende methodes. Er worden 2 methodes besproken.

Standaard deviatie

Een methode voor het detecteren en verwijderen van NA’s is met de standaard deviatie. Met de aanname dat de steekproef uit een normaal verdeelde populatie komt verwacht je dat 99,7% van de steekproef binnen +/- 3 standaard deviatie’s (SD) valt.

# Verander alle waarden uit een variabele als deze >3SD afwijkt van het gemiddelde
data <- mydata$age  # defineer variabele om te recoderen
data <- data.frame(data)
data <- data.frame(data, abs((data[,1] - mean(data[,1])))/sd(data[,1]))       # voeg tau waarde toe
data <- data.frame(data, abs((data[,1] - mean(data[,1])))/sd(data[,1]) > 3.0) # voeg logical TRUE/FALSE toe als de tau waarde > 3 is
colnames(data) <- c("x", "tau", "outlier") 

# Wat moet er gebeuren met de outliers
data$x[data$outlier == TRUE] <- NA     # verander outliers in NA
                 
# Defineer nieuwe variabele met de outliers +/- 3 SD verwijderd.
mydata$ageNAoutlier <- data$x

Lees meer Over outlier detectie.

Smirnov-Grubbs Test

Een andere methode van outlier verwijdering is met de Smirnov-Grubb test. Dit is een null hypothese significantie test met als H0: er zijn geen outliers. De alternatieve hypothese H1: er is minstens 1 outlier. Deze test vertelt niet hoeveel outliers er zijn.

library(outliers)
# gebruik de functie grubbs.test() om de 'Smirnov-Grubb test' uit te voeren.
grubbs.test(mydata$age)
## 
##  Grubbs test for one outlier
## 
## data:  mydata$age
## G = 3.2943, U = 0.9816, p-value = 0.2784
## alternative hypothesis: highest value 98 is an outlier

Lees meer.

Visualiseren

Het is mogelijk om de grubbs-test te gebruiken in een loop-functie om zo alle statistische outliers te detecteren, visualiseren en verwijderen. Onderstaande functie is gebaseerd op deze discussie.

# 1 bij 2 grafiek tabel
par(mfrow=c(2,1))

# Maak een boxplot met outliers
boxplot(mydata$mmse, horizontal=TRUE,
        main="Boxplot van de MMSE scores Outliers",             # De titel van de grafiek
        sub="punten vallen buiten 1,5 IQR"             # De subtitel van de grafiek
       )
## Maak een boxplot met extremes 
# bron: http://stackoverflow.com/questions/4994313/changing-the-outlier-rule-in-a-boxplot
boxplot(mydata$mmse, horizontal=TRUE,
        range = 3,                                     # Gebruik range = 3 voor het veranderen van de boxplot parameters
        main="Boxplot van de MMSE scores Extremes",             # De titel van de grafiek
        sub="punten vallen buiten 3 IQR"               # De subtitel van de grafiek
       )

## Dit is de functie
# bron: http://stackoverflow.com/questions/22837099/how-to-repeat-the-grubbs-test-and-flag-the-outliers
library(outliers)
library(ggplot2)

outliervisual <- function (X) {
  grubbs.flag <- function(x) {
    outliers <- NULL
    test <- x
    grubbs.result <- grubbs.test(test)
    pv <- grubbs.result$p.value
    while(pv < 0.05) {
      outliers <- c(outliers,as.numeric(strsplit(grubbs.result$alternative," ")[[1]][3]))
      test <- x[!x %in% outliers]
      grubbs.result <- grubbs.test(test)
      pv <- grubbs.result$p.value
    }
    return(data.frame(X=x,Outlier=(x %in% outliers)))
  }
  grubbs.flag(X)
  ggplot(grubbs.flag(X),aes(x=X,color=Outlier,fill=Outlier))+
    geom_histogram(binwidth=diff(range(X))/30)+
    theme_bw()
}

# Pas de nieuw gemaakte functie toe
outliervisual(mydata$mmse) # het visualiseert de extremen (>+/-3 IQR)


Missing values

Visualiseren

Een goede manier voor het visualiseren van de missing values is met de ’VIM’package.

library(VIM)
aggr(mydata, plot = TRUE)  

Imputeren

Een statistische methode voor het omgaan met missing values is imputeren. Met imputeren vervang je de NA door een “voorspelling”. Dit kan met verschillende methodes. Als je imputeert met een schatting zoals via regressie, een gemiddelde of een groepsgemiddelde dan zal de verdeling van de desbetreffende variabele gepiekter worden. Hoe meer je imputeert, hoe meer waarden in het centrum van de verdeling komen te liggen. De variantie van de variabele zal daardoor afnemen. Dit kan gevolgen hebben voor de statistische analyse. Bron.

## Voorbeeld van een imputatie met Hmisc package
library(Hmisc)
DF <- data.frame(age = c(10, 20, NA, 40), sex = c('male','female'))

# imputeer met afgeronde gemiddelde
with(DF, round(impute(age, mean), digits=2))
##      1      2      3      4 
##  10.00  20.00 23.33*  40.00
# impute met een random waarde
with(DF, impute(age, 'random'))
##   1   2   3   4 
##  10  20 40*  40
# imputeer met de mediaan
with(DF, impute(age, median))
##   1   2   3   4 
##  10  20 20*  40
# imputeer met het minimum
with(DF, impute(age, min))
##   1   2   3   4 
##  10  20 10*  40
# imputeer met het maximum
with(DF, impute(age, max))
##   1   2   3   4 
##  10  20 40*  40

Lees meer over omgaan met missing values.

Er zijn voor R vele packages beschikbaar voor imputeren, kies een package op basis van de methodes die het meest geschikt is voor jou data:

Imputatietabel

Lees meer over multipele imputaties of bekijk de presentatie.


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