Язык программирования R для анализа данных: лекция 9

Ограничения ANOVA, линейная регрессия

Elena U

План лекции

  • Ограничения ANOVA

  • Ковариация и корреляция

  • Простая линейная регрессия: формула, реализация в R, интерпретация

  • Ограничения линейной регрессии (в т.ч. нормальность распределения остатков)

  • Множественная линейная регрессия

  • Ограничения ANOVA

Ковариация и корреляция

Данные для работы

Датасет с каггла о связи ВВП (GDP) на душу населения и счастья.

library(tidyverse)
df <- read_csv('https://raw.githubusercontent.com/ubogoeva/Rcourse_NSU/master/posts/lectures/data/happyscore_income.csv')
head(df)
# A tibble: 6 × 11
  country...1 adjusted_satisfaction avg_satisfaction std_satisfaction avg_income
  <chr>                       <dbl>            <dbl>            <dbl>      <dbl>
1 Armenia                        37              4.9             2.42      2097.
2 Angola                         26              4.3             3.19      1449.
3 Argentina                      60              7.1             1.91      7101.
4 Austria                        59              7.2             2.11     19457.
5 Australia                      65              7.6             1.8      19917 
6 Azerbaijan                     46              5.8             2.27      3382.
# ℹ 6 more variables: median_income <dbl>, income_inequality <dbl>,
#   region <chr>, happyScore <dbl>, GDP <dbl>, country...11 <chr>
df <- df %>% 
  rename(country = country...1) %>% 
  select(!country...11)

Посмотрим на данные

Используя базовую generic-функцию summary().

summary(df)
   country          adjusted_satisfaction avg_satisfaction std_satisfaction
 Length:111         Min.   :19.00         Min.   :2.500    Min.   :1.380   
 Class :character   1st Qu.:40.00         1st Qu.:5.100    1st Qu.:1.910   
 Mode  :character   Median :48.00         Median :6.000    Median :2.130   
                    Mean   :48.73         Mean   :5.937    Mean   :2.125   
                    3rd Qu.:57.00         3rd Qu.:7.000    3rd Qu.:2.330   
                    Max.   :74.00         Max.   :8.500    Max.   :3.190   
   avg_income      median_income     income_inequality    region         
 Min.   :  572.9   Min.   :  415.5   Min.   :24.21     Length:111        
 1st Qu.: 1519.4   1st Qu.: 1167.7   1st Qu.:32.18     Class :character  
 Median : 3889.3   Median : 2647.0   Median :36.48     Mode  :character  
 Mean   : 6442.8   Mean   : 5186.0   Mean   :38.42                       
 3rd Qu.: 7867.4   3rd Qu.: 6581.1   3rd Qu.:43.38                       
 Max.   :26182.3   Max.   :22240.2   Max.   :63.73                       
   happyScore         GDP        
 Min.   :2.839   Min.   :0.0153  
 1st Qu.:4.568   1st Qu.:0.5099  
 Median :5.268   Median :0.9186  
 Mean   :5.422   Mean   :0.8419  
 3rd Qu.:6.392   3rd Qu.:1.1495  
 Max.   :7.587   Max.   :1.5639

Больше всего нас будут интересовать переменные happyScore и GDP.

Посмотрим на данные

С помощью функции skim() из пакета skimr.

skimr::skim(df)

(Аутпут очень большой, поэтому прикрепила скрином).

Ковариация

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

Формула ковариации:

\[ \sigma_{xy} = cov(x, y) = \frac{\sum_{i = 1}^n(x_i - \overline{x})(y_i - \overline{y})}{n} \]

Оценка ковариации по выборке:

\[ \hat{\sigma}_{xy} = \frac{\sum_{i = 1}^n(x_i - \overline{x})(y_i - \overline{y})}{n-1} \]

Если оценить ковариацию переменной с самой собой, какая будет формула?

Ковариация

Ковариация переменной с самой собой - это дисперсия.

Давайте посчитаем ковариацию: функция cov().

df %>% 
  select(happyScore, GDP) %>% 
  cov()
           happyScore       GDP
happyScore  1.3942899 0.3615849
GDP         0.3615849 0.1502264

По главной диагонали ковариация переменной с самой собой - дисперсия.

df %>% 
  select(happyScore, GDP) %>% 
  var() # тот же самый результат с помощью функции var
           happyScore       GDP
happyScore  1.3942899 0.3615849
GDP         0.3615849 0.1502264

Можем ли мы сделать вывод, это относительно большая ковариация между переменными или нет?

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

