суббота, 10 октября 2020 г.

Статистика в R : корреляционный анализ

 Изучение регрессии в R будем на примере работы с встроенными в R данными mtcars. Посмотрим help по этому датафрейму :

Эти данные были извлечены из журнала Motor Trend США в 1974 году, и включает в себя потребление топлива и 10 аспектов автомобильного дизайна и производительности для 32 автомобилей (1973-74 моделей).


Data frame с 32 строками по 11 наблюдений (числовым) переменным. 

[, 1] mpg Miles/(US) gallon (расход топлива мили/галлон)

[, 2] cyl Number of cylinders (число цилиндров)

[, 3] disp Displacement (cu.in.) (объем двигателя)

[, 4] hp Gross horsepower (мощность в л.с.) 

[, 5] drat Rear axle ratio (передаточное число ведущего моста, передаточное число главной передачи)

[, 6] wt Weight (1000 lbs) (вес в 1000 фунтах)

[, 7] qsec 1/4 mile time (время прохождения 1/4 мили 402м в сек.)

[, 8] vs Engine (0 = V-shaped, 1 = straight) (форма двигателя)

[, 9] am Transmission (0 = automatic, 1 = manual) (вид коробки передач)

[,10] gear Number of forward gears  (количество передач)

[,11] carb Number of carburetors (количество карбюраторов)


Перед началом работы устанавливаем опции для отображения десятичных дробей :

options(scipen = 6, digits = 3)

df <- mtcars

str(df)

data.frame': 32 obs. of  14 variables:
 $ mpg     : num  10.4 10.4 13.3 14.3 14.7 15 15.2 15.2 15.5 15.8 ...
 $ cyl     : num  8 8 8 8 8 8 8 8 8 8 ...
 $ disp    : num  472 460 350 360 440 ...
 $ hp      : num  205 215 245 245 230 335 180 150 150 264 ...
 $ drat    : num  2.93 3 3.73 3.21 3.23 3.54 3.07 3.15 2.76 4.22 ...
 $ wt      : num  5.25 5.42 3.84 3.57 5.34 ...
 $ qsec    : num  18 17.8 15.4 15.8 17.4 ...
 $ vs      : num  0 0 0 0 0 0 0 0 0 0 ...
 $ am      : num  0 0 0 0 0 1 0 0 0 1 ...
 $ gear    : num  3 3 3 3 3 5 3 3 3 5 ...
 $ carb    : num  4 4 4 4 4 8 3 2 2 4 ...
 $ car name: Factor w/ 32 levels "Cadillac Fleetwood",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ mpg_z   : num  -1.61 -1.61 -1.13 -0.96 -0.89 -0.84 -0.81 -0.81 -0.76 -0.71 ...
 $ mpg_type: chr  "below" "below" "below" "below" ...

Как видно, help выдал не всю информацию, переменных 14. Что такое "car name" понятно, а вот что означает mpg_z и mpg_type не ясно и потому не будем их использовать. 
В help есть полезные преобразования показателей и мы ими воспользуемся :

df <- within(mtcars, {
  vs <- factor(vs, labels = c("V", "S"))
  am <- factor(am, labels = c("automatic", "manual"))
  cyl  <- ordered(cyl)
  gear <- ordered(gear)
  carb <- ordered(carb)
})

str(df)

'data.frame': 32 obs. of  14 variables:
 $ mpg     : num  10.4 10.4 13.3 14.3 14.7 15 15.2 15.2 15.5 15.8 ...
 $ cyl     : Ord.factor w/ 3 levels "4"<"6"<"8": 3 3 3 3 3 3 3 3 3 3 ...
 $ disp    : num  472 460 350 360 440 ...
 $ hp      : num  205 215 245 245 230 335 180 150 150 264 ...
 $ drat    : num  2.93 3 3.73 3.21 3.23 3.54 3.07 3.15 2.76 4.22 ...
 $ wt      : num  5.25 5.42 3.84 3.57 5.34 ...
 $ qsec    : num  18 17.8 15.4 15.8 17.4 ...
 $ vs      : Factor w/ 2 levels "V","S": 1 1 1 1 1 1 1 1 1 1 ...
 $ am      : Factor w/ 2 levels "automatic","manual": 1 1 1 1 1 2 1 1 1 2 ...
 $ gear    : Ord.factor w/ 3 levels "3"<"4"<"5": 1 1 1 1 1 3 1 1 1 3 ...
 $ carb    : Ord.factor w/ 6 levels "1"<"2"<"3"<"4"<..: 4 4 4 4 4 6 3 2 2 4 ...
 $ car name: Factor w/ 32 levels "Cadillac Fleetwood",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ mpg_z   : num  -1.61 -1.61 -1.13 -0.96 -0.89 -0.84 -0.81 -0.81 -0.76 -0.71 ...
 $ mpg_type: chr  "below" "below" "below" "below" ...

