Ga naar hoofdinhoud

R-tutorial kruistabellen

In deze tutorial leer je hoe je kruistabellen maakt en hoe je relaties die daarin zichtbaar zijn test en kwantificeert.
· 10 min lezen

Verken met AI

Hoe maak je een tabel

Laten we eerst wat data pakken. Het pakket MASS bevat gegevens over 93 auto's die in 1993 in de VS te koop waren. Ze zijn opgeslagen in het object Cars93 en bevatten 27 kenmerken per auto, waarvan sommige categorisch zijn. Laten we dus het pakket MASS laden en kijken naar het type voertuigen in cars93:

library(MASS)
Cars93$Type

##  [1] Small   Midsize Compact Midsize Midsize Midsize Large   Large  
##  [9] Midsize Large   Midsize Compact Compact Sporty  Midsize Van    
## [17] Van     Large   Sporty  Large   Compact Large   Small   Small  
## [25] Compact Van     Midsize Sporty  Small   Large   Small   Small  
## [33] Compact Sporty  Sporty  Van     Midsize Large   Small   Sporty
## [41] Sporty  Small   Compact Small   Small   Sporty  Midsize Midsize
## [49] Midsize Midsize Midsize Large   Small   Small   Compact Van    
## [57] Sporty  Compact Midsize Sporty  Midsize Small   Midsize Small  
## [65] Compact Van     Midsize Compact Midsize Van     Large   Sporty
## [73] Small   Compact Sporty  Midsize Large   Compact Small   Small  
## [81] Small   Compact Small   Small   Sporty  Midsize Van     Small  
## [89] Van     Compact Sporty  Compact Midsize
## Levels: Compact Large Midsize Small Sporty Van

We hebben daar 6 autotypen. De functie table vertelt hoeveel we er van elk type hebben:

table(Cars93$Type)

##
## Compact   Large Midsize   Small  Sporty     Van
##      16      11      22      21      14       9

prop.table zet dit om naar fracties:

prop.table(table(Cars93$Type))

##
##    Compact      Large    Midsize      Small     Sporty        Van
## 0.17204301 0.11827957 0.23655914 0.22580645 0.15053763 0.09677419

Hetzelfde met de herkomst van de auto's:

table(Cars93$Origin)

##
##     USA non-USA
##      48      45

prop.table(table(Cars93$Origin))

##
##      USA  non-USA
## 0.516129 0.483871

Hoe maak je een kruistabel

Mooi, we zagen dat onze dataset een vergelijkbaar aantal Amerikaanse en niet-Amerikaanse auto's bevat en dat de meest voorkomende typen Midsize en Small zijn. Maar misschien verschillen de Amerikaanse en niet-Amerikaanse auto's in type?

Laten we naar de autotypen kijken in relatie tot hun herkomst. We kunnen table opnieuw gebruiken, maar nu met twee argumenten. Het eerste wordt de rijvariabele en het tweede de kolomvariabele:

table(Cars93$Type, Cars93$Origin)

##          
##           USA non-USA
##   Compact   7       9
##   Large    11       0
##   Midsize  10      12
##   Small     7      14
##   Sporty    8       6
##   Van       5       4

Nu zagen we wat iedereen weet: Amerikanen houden van grote voertuigen! De bovenstaande tabel laat de gezamenlijke verdeling zien van twee categorische variabelen (Type en Origin). Zulke tabellen noemen we kruistabellen.

Hoe krijg je marginalen uit een kruistabel

De functies rowSums en colSums spreken eigenlijk voor zich

(tab1<-table(Cars93$Type, Cars93$Origin))

##          
##           USA non-USA
##   Compact   7       9
##   Large    11       0
##   Midsize  10      12
##   Small     7      14
##   Sporty    8       6
##   Van       5       4

rowSums(tab1)

## Compact   Large Midsize   Small  Sporty     Van
##      16      11      22      21      14       9

colSums(tab1)

##     USA non-USA
##      48      45

Hoe krijg je percentages uit een kruistabel

prop.table genest met table geeft frequenties

prop.table(table(Cars93$Type, Cars93$Origin))

##          
##                  USA    non-USA
##   Compact 0.07526882 0.09677419
##   Large   0.11827957 0.00000000
##   Midsize 0.10752688 0.12903226
##   Small   0.07526882 0.15053763
##   Sporty  0.08602151 0.06451613
##   Van     0.05376344 0.04301075

