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

Статистика в R : простая линейная регрессия

Регрессионный анализ является развитием идеи корреляционного анализа. Корреляционный анализ позволяет выявить силу и направление связи, а также оценить статистическую значимость этой связи. Регрессионный анализ дает возможность построить модель, описывающую эту связь. То есть с его помощью можно вывести формулу модели и построить график, визуализирующий данную модель. 


 Термин "регрессионный анализ" был введен Фрэнсисом Гальтоном, кузеном Чарльза Дарвина, который показал, что рост потомков очень высоких людей оказывается меньше, чем их родителей, то есть рост регрессирует к популяционной средней. Коэффициент, описывающий связь роста родителей и роста потомков, Гальтон назвал коэффициентом регрессии.

Регрессионный анализ позволяет построить стохастическую линейную модель, описывающую взаимосвязь между величинами в генеральной совокупности. Но  то, что происходит в генеральной совокупности, нам неизвестно. Мы можем лишь оценить то, что там происходит по выборке. Поэтому в регрессионном анализе  мы только предполагаем, что в генеральной совокупности существует некоторая взаимосвязь между величиной x и величиной y. Эта взаимосвязь определяется формулой с коэффициентами β0 и β1. И поскольку в природе все величины, которые существуют объективно, подвержены изменчивости, поэтому данная модель обязательно будет иметь случайную часть ε, вот она. Соответственно, когда мы пытаемся оценить эту связь, то мы должны взять некоторую выборку, и данная выборка сформирует нам два сопряженных вектора x и y, в каждом из которых будет пара чисел. То есть у некоторого объекта мы изучили объект i и изучили величину x и величину y. Соответственно, дальше полученные два ряда чисел ложатся в основу регрессионного анализа, который позволяет построить модель, которая будет являться оценкой той модели, которая описывает связь между xи y в генеральной совокупности. Такая модель называется регрессионной моделью. Так как речь идет об оценках, то должны появиться стандартные ошибки и оценки статистической значимости. 

Компоненты регрессионной модели :

$$y_{i}=b_{0}+b_{1}x_{i}+e_{i}$$

Переменные :
y- зависимая (отклик, response)
x - независимая (предиктор, predictor)

Параметры модели :
$b_{0}$ - свободный член (отрезок, intercept)
$b_{1}$ - угловой коэффициент (slope)

Простая регрессионная модель — это простая стохастическая линейная модель. Формула включает в себя, собственно, следующие величины: y, это зависимая переменная, ее еще называют переменной отклика, или response variable.  Это та величина, поведение которой нас интересует. x — это независимая переменная, предиктор, которая, как мы предполагаем, определяет поведение нашей зависимой переменной. Параметры модели b0 и b1имеют специальное название. Параметры b0 — это свободный член, или его еще называют отрезком, или intercept. И второй параметр b1 — его называют угловым коэффициентом, или slope. Случайная часть модели —ее называют остатками, она обозначается латинской буквой e, и она определяет варьирование, связанное с неучтенными фиксированной частью модели, то есть вот этой частью модели, какими-то факторами. Данные, вошедшие в выборку, на основе которой проведен данный анализ, можно изобразить в виде графика, который имеет вид облака точек. Каждая точка — это отдельная пара измерений.

Разберем составляющие модели отдельно :

$$y_{i}=b_{0}+b_{1}x_{i}+e_{i}$$

$$\widehat{y}_{i}=b_{0}+b_{1}x_{i}$$

предсказанные значения лежат на линии регрессии

$$e_{i}=y_{i}-\widehat{y}_{i}$$
остатки, отклонение наблюдаемых значений зависимой переменной от предсказанных моделью

$b_{1}$ - угловой коэффициент (slope) , показывает на сколько единиц изменится предсказанное значение $\widehat{y}$ при изменении предиктора x на одну единицу

$b_{0}$ - свободный член (intercept) , определяет предсказанное значение $\widehat{y}$ при значении предиктора x=0. Иногда не имеет "смысла", это просто поправочный коэффициент.

Коэффициенты модели рассчитываются по исходным данным методом наименьших квадратов.

