Скажем knitr'у, что надо кэшировать куски для ускорения работы…
library(knitr)
opts_chunk$set(cache = TRUE)
Загрузим данные по стоимости квартир в Москве. Проверим, что они загрузились…
setwd("~/science/econometrix/em301/datasets/")
h <- read.table("flats_moscow.txt", header = TRUE)
str(h)
## '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 ...
Загрузим пакет ggplot2 и установим чёрно-белую тему для всех последующих графиков
library(ggplot2)
# theme_set(theme_gray())
theme_set(theme_bw())
Оценим logit модель, полином 5-ой степени, \( P(y_i=1)=F(\beta_1+\beta_2 x_i+\ldots+\beta_6 x_i^5) \)
model <- glm(brick ~ poly(price, 5), data = h, family = binomial)
summary(model)
##
## Call:
## glm(formula = brick ~ poly(price, 5), family = binomial, data = h)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.978 -0.777 -0.693 1.078 1.777
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.7877 0.0506 -15.56 < 2e-16 ***
## poly(price, 5)1 25.9730 3.1163 8.33 < 2e-16 ***
## poly(price, 5)2 -9.2967 5.6651 -1.64 0.101
## poly(price, 5)3 5.7183 5.7152 1.00 0.317
## poly(price, 5)4 6.0555 3.2104 1.89 0.059 .
## poly(price, 5)5 -13.7489 2.6795 -5.13 2.9e-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
##
## Number of Fisher Scoring iterations: 6
Строим графики…
ggplot(data = h, aes(x = price, y = brick)) + geom_point(position = position_jitter(height = 0.05),
alpha = 0.2) + stat_smooth(method = "glm", family = binomial, formula = y ~
poly(x, 4))
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
Две гистограммки
ggplot(data = h, aes(x = price, fill = factor(brick))) + geom_density(alpha = 0.5)
Виолончельки
ggplot(data = h, aes(y = price, x = factor(brick))) + geom_violin()
Прогнозируем вероятности по нашей модели. Сравниваем их с настоящими \( y_i \)
h$prob.brick <- fitted(model)
head(h[c("brick", "prob.brick")])
## brick prob.brick
## 1 1 0.2225
## 2 0 0.2445
## 3 1 0.2961
## 4 0 0.2065
## 5 1 0.5003
## 6 1 0.3451
Создаём функцию, которая по порогу считает количество 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, h$brick, h$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], h$brick, h$prob.brick)
}
t$spec <- t$tn/(t$tn + t$fp)
t$sens <- t$tp/(t$tp + t$fn)
head(t)
## cuts tp tn fp fn spec sens
## 1 0.00000 659 0 1381 0 0 1
## 2 0.03125 659 0 1381 0 0 1
## 3 0.06250 659 0 1381 0 0 1
## 4 0.09375 659 0 1381 0 0 1
## 5 0.12500 659 0 1381 0 0 1
## 6 0.15625 659 0 1381 0 0 1
Зависимость специфичности от порога
ggplot(data = t, aes(x = cuts, y = sens)) + geom_line()
ggplot(data = t, aes(x = 1 - spec, y = sens)) + geom_line()
Еще один график симпатяшка для качественных переменных
library(vcd)
## Loading required package: MASS
## Loading required package: grid
## Loading required package: colorspace
mosaic(data = h, ~brick + walk + floor, shade = TRUE)
Остальные заметочки, https://github.com/bdemeshev/em301/wiki/r_cycle
xxx