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

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

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

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

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

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

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

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

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

з

• , тобто стандартне відхилення від
• , тобто стандартне відхилення від

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

• Значення, що дорівнює або наближається до 0, передбачає невеликий або відсутність лінійного зв’язку між і.
• Навпаки, чим ближче до 1 або -1, тим сильніший лінійний зв’язок.

Ми можемо обчислити t-тест наступним чином і перевірити таблицю розподілу зі ступенем свободи, що дорівнює:

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

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

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

із зазначеними коваріаціями між рангом і . Знаменник обчислює стандартні відхилення.

У R ми можемо використовувати функцію cor(). Він приймає три аргументи, і метод.

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)

Пояснення коду

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

вихід:

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

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

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

Функція 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

Пояснення коду

• cor(дані): відображення кореляційної матриці
• раунд(дані, 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

Пояснення коду

• 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. Наразі він недоступний у бібліотеці 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)

Пояснення коду

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

вихід:

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

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

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

Пояснення коду

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

вихід:

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

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

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

Пояснення коду

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

вихід:

Функція ggpairs

Нарешті, ми представляємо ще одну функцію з бібліотеки GGaly. Ggpair. Він створює графік у матричному форматі. Ми можемо відобразити три види обчислень на одному графіку. Матриця – це розмірність, що дорівнює кількості спостережень. У верхній/нижній частині відображаються вікна та по діагоналі. Ми можемо контролювати, яку інформацію ми хочемо показати в кожній частині матриці. Формула для ggpair така:

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

аргументи:

• df: використаний набір даних
• стовпців: виберіть стовпці, щоб намалювати графік
• назву: Додайте назву
• верхній: контролюйте прямокутники над діагоналлю ділянки. Необхідно вказати тип обчислень або графіка для повернення. Якщо безперервний = “cor”, ми просимо R обчислити кореляцію. Зауважте, що аргумент має бути списком. Можна використовувати інші аргументи, див. [vignette](“http://ggobi.github.io/ggally/#custom_functions”) для отримання додаткової інформації.
• Опустіть: Керуйте прямокутниками під діагоналлю.
• Карт: вказує на естетику графіка. Наприклад, ми можемо обчислити графік для різних груп.

Двомірний аналіз за допомогою ggpair з групуванням

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

• Кореляційна матриця між 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))

Пояснення коду

• 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)
• відображення = aes(колір = children_fac): ми хочемо, щоб кожна частина графіка була складена за змінною children_fac, яка є категоріальною змінною, яка приймає значення 1, якщо в домогосподарстві немає дітей, і 2 в іншому випадку

вихід:

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

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

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

Пояснення коду

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

вихід:

Підсумки

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

Ми можемо узагальнити всі кореляційні функції в R у таблиці нижче:

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)

cf code below