Загружаем нужные библиотеки:
library(knitr)
opts_chunk$set(cache = TRUE) # кэшируем куски Rmd файла для скоростной компиляции после
library(plyr) # ddply
library(reshape) # melt - cast
library(kernlab) # SVM
library(ggplot2) # графики
Многие преобразования могут быть описаны с помощью принципа “Разделяй и властвуй! Затем соединяй!”:
Загрузим данные по стоимости квартир в Москве:
file <- "~/science/econometrix/em301/datasets/flats_moscow.txt"
h <- read.table(file, header = TRUE)
Пример. Подсчёт описательных статистик внутри групп
Мы разбиваем большую таблицу на маленькие таблички по переменным code (географический район квартиры) и brick (кирпичность). По каждой маленькой табличке считаем число наблюдений и среднюю цену квартиры. Затем сводим в итоговую таблицу:
h.summ <- ddply(h, ~code + brick, summarize, n = length(price), mean = mean(price))
head(h.summ)
## code brick n mean
## 1 1 0 180 123.0
## 2 1 1 92 156.7
## 3 2 0 186 106.0
## 4 2 1 29 138.3
## 5 3 0 169 134.7
## 6 3 1 176 161.2
Пример. Корректировка наблюдения на групповое среднее. Из цены каждой квартиры будет вычтена средняя цена по квартирам данного района и данной кирпичности.
h.new <- ddply(h, ~code + brick, transform, delta = price - mean(price))
head(h.new)
## n price totsp livesp kitsp dist metrdist walk brick floor code delta
## 1 4 95 61 37 6 13.5 7 1 0 1 1 -28
## 2 24 95 58 38 6 12.0 8 1 0 1 1 -28
## 3 38 120 77 47 12 13.5 10 1 0 1 1 -3
## 4 47 129 79 46 13 12.0 15 1 0 1 1 6
## 5 67 135 77 47 10 13.5 10 0 0 1 1 12
## 6 68 100 78 44 8 13.5 20 0 0 1 1 -23
Итого:
transform в финальной табличке создаёт столько строк, сколько было в изначальной табличкеsummarize в финальной табличке создаёт столько строк, на сколько мелких табличек разрезали изначальную Из опыта. Для подсчета числа наблюдений бывает полезна функция nrow(). Функция nrow() работает только без имени переменной и без summarize.
Доказательство.
Шаг 1. Крокодил более широкий, чем зеленый. Зеленый он только снаружи, а широкий он и снаружи и внутри.
Шаг 2. Крокодил более зеленый, чем длинный. Зеленый он и в длину, и в ширину, а длинный - только в длину.
Шаг 3. Отношение “больше” транзитивно :)
Пример “широкой” таблицы
tab.wide <- smiths
head(tab.wide)
## subject time age weight height
## 1 John Smith 1 33 90 1.87
## 2 Mary Smith 1 NA NA 1.54
Пример “длинной” таблицы
names(airquality) <- tolower(names(airquality))
tab.long <- melt(airquality, id = c("month", "day"), na.rm = TRUE)
head(tab.long)
## month day variable value
## 1 5 1 ozone 41
## 2 5 2 ozone 36
## 3 5 3 ozone 12
## 4 5 4 ozone 18
## 5 5 6 ozone 28
## 6 5 7 ozone 23
Данные удобнее обрабатывать в “длинном” формате. В “длинном” формате при получении новых наблюдений не нужно добавлять столбцы в data.frame. “Широкий” формат изредка бывает удобен для представления результатов. Или клиент так хочет.
Чтобы табличка растаяла из “широкого” формата и стекла вниз в “длинный” её можно растопить командой melt из пакета reshape.
tab.long2 <- melt(tab.wide, id = "subject")
head(tab.long2)
## subject variable value
## 1 John Smith time 1
## 2 Mary Smith time 1
## 3 John Smith age 33
## 4 Mary Smith age NA
## 5 John Smith weight 90
## 6 Mary Smith weight NA
Обратное действие выполняется командой cast.
tab.wide2 <- cast(tab.long, day + month ~ variable, mean)
head(tab.wide2)
## day month ozone solar.r wind temp
## 1 1 5 41 190 7.4 67
## 2 1 6 NaN 286 8.6 78
## 3 1 7 135 269 4.1 84
## 4 1 8 39 83 6.9 81
## 5 1 9 96 167 6.9 91
## 6 2 5 36 118 8.0 72
Во многих моделях есть параметр, отвечающий за простоту модели. Чем проще модель, тем хуже модель описывает выборку, по которой она оценивалась.
В линейной регрессии параметр сложности модели — это количество регрессоров, \( k \). Чем больше регрессоров, тем ниже сумма квадратов остатков, \( RSS \).
На гистограмме параметр сложности — число столбцов.
В ridge regression и LASSO параметр простоты — это \( \lambda \). Действительно, LASSO минимизирует \[ \sum_{i=1}^n (y_i-\hat{y}_i)^2+\lambda \sum_{j=1}^k |\hat{\beta}_j| \] Значит, чем больше параметр \( \lambda \), тем больше стремление алгоритма LASSO занулить некоторые \( \hat{\beta}_j \).
Этот параметр сложности не так легко оценить, как коэффициенты модели. Если наивно попытаться выбрать параметр сложности так, чтобы модель как можно лучше описывала бы выборку, по которой она оценивалась… То ничего хорошего не выйдет. Окажется, что оптимальное количество регрессоров равно плюс бесконечности, а оптимальное \( \lambda \) в LASSO равно нулю. В регрессии одно из решений этой проблемы — это проверка гипотез о значимости коэффициентов.
Есть и универсальный способ выбора сложности модели — кросс валидация (перекрёстная проверка, cross validation). Её идея состоит в том, что надо оценивать качество прогнозов не по той же выборке, на которой оценивалась модель, а на новых наблюдениях.
set.seed() для воспроизводимости эксперимента.Популярные значения k:
Создаём табличку перебираемых \( C \) и \( \sigma \):
C <- c(1, 10, 100)
sigma <- c(0.1, 1, 10)
d <- expand.grid(C, sigma) # Случайно не декартово ли произведение это? :)
colnames(d) <- c("C", "sigma")
head(d)
## C sigma
## 1 1 0.1
## 2 10 0.1
## 3 100 0.1
## 4 1 1.0
## 5 10 1.0
## 6 100 1.0
Пробуем 10-кратную кросс-валидацию для конкретных \( C \) и \( \sigma \):
set.seed(33222233) # любимый seed Фрекен Бок
m1 <- ksvm(data = h, brick ~ ., kernel = rbfdot, kpar = list(sigma = 1), C = 1,
cross = 10)
cross(m1) # cross validation error
## [1] 0.1669
Для удобства оформляем это в функцию двух переменных
f_cross <- function(sig, C) {
model <- ksvm(data = h, brick ~ ., kernel = rbfdot, kpar = list(sigma = sig),
C = C, cross = 10)
return(cross(model))
}
f_cross(1, 1) # тестируем сделанную функцию
## [1] 0.1655
Применяем её ко всем возможным \( C \) и \( \sigma \)
d <- ddply(d, ~C + sigma, transform, cr = f_cross(sigma, C))
head(d)
## C sigma cr
## 1 1 0.1 0.1691
## 2 1 1.0 0.1644
## 3 1 10.0 0.2101
## 4 10 0.1 0.1792
## 5 10 1.0 0.1763
## 6 10 10.0 0.2065
Ссылки на дополнительные материалы:
ps Придумайте хороший перевод названия функции cast на русский язык.