Введение

Шла 7-ая неделя, придумывать новые интересные штуки становилось все сложнее … Поэтому от мира чудесных, но абстрактных, моделей прогнозирования рисков будем плавно перетекать к сколько-нибудь полезным вещам. Лекция будет обо всем по чуть-чуть, постараюсь рассказать про интересные модели из мира анализа данных и эконометрики. А именно то, что в свое время стало открытием для меня.

В этом файле в основном картинки, основное постараюсь рассказать словами.

Проклятие размерности

Серьезно, это официальный термин. Один из самых крутых, что я знаю. Проклятие размерности – это проблема всех аналитиков данных. Хотя это не так плохо как быть проклятым на вечное хождение по морю на Летучем Голландце.

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

Метод главных компонент

Я думаю все помнят, что это такое. А тем кто не помнит, я сейчас расскажу у доски (не зря же Вы пришли на лекцию). Возьмем нашего старого знакомого – mtcars. И вспомним корреляционную матрицу:

library(corrplot)
corrplot(cor(mtcars))

Датасет сильно скоррелирован, а значит можно много выкинуть, без (сильной) потери качества. Выделим главные компоненты:

PCA <- princomp(scale(mtcars))
summary(PCA)
## Importance of components:
##                           Comp.1    Comp.2     Comp.3     Comp.4
## Standard deviation     2.5301952 1.6023860 0.77948532 0.51105040
## Proportion of Variance 0.6007637 0.2409516 0.05701793 0.02450886
## Cumulative Proportion  0.6007637 0.8417153 0.89873322 0.92324208
##                            Comp.5     Comp.6     Comp.7     Comp.8
## Standard deviation     0.46526149 0.45275130 0.36198764 0.34505183
## Proportion of Variance 0.02031374 0.01923601 0.01229654 0.01117286
## Cumulative Proportion  0.94355581 0.96279183 0.97508837 0.98626123
##                             Comp.9     Comp.10     Comp.11
## Standard deviation     0.273201294 0.224520229 0.146135274
## Proportion of Variance 0.007004241 0.004730495 0.002004037
## Cumulative Proportion  0.993265468 0.997995963 1.000000000
PCA$loadings
## 
## Loadings:
##      Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9
## mpg   0.363        -0.226        -0.103 -0.109  0.368  0.754  0.236
## cyl  -0.374        -0.175                0.169         0.231       
## disp -0.368                0.257 -0.394 -0.336  0.214         0.198
## hp   -0.330 -0.249  0.140        -0.540                0.222 -0.576
## drat  0.294 -0.275  0.161  0.855         0.244                     
## wt   -0.346  0.143  0.342  0.246        -0.465                0.359
## qsec  0.200  0.463  0.403         0.165 -0.330         0.232 -0.528
## vs    0.307  0.232  0.429 -0.215 -0.600  0.194 -0.266         0.359
## am    0.235 -0.429 -0.206               -0.571 -0.587              
## gear  0.207 -0.462  0.290 -0.265        -0.244  0.605 -0.336       
## carb -0.214 -0.414  0.529 -0.127  0.361  0.184 -0.175  0.396  0.171
##      Comp.10 Comp.11
## mpg   0.139  -0.125 
## cyl  -0.846  -0.141 
## disp          0.661 
## hp    0.248  -0.256 
## drat -0.101         
## wt           -0.567 
## qsec -0.271   0.181 
## vs   -0.159         
## am   -0.178         
## gear -0.214         
## carb          0.320 
## 
##                Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
## SS loadings     1.000  1.000  1.000  1.000  1.000  1.000  1.000  1.000
## Proportion Var  0.091  0.091  0.091  0.091  0.091  0.091  0.091  0.091
## Cumulative Var  0.091  0.182  0.273  0.364  0.455  0.545  0.636  0.727
##                Comp.9 Comp.10 Comp.11
## SS loadings     1.000   1.000   1.000
## Proportion Var  0.091   0.091   0.091
## Cumulative Var  0.818   0.909   1.000

И зачем? Во-первых, теперь можно использовать меньше переменных:

plot(cumsum(PCA$sdev)/sum(PCA$sdev), ylab = "Доля объясненной дисперсии")

plot(PCA)

Во-вторых, можно визуализировать.

plot(PCA$scores)

Или даже так:

biplot(PCA)

