Логит-регрессия в R

Общее описание логичестической регрессии

Логистическая регрессия - регрессия, в которой объясняемая переменная - факторная, а объясняющие переменные могут быть и количественными, и факторными.

Возможные примеры бинарных факторных признаков:

мужчина или женщина
наличие или отсутствие заболевания
обнаружен или нет биологический вид в пробе
голосовал или нет избиратель
произошел ли дефолт по кредиту или нет

Бинарная логистическая модель - только два признака, многомерная регрессия - несколько признаков. Логистическая регрессия - один из видов обобщенных линейных моделей (GLM). Допущение - ошибки (и Y) имеют биномиальное распределение.

\[ log_e\frac{\pi}{1-\pi} = \beta_0 + \sum_{j=1}^{p} \beta_j X_j \]

где \(\pi = \mu_Y\) - условное мат ожидание Y (то есть вероятность того, что Y = 1 при имеющихся значениях X), \(\pi/(1-\pi)\) - шансы того, что Y = 1, \(log(\pi/(1-\pi))\) - логарифм шансов (или логит).
\(log(\pi/(1-\pi))\) - связующая функция (link function).

Почему нельзя использовать МНК для факторных Y?

В принципе это можно делать, но результаты будут плохими

Регрессия может “предсказывать” значения, выходящие за пределы возможных значений?
Ошибки обязательно не имеют нормального распределения
Гетерокседаксичность

Каковы общие принципы логистической регрессии?

В логистической регрессии, вместо того, чтобы предсказывать значения \(Y\) в зависимости от значений \(X_i\), предсказываются вероятность Y при реализации тех или значений \(X_i\). В простейшем виде (одна объясняющая переменная) логистическая регрессия имеет вид:

\[ P(Y) = \frac{1}{1+e^{-(b_0 + b_1* X_1)}} \]

Оценка логистической модель - статистика log-likehood

P(Y) может принимать значения между 0 (событие не произошло) и 1 (событие произошло). Общий подход - сравнивать предсказанные моделью и фактические значения сохранятся, но несколько модифицируется. Log-likehood - сумма верятностей, связанных с предсказанными и фактическими значениями Y. Статистика аналогична RSS (residiual sum squares) в линейной регрессии - то есть, показывают сколько сколько не-объясненной информации моделью осталось. Большие значения log-likehood указывают на плохую предсказательную способность модель. Чем меньше log-likehood - тем лучше.

\[ log-likehood = \sum_{i=1}^{N}[Y_i ln(P(Y_i)) + (1-Y_i) ln(1-P(Y_i))] \]

Девианс определяется как и имеет распределение хи-квадрат. \[ deviance = -2 * log-likehood \]

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

Интепретация логистической регрессии - шансы (odds)

Шансы события - отношение вероятности выполнения события к вероятности его не выполнения. К примеру, шансы забеременеть - вероятность забеременеть делить на вероятность не забеременеть.

\[ odds = \frac{P(event)}{P(no event)} \]

Проблемы, возникающие в расчетах логистической регрессии

полное разделение (complete separation)
неполная информация по объясняющим факторам

Пример расчета логистической регрессиии

Данные по некоторой группе пациентов:

Cured - вылечился/нет
Intervention - лечение/нет
Duration - количество дней до того, как началось лечение

Загрузим и немного трансформируем данные

eelData <- read.delim('datasets/eel.dat', header = TRUE)
eelData$Cured <- relevel(eelData$Cured, "Not Cured")
eelData$Intervention <- relevel(eelData$Intervention, "No Treatment")
head(eelData)

##       Cured Intervention Duration
## 1 Not Cured No Treatment        7
## 2 Not Cured No Treatment        7
## 3 Not Cured No Treatment        6
## 4     Cured No Treatment        8
## 5     Cured Intervention        7
## 6     Cured No Treatment        6

Проведение расчетов

eelModel1 <- glm( Cured ~ Intervention, data = eelData, family = binomial())
summary(eelModel1)

## 
## Call:
## glm(formula = Cured ~ Intervention, family = binomial(), data = eelData)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5940  -1.0579   0.8118   0.8118   1.3018  
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)   
## (Intercept)               -0.2877     0.2700  -1.065  0.28671   
## InterventionIntervention   1.2287     0.3998   3.074  0.00212 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 154.08  on 112  degrees of freedom
## Residual deviance: 144.16  on 111  degrees of freedom
## AIC: 148.16
## 
## Number of Fisher Scoring iterations: 4

