library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.0     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.2     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(haven)
library(dplyr)

df <- read_sav("~/files/База_КлимРиск_2023.sav")

Создаем бинарную зависимую переменную (V19 - оценка опасности). Предположим, что 1 - “опасно”, 0 - “не опасно”

df <- df %>%
  mutate( 
    danger = case_when(
      V19 == 1 ~ 1,  # считает опасным
      V19 == 2 ~ 0,  # не считает опасным
      TRUE ~ NA_real_
    ),
    sex = factor(V1, levels = c(1, 2), labels = c("Мужской", "Женский")),
    type = factor(type, levels = c(1, 2, 3, 4, 5, 6, 7), 
                  labels = c("тип 1", "тип 2", "тип 3", "тип 4", "тип 5", "тип 6", "тип 7")) 
  )

Проанализируем зависимую переменную

table(df$danger)
## 
##   0   1 
## 396 476
prop.table(table(df$danger))
## 
##         0         1 
## 0.4541284 0.5458716

Переменная danger принимает значения 0 и 1, можем приступить к построению логистической регрессии

model_climate <- glm(danger ~ sex + age + type, 
                     data = df, 
                     family = binomial,
                     na.action = na.omit)

summary(model_climate)
## 
## Call:
## glm(formula = danger ~ sex + age + type, family = binomial, data = df, 
##     na.action = na.omit)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.223006   0.323422  -3.781 0.000156 ***
## sexЖенский   0.746804   0.154576   4.831 1.36e-06 ***
## age          0.019764   0.005167   3.825 0.000131 ***
## typeтип 2   -0.594708   0.485797  -1.224 0.220880    
## typeтип 3    0.323568   0.525317   0.616 0.537929    
## typeтип 4   -0.271920   0.263149  -1.033 0.301449    
## typeтип 5   -0.305343   0.294134  -1.038 0.299219    
## typeтип 6   -0.069167   0.409702  -0.169 0.865936    
## typeтип 7    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
##   (78 пропущенных наблюдений удалены)
## AIC: 1107.7
## 
## Number of Fisher Scoring iterations: 4

Из полученных данных отметим, что все типы поселений (тип 2–7) — не оказывают статистически значимого влияния. Женщины значимо чаще считают проживание вблизи ледников опасным (0,747), также С возрастом оценка опасности значимо увеличивается (0,02). Пол и возраст являются значимыми факторами, влияющими на восприятие опасности проживания вблизи ледников, тогда как тип поселения значимой роли не играет. Женщины и люди старшего возраста склонны оценивать риски выше.

exp(cbind(OR = coef(model_climate), confint(model_climate)))
## Выполняю профилирование...
##                    OR     2.5 %   97.5 %
## (Intercept) 0.2943441 0.1553194 0.552889
## sexЖенский  2.1102447 1.5605409 2.861582
## age         1.0199606 1.0097583 1.030437
## typeтип 2   0.5517235 0.2067296 1.413412
## typeтип 3   1.3820495 0.5026969 4.049814
## typeтип 4   0.7619153 0.4537887 1.275211
## typeтип 5   0.7368704 0.4127697 1.309958
## typeтип 6   0.9331710 0.4170207 2.093803
## typeтип 7   1.3981061 0.8713250 2.237323

Данная таблица поддвеждает заявления, сделанные ранее: пол — самый сильный предиктор: женщины значительно чаще воспринимают ледники как опасные. Женщины: OR = 2.11 → вероятность выше на 111%. озраст: OR = 1.02 → с каждые 10 лет шансы растут на 22%.

pred_probs <- predict(model_climate, type = "response")

pred_class <- ifelse(pred_probs > 0.5, 1, 0)

confusion <- table(Предсказано = pred_class[!is.na(df$danger)], 
                   Фактически = df$danger[!is.na(df$danger)])
print(confusion)
##            Фактически
## Предсказано   0   1
##           0 146 155
##           1 218 282

Модель правильно предсказала 146 человек, которые НЕ считают ледники опасными. Модель ошиблась: предсказала “не опасно”, но на самом деле 155 человек считают опасно. Модель лучше предсказывает тех, кто считает ледники опасными (64.5%), чем тех, кто не считает (40.1%). Общая точность невысокая (53.4%), что говорит о необходимости добавления других предикторов.

# McFadden's R²
with(summary(model_climate), 1 - deviance/null.deviance)
## [1] 0.05342602
library(pscl)
## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.
pR2(model_climate)
## fitting null model for pseudo-r2
##           llh       llhNull            G2      McFadden          r2ML 
## -544.83170413 -575.58280042   61.50219257    0.05342602    0.07100816 
##          r2CU 
##    0.09492061
# Nagelkerke R²

library(fmsb)
NagelkerkeR2(model_climate)
## $N
## [1] 835
## 
## $R2
## [1] 0.09492061

McFadden R² = 0.053 — модель объясняет около 5.3% вариации зависимой переменной. В социальных науках значения 0.2–0.4 считаются очень хорошими, 0.05 — скромный результат. Nagelkerke R² = 0.095 — скорректированный псевдо R² показывает, что модель объясняет менее 10% дисперсии.Значения согласованы — все три показателя указывают на невысокую объясняющую способность модели. Модель с предикторами (пол, возраст, тип поселения) объясняет лишь около 9.5% вариации в восприятии опасности ледников. Это означает, что основные факторы, влияющие на мнение респондентов, не вошли в модель.