Кореляційна матриця Пірсона та Спірмена в R з прикладом

Двомірна кореляція в R

Двомірний зв’язок описує зв’язок або кореляцію між двома змінними в R. У цьому посібнику ми обговоримо концепцію кореляції та покажемо, як її можна використовувати для вимірювання зв’язку між будь-якими двома змінними в R.

Кореляція в програмуванні R

Існує два основні методи обчислення кореляції між двома змінними в програмуванні R:

• Pearson: Параметрична кореляція
• Сперман: Непараметрична кореляція

Кореляційна матриця Пірсона в R

Метод кореляції Пірсона зазвичай використовується як первинна перевірка зв'язку між двома змінними.

Команда коефіцієнт кореляції, що позначається як r, вимірює силу лінійний зв'язок між двома змінними x та y. Він обчислюється наступним чином:

з

• стандартне відхилення x
• — стандартне відхилення y

Кореляція коливається від -1 до 1.

• Значення r близьке або рівне 0 означає незначний або відсутній лінійний зв'язок між x та y.
• Чим ближче r до 1 або -1, тим сильніший лінійний зв'язок.

Ви можете перевірити, чи r відрізняється від нуля, за допомогою t-статистики нижче, порівнявши її з розподілом Стьюдента з n – 2 ступенями свободи:

Рангова кореляція Спірмена в R

Рангова кореляція сортує спостереження за рангом та обчислює рівень подібності між рангами. Рангова кореляція має перевагу в тому, що вона стійка до викидів і не пов'язана з розподілом даних. Рангова кореляція також є правильним вибором для порядкових змінних.

Рангова кореляція Спірмена, що позначається як rho, також проходить від -1 до 1, а значення, близькі до будь-якого з крайніх значень, вказують на сильний монотонний зв'язок. Вона обчислюється наступним чином:

Чисельник — це коваріація між рангами x та y, а знаменник — добуток їхніх стандартних відхилень.

У R обидва обчислюються за допомогою функції cor(), яка приймає три аргументи: x, y та метод.

cor(x, y, method)

Аргументи:

• x: Перший вектор
• y: Другий вектор
• метод: формула, яка використовується для обчислення кореляції. Три значення рядка:
  • «пірсон»
  • «кендалл»
  • «списоносець»

Можна додати необов’язковий аргумент, якщо вектори містять відсутнє значення: use = “complete.obs”

Ми будемо використовувати набір даних BudgetUK. Цей набір даних повідомляє про розподіл бюджету британських домогосподарств між 1980 і 1982 роками. Існує 1519 спостережень із десятьма особливостями, серед яких:

• wfood: ділитися їжею ділитися витрачати
• wfuel: розділити витрату палива
• тканина: частка бюджету на витрати на одяг
• Walc: ділити алкоголь, витрачати
• wtrans: розділити витрати на транспорт
• мати: частка витрат на інші товари
• totexp: загальні витрати домогосподарств у фунтах
• дохід: загальний чистий дохід домогосподарства
• вік: вік дом
• діти: Кількість дітей

Приклад

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 Пояснення

• Спочатку ми імпортуємо дані та переглядаємо їх за допомогою функції glimpse() із бібліотеки dplyr.
• Три домогосподарства повідомляють про дохід 500 або більше, тому фільтр (дохід < 500) видаляє їх, а кількість рядків зменшується з 1,519 до 1,516.
• Перетворення монетарної змінної в журнал є звичайною практикою. Це допомагає зменшити вплив викидів і зменшує нерівність у наборі даних.

вихід:

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

Ми можемо обчислити коефіцієнт кореляції між змінними доходу та wfood за допомогою методів «Пірсона» та «Спірмена».

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

вихід:

## [1] -0.2466986

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

вихід:

## [1] -0.2501252

Перш ніж поширювати це на кожну пару змінних, варто визначити, як слід читати окремий коефіцієнт.

Як інтерпретувати коефіцієнт кореляції

Коефіцієнт корисний лише тоді, коли ви можете сказати, що він означає. Смуги нижче є загальноприйнятими показниками, а знак зчитується окремо від сили.

Значення -0.2467, обчислене раніше між log_income та wfood, таким чином, є слабким негативним зв'язком: багатші домогосподарства витрачають дещо меншу частку свого бюджету на їжу.

До кожного коефіцієнта застосовуються три застереження.

• Кореляція не є причинно-наслідковим зв’язком. Сильне r означає, що дві змінні рухаються разом, а не те, що одна з них спричиняє іншу. Третя, невимірювана змінна часто керує обома.
• Пірсон бачить лише прямі лінії. Ідеальне U-подібне співвідношення повертає r, близьке до нуля. Завжди побудуйте графік даних, перш ніж довіряти числу.
• Розмір має більше значення, ніж значення. З 1,516 спостереженнями коефіцієнт 0.06 може бути статистично значущим і все ще практично безглуздим.

