Pearson- en Spearman-correlatiematrix in R met voorbeeld

Bivariate correlatie in R

Een bivariate relatie beschrijft een relatie (of correlatie) tussen twee variabelen in R. In deze tutorial bespreken we het concept van correlatie en laten we zien hoe het kan worden gebruikt om de relatie tussen twee willekeurige variabelen in R te meten.

Correlatie in R-programmering

Er zijn twee primaire methoden om de correlatie tussen twee variabelen in R-programmering te berekenen:

• Pearson: Parametrische correlatie
• Spearman: Niet-parametrische correlatie

Pearson-correlatiematrix in R

De Pearson-correlatiemethode wordt meestal gebruikt als primaire controle voor de relatie tussen twee variabelen.

De Correlatiecoëfficiënt, geschreven als r, meet de sterkte van de lineair De relatie tussen twee variabelen x en y wordt als volgt berekend:

with

• is de standaardafwijking van x
• is de standaardafwijking van y

De correlatie varieert tussen -1 en 1.

• Een waarde van r die dicht bij of gelijk aan 0 ligt, impliceert een gering of afwezig lineair verband tussen x en y.
• Hoe dichter r bij 1 of -1 komt, hoe sterker het lineaire verband.

Je kunt met behulp van de onderstaande t-statistiek testen of r verschilt van nul, door deze te vergelijken met de Student-verdeling met n – 2 vrijheidsgraden:

Spearman-rangcorrelatie in R

Een rangcorrelatie sorteert de waarnemingen op rang en berekent de mate van gelijkenis tussen de rangen. Een rangcorrelatie heeft als voordeel dat deze robuust is tegen uitschieters en niet afhankelijk is van de verdeling van de gegevens. Een rangcorrelatie is ook de juiste keuze voor ordinale variabelen.

De rangcorrelatiecoëfficiënt van Spearman, afgekort rho, loopt ook van -1 tot 1, en waarden dicht bij een van beide uitersten duiden op een sterke monotone relatie. Deze wordt als volgt berekend:

De teller is de covariantie tussen de rangen van x en y, en de noemer is het product van hun standaarddeviaties.

In R worden beide berekend met de functie cor(), die drie argumenten accepteert: x, y en methode.

cor(x, y, method)

argumenten:

• x: Eerste vector
• y: Tweede vector
• methode: De formule die wordt gebruikt om de correlatie te berekenen. Drie stringwaarden:
  • “pearson”
  • “kendall”
  • “speerman”

Er kan een optioneel argument worden toegevoegd als de vectoren een ontbrekende waarde bevatten: use = “complete.obs”

We zullen de BudgetUK-dataset gebruiken. Deze dataset rapporteert de budgettoewijzing van Britse huishoudens tussen 1980 en 1982. Er zijn 1519 observaties met tien kenmerken, waaronder:

• w eten: deel de uitgaven voor voedsel
• brandstof: brandstofuitgaven delen
• kleding: budgetaandeel voor kledinguitgaven
• loop: deel alcoholuitgaven
• wtrans: deel transportuitgaven
• woor: aandeel van de overige goederenuitgaven
• totexp: totale uitgaven van huishoudens in ponden
• inkomen: totaal netto huishoudinkomen
• leeftijd: leeftijd van het huishouden
• kinderen: Aantal kinderen

Voorbeeld

library(dplyr)
PATH <- "https://raw.githubusercontent.com/guru99-edu/R-Programming/master/british_household.csv"
data <- read.csv(PATH) %>%
    filter(income < 500) %>%
    mutate(log_income = log(income),
           log_totexp = log(totexp),
           children_fac = factor(children, order = TRUE, labels = c("No", "Yes"))) %>%
    select(-c(X, X.1, children, totexp, income))
glimpse(data)

Code Uitleg

• We importeren eerst de gegevens en kijken met de functie glimp() uit de dplyr-bibliotheek.
• Drie huishoudens geven een inkomen van 500 of meer op, dus filter(inkomen < 500) verwijdert deze en het aantal rijen daalt van 1,519 naar 1,516.
• Het is gebruikelijk om een ​​monetaire variabele in log om te zetten. Het helpt de impact van uitschieters te verminderen en de scheefheid in de dataset te verminderen.

Output:

