Бинарная логистическая регрессия: оценка опасности проживания вблизи ледников
2026-04-12
1. Постановка задачи
Методом бинарной логистической регрессии анализируется связь между оценкой опасности проживания вблизи ледников (вопрос V19) и тремя предикторами: полом (V1), возрастом (age) и типом района проживания (type).
2. Загрузка и подготовка данных
if (!require("readxl")) install.packages("readxl", repos = "https://cran.r-project.org")
if (!require("knitr")) install.packages("knitr", repos = "https://cran.r-project.org")
if (!require("kableExtra")) install.packages("kableExtra", repos = "https://cran.r-project.org")
if (!require("ggplot2")) install.packages("ggplot2", repos = "https://cran.r-project.org")
if (!require("pROC")) install.packages("pROC", repos = "https://cran.r-project.org")
if (!require("ResourceSelection")) install.packages("ResourceSelection", repos = "https://cran.r-project.org")
library(readxl); library(knitr); library(kableExtra); library(ggplot2); library(pROC); library(ResourceSelection)data <- read_excel("База_КлимРиск_2023 (3).xlsx")
df <- data.frame(V19 = as.numeric(data$V19), sex = as.factor(data$V1), age = as.numeric(data$age), type = as.factor(data$type))
df <- df[complete.cases(df), ]
df$danger <- ifelse(df$V19 == 1, 1, 0)
n_total <- nrow(df); n_danger <- sum(df$danger == 1); n_safe <- sum(df$danger == 0); pct_danger <- round(n_danger / n_total * 100, 1)Комментарий: Загружено 835 наблюдений. Из них 454 (54.4%) считают проживание вблизи ледников опасным, 381 — не считают.
3. Описательная статистика
sex_table <- table(df$sex)
age_mean <- round(mean(df$age), 1); age_sd <- round(sd(df$age), 1)
danger_by_sex <- round(prop.table(table(df$sex, df$danger), margin = 1) * 100, 1)
kable(danger_by_sex, caption = "Таблица 1 — Доля оценки «опасно» по полу (%)") %>% kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| 0 | 1 |
|---|---|
| 58.5 | 41.5 |
| 39.3 | 60.7 |
Вывод: В выборке 275 мужчин и 560 женщин. Средний возраст — 43.9 лет (SD = 14.3). Среди мужчин «опасно» выбрали 41.5%, среди женщин — 60.7%.
4. Построение модели
model <- glm(danger ~ sex + age + type, data = df, family = binomial(link = "logit"))
model_summary <- summary(model)
model_summary##
## Call:
## glm(formula = danger ~ sex + age + type, family = binomial(link = "logit"),
## data = df)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.223006 0.323422 -3.781 0.000156 ***
## sex2 0.746804 0.154576 4.831 1.36e-06 ***
## age 0.019764 0.005167 3.825 0.000131 ***
## type2 -0.594708 0.485797 -1.224 0.220880
## type3 0.323568 0.525317 0.616 0.537929
## type4 -0.271920 0.263149 -1.033 0.301449
## type5 -0.305343 0.294134 -1.038 0.299219
## type6 -0.069167 0.409702 -0.169 0.865936
## type7 0.335119 0.240062 1.396 0.162724
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1151.2 on 834 degrees of freedom
## Residual deviance: 1089.7 on 826 degrees of freedom
## AIC: 1107.7
##
## Number of Fisher Scoring iterations: 4coef_table <- model_summary$coefficients
sig_vars <- rownames(coef_table)[coef_table[, 4] < 0.05 & rownames(coef_table) != "(Intercept)"]Вывод: Девианс снизился с 1151.2 (нулевая модель) до 1089.7 (полная), что указывает на улучшение. Значимые предикторы (p < 0.05): sex2, age.
5. Отношения шансов
ci <- confint(model)
or_table <- data.frame(Переменная = names(coef(model)), Коэффициент = round(coef(model), 4), OR = round(exp(coef(model)), 4), CI_lower = round(exp(ci)[, 1], 4), CI_upper = round(exp(ci)[, 2], 4), p_value = round(model_summary$coefficients[, 4], 4))
kable(or_table, caption = "Таблица 2 — Коэффициенты и отношения шансов (OR)", row.names = FALSE) %>% kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Переменная | Коэффициент | OR | CI_lower | CI_upper | p_value |
|---|---|---|---|---|---|
| (Intercept) | -1.2230 | 0.2943 | 0.1553 | 0.5529 | 0.0002 |
| sex2 | 0.7468 | 2.1102 | 1.5605 | 2.8616 | 0.0000 |
| age | 0.0198 | 1.0200 | 1.0098 | 1.0304 | 0.0001 |
| type2 | -0.5947 | 0.5517 | 0.2067 | 1.4134 | 0.2209 |
| type3 | 0.3236 | 1.3820 | 0.5027 | 4.0498 | 0.5379 |
| type4 | -0.2719 | 0.7619 | 0.4538 | 1.2752 | 0.3014 |
| type5 | -0.3053 | 0.7369 | 0.4128 | 1.3100 | 0.2992 |
| type6 | -0.0692 | 0.9332 | 0.4170 | 2.0938 | 0.8659 |
| type7 | 0.3351 | 1.3981 | 0.8713 | 2.2373 | 0.1627 |
or_sex <- round(exp(coef(model)["sex2"]), 2); p_sex <- round(coef_table["sex2", 4], 4)
or_age <- round(exp(coef(model)["age"]), 4); p_age <- round(coef_table["age", 4], 4)Вывод: OR для пола (жен. vs муж.) = 2.11 (p = 0): шансы женщин оценить проживание как опасное выше в 2.11 раза. OR для возраста = 1.02 (p = 10^{-4}): каждый дополнительный год увеличивает шансы на 2%.
6. Качество модели
hl_test <- hoslem.test(model$y, fitted(model), g = 10)
pseudo_r2 <- round(1 - (model$deviance / model$null.deviance), 4)
aic_val <- round(AIC(model), 2)
roc_obj <- roc(df$danger, fitted(model)); auc_val <- round(auc(roc_obj), 3)| Показатель | Значение | Интерпретация |
|---|---|---|
| Тест Хосмера-Лемешова (χ²) | 9.703 (p = 0.2865) | Модель адекватна |
| Псевдо-R² McFadden | 0.0534 | Слабая объяснительная сила |
| AIC | 1107.66 | Чем ниже, тем лучше |
| AUC | 0.651 | Слабое качество классификации |
7. ROC-кривая
plot(roc_obj, col = "#5C6B3C", lwd = 2, main = paste0("ROC-кривая (AUC = ", auc_val, ")"), print.auc = TRUE, auc.polygon = TRUE, auc.polygon.col = "#5C6B3C20")
Рис. 1 — ROC-кривая
Вывод: AUC = 0.651 — модель слабо различает группы, что типично для субъективных оценок в социологии.
8. Классификационная таблица
predicted <- ifelse(fitted(model) > 0.5, 1, 0)
conf_matrix <- table(Факт = df$danger, Прогноз = predicted)
kable(conf_matrix, caption = "Таблица 3 — Матрица классификации (порог 0.5)") %>% kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| 0 | 1 | |
|---|---|---|
| 0 | 191 | 190 |
| 1 | 118 | 336 |
accuracy <- round(sum(diag(conf_matrix)) / sum(conf_matrix) * 100, 1)
sensitivity <- round(conf_matrix[2,2] / sum(conf_matrix[2,]) * 100, 1)
specificity <- round(conf_matrix[1,1] / sum(conf_matrix[1,]) * 100, 1)Вывод: Точность — 63.1%, чувствительность — 74%, специфичность — 50.1%. Модель лучше предсказывает один класс.
9. Визуализация эффектов
df$predicted_prob <- fitted(model)
ggplot(df, aes(x = age, y = predicted_prob, color = sex)) + geom_smooth(method = "loess", se = TRUE) +
labs(title = "Предсказанная вероятность оценки «опасно»", x = "Возраст", y = "P(опасно)", color = "Пол") +
scale_color_manual(values = c("1" = "#5C6B3C", "2" = "#CC3333"), labels = c("1" = "Мужской", "2" = "Женский")) + theme_minimal()
Рис. 2 — Предсказанная вероятность по возрасту и полу
prob_male <- round(mean(df$predicted_prob[df$sex == "1"]), 3); prob_female <- round(mean(df$predicted_prob[df$sex == "2"]), 3)
prob_young <- round(mean(df$predicted_prob[df$age < 35]), 3); prob_old <- round(mean(df$predicted_prob[df$age >= 55]), 3)Вывод: Средняя P(опасно) для мужчин = 0.415, для женщин = 0.607 (разница 19.2 п.п.). Для молодых (< 35) = 0.469, для старших (≥ 55) = 0.63.
10. Итоговые выводы
- Модель построена на 835 наблюдениях. Значимые предикторы: sex2, age.
- Псевдо-R² = 0.0534 — слабая объяснительная сила.
- AUC = 0.651 — слабое качество классификации.
- Точность = 63.1%, что превышает базовую частоту наиболее частого класса.
- Тип района и пол могут объяснять часть различий в восприятии опасности проживания вблизи ледников.
Анализ выполнен в R с использованием пакетов pROC и ResourceSelection.