Загружаем необходимые пакеты, говорим knitr’у, что надо кэшировать куски для ускорения работы, и устанавливаем чёрно-белую тему оформления графиков:

library("knitr") # создание отчётов
library("reshape2") # структуризация данных
library("ggplot2") # графики
library("vcd") # графики для качественных переменных

# opts_chunk$set(cache=TRUE)

# theme_set(theme_gray())
theme_set(theme_bw())

Загружаем данные по стоимости квартир в Москве и проверяем, что они загрузились:

filename <- "../datasets/flats_moscow.txt"
flats <- read.table(filename, header = TRUE)

str(flats)
## 'data.frame':    2040 obs. of  11 variables:
##  $ n       : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ price   : int  81 75 128 95 330 137 98 88 225 140 ...
##  $ totsp   : int  58 44 70 61 104 76 59 55 80 86 ...
##  $ livesp  : int  40 28 42 37 60 50 39 36 56 51 ...
##  $ kitsp   : num  6 6 6 6 11 9 6 6 9 10 ...
##  $ dist    : num  12.5 13.5 14.5 13.5 10.5 11 7.5 9 9 12.7 ...
##  $ metrdist: int  7 7 3 7 7 7 10 5 5 10 ...
##  $ walk    : int  1 1 1 1 0 1 0 1 1 1 ...
##  $ brick   : int  1 0 1 0 1 1 0 1 1 0 ...
##  $ floor   : int  1 1 1 1 1 1 1 0 1 1 ...
##  $ code    : int  3 6 3 1 3 8 8 4 3 5 ...

Оцениваем logit модель, полином 5-ой степени, \(P(y_i=1)=F(\beta_1+\beta_2 x_i+\ldots+\beta_6 x_i^5)\)

m1 <- glm(brick ~ poly(price, 5), data = flats, family = binomial)

summary(m1)
## 
## Call:
## glm(formula = brick ~ poly(price, 5), family = binomial, data = flats)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9784  -0.7769  -0.6932   1.0781   1.7766  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -0.78768    0.05063 -15.558  < 2e-16 ***
## poly(price, 5)1  25.97300    3.11629   8.335  < 2e-16 ***
## poly(price, 5)2  -9.29670    5.66507  -1.641   0.1008    
## poly(price, 5)3   5.71829    5.71520   1.001   0.3170    
## poly(price, 5)4   6.05547    3.21044   1.886   0.0593 .  
## poly(price, 5)5 -13.74891    2.67952  -5.131 2.88e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2566.9  on 2039  degrees of freedom
## Residual deviance: 2354.5  on 2034  degrees of freedom
## AIC: 2366.5
## 
## Number of Fisher Scoring iterations: 6

Строим графики:

ggplot(data = flats, aes(x = price, y = brick)) +
  geom_point(position = position_jitter(height = 0.05), alpha = 0.2) +
  stat_smooth(method.args = list(family = "binomial"), 
              method = "glm", 
              formula = y ~ poly(x, 4)) +
  labs(x = "Цена квартиры, $1000", 
       y = "Кирпичный дом или нет", 
       title = "График logit модели")

Две гистограммки:

ggplot(data = flats, aes(x = price, fill = factor(brick))) + geom_density(alpha = 0.5) + 
  scale_fill_discrete(name = "Материал", 
                      breaks = c("1", "0"),  
                      labels = c("Кирпич", "Другой")) +
  labs(x = "Цена квартиры, $1000", 
       y = "Плотность", 
       title = "Цены квартир в кирпичных и не кирпичных домах")

Виолончельки:

ggplot(data = flats, aes(y = price, x = factor(brick))) + geom_violin() + 
  scale_x_discrete(labels = c("Другой", "Кирпич")) + 
  labs(x = "Материал", 
       y = "Цена квартиры, $1000", 
       title = "Цены квартир в кирпичных и не кирпичных домах")

Прогнозируем вероятности по нашей модели. Сравниваем их с настоящими \(y_i\)

flats$prob.brick <- fitted(m1)

head(flats[c("brick", "prob.brick")])
##   brick prob.brick
## 1     1  0.2225129
## 2     0  0.2444560
## 3     1  0.2960600
## 4     0  0.2065029
## 5     1  0.5002794
## 6     1  0.3450575

Создаём функцию, которая по порогу считает количество true positive, true negative и т.д.

tfnp <- function(cut = 0.5, y.true, prob) {
  y.pred <- (prob > cut)
  tp <- sum((y.true == 1) & (y.pred == 1))
  tn <- sum((y.true == 0) & (y.pred == 0))
  fp <- sum((y.true == 0) & (y.pred == 1))
  fn <- sum((y.true == 1) & (y.pred == 0))
  return(c(tp, tn, fp, fn))
}

tfnp(0.5, flats$brick, flats$prob.brick) # проверяем
## [1]  185 1275  106  474

Для разных порогов считаем специфичность и чувствительность

t <- data.frame(cuts = seq(0,1, len = 33), tp = 0, tn = 0, fp = 0, fn = 0)

for (i in 1:nrow(t)) {
  t[i, 2:5] <- tfnp(t$cuts[i], flats$brick, flats$prob.brick)
}

t$spec <- t$tn / (t$tn + t$fp)
t$sens <- t$tp / (t$tp + t$fn)
t$correct <- (t$tp + t$tn) / (t$tp + t$tn + t$fn + t$fp)

head(t)
##      cuts  tp tn   fp fn spec sens   correct
## 1 0.00000 659  0 1381  0    0    1 0.3230392
## 2 0.03125 659  0 1381  0    0    1 0.3230392
## 3 0.06250 659  0 1381  0    0    1 0.3230392
## 4 0.09375 659  0 1381  0    0    1 0.3230392
## 5 0.12500 659  0 1381  0    0    1 0.3230392
## 6 0.15625 659  0 1381  0    0    1 0.3230392

Зависимость специфичности от порога

matplot(t$cuts, cbind(t$spec, t$sens, t$correct), type = "l")

molten.t <- melt(t, id.vars = "cuts", measure.vars = c("spec", "sens", "correct"))

ggplot(data = molten.t, aes(x = cuts, y = value, color = variable)) + geom_line() +
  labs(title = "График зависимости специфичности от порога")

Кривая ROC

ggplot(data = t, aes(x = 1 - spec, y = sens)) + geom_line()

Красивый график для качественных переменных:

mosaic(data = flats, ~ brick + walk + floor, shade = TRUE)