R_HW_2.rnw

\documentclass[a4paper, 11pt]{article}
\usepackage[T2A]{fontenc}			% кодировка
\usepackage[utf8]{inputenc}			% кодировка исходного текста
\usepackage[english,russian]{babel}	% локализация и переносы
\usepackage{indentfirst}
\frenchspacing
\usepackage{geometry}
\geometry{left=2cm, right=1.5cm, top=2cm, bottom=1.5cm}
\usepackage{soul}
%%% Дополнительная работа с математикой
\usepackage{amsmath,amsfonts,amssymb,amsthm,mathtools,bm} %
\usepackage{icomma} % "Умная" запятая: $0,2$ --- число, $0, 2$ --- перечисление
\usepackage{etoolbox}
\usepackage{enumitem}
\usepackage{array,tabularx,tabulary,booktabs} % Дополнительная работа с таблицами
\usepackage{longtable}  % Длинные таблицы
\usepackage{multirow} % Слияние строк в таблице

\usepackage{lastpage}

\usepackage{url}
\begin{document}

% learn more about knitr: https://yihui.name/knitr/

<<setup, include=FALSE, cache=FALSE>>=
library(knitr)
opts_chunk$set(cache=TRUE, autodep=TRUE)
options(formatR.arrow=TRUE, width=90)
@

\title{Линейная регрессия в R}

\author{Варвара Кожухова}

\maketitle

\section*{Задание}
\begin{enumerate}
\item Загрузить набор данных для своего варианта, ознакомиться с его содержимым.
\item Построить график корреляционного поля для каждого фактора.
\item Построить уравнение парной линейной регрессии для каждого фактора.
\item Проверить значимость каждого из полученных уравнений регрессии. Показать уравнения регрессии с заданным в варианте доверительным интервалом на графиках.
\item Построить прогнозы по каждому из уравнений парной регрессии для заданных в варианте значений факторов.
\item Построить уравнение множественной линейной регрессии и получить корреляционную матрицу.
\item Построить прогноз по уравнению множественной регрессии для заданных в варианте значений факторов.

\end{enumerate}

\section*{Варианты}
\noindent \begin{tabularx}{\textwidth}{|X|X|X|X|X|X|} \hline
№ варианта & 1 & 2 & 3 & 4 & 5 \\ \hline
Набор данных & airquality & state.x77 & Cars93\footnote{необходимо подключить библиотеку MASS} & Cars93 & Cars93 \\ \hline
y & Ozone & Life Exp\footnote{поскольку в названии переменной нельзя использовать символ пробела, данную переменную необходимо переименовать командой colnames, например, в $Life\_Exp$, либо обращаться к ней по номеру столбца [,4].} & Price & Min.Price & Max.Price \\ \hline
Факторы и их значения для прогноза & Solar.R = 350, Wind = 8,3, Temp = 80 & Illiteracy = 1,0, Murder = 12, Income = 4000 & Horsepower = 200, RPM = 5200, Passengers = 4 & Horsepower = 210, RPM = 5500, Passengers = 4 & Horsepower = 220, RPM = 6000, Passengers = 6 \\ \hline
№ варианта & 6 & 7 & 8 & 9 & 10 \\ \hline
Набор данных & stackloss & longley & longley & LifeCycleSavings & Anscombe\footnote{необходимо подключить библиотеку car} \\ \hline
y & stack.loss & GNP & Employed & sr & education \\ \hline
Факторы и их значения для прогноза & Air.Flow = 55, Water.Temp = 20, Acid.Conc = 89 & Unemployed = 221, Armed.Forces = 180, Population = 125 & GNP.deflator = 102, Armed.Forces = 170, Population = 110 & pop15 = 35,5, pop75 = 1,5, dpi = 2500, ddpi = 2,15 & income = 3200, young = 347,8, urban = 425 \\ \hline
\end{tabularx}

\section*{Парная линейная регрессия}

Рассмотрим построение парной линейной регрессии на встроенном наборе данных cars. Будем рассматривать зависимость длины тормозного пути (переменная dist) от скорости (переменная speed). Построим график зависимости длины тормозного пути от скорости автомобиля.

<<libs, message=FALSE, echo=FALSE>>=
library("psych")# описательные статистики
library("lmtest") # тестирование гипотез в линейных моделях
library("ggplot2")# графики
library("dplyr") # манипуляции с данными
library("MASS") # подгонка распределений
@

<<df1, fig.width=4, fig.height=4>>=
d <- cars
qplot(data=d, speed, dist) +
  geom_point(aes(x=d$speed, y=d$dist), size = 1) + theme_bw(base_size = 12) +
  xlab("speed, miles/h") + ylab("distance, feet") +
  labs(title = "scatter plot")
@

Оценим модель линейной регрессии длины тормозного пути на скорость автомобиля. 
Для этого командой lm поместим в переменную model модель линейной регрессии, указав dist в качестве зависимой переменной, и через значок ~ переменную speed в качестве регрессора.

Тип lm представляет собой список из 12 элементов. Посмотрим на коэффициенты уравнения линейной регрессии.

<<df2>>=
model <- lm(data=d, dist~speed) # базовый пакет stats
model$coefficients
@

