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.


Karakteristieken

De eerste stap in het analyseren van de data is het beschrijven van de variabelen in statistische termen. Descriptieve statistiek heeft als doel het samenvatten van de steekproef. Het gebruiken van de steekproef om uitspraken te doen over de gehele populatie valt onder inferentiële statistiek. Karakteristieken valt onder de univariate data analyse.

Voor het beschrijven van van de karakteristieken zijn er twee vormen: het beschrijven van de Central Tendancy en de Variance. Lees meer.

Centrum maten

De maten voor het beschrijven van het centrum van de variabele zijn:

  • Continue data:
    • gemiddelde: mean()
    • mediaan: median()
  • Categoriale data:
    • mode: de meest voorkomende waarde
## continue data
# gemiddelde leeftijd afgerond op gehele getallen
round(mean(mydata$age),0)
## [1] 78
# gemiddelde visusscore
mean(mydata$visus)
## [1] 0.4174283
# mediane leeftijd: helft van de patiënten is jonger, de helft is ouder
median(mydata$age)
## [1] 78
# mediane MMSE score
median(mydata$mmse)
## [1] 26
## categoriale data
# mode: de meest voorkomende geslacht
mode <- table(as.vector(mydata$sexe))
names(mode)[mode == max(mode)]
## [1] "vrouw"

Lees meer over Mode. Het beschrijven van categoriale data kan beter met een grafische weergave zoals bijvoorbeeld een histogram of een stem & leaf.

Spreidingsmaten

De spreiding van continue data wordt beschreven met de volgende maten:

  • variantie: var()
  • minimum en maximum: min() en max()
  • kwantielen: quantiles()
    • interquantile range: IQR()
  • standaard deviatie: sd()
  • standaardfout: sd(x)/sqrt(length(x))
# variatie in leeftijd
var(mydata$age, na.rm = FALSE, use = "everything")
## [1] 35.35098
# De leeftijd is dus gemiddeld 78 met een variatie (min tot max) van ~35 jaar.

# De minimum en maximum van de apache score
min(mydata$apache)
## [1] 7
max(mydata$apache)
## [1] 31
# kwantielen van visus
quantile(mydata$visus)
##    0%   25%   50%   75%  100% 
## 0.025 0.300 0.400 0.500 0.700
# kwantielen zijn net als de mediaan een afkappunt. De visus loopt dus van vrijwel niks naar redelijke zicht (70%), er zijn dus geen ouderen in deze patiëntpopulatie met perfecte visus (visus=1). Quantielen zijn erg geschikt om OUTLIERS te vinden.

# Interquantile Range
# De observaties tussen de tweede en vierde kwantiel. Met andere woorden: 50% van de observaties vallen tussen de IQR. zie ?IQR() voor meer opties.
IQR(mydata$visus)
## [1] 0.2
# Standaard deviatie
mean(mydata$mmse); sd(mydata$mmse)
## [1] 25.27993
## [1] 4.160937
# De gemiddelde MMSE score is ~25 met een SD van 4.

# Standaard error of mean is de standaard deviatie gedeeld door wortel van N.
x <- mydata$mmse
std <- sd(x)/sqrt(length(x)) # length(x) is de N (de hoeveelheid patiënten)
print(std)
## [1] 0.1708692

Wat heb ik aan de standaarddeviatie: De combinatie van een gemiddelde en de standaard deviatie vertelt wat over de distributie van de data. De empirische wet van Bell stelt dat de ca. 68% van de data binnen +-1 SD valt. Ca. 95% van de data valt binnen +-2SD en 99.7% van de data valt binnen +-3SD indien er sprake is van een normaalverdeling. Lees meer.

Standaard Deviatie

Wat heb ik aan een standaardfout: De standaardfout is in principe kleiner naarmate de steekproef groter is. De standaard fout in het steekproefgemiddelde is recht evenredig met de standaarddeviatie van de populatie waaruit de steekproef is getrokken en omgekeerd evenredig met de wortel van het aantal onafhankelijke waarnemingen in de steekproef. Lees meer.