Проблема ковариации в том, что она привязана к исходной шкале и измеряется в пределах \([-\infty; \infty]\). Неплохо было бы иметь возможность нормировать. Для этого используется коэффициент корреляции.

Формула коэффициента корреляции Пирсона:

\[ \rho_{xy} = \frac{\sigma_{xy}}{\sigma_x \sigma_y} = \frac{\sum_{i = 1}^n(x_i - \overline{x})(y_i - \overline{y})}{\sqrt{\sum_{i = 1}^n(x_i - \overline{x})^2}\sqrt{\sum_{i = 1}^n(y_i - \overline{y})^2}} = \frac{1}{n}\sum_{i = 1}^n z_{x,i} z_{y, i} \]

Оценка коэффициента корреляции Пирсона по выборке:

\[ r_{xy} = \frac{\hat{\sigma}_{xy}}{\hat{\sigma}_x \hat{\sigma}_y} = \frac{\sum_{i = 1}^n(x_i - \overline{x})(y_i - \overline{y})}{\sqrt{\sum_{i = 1}^n(x_i - \overline{x})^2}\sqrt{\sum_{i = 1}^n(y_i - \overline{y})^2}} = \frac{1}{n - 1}\sum_{i = 1}^n z_{x,i} z_{y, i} \]

По сути это ковариация, деленная на стандартное отклонение обеих переменных.

Коэффициент корреляции измеряется в пределах \([-1; 1]\).

Вычисляем корреляцию

Используем функцию cor() базового R.

df %>% 
  select(happyScore, GDP) %>% 
  cor() # много это или мало?
           happyScore       GDP
happyScore  1.0000000 0.7900609
GDP         0.7900609 1.0000000

Посмотрим визуально на распределение данных.

plot(df$GDP, df$happyScore)

Значим ли этот коэффициент корреляции?

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

Для оценки статистической значимости коэффициента корреляции используется функция cor.test().

Нулевая гипотеза - коэффициент корреляции равен нулю, альтернативная - не равен.

cor.test(df$happyScore, df$GDP)

    Pearson's product-moment correlation

data:  df$happyScore and df$GDP
t = 13.455, df = 109, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.7079171 0.8511169
sample estimates:
      cor 
0.7900609

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

Непараметрические аналоги коэффициента корреляции Пирсона

Часто используется коэффициент корреляции Спирмена (Spearman), если в данных есть выбросы. Математика метода точно такая же как у коэффициента Пирсона, только вместо оригинальных значений используются ранги.

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

cor.test(df$happyScore, df$GDP, method = 'spearman')

    Spearman's rank correlation rho

data:  df$happyScore and df$GDP
S = 47026, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.7936732

Наличие значимого коэффициента корреляции ничего не говорит о причинно-следственной связи между переменными!

Непараметрические аналоги коэффициента корреляции Пирсона

Также есть коэффициент корреляции Кендалла.

cor.test(df$happyScore, df$GDP, method = 'kendall')

    Kendall's rank correlation tau

data:  df$happyScore and df$GDP
z = 9.3551, p-value < 2.2e-16
alternative hypothesis: true tau is not equal to 0
sample estimates:
      tau 
0.6013104

Сравнение различных коэффициентов корреляции

Пирсона Спирмена Кендалла
Выявляет линейную зависимость Выявляет любую монотонную зависимость Выявляет любую монотонную зависимость
Количественная или интервальная шкала шкала >= ранговая (то есть ранговая, интервальная, количественная) шкала >= ранговая
Неустойчивость к выбросам Устойчивость к выбросам Устойчивость к выбросам

Пример ситуации, когда коэффициент Спирмена более применим:

Проведя корреляционный анализ, мы лишь ответили на вопрос “Существует ли статистически значимая связь между величинами?”

Сможем ли мы, используя это знание, предсказать значения одной величины, исходя из знаний другой?

Линейная регрессия

Какие бывают регрессионные модели?

  • Линейные и нелинейные

    • Линейные

      \[ y = b_0 + b_1x \] \[ y = b_0 + b_1x_1 + b_2x_2 \]

      \[ y = b_0 + b_1x_1^2 + b_2x_2^3 \]

    • Нелинейные

      \[ y = b_0 + b_1^x \] \[ y = b_0^{b_1x_1+b_2x_2} \]

  • Простые и множественные

Линейная регрессия: формула

Линейная регрессия позволяет предсказать значение одной переменной на основании значений другой. Обе переменные должны быть количественными.

Формула для простой линейной регрессии - по сути это формула прямой из алгебры + остатки:

\[ Y = ax + b + \epsilon \]

