\documentclass[a4paper, 11pt]{article}\usepackage[]{graphicx}\usepackage[]{color}
\makeatletter
\def\maxwidth{
\ifdim\Gin@nat@width>\linewidth
\linewidth
\else
\Gin@nat@width
\fi
}
\makeatother
\definecolor{fgcolor}{rgb}{0.345, 0.345, 0.345}
\newcommand{\hlnum}[1]{\textcolor[rgb]{0.686,0.059,0.569}{#1}}
\newcommand{\hlstr}[1]{\textcolor[rgb]{0.192,0.494,0.8}{#1}}
\newcommand{\hlcom}[1]{\textcolor[rgb]{0.678,0.584,0.686}{\textit{#1}}}
\newcommand{\hlopt}[1]{\textcolor[rgb]{0,0,0}{#1}}
\newcommand{\hlstd}[1]{\textcolor[rgb]{0.345,0.345,0.345}{#1}}
\newcommand{\hlkwa}[1]{\textcolor[rgb]{0.161,0.373,0.58}{\textbf{#1}}}
\newcommand{\hlkwb}[1]{\textcolor[rgb]{0.69,0.353,0.396}{#1}}
\newcommand{\hlkwc}[1]{\textcolor[rgb]{0.333,0.667,0.333}{#1}}
\newcommand{\hlkwd}[1]{\textcolor[rgb]{0.737,0.353,0.396}{\textbf{#1}}}
\let\hlipl\hlkwb
\usepackage{framed}
\makeatletter
\newenvironment{kframe}{
\def\at@end@of@kframe{}
\ifinner\ifhmode
\def\at@end@of@kframe{\end{minipage}}
\begin{minipage}{\columnwidth}
\fi\fi
\def\FrameCommand##1{\hskip\@totalleftmargin \hskip-\fboxsep
\colorbox{shadecolor}{##1}\hskip-\fboxsep
\hskip-\linewidth \hskip-\@totalleftmargin \hskip\columnwidth}
\MakeFramed {\advance\hsize-\width
\@totalleftmargin\z@ \linewidth\hsize
\@setminipage}}
{\par\unskip\endMakeFramed
\at@end@of@kframe}
\makeatother
\definecolor{shadecolor}{rgb}{.97, .97, .97}
\definecolor{messagecolor}{rgb}{0, 0, 0}
\definecolor{warningcolor}{rgb}{1, 0, 1}
\definecolor{errorcolor}{rgb}{1, 0, 0}
\newenvironment{knitrout}{}{}
\usepackage{alltt}
\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}
\usepackage{etoolbox}
\usepackage{enumitem}
\usepackage{array,tabularx,tabulary,booktabs}
\usepackage{longtable}
\usepackage{multirow}
\usepackage{lastpage}
\usepackage{url}
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
\begin{document}
\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). Построим график зависимости длины тормозного пути от скорости автомобиля.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{d} \hlkwb{<-} \hlstd{cars}
\hlkwd{qplot}\hlstd{(}\hlkwc{data}\hlstd{=d, speed, dist)} \hlopt{+}
\hlkwd{geom_point}\hlstd{(}\hlkwd{aes}\hlstd{(}\hlkwc{x}\hlstd{=d}\hlopt{$}\hlstd{speed,} \hlkwc{y}\hlstd{=d}\hlopt{$}\hlstd{dist),} \hlkwc{size} \hlstd{=} \hlnum{1}\hlstd{)} \hlopt{+} \hlkwd{theme_bw}\hlstd{(}\hlkwc{base_size} \hlstd{=} \hlnum{12}\hlstd{)} \hlopt{+}
\hlkwd{xlab}\hlstd{(}\hlstr{"speed, miles/h"}\hlstd{)} \hlopt{+} \hlkwd{ylab}\hlstd{(}\hlstr{"distance, feet"}\hlstd{)} \hlopt{+}
\hlkwd{labs}\hlstd{(}\hlkwc{title} \hlstd{=} \hlstr{"scatter plot"}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/df1-1}
\end{knitrout}
Оценим модель линейной регрессии длины тормозного пути на скорость автомобиля.
Для этого командой lm поместим в переменную model модель линейной регрессии, указав dist в качестве зависимой переменной, и через значок ~ переменную speed в качестве регрессора.
Тип lm представляет собой список из 12 элементов. Посмотрим на коэффициенты уравнения линейной регрессии.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{model} \hlkwb{<-} \hlkwd{lm}\hlstd{(}\hlkwc{data}\hlstd{=d, dist}\hlopt{~}\hlstd{speed)} \hlcom{# базовый пакет stats}
\hlstd{model}\hlopt{$}\hlstd{coefficients}
\end{alltt}
\begin{verbatim}
## (Intercept) speed
## -17.579095 3.932409
\end{verbatim}
\end{kframe}
\end{knitrout}
(Intercept) – это константа в уравнении регрессии, speed – коэффициент регрессии.
Таким образом, уравнение регрессии имеет вид:
\[dist_i^m = -17.579 + 3.9324 \cdot speed_i\]
Так же можно посмотреть значения вектора ошибок модели – разницу между реальной длиной тормозного пути dist и полученной по модели $dist_i^m$. Выведем первые 10 значений этого вектора с точностью две цифры после запятой. Более полный набор расчетов по модели можно получить командой summary.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{options}\hlstd{(}\hlkwc{digits} \hlstd{=} \hlnum{3}\hlstd{)}
\hlstd{model}\hlopt{$}\hlstd{residuals[}\hlnum{1}\hlopt{:}\hlnum{10}\hlstd{]}
\end{alltt}
\begin{verbatim}
## 1 2 3 4 5 6 7 8 9 10
## 3.85 11.85 -5.95 12.05 2.12 -7.81 -3.74 4.26 12.26 -8.68
\end{verbatim}
\begin{alltt}
\hlkwd{summary}\hlstd{(model)} \hlcom{# базовый пакет base}
\end{alltt}
\begin{verbatim}
##
## Call:
## lm(formula = dist ~ speed, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.07 -9.53 -2.27 9.21 43.20
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17.579 6.758 -2.60 0.012 *
## speed 3.932 0.416 9.46 1.5e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.4 on 48 degrees of freedom
## Multiple R-squared: 0.651, Adjusted R-squared: 0.644
## F-statistic: 89.6 on 1 and 48 DF, p-value: 1.49e-12
\end{verbatim}
\end{kframe}
\end{knitrout}
Помимо коэффициентов регрессии, 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\% доверительными интервалами.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{qplot}\hlstd{(}\hlkwc{data} \hlstd{= d, speed, dist)} \hlopt{+} \hlkwd{stat_smooth}\hlstd{(}\hlkwc{method}\hlstd{=}\hlstr{"lm"}\hlstd{,} \hlkwc{level} \hlstd{=} \hlnum{0.95}\hlstd{)} \hlopt{+}
\hlkwd{theme_bw}\hlstd{(}\hlkwc{base_size} \hlstd{=} \hlnum{12}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/df4-1}
\end{knitrout}
Рассчитаем также 95\% доверительные интервалы для параметров линейной регрессии. Полученные по модели значения можно вывести командой fitted.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{confint}\hlstd{(model,} \hlkwc{level} \hlstd{=} \hlnum{0.95}\hlstd{)} \hlcom{# базовый пакет stats}
\end{alltt}
\begin{verbatim}
## 2.5
## (Intercept) -31.2 -3.99
## speed 3.1 4.77
\end{verbatim}
\begin{alltt}
\hlkwd{options}\hlstd{(}\hlkwc{digits}\hlstd{=}\hlnum{4}\hlstd{)}
\hlkwd{fitted}\hlstd{(model)} \hlcom{# базовый пакет stats}
\end{alltt}
\begin{verbatim}
## 1 2 3 4 5 6 7 8 9 10 11 12
## -1.849 -1.849 9.948 9.948 13.880 17.813 21.745 21.745 21.745 25.677 25.677 29.610
## 13 14 15 16 17 18 19 20 21 22 23 24
## 29.610 29.610 29.610 33.542 33.542 33.542 33.542 37.475 37.475 37.475 37.475 41.407
## 25 26 27 28 29 30 31 32 33 34 35 36
## 41.407 41.407 45.339 45.339 49.272 49.272 49.272 53.204 53.204 53.204 53.204 57.137
## 37 38 39 40 41 42 43 44 45 46 47 48
## 57.137 57.137 61.069 61.069 61.069 61.069 61.069 68.934 72.866 76.799 76.799 76.799
## 49 50
## 76.799 80.731
\end{verbatim}
\end{kframe}
\end{knitrout}
Рассчитать необъясненную сумму квадратов отклонений и полную сумму квадратов можно, воспользовавшись функцией deviance и уже известными функциями sum и mean. Для того, чтобы построить прогноз по полученной модели, нужно задать значения регрессора и поместить их в новый data.frame.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{RSS} \hlkwb{<-} \hlkwd{deviance}\hlstd{(model)} \hlcom{# базовый пакет stats}
\hlstd{TSS} \hlkwb{<-} \hlkwd{sum}\hlstd{((d}\hlopt{$}\hlstd{dist}\hlopt{-}\hlkwd{mean}\hlstd{(d}\hlopt{$}\hlstd{dist))}\hlopt{^}\hlnum{2}\hlstd{)}
\hlstd{RSS; TSS}
\end{alltt}
\begin{verbatim}
## [1] 11354
## [1] 32539
\end{verbatim}
\begin{alltt}
\hlcom{# создаем новый набор данных}
\hlstd{nd} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(}\hlkwc{speed}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{40}\hlstd{,}\hlnum{60}\hlstd{))}
\hlcom{# Строим прогноз функцией predict}
\hlkwd{predict}\hlstd{(model,nd)}
\end{alltt}
\begin{verbatim}
## 1 2
## 139.7 218.4
\end{verbatim}
\end{kframe}
\end{knitrout}
\section*{Множественная линейная регрессия}
Рассмотрим встроенный набор данных по социально-экономическим показателям в 47 провинциях Швейцарии в 1888 г. Этот набор данных содержит 6 переменных по 47 наблюдений, каждая из которых измеряется в процентах (help(swiss)):
Fertility – рождаемость;
Agriculture –
Examination –
Education –
Catholic –
Infant.Mortality –
Посмотрим на этот набор данных.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{t} \hlkwb{<-} \hlstd{swiss} \hlcom{# встроенный набор данных по Швейцарии}
\hlkwd{glimpse}\hlstd{(t)}
\end{alltt}
\begin{verbatim}
## Observations: 47
## Variables: 6
## $ Fertility <dbl> 80.2, 83.1, 92.5, 85.8, 76.9, 76.1, 83.8, 92.4, 82.4, 82.9, ...
## $ Agriculture <dbl> 17.0, 45.1, 39.7, 36.5, 43.5, 35.3, 70.2, 67.8, 53.3, 45.2, ...
## $ Examination <int> 15, 6, 5, 12, 17, 9, 16, 14, 12, 16, 14, 21, 14, 19, 22, 18,...
## $ Education <int> 12, 9, 5, 7, 15, 7, 7, 8, 7, 13, 6, 12, 7, 12, 5, 2, 8, 28, ...
## $ Catholic <dbl> 9.96, 84.84, 93.40, 33.77, 5.16, 90.57, 92.85, 97.16, 97.67,...
## $ Infant.Mortality <dbl> 22.2, 22.2, 20.2, 20.3, 20.6, 26.6, 23.6, 24.9, 21.0, 24.4, ...
\end{verbatim}
\end{kframe}
\end{knitrout}
Встроенный пакет graphics содержит функцию pairs, позволяющую получить все возможные диаграммы рассеяния на одном графике, а также выполнить их сглаживание с помощью опции panel.smooth:
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{pairs}\hlstd{(swiss,} \hlkwc{panel} \hlstd{= panel.smooth)}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/df1_2-1}
\end{knitrout}
Функция cor позволяет как вычислить корреляцию между двумя выборками, так и получить корреляционную матрицу для всех переменных из набора данных.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{options}\hlstd{(}\hlkwc{digits} \hlstd{=} \hlnum{5}\hlstd{)}
\hlkwd{cor}\hlstd{(swiss)}
\end{alltt}
\begin{verbatim}
## Fertility Agriculture Examination Education Catholic Infant.Mortality
## Fertility 1.00000 0.353079 -0.64588 -0.663789 0.46368 0.416556
## Agriculture 0.35308 1.000000 -0.68654 -0.639523 0.40110 -0.060859
## Examination -0.64588 -0.686542 1.00000 0.698415 -0.57274 -0.114022
## Education -0.66379 -0.639523 0.69842 1.000000 -0.15386 -0.099322
## Catholic 0.46368 0.401095 -0.57274 -0.153859 1.00000 0.175496
## Infant.Mortality 0.41656 -0.060859 -0.11402 -0.099322 0.17550 1.000000
\end{verbatim}
\end{kframe}
\end{knitrout}
Существует еще одна функция, позволяющая получить корреляционную матрицу, диаграммы рассеяния и сглаженные распределения одновременно.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{library}\hlstd{(}\hlstr{"GGally"}\hlstd{)}
\hlkwd{ggpairs}\hlstd{(t)} \hlcom{# функция из пакета GGally}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/df1_5-1}
\end{knitrout}
Чтобы оценить регрессию рождаемости на остальные переменные, можно воспользоваться уже знакомой функцией lm, а регрессоры перечислить через знак «плюс»:
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{model2} \hlkwb{<-} \hlkwd{lm}\hlstd{(}\hlkwc{data}\hlstd{=t, Fertility}\hlopt{~}\hlstd{Agriculture}\hlopt{+}\hlstd{Education}\hlopt{+}\hlstd{Catholic)}
\end{alltt}
\end{kframe}
\end{knitrout}
В данном случае регрессорами стали \% занятых в с/х, \% католического населения и \% имеющих образование выше начального.
Получить оценки коэффициентов уравнения регрессии, а также проверить основные гипотезы поможет функция summary:
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{summary}\hlstd{(model2)}
\end{alltt}
\begin{verbatim}
##
## Call:
## lm(formula = Fertility ~ Agriculture + Education + Catholic,
## data = t)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.178 -6.548 1.379 5.822 14.840
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 86.22502 4.73472 18.211 < 2e-16 ***
## Agriculture -0.20304 0.07115 -2.854 0.00662 **
## Education -1.07215 0.15580 -6.881 1.91e-08 ***
## Catholic 0.14520 0.03015 4.817 1.84e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.728 on 43 degrees of freedom
## Multiple R-squared: 0.6423, Adjusted R-squared: 0.6173
## F-statistic: 25.73 on 3 and 43 DF, p-value: 1.089e-09
\end{verbatim}
\end{kframe}
\end{knitrout}
Построим прогноз по аналогии с парной линейной регрессией. Отличие заключается лишь в том, что в наборе данных необходимо указать значения каждого фактора.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{# создаем новый набор данных}
\hlstd{nd2} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(}\hlkwc{Agriculture}\hlstd{=}\hlnum{0.5}\hlstd{,} \hlkwc{Catholic}\hlstd{=}\hlnum{0.5}\hlstd{,} \hlkwc{Education}\hlstd{=}\hlnum{20}\hlstd{)}
\hlkwd{predict}\hlstd{(model2,nd2)}
\end{alltt}
\begin{verbatim}
## 1
## 64.8
\end{verbatim}
\end{kframe}
\end{knitrout}
Построение прогноза по нескольким точкам выполняется с помощью векторов значений.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{# создаем новый набор данных}
\hlstd{nd2} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(}\hlkwc{Agriculture}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{0.5}\hlstd{,}\hlnum{0.8}\hlstd{),} \hlkwc{Catholic}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{0.5}\hlstd{,} \hlnum{0.65}\hlstd{),} \hlkwc{Education}\hlstd{=}\hlkwd{c}\hlstd{(}\hlnum{20}\hlstd{,} \hlnum{25}\hlstd{))}
\hlcom{# прогнозируем}
\hlkwd{predict}\hlstd{(model2,nd2)}
\end{alltt}
\begin{verbatim}
## 1 2
## 64.8 59.4
\end{verbatim}
\end{kframe}
\end{knitrout}
\end{document}