Verschil standaarddeviatie en standaardfout: De standaarddeviatie is een beschrijvende maat waarmee je de positie ten opzichte van het gemiddelde aangeeft en de standaardfout is een maat voor de betrouwbaarheid waarmee je een significant verschil met het gemiddelde kunt berekenen. Lees meer.

Waarom kwantielen berekenen als ik al een standaarddeviatie heb: Bij een kleine patiëntpopulatie is de standaarddeviatie onbetrouwbaar. Dan kan beter worden gekozen van de mediaan en de kwantielen. Deze worden gebruikt voor het maken van een boxplot en zijn daardoor geschikt voor het visualiseren van outliers. Lees meer.

Distributie

De distributie van de data wordt beschreven met twee maten:

  • Kurtosis: gewelfdheid
  • Skewness: scheefheid (“tendancy to assymetry”)

Het beschrijven van de kurtosis en skewness kan het beste met behulp van een Package. Gebruik bijvoorbeeld de package e1071.

# Laad de package e1071
library(e1071)

# Bereken de Kurtosis van de leeftijd
kurtosis(mydata$age)
## [1] -0.07553439
# Bereken de Skewness van de leeftijd
skewness(mydata$age)
## [1] 0.6701232

Waarom Kurtosis en Skewness bepalen: Omdat deze twee iets zeggen over de normaalverdeling. Het is voor het kiezen van de juiste statische test noodzakelijk om te bepalen of de data normaal verdeeld is. Is de data niet normaal verdeeld dan moet men niet-parametrische testen gebruiken. Het testen van de normaalverdeling gebeurt met de Kolmogorov-Smirnoff toets en de Shapiro-Wilk toets.

Interpretatie:

  • Kurtosis
    • Kurtosis = 3: Mesokurtic distribution - normaalverdeling.
    • Kurtosis > 3: Leptokurtic distribution - scherper dan een normaalverdeling met waarden geconcentreerd rond het gemiddelde met een dikke staart. Er is een hogere kans op extreme outliers dan bij een normaalverdeling.
    • Kurtosis < 3: Platykurtic distribution - platter dan een normaalverdeling met een wijdere piek. De kans op extreme outliers is kleiner dan bij een normaalverdeling en de data is ver verspreid van het gemiddelde.

kurtosis

  • Skewness
    • positieve skew: de data is assymetrisch naar links afgebogen
    • negatieve skew: de data is assymetrisch naar rechts afgebogen

Skewness Skewness

Lees meer over distributie, Skweness en Kurtosis.

Samenvatting

R heeft ingebouwde functies om alle karakteristieken die hierboven zijn besproken snel weer te geven per variabele (kolom) in een data.frame: summary(). Daarnaast zijn er Packages beschikbaar die extra functionaliteit bieden, zoals bijvoorbeeld de Package Hmisc met de functie describe() of de package psych met de functie describe().

