<jupyter_start><jupyter_text>Регрессия: [датасет Yacht Hydrodynamics](http://archive.ics.uci.edu/ml/datasets/Yacht+Hydrodynamics)Для парусных яхт нужно предсказать остаточное сопротивление на единицу массы смещения от размеров яхты и ее скорости.  Загрузка данных<jupyter_code>data <- read.table('yacht_hydrodynamics.data', header = TRUE, sep = ',')
data[1:5,]<jupyter_output><empty_output><jupyter_text>Разбиение данных на обучающую и тестовую выборки[createDataPartition](https://www.rdocumentation.org/packages/caret/versions/6.0-86/topics/createDataPartition) `(y,  times = 1,  p = 0.5,  list = TRUE,  groups = min(5, length(y)))` — функция для разделения данных. Требуется пакет 'caret'.Возвращает индексы объектов для обучения.* `y` — данные;* `times` — колличество разделений;* `p` — процент данных для train'а;* `list` — если TRUE, результат в виде списка, иначе матрица с floor(p * length(y)) строками и times столбцами;* `groups` — для числовых данных количество разбиений данных по квантилям для дальнейшего сэмплирования в каждом промежутке.<jupyter_code># установка пакета
install.packages('caret')
library(caret)
# разбиение данных
partition <- createDataPartition(data$Residuary_resistance, p = 0.7, list = FALSE)
train <- data[partition,]
test <- data[-partition,]<jupyter_output><empty_output><jupyter_text>Обучение[lm](https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/lm)`(formula, data, subset, weights, na.action,   method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,   singular.ok = TRUE, contrasts = NULL, offset, …)` — функция для построения линейных моделей.   * `formula` — формула;   * `data` — данные;   * `subset` — указывает на подмножество наблюдений, которые нужно использоваться для обучения;   * `weights` — веса для взвешенного МНК;   * `na.action` — функция, указывающая, что делать с пропусками.Возврвщает объект с полями:* `coefficients` — вектор коэффициентов* `residuals` — остатки моделиПолная справка:<jupyter_code>#?lm<jupyter_output><empty_output><jupyter_text>Обучаем на train модель вида $Residuary\_resistance  \sim \theta_1 + \theta_1 \cdot Froude_number$.<jupyter_code># обучение модели
model <- lm(formula = Residuary_resistance ~ Froude_number, data = train)
model<jupyter_output><empty_output><jupyter_text>Оценки параметров линейной регрессии:<jupyter_code># все коэффициенты
model$coefficients 

# взять первый коэффициент
model$coefficients[1]<jupyter_output><empty_output><jupyter_text>Ковариационная матрица вектора $\widehat \theta$ в условиях гомоскедастичности:<jupyter_code>vcov(model)<jupyter_output><empty_output><jupyter_text>Для записи формулы, с зависимостью от нескольких параметров, например $Residuary\_resistance  \sim \theta_1 + \theta_1 \cdot Froude_number + \theta_2 \cdot Prismatic\_coefficient$, можно записать названия признаков через `+`.<jupyter_code>model2 <- lm(formula = Residuary_resistance ~ Froude_number + Prismatic_coefficient, data = train)
model2<jupyter_output><empty_output><jupyter_text>Свойства (в гауссовской линейной модели)Функция [summary](https://www.rdocumentation.org/packages/base/versions/3.6.2/topics/summary)`(object, …)` печатает:* Остатки: минимум, 0.25-квантиль, медиана, 0.75-квантиль, максимум* Для каждого коэффициента: его оценка, стандартная ошибка, значение t-статистики гипотезы онезначимости коэффициента, pvalue этой гипотезы, звездочки значимости (чем больше, тем болеезначим коэффициент)* RSS и число степеней свободы, R^2 и его поправленная версия, значение F-статистики критерияФишера о значимости регрессии вообще, число степеней свободы распределения Фишера, pvalueэтой гипотезы.Для данной модели оба признака значимы, поскольку pvalue мало. То есть отвергаются гипотезы $\theta_1 = 0$ и $\theta_2 = 0$.<jupyter_code>summary(model)<jupyter_output><empty_output><jupyter_text>Можно вытащить отдельные числа:<jupyter_code>summary(model)$r.squared<jupyter_output><empty_output><jupyter_text>Доверительные интервалы для коэффициентов:<jupyter_code>confint(model, level = 0.95)<jupyter_output><empty_output><jupyter_text>ПредсказанияПредсказания строятся с помощью универсальной функции predict . Для линейной регрессии онаэквивалентна функции predict.lm .[predict](https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/predict.lm)`(object, newdata, se.fit = FALSE, scale = NULL, df = Inf,        interval = c("none", "confidence", "prediction"),        level = 0.95, type = c("response", "terms"),        terms = NULL, na.action = na.pass,        pred.var = res.var/weights, weights = 1, …)`* object — обученная ранее модель;* se.fit — индикатор того, что требуется подсчет стандартных ошибок;* scale — параметр для подсчета std.err;* df — степени свободы scale;* level — уровень доверия;* interval — тип интервала.По ней можно получить справку:<jupyter_code># ?predict.lm<jupyter_output><empty_output><jupyter_text>В предположениях гауссовской линейной модели можно построить два типа интервалов —доверительный ( confidence ) и предсказательный ( prediction ). Первый является доверительнымиинтервалом в обычном смысле для среднего значения отклика. Второй является интервалом, в которомс большой вероятностью лежит само значение отклика. Второй итервал всегда шире первого.Доверительный интервал для  $x_{new}\theta$ :$$x_{new}^T\hat\theta \pm T_{n -k - 1, 1 -  \alpha/2} \hat \sigma \sqrt{x_{new}^T(X^TX)^{-1}x_{new}}$$Предсказательный интервал для $x_{new} \theta + \varepsilon$:$$x_{new}^T\hat\theta \pm T_{n -k - 1, 1 - \alpha / 2} \hat \sigma \sqrt{1 + x_{new}^T(X^TX)^{-1}x_{new}}$$Предсказание значений на новых объектах вместе с доверительным интервалом:<jupyter_code>predicted <- predict(model, test, level = 0.95, interval = 'confidence')
predicted[1:3,]<jupyter_output><empty_output><jupyter_text>Предсказание значений на новых объектах вместе с предсказательным интервалом:<jupyter_code>predicted <- predict(model, test, level = 0.95, interval = 'prediction')
predicted[1:3,]<jupyter_output><empty_output><jupyter_text>MSE посчитаем ручками.<jupyter_code>mean((predicted[,1] - test$Residuary_resistance) ^ 2)<jupyter_output><empty_output><jupyter_text>Посмотрим на график предсказания.<jupyter_code>plot(data$Residuary_resistance ~ data$Froude_number)
x <- seq(from = 0, to = 0.5, by = 0.1)
lines(x, model$coefficients[1] + model$coefficients[2] * x, col = "red", lwd = 2)<jupyter_output><empty_output><jupyter_text>Еще примерыОбучаем на train модель вида $Residuary\_resistance = \theta_1 + \theta_2 \cdot Froude\_number + \theta_3 \cdot Froude\_number^2 + \theta_4 \cdot Froude\_number^3$*Обозначения в формуле:*`(x+y)^2` эквивалентно `x^2 + y^2 + xy`, что означает взять признаки `x^2`, `y^2`, `xy``I((x+y)^2)` означает взять признак `(x+y)^2`.<jupyter_code>model_2 <- lm(
    formula = Residuary_resistance ~ Froude_number + I(Froude_number^2) + I(Froude_number^3), 
    data = train
)
model_2$coefficients
summary(model_2)<jupyter_output><empty_output><jupyter_text>Сравнить две модели можно в одной таблице при помощи функции mtable.<jupyter_code>install.packages('memisc')
library("memisc")
print(mtable(model, model_2))<jupyter_output>Calls:
model: lm(formula = Residuary_resistance ~ Froude_number, data = train)
model_2: lm(formula = Residuary_resistance ~ Froude_number + I(Froude_number^2) + 
    I(Froude_number^3), data = train)

================================================
                        model       model_2     
------------------------------------------------
  (Intercept)         -25.098***    -47.514***  
                       (1.861)       (3.253)    
  Froude_number       123.299***    684.900***  
                       (6.091)      (38.674)    
  I(Froude_number^2)              -3105.799***  
                                   (141.931)    
  I(Froude_number^3)               4592.528***  
                                   (163.041)    
------------------------------------------------
  R-squared             0.657         0.985     
  N                   216           216         
================================================
  Significance: *** = p < 0.001;   
                ** = p < 0.[...]<jupyter_text>Значение MSE куда меньше во второй модели.<jupyter_code>predicted <- predict(model_2, test)
mean((predicted - test$Residuary_resistance) ^ 2)<jupyter_output><empty_output><jupyter_text>Посмотрим на график:<jupyter_code>plot(data$Residuary_resistance ~ data$Froude_number)
x <- seq(from = 0, to = 0.5, by = 0.01)
lines(x, 
      model_2$coefficients[1] + model_2$coefficients[2] * x + model_2$coefficients[3] * x^2 + model_2$coefficients[4] * x^3, 
      col = "red", lwd = 2)<jupyter_output><empty_output><jupyter_text>Обучим регрессию на всех фичах. Как видим, все остальные фичи незначимы - pvalue мало.*Чтобы взять в формулу все признаки, можно поставить просто точку*<jupyter_code>model <- lm(formula = Residuary_resistance ~ . + I(Froude_number^2) + I(Froude_number^3), 
            data = train)
model$coefficients
summary(model)<jupyter_output><empty_output><jupyter_text>MSE почти такое же.<jupyter_code>predicted <- predict(model, test)
mean((predicted - test$Residuary_resistance) ^ 2)<jupyter_output><empty_output><jupyter_text>Отбор признаковПроцедура отбора фичей [bestglm](https://cran.r-project.org/web/packages/bestglm/bestglm.pdf)`(Xy, family = gaussian, IC = "BIC", t = "default",CVArgs = "default", qLevel = 0.99, TopModels = 5,method = "exhaustive", intercept = TRUE, weights = NULL,nvmax = "default", RequireFullEnumerationQ = FALSE, ...)` использует алгоритм ’leaps’ для гауссовского случая.* `Xy` — данные;* `IC` — критерий информативности;<jupyter_code>install.packages('bestglm')
library('bestglm')

bestglm(train, family = gaussian, IC = "BIC")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14187384425724264, 'fraction_numerical': 0.021676597493322375, 'mean_word_length': 4.896426408237432, 'pattern_counts': {'":': 0, '<': 104, '<?xml version=': 0, '>': 88, 'https://': 5, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '232196', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3889, 'n_tokens_neox': 4035, 'n_words': 745, 'path': 'notebooks/MatStat6sem/R/3.linejjnajaregressijavr.ipynb', 'repo_name': 'TopCoder2K/DS-introduction', 'size': '9734'}