В этом датафрейме нет пропущенных значений, но для того, чтобы проверить методы работы с ними, введем их искусственно

df$mpg_type[5] <- NA

df$mpg_z[10] <- NA

Проверяем, есть ли пропущенные значения?

sum(!complete.cases(df))

[1] 2

Где именно пропущенные значения?

colSums(is.na(df))

  mpg  cyl   disp    hp    drat    wt     qsec       vs       am     gear     carb car name  mpg_z mpg_type
       0     0      0        0        0        0        0        0        0        0        0        0                         1            1

Выведем строки с пропущенными значениями

df[!complete.cases(df),]

                               mpg cyl disp  hp drat    wt  qsec vs        am gear carb          car name         mpg_z  mpg_type
Chrysler Imperial 14.7   8  440 230 3.23 5.345 17.42  V automatic    3    4 Chrysler Imperial  -0.89        NA
Ford Pantera L    15.8   8  351 264 4.22 3.170 14.50  V    manual    5    4    Ford Pantera L       NA          below

Цель практически любого исследования заключается в поиске взаимосвязей величин и создании основы для предсказания чего-то неизвестного на основе имеющихся данных. Решением такой задачи занимается корреляционный анализ. 

Имея дело с корреляционным анализом, мы должны принять два базовых утверждения. Во-первых, это то что наличие связи между явлениями не означает, что между ними существует причинно-следственная связь. То есть наличие корреляции ничего не говорит о том, что является причиной, а что — следствием. 
И вторая важная посылка заключается в том, что связь между явлениями можно измерить количественно.
 Однако для начала эту связь необходимо выявить, обнаружить. Если мы умеем изучаемые явления описывать количественно, то оценить наличие связей между двумя явлениями довольно просто. 
Начнем анализ с визуального представления. Лучше привести один график, чем много чисел в каких-то таблицах. 
Для визуализации связей между двумя численно выраженными величинами можно использовать простейший прием, который называется построение точечной диаграммы. 

library (ggplot2)
theme_set(theme_bw())
g_mpg_hp <- ggplot(data = df, aes(x = hp, y = mpg)) +
  geom_point() +
  labs(x = "Мощьность двигателя, л.с.", y = "Потребление топлива, мили/галлон")
g_mpg_hp



Как видно по этому графику, связь безусловно есть и она нелинейная.

В качестве меры корреляционной связи обычно используют две величины — это коэффициент корреляции и коэффициент ковариации. 

Коэффициент ковариации :

$$cov_{x,y}=\frac{\sum(x_{i}-\bar{x})*(y_{i}-\bar{y}) }{n-1}$$

Rоэффициент ковариации, оценивает степень сонаправленности отклонения двух величин от своих средних значений. Отклонение от средних принято называть вариацией а ковариация означает совместную вариацию двух сопряженных величин. В формуле x с индексом i обозначается отдельное i-тое измерение величины x. x с чертой, вот она, — это среднее значение для ряда, которое составляет величину x. То же самое в отношении y. n — это объем выборки, это количество сопряженных пар, которые составляют x и y. Этот коэффициент может принимать значение от −∞ до +∞. 

Коэффициент корреляции :

$$r_{x,y}=\frac{\sum(x_{i}-\bar{x})*(y_{i}-\bar{y}) }{\sqrt{\sum(x_{i}-\bar{x})^2}*\sqrt{\sum(y_{i}-\bar{y})^2} }=\frac{cov_{x,y}}{\sigma _{x}\sigma _{y}}$$

Коэффициент корреляции это не что иное, как стандартизованный коэффициент ковариации. Мы можем в числителе поставить коэффициент ковариации и отнести его, стандартизовать по произведению двух среднеквадратичных отклонений. Коэффициент корреляции варьирует от −1 до +1, что намного удобнее, чем неограниченное распределение коэффициента ковариации. Но суть этого коэффициента в принципе та же самая: отрицательные значения говорят об отрицательной связи, а положительные — о положительной. А значения, близкие к нулю,— об отсутствии связи. Этот коэффициент  был введен Карлом Пирсоном и параллельно с ним Браве. Поэтому этот коэффициент обычно называют коэффициентом Пирсона или коэффициентом Браве — Пирсона. Но эту корреляцию, которую еще называют корреляцией Пирсона, у нее есть одна особенность. Этот коэффициент учитывает только линейную составляющую ковариации величин. Например, при криволинейной зависимости между величинами, коэффициент корреляции Пирсона выдаст только линейную составляющую этой взаимосвязи. Соответственно, мы немножко занизим силу связи между такими величинами. Эту особенность надо принимать в расчет. Если при визуальном рассмотрении облака точек видны признаки нелинейной связи, то для оценки корреляции между величинами лучше подойдут другие коэффициенты.

