Pearson & Spearman korrelationsmatris i R med exempel

Bivariat korrelation i R

En bivariat relation beskriver ett samband -eller korrelation- mellan två variabler i R. I denna handledning kommer vi att diskutera begreppet korrelation och visa hur det kan användas för att mäta sambandet mellan två valfria variabler i R.

Korrelation i R-programmering

Det finns två primära metoder för att beräkna korrelationen mellan två variabler i R-programmering:

• Pearson: Parametrisk korrelation
• Spearman: Icke-parametrisk korrelation

Pearson Correlation Matrix i R

Pearson-korrelationsmetoden används vanligtvis som en primär kontroll för sambandet mellan två variabler.

Ocuco-landskapet korrelationskoefficient, skriven r, mäter styrkan hos linjär sambandet mellan två variabler x och y. Det beräknas enligt följande:

med

• är standardavvikelsen för x
• är standardavvikelsen för y

Korrelationen sträcker sig mellan -1 och 1.

• Ett värde på r nära eller lika med 0 antyder litet eller inget linjärt samband mellan x och y.
• Ju närmare r kommer 1 eller -1, desto starkare är det linjära sambandet.

Du kan testa om r skiljer sig från noll med t-statistiken nedan, och jämföra den med Studentfördelningen med n – 2 frihetsgrader:

Spearman Rank Correlation i R

En rangkorrelation sorterar observationerna efter rang och beräknar likhetsnivån mellan rangordningarna. En rangkorrelation har fördelen att den är robust mot extremvärden och inte är kopplad till datafördelningen. En rangkorrelation är också rätt val för ordinala variabler.

Spearmans rangkorrelation, skriven rho, går också från -1 till 1, och värden nära endera extremen indikerar ett starkt monotont samband. Den beräknas enligt följande:

Täljaren är kovariansen mellan raderna x och y, och nämnaren är produkten av deras standardavvikelser.

I R beräknas båda med cor()-funktionen, som tar tre argument: x, y och method.

cor(x, y, method)

Argument:

• x: Första vektorn
• y: Andra vektorn
• metod: Formeln som används för att beräkna korrelationen. Tre strängvärden:
  • "pearson"
  • "kendall"
  • "spearman"

Ett valfritt argument kan läggas till om vektorerna innehåller saknat värde: use = "complete.obs"

Vi kommer att använda datauppsättningen BudgetUK. Denna datauppsättning rapporterar budgettilldelningen för brittiska hushåll mellan 1980 och 1982. Det finns 1519 observationer med tio funktioner, bland dem:

• wfood: dela mat dela spendera
• wbränsle: dela bränsleförbrukning
• wcloth: budgetandel för utgifter för kläder
• valk: dela alkoholutgifter
• wtrans: dela transportkostnader
• wother: andel av andra varor som spenderas
• totexp: hushållens totala utgifter i pund
• intäkter: hushållens totala nettoinkomst
• ålder: hushållets ålder
• barn: Antal barn

Exempelvis

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 Förklaring

• Vi importerar först data och tar en titt med funktionen glimpse() från dplyr-biblioteket.
• Tre hushåll rapporterar en inkomst på 500 eller mer, så filter (inkomst < 500) tar bort dem och radantalet sjunker från 1 519 till 1 516.
• Det är vanligt att konvertera en monetär variabel i log. Det hjälper till att minska effekten av extremvärden och minskar skevheten i datamängden.

Produktion:

## 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, ...

Vi kan beräkna korrelationskoefficienten mellan inkomst- och wfood-variabler med metoderna "pearson" och "spearman".

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

Produktion:

## [1] -0.2466986

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

Produktion:

## [1] -0.2501252

Innan detta utvidgas till varje variabelpar är det värt att fastställa hur en enda koefficient ska läsas.

Hur man tolkar en korrelationskoefficient

En koefficient är bara användbar när du kan säga vad den betyder. Banden nedan är den konventionella avläsningen, och tecknet avläses separat från styrkan.

Värdet -0.2467 som beräknades tidigare mellan log_income och wfood är därför ett svagt negativt samband: rikare hushåll spenderar en något mindre andel av sin budget på mat.

Tre försiktighetsåtgärder gäller för varje koefficient.

• Korrelation är inte orsakssamband. Ett starkt r säger att de två variablerna rör sig tillsammans, aldrig att den ena orsakar den andra. En tredje, ouppmätt variabel driver ofta båda.
• Pearson ser bara raka linjer. Ett perfekt U-format förhållande returnerar ett r nära noll. Rita alltid in data innan du litar på talet.
• Storleken är viktigare än betydelsen. Med 1 516 observationer kan en koefficient på 0.06 vara statistiskt signifikant och fortfarande praktiskt taget meningslös.

Hur man testar korrelationssignifikans med cor.test()