Метод наименьших квадратов (МНК / OLS)

МНК - один из способов поиска коэффициентов регрессии. В основе метода лежит поиск линии регрессии, имеющую минимальную сумму квадратов остатков


$$\sum e_{i}^2\rightarrow min$$


Подбор коэффициентов линейной регресcии в R


В R функция для работы с линейным моделями называется lm, от английского linear model.

fit <- lm(formula, data)

data - исходные данные

Модель :
    - простая линейная регрессия
        - уравнение :

$$\widehat{y}_{i}=b_{0}+b_{1}x_{i}$$

        - formula :
                Y~X
                Y~1+X
                Y~X+1
    - простая регрессия, $b_{0}=0$
        -уравнение :

$$\widehat{y}_{i}=b_{1}x_{i}$$

        - formula :
                Y~-1+X
                Y~X-1
    - уменьшенная простая регрессия
        -уравнение :
 $$\widehat{y}_{i}=b_{0}$$
        - formula :
                Y~1
                Y~1-X
    - множественная линейная регрессия
        - уравнение :
 $$\widehat{y}_{i}=b_{0}+b_{1}x_{1i}+b_{2}x_{2i}$$
        - formula :
                Y~X1+X2

Модели линейной регрессии можно записать на специальном языке формул. Модель простой линейной регрессии с одним единственным предиктором, ее можно записать на языке формул тремя разными способами. Во всех формула единообразно, в середине стоит знак тильда, а слева указывается зависимый аргумент прямо по имени переменной и как она хранится в названии датафрейма из аргумента data. Справа в самом минимальном варианте надо указать название этого единственного предиктора, в примере это X. Коэффициент b нулевое никаким образом не обозначается, по умолчанию он будет добавлен в модель, но если вы хотите его упомянуть в явном виде, нужно записать единицу. Неважно, в каком месте вы ее запишите, вы можете записать перед именем предиктора X или после. Модель все равно будет подобрана одна и та же, модель простой линейной регрессии. Что будет если вы хотите подобрать такую зависимость, которая проходит через ноль? У такой линейный регрессии коэффициент нулевое будет равен нулю, то есть intercept-а в модели не будет таким образом. Очень просто можно записать формулу такой модели. Вам нужно перед знаком intercept-а написать минус. Это значит, что этот элемент формулы будет убран из модели, и опять же неважно в начале это записывать или в конце. Модель линейной регрессии, которая проходит через ноль без intercept-а. Бывают случаи, когда нам нужна самая простая модель простой линейной регрессии, в которой вообще нет никаких предикторов, в которых все описывается одним единственным коэффициентом, b нулевое, т.е. мы описываем среднее значение нашей переменной отклика, среднее значение Y. Такую модель тоже очень легко записать в R. Вы можете написать справа от тильды просто единицу, то есть в модели есть единственный коэффициент, это intercept, или вы можете в явном виде перечислить те предикторы, которые вы хотите из нее исключить. Эта запись будет синонимична, так обычно никто не пишет, конечно, пишут вот первым способом.  Если у вас несколько предикторов в модели, вы их просто перечисляете, соединяя знаком плюс. 

Подберем модель в R для наших данных и рассмотрим результаты

model1 <- lm(mpg ~ hp, data = df)

summary(model1)

Это можно сделать при помощи функции lm и summary. Функция summary, это функция generic, по английский называется функция обертка, которая смотрит к какому классу относится ее аргумент и пытается найти метод, который подходит для этого класса. Если мы передадим функции summary линейную модель, то тогда она посмотрит к какому классу относится это линейная модель, потому что линейных моделей тоже бывает много разных, найдет подходящий метод и передаст эту модель дальше этому методу, который сделает с ней что-то, покажет нам результаты в человеко-читаемом виде. Передадим нашу модель model1 в функцию summary и нам сразу покажут в консоли, как мы модель создавали.  Последний раздел, это таблицы с коэффициентами линейной регрессии. Это оценки коэффициентов intercept-а, и коэффициента угла наклона. Здесь записаны стандартные ошибки этих коэффициентов, а вся правая часть таблицы связана косвенным образом с этими стандартными ошибками, потому что это тест о значимости коэффициентов.