Як перевірити значущість кореляції за допомогою cor.test()

cor() повертає коефіцієнт і нічого більше. Для однієї пари, cor.test() додає p-значення та довірчий інтервал за один виклик.

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

Вихідний матеріал має чотири частини, які варто прочитати.

1. t та df: тестова статистика та її ступені свободи, n – 2.
2. р-значення: ймовірність побачити такий великий коефіцієнт, якби справжня кореляція була нульовою.
3. 95-відсотковий довірчий інтервал: правдоподібний діапазон для справжньої кореляції. Якщо він виключає нуль, зв'язок є значущим на цьому рівні.
4. зразок кошторису: сам коефіцієнт, ідентичний тому, що повертає cor().

Та сама функція запускає тести на основі рангу, змінюючи один аргумент:

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

Коли який використовувати. Використовуйте cor.test(), коли досліджуєте одну конкретну пару, оскільки він дає довірчий інтервал, який rcorr() пропускає. Використовуйте rcorr() з Hmisc, як показано вище, коли вам потрібні p-значення для всієї матриці одночасно. Зверніть увагу, що тестування багатьох пар завищує рівень хибнопозитивних результатів, тому відкоригуйте p-значення за допомогою p.adjust(p_value, method = “BH”), перш ніж робити висновки з великої матриці.

Кореляційна матриця в R

Двовимірна кореляція – це гарний початок, але багатовимірний підхід дає ширшу картину. кореляційна матриця являє собою квадратну таблицю, що містить попарну кореляцію кожної змінної відносно кожної іншої.

Функція cor() повертає кореляційну матрицю. Єдина відмінність від двофакторної кореляції полягає в тому, що нам не потрібно вказувати, які саме змінні. За замовчуванням R обчислює кореляцію між усіма змінними.

Кореляцію неможливо обчислити для фактора, тому перед передачею фрейму даних до cor() слід видалити кожен категоріальний стовпець.

Кореляційна матриця є симетричною, що означає, що значення над діагоналлю мають ті самі значення, що й значення нижче. Більш наочно показати половину матриці.

children_fac виключається, оскільки cor() не може працювати з фактором.

# 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 Пояснення

• cor(дані[, 1:9])Обчисліть матрицю кореляції для дев'яти числових стовпців
• круглий(…, 2)Округліть кожен коефіцієнт до двох знаків під комою
• as.dist()Вивести лише нижній трикутник, оскільки матриця симетрична

вихід:

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

Рівень значущості

Сам по собі коефіцієнт не говорить про те, чи є зв'язок статистично достовірним. Функція rcorr() з бібліотеки Hmisc повертає p-значення для кожної пари. Ми можемо завантажити бібліотеку з Конда і скопіюйте код, щоб вставити його в термінал:

conda install -c r r-hmisc

Для rcorr() потрібно, щоб кадр даних зберігався як матриця. Ми можемо перетворити наші дані в матрицю перед обчисленням кореляційної матриці з p-значенням.

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

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

Об’єкт списку mat_2 містить три елементи:

• r: Вихід кореляційної матриці
• n: Номер спостереження
• P: p-значення

Нас цікавить третій елемент, p-value. Зазвичай кореляційну матрицю показують із значенням p замість коефіцієнта кореляції.

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

Code Пояснення

• mat_2[[“P”]]: p-значення зберігаються в елементі під назвою P
• round(mat_2[[“P”]], 3): Округліть елементи з трьома цифрами

вихід:

           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

Візуалізація кореляційної матриці в R

Теплова карта – це ще один спосіб зчитування матриці кореляції. Бібліотека GGally розширює ggplot2 та встановлюється з CRAN, а не з conda:

install.packages("GGally")

Бібліотека містить різні функції для показу зведеної статистики, наприклад кореляції та розподілу всіх змінних у матриця.

Функція ggcorr() має багато аргументів. Ми представимо лише ті аргументи, які використовуватимемо в підручнику:

Функція ggcorr

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

аргументи:

• df: використаний набір даних
• метод: Формула для обчислення кореляції. За замовчуванням обчислюються попарно та Пірсона
• nbreaks: повертає категорійний діапазон для забарвлення коефіцієнтів. За замовчуванням немає розривів і градієнт кольору безперервний
• цифр: Округліть коефіцієнт кореляції. За замовчуванням встановлено 2
• низький: керування нижнім рівнем забарвлення
• в середині: контроль середнього рівня кольору
• висока: контроль високого рівня забарвлення
• геом: керування формою геометричного аргументу. За замовчуванням «плитка»
• етикетка: логічне значення. Відображати чи ні мітку. За замовчуванням встановлено `FALSE`

Базова теплова карта

Найпростіший малюнок пакету – теплова карта. Легенда графіка показує колір градієнта від – 1 до 1, де гарячий колір вказує на сильну позитивну кореляцію, а холодний колір – на негативну кореляцію.