Omzetten naar percentages is simpelweg vermenigvuldigen met 100:

prop.table(table(Cars93$Type, Cars93$Origin))*100

##          
##                 USA   non-USA
##   Compact  7.526882  9.677419
##   Large   11.827957  0.000000
##   Midsize 10.752688 12.903226
##   Small    7.526882 15.053763
##   Sporty   8.602151  6.451613
##   Van      5.376344  4.301075

Let op: dit is een gezamenlijke kansverdeling, waaruit we bijvoorbeeld zien dat zo'n 7,5% van de auto's klein is en van Amerikaanse herkomst.

Vaker zijn we geïnteresseerd in de verdeling van één variabele binnen groepen die door een andere variabele zijn gevormd. Hier lijkt de verdeling van autotypen onder de Amerikaanse en (afzonderlijk) niet-Amerikaanse auto's interessant. Hiervoor gebruiken we het argument margin van de functie prop.table. Het geeft aan waar de groeperingsvariabele staat: in rijen (margin=1) of in kolommen (margin=2):

prop.table(table(Cars93$Type, Cars93$Origin), margin=2)*100

##          
##                 USA   non-USA
##   Compact 14.583333 20.000000
##   Large   22.916667  0.000000
##   Midsize 20.833333 26.666667
##   Small   14.583333 31.111111
##   Sporty  16.666667 13.333333
##   Van     10.416667  8.888889

Nu zien we makkelijk dat kleine auto's twee keer zo vaak voorkomen in het niet-Amerikaanse deel van onze dataset als in het Amerikaanse.

Let ook op dat kolommen tot 100 optellen, terwijl in de gezamenlijke verdelingstabel (degene zonder het argument margin) 100 de som van de hele tabel was.

(tab2<-prop.table(table(Cars93$Type, Cars93$Origin), margin=2)*100)

##          
##                 USA   non-USA
##   Compact 14.583333 20.000000
##   Large   22.916667  0.000000
##   Midsize 20.833333 26.666667
##   Small   14.583333 31.111111
##   Sporty  16.666667 13.333333
##   Van     10.416667  8.888889

colSums(tab2)

##     USA non-USA
##     100     100

Chi-kwadraattoets

De meest voorkomende vraag bij kruistabellen is of de rij- en kolomvariabelen onafhankelijk zijn. De meest basale manier om dit te beantwoorden is een chi-kwadraattoets uitvoeren. Die wordt uitgebreid behandeld in deze tutorial. Laten we controleren of Type en Origin onafhankelijk zijn:

chisq.test(Cars93$Type, Cars93$Origin)

## Warning in chisq.test(Cars93$Type, Cars93$Origin): Chi-squared
## approximation may be incorrect

##
##  Pearson's Chi-squared test
##
## data:  Cars93$Type and Cars93$Origin
## X-squared = 14.08, df = 5, p-value = 0.01511

Blijkbaar niet, maar we kregen ook de waarschuwing Chi-squared approximation may be incorrect. Dit komt doordat de chi-kwadraatstatistiek slechts bij benadering een chi-kwadraatverdeling volgt. Hoe meer observaties we hebben, hoe beter de benadering. De functie chisq.test geeft bovenstaande waarschuwing zodra een van de verwachte aantallen lager is dan 5 (voor wat een "expected count" is, zie de gelinkte tutorial).

Exacte toets van Fisher

De exacte toets van Fisher is een alternatief voor de chi-kwadraattoets, dat vooral wordt gebruikt wanneer de chi-kwadraatbenadering niet toereikend is. Laten we die uitvoeren:

fisher.test(Cars93$Type, Cars93$Origin)

##
##  Fisher's Exact Test for Count Data
##
## data:  Cars93$Type and Cars93$Origin
## p-value = 0.007248
## alternative hypothesis: two.sided

De resultaten lijken behoorlijk op die van de chi-kwadraat, maar dat is niet altijd zo. Een belangrijk nadeel van de exacte toets van Fisher is dat deze voor grote tabellen (of grote steekproeven) computationeel inefficiënt wordt.

G-toets