Вот если зависимость нелинейна, но при этом монотонна, то есть она либо убывает, либо постоянно возрастает, то лучше для таких целей, для оценки связи между величинами в такой ситуации лучше подходят ранговые коэффициенты корреляции. Обычно используют два близких по смыслу коэффициента: это коэффициент Спирмена и коэффициент Кендалла. Поэтому, если возникнет необходимость в оценке корреляции между величинами, которые связаны друг с другом нелинейно, в этой ситуации мы должны перейти к ранговым коэффициентам корреляции.

Тестирование статистической значимости коэффициента корреляции

Коэффициент корреляции, это статистика, описывающая степень взаимосвязи двух сопряженных переменных. Следовательно здесь применима логика статистического критерия. Как и в случае, например, t-критерия,  для коэффициента корреляции есть нулевая гипотеза, а именно, что в генеральной совокупности связь отсутствует, то есть корреляция равна нулю. Важно, что альтернативная гипотеза должна быть сформулирована до начала тестирования, поскольку здесь могут быть как двусторонние так и односторонние альтернативные гипотезы. Тестирование этих двух типов гипотез зависит от того, какую альтернативу мы рассматриваем. Разными будут при этом пороги отвержения нулевой гипотезы, и при одном и том же уровне альфа мы будем отвергать нулевую гипотезу при одном и том же, при разных значениях. Соответственно, двусторонняя альтернативная гипотеза будет говорить о том, что связь существует, то есть корреляция не равно нулю, то есть связь есть, неважно какая, положительная - отрицательная, но связь есть. Одностороннее же гипотеза подразумевает, что мы утверждаем, что корреляция строго положительная или строго отрицательная, и именно это мы должны проверить нашим тестом. То есть мы должны заранее утверждать какой тип связи мы считаем присутствует в генеральной совокупности. Это важно потому что односторонние тесты более мягкие, и мы можем прийти к выводу что связь есть там где ее нет. То есть мы можем совершить статистическую ошибку и найти корреляцию там, где ее нет. 

Аналогично t-критерию, если нулевая гипотеза не может быть отвергнута, то никаких суждений о существовании взаимосвязи мы вынести не можем. Если мы отвергаем нулевую гипотезу при определенном значении альфа, то мы сохраняем альтернативную гипотезу. Собственно, вот такова логика тестирования этой гипотезы. Поскольку вычисленный коэффициент корреляции, это лишь выборочная оценка некоторого генерального коэффициента, то для этой оценки необходимо вычислить стандартную ошибкy.

Стандартная ошибка коэффициента корреляции Пирсона вычисляется по  формуле : 

$$SE_{r}=\sqrt{\frac{1-r^2}{n-2}}$$

Согласно $H_{0}:\rho =0$ , поэтому t-статистика выглядит так :

$$t=\frac{r-\rho }{SE_{r}}=\frac{r}{SE_{r}}$$

Как и в случае с t-критерием, нам понадобится стандартизовать. Стандартизированная величина t, которая будет равна отношению самого коэффициента корреляции к его стандартной ошибке, и далее, эта величина, она тоже подчиняется тому же t-распределению как и в случае с разностью средних значений, и единственный параметр t-распределение, который называется числом степеней свободы, в данном случае будет равен n-2. Далее, точно так же надо сравнить полученные значение t с пороговым значением, если модуль этой величины превышает пороговое значение, то можно утверждать, что корреляция статистически значима.

Корреляционный анализ в R

Для расчета и проверки значимости различных типов коэффициентов корреляции (Пирсона, Спирмена или Кендала) в R используется функция

cor.test(x,y,method)

x,y - данные в виде векторов одинаковой длины
method - тип коэффициента корреляции

Оценим корреляцию ранее мощности двигателя и расхода топлива двумя методами :

Коэффициент корреляции Пирсона

cor.test(x = df$hp, y = df$mpg, method = "pearson")

Pearson's product-moment correlation

data:  df$hp and df$mpg
t = -6.742, df = 30, p-value = 1.788e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.885 -0.586
sample estimates:
       cor 
-0.776 

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

Ранговый коэффициент корреляции Спирмена

cor.test(x = df$hp+rnorm(32)*0.0001, y = df$mpg+rnorm(32)*0.0001, method = "spearman")

Spearman's rank correlation rho

data:  df$hp + rnorm(32) * 1e-04 and df$mpg + rnorm(32) * 1e-04
S = 10344, p-value = 7.386e-08
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.896 

Как и ожидалось, из-за нелинейности связи, Пирсон показал заниженное значение.


Комментариев нет:

Отправить комментарий