## Observations: 1,516
## Variables: 10
## $ wfood        <dbl> 0.4272, 0.3739, 0.1941, 0.4438, 0.3331, 0.3752, 0...
## $ wfuel        <dbl> 0.1342, 0.1686, 0.4056, 0.1258, 0.0824, 0.0481, 0...
## $ wcloth       <dbl> 0.0000, 0.0091, 0.0012, 0.0539, 0.0399, 0.1170, 0...
## $ walc         <dbl> 0.0106, 0.0825, 0.0513, 0.0397, 0.1571, 0.0210, 0...
## $ wtrans       <dbl> 0.1458, 0.1215, 0.2063, 0.0652, 0.2403, 0.0955, 0...
## $ wother       <dbl> 0.2822, 0.2444, 0.1415, 0.2716, 0.1473, 0.3431, 0...
## $ age          <int> 25, 39, 47, 33, 31, 24, 46, 25, 30, 41, 48, 24, 2...
## $ log_income   <dbl> 4.867534, 5.010635, 5.438079, 4.605170, 4.605170,...
## $ log_totexp   <dbl> 3.912023, 4.499810, 5.192957, 4.382027, 4.499810,...
## $ children_fac <ord> Yes, Yes, Yes, Yes, No, No, No, No, No, No, Yes, ...

We kunnen de correlatiecoëfficiënt tussen inkomens- en wfood-variabelen berekenen met de ‘pearson’- en ‘spearman’-methoden.

cor(data$log_income, data$wfood, method = "pearson")

Output:

## [1] -0.2466986

cor(data$log_income, data$wfood, method = "spearman")

Output:

## [1] -0.2501252

Voordat we dit uitbreiden naar elk paar variabelen, is het de moeite waard om vast te stellen hoe een enkele coëfficiënt moet worden geïnterpreteerd.

Hoe interpreteer je een correlatiecoëfficiënt?

Een coëfficiënt is pas nuttig als je kunt uitleggen wat hij betekent. De onderstaande banden geven de conventionele aflezing weer, waarbij het teken los van de sterkte wordt afgelezen.

De eerder berekende waarde van -0.2467 tussen logaritme van het inkomen en voedseluitgaven duidt dus op een zwakke negatieve correlatie: rijkere huishoudens besteden een iets kleiner deel van hun budget aan voedsel.

Voor elke coëfficiënt gelden drie waarschuwingen.

• Correlatie is geen causaliteit. Een sterke correlatiecoëfficiënt r geeft aan dat de twee variabelen samen bewegen, niet dat de ene de andere veroorzaakt. Vaak is er een derde, niet-gemeten variabele die beide beïnvloedt.
• Pearson ziet alleen rechte lijnen. Een perfecte U-vormige relatie levert een r-waarde op die dicht bij nul ligt. Plot de gegevens altijd voordat je op het getal vertrouwt.
• Grootte is belangrijker dan betekenis. Met 1,516 waarnemingen kan een coëfficiënt van 0.06 statistisch significant zijn, maar in de praktijk toch betekenisloos.

Hoe test je de significantie van een correlatie met cor.test()?

`cor()` retourneert de coëfficiënt en verder niets. Voor een enkel paar voegt `cor.test()` de p-waarde en een betrouwbaarheidsinterval in één aanroep toe.

cor.test(data$log_income, data$wfood, method = "pearson")

De output bestaat uit vier delen die de moeite waard zijn om te lezen.

1. t en df: de toetsingsstatistiek en het aantal vrijheidsgraden, n – 2.
2. p-waarde: de kans om een ​​coëfficiënt van deze omvang te zien als de werkelijke correlatie nul zou zijn.
3. 95 procent betrouwbaarheidsinterval: het plausibele bereik voor de werkelijke correlatie. Als het nul uitsluit, is de relatie significant op dat niveau.
4. steekproefschatting: de coëfficiënt zelf, identiek aan wat cor() retourneert.

Dezelfde functie voert de rangorde-gebaseerde tests uit door één argument te wijzigen:

# Spearman rank correlation with a p-value
cor.test(data$log_income, data$wfood, method = "spearman")

# One-sided test: is the correlation greater than zero?
cor.test(data$log_income, data$wfood, alternative = "greater")

Wanneer welke te gebruiken. Gebruik `cor.test()` wanneer u één specifiek paar onderzoekt, omdat deze functie het betrouwbaarheidsinterval geeft dat `rcorr()` niet weergeeft. Gebruik `rcorr()` uit Hmisc, zoals hierboven weergegeven, wanneer u p-waarden voor een hele matrix tegelijk nodig hebt. Houd er rekening mee dat het testen van veel paren de kans op valse positieven vergroot, dus corrigeer de p-waarden met `p.adjust(p_value, method = “BH”)` voordat u conclusies trekt uit een grote matrix.