Стандартные ошибки коэффициентов

Параметр :

$$\beta _{0}$$

    - Оценка

 $$b_{0}=\bar{y}-b_{1}\bar{x}$$

    - Стандартная ошибка

$$SE_{b0}=\sqrt{MS_{e}\left [ \frac{1}{n}+\frac{\bar{x}}{\sum \left ( x_{i} - \bar{x}\right )^2} \right ]}$$   

$$\beta _{1}$$ 

    - Оценка

$$b_{1}=\frac{\sum \left [ \left ( x_{i}-\bar{x} \right ) \left ( y_{i}-\bar{y} \right )\right ]}{\sum \left ( x_{i}-\bar{x} \right )^2}$$

    - Стандартная ошибка

$$SE_{b1}=\sqrt{\frac{MS_{e}}{\sum \left ( x_{i}-\bar{x} \right )^2}}$$

 $$\varepsilon_{i}$$

    - Оценка

$$e_{i}=y_{i}-\hat{y}_{i}$$

    - Стандартная ошибка

$$\approx \sqrt{MS_{e}}$$

 Какова природа стандартных ошибок? Дело в том, что при помощи линейной регрессии мы пытаемся описать зависимость, которая есть в генеральной совокупности используя данные о выборке. В генеральной совокупности, в зависимости, которые там существует, помимо нас и наших знаний о ней, есть параметры. Эти параметры обозначаются греческими буквами бета, бета нулевое, бета первое и остатки обозначаются эпсилон. Но мы работаем с выборкой, и нам не доступны эти значения параметров, мы их пытаемся оценить по выборке и получаем оценки коэффициентов, b нулевое, b первое и к каждому наблюдению можно посчитать величину остатка, то есть насколько от предсказанного значения отличается реальное значение наблюдения. В связи с тем, что мы работаем с выборкой, и выборка несовершенна, одна выборка может отличаться от другой и никогда выборка полностью не идентична генеральной совокупности, у нас возникает неопределенность, и эту неопределенность мы можем оценить, используя стандартные ошибки. По своей формуле можно посчитать стандартную ошибку для каждого из коэффициентов, а в качестве оценки стандартной ошибки остатков можно использовать корень из остаточной суммы квадратов. В нашем случае мы имеем дело с выборкой из 32 автомобилей, но они все не идентичны, они все отличаются друг от друга и они явно отличаются от той генеральной совокупности, на которую бы мы хотели экстраполировать наши результаты. Поэтому возникает некоторая неопределенность, наша неуверенность в полученных коэффициентах регрессии, которая выражается в определении доверительных интервалов.

Доверительные интервалы коэффициентов и доверительная зона регрессии

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

 Если таких выборок сделать достаточно много, мы можем построить распределение для каждого из этих параметров, распределение для интерсепт и распределение для коэффициента угла наклона. И эти распределения выглядят достаточно колоколообразно. Если объем выборки у нас достаточно велик, если выполняются условия, при которых справедлива центральная предельная теорема, то распределение значений параметров будет приближаться к нормальному, поэтому мы можем использовать характеристики этого распределения для того, чтобы построить доверительный интервал — зону, куда попадает 95 % всех значений параметров в повторных выборках. 

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

 $$b_{k}\pm t_{\alpha ,df}*SE_{bk}$$

$t_{\alpha ,df}$ - квантиль t-распределения

df=n-p - число степеней свободы

p - число параметров (для простой линейной регрессии p=2)

k=0,...,p - конкретный коэффициент линейной регрессии

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

Доверительные интервалы к параметрам регрессии в R

Доверительные интервалы к параметрам регрессии в R можно получить при помощи функции confint У этой функции один аргумент, мы должны туда передать название переменной, в которую мы сохранили модель. Функция confint выдаст нам нижнюю и верхнюю границы доверительных интервалов. Таким образом, если нас интересует зависимость величины mpg от мощности двигателя и коэффициент угла наклона, соответственно, мы получим информацию, что в 95 %повторных выборок из генеральной совокупности значение этого коэффициента будет лежать вот в таких пределах. 