exp(eelModel1$coefficients)

##              (Intercept) InterventionIntervention 
##                 0.750000                 3.416667

exp(confint(eelModel1))

## Waiting for profiling to be done...

##                              2.5 %   97.5 %
## (Intercept)              0.4374531 1.268674
## InterventionIntervention 1.5820127 7.625545

Интерпретация полученных коэфициентов.

Шансы выздороветь увеличились в 3,42 раза после использования лечения. Проверить с помощью расчетов!
доверительные интервалы для коэффициентов в логистической регрессии не должны пересекать 0 (или 1 для логитов). Пересечение означает “статистическую незначимость” значений коэффициентов. Значение логит-коэффициента больше 1 означает, что увеличение шансов выздороветь в результате изменения X, значения меньше 1 указывает на снижение шансов выздороветь.

Модель 2 = Модель 1 + продолжительность

eelModel2 <- glm( Cured ~ Intervention + Duration, data = eelData, family = binomial())
eelModel_ols <- lm( Cured ~ Intervention, data = eelData)

## Warning in model.response(mf, "numeric"): using type = "numeric" with a
## factor response will be ignored

## Warning in Ops.factor(y, z$residuals): '-' not meaningful for factors

summary(eelModel2)

## 
## Call:
## glm(formula = Cured ~ Intervention + Duration, family = binomial(), 
##     data = eelData)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6025  -1.0572   0.8107   0.8161   1.3095  
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)   
## (Intercept)              -0.234660   1.220563  -0.192  0.84754   
## InterventionIntervention  1.233532   0.414565   2.975  0.00293 **
## Duration                 -0.007835   0.175913  -0.045  0.96447   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 154.08  on 112  degrees of freedom
## Residual deviance: 144.16  on 110  degrees of freedom
## AIC: 150.16
## 
## Number of Fisher Scoring iterations: 4

anova(eelModel1, eelModel2)

## Analysis of Deviance Table
## 
## Model 1: Cured ~ Intervention
## Model 2: Cured ~ Intervention + Duration
##   Resid. Df Resid. Dev Df  Deviance
## 1       111     144.16             
## 2       110     144.16  1 0.0019835

На основе сравнения модели 1 и модели 2 можно сделать, вывод, что использование модели 1 предпочтительнее. Модель 2 не улучшает существенным расчетом “качество” модели

Графическое отображение результатов логистической регрессии

##    bodysize survive
## 1  25.19093       0
## 2  26.38245       0
## 3  26.57807       0
## 4  27.90968       0
## 5  28.33787       0
## 6  28.59298       1
## 7  28.93584       0
## 8  29.06218       1
## 9  30.29331       0
## 10 30.52695       0
## 11 30.61701       1
## 12 30.77641       1
## 13 30.92684       0
## 14 31.05090       1
## 15 31.86503       1
## 16 32.57736       1
## 17 32.59838       0
## 18 32.65738       1
## 19 33.78504       1
## 20 33.84042       1

Модель измен - пример логистической регрессии с несколькими факторами

загрузить данные

data(Affairs, package = "AER")
summary(Affairs)

##     affairs          gender         age         yearsmarried    children 
##  Min.   : 0.000   female:315   Min.   :17.50   Min.   : 0.125   no :171  
##  1st Qu.: 0.000   male  :286   1st Qu.:27.00   1st Qu.: 4.000   yes:430  
##  Median : 0.000                Median :32.00   Median : 7.000            
##  Mean   : 1.456                Mean   :32.49   Mean   : 8.178            
##  3rd Qu.: 0.000                3rd Qu.:37.00   3rd Qu.:15.000            
##  Max.   :12.000                Max.   :57.00   Max.   :15.000            
##  religiousness     education       occupation        rating     
##  Min.   :1.000   Min.   : 9.00   Min.   :1.000   Min.   :1.000  
##  1st Qu.:2.000   1st Qu.:14.00   1st Qu.:3.000   1st Qu.:3.000  
##  Median :3.000   Median :16.00   Median :5.000   Median :4.000  
##  Mean   :3.116   Mean   :16.17   Mean   :4.195   Mean   :3.932  
##  3rd Qu.:4.000   3rd Qu.:18.00   3rd Qu.:6.000   3rd Qu.:5.000  
##  Max.   :5.000   Max.   :20.00   Max.   :7.000   Max.   :5.000