# Voorbeeld van summary
summary(mydata)
##      gebdat              sexe          medicijn       apache     
##  Min.   :1903-01-11   man  :137   geen med.:172   Min.   : 7.00  
##  1st Qu.:1919-11-04   vrouw:456   haldol   :208   1st Qu.:11.00  
##  Median :1923-09-17               placebo  :213   Median :12.00  
##  Mean   :1923-02-24                               Mean   :13.04  
##  3rd Qu.:1927-12-16                               3rd Qu.:15.00  
##  Max.   :1935-09-11                               Max.   :31.00  
##       mmse           visus           creaurer             delier   
##  Min.   : 0.00   Min.   :0.0250   Min.   : 0.00   delier     : 71  
##  1st Qu.:24.00   1st Qu.:0.3000   1st Qu.:10.16   geen delier:522  
##  Median :26.00   Median :0.4000   Median :12.20                    
##  Mean   :25.28   Mean   :0.4174   Mean   :12.69                    
##  3rd Qu.:28.00   3rd Qu.:0.5000   3rd Qu.:14.31                    
##  Max.   :30.00   Max.   :0.7000   Max.   :39.51                    
##      opndat               ontdat               electief  
##  Min.   :2000-08-07   Min.   :2000-08-25   acuut   :133  
##  1st Qu.:2001-02-07   1st Qu.:2001-03-08   electief:460  
##  Median :2001-07-25   Median :2001-08-17                 
##  Mean   :2001-07-24   Mean   :2001-08-11                 
##  3rd Qu.:2001-12-24   3rd Qu.:2002-01-13                 
##  Max.   :2002-11-27   Max.   :2002-08-15                 
##     datdood                       dood          age       
##  Min.   :2000-09-09   niet overleden:401   Min.   :66.00  
##  1st Qu.:2002-05-27   overleden     :192   1st Qu.:74.00  
##  Median :2004-01-01                        Median :78.00  
##  Mean   :2004-01-04                        Mean   :78.41  
##  3rd Qu.:2005-08-15                        3rd Qu.:82.00  
##  Max.   :2007-05-03                        Max.   :98.00
## Probeer deze functie's zelf uit. De output wordt hier NIET gegenereerd. In plaats van mydata kan ook de dataset 'mtcars' worden gebruikt voor een snelle overzicht, deze zit standaard in R.
# Voorbeeld van describe met Hmisc
library(Hmisc)
describe(mydata)

# Voorbeeld van describe met psych
library(psych)
describe(mydata)

Samenvatting per groep

Naast een samenvatting van alle data zoals hierboven is het ook mogelijk om een samenvatting te geven van de data per groep. Dit kan met de functie summaryBy() uit de Package doBy. De Package psych heeft hier ook een functie voor: describeBy.

structuur summaryBy(): summaryBy(formule, data=dataframe, FUN=function)
De formule heeft de vorm: var1 + var2 + var3 + … + varN ~ groupvar1 + groupvar2 + … + groupvarN

## Probeer deze functie's zelf uit. De output wordt hier NIET gegenereerd.

# Voorbeeld van summary.by() met doBY. Data: leeftijd, groep: medicijn, functie: summary()
library(doBy)
summaryBy(age~medicijn, data=mydata, FUN=summary)
##    medicijn age.Min. age.1st Qu. age.Median age.Mean age.3rd Qu. age.Max.
## 1 geen med.       66          72         74    75.48          78       89
## 2    haldol       70          74         79    79.06          83       94
## 3   placebo       70          75         79    80.15          84       98
# Voorbeeld van describe() met psych. Data: leeftijd, groep: medicijn.
library(psych)
describeBy(mydata$age, mydata$medicijn)
## group: geen med.
##   vars   n  mean   sd median trimmed  mad min max range skew kurtosis   se
## 1    1 172 75.48 4.24     74   75.11 4.45  66  89    23 0.75     0.15 0.32
## -------------------------------------------------------- 
## group: haldol
##   vars   n  mean   sd median trimmed  mad min max range skew kurtosis   se
## 1    1 208 79.06 5.85     79   78.82 7.41  70  94    24 0.29    -0.74 0.41
## -------------------------------------------------------- 
## group: placebo
##   vars   n  mean   sd median trimmed  mad min max range skew kurtosis   se
## 1    1 213 80.15 6.35     79    79.7 5.93  70  98    28 0.61    -0.27 0.44

Lees meer over Summary Statistics.


Tabellen

R heeft verschillende ingebouwde functies voor het creëren van tabellen. Hieronder volgt een kort overzicht van de functies voor het maken van tabellen:

Functie Beschrijving
table(var1, var2, ., varN) N-way kruistabel voor N categoriale variabele
xtabs(formule, data) N-way kruistabel op basis van een formule (bv gemiddelde)
prop.table(table, margins) Maakt van een absolute frequentietabel een relatieve kruistabel
margin.table(table, margins) Maakt een som per rij (1) en/of kolommen (2)
addmargins(table, margins) Voegt een samenvattingsmarge toe aan een tabel
ftable(table) Maakt een compacte kruistabel van een tabel