(Intercept) – это константа в уравнении регрессии, speed – коэффициент регрессии. 
Таким образом, уравнение регрессии имеет вид: 
\[dist_i^m = -17.579 + 3.9324 \cdot speed_i\]

Так же можно посмотреть значения вектора ошибок модели – разницу между реальной длиной тормозного пути dist и полученной по модели $dist_i^m$. Выведем первые 10 значений этого вектора с точностью две цифры после запятой. Более полный набор расчетов по модели можно получить командой summary.

<<df3>>=
options(digits = 3)
model$residuals[1:10]
summary(model) # базовый пакет base
@

Помимо коэффициентов регрессии, R выводит:
\begin{itemize}
\item стандартные ошибки коэффициентов (Std. Error); 
\item наблюдаемые значения t-критерия при проверке значимости коэффициентов регрессии (t value); 
\item P-значения для коэффициентов регрессии (P-value). 
\end{itemize}
Звездочками или точками в столбце справа R показывает значимость или незначимость коэффициентов: *** – значимы на уровне значимости менее 0.001; ** – значимы на уровне значимости 0.001; * – значимы на уровне значимости 0.01; . – значимы на уровне значимости 0.05 и т.д. Эти обозначения приведены в разделе Signif. codes.
Коэффициент детерминации (Multiple R-squared) равен 0.6511; скорректированный коэффициент детерминации (Adjusted R-squared) равен 0.6438. Наблюдаемое значения F-критерия проверки значимости уравнения в целом и P-значение:

F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

Таким образом, уравнение регрессии получилось значимым.
Проведем на графике полученную линию регрессии с 95\% доверительными интервалами.

<<df4, fig.width=4, fig.height=4>>=
qplot(data = d, speed, dist) + stat_smooth(method="lm", level = 0.95) + 
  theme_bw(base_size = 12)
@

Рассчитаем также 95\% доверительные интервалы для параметров линейной регрессии. Полученные по модели значения можно вывести командой fitted.
<<df5>>=
confint(model, level = 0.95) # базовый пакет stats
options(digits=4)
fitted(model) # базовый пакет stats
@

Рассчитать необъясненную сумму квадратов отклонений и полную сумму квадратов можно, воспользовавшись функцией deviance и уже известными функциями sum и mean. Для того, чтобы построить прогноз по полученной модели, нужно задать значения регрессора и поместить их в новый data.frame.

<<df6>>=
RSS <- deviance(model) # базовый пакет stats
TSS <- sum((d$dist-mean(d$dist))^2)
RSS; TSS
# создаем новый набор данных
nd <- data.frame(speed=c(40,60))
# Строим прогноз функцией predict
predict(model,nd)
@

\section*{Множественная линейная регрессия}

Рассмотрим встроенный набор данных по социально-экономическим показателям в 47 провинциях Швейцарии в 1888 г. Этот набор данных содержит 6 переменных по 47 наблюдений, каждая из которых измеряется в процентах (help(swiss)):
Fertility – рождаемость;
Agriculture – % мужчин, занятых в сельском хозяйстве;
Examination – % призывников, получивших высшую оценку на экзамене в армии;
Education – % призывников, имеющих образование помимо начального;
Catholic – % католиков среди населения;
Infant.Mortality – % детей, умерших до года.
Посмотрим на этот набор данных.

<<df1_1>>=
t <- swiss # встроенный набор данных по Швейцарии
glimpse(t)
@

Встроенный пакет graphics содержит функцию pairs, позволяющую получить все возможные диаграммы рассеяния на одном графике, а также выполнить их сглаживание с помощью опции panel.smooth:
<<df1_2, fig.width=7, fig.height=7>>=
pairs(swiss, panel = panel.smooth)
@

Функция cor позволяет как вычислить корреляцию между двумя выборками, так и получить корреляционную матрицу для всех переменных из набора данных.

<<df1_3>>=
options(digits = 5)
cor(swiss)
@

Существует еще одна функция, позволяющая получить корреляционную матрицу, диаграммы рассеяния и сглаженные распределения одновременно.
<<df1_5, fig.width=7, fig.height=7, message=FALSE, warning=FALSE>>=
library("GGally")
ggpairs(t) # функция из пакета GGally
@

Чтобы оценить регрессию рождаемости на остальные переменные, можно воспользоваться уже знакомой функцией lm, а регрессоры перечислить через знак «плюс»:
<<df1_6>>=
model2 <- lm(data=t, Fertility~Agriculture+Education+Catholic)
@
В данном случае регрессорами стали \% занятых в с/х, \% католического населения и \% имеющих образование выше начального.

Получить оценки коэффициентов уравнения регрессии, а также проверить основные гипотезы поможет функция summary:
<<df1_7>>=
summary(model2)
@
Построим прогноз по аналогии с парной линейной регрессией. Отличие заключается лишь в том, что в наборе данных необходимо указать значения каждого фактора.
<<df1_8>>=
# создаем новый набор данных
nd2 <- data.frame(Agriculture=0.5, Catholic=0.5, Education=20)
predict(model2,nd2)
@

Построение прогноза по нескольким точкам выполняется с помощью векторов значений.
<<df1_9>>=
# создаем новый набор данных
nd2 <- data.frame(Agriculture=c(0.5,0.8), Catholic=c(0.5, 0.65), Education=c(20, 25))
# прогнозируем
predict(model2,nd2)
@

\end{document}