Affairs$ynaffair[Affairs$affairs >0] <- 1
Affairs$ynaffair[Affairs$affairs == 0] <- 0
Affairs$ynaffair <- factor(Affairs$ynaffair, levels = c(0,1), labels = c("No", "Yes"))

table(Affairs$ynaffair)

## 
##  No Yes 
## 451 150

посчитать модель

fit.full <- glm(ynaffair ~ gender + age + yearsmarried + children + religiousness + education + occupation + rating, data = Affairs, family = binomial())
summary(fit.full)

## 
## Call:
## glm(formula = ynaffair ~ gender + age + yearsmarried + children + 
##     religiousness + education + occupation + rating, family = binomial(), 
##     data = Affairs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5713  -0.7499  -0.5690  -0.2539   2.5191  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    1.37726    0.88776   1.551 0.120807    
## gendermale     0.28029    0.23909   1.172 0.241083    
## age           -0.04426    0.01825  -2.425 0.015301 *  
## yearsmarried   0.09477    0.03221   2.942 0.003262 ** 
## childrenyes    0.39767    0.29151   1.364 0.172508    
## religiousness -0.32472    0.08975  -3.618 0.000297 ***
## education      0.02105    0.05051   0.417 0.676851    
## occupation     0.03092    0.07178   0.431 0.666630    
## rating        -0.46845    0.09091  -5.153 2.56e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 675.38  on 600  degrees of freedom
## Residual deviance: 609.51  on 592  degrees of freedom
## AIC: 627.51
## 
## Number of Fisher Scoring iterations: 4

exp(fit.full$coefficients)

##   (Intercept)    gendermale           age  yearsmarried   childrenyes 
##     3.9640180     1.3235091     0.9567099     1.0994093     1.4883560 
## religiousness     education    occupation        rating 
##     0.7227292     1.0212740     1.0314027     0.6259691

уберем незначимые переменные (пол, дети, образование, работа)

fit.reduced <- glm(ynaffair ~ age + yearsmarried + religiousness + rating, data = Affairs, family = binomial())
summary(fit.reduced)

## 
## Call:
## glm(formula = ynaffair ~ age + yearsmarried + religiousness + 
##     rating, family = binomial(), data = Affairs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6278  -0.7550  -0.5701  -0.2624   2.3998  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    1.93083    0.61032   3.164 0.001558 ** 
## age           -0.03527    0.01736  -2.032 0.042127 *  
## yearsmarried   0.10062    0.02921   3.445 0.000571 ***
## religiousness -0.32902    0.08945  -3.678 0.000235 ***
## rating        -0.46136    0.08884  -5.193 2.06e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 675.38  on 600  degrees of freedom
## Residual deviance: 615.36  on 596  degrees of freedom
## AIC: 625.36
## 
## Number of Fisher Scoring iterations: 4

так модели вложенные, их можно сравнить с помощью anova

anova(fit.reduced, fit.full, test = 'Chisq')

## Analysis of Deviance Table
## 
## Model 1: ynaffair ~ age + yearsmarried + religiousness + rating
## Model 2: ynaffair ~ gender + age + yearsmarried + children + religiousness + 
##     education + occupation + rating
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       596     615.36                     
## 2       592     609.51  4   5.8474   0.2108

Так значение хи-квадрат незначимо (p = 0,21), можно сделать вывод, что сокращенная модель работает не хуже полной. Значит, можно пользоваться “сокращенной” моделью.

exp(coef(fit.reduced))

##   (Intercept)           age  yearsmarried religiousness        rating 
##     6.8952321     0.9653437     1.1058594     0.7196258     0.6304248

Стандартные проблемы в логит-моделях

Overdispersion - ситуация, когда наблюдаемый разброс наблюдаемой величины больше того, который можно было бы ожидать из биномиального распределения. Это приводит к завышению стандартных ошибок и оценок статистической значимости коээффициентов (сами значения коээфициентов не изменяются при оverdispersion.

Критерий: Residual deviance / Residual df (должно быть не сильно больше 1).

При overdispersion необходимо задавать другое распределение ошибок (quasibinomial), а не просто binomial.