library(GGally)
ggcorr(data)

Code Пояснення

• ggcorr(дані): потрібен лише один аргумент, яким є ім’я кадру даних. Змінні рівня факторів не включені в графік.

вихід:

Додавання контролю до теплової карти

Ми можемо додати більше елементів керування до графіка:

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

Code Пояснення

• nbreaks=6: розбийте легенду за допомогою 6 рангів.
• низький = "сталевий синій": використовуйте світліші кольори для негативної кореляції
• mid = «білий»: використовуйте білі кольори для кореляції середніх діапазонів
• високий = «темний»: використовуйте темні кольори для позитивної кореляції
• geom = «коло»: використовуйте коло як форму вікон на тепловій карті. Розмір кола пропорційний абсолютному значенню кореляції.

вихід:

Додавання мітки до теплової карти

GGally дозволяє нам додавати мітку всередині вікон:

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

Code Пояснення

• мітка = TRUE: Додайте значення коефіцієнтів кореляції всередину теплової карти.
• колір = "сірий50": Виберіть колір, тобто сірий
• розмір_мітки = 3: Встановіть розмір етикетки рівним 3

вихід:

Функція ggpairs

Бібліотека GGally також надає функцію ggpairs(), яка повертає матрицю графіків. Для k вибраних змінних результатом є сітка розміром ak x k: діагональ показує розподіл кожної змінної, тоді як панелі над та під діагоналлю можуть містити різні обчислення. Синтаксис такий:

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

аргументи:

• df: використаний набір даних
• стовпців: виберіть стовпці, щоб намалювати графік
• назву: Додайте назву
• верхній: Керуйте полями над діагоналлю графіка. Потрібно вказати тип обчислень або графіка, який потрібно повернути. Якщо continuous = “cor”, ми просимо R обчислити кореляцію. Зверніть увагу, що аргумент має бути списком. Доступні й інші аргументи; див. Документація GGally для отримання додаткової інформації.
• знизити: Керуйте прямокутниками під діагоналлю.
• картаping: вказує на естетику графіка. Наприклад, ми можемо обчислити графік для різних груп.

Двовимірний аналіз з ggpair з grouping

Наступний графік відображає три відомості:

• Кореляційна матриця між log_totex, log_income, віком і змінною wtrans, згрупована за тим, чи є в домогосподарстві дитина чи ні.
• Побудуйте графік розподілу кожної змінної за групами
• Відобразіть діаграму розсіювання з тенденцією за групами

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 Пояснення

• columns = c(“log_totexp”, “log_income”, “age”, “wtrans”): виберіть змінні для відображення на графіку
• title = «Двофакторний аналіз витрат доходів британського домогосподарства»: додайте назву
• верхній = список(): контроль верхньої частини графіка. Тобто вище діагоналі
• безперервний = wrap(“cor”, size = 3)): обчисліть коефіцієнт кореляції. Ми загортаємо аргумент безперервно всередину функції wrap(), щоб контролювати естетику графіка (тобто розмір = 3) -lower = list(): контролювати нижню частину графіка. Тобто нижче діагоналі.
• безперервний = обертання (“гладкий”, альфа = 0.3, розмір = 0.1): додайте діаграму розсіювання з лінійним трендом. Ми загортаємо аргумент безперервно всередину функції wrap(), щоб контролювати естетику графіка (тобто розмір=0.1, альфа=0.3)
• картаping = aes(колір = дитячий_факт)Розділити кожну панель за параметром children_fac, упорядкованим фактором з позначкою «Ні» для домогосподарств без дітей та «Так» для домогосподарств з дітьми.

вихід:

Двовимірний аналіз з ggpair з частковою групоюping

Графік нижче трохи відрізняється. Ми змінюємо положення картиping всередині верхнього аргументу.

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 Пояснення

• Точно такий самий код, як у попередньому прикладі, за винятком:
• картаping = aes(color = children_fac): Перемістити список у upper = list(). Нам потрібно, щоб обчислення було укладено за групами лише у верхній частині графіка.

вихід:

Кореляція в R: ключові висновки та довідник функцій

• Двомірний зв’язок описує зв’язок або кореляцію між двома змінними в R.
• Існує два основні методи обчислення кореляції між двома змінними R Програмування: Pearson & Spearman.
• Метод кореляції Пірсона зазвичай використовується як первинна перевірка зв'язку між двома змінними.
• Рангова кореляція сортує спостереження за рангом і обчислює рівень подібності між рангом.
• Рангова кореляція Спірмена коливається від -1 до 1, а значення поблизу будь-якого з екстремумів вказують на сильний монотонний зв'язок.
• Кореляційна матриця — це квадратна таблиця, яка містить попарну кореляцію кожної змінної.
• p-значення показує, чи спостережувана кореляція статистично відрізняється від нуля.

Кожна функція кореляції, що використовується в цьому посібнику, перелічена нижче:

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

Поширені запитання