В R есть куча встроенных наборов данных. Многие пакеты в себе помимо функций содержат также еще и примеры данных. Например, пакет ggplot2 содержит в себе данные по стоимости бриллиантов diamonds.
Список всех наборов данных можно узнать с помощью команды data()
data()
Мы воспользуемся набором про апельсиновые деревья.
h <- Orange
summary(h)
## Tree age circumference
## 3:7 Min. : 118 Min. : 30.0
## 1:7 1st Qu.: 484 1st Qu.: 65.5
## 5:7 Median :1004 Median :115.0
## 2:7 Mean : 922 Mean :115.9
## 4:7 3rd Qu.:1372 3rd Qu.:161.5
## Max. :1582 Max. :214.0
plot(h$age, h$circumference)
Более подробную информацию о наборе данных про Апельсины ищем у Чебурашки. Чебурашку можно позвать командой help(Orange).
Оценим обыкновенную регрессию
model <- lm(circumference ~ age, data = h) # оцениваем модель
report <- summary(model) # создаём отчет по модели
report # выводим отчет на экран
##
## Call:
## lm(formula = circumference ~ age, data = h)
##
## Residuals:
## Min 1Q Median 3Q Max
## -46.31 -14.95 -0.08 19.70 45.11
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.39965 8.62266 2.02 0.052 .
## age 0.10677 0.00828 12.90 1.9e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 23.7 on 33 degrees of freedom
## Multiple R-squared: 0.835, Adjusted R-squared: 0.83
## F-statistic: 166 on 1 and 33 DF, p-value: 1.93e-14
Из результатов оценивания можно вытащить
Ковариационную матрицу, \( \hat{\sigma}^2(X'X)^{-1} \):
vcov(model)
## (Intercept) age
## (Intercept) 74.35026 -0.0631691
## age -0.06317 0.0000685
Любимый всеми \( R^2 \):
report$r.squared
## [1] 0.8345
А также скорректированный \( R^2_{adj} \):
report$adj.r.squared
## [1] 0.8295
Коэффициенты:
model$coefficients
## (Intercept) age
## 17.3997 0.1068
Доверительные интервалы для коэффициентов:
confint(model, level = 0.9)
## 5 % 95 %
## (Intercept) 2.80700 31.9923
## age 0.09276 0.1208
Построение прогнозов с доверительным интервалом для среднего:
new.data <- data.frame(age = c(100, 200, 300))
predict(model, new.data, interval = "confidence")
## fit lwr upr
## 1 28.08 12.01 44.15
## 2 38.75 24.11 53.40
## 3 49.43 36.15 62.71
И с предиктивным интервалом:
predict(model, new.data, interval = "prediction")
## fit lwr upr
## 1 28.08 -22.8219 78.98
## 2 38.75 -11.7129 89.22
## 3 49.43 -0.6568 99.52
Что еще можно вытащить из оцененной модели, можно узнать несколькими способами. Например, стоит почитать help(summary.lm). Можно еще набрать в консоли model$ и пару раз кликнуть <Tab>. Еще имеет смысл попробовать нажать <Tab> после summ$. Да еще гугл есть.
А вообще бывает полезно изобрести велосипед. Пусть где-то в R уже есть функция, считающая доверительный интервал для прогноза. Если вместо этого написать эту функцию своими руками, то можно получить плюс один к экспириенсу и перейти на следующий уровень.