Een ander alternatief is de zogeheten G-toets. De statistiek volgt ook bij benadering een chi-kwadraatverdeling, maar voor kleine steekproeven ligt deze benadering dichterbij dan die van de chi-kwadraattoets. Voor de G-toets kunnen we de functie GTest uit het pakket DescTools gebruiken. De resultaten zijn opnieuw vergelijkbaar met de twee vorige toetsen: Type en Origin zijn niet onafhankelijk.

library(DescTools)
GTest(Cars93$Type, Cars93$Origin)

##
##  Log likelihood ratio (G-test) test of independence without
##  correction
##
## data:  Cars93$Type and Cars93$Origin
## G = 18.362, X-squared df = 5, p-value = 0.002526

Correctie van Yates

In 2x2-kruistabellen kan de chi-kwadraattoets worden verbeterd met de continuïteitscorrectie van Yates. Deze trekt simpelweg 0,5 af van elke term | Observed - Expected | in de chi-kwadraatstatistiek. Raadpleeg opnieuw deze tutorial als je nu de draad kwijt bent. Bovendien past R de correctie van Yates automatisch toe wanneer dat nodig is. Laten we kijken naar de beschikbaarheid van uitvoeringen met handgeschakelde transmissie bij Amerikaanse en niet-Amerikaanse auto's:

(tab3<-table(Cars93$Man.trans.avail, Cars93$Origin))

##      
##       USA non-USA
##   No   26       6
##   Yes  22      39

chisq.test(tab3)

##
##  Pearson's Chi-squared test with Yates' continuity correction
##
## data:  tab3
## X-squared = 15.397, df = 1, p-value = 8.712e-05

De correctie van Yates werd automatisch toegepast, en R gaf dat aan.

3 (of meer) dimensionale tabel

Laten we tegelijk kijken naar Man.trans.avail, Origin en Type. table splitst de dataset op volgens de als derde opgegeven variabele:

table(Cars93$Man.trans.avail, Cars93$Origin, Cars93$Type)

## , ,  = Compact
##
##      
##       USA non-USA
##   No    2       0
##   Yes   5       9
##
## , ,  = Large
##
##      
##       USA non-USA
##   No   11       0
##   Yes   0       0
##
## , ,  = Midsize
##
##      
##       USA non-USA
##   No    9       4
##   Yes   1       8
##
## , ,  = Small
##
##      
##       USA non-USA
##   No    0       0
##   Yes   7      14
##
## , ,  = Sporty
##
##      
##       USA non-USA
##   No    0       0
##   Yes   8       6
##
## , ,  = Van
##
##      
##       USA non-USA
##   No    4       2
##   Yes   1       2

ftable geeft een compacter overzicht:

ftable(Cars93$Man.trans.avail, Cars93$Origin, Cars93$Type)

##              Compact Large Midsize Small Sporty Van
##                                                    
## No  USA            2    11       9     0      0   4
##     non-USA        0     0       4     0      0   2
## Yes USA            5     0       1     7      8   1
##     non-USA        9     0       8    14      6   2

Cochran-Mantel-Haenszel-toets

Uit de bovenstaande table-resultaten zou duidelijk moeten zijn dat de relatie tussen Origin en Man.trans.avail verschilt per Type. Voor kleine en sportieve auto's is er bijvoorbeeld geen verband: elke Amerikaanse, evenals elke niet-Amerikaanse kleine of sportieve auto, is verkrijgbaar in een handgeschakelde versie. Aan de andere kant hebben de meeste middelgrote Amerikaanse auto's geen handgeschakelde versie, terwijl de meeste niet-Amerikaanse auto's van dit type die wel hebben. Om rekening te houden met mogelijk verschillende relaties per laag (stratum), kan de Cochran-Mantel-Haenszel-toets worden gebruikt:

mantelhaen.test(Cars93$Man.trans.avail, Cars93$Origin, Cars93$Type)

##
##  Mantel-Haenszel chi-squared test with continuity correction
##
## data:  Cars93$Man.trans.avail and Cars93$Origin and Cars93$Type
## Mantel-Haenszel X-squared = 8.0153, df = 1, p-value = 0.004638
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##   2.226531 76.891307
## sample estimates:
## common odds ratio
##          13.08438

Het derde argument van mantelhaen.test identificeert de strata. Vergelijk bovenstaande resultaten met die zonder stratificatie (het voorbeeld met de correctie van Yates hierboven). Er is nog steeds een verband, maar het bewijs ervoor is een stuk zwakker.