Доверительные интервалы коэффициентов и доверительная зона регрессии  

Доверительные интервалы к параметрам регрессии в R

coef(model1)

(Intercept)          hp 
    30.0989     -0.0682 

confint(model1)

              2.5 %  97.5 %
(Intercept) 26.7619 33.4358
hp          -0.0889 -0.0476

Доверительная зона регрессии в R

Построим доверительную зону регрессии при помощи пакета ggplot. За это построение отвечает geom_smooth. Smooth — сглаживать. Это geom, который строит разного вида зависимости. В данном случае мы используем метод «Линейная модель». geom_smooth сразу построит нам и линию регрессии, и ее доверительную зону.

ggplot(df, aes(x = hp, y = mpg)) +
  geom_point() +
  labs(x="Мощность двигателя, л.с.", y = "Потребление топлива, мили/галлон") +
  geom_smooth(method = "lm")




По умолчанию geom_smooth строит доверительную зону для 95 %-го интервала. Но на самом деле мы можем построить доверительную зону для любых других интервалов, если изменим аргумент level,который задает уровень значимости α. Соответственно, если вот в эту узкую зону линия регрессии попадает в 95 % случаев при повторных выборках, то вот в эту более широкую зону она попадет уже в 99,9 % повторных выборок, то есть почти всегда. Естественно, это зона будет гораздо шире. 

Использование регрессии для предсказаний

Можно использовать уравнение регрессии для предсказаний. Для этого нам нужно просто подставить конкретное значение в уравнение. Например, мы можем захотеть узнать, какой сколь миль можно проехать на одном галлоне топлива с двигателем мощностью 250 л.с. Но правда ли, что у всех двигателей с такой мощностью будет такой же пробег. Нет, вряд ли, потому что двигатели отличаются друг от друга. Соответственно, у нас есть некоторая неопределенность, когда мы оцениваем по выборке значения параметров линейной регрессии, и значит должна быть какая-то неопределенность в предсказаниях, когда мы используем уравнение регрессии для предсказания. 

Эту неопределенность можно оценить двумя способами, можно получить групповые предсказания или предсказания на индивидуальном уровне. Групповые предсказания - это значит, что мы показываем предсказывает положение среднего значения в выборки. Если мы говорим об индивидуальных предсказаниях, мы предсказываем зону, в которую попадут 95% всех наблюдений в повторной выборке, то есть уже не среднего значения, а конкретных наблюдений. Естественно, что зона групповых предсказаний, она уже, чем зона предсказаний индивидуальных. 

Функция predict()

predict(object, newdata, interval, level)

object - модель

newdata - новые данные для предсказания

interval - тип предсказаний

level = 1 - alf

В R можно получить предсказание и групповые и индивидуальные при помощи функции predict. Функция predict - это функция generic, на самом деле она скрывает под одной оберткой много разных функций. Функция predict смотрит к какому классу относится ее аргумент и выбирает подходящий метод работы с ним. Например, если мы ей дадим объект, который создан при помощи функции lm, она на самом деле вызовет функцию predict.lm со своими параметрами. Это функция уже умеет работать с линейными моделями и она нам даст предсказание, поэтому, если вы хотите посмотреть какие аргументы мы можем использовать в методе predict, вам нужно смотреть help именно к predict lm. Какие аргументы нам потребуется? Ну во-первых, конечно же первый аргумент, это модель, которую мы создали при помощи функции lm. Мы должны этой модели будем передать новые данные в формате датафрейма, хорошо бы ей указать, какой тип предсказаний мы хотим получить, групповые или индивидуальные. По умолчанию это будут групповые предсказания. И последний аргумент, который нам сейчас важен, это level, это вероятность, которая описывает ширину получившегося интервала. 

 Подготовим данные, которые мы будем использовать для предсказаний функции predict. Нам нужно будет вручную создать датафрейм при помощи конструктора data.frame, и в этом датафрейме будет единственная переменная. Дело в том, что если мы хотим использовать регрессию для предсказаний, в датафрейме обязательно должны содержаться все переменные, которые участвовали в построении модели. Здесь у нас одна переменная, и соответственно, в датафрейме будет тоже одна. Если мы хотим получить предсказание для одного наблюдения, в датафрейме будет всего одна строка. В данном случае там будет записано 250