Indien het nodig is om data te transformeren naar een tabel voor verdere data manipulatie is het aangeraden gebruik te maken van de functie count() in plaats van table() omdat zo data class data.frame te behouden. Waarom count?

Keuze tabel

In de statistiek worden tabellen gebruikt voor het weergeven van de frequentie distributie. Dit kan op verschillende manieren:

  • Numerieke data:
    • Absolute frequentie distributie
    • Relatieve frequentie distributie
    • Cummulatieve frequentie distributie
  • categoriale data:
    • Absolute frequentie distributie
    • Relatieve frequentie distributie
# Voorbeeld van een absolute frequentietabel voor geslacht
table(mydata$sexe)
## 
##   man vrouw 
##   137   456
# Voorbeeld van een relatieve frequentietabel voor geslacht
prop.table(table(mydata$sexe))
## 
##       man     vrouw 
## 0.2310287 0.7689713

Data transformatie

Een frequentietabel van 1 variabele is alleen geschikt voor categoriale of interval data. Een numerieke variabele (bijvoorbeeld leeftijd) is te recoderen naar een interval met de functie cut().

## Voorbeeld van een recodering ##

head(mydata$age)
## [1] 87 92 83 82 95 94
# Recodeer leeftijd naar 5 groepen (5 groepen van 20% van de data)
newdata <- cut(mydata$age, 5)
# Laat de nieuwe variabele zien
head(newdata)
## [1] (85.2,91.6] (91.6,98]   (78.8,85.2] (78.8,85.2] (91.6,98]   (91.6,98]  
## Levels: (66,72.4] (72.4,78.8] (78.8,85.2] (85.2,91.6] (91.6,98]
# Maak er een tabel van
table(newdata)
## newdata
##   (66,72.4] (72.4,78.8] (78.8,85.2] (85.2,91.6]   (91.6,98] 
##         103         229         176          66          19
# Maak er een relatieve frequentietabel van
prop.table(table(newdata))
## newdata
##   (66,72.4] (72.4,78.8] (78.8,85.2] (85.2,91.6]   (91.6,98] 
##  0.17369309  0.38617201  0.29679595  0.11129848  0.03204047

2x2 Kruistabellen

2x2 kruistabellen worden gebruikt voor het uitdrukken van de variatie van een variabele (A) ten opzichte van een variabele (B). Dit kunnen absolute en relatieve frequentietabellen zijn.

# Voorbeeld van een absolute frequentietabel(A,B) voor (A) geslacht (rijen)   en (B) medicijn (kolomen)
mytable <- table(mydata$sexe,mydata$medicijn)
print(mytable)
##        
##         geen med. haldol placebo
##   man          52     40      45
##   vrouw       120    168     168
# Voorbeeld van een (getransformeerde) relatieve frequentietabel
prop.table(mytable) # fractie per cel (alle cellen bij elkaar opgeteld = 1)
##        
##          geen med.     haldol    placebo
##   man   0.08768971 0.06745363 0.07588533
##   vrouw 0.20236088 0.28330523 0.28330523
prop.table(mytable, 1) # fractie per rij (som per rij = 1)
##        
##         geen med.    haldol   placebo
##   man   0.3795620 0.2919708 0.3284672
##   vrouw 0.2631579 0.3684211 0.3684211
prop.table(mytable, 2) # fractie per kollom (som per kolom = 1)
##        
##         geen med.    haldol   placebo
##   man   0.3023256 0.1923077 0.2112676
##   vrouw 0.6976744 0.8076923 0.7887324

Esthetisch aantrekkelijke en informatieve (zowel absoluut als relatieve) 2x2 kruistabellen kunnen worden gemaakt met de functie CrossTable() uit de Package gmodels. Lees meer over deze functie. Bijvoorbeeld:

# Kruistabel van geslacht en medicijn.
library(gmodels)