Associatiematen

Nu we enkele verbanden tussen variabelen hebben ontdekt, is het tijd om hun sterkte te meten. Er zijn tal van mogelijke maten. Veel daarvan worden hier beschreven. Ik focus nu op de twee meest gebruikte.

Cramers V

V is gebaseerd op de chi-kwadraatstatistiek:

$$ V = \sqrt{\frac{\chi^2/N}{\min(C-1, R-1)}}, $$ waarbij:

  • N de totale som in de kruistabel is (som van alle cellen),
  • C het aantal kolommen is
  • R het aantal rijen is.

V ∈ [0; 1]. Hoe groter V, hoe sterker de relatie tussen variabelen. V = 0 kan worden geïnterpreteerd als onafhankelijkheid (omdat V = 0 dan en slechts dan wanneer χ2 = 0). Het belangrijkste nadeel van V is het gebrek aan precieze interpretatie. Is V = 0,6 een sterke, middelmatige of zwakke associatie?

De functie CramerV uit DescTools kan dit voor ons berekenen:

CramerV(Cars93$Type, Cars93$Origin)

## [1] 0.3890967

Nogmaals: is dit een sterke, middelmatige of zwakke associatie?

Goodman en Kruskal's lambda

Goodman en Kruskal's lambda is een voorbeeld van een maat gebaseerd op de proportionele reductie in variatie. Dit soort maten probeert R2 — de determinatiecoëfficiënt uit de lineaire regressie — na te bootsen:

  • ze nemen waarden uit [0, 1]
  • ze drukken een fractie uit van de variatie die door de onafhankelijke variabele wordt verklaard.

Laten we een kolomvariabele gebruiken als de onafhankelijke. De formule van Goodman en Kruskal's lambda is dan: $$\lambda = \frac{L - \sum_j L_j}{L},$$ waarbij:

  • Lj de som is van niet-modale frequenties in de j-de kolom,
  • L de som is van niet-modale frequenties in de "totaal"-kolom

Dit is het best uit te leggen met een voorbeeld. Laten we kijken naar Type en Origin, waarbij de laatste de onafhankelijke is:

table(Cars93$Type, Cars93$Origin)

##          
##           USA non-USA
##   Compact   7       9
##   Large    11       0
##   Midsize  10      12
##   Small     7      14
##   Sporty    8       6
##   Van       5       4

De modus van de VS-kolom is 11, dus L1 = 7 + 10 + 7 + 8 + 5 = 37. De modus van de niet-VS-kolom is 14, dus L2 = 9 + 0 + 12 + 6 + 4 = 31. De "totaal"-kolom is

rowSums(table(Cars93$Type, Cars93$Origin))

## Compact   Large Midsize   Small  Sporty     Van
##      16      11      22      21      14       9

De modus is 22, dus L = 16 + 11 + 21 + 14 + 9 = 71 en $$\lambda = \frac{71 - (37+31)}{71}=0.042$$

We kunnen in plaats daarvan de functie Lambda uit het pakket DescTools gebruiken:

Lambda(Cars93$Type, Cars93$Origin, direction='row')

## [1] 0.04225352

De parameter direction bepaalt waar (in row of in column) de afhankelijke variabele staat.

Origin verklaart slechts zo'n 4% van de variatie in Type. Let op dat lambda variatie definieert als een dichotomie tussen wel of niet behoren tot de grootste groep.

Het is het vermelden waard dat lambda nul is wanneer de modale categorie van elke kolom dezelfde is. Kijk naar deze tabel:

(lambdaTab<-cbind(c(0,100), c(49,51)))

##      [,1] [,2]
## [1,]    0   49
## [2,]  100   51

chisq.test(lambdaTab)

##
##  Pearson's Chi-squared test with Yates' continuity correction
##
## data:  lambdaTab
## X-squared = 62.279, df = 1, p-value = 2.981e-15

De modale categorie is in elke kolom dezelfde (dat is de tweede rij). Dus moet lambda nul zijn ondanks een significant en zichtbaar verband:

Lambda(lambdaTab, direction='row')

## [1] 0

Wil je meer leren over statistiek in R, volg dan de DataCamp-cursus Statistical Modeling in R (Part 1).

Onderwerpen

Gerelateerd