Correlatiematrix in R

Een bivariate correlatieanalyse is een goed begin, maar een multivariate analyse geeft een completer beeld. correlatiematrix is een vierkante tabel die de paarsgewijze correlatie van elke variabele ten opzichte van elke andere variabele weergeeft.

De functie cor() retourneert een correlatiematrix. Het enige verschil met de bivariate correlatie is dat we niet hoeven te specificeren welke variabelen. Standaard berekent R de correlatie tussen alle variabelen.

Voor een factor kan geen correlatie worden berekend, dus verwijder alle categorische kolommen voordat je het dataframe aan cor() doorgeeft.

Een correlatiematrix is ​​symmetrisch, wat betekent dat de waarden boven de diagonaal dezelfde waarden hebben als die eronder. Het is visueler om de helft van de matrix weer te geven.

children_fac is uitgesloten omdat cor() niet op een factor kan werken.

# the last column of data is a factor level. We don't include it in the code
mat_1 <-as.dist(round(cor(data[,1:9]),2))
mat_1

Code Uitleg

• cor(data[, 1:9])Bereken de correlatiematrix op basis van de negen numerieke kolommen.
• ronde(…, 2)Rond elke coëfficiënt af op twee decimalen.
• as.dist(): Print alleen de onderste driehoek, aangezien de matrix symmetrisch is.

Output:

##            wfood wfuel wcloth  walc wtrans wother   age log_income
## wfuel       0.11                                                  
## wcloth     -0.33 -0.25                                            
## walc       -0.12 -0.13  -0.09                                     
## wtrans     -0.34 -0.16  -0.19 -0.22                               
## wother     -0.35 -0.14  -0.22 -0.12  -0.29                        
## age         0.02 -0.05   0.04 -0.14   0.03   0.02                 
## log_income -0.25 -0.12   0.10  0.04   0.06   0.13  0.23           
## log_totexp -0.50 -0.36   0.34  0.12   0.15   0.15  0.21       0.49

Mate van belangrijkheid

Een coëfficiënt op zich zegt niets over de statistische betrouwbaarheid van de relatie. De functie rcorr() uit de Hmisc-bibliotheek geeft de p-waarde voor elk paar terug. We kunnen de bibliotheek downloaden via Conda en kopieer de code om deze in de terminal te plakken:

conda install -c r r-hmisc

De rcorr() vereist dat een dataframe wordt opgeslagen als een matrix. We kunnen onze gegevens eerst omzetten in een matrix om de correlatiematrix met de p-waarde te berekenen.

library("Hmisc")
data_rcorr <-as.matrix(data[, 1: 9])

mat_2 <-rcorr(data_rcorr)
# mat_2 <-rcorr(as.matrix(data)) returns the same output

Het lijstobject mat_2 bevat drie elementen:

• r: Uitvoer van de correlatiematrix
• n: Aantal observaties
• P: p-waarde

We zijn geïnteresseerd in het derde element, de p-waarde. Het is gebruikelijk om de correlatiematrix weer te geven met de p-waarde in plaats van de correlatiecoëfficiënt.

p_value <-round(mat_2[["P"]], 3)
p_value

Code Uitleg

• mat_2[[“P”]]: De p-waarden worden opgeslagen in het element genaamd P
• round(mat_2[[“P”]], 3): Rond de elementen af ​​met drie cijfers

Output:

           wfood wfuel wcloth  walc wtrans wother   age log_income log_totexp
wfood         NA 0.000  0.000 0.000  0.000  0.000 0.365      0.000          0
wfuel      0.000    NA  0.000 0.000  0.000  0.000 0.076      0.000          0
wcloth     0.000 0.000     NA 0.001  0.000  0.000 0.160      0.000          0
walc       0.000 0.000  0.001    NA  0.000  0.000 0.000      0.105          0
wtrans     0.000 0.000  0.000 0.000     NA  0.000 0.259      0.020          0
wother     0.000 0.000  0.000 0.000  0.000     NA 0.355      0.000          0
age        0.365 0.076  0.160 0.000  0.259  0.355    NA      0.000          0
log_income 0.000 0.000  0.000 0.105  0.020  0.000 0.000         NA          0
log_totexp 0.000 0.000  0.000 0.000  0.000  0.000 0.000      0.000         NA