# Maak een kruistabel, gebruik expected=TRUE om een Pearson's Chi-squared test toe te voegen.
CrossTable(mydata$sexe, mydata$medicijn, expected = TRUE) 
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |              Expected N |
## | Chi-square contribution |
## |           N / Row Total |
## |           N / Col Total |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  593 
## 
##  
##              | mydata$medicijn 
##  mydata$sexe | geen med. |    haldol |   placebo | Row Total | 
## -------------|-----------|-----------|-----------|-----------|
##          man |        52 |        40 |        45 |       137 | 
##              |    39.737 |    48.054 |    49.209 |           | 
##              |     3.784 |     1.350 |     0.360 |           | 
##              |     0.380 |     0.292 |     0.328 |     0.231 | 
##              |     0.302 |     0.192 |     0.211 |           | 
##              |     0.088 |     0.067 |     0.076 |           | 
## -------------|-----------|-----------|-----------|-----------|
##        vrouw |       120 |       168 |       168 |       456 | 
##              |   132.263 |   159.946 |   163.791 |           | 
##              |     1.137 |     0.406 |     0.108 |           | 
##              |     0.263 |     0.368 |     0.368 |     0.769 | 
##              |     0.698 |     0.808 |     0.789 |           | 
##              |     0.202 |     0.283 |     0.283 |           | 
## -------------|-----------|-----------|-----------|-----------|
## Column Total |       172 |       208 |       213 |       593 | 
##              |     0.290 |     0.351 |     0.359 |           | 
## -------------|-----------|-----------|-----------|-----------|
## 
##  
## Statistics for All Table Factors
## 
## 
## Pearson's Chi-squared test 
## ------------------------------------------------------------
## Chi^2 =  7.145067     d.f. =  2     p =  0.02808462 
## 
## 
## 

Ook kan gebruik worden gemaakt van de functie crosstab(). Deze functie heeft geen package maar is wel hier te vinden inclusief een uitleg.

2+ groepen

Voor het maken van een frequentietabel van meer dan 2 variabelen kan dat met de functie ftable() voor het transformeren van een tabel. Ook kan gebruik worden gemaakt van de functie xtabs() voor het maken van een tabel.

# tussenstap: Maak een tabel
mytable <- table(mydata$sexe, mydata$medicijn, mydata$delier)
print(mytable)
## , ,  = delier
## 
##        
##         geen med. haldol placebo
##   man           0      8      12
##   vrouw         6     23      22
## 
## , ,  = geen delier
## 
##        
##         geen med. haldol placebo
##   man          52     32      33
##   vrouw       114    145     146
# Transformeer naar frequentietabel met ftable()
ftable(mytable) 
##                  delier geen delier
##                                    
## man   geen med.       0          52
##       haldol          8          32
##       placebo        12          33
## vrouw geen med.       6         114
##       haldol         23         145
##       placebo        22         146
# frequentietabel met xtabs()
mytable <- xtabs(~sexe+medicijn+delier, data=mydata)
# Transformeer met ftable()
ftable(mytable) # maak er nu een 3-staps tabel van
##                 delier delier geen delier
## sexe  medicijn                           
## man   geen med.             0          52
##       haldol                8          32
##       placebo              12          33
## vrouw geen med.             6         114
##       haldol               23         145
##       placebo              22         146

Analyse

Met het visueel inspecteren van kruistabellen kan worden geïnspecteerd of er een verband is tussen de variabelen. De statistische analyse van een kruistabel kan met verschillende testen:

  • Associatie
    • Cramers V: associatiemaat tussen twee categorische variabelen
    • Yules Q: een associatiemaat tussen twee dichotome variabelen
    • Yules Y: een associatiemaat tussen twee dichotome variabelen
    • Kendalls Tau*: associatiemaat gebaseerd op de rangnummers (non-parametrisch) van de data
  • Verschil
    • Chi-kwadraattoets*
    • Fishers exacte toets*

* Deze analyses komen aan bod bij de Statistische Analyse


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