cor() returnerar koefficienten och inget annat. För ett enda par adderar cor.test() p-värdet och ett konfidensintervall i ett anrop.

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

Resultatet har fyra läsvärda delar.

1. t och df: teststatistiken och dess frihetsgrader, n – 2.
2. p-värdesannolikheten att se en så stor koefficient om den verkliga korrelationen var noll.
3. 95 procents konfidensintervall: det rimliga intervallet för den sanna korrelationen. Om den utesluter noll är sambandet signifikant på den nivån.
4. exempeluppskattningsjälva koefficienten, identisk med vad cor() returnerar.

Samma funktion kör de rangbaserade testerna genom att ändra ett argument:

# 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")

När man ska använda vilket. Använd cor.test() när du undersöker ett specifikt par, eftersom det ger det konfidensintervall som rcorr() utelämnar. Använd rcorr() från Hmisc, som visas ovan, när du behöver p-värden för en hel matris samtidigt. Observera att testning av många par ökar andelen falskt positiva resultat, så justera p-värdena med p.adjust(p_value, method = "BH") innan du drar slutsatser från en stor matris.

Korrelationsmatris i R

En bivariat korrelation är en bra början, men en multivariat synvinkel ger en bredare bild. korrelationsmatris är en kvadratisk tabell som håller den parvisa korrelationen för varje variabel mot varandra.

Cor()-funktionen returnerar en korrelationsmatris. Den enda skillnaden med den bivariata korrelationen är att vi inte behöver specificera vilka variabler. Som standard beräknar R korrelationen mellan alla variabler.

En korrelation kan inte beräknas för en faktor, så ta bort varje kategorisk kolumn innan dataframen skickas till cor().

En korrelationsmatris är symmetrisk vilket innebär att värdena ovanför diagonalen har samma värden som den nedanför. Det är mer visuellt att visa hälften av matrisen.

children_fac är uteslutet eftersom cor() inte kan operera på en faktor.

# 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 Förklaring

• kor(data[, 1:9])Beräkna korrelationsmatrisen på de nio numeriska kolumnerna
• runda(…, 2)Avrunda varje koefficient till två decimaler
• as.dist()Skriv endast ut den nedre triangeln, eftersom matrisen är symmetrisk

Produktion:

##            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

Signifikansnivå

En koefficient i sig säger inte om sambandet är statistiskt tillförlitligt. Funktionen rcorr() från Hmisc-biblioteket returnerar p-värdet för varje par. Vi kan ladda ner biblioteket från Conda och kopiera koden för att klistra in den i terminalen:

conda install -c r r-hmisc

rcorr() kräver att en dataram lagras som en matris. Vi kan konvertera våra data till en matris innan vi beräknar korrelationsmatrisen med p-värdet.

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

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

Listobjektet mat_2 innehåller tre element:

• r: Utdata från korrelationsmatrisen
• n: Antal observationer
• P: p-värde

Vi är intresserade av det tredje elementet, p-värdet. Det är vanligt att visa korrelationsmatrisen med p-värdet istället för korrelationskoefficienten.

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

Code Förklaring

• mat_2[[“P”]]: P-värdena lagras i elementet som kallas P
• round(mat_2[[“P”]], 3): Avrunda elementen med tre siffror

Produktion:

           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

Visualisera korrelationsmatris i R

En värmekarta är ett annat sätt att läsa en korrelationsmatris. GGally-biblioteket utökar ggplot2 och installeras från CRAN snarare än conda:

install.packages("GGally")

Biblioteket innehåller olika funktioner för att visa den sammanfattande statistiken såsom korrelationen och fördelningen av alla variabler i en matris.

Funktionen ggcorr() har många argument. Vi kommer endast att introducera de argument vi kommer att använda i handledningen:

ggcorr-funktionen

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

argument:

• df: Dataset används
• metod: Formel för att beräkna korrelationen. Som standard beräknas parvis och Pearson
• nbreaks: Returnera ett kategoriskt område för färgen på koefficienterna. Som standard är det inget avbrott och färggradienten är kontinuerlig
• siffror: Avrunda korrelationskoefficienten. Som standard ställs in på 2
• låg: Styr den lägre nivån av färgningen
• mitten: Styr mellannivån för färgningen
• hög: Kontrollera den höga färgnivån
• geom: Styr formen på det geometriska argumentet. Som standard, "bricka"
• etikett: Booleskt värde. Visa eller inte etiketten. Som standard ställs in på "FALSE".

Grundläggande värmekarta

Den mest grundläggande handlingen i paketet är en värmekarta. Förklaringen av grafen visar en gradientfärg från – 1 till 1, med varm färg som indikerar stark positiv korrelation och kall färg, en negativ korrelation.

library(GGally)
ggcorr(data)

Code Förklaring

