Отчет по домашней работе №3
# install.packages("pacman")
pacman::p_load(ISwR, ggplot2, dplyr, nortest, car, tidyverse, scales, emmeans, multcomp, multcompView, gtsummary, knitr)- В пакете ISwR содержится набор данных juul, который содержит концентрации инсулин-подобного фактора роста (igf1) и описательные данные пациентов (sex – пол, age – возраст). Загрузите этот набор данных.
data(juul)kable(head(juul))| age | menarche | sex | igf1 | tanner | testvol |
|---|---|---|---|---|---|
| NA | NA | NA | 90 | NA | NA |
| NA | NA | NA | 88 | NA | NA |
| NA | NA | NA | 164 | NA | NA |
| NA | NA | NA | 166 | NA | NA |
| NA | NA | NA | 131 | NA | NA |
| 0.17 | NA | 1 | 101 | 1 | NA |
str(juul)## 'data.frame': 1339 obs. of 6 variables:
## $ age : num NA NA NA NA NA 0.17 0.17 0.17 0.17 0.17 ...
## $ menarche: int NA NA NA NA NA NA NA NA NA NA ...
## $ sex : num NA NA NA NA NA 1 1 1 1 1 ...
## $ igf1 : num 90 88 164 166 131 101 97 106 111 79 ...
## $ tanner : int NA NA NA NA NA 1 1 1 1 1 ...
## $ testvol : int NA NA NA NA NA NA NA NA NA NA ...kable(summary(juul))| age | menarche | sex | igf1 | tanner | testvol | |
|---|---|---|---|---|---|---|
| Min. : 0.170 | Min. :1.000 | Min. :1.000 | Min. : 25.0 | Min. :1.00 | Min. : 1.000 | |
| 1st Qu.: 9.053 | 1st Qu.:1.000 | 1st Qu.:1.000 | 1st Qu.:202.2 | 1st Qu.:1.00 | 1st Qu.: 1.000 | |
| Median :12.560 | Median :1.000 | Median :2.000 | Median :313.5 | Median :2.00 | Median : 3.000 | |
| Mean :15.095 | Mean :1.476 | Mean :1.534 | Mean :340.2 | Mean :2.64 | Mean : 7.896 | |
| 3rd Qu.:16.855 | 3rd Qu.:2.000 | 3rd Qu.:2.000 | 3rd Qu.:462.8 | 3rd Qu.:5.00 | 3rd Qu.:15.000 | |
| Max. :83.000 | Max. :2.000 | Max. :2.000 | Max. :915.0 | Max. :5.00 | Max. :30.000 | |
| NA’s :5 | NA’s :635 | NA’s :5 | NA’s :321 | NA’s :240 | NA’s :859 |
- Создайте во фрейме новую переменную, содержащую квадратный корень из концентрации инсулин-подобного фактора роста.
juul$sqrt_igf1 <- sqrt(juul$igf1)
kable(tail(juul))| age | menarche | sex | igf1 | tanner | testvol | sqrt_igf1 | |
|---|---|---|---|---|---|---|---|
| 1334 | 58.95 | 2 | 2 | 218 | 5 | NA | 14.76482 |
| 1335 | 60.99 | 2 | 2 | 226 | 5 | NA | 15.03330 |
| 1336 | 62.73 | 2 | 2 | NA | NA | NA | NA |
| 1337 | 65.00 | 2 | 2 | 106 | NA | NA | 10.29563 |
| 1338 | 67.88 | 2 | 2 | 217 | NA | NA | 14.73092 |
| 1339 | 75.12 | 2 | 2 | 135 | NA | NA | 11.61895 |
- Постройте регрессионную модель зависимости квадратного корня концентрации инсулин-подобного фактора роста от возраста для пациентов старше 25 лет.
juul_over_25 <- subset(juul, age > 25)
fit <- lm(sqrt_igf1 ~ age, data = juul_over_25)
summary(fit)##
## Call:
## lm(formula = sqrt_igf1 ~ age, data = juul_over_25)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.8642 -1.1661 0.1018 0.9450 4.1136
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.71025 0.49462 37.828 <2e-16 ***
## age -0.10533 0.01072 -9.829 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.741 on 120 degrees of freedom
## (4 пропущенных наблюдений удалены)
## Multiple R-squared: 0.446, Adjusted R-squared: 0.4414
## F-statistic: 96.6 on 1 and 120 DF, p-value: < 2.2e-16- Проанализируйте характер распределения остатков регрессионной модели. Постройте гистограмму их распределения, примените три критерия согласия. В комментарии приведите p-значения и сделайте итоговый вывод.
residuals <- residuals(fit)
# 1. Критерий Шапиро-Уилка
shapiro_test <- shapiro.test(residuals)
# 2. Критерий Лиллиефорса
lillie <- lillie.test(residuals)
# 3. Критерий Андерсона-Дарлинга
ad_test <- ad.test(residuals)p_values <- data.frame(
Test = c("Шапиро-Уилка", "Лиллиефорса", "Андерсона-Дарлинга"),
`p-value` = c(shapiro_test$p.value, lillie$p.value, ad_test$p.value)
)
kable(p_values, col.names = c("Критерий", "p-value"))| Критерий | p-value |
|---|---|
| Шапиро-Уилка | 0.5945936 |
| Лиллиефорса | 0.7321537 |
| Андерсона-Дарлинга | 0.6734449 |
Поскольку все три теста показывают p > 0.05 (0.59, 0.73 и 0.67 соответственно), распределение остатков можно считать нормальным.
# Гистограмма остатков
hist(residuals, freq = FALSE)
Наблюдается распределение, близкое к нормальному.
- Проведите анализ гомоскедастичности регрессионной модели. В комментарии приведите название процедуры, p-значение и сделайте вывод.
Тест Бройша-Пагана:
ncvRes <- ncvTest(fit)
ncvRes## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.8506106, Df = 1, p = 0.35638bp_p_value <- ncvRes$p
cat("p-значение теста Бройша-Пагана:", bp_p_value, "\n")## p-значение теста Бройша-Пагана: 0.3563796if (bp_p_value > 0.05) {
cat("Гомоскедастичность подтверждена (p > 0.05). Модель удовлетворяет условию равенства дисперсий остатков.\n")
} else {
cat("Гетероскедастичность обнаружена (p < 0.05). Модель не удовлетворяет условию равенства дисперсий остатков.\n")
}## Гомоскедастичность подтверждена (p > 0.05). Модель удовлетворяет условию равенства дисперсий остатков.- Проведите анализ статистической значимости коэффициентов регрессионной модели. В комментарии приведите название процедуры, p-значение и сделайте вывод.
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.71025 0.49462 37.828 <2e-16 ***
## age -0.10533 0.01072 -9.829 <2e-16 ***Для каждого коэффициента модели проводится t-тест на статистическую значимость. Интерпретация:
p-value для обоих коэффициентов (интерсепта и переменной
age) < 2e-16, что намного меньше уровня значимости 0.05:
оба коэффициента (интерсепт и возраст) статистически
значимы. Это означает, что возраст оказывает значительное
влияние на уровень IGF-1 (с учетом проведенной трансформации), и
модель линейной регрессии значима.
- Как изменится ожидаемое значение переменной отклика в регрессионной модели при увеличении возраста на 5 месяцев?
# Извлекаем коэффициенты из модели
coefficients_model <- coef(fit)
# Коэффициент при переменной age
beta_1 <- coefficients_model['age']
# Рассчитываем изменение ожидаемого значения sqrt_igf1 при увеличении возраста на 5 месяцев
change_in_response <- beta_1 * (5 / 12)
# Выводим результат
cat("Изменение ожидаемого значения sqrt_igf1 при увеличении возраста на 5 месяцев:", change_in_response, "\n")## Изменение ожидаемого значения sqrt_igf1 при увеличении возраста на 5 месяцев: -0.04388695- Каково ожидаемое значение переменной отклика для следующих значений возраста: 41.77, 41.88, 43.84, 46.08, 54.91?
ages <- c(41.77, 41.88, 43.84, 46.08, 54.91)
new_data <- data.frame(age = ages)
prediction <- predict(fit, new_data)
kable(data.frame(age=ages, predicted=prediction))| age | predicted |
|---|---|
| 41.77 | 14.31067 |
| 41.88 | 14.29909 |
| 43.84 | 14.09264 |
| 46.08 | 13.85670 |
| 54.91 | 12.92665 |
- Рассчитайте и приведите в комментарии 95 %-ный доверительный интервал для коэффициента наклона.
conf_intervals <- confint(fit)
kable(conf_intervals['age', ], col.names = c('age'), label = '95%-ный доверительный интервал для коэффициента наклона (age)')| age | |
|---|---|
| 2.5 % | -0.1265467 |
| 97.5 % | -0.0841107 |
Интерпретация:
Если доверительный интервал не включает 0 (например, если оба конца интервала строго положительные или строго отрицательные), это указывает на статистически значимый эффект возраста на переменную отклика.
Если доверительный интервал включает 0, это означает, что влияние возраста на переменную отклика может быть статистически незначимым на уровне 95%-ного доверия.
В данном случае оба конца интервала строго отрицательные, значит есть статистически значимый эффект возраста на переменную отклика.
- Постройте диаграмму рассеяния для данных, на основе которых построена анализируемая регрессионная модель. Используйте залитые полупрозрачные круги зеленого цвета. Подпишите оси. Добавьте на график линию регрессии.
# Построение диаграммы рассеяния
plot(juul_over_25$age, juul_over_25$sqrt_igf1,
col = rgb(0, 1, 0, 0.5),
pch = 16,
xlab = "Возраст",
ylab = "IGF-1",
main = "Диаграмма рассеяния с линией регрессии")
# Добавление линии регрессии
abline(fit, col = "blue", lwd = 2) # Линия регрессии синего цвета
- Добавьте на график доверительные интервалы для положения линии регресии и для прогнозируемых данных. Границы интервалов должны отличаться цветом от линии регрессии и друг от друга.
age_range <- data.frame(age = seq(min(juul_over_25$age), max(juul_over_25$age), length.out = 100))
confidence_intervals <- predict(fit, newdata = age_range, interval = "confidence")
prediction_intervals <- predict(fit, newdata = age_range, interval = "prediction")
plot(juul_over_25$age, juul_over_25$sqrt_igf1,
col = rgb(0, 1, 0, 0.5), # Полупрозрачный зелёный цвет
pch = 16, # Залитые круги
xlab = "Возраст", # Подпись оси X
ylab = "IGF-1", # Подпись оси Y
main = "Диаграмма с линией регрессии и доверительными интервалами")
abline(fit, col = "blue", lwd = 2) # Линия регрессии
lines(age_range$age, confidence_intervals[, "lwr"], col = "red", lty = 2, lwd = 2) # Нижняя граница доверительного интервала
lines(age_range$age, confidence_intervals[, "upr"], col = "red", lty = 2, lwd = 2) # Верхняя граница доверительного интервала
lines(age_range$age, prediction_intervals[, "lwr"], col = "purple", lty = 3, lwd = 2) # Нижняя граница прогноза
lines(age_range$age, prediction_intervals[, "upr"], col = "purple", lty = 3, lwd = 2) # Верхняя граница прогноза