Данные для предсказаний

new_data <- data.frame(hp = 250)
new_data

  hp
1 250

Если мы хотим получить групповые предсказания, то функцией predict можно сказать, а можно и не говорить, что interval = 'confidance'. По умолчанию параметр interval как раз принимает это значение. И по умолчанию уровень значимости, это 95%. Результаты работы функции predict, это предсказанное значение и нижняя и верхняя граница доверительного интервала. Таким образом мы узнаем, что машин с двигателем мощностью 250 л.с.  mpg в 95% повторных выборок будет заключаться вот в этом интервале.

Групповые предсказания

predict(model1, newdata = new_data, interval = "confidence", level = 0.95)

  fit  lwr  upr

1  13 10.5 15.6

Интервал от 10.5 до 15.6 с вероятностью 95% будет включать в себя среднее значение расстояния у машины с двигателем мощностью 250 л.с.

Индивидуальные предсказания

Какое значение mpg можно ожидать у машины с двигателем в 250 л.с.

predict(model1, newdata = new_data, interval = "prediction", level = 0.95)

  fit  lwr  upr
1  13 4.75 21.3

Чтобы получить индивидуальное предсказание при помощи функции predict, нам нужно изменить значение только одного аргумента, interval, мы задаем аргумент interval равный prediction. Предсказательный интервал, он шире, и значение mpg у всех машин с двигателем в 250 л.с.  будут с 95% вероятностью содержаться в этом интервале, от 4.75 до 21.3. 

Чтобы построить график, на котором мы одновременно изобразим линию зависимости и доверительные зоны предсказаний, доверительные зоны регрессий и точки исходных наблюдений, нам понадобится извлечь предсказание из результатов функции predict. Мы сохраним эти результаты в переменную model1, а чтобы нам в том же датафрейме были доступны исходные точки наблюдений, давайте мы их просто добавим при помощи функции data.frame. Теперь у нас есть все данные для того, чтобы построить график. 

Данные для графика доверительных зон регрессии и ее предсказаний

df_predicted <- predict(model1, interval = "prediction") # Подготовка данных
head(df_predicted, 3)

                     fit  lwr  upr
Cadillac Fleetwood  16.1 8.01 24.2
Lincoln Continental 15.4 7.29 23.6
Camaro Z28          13.4 5.12 21.6

df_predicted <- data.frame(df, df_predicted)          # Добавим исходные данные
head(df_predicted, 3)

                     mpg cyl disp  hp drat   wt qsec vs        am gear carb            car.name mpg_z
Cadillac Fleetwood  10.4   8  472 205 2.93 5.25 18.0  V automatic    3    4  Cadillac Fleetwood -1.61
Lincoln Continental 10.4   8  460 215 3.00 5.42 17.8  V automatic    3    4 Lincoln Continental -1.61
Camaro Z28          13.3   8  350 245 3.73 3.84 15.4  V automatic    3    4          Camaro Z28 -1.13
                    mpg_type  fit  lwr  upr
Cadillac Fleetwood     below 16.1 8.01 24.2
Lincoln Continental    below 15.4 7.29 23.6
Camaro Z28             below 13.4 5.12 21.6


pl_limits <- ggplot(df_predicted, aes(x = hp, y = mpg)) +
  # Доверительная зона предсказаний
  geom_ribbon(aes(ymin = lwr, ymax = upr, fill = "Для предсказаний"), alpha = 0.5) +
  # Линия регрессии и ее доверительная зона
  geom_smooth(method = "lm", aes(fill = "Для регрессии")) +
  # Вручную настраиваем цвета заливки
  scale_fill_manual(name = "Доверительные \nзоны", values = c("powderblue", "gray")) +
  # Исходные наблюдения
  geom_point() +
  # Подписи
  labs(x = "Мощьность двигателя, л.с.", y = "Потребление топлива, мили/галлон",
       title = "Доверительные зоны регрессии и ее предсказаний")
pl_limits








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

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