\(Y\) - зависимая переменная, что пытаемся предсказать, \(\epsilon\) - это ошибки или остатки модели (residuals) - то есть то, что наша модель не смогла предсказать.

Для более сложных моделей нотация немного меняется:

\[ Y = \beta_0 + \beta_1x_1 + \beta_2x_2 + ... + \beta_nx_n + \epsilon \]

Линейная регрессия

В R функция для линейной регрессии - lm(). Давайте попробуем предсказать значение happyScore на основании ВВП на душу населения (GDP).

model <- lm(happyScore ~ GDP, data = df)
model

Call:
lm(formula = happyScore ~ GDP, data = df)

Coefficients:
(Intercept)          GDP  
      3.395        2.407

Мы получили коэффициенты: Intercept, GDP. Что это значит?

Как эти коэффициенты подбираются?

  • Метод наименьших квадратов (для простых моделей) -> нарисовать

  • Метод максимального правдоподобия (для более сложных)

График регрессионной прямой на основании полученных коэффициентов

happyScore= Intercept + GDP * Slope = 3.395 + GDP*2.407

Попробуем нарисовать регрессионную прямую на нашем графике.

ggplot(df, aes(GDP, happyScore))+
  geom_point(alpha = 0.8)+
  geom_abline(slope = model$coefficients[2], 
              intercept = model$coefficients[1], color = 'blue')

Немного приведем код в порядок и разобьем еще по региону, чтобы было интереснее анализировать.

График зависимости happyScore от ВВП

ggplot(df, aes(GDP, happyScore))+
  geom_point(aes(color = region), alpha = 0.8)+
  geom_abline(slope = model$coefficients[2], 
              intercept = model$coefficients[1], color = 'blue')+
  scale_x_continuous(limits = c(0, 1.5), breaks = seq(0, 1.5, 0.25),
                     expand = c(0,0))+
  theme(axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))+
  theme_minimal()

С помощью geom_abline() задаем регрессионную прямую: подаем значения ax+b. Еще можно это сделать с помощью geom_smooth(), но тут задача показать, как прямую задать самостоятельно.

График зависимости happyScore от ВВП

Предсказание значений на основе регрессионной прямой

Функция predict() принимает на вход независимую переменную и предсказывает значение зависимой на основании нашей модели (подставляет в формулу).

Например, предскажем happyScore для GDP = 0.8. (Нужно подать именно датафрейм).

predict(model, data.frame(GDP = 0.8))
       1 
5.321038

Функция predict() ничего не знает о физическом смысле наших данных, поэтому мы можем предсказывать любые, даже неадекватные значения, например

predict(model, data.frame(GDP = -5))
        1 
-8.639174

Предсказания регрессионной модели имеют смысл только “внутри” графика!

Это называется интерполяция, а экстраполяция почти всегда будет неверной.

Экстраполяция и интерполяция

Более внимательно посмотрим на вывод регрессии

summary(model)