• ggcorr(data): Endast ett argument behövs, vilket är dataramens namn. Faktornivåvariabler ingår inte i diagrammet.

Produktion:

Lägga till kontroll till värmekartan

Vi kan lägga till fler kontroller i diagrammet:

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

Code Förklaring

• nbreaks=6: bryt legenden med 6 led.
• låg = "stålblå": Använd ljusare färger för negativ korrelation
• mitt = "vit": Använd vita färger för mellanintervallkorrelation
• hög = "mörkröd": Använd mörka färger för positiv korrelation
• geom = "cirkel": Använd cirkel som formen på fönstren i värmekartan. Cirkelns storlek är proportionell mot korrelationens absoluta värde.

Produktion:

Lägger till etikett till värmekartan

GGally låter oss lägga till en etikett i fönstren:

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

Code Förklaring

• etikett = SANT: Lägg till värdena för korrelationskoefficienterna inuti värmekartan.
• färg = "grå50": Välj färg, dvs grått
• label_size = 3: Ställ in storleken på etiketten lika med 3

Produktion:

ggpairs-funktionen

GGally-biblioteket tillhandahåller även ggpairs(), som returnerar en matris av diagram. För k valda variabler blir resultatet ak gånger k rutnät: diagonalen visar fördelningen av varje variabel, medan panelerna ovanför och under diagonalen kan utföra olika beräkningar. Syntaxen är:

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

argument:

• df: Dataset används
• kolonner: Välj kolumnerna för att rita plotten
• rubricerade: Inkludera en titel
• övreKontrollera rutorna ovanför diagrammets diagonal. Behöver ange vilken typ av beräkningar eller graf som ska returneras. Om kontinuerlig = "kor" ber vi R att beräkna korrelationen. Observera att argumentet måste vara en lista. Andra argument är tillgängliga; se GGally-dokumentation för mer information.
• lägre: Styr rutorna under diagonalen.
• kartaping: Indikerar grafens estetik. Vi kan till exempel beräkna grafen för olika grupper.

Bivariat analys med ggpair med grouping

Nästa graf visar tre information:

• Korrelationsmatrisen mellan log_totexp, log_income, age och wtrans variabel grupperad efter om hushållet har ett barn eller inte.
• Rita fördelningen av varje variabel per grupp
• Visa spridningsdiagrammet med trenden för grupp

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 Förklaring

• kolumner = c(“log_totexp”, “log_income”, “age”, “wtrans”): Välj de variabler som ska visas i grafen
• title = "Bivariat analys av det brittiska hushållets inkomstutgifter": Lägg till en titel
• övre = lista(): Styr den övre delen av grafen. Dvs ovanför diagonalen
• kontinuerlig = wrap(“cor”, storlek = 3)): Beräkna korrelationskoefficienten. Vi lindar argumentet kontinuerligt inuti wrap()-funktionen för att kontrollera grafens estetik (dvs storlek = 3) -lower = list(): Styr den nedre delen av grafen. Dvs. Under diagonalen.
• kontinuerlig = wrap(“slät”,alfa = 0.3,storlek=0.1): Lägg till ett spridningsdiagram med en linjär trend. Vi lindar argumentet kontinuerligt inuti wrap()-funktionen för att kontrollera grafens estetik (dvs storlek=0.1, alfa=0.3)
• kartaping = aes(färg = barn_fac)Dela upp varje panel med children_fac, den ordnade faktorn märkt "Nej" för hushåll utan barn och "Ja" för hushåll med barn

Produktion:

Bivariat analys med ggpair med partiell grouping

Diagrammet nedan är lite annorlunda. Vi ändrar kartans position.ping inuti det övre argumentet.

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 Förklaring

• Exakt samma kod som föregående exempel förutom:
• kartaping = aes(color = children_fac): Flyttar listan i upper = list(). Vi vill bara ha beräkningen staplad gruppvis i den övre delen av grafen.

Produktion:

Korrelation i R: Viktiga slutsatser och funktionsreferens

• Ett bivariat samband beskriver ett samband -eller korrelation- mellan två variabler i R.
• Det finns två primära metoder för att beräkna korrelationen mellan två variabler i R-programmering: Pearson & Spearman.
• Pearson-korrelationsmetoden används vanligtvis som en primär kontroll för sambandet mellan två variabler.
• En rangkorrelation sorterar observationerna efter rangordning och beräknar graden av likhet mellan rangordningen.
• Spearmans rangkorrelation går från -1 till 1, och värden nära endera extremen indikerar ett starkt monotont samband.
• En korrelationsmatris är en kvadratisk tabell som innehåller den parvisa korrelationen för varje variabel.
• Ett p-värde visar om en observerad korrelation är statistiskt särskiljbar från noll.

Varje korrelationsfunktion som används i den här handledningen listas nedan:

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"))

Vanliga frågor