Проверим условия для получения состоятельных, несмещенных, эффективных оценок
Проверка случайности остаточной компоненты. Построим график зависимости остатков от спрогнозированных значений модели
qplot(y = model_1$fitted.values, x = model_1$residuals, xlab = "Остаточные степени свободы, Ei", ylab="Спрогнозированные значения, y^")
Проверим условие M(Ei) = 0.
H0: M(Ei) = 0(математическое ожидание остатков равно нулю); H1: M(Ei) != 0(математическое ожидание остатков отлично от нуля).
mean(model_1$residuals)
## [1] 18.66859
Вычисляем t расч.:
a<-mean(model_1$residuals) b<-sd(model_1$residuals, na.rm = FALSE) n <- sqrt(nrow(train)) tm <- a/b*n str(tm)
## num 1.1
Вычисляем t табл. при p=0.95 и степени свободы N-1:
df <- nrow(train) qt(0.95, df - 1)
## [1] 1.671553
|t расч.| > t табл. => отвергаем H0 с вероятностью p=0.95.
Протестируем на наличие гетероскедастичности в остатках регрессионное уравнение с помощью теста Бройша-Погана: H0: D(Ei) = sig^2 (отсутствие гетероскедастичности, наличие гомоскедастичности) H1: D(Ei) != sig^2 (наличие гетероскедастичности)
#bptest(model_1)
Так как значение p-value больше 0.05, нулевая ги потеза об отсутствие гетероскедастичности остатков принимается.
Наличие гетероскедастичности в остатках регрессии можно проверить с помощью теста ранговой корреляции Спирмена.
Суть теста заключается в определении наличия связи между ростом остаточной компоненты и ростом независимого фактора, то есть определение роста дисперсии остатков. Такая зависимость проверяется на основе расчета коэффициента ранговой корреляции Спирмена ρ между остатками модели E и независимым фактором liquidated_2020_01.
cor.test (train$liquidated_2020_01, model_1$residuals, method = "spearman")
## Warning in cor.test.default(train$liquidated_2020_01, model_1$residuals, : ## Cannot compute exact p-value with ties
## ## Spearman's rank correlation rho ## ## data: train$liquidated_2020_01 and model_1$residuals ## S = 23322, p-value = 0.01396 ## alternative hypothesis: true rho is not equal to 0 ## sample estimates: ## rho ## 0.318478
Так как значение p-value меньше 0.05, нулевая гипотеза об отсутствие гетероскедастичности остатков отвергается.
Проверим условие cov(Ei,Ej)=0, i!=j
Н0: rho=0 (т.е. автокорреляция остатков отсутствует); H1: rho>0 или rho<0 (наличие положительной или отрицательной автокорреляции остатков).
# тест Дарбина-Уотсона bgtest(model_1)
## ## Breusch-Godfrey test for serial correlation of order up to 1 ## ## data: model_1 ## LM test = 0.034515, df = 1, p-value = 0.8526
# тест Бройша- Годфри dwtest(model_1)
## ## Durbin-Watson test ## ## data: model_1 ## DW = 1.9223, p-value = 0.4048 ## alternative hypothesis: true autocorrelation is greater than 0
Так как значение p-value больше 0.05 => гипотеза H0 об отсутствии автокорреляции остатков принимается
Проверим условие Ei ~ N(0, sig^2), т.е. согласуются ли остатки регрессии с нормальным законом:
H0: остатки регрессии согласуются с нормальным законом распределения H1: остатки регрессии не согласуются с нормальным законом распределения
Рассмотрим графические тесты:
qqnorm(model_1$residuals)
library(sm); Z <- model_1$residuals hist(Z, main = "Гистограмма остатков")
sm.density(model_1$residuals, model = "Normal",xlab = "Остатки", ylab = "Функция плотности распределения")
Рассмотрим параметрические тесты:
library(nortest) # тест Lilliefors (Kolmogorov-Smirnov) lillie.test(model_1$residuals)
## ## Lilliefors (Kolmogorov-Smirnov) normality test ## ## data: model_1$residuals ## D = 0.17603, p-value = 9.643e-05
Так как значение p-value меньше 0.05, нулевая гипотеза о согласии распределения остатков с нормальным законом распределения отвергается.
Рассмотрим также тест Шапиро-Уилка.
# тест Шапиро-Уилка sh <- model_1$residuals #View(sh) shapiro.test(sh)
## ## Shapiro-Wilk normality test ## ## data: sh ## W = 0.90312, p-value = 0.0001938
Так как значение p-value меньше 0.05, нулевая гипотеза о согласии распределения остатков с нормальным законом распределения отвергается.
Прогнозирование значений по полученной модели на тестовой выборке:
plot(data$registered_2020_01, data$liquidated_2020_01, ylab="Ликвидированые организации", xlab="Зарегистрированные организации") lines(data$registered_2020_01, predict(model_1, data), col = 'red')
Цель
Провести множественный регрессионный анализ между максимальным числом сердечных сокращений (thalach), который зависит от следующих параметров:
age - возраст
sex - пол
cp - тип боли в груди
trestbps - артериальное давление в состоянии покоя
chol - уровень холестерола в сыворотке крови
restecg - результаты электрокардиографии в состоянии покоя
exang - наличие стенокардии
oldpeak - депрессия сегмента ST сердца
Загрузим необходимые пакеты
library(openxlsx) library(caTools) library(DescTools) library(dplyr) library(ggplot2) library(glmnet)
## Loading required package: Matrix
## Loaded glmnet 4.1-1
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v tibble 3.1.1 v purrr 0.3.4 ## v tidyr 1.1.3 v stringr 1.4.0 ## v readr 1.4.0 v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() -- ## x purrr::%@%() masks memisc::%@%() ## x dplyr::collect() masks memisc::collect() ## x tidyr::expand() masks Matrix::expand() ## x dplyr::filter() masks stats::filter() ## x dplyr::lag() masks stats::lag() ## x tidyr::pack() masks Matrix::pack() ## x dplyr::recode() masks memisc::recode() ## x dplyr::rename() masks memisc::rename() ## x dplyr::select() masks MASS::select() ## x dplyr::syms() masks memisc::syms(), ggplot2::syms() ## x tidyr::unpack() masks Matrix::unpack() ## x tibble::view() masks memisc::view()
Загрузим набор данных. Источник - https://www.kaggle.com/imnikhilanand/heart-attack-prediction
data <- read.csv("data.csv", sep = ";")
Для оценки качества регресионной модели разделим выборку на тестовую и обучающую
set.seed(1821) split <- sample.split(data$thalach, SplitRatio = 0.75) train <- subset(data, split == TRUE) test <- subset(data, split == FALSE)
Построим модель линейной регрессии (для оценки параметров применяется метод наименьших квадратов)
m_lm <- lm(thalach ~ ., data = train) summary(m_lm)
## ## Call: ## lm(formula = thalach ~ ., data = train) ## ## Residuals: ## Min 1Q Median 3Q Max ## -64.279 -12.048 0.885 14.177 39.009 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 2.191e+02 1.351e+01 16.215 < 2e-16 *** ## age -1.061e+00 1.791e-01 -5.921 1.38e-08 *** ## sex 1.311e+00 3.026e+00 0.433 0.6652 ## cp -5.889e+00 1.635e+00 -3.601 0.0004 *** ## trestbps -6.435e-02 7.930e-02 -0.811 0.4180 ## chol -6.939e-05 2.102e-02 -0.003 0.9974 ## restecg 1.309e+00 2.728e+00 0.480 0.6318 ## exang -8.930e+00 4.465e+00 -2.000 0.0468 * ## oldpeak -2.502e+00 2.216e+00 -1.129 0.2602 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 19.2 on 200 degrees of freedom ## Multiple R-squared: 0.3725, Adjusted R-squared: 0.3474 ## F-statistic: 14.84 on 8 and 200 DF, p-value: < 2.2e-16
Проверим линейную связь F-критерий H0: a=b 0 , (т.е. линейная связь между x и y отсутствует); H1: a>0 или b>0 (т.е. наличие линейной связи).
mean_model <- mean(train$thalach) f <- sum((m_lm$fitted.values - mean_model)^2) / summary(m_lm)$sigma^2 f
## [1] 118.7308
f_tabl <- qf(0.95, 1, nrow(train)-2) f_tabl
## [1] 3.886774
F расч | > F табл => отвергаем H0 в пользу наличия линейной связи с вероятностью 0.95.
Находим коэффициенты детерминации:
#коэффициент множественной корреляции summary(m_lm)$r.squared
## [1] 0.3725113
#квадрат коэфициента скоректированной корреляции (доля дисперсии зависимой переменной price) summary(m_lm)$adj.r.squared
## [1] 0.3474117
Значения коэффициентов детерминации довольно низки (ниже 0.5), подгонку построенной модели множественной регресии можно считать некачественной
Находим остаточную дисперсию:
summary(m_lm)$sigma^2
## [1] 368.6646
Средняя ошибка аппроксимации
A <- sum(abs((train$thalach - m_lm$fitted.values)/train$thalach)) / nrow(train) * 100 str(A)
## num 11.7
А > 5-7% => качество модели плохое
Построим прогноз по полученной модели на тестовой выборке
test$predict <- predict(m_lm, test) ggplot(test, aes(x = thalach, y = predict)) + geom_point() + geom_abline() + scale_x_continuous(limits = c(75, 200), expand = c(0, 0)) + scale_y_continuous(limits = c(75, 200), expand = c(0, 0))
## Warning: Removed 1 rows containing missing values (geom_abline).
Построим модель линейной регрессии, исключив незначительный параметр chol
m_lm <- lm(thalach ~ age + sex + cp + trestbps + restecg + exang + oldpeak, data = train) summary(m_lm)
## ## Call: ## lm(formula = thalach ~ age + sex + cp + trestbps + restecg + ## exang + oldpeak, data = train) ## ## Residuals: ## Min 1Q Median 3Q Max ## -64.282 -12.042 0.884 14.175 39.004 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 219.05567 13.07379 16.755 < 2e-16 *** ## age -1.06073 0.17781 -5.965 1.08e-08 *** ## sex 1.31167 3.01820 0.435 0.664328 ## cp -5.88917 1.63036 -3.612 0.000384 *** ## trestbps -0.06437 0.07888 -0.816 0.415467 ## restecg 1.30835 2.71386 0.482 0.630259 ## exang -8.93161 4.40779 -2.026 0.044053 * ## oldpeak -2.50171 2.20882 -1.133 0.258731 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 19.15 on 201 degrees of freedom ## Multiple R-squared: 0.3725, Adjusted R-squared: 0.3507 ## F-statistic: 17.05 on 7 and 201 DF, p-value: < 2.2e-16
Проверим линейную связь F-критерий H0: a=b 0 , (т.е. линейная связь между x и y отсутствует); H1: a>0 или b>0 (т.е. наличие линейной связи).
mean_model <- mean(train$thalach) f <- sum((m_lm$fitted.values - mean_model)^2) / summary(m_lm)$sigma^2 f
## [1] 119.3245
f_tabl <- qf(0.95, 1, nrow(train)-2) f_tabl
## [1] 3.886774
F расч | > F табл => отвергаем H0 в пользу наличия линейной связи с вероятностью 0.95.
Находим коэффициенты детерминации:
#коэффициент множественной корреляции summary(m_lm)$r.squared
## [1] 0.3725112
#квадрат коэфициента скоректированной корреляции (доля дисперсии зависимой переменной price) summary(m_lm)$adj.r.squared
## [1] 0.3506584
Значения коэффициентов детерминации довольно низки (ниже 0.5), подгонку построенной модели множественной регресии можно считать некачественной
Находим остаточную дисперсию:
summary(m_lm)$sigma^2
## [1] 366.8305
Средняя ошибка аппроксимации
A <- sum(abs((train$thalach - m_lm$fitted.values)/train$thalach)) / nrow(train) * 100 str(A)
## num 11.7
А > 5-7% => качество модели плохое
Построим прогноз по полученной модели на тестовой выборке
test$predict <- predict(m_lm, test) ggplot(test, aes(x = thalach, y = predict)) + geom_point() + geom_abline() + scale_x_continuous(limits = c(75, 200), expand = c(0, 0)) + scale_y_continuous(limits = c(75, 200), expand = c(0, 0))
## Warning: Removed 1 rows containing missing values (geom_abline).
Построим модель линейной регрессии, исключив следующий незначительный параметр sex
m_lm <- lm(thalach ~ age + cp + trestbps + restecg + exang + oldpeak, data = train) summary(m_lm)
## ## Call: ## lm(formula = thalach ~ age + cp + trestbps + restecg + exang + ## oldpeak, data = train) ## ## Residuals: ## Min 1Q Median 3Q Max ## -63.727 -12.742 1.106 14.423 39.062 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 219.47518 13.01191 16.867 < 2e-16 *** ## age -1.06531 0.17714 -6.014 8.36e-09 *** ## cp -5.77391 1.60541 -3.597 0.000405 *** ## trestbps -0.06150 0.07845 -0.784 0.433971 ## restecg 1.21115 2.69919 0.449 0.654123 ## exang -9.01129 4.39512 -2.050 0.041627 * ## oldpeak -2.41331 2.19501 -1.099 0.272880 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 19.11 on 202 degrees of freedom ## Multiple R-squared: 0.3719, Adjusted R-squared: 0.3533 ## F-statistic: 19.94 on 6 and 202 DF, p-value: < 2.2e-16
Проверим линейную связь F-критерий H0: a=b 0 , (т.е. линейная связь между x и y отсутствует); H1: a>0 или b>0 (т.е. наличие линейной связи).
mean_model <- mean(train$thalach) f <- sum((m_lm$fitted.values - mean_model)^2) / summary(m_lm)$sigma^2 f
## [1] 119.6159
f_tabl <- qf(0.95, 1, nrow(train)-2) f_tabl
## [1] 3.886774
F расч | > F табл => отвергаем H0 в пользу наличия линейной связи с вероятностью 0.95.
Находим коэффициенты детерминации:
#коэффициент множественной корреляции summary(m_lm)$r.squared
## [1] 0.3719216
#квадрат коэфициента скоректированной корреляции (доля дисперсии зависимой переменной price) summary(m_lm)$adj.r.squared
## [1] 0.3532658
Значения коэффициентов детерминации довольно низки (ниже 0.5), подгонку построенной модели множественной регресии можно считать некачественной
Находим остаточную дисперсию:
summary(m_lm)$sigma^2
## [1] 365.3574
Средняя ошибка аппроксимации
A <- sum(abs((train$thalach - m_lm$fitted.values)/train$thalach)) / nrow(train) * 100 str(A)
## num 11.7
А > 5-7% => качество модели плохое
Построим прогноз по полученной модели на тестовой выборке
test$predict <- predict(m_lm, test) ggplot(test, aes(x = thalach, y = predict)) + geom_point() + geom_abline() + scale_x_continuous(limits = c(75, 200), expand = c(0, 0)) + scale_y_continuous(limits = c(75, 200), expand = c(0, 0))
## Warning: Removed 1 rows containing missing values (geom_abline).
Построим модель линейной регрессии, исключив следующий незначительный параметр restecg
m_lm <- lm(thalach ~ age + cp + trestbps + exang + oldpeak, data = train) summary(m_lm)
## ## Call: ## lm(formula = thalach ~ age + cp + trestbps + exang + oldpeak, ## data = train) ## ## Residuals: ## Min 1Q Median 3Q Max ## -62.764 -12.974 1.058 14.696 38.788 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 219.64576 12.98074 16.921 < 2e-16 *** ## age -1.06003 0.17640 -6.009 8.51e-09 *** ## cp -5.78974 1.60186 -3.614 0.00038 *** ## trestbps -0.06216 0.07828 -0.794 0.42806 ## exang -8.90686 4.38032 -2.033 0.04332 * ## oldpeak -2.45687 2.18855 -1.123 0.26293 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 19.08 on 203 degrees of freedom ## Multiple R-squared: 0.3713, Adjusted R-squared: 0.3558 ## F-statistic: 23.98 on 5 and 203 DF, p-value: < 2.2e-16
Проверим линейную связь F-критерий H0: a=b 0 , (т.е. линейная связь между x и y отсутствует); H1: a>0 или b>0 (т.е. наличие линейной связи).
mean_model <- mean(train$thalach) f <- sum((m_lm$fitted.values - mean_model)^2) / summary(m_lm)$sigma^2 f
## [1] 119.8862
f_tabl <- qf(0.95, 1, nrow(train)-2) f_tabl
## [1] 3.886774
F расч | > F табл => отвергаем H0 в пользу наличия линейной связи с вероятностью 0.95.
Находим коэффициенты детерминации:
#коэффициент множественной корреляции summary(m_lm)$r.squared
## [1] 0.3712956
#квадрат коэфициента скоректированной корреляции (доля дисперсии зависимой переменной price) summary(m_lm)$adj.r.squared
## [1] 0.3558103
Значения коэффициентов детерминации довольно низки (ниже 0.5), подгонку построенной модели множественной регресии можно считать некачественной
Находим остаточную дисперсию:
summary(m_lm)$sigma^2
## [1] 363.92
Средняя ошибка аппроксимации
A <- sum(abs((train$thalach - m_lm$fitted.values)/train$thalach)) / nrow(train) * 100 str(A)
## num 11.7
А > 5-7% => качество модели плохое
Построим прогноз по полученной модели на тестовой выборке
test$predict <- predict(m_lm, test) ggplot(test, aes(x = thalach, y = predict)) + geom_point() + geom_abline() + scale_x_continuous(limits = c(75, 200), expand = c(0, 0)) + scale_y_continuous(limits = c(75, 200), expand = c(0, 0))
## Warning: Removed 1 rows containing missing values (geom_abline).
Построим модель линейной регрессии, исключив следующий незначительный параметр trestbps
m_lm <- lm(thalach ~ age + cp + exang + oldpeak, data = train) summary(m_lm)
## ## Call: ## lm(formula = thalach ~ age + cp + exang + oldpeak, data = train) ## ## Residuals: ## Min 1Q Median 3Q Max ## -62.948 -12.818 1.443 14.557 38.257 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 212.4188 9.2477 22.970 < 2e-16 *** ## age -1.0819 0.1741 -6.215 2.84e-09 *** ## cp -5.7203 1.5980 -3.580 0.00043 *** ## exang -9.1658 4.3642 -2.100 0.03694 * ## oldpeak -2.6810 2.1683 -1.236 0.21771 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 19.06 on 204 degrees of freedom ## Multiple R-squared: 0.3693, Adjusted R-squared: 0.357 ## F-statistic: 29.87 on 4 and 204 DF, p-value: < 2.2e-16
Проверим линейную связь F-критерий H0: a=b 0 , (т.е. линейная связь между x и y отсутствует); H1: a>0 или b>0 (т.е. наличие линейной связи).
mean_model <- mean(train$thalach) f <- sum((m_lm$fitted.values - mean_model)^2) / summary(m_lm)$sigma^2 f
## [1] 119.472
f_tabl <- qf(0.95, 1, nrow(train)-2) f_tabl
## [1] 3.886774
F расч | > F табл => отвергаем H0 в пользу наличия линейной связи с вероятностью 0.95.
Находим коэффициенты детерминации:
#коэффициент множественной корреляции summary(m_lm)$r.squared
## [1] 0.3693425
#квадрат коэфициента скоректированной корреляции (доля дисперсии зависимой переменной price) summary(m_lm)$adj.r.squared
## [1] 0.3569767
Значения коэффициентов детерминации довольно низки (ниже 0.5), подгонку построенной модели множественной регресии можно считать некачественной
Находим остаточную дисперсию:
summary(m_lm)$sigma^2
## [1] 363.2611
Средняя ошибка аппроксимации
A <- sum(abs((train$thalach - m_lm$fitted.values)/train$thalach)) / nrow(train) * 100 str(A)
## num 11.7
А > 5-7% => качество модели плохое
Построим прогноз по полученной модели на тестовой выборке
test$predict <- predict(m_lm, test) ggplot(test, aes(x = thalach, y = predict)) + geom_point() + geom_abline() + scale_x_continuous(limits = c(75, 200), expand = c(0, 0)) + scale_y_continuous(limits = c(75, 200), expand = c(0, 0))
## Warning: Removed 1 rows containing missing values (geom_abline).
Построим модель линейной регрессии, исключив следующий незначительный параметр oldpeak
m_lm <- lm(thalach ~ age + cp + exang, data = train) summary(m_lm)
## ## Call: ## lm(formula = thalach ~ age + cp + exang, data = train) ## ## Residuals: ## Min 1Q Median 3Q Max ## -62.72 -13.53 1.08 14.21 34.90 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 212.7708 9.2553 22.989 < 2e-16 *** ## age -1.0929 0.1741 -6.279 2.01e-09 *** ## cp -5.8081 1.5985 -3.633 0.000353 *** ## exang -12.6927 3.3073 -3.838 0.000165 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 19.08 on 205 degrees of freedom ## Multiple R-squared: 0.3646, Adjusted R-squared: 0.3553 ## F-statistic: 39.21 on 3 and 205 DF, p-value: < 2.2e-16
Проверим линейную связь F-критерий H0: a=b 0 , (т.е. линейная связь между x и y отсутствует); H1: a>0 или b>0 (т.е. наличие линейной связи).
mean_model <- mean(train$thalach) f <- sum((m_lm$fitted.values - mean_model)^2) / summary(m_lm)$sigma^2 f
## [1] 117.6397
f_tabl <- qf(0.95, 1, nrow(train)-2) f_tabl
## [1] 3.886774
F расч | > F табл => отвергаем H0 в пользу наличия линейной связи с вероятностью 0.95.
Находим коэффициенты детерминации:
#коэффициент множественной корреляции summary(m_lm)$r.squared
## [1] 0.3646163
#квадрат коэфициента скоректированной корреляции (доля дисперсии зависимой переменной price) summary(m_lm)$adj.r.squared
## [1] 0.355318
Значения коэффициентов детерминации довольно низки (ниже 0.5), подгонку построенной модели множественной регресии можно считать некачественной
Находим остаточную дисперсию:
summary(m_lm)$sigma^2
## [1] 364.1981
Средняя ошибка аппроксимации
A <- sum(abs((train$thalach - m_lm$fitted.values)/train$thalach)) / nrow(train) * 100 str(A)
## num 11.9
А > 5-7% => качество модели плохое
Построим прогноз по полученной модели на тестовой выборке
test$predict <- predict(m_lm, test) ggplot(test, aes(x = thalach, y = predict)) + geom_point() + geom_abline() + scale_x_continuous(limits = c(75, 200), expand = c(0, 0)) + scale_y_continuous(limits = c(75, 200), expand = c(0, 0))
## Warning: Removed 1 rows containing missing values (geom_abline).
Удалив все незначительные факторы, качество модели осталось плохим, даже понемного ухудшалось.