Call:
lm(formula = happyScore ~ GDP, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.05142 -0.47761 -0.02728  0.60008  1.53001 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.3955     0.1657   20.50   <2e-16 ***
GDP           2.4069     0.1789   13.46   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7272 on 109 degrees of freedom
Multiple R-squared:  0.6242,    Adjusted R-squared:  0.6207 
F-statistic:   181 on 1 and 109 DF,  p-value: < 2.2e-16

Что все это значит?

Самое важное из вывода регрессии

summary(model)$coefficients
            Estimate Std. Error  t value     Pr(>|t|)
(Intercept) 3.395491  0.1656664 20.49596 3.459016e-39
GDP         2.406933  0.1788837 13.45530 6.570905e-25

p-value есть для каждого коэффициента,

\(H_0: Intercept = 0, Slope = 0\).

\(H_1: Intercept \neq 0, Slope \neq 0\).

Какое значение p-value для нас важнее?

Коэффициент наклона.

Как рассчитывается значимость коэффициента наклона?

Рассчитывается относительно “нулевой” модели - то есть модели, в которую входит только интерсепт.

Коэффициент детерминации - мера качества модели

\(R^2\) - коэффициент детерминации. Описывает какую долю дисперсии зависимой переменной объясняет модель.

\[ R^2 = \frac{\color{blue}{SS_{r}}}{SS_{t}} \]

  • Измеряется от \([0, 1]\).

  • \(R^2 = r^2\) - для простой линейной регрессии коэффициент детерминации - квадрат коэффициента корреляции Пирсона.

Допущения линейной регрессии

  • Линейная связь

  • Независимость

  • Гомогенность дисперсий

  • Отсутствие коллинеарности предикторов (для множественной регрессии)

  • Нормальное распределение ошибок -> проиллюстрировать другие графики в R

    plot(model)

Линейность связи

Нелинейность связи видно на графиках остатков.

Проверка на линейность связи:

  • График зависимости y от x (и от других переменных, не включенных в модель).

  • График остатков от предсказанных значений.

Что делать с нелинейностью данных?

  • Добавить неучтенные переменные или взаимодействия

  • Применить линеаризующее преобразование (Осторожно!)

  • Применить обобщенную линейную модель с другой функцией связи (GLM)

График остатков для наших данных

ggplot(data.frame(fit = fitted(model), 
                             res = residuals(model)), 
                  aes(x = fit, y = res)) + 
  geom_point() + 
  geom_hline(yintercept = 0, linetype = 2) + 
  xlab("Fitted") + 
  ylab("Residuals")+
  theme_bw()

График остатков для наших данных

Независимость данных

Каждое значение \(y_i\) должно быть независимо от любого другого \(y_j\).

Это возможно проконтролировать на этапе сбора данных

Наиболее частые источники зависимостей:

  • Псевдоповторности (когда измерили один и тот же объект несколько раз)

  • Временные и пространственные автокорреляции

  • Неучтенные переменные

Нормальное распределение ошибок

Можно проверить с помощью qqPlot().

car::qqPlot(model$residuals, id = FALSE)

Проверка на гомогенность дисперсий (гомоскедастичность)

Многие тесты чувствительны к гетероскедастичности.

Лучший способ проверки на гомогенность дисперсий — график остатков от предсказанных значений.

Проверка на гомогенность дисперсий

взято из презентации Марины Варфоломеевой

Множественная линейная регрессия

Можно использовать несколько предикторов для объяснения зависимой переменной.

Больше предикторов - лучше модель?

model_income <- lm(happyScore ~ income_inequality + GDP, data = df)
summary(model_income)

Call:
lm(formula = happyScore ~ income_inequality + GDP, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.09360 -0.45110  0.01026  0.54468  1.43782 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       3.038567   0.416302   7.299 5.17e-11 ***
income_inequality 0.008124   0.008692   0.935    0.352    
GDP               2.460162   0.187829  13.098  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7276 on 108 degrees of freedom
Multiple R-squared:  0.6272,    Adjusted R-squared:  0.6203 
F-statistic: 90.85 on 2 and 108 DF,  p-value: < 2.2e-16

Множественная линейная регрессия

Переменная income_inequality оказалась незначимым предиктором. Возможно, есть смысл включить все возможные переменные в модель?

Проблема мультиколлинеарности во множественной линейной регрессии

Для множественной линейной регрессии очень плохо, когда независимые переменные коррелируют друг с другом. Чтобы оценить, насколько у нас выражена автокорреляция, построим плот корреляций (пакет corrplot, функция corrplot()).

corrplot::corrplot(cor(df %>% select(where(is.numeric))))

Большинство переменных так или иначе коррелируют с GDP, следовательно, включать их в модель нет особого смысла.

Линейная регрессия - частный случай общей линейной модели

И вообще большинство статистических тестов – это частный случай общей или обобщенной линейной модели.

Но про это подробнее в другой раз.

Вернемся к ANOVA

Ограничения ANOVA

  • Нормальность распределения - под вопросом (тоже самое, что и про t-test).

  • Равенство дисперсий - необязательно, если дизайн сбалансированный.

  • Нормальность распределения остатков - очень важно.

  • Независимость наблюдений.

Что такое остатки в ANOVA?

Дисперсионный анализ можно воспринимать как частный случай линейной модели, когда зависимая переменная количественная, а независимая номинативная.

Объяснение на рисунке.

Проверка на нормальность распределения остатков

Вспомним, что были за данные и проведем однофакторную анову:

diet <- readr::read_csv("https://raw.githubusercontent.com/Pozdniakov/tidy_stats/master/data/stcp-Rdataset-Diet.csv")
diet <- diet %>%
  mutate(weight_loss = weight6weeks - pre.weight,
         Dietf = factor(Diet, labels = LETTERS[1:3]),
         Person = factor(Person)) %>%
  drop_na()
fit_diet <- aov(weight_loss ~ Dietf, data = diet)

Проверка на нормальность распределения остатков

hist(residuals(fit_diet))

Распределение похоже на нормальное.

car::qqPlot(residuals(fit_diet))
[1] 48 15

Независимость наблюдений

Актуальны все те же правила, что были для линейной регрессии и т-теста.

Спасибо за внимание!

Если понравилось, переходите по ссылке: