library(readr)
## Warning: пакет 'readr' был собран под R версии 4.5.2
library(haven)
## Warning: пакет 'haven' был собран под R версии 4.5.2
climat = read_sav(file.choose())
table(climat$V19)
## 
##   1   2 
## 476 396
unique(climat$V19)
## <labelled<double>[3]>: Проживание вблизи ледников, в зоне вечной мерзлоты в условиях изменения климата представляет какую-то особую опасность, по Вашему мнению, или нет?
## [1]  2 NA  1
## 
## Labels:
##  value                label
##      1     Да, представляет
##      2 Нет, не представляет
summary(climat$V19)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   1.000   1.000   1.000   1.454   2.000   2.000      41
climat$V19_binary <- ifelse(climat$V19 == 1, 1, 0)
table(climat$V19_binary, climat$V19)
##    
##       1   2
##   0   0 396
##   1 476   0
model_glm <- glm(V19_binary ~ V1 + age + type, 
                 data = climat, 
                 family = "binomial")

summary(model_glm)
## 
## Call:
## glm(formula = V19_binary ~ V1 + age + type, family = "binomial", 
##     data = climat)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.484501   0.392552  -6.329 2.47e-10 ***
## V1           0.803170   0.152273   5.275 1.33e-07 ***
## age          0.019434   0.005111   3.802 0.000143 ***
## type         0.092082   0.035109   2.623 0.008722 ** 
## ---
## 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: 1100.0  on 831  degrees of freedom
##   (78 пропущенных наблюдений удалены)
## AIC: 1108
## 
## Number of Fisher Scoring iterations: 4

V1=1 к V1=2 ( от мужчины к женщине) логарифм шансов V19_binary = 1 увеличивается на 0,803. У второй группы (женщин, V1=2) шансы ответить “да,представляет” в 2,23 раза выше, чем у первой группы (мужчин). При увеличении возраста на 1 год логарифм шансов увеличивается на 0,019. Все предикторы (V1, age, type) статистически значимы и влияют на V19_binary. V1 — самый сильный предиктор (наибольший z-score: 5.275).

exp(cbind(OR = coef(model_glm), confint(model_glm)))
## Выполняю профилирование...
##                     OR      2.5 %    97.5 %
## (Intercept) 0.08336714 0.03822542 0.1783554
## V1          2.23260658 1.65887460 3.0144790
## age         1.01962401 1.00953291 1.0299814
## type        1.09645438 1.02367750 1.1748684

отношения шансов (Odds Ratios, OR) с 95% доверительными интервалами для логистической регрессионной модели.

model_glm_pred = ifelse(predict(model_glm, type = "response") > 0.5, 1, 0)
head(model_glm_pred)
##  5  6  8  9 12 14 
##  1  1  1  1  1  1

Вектор предсказанных классов на основе модели логистической регрессии с порогом 0.5. Все указанные наблюдения предсказаны как класс 1 “да, представляет”.

climat_clean <- na.omit(climat[, c("V1", "age", "type", "V19_binary")])

model_glm <- glm(V19_binary ~ V1 + age + type, 
                 data = climat_clean, 
                 family = "binomial")

# Предсказания
model_glm_pred <- ifelse(predict(model_glm, type = "response") >= 0.5, 1, 0)

# Теперь длины совпадают
tab = table(predicted = model_glm_pred, actual = climat_clean$V19_binary)
tab
##          actual
## predicted   0   1
##         0 173 117
##         1 208 337

Модель завышает оценку обеспокоенности: 208 респондентов, которые реально не считают ледники опасными, были ошибочно отнесены к группе “обеспокоенных”. (True Positives = 337) Модель успешно идентифицирует респондентов, которые действительно считают проживание вблизи ледников опасным. False Negative (FN) = 117 — модель ошибочно предсказала «0», хотя на самом деле ответ был «1» (опасность есть). True Negative (TN) = 173 — модель верно предсказала «0» (опасности нет), и это соответствует реальным данным.

library(caret)
## Warning: пакет 'caret' был собран под R версии 4.5.3
## Загрузка требуемого пакета: ggplot2
## Warning: пакет 'ggplot2' был собран под R версии 4.5.2
## Загрузка требуемого пакета: lattice
confusionMatrix = confusionMatrix(tab, positive = "1")
c(confusionMatrix$overall["Accuracy"], 
  confusionMatrix$byClass["Sensitivity"], 
  confusionMatrix$byClass["Specificity"])
##    Accuracy Sensitivity Specificity 
##   0.6107784   0.7422907   0.4540682

модель правильно классифицирует 61,08 % всех наблюдений из всех людей, которые считают проживание вблизи ледников опасным (фактические «1»), модель правильно идентифицировала 74,23 % из всех людей, которые не считают проживание опасным (фактические «0»), модель правильно идентифицировала только 45,41 %.

library(pscl)
## Warning: пакет 'pscl' был собран под R версии 4.5.3
## 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_glm)['McFadden']
## fitting null model for pseudo-r2
## McFadden 
## 0.044471

0,0 – 0,1 Очень слабая объяснительная способность

library(fmsb)
## Warning: пакет 'fmsb' был собран под R версии 4.5.3
## Registered S3 methods overwritten by 'fmsb':
##   method    from
##   print.roc pROC
##   plot.roc  pROC
NagelkerkeR2(model_glm)
## $N
## [1] 835
## 
## $R2
## [1] 0.07949419

Модель, включающая пол, возраст и тип населенного пункта, объясняет около 8% дисперсии в восприятии опасности

P_hat <- predict(model_glm,
                 type = "response")

pred_class <- ifelse(P_hat >= 0.5, 1, 0)

y_true <- climat_clean$V19_binary

mean(pred_class == y_true)
## [1] 0.6107784

модель правильно предсказывает класс (0 или 1) для 61,1% респондентов.

AIC(model_glm)
## [1] 1107.972

значение 1108 само по себе не говорит “хорошо” или “плохо”. Оно имеет смысл только при сравнении с другими моделями на тех же данных.