Visualisatie van de correlatiematrix in R

Een heatmap is een andere manier om een ​​correlatiematrix te interpreteren. De GGally-bibliotheek is een uitbreiding van ggplot2 en wordt geïnstalleerd via CRAN in plaats van conda:

install.packages("GGally")

De bibliotheek bevat verschillende functies om de samenvattende statistieken weer te geven, zoals de correlatie en verdeling van alle variabelen in a Matrix.

De functie ggcorr() heeft veel argumenten. We introduceren alleen de argumenten die we in de tutorial zullen gebruiken:

De ggcorr-functie

ggcorr(df, method = c("pairwise", "pearson"),
  nbreaks = NULL, digits = 2, low = "#3B9AB2",
  mid = "#EEEEEE", high = "#F21A00",
  geom = "tile", label = FALSE,
  label_alpha = FALSE)

Argumenten:

• df: Dataset gebruikt
• methode: Formule om de correlatie te berekenen. Standaard worden pairwise en Pearson berekend
• npauzes: Retourneert een categorisch bereik voor de kleuring van de coëfficiënten. Standaard is er geen onderbreking en is het kleurverloop continu
• cijfers: Rond de correlatiecoëfficiënt af. Standaard ingesteld op 2
• lage: Regel het lagere niveau van de kleuring
• mid: Beheer het middelste niveau van de kleuring
• hoog: Controleer het hoge niveau van de kleuring
• geom: Bepaal de vorm van het geometrische argument. Standaard “tegel”
• label: Booleaanse waarde. Geef het label wel of niet weer. Standaard ingesteld op `FALSE`

Basis hittekaart

De meest basale plot van het pakket is een heatmap. De legenda van de grafiek toont een gradiëntkleur van – 1 tot 1, waarbij warme kleur een sterke positieve correlatie aangeeft en koude kleur een negatieve correlatie.

library(GGally)
ggcorr(data)

Code Uitleg

• ggcorr(gegevens): Er is slechts één argument nodig, namelijk de dataframenaam. Variabelen op factorniveau zijn niet in de grafiek opgenomen.

Output:

Controle toevoegen aan de hittekaart

We kunnen meer besturingselementen aan de grafiek toevoegen:

ggcorr(data,
    nbreaks = 6,
    low = "steelblue",
    mid = "white",
    high = "darkred",
    geom = "circle")

Code Uitleg

• npauzes=6: doorbreek de legende met 6 rangen.
• laag = “staalblauw”: Gebruik lichtere kleuren voor negatieve correlatie
• midden = “wit”: Gebruik witte kleuren voor correlatie in het middenbereik
• hoog = “donkerrood”: Gebruik donkere kleuren voor positieve correlatie
• geom = “cirkel”: Gebruik cirkel als de vorm van de vensters in de heatmap. De grootte van de cirkel is proportioneel aan de absolute waarde van de correlatie.

Output:

Label toevoegen aan de hittekaart

Met GGally kunnen we een label aan de vensters toevoegen:

ggcorr(data,
    nbreaks = 6,
    label = TRUE,
    label_size = 3,
    color = "grey50")

Code Uitleg

• etiket = WAAR: voeg de waarden van de correlatiecoëfficiënten toe binnen de hittekaart.
• kleur = “grijs50”: Kies de kleur, bijvoorbeeld grijs
• label_grootte = 3: Stel de grootte van het label in op 3

Output:

De ggpairs-functie

De GGally-bibliotheek biedt ook de functie ggpairs(), die een matrix van grafieken retourneert. Voor k geselecteerde variabelen is het resultaat een raster van ak bij k: de diagonaal toont de verdeling van elke variabele, terwijl de panelen boven en onder de diagonaal elk een andere berekening kunnen bevatten. De syntax is:

ggpairs(df, columns = 1:ncol(df), title = NULL,
    upper = list(continuous = "cor"),
    lower = list(continuous = "smooth"),
    mapping = NULL)		

Argumenten:

• df: Dataset gebruikt
• kolommen: Selecteer de kolommen om de plot te tekenen
• titel: Voeg een titel toe
• bovenste: Hiermee kunt u de vakjes boven de diagonaal van de grafiek beheren. U moet het type berekening of grafiek opgeven dat moet worden geretourneerd. Als `continue = "cor"`, vragen we R om de correlatie te berekenen. Houd er rekening mee dat het argument een lijst moet zijn. Andere argumenten zijn beschikbaar; zie de documentatie. GGally-documentatie voor meer informatie.
• te verlagen: Bestuur de vakken onder de diagonaal.
• kaartping: Geeft de esthetiek van de grafiek aan. We kunnen bijvoorbeeld de grafiek voor verschillende groepen berekenen.

Bivariate analyse met ggpair met grouping

De volgende grafiek toont drie informatie:

• De correlatiematrix tussen log_totexp, log_income, leeftijd en wtrans-variabele, gegroepeerd op basis van het feit of het huishouden een kind heeft of niet.
• Teken de verdeling van elke variabele per groep
• Geef het spreidingsdiagram weer met de trend per groep

library(ggplot2)
ggpairs(data, columns = c("log_totexp", "log_income", "age", "wtrans"), title = "Bivariate analysis of revenue expenditure by the British household", upper = list(continuous = wrap("cor",
        size = 3)),
    lower = list(continuous = wrap("smooth",
        alpha = 0.3,
        size = 0.1)),
    mapping = aes(color = children_fac))

Code Uitleg

• columns = c(“log_totexp”, “log_income”, “leeftijd”, “wtrans”): Kies de variabelen die u in de grafiek wilt weergeven
• title = “Bivariate analyse van de inkomstenuitgaven door het Britse huishouden”: Voeg een titel toe
• bovenste = lijst(): beheer het bovenste deel van de grafiek. Dat wil zeggen: boven de diagonaal
• continu = wrap(“cor”, maat = 3)): Bereken de correlatiecoëfficiënt. We plaatsen het argument continu in de functie wrap() om de esthetiek van de grafiek te controleren (dwz grootte = 3) -lower = list(): Beheer het onderste deel van de grafiek. Dat wil zeggen onder de diagonaal.
• continu = wrap(“glad”,alpha = 0.3,grootte=0.1): Voeg een spreidingsdiagram toe met een lineaire trend. We plaatsen het argument continu in de functie wrap() om de esthetiek van de grafiek te controleren (dwz grootte = 0.1, alpha = 0.3)
• kaartping = aes(kleur = children_fac)Splits elk paneel op basis van children_fac, de geordende factor met het label "Nee" voor huishoudens zonder kinderen en "Ja" voor huishoudens met kinderen.

Output:

Bivariate analyse met ggpair met gedeeltelijke groeperingping

De onderstaande grafiek is iets anders. We veranderen de positie van de kaart.ping binnen het bovenste argument.

ggpairs(data, columns = c("log_totexp", "log_income", "age", "wtrans"),
    title = "Bivariate analysis of revenue expenditure by the British household",
    upper = list(continuous = wrap("cor",
            size = 3),
        mapping = aes(color = children_fac)),
    lower = list(
        continuous = wrap("smooth",
            alpha = 0.3,
            size = 0.1))
)

Code Uitleg

• Exact dezelfde code als het vorige voorbeeld, behalve:
• kaartping = aes(color = children_fac): Verplaats de lijst naar boven = list(). We willen alleen de berekening gestapeld per groep in het bovenste deel van de grafiek.

Output:

Correlatie in R: Belangrijkste conclusies en functiereferentie

• Een bivariate relatie beschrijft een relatie – of correlatie – tussen twee variabelen in R.
• Er zijn twee primaire methoden om de correlatie tussen twee variabelen te berekenen R Programmeren: Pearson & Speerman.
• De Pearson-correlatiemethode wordt meestal gebruikt als primaire controle voor de relatie tussen twee variabelen.
• Een rangcorrelatie sorteert de waarnemingen op rang en berekent de mate van overeenkomst tussen de rangschikkingen.
• De rangcorrelatiecoëfficiënt van Spearman loopt van -1 tot 1, en waarden dicht bij een van beide uitersten duiden op een sterke monotone relatie.
• Een correlatiematrix is ​​een vierkante tabel die de paarsgewijze correlatie van elke variabele weergeeft.
• Een p-waarde geeft aan of een waargenomen correlatie statistisch significant verschilt van nul.

Alle correlatiefuncties die in deze handleiding worden gebruikt, staan ​​hieronder vermeld:

cor(dfx2, method = "pearson")

cor(dfx2, method = "spearman")

cor(df, method = "pearson")

cor(df, method = "spearman")

rcorr(as.matrix(data[,1:9]))[["P"]]

ggcorr(df)

ggpairs(df, columns = c("x1", "x2"))

Veelgestelde vragen