Это ли не наглядная визуализация?

Но это можно было и раньше. Но зато теперь в этих компонентах максимум дисперсии. А значит эта визуализация хороша. Я вот о чем:

plot(PCA$scores, col=factor(mtcars$cyl))

plot(PCA$scores, col=factor(mtcars$am))

plot(PCA$scores, col=factor(mtcars$vs))

plot(PCA$scores, col=factor(mtcars$carb))

LDA – Дискриминантный анализ

Я не буду притворяться, что хорошо понимаю, как работает LDA, потому что PCA всегда хватает. Но тем не менее:

library(MASS)
lda(cyl ~., data=mtcars)
## Call:
## lda(cyl ~ ., data = mtcars)
## 
## Prior probabilities of groups:
##       4       6       8 
## 0.34375 0.21875 0.43750 
## 
## Group means:
##        mpg     disp        hp     drat       wt     qsec        vs
## 4 26.66364 105.1364  82.63636 4.070909 2.285727 19.13727 0.9090909
## 6 19.74286 183.3143 122.28571 3.585714 3.117143 17.97714 0.5714286
## 8 15.10000 353.1000 209.21429 3.229286 3.999214 16.77214 0.0000000
##          am     gear     carb
## 4 0.7272727 4.090909 1.545455
## 6 0.4285714 3.857143 3.428571
## 8 0.1428571 3.285714 3.500000
## 
## Coefficients of linear discriminants:
##               LD1          LD2
## mpg   0.051343477  0.183229085
## disp  0.008008717 -0.008843922
## hp    0.024384632  0.044985875
## drat -0.416754228  1.767536660
## wt   -0.158957360  1.088695168
## qsec -0.314124831  0.396446017
## vs   -2.211101707 -1.942256185
## am   -1.213501776 -0.039106559
## gear -1.106168281  0.270037486
## carb -0.079499680 -1.485308396
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9327 0.0673
plot(lda(cyl ~., data=mtcars))

plot(lda(gear ~., data=mtcars))

plot(lda(carb ~., data=mtcars))

plot(lda(am ~., data=mtcars))

Теперь даже совсем деревянный классификатор все сумеет отличить.

Факторный анализ

Помогает вскрывать интересные зависимости там, где их не было видно:

factanal(mtcars, 3)
## 
## Call:
## factanal(x = mtcars, factors = 3)
## 
## Uniquenesses:
##   mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb 
## 0.135 0.055 0.090 0.127 0.290 0.060 0.051 0.223 0.208 0.125 0.158 
## 
## Loadings:
##      Factor1 Factor2 Factor3
## mpg   0.643  -0.478  -0.473 
## cyl  -0.618   0.703   0.261 
## disp -0.719   0.537   0.323 
## hp   -0.291   0.725   0.513 
## drat  0.804  -0.241         
## wt   -0.778   0.248   0.524 
## qsec -0.177  -0.946  -0.151 
## vs    0.295  -0.805  -0.204 
## am    0.880                 
## gear  0.908           0.224 
## carb  0.114   0.559   0.719 
## 
##                Factor1 Factor2 Factor3
## SS loadings      4.380   3.520   1.578
## Proportion Var   0.398   0.320   0.143
## Cumulative Var   0.398   0.718   0.862
## 
## Test of the hypothesis that 3 factors are sufficient.
## The chi square statistic is 30.53 on 25 degrees of freedom.
## The p-value is 0.205

Найдем количество факторов по p-value:

pvals <- c()
for (f in 1:6){
  pvals <- c(pvals, factanal(mtcars, f)$PVAL)
}
barplot(pvals, names.arg=1:6)

Остановимся на 5:

factanal(mtcars, 5)
## 
## Call:
## factanal(x = mtcars, factors = 5)
## 
## Uniquenesses:
##   mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb 
## 0.132 0.039 0.014 0.062 0.277 0.005 0.093 0.125 0.154 0.138 0.053 
## 
## Loadings:
##      Factor1 Factor2 Factor3 Factor4 Factor5
## mpg   0.639  -0.438  -0.474  -0.207         
## cyl  -0.612   0.691   0.297   0.111         
## disp -0.654   0.543   0.209   0.429   0.185 
## hp   -0.272   0.663   0.537   0.191   0.315 
## drat  0.807  -0.244                         
## wt   -0.730   0.231   0.429   0.473         
## qsec -0.178  -0.895  -0.257                 
## vs    0.260  -0.844  -0.215  -0.151   0.159 
## am    0.908                                 
## gear  0.898           0.228                 
## carb          0.476   0.841                 
## 
##                Factor1 Factor2 Factor3 Factor4 Factor5
## SS loadings      4.208   3.267   1.712   0.536   0.183
## Proportion Var   0.383   0.297   0.156   0.049   0.017
## Cumulative Var   0.383   0.680   0.835   0.884   0.901
## 
## Test of the hypothesis that 5 factors are sufficient.
## The chi square statistic is 3.11 on 10 degrees of freedom.
## The p-value is 0.979

А теперь можно рассуждать о факторах во все тяжкие.

Красиво можно интерпретировать факторы в movies:

library(ggplot2)
mov <- na.omit(movies[c('year', 'length', 'budget', 'rating', 'Action', 'Animation', 'Comedy', 'Drama', 'Documentary', 'Romance', 'Short')])
factanal(mov, 5)
## 
## Call:
## factanal(x = mov, factors = 5)
## 
## Uniquenesses:
##        year      length      budget      rating      Action   Animation 
##       0.850       0.005       0.060       0.835       0.851       0.949 
##      Comedy       Drama Documentary     Romance       Short 
##       0.005       0.005       0.964       0.928       0.005 
## 
## Loadings:
##             Factor1 Factor2 Factor3 Factor4 Factor5
## year        -0.112   0.342                   0.113 
## length       0.967           0.203                 
## budget       0.219   0.932   0.134                 
## rating                       0.152           0.375 
## Action       0.114   0.336          -0.135         
## Animation   -0.120   0.171                         
## Comedy      -0.116          -0.103   0.983         
## Drama        0.123  -0.150   0.954  -0.116   0.186 
## Documentary                 -0.140                 
## Romance      0.115           0.135   0.197         
## Short       -0.625          -0.191           0.749 
## 
##                Factor1 Factor2 Factor3 Factor4 Factor5
## SS loadings      1.461   1.169   1.082   1.063   0.768
## Proportion Var   0.133   0.106   0.098   0.097   0.070
## Cumulative Var   0.133   0.239   0.337   0.434   0.504
## 
## Test of the hypothesis that 5 factors are sufficient.
## The chi square statistic is 434.34 on 10 degrees of freedom.
## The p-value is 4.56e-87

Какие красивые истории можно придумать тут?

Классификация

Иерархическая классификация

Эта штука умеет превращать матрицу расстояний в классификационное дерево. Получить матрицу расстояния можно командой dist:

#dist(mtcars) # Она занимает много места

Расстояния бывают разные, лишь бы выполнялись три условия: Какие?

Остановимся на Евклиде. Теперь подаем эту штуку в hclust.

hcl <- hclust(dist(mtcars))
hcl
## 
## Call:
## hclust(d = dist(mtcars))
## 
## Cluster method   : complete 
## Distance         : euclidean 
## Number of objects: 32
plot(hcl)

Теперь мы просто обязаны проверить классификатор. Вот эти две машины оказались рядом:

И ее сосед:

Похоже, что все хорошо. Но где мы могли сильно ошибиться? А именно: почему Мазератти и Феррари упали в такие разные кластеры?

Шкалирование данных:

mtcars.scaled <- scale(mtcars, center = FALSE, scale = TRUE)
plot(hclust(dist(mtcars.scaled)))

Совсем другое дело.

Мораль: всегда будьте аккуратнее, когда имеете дело с расстояниями.

K-средних

Про алгоритм я тоже расскажу у доски.

km <- kmeans(mtcars.scaled, 3)
km
## K-means clustering with 3 clusters of sizes 14, 11, 7
## 
## Cluster means:
##         mpg       cyl      disp        hp      drat        wt      qsec
## 1 0.7094773 1.2241280 1.3316418 1.2753130 0.8744291 1.1720913 0.9204247
## 2 1.2228980 0.6955273 0.4161788 0.5818648 1.1050325 0.6555143 0.9807230
## 3 0.9746084 0.7869394 0.6604064 0.6226349 0.9666881 0.9361825 1.0957605
##          vs        am      gear      carb
## 1 0.0000000 0.2206029 0.8604813 1.0662908
## 2 0.9469394 1.5442200 1.1189658 0.7200925
## 3 1.4880476 0.0000000 0.9353057 0.6528311
## 
## Clustering vector:
##           Mazda RX4       Mazda RX4 Wag          Datsun 710 
##                   2                   2                   2 
##      Hornet 4 Drive   Hornet Sportabout             Valiant 
##                   3                   1                   3 
##          Duster 360           Merc 240D            Merc 230 
##                   1                   3                   3 
##            Merc 280           Merc 280C          Merc 450SE 
##                   3                   3                   1 
##          Merc 450SL         Merc 450SLC  Cadillac Fleetwood 
##                   1                   1                   1 
## Lincoln Continental   Chrysler Imperial            Fiat 128 
##                   1                   1                   2 
##         Honda Civic      Toyota Corolla       Toyota Corona 
##                   2                   2                   3 
##    Dodge Challenger         AMC Javelin          Camaro Z28 
##                   1                   1                   1 
##    Pontiac Firebird           Fiat X1-9       Porsche 914-2 
##                   1                   2                   2 
##        Lotus Europa      Ford Pantera L        Ferrari Dino 
##                   2                   1                   2 
##       Maserati Bora          Volvo 142E 
##                   1                   2 
## 
## Within cluster sum of squares by cluster:
## [1] 10.606567  9.992138  1.866397
##  (between_SS / total_SS =  65.4 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"    
## [5] "tot.withinss" "betweenss"    "size"         "iter"        
## [9] "ifault"

Как выбрать оптимальное число кластеров? Визуально!

ss <- c()
for (k in 2:10) {
  km <- kmeans(mtcars.scaled, k, nstart = 25) # Попробуйте уменьшить nstart
  ss <- c(ss, km$betweenss/km$totss)
}
plot(2:10, ss, type="l")

Остановимся на четырех кластерах:

km <- kmeans(mtcars, 4, nstart = 25)
km$cluster
##           Mazda RX4       Mazda RX4 Wag          Datsun 710 
##                   2                   2                   2 
##      Hornet 4 Drive   Hornet Sportabout             Valiant 
##                   1                   4                   1 
##          Duster 360           Merc 240D            Merc 230 
##                   3                   2                   2 
##            Merc 280           Merc 280C          Merc 450SE 
##                   2                   2                   1 
##          Merc 450SL         Merc 450SLC  Cadillac Fleetwood 
##                   1                   1                   4 
## Lincoln Continental   Chrysler Imperial            Fiat 128 
##                   4                   4                   2 
##         Honda Civic      Toyota Corolla       Toyota Corona 
##                   2                   2                   2 
##    Dodge Challenger         AMC Javelin          Camaro Z28 
##                   1                   1                   3 
##    Pontiac Firebird           Fiat X1-9       Porsche 914-2 
##                   4                   2                   2 
##        Lotus Europa      Ford Pantera L        Ferrari Dino 
##                   2                   3                   2 
##       Maserati Bora          Volvo 142E 
##                   3                   2
table(km$cluster)
## 
##  1  2  3  4 
##  7 16  4  5

Разукрасим график кластерами:

plot(mpg ~ disp, data=mtcars, col=factor(km$cluster))

plot(drat ~ wt, data=mtcars, col=factor(km$cluster))

А знаете, где будет хорошее разделение?

plot(PCA$scores, col=factor(km$cluster))

Немного эконометрики

Прям вот совсем чуть-чуть. Будем для начала исследовать деревья.

Линейная модель

Будем исследовать деревья:

head(trees)
##   Girth Height Volume
## 1   8.3     70   10.3
## 2   8.6     65   10.3
## 3   8.8     63   10.2
## 4  10.5     72   16.4
## 5  10.7     81   18.8
## 6  10.8     83   19.7
dim(trees)
## [1] 31  3
hist(trees$Volume/(trees$Girth^2*trees$Height), main="\"Эффективность каждого дерева\"")

Обычная линейная модель:

model <- lm(data=trees, Volume ~ Girth + Height)
summary(model)
## 
## Call:
## lm(formula = Volume ~ Girth + Height, data = trees)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.4065 -2.6493 -0.2876  2.2003  8.4847 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -57.9877     8.6382  -6.713 2.75e-07 ***
## Girth         4.7082     0.2643  17.816  < 2e-16 ***
## Height        0.3393     0.1302   2.607   0.0145 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.882 on 28 degrees of freedom
## Multiple R-squared:  0.948,  Adjusted R-squared:  0.9442 
## F-statistic:   255 on 2 and 28 DF,  p-value: < 2.2e-16
vcov(model)
##             (Intercept)       Girth      Height
## (Intercept)  74.6189461  0.43217138 -1.05076889
## Girth         0.4321714  0.06983578 -0.01786030
## Height       -1.0507689 -0.01786030  0.01693933
plot(model, which = 1)

Проверим ручками (%*% – умножение матриц, t – транспонирование, solve – обратная):

y <- trees$Volume
X <- cbind(1, trees$Girth, trees$Height)
solve(t(X) %*% X) %*% t(X) %*% y
##             [,1]
## [1,] -57.9876589
## [2,]   4.7081605
## [3,]   0.3392512

Это тот случай, где уместны логарифмы:

model2 <- lm(data=trees, log(Volume) ~ log(Girth) + log(Height))
summary(model2)
## 
## Call:
## lm(formula = log(Volume) ~ log(Girth) + log(Height), data = trees)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.168561 -0.048488  0.002431  0.063637  0.129223 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.63162    0.79979  -8.292 5.06e-09 ***
## log(Girth)   1.98265    0.07501  26.432  < 2e-16 ***
## log(Height)  1.11712    0.20444   5.464 7.81e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08139 on 28 degrees of freedom
## Multiple R-squared:  0.9777, Adjusted R-squared:  0.9761 
## F-statistic: 613.2 on 2 and 28 DF,  p-value: < 2.2e-16
plot(model2, which = 1)

Что реальный мир думает про МНК

Допустим количество данных увеличилось. Теперь матрица X стала размером 10К на 100. 10К наблюдений и 100 регрессоров.

M <- 10
N <- 1e4
X <- cbind(1, matrix(rnorm(M*N), ncol = M))
y <- runif(N)
system.time(solve(t(X) %*% X) %*% t(X) %*% y)
##    user  system elapsed 
##   0.004   0.001   0.005

А теперь данных стало еще больше: N=10K, M=100.

M <- 500
N <- 1e4
X <- cbind(1, matrix(rnorm(M*N), ncol = M))
y <- runif(N)
system.time(solve(t(X) %*% X) %*% t(X) %*% y)
##    user  system elapsed 
##   6.538   0.196   7.464

Вот это оно и есть: проклятие размерности. А что делать, если N>1e6, M>1e4?

Ответ: численная оптимизация. Расскажу только на уровне идеи про градиентный спуск. Опять же, мне лень печатать, поэтому все у доски. Задача у метода простая: оптимизировать функцию потерь:

\[J(b) = \sum \left(\hat{y} - y\right)^2 = \sum \left(X\cdot b - y\right)^2\]

M <- 3
N <- 10
X <- cbind(1, matrix(rnorm(M*N), ncol = M))
y <- runif(N)
b0 <- rep(0, M + 1)
J <- function(b) {
  sum((X %*% b - y)^2)
}

lm(y ~ 0 + X)
## 
## Call:
## lm(formula = y ~ 0 + X)
## 
## Coefficients:
##      X1       X2       X3       X4  
## 0.59548  0.05003  0.13321  0.16426
optim(b0, J) # Это не градиентный спуск! Так как мы не передаем третий аргумент - функцию gr
## $par
## [1] 0.5954962 0.0500215 0.1331427 0.1642000
## 
## $value
## [1] 0.2693185
## 
## $counts
## function gradient 
##      197       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Регуляризация: LASSO и Ridge

Зачем нужна регуризация?

Ridge нормального человека:

Ridge курильщика: \[J(b) = \sum \left(X\cdot b - y\right)^2 + \lambda \sum b^2\]

LASSO: \[J(b) = \sum \left(X\cdot b - y\right)^2 + \lambda \sum |b|\]

В чем принципиальное отличие? В чем преимущество каждого? (Хочется это обсудить, поэтому тут ответов не будет)

Вероятностные модели

Логит, пробит, tanh – все суть об одном и том же. И главное, все можно решить градиентным спуском.

data <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv") # Ссылка на данные: http://www.ats.ucla.edu/stat/r/dae/probit.htm
data$rank <- factor(data$rank)
head(data)
##   admit gre  gpa rank
## 1     0 380 3.61    3
## 2     1 660 3.67    3
## 3     1 800 4.00    1
## 4     1 640 3.19    4
## 5     0 520 2.93    4
## 6     1 760 3.00    2

Разные модели:

myprobit <- glm(admit ~ gre + gpa + rank, family=binomial(link="probit"), data=data)
summary(myprobit)
## 
## Call:
## glm(formula = admit ~ gre + gpa + rank, family = binomial(link = "probit"), 
##     data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6163  -0.8710  -0.6389   1.1560   2.1035  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.386836   0.673946  -3.542 0.000398 ***
## gre          0.001376   0.000650   2.116 0.034329 *  
## gpa          0.477730   0.197197   2.423 0.015410 *  
## rank2       -0.415399   0.194977  -2.131 0.033130 *  
## rank3       -0.812138   0.208358  -3.898 9.71e-05 ***
## rank4       -0.935899   0.245272  -3.816 0.000136 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 499.98  on 399  degrees of freedom
## Residual deviance: 458.41  on 394  degrees of freedom
## AIC: 470.41
## 
## Number of Fisher Scoring iterations: 4
confint(myprobit)
## Waiting for profiling to be done...
##                     2.5 %       97.5 %
## (Intercept) -3.7201050682 -1.076327713
## gre          0.0001104101  0.002655157
## gpa          0.0960654793  0.862610221
## rank2       -0.7992113929 -0.032995019
## rank3       -1.2230955861 -0.405008112
## rank4       -1.4234218227 -0.459538829
plot(myprobit, which = 1)

mylogit <- glm(admit ~ ., family=binomial(link="logit"), data=data)
summary(mylogit)
## 
## Call:
## glm(formula = admit ~ ., family = binomial(link = "logit"), data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6268  -0.8662  -0.6388   1.1490   2.0790  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.989979   1.139951  -3.500 0.000465 ***
## gre          0.002264   0.001094   2.070 0.038465 *  
## gpa          0.804038   0.331819   2.423 0.015388 *  
## rank2       -0.675443   0.316490  -2.134 0.032829 *  
## rank3       -1.340204   0.345306  -3.881 0.000104 ***
## rank4       -1.551464   0.417832  -3.713 0.000205 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 499.98  on 399  degrees of freedom
## Residual deviance: 458.52  on 394  degrees of freedom
## AIC: 470.52
## 
## Number of Fisher Scoring iterations: 4

Что делать, если классов много?

Измерение точности классификации

На выходи из logit/probit и т.д. мы получаем вероятность принадлежать классу 1. Определим точность нашего классификатора.

library(ROCR)
## Loading required package: gplots
## 
## Attaching package: 'gplots'
## 
## The following object is masked from 'package:stats':
## 
##     lowess
model <- glm(admit ~ ., family=binomial(link="probit"), data=data)
pred <- prediction(model$fitted.values, data$admit)
plot(performance(pred, measure = "acc"), main="Точность")

Какие бывают проблемы:

  1. Скошенные классы
  2. Цена ошибки
table(data$admit)
## 
##   0   1 
## 273 127
plot(performance(pred, measure = "fpr"), main="Ложно-положительные")

plot(performance(pred, measure = "fnr"), main="Ложно-отрицательные")

plot(performance(pred, measure = "tpr", x.measure = "fpr"), main="ROC")

# AUC
performance(prediction(model$fitted.values, data$admit), measure = "auc")@y.values[[1]]
## [1] 0.6925817

Что выбрать, если вы предсказываете наличие рака? А если спам на почтовом севере?

А какую метрику взять для качества регрессии?

Задание

Есть два пути:

  1. Более простой (для экономистов-математиков): выбрать многомерный датасет (да, можно из ДЗ по визуализации) и провести уменьшение размерности (PCA, LDA, факторы) и классификацию (иерархия, k-means).
  2. Более сложный (для экономистов-программистов): научиться находить коэффициенты обычного МНК, Ridge и logit через оптимизацию функции потерь. Можно писать свой градиентный спуск (что круто), либо использовать готовую функцию optim, передавая в нее градиент (что тоже хорошо). Сравнить время работы на малых данных (N<100, M<10) и средних (N>1e4, M>1e2, тут аккуратнее, стандартные методы могут повиснуть).