Семинар 9. Карты, вопросы :)
Июнь, 28, 2018
Новые пакеты этого семинара: * leaflet, units, geojsonio, maptools, rgeos
- пакет со статичными данными для регионов России с гитхаба:
devtools::install_github('akondrashov96/rusmaps')library(tidyverse) # обработка данных, графики...
library(skimr) # описательные статистики
library(rio) # импорт фантастического количества форматов данных
library(ISLR) # ещё данные
library(caret) # пакет для подбора параметров разных моделей
library(elasticnet) # LASSO
library(latex2exp) # формулы в подписях к графику
library(rusmaps) # подборка карт России
library(geojsonio) # чтение карт в формате geojson
library(leaflet) # рисование динамических карт в html
library(maptools) # утилиты для работы с картами
library(rgeos) # пересчёт картографических проекций
library(plm) # анализ панельных данных
library(texreg) # таблички со сравнением моделейСнова импортируем набор данных по успеваемости студентов и смотрим на него :)
educ <- import('data/xAPI-Edu-Data.csv')
skim(educ)## Skim summary statistics
## n obs: 480
## n variables: 17
##
## Variable type: character
## variable missing complete n min max empty n_unique
## Class 0 480 480 1 1 0 3
## gender 0 480 480 1 1 0 2
## GradeID 0 480 480 4 4 0 10
## NationalITy 0 480 480 2 11 0 14
## ParentAnsweringSurvey 0 480 480 2 3 0 2
## ParentschoolSatisfaction 0 480 480 3 4 0 2
## PlaceofBirth 0 480 480 3 11 0 14
## Relation 0 480 480 3 6 0 2
## SectionID 0 480 480 1 1 0 3
## Semester 0 480 480 1 1 0 2
## StageID 0 480 480 10 12 0 3
## StudentAbsenceDays 0 480 480 7 7 0 2
## Topic 0 480 480 2 9 0 12
##
## Variable type: integer
## variable missing complete n mean sd p0 p25 p50 p75 p100
## AnnouncementsView 0 480 480 37.92 26.61 0 14 33 58 98
## Discussion 0 480 480 43.28 27.64 1 20 39 70 99
## raisedhands 0 480 480 46.77 30.78 0 15.75 50 75 100
## VisITedResources 0 480 480 54.8 33.08 0 20 65 84 99
## hist
## ▇▆▅▅▅▃▃▂
## ▇▇▇▇▃▅▅▃
## ▇▆▃▃▃▆▆▃
## ▆▃▂▁▃▃▇▇Данные взяты с kaggle.com. Целевая переменная — успеваемость студента, Class, принимает три значения, ‘H’ (high), ‘M’ (mid), ‘L’ (low). Остальные переменные — характеристики студента.
Повторяем шаги предыдущего семинара:
- Для целей бинарной классификации объединяем две верхних категории в одну
- Уберём старую переменную ‘Class’ с тремя категориями
- Объявляем все текстовые переменные факторными:
educ_logit <- mutate(educ, y = fct_collapse(Class, H = c('M', 'H'))) %>%
select(-Class)
educ_fct <- mutate_if(educ_logit, is.character, factor)- Разбиваем выборку на две части: 80% в обучающей и 20% в тестовой
set.seed(777)
train_rows <- createDataPartition(educ_fct$y, p = 0.8, list = FALSE)
educ_train <- educ_fct[train_rows, ]
educ_test <- educ_fct[-train_rows, ]1 Кросс-валидация для k ближайших соседей
Метод классификации k-ближайших соседей реализует пословицу «Скажи мне, кто твой друг, и я скажу, кто ты». Мы прогнозируем класс нового наблюдения, как наиболее часто встречающийся класс у ближайших соседей данного наблюдения.
Реализуем метод k ближайших соседей для набора данных про студентов. Подбирать число соседей будем с помощью кросс-валидации. Для этого у объекта trainControl зададим опции method = 'cv' и number = 5. Первая отвечает за тип кросс-валидации, а вторая за колиество разбиений выборки. В настройках обучения train укажем метод knn, добавим параметры кросс-валидции из объекта trxrl, а также попросим центрировать и отмасштабировать данные.
trctrl <- trainControl(method = 'cv', number = 5)
set.seed(777)
educ_knn_fit <- train(y ~ ., data = educ_train, method = 'knn',
trControl = trctrl,
preProcess = c('center', 'scale'))
educ_knn_fit## k-Nearest Neighbors
##
## 385 samples
## 16 predictor
## 2 classes: 'H', 'L'
##
## Pre-processing: centered (60), scaled (60)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 308, 308, 308, 309, 307
## Resampling results across tuning parameters:
##
## k Accuracy Kappa
## 5 0.8285960 0.5182881
## 7 0.8129081 0.4794467
## 9 0.8259644 0.5090629
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 5.Визуализируем долю правильных ответов в зависимости от числа соседей.
plot(educ_knn_fit, xlab = 'Число соседей',
ylab = 'Доля правильных ответов',
main = 'Метод k ближаших соседей')
Вспомним данные о кредитах :)
def <- Default
skim(def)## Skim summary statistics
## n obs: 10000
## n variables: 4
##
## Variable type: factor
## variable missing complete n n_unique top_counts
## default 0 10000 10000 2 No: 9667, Yes: 333, NA: 0
## student 0 10000 10000 2 No: 7056, Yes: 2944, NA: 0
## ordered
## FALSE
## FALSE
##
## Variable type: numeric
## variable missing complete n mean sd p0 p25
## balance 0 10000 10000 835.37 483.71 0 481.73
## income 0 10000 10000 33516.98 13336.64 771.97 21340.46
## p50 p75 p100 hist
## 823.64 1166.31 2654.32 ▅▆▇▆▃▁▁▁
## 34552.64 43807.73 73554.23 ▁▆▆▆▇▅▁▁Разбиваем данные о кредитах на две части:
set.seed(777)
train_rows <- createDataPartition(def$default, p = 0.7, list = FALSE)
def_train <- def[train_rows, ]
def_test <- def[-train_rows, ]- Упражнение 1.
Реализуйте метод k ближайших соседей на данных о кредитах. Не забудьте центрировать и отмасштабировать их! Затем визуализируйте полученный результат.
def_knn_fit <- train(default ~ ., data = def_train, method = 'knn',
trControl = trctrl,
preProcess = c('center', 'scale'))
plot(def_knn_fit, xlab = 'Число соседей',
ylab = 'Доля правильных ответов',
main = 'Метод k ближаших соседей')
Теперь получим прогнозы для тестовой выборки и оценим их качество.
educ_knn_predict <- predict(educ_knn_fit, newdata = educ_test)
confusionMatrix(educ_knn_predict, educ_test$y)## Confusion Matrix and Statistics
##
## Reference
## Prediction H L
## H 65 7
## L 5 18
##
## Accuracy : 0.8737
## 95% CI : (0.7897, 0.933)
## No Information Rate : 0.7368
## P-Value [Acc > NIR] : 0.0009533
##
## Kappa : 0.6657
## Mcnemar's Test P-Value : 0.7728300
##
## Sensitivity : 0.9286
## Specificity : 0.7200
## Pos Pred Value : 0.9028
## Neg Pred Value : 0.7826
## Prevalence : 0.7368
## Detection Rate : 0.6842
## Detection Prevalence : 0.7579
## Balanced Accuracy : 0.8243
##
## 'Positive' Class : H
## - Упражнение 2.
Постройте прогнозы метода k ближайших соседей для данных по кредитам и оцените их качество.
def_knn_pred <- predict(def_knn_fit, newdata = def_test)
confusionMatrix(def_knn_pred, def_test$default)## Confusion Matrix and Statistics
##
## Reference
## Prediction No Yes
## No 2878 63
## Yes 22 36
##
## Accuracy : 0.9717
## 95% CI : (0.9651, 0.9773)
## No Information Rate : 0.967
## P-Value [Acc > NIR] : 0.08142
##
## Kappa : 0.4451
## Mcnemar's Test P-Value : 1.434e-05
##
## Sensitivity : 0.9924
## Specificity : 0.3636
## Pos Pred Value : 0.9786
## Neg Pred Value : 0.6207
## Prevalence : 0.9670
## Detection Rate : 0.9597
## Detection Prevalence : 0.9807
## Balanced Accuracy : 0.6780
##
## 'Positive' Class : No
## 2 Кросс-валидация для LASSO
Вернёмся к построению линейной модели регрессии.
Чтобы получить простую модель с хорошими предсказаниями, можно использовать LASSO-регрессию. От обычной она отличается тем, что в целевую функцию мы добавляем дополнительное слагаемое — сумму модулей всех коэффициентов кроме первого: \[ Q = \sum_{i=1}^n (y_i - \hat y_i)^2 + \lambda \sum_{i=2}^k |\beta_i|, \] где \(\sum_{i=1}^n (y_i - \hat y_i)^2\) — сумма квадратов остатков, \(\lambda\) — параметр регуляризации.
Добавляя сумму модулей коэффициентов, мы штрафуем модель за слишком большие веса. При значении параметра \(\lambda = 0\) эта модель совпадает с обычным МНК, а при росте \(\lambda\) всё большее число коэффициентов зануляется.
Поcтроим LASSO-регрессию для данных о студентах. В качестве зависимой переменной возьмём количество раз, когда студент поднимал руку, а в качестве объясняющих — все остальные. В параметрах обучения укажем method = 'lasso' и зададим явно значения параметров, среди которых мы хотим найти наилучший, аргументом tuneGrid.
the_grid <- expand.grid(.fraction = c(0.1, 0.2, 0.3, 0.5, 0.7))
the_control <- trainControl(method = 'cv', number = 5)
educ_lasso_fit <- train(raisedhands ~ ., data = educ_train,
method = 'lasso', trControl = the_control,
tuneGrid = the_grid)
educ_lasso_fit## The lasso
##
## 385 samples
## 16 predictor
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 307, 309, 307, 308, 309
## Resampling results across tuning parameters:
##
## fraction RMSE Rsquared MAE
## 0.1 20.39864 0.5801659 16.46697
## 0.2 19.71801 0.5944004 15.41892
## 0.3 20.09979 0.5813415 15.79221
## 0.5 20.43065 0.5735971 16.01194
## 0.7 20.66105 0.5662569 16.13775
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was fraction = 0.2.plot(educ_lasso_fit,
xlab = TeX('Доля $R^2$ от МНК'), ylab = 'RMSE',
main = 'LASSO регрессия')
- Упражнение 3.
Постройте прогнозы для тестовой выборки.
educ_lasso_pred <- predict(educ_lasso_fit, newdata = educ_test)3 Карты
3.1 Динамичные
Динамичные карты будем рисовать с помощью пакета leaflet. Для этого нам прежде всего понадобятся данные для карты в хитром формате .geojson, которые мы взяли у Code for America.
rus <- geojson_read('data/russia.geojson', what = 'sp')
class(rus)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"Добавим к каждому региону вектор из случайных чисел:
set.seed(777)
rus@data <- mutate(rus@data, num = rnorm(n = 83, mean = 100, sd = 30))Укажем цвет для каждого региона и подпись, которая будет появляться при наведении мышкой:
pal <- colorNumeric('BuPu', NULL) # задаём палитру цветов
rus@data <- mutate(rus@data,
reg_color = pal(num),
reg_label = paste0(name, ': ', formatC(num), ' попугаев'))
glimpse(rus@data)## Observations: 83
## Variables: 8
## $ name <fct> Бурятия, Карачаево-Черкесская республика, Сахалинск...
## $ cartodb_id <int> 10, 13, 11, 16, 18, 19, 3, 7, 8, 12, 15, 17, 42, 76...
## $ created_at <fct> 2013/12/04 04:23:51+01, 2013/12/04 04:23:51+01, 201...
## $ updated_at <fct> 2013/12/04 08:04:02+01, 2013/12/04 08:04:37+01, 201...
## $ name_latin <fct> Republic of Buryatia, Karachay-Cherkess Republic, S...
## $ num <dbl> 114.69359, 88.04376, 115.32509, 88.03564, 149.16058...
## $ reg_color <chr> "#8C6FB3", "#97ACD2", "#8C6DB2", "#97ACD2", "#821B8...
## $ reg_label <chr> "Бурятия: 114.7 попугаев", "Карачаево-Черкесская ре...- Что находится в переменных
reg_colorиreg_label?
И теперь можно рисовать карту! Мы раскрасим все регионы России в зависимости от случайного числа num, добавим подпись при наведении курсора и легенду, а ещё попросим показывать не весь мир целиком, а только Россию.
rus %>%
leaflet() %>%
addTiles() %>% # инициализируем карту
addPolygons(stroke = FALSE, fillColor = ~reg_color, # раскрашиваем области
label = ~reg_label) %>% # добавляем подписи
addLegend(pal = pal, values = ~num, title = 'Попугаи') %>% # добавляем легенду
setView(lng = 100, lat = 66, zoom = 2) # просим показывать РоссиюЗоркий чукча, глядя на эту карту, обязательно скажет “Однако!”
Проблема состоит в том, что Россия — одна из двух стран, пересекающих 180-ый меридиан. Примерно по 180-му меридиану проходит линия перемены дат. Слева — “вчера”, справа — “сегодня”:
Левый остров — Ратманова (Россия), правый — Крузенштерна (США).
Здесь нам пришлось вручную поправить скачанный geojson. Технически Чукотка реализована в виде набора полигонов. Мы прибавили к отрицательным долготам некоторых полигонов 360 градусов:
chukotka <- rus@polygons[[18]]
for (i in 1:length(chukotka@Polygons)) {
polygon_long <- chukotka@Polygons[[i]]@coords[, 1]
if (mean(polygon_long) < 0) {
polygon_long <- 360 + polygon_long
}
chukotka@Polygons[[i]]@coords[, 1] <- polygon_long
}
rus_corrected <- rus
rus_corrected@polygons[[18]] <- chukotkaИ заново строим карту:
rus_corrected %>%
leaflet() %>%
addTiles() %>% # инициализируем карту
addPolygons(stroke = FALSE, fillColor = ~reg_color, # раскрашиваем области
label = ~reg_label) %>% # добавляем подписи
addLegend(pal = pal, values = ~num, title = 'Попугаи') %>% # добавляем легенду
setView(lng = 100, lat = 66, zoom = 2) # просим показывать РоссиюСохранить карту в статичный png можно кликнув Export над картинкой в Rstudio.
Обновлённую карту со скорректированной Чукоткой и данными по попугаям можно сохранить командой geojson_write() из пакета geojsonio.
geojson_write(rus, file = 'russia_new.geojson')## <geojson-file>
## Path: russia_new.geojson
## From class: SpatialPolygonsDataFrameПро динамичные карты можно прочитать подробнее в документации, а попрактиковаться — на datacamp.
- Упражнение 4.
Отберите из старой таблички rus_corrected@data только содержательные колонки:
data_old <- rus_corrected@data %>% select(name, num, cartodb_id)
glimpse(data_old)## Observations: 83
## Variables: 3
## $ name <fct> Бурятия, Карачаево-Черкесская республика, Сахалинск...
## $ num <dbl> 114.69359, 88.04376, 115.32509, 88.03564, 149.16058...
## $ cartodb_id <int> 10, 13, 11, 16, 18, 19, 3, 7, 8, 12, 15, 17, 42, 76...Экспортируйте таблицу data_old c названиями регионов и их идентификаторами в формат csv:
export(data_old, 'data_old.csv')В экселе руками поменяйте столбец num и импортируйте эту таблицу обратно в R.
new_data <- import('data_old_corrected.csv')Добавьте в таблицу new_data код цвета и подпись для каждого региона:
pal <- colorNumeric('BuPu', NULL) # задаём палитру цветов
new_data2 <- mutate(new_data,
reg_color = pal(num),
reg_label = paste0(name, ': ', formatC(num)))Подменяем старые данные внутри нашей карты отредактированными new_data2:
rus_corrected@data <- new_data2И постройте карту России, раскрашенную по регионам в соответсвии с вашими данными. Добавьте легенду и подписи.
rus_corrected %>%
leaflet() %>%
addTiles() %>% # инициализируем карту
addPolygons(stroke = FALSE, fillColor = ~reg_color, # раскрашиваем области
label = ~reg_label) %>% # добавляем подписи
addLegend(pal = pal, values = ~num, title = 'Попугаи') %>% # добавляем легенду
setView(lng = 100, lat = 66, zoom = 2) # просим показывать Россию3.2 Статические
Куча готовых и очень подробных, даже слишком, карт по России в пакете Артёма Кондрашова rusmaps. К пакету есть шикарная документация.
Смотрим богатства пакета rusmaps. Cписок карт, которые содержит пакет, в студию:
glimpse(rusmaps.dataframe)## Observations: 30
## Variables: 4
## $ map_name <fct> cities, Voronezh, Novgorod, Novosibirsk, Omsk, Per...
## $ description <fct> Cities of Russia with 1mil population, Voronezh ci...
## $ format <chr> "SpatialPolygonDataFrame", "SpatialPolygonDataFram...
## $ date <chr> "27.06.2016", "27.06.2016", "27.06.2016", "27.06.2...Для просмотра таблицы с доступными картами внутри Rstudio можно набрать View(rusmaps.dataframe).
Карты хранятся в старом стандарте, который можно изображать сразу:
plot(rus_fd)
Мы переведём карты в формат, удобный для рисовки с помощью ggplot2:
Извлекаем числовые данные о регионах, встроенные в старую карту:
info <- rus_fd@data
glimpse(info)## Observations: 9
## Variables: 5
## $ id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9
## $ name <fct> Северо-Западный федеральный округ, Сибирский федерал...
## $ admin_lvl <int> 3, 3, 3, 3, 3, 3, 3, 3, 3
## $ pop_2015 <int> 13843600, 19312200, 12275800, 6211000, 14003800, 297...
## $ pop_2016 <int> 13853694, 19324031, 12308103, 6194969, 14044580, 296...- Какая информация встроена в карту
rus_fdв слотdata? - Как называется переменная с названием региона?
Переводим старую карту в простую таблицу с потерей числовых данных:
rus_nodata <- fortify(rus_fd, region = 'name')
glimpse(rus_nodata)## Observations: 16,381
## Variables: 7
## $ long <dbl> 105.5083, 105.9485, 106.0648, 106.0713, 106.3421, 106.54...
## $ lat <dbl> 67.01577, 67.04022, 67.11533, 67.16792, 67.20077, 67.291...
## $ order <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 1...
## $ hole <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, ...
## $ piece <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ id <chr> "Дальневосточный федеральный округ", "Дальневосточный фе...
## $ group <fct> Дальневосточный федеральный округ.1, Дальневосточный фед...- Какие данные содержатся в таблице
rus_nodata? - Как называется переменная с названием региона?
Подклеиваем данные к нашей карте:
rus_final <- left_join(rus_nodata, info, by = c('id' = 'name'))
glimpse(rus_final)## Observations: 16,381
## Variables: 11
## $ long <dbl> 105.5083, 105.9485, 106.0648, 106.0713, 106.3421, 10...
## $ lat <dbl> 67.01577, 67.04022, 67.11533, 67.16792, 67.20077, 67...
## $ order <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1...
## $ hole <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FAL...
## $ piece <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ id <chr> "Дальневосточный федеральный округ", "Дальневосточны...
## $ group <fct> Дальневосточный федеральный округ.1, Дальневосточный...
## $ id.y <int> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4...
## $ admin_lvl <int> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3...
## $ pop_2015 <int> 6211000, 6211000, 6211000, 6211000, 6211000, 6211000...
## $ pop_2016 <int> 6194969, 6194969, 6194969, 6194969, 6194969, 6194969...Базовый вариант карты России:
base <- ggplot(rus_final) +
geom_polygon(aes(x = long, y = lat,
fill = pop_2016, group = group), color = 'white') +
labs(x = '', y = '', fill = 'Население 2016', title = 'Население России')
base
Небольшое преобразование координат:
base + coord_quickmap()
Убираем легенду полностью:
base + coord_quickmap() + theme(legend.position = 'none')
Упражнение 5.
Сохраните табличку
infoв форматеcsv
export(info, 'info.csv')- Добавьте в экселе к ней новый столбец с произвольными данными
- Импотрируйте обновлённую табличку
new_info <- import('info_new.csv') # создали стобец hop_2016
rus_new <- left_join(rus_nodata, new_info, by = c('id' = 'name'))
glimpse(rus_new)## Observations: 16,381
## Variables: 12
## $ long <dbl> 105.5083, 105.9485, 106.0648, 106.0713, 106.3421, 10...
## $ lat <dbl> 67.01577, 67.04022, 67.11533, 67.16792, 67.20077, 67...
## $ order <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1...
## $ hole <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FAL...
## $ piece <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ id <chr> "Дальневосточный федеральный округ", "Дальневосточны...
## $ group <fct> Дальневосточный федеральный округ.1, Дальневосточный...
## $ id.y <int> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4...
## $ admin_lvl <int> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3...
## $ pop_2015 <int> 6211000, 6211000, 6211000, 6211000, 6211000, 6211000...
## $ pop_2016 <int> 6194969, 6194969, 6194969, 6194969, 6194969, 6194969...
## $ hop_2016 <int> 754, 754, 754, 754, 754, 754, 754, 754, 754, 754, 75...- Постройте статичную карту России с обновлёнными данными
new_base <- ggplot(rus_new) +
geom_polygon(aes(x = long, y = lat,
fill = hop_2016, group = group), color = 'white') +
labs(x = '', y = '', fill = 'Какой-то показатель', title = 'Какой-то показатель по России')
new_base
4 Немного панельных данных
Замечательная инструкция по анализу панельных данных в R.
Активируем из спячки встроенный в пакет plm набор данных по предприятиям:
data('Grunfeld', package='plm')
skim(Grunfeld)## Skim summary statistics
## n obs: 200
## n variables: 5
##
## Variable type: integer
## variable missing complete n mean sd p0 p25 p50 p75
## firm 0 200 200 5.5 2.88 1 3 5.5 8
## year 0 200 200 1944.5 5.78 1935 1939.75 1944.5 1949.25
## p100 hist
## 10 ▇▃▃▃▃▃▃▇
## 1954 ▇▅▇▅▅▇▅▇
##
## Variable type: numeric
## variable missing complete n mean sd p0 p25 p50 p75
## capital 0 200 200 276.02 301.1 0.8 79.17 205.6 358.1
## inv 0 200 200 145.96 216.88 0.93 33.56 57.48 138.04
## value 0 200 200 1081.68 1314.47 58.12 199.97 517.95 1679.85
## p100 hist
## 2226.3 ▇▃▁▁▁▁▁▁
## 1486.7 ▇▁▁▁▁▁▁▁
## 6241.7 ▇▁▂▁▁▁▁▁glimpse(Grunfeld)## Observations: 200
## Variables: 5
## $ firm <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
## $ year <int> 1935, 1936, 1937, 1938, 1939, 1940, 1941, 1942, 1943, ...
## $ inv <dbl> 317.6, 391.8, 410.6, 257.7, 330.8, 461.2, 512.0, 448.0...
## $ value <dbl> 3078.5, 4661.7, 5387.1, 2792.2, 4313.2, 4643.9, 4551.2...
## $ capital <dbl> 2.8, 52.6, 156.9, 209.2, 203.4, 207.2, 255.2, 303.7, 2...В панельных данных важно объявить переменные отвечающие за время и номер объекта:
grun <- pdata.frame(Grunfeld, index=c('firm', 'year'))
skim(grun)## Skim summary statistics
## n obs: 200
## n variables: 5
##
## Variable type: factor
## variable missing complete n n_unique top_counts
## firm 0 200 200 10 1: 20, 2: 20, 3: 20, 4: 20
## year 0 200 200 20 193: 10, 193: 10, 193: 10, 193: 10
## ordered
## FALSE
## FALSE
##
## Variable type: numeric
## variable missing complete n mean sd p0 p25 p50 p75
## capital 0 200 200 276.02 301.1 0.8 79.17 205.6 358.1
## inv 0 200 200 145.96 216.88 0.93 33.56 57.48 138.04
## value 0 200 200 1081.68 1314.47 58.12 199.97 517.95 1679.85
## p100 hist
## 2226.3 ▇▃▁▁▁▁▁▁
## 1486.7 ▇▁▁▁▁▁▁▁
## 6241.7 ▇▁▂▁▁▁▁▁glimpse(grun)## Observations: 200
## Variables: 5
## $ firm <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
## $ year <fct> 1935, 1936, 1937, 1938, 1939, 1940, 1941, 1942, 1943, ...
## $ inv <S3: pseries> 317.6, 391.8, 410.6, 257.7, 330.8, 461.2, 512....
## $ value <S3: pseries> 3078.5, 4661.7, 5387.1, 2792.2, 4313.2, 4643.9...
## $ capital <S3: pseries> 2.8, 52.6, 156.9, 209.2, 203.4, 207.2, 255.2, ...Теперь можно оценивать FE модель :)
grun_fe <- plm(inv ~ value + capital, data = Grunfeld, model = "within")
summary(grun_fe)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = inv ~ value + capital, data = Grunfeld, model = "within")
##
## Balanced Panel: n = 10, T = 20, N = 200
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -184.00857 -17.64316 0.56337 19.19222 250.70974
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## value 0.110124 0.011857 9.2879 < 2.2e-16 ***
## capital 0.310065 0.017355 17.8666 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 2244400
## Residual Sum of Squares: 523480
## R-Squared: 0.76676
## Adj. R-Squared: 0.75311
## F-statistic: 309.014 on 2 and 188 DF, p-value: < 2.22e-16И RE модель:
grun_re <- plm(inv ~ value + capital, data = Grunfeld, model = "random")
summary(grun_re)## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = inv ~ value + capital, data = Grunfeld, model = "random")
##
## Balanced Panel: n = 10, T = 20, N = 200
##
## Effects:
## var std.dev share
## idiosyncratic 2784.46 52.77 0.282
## individual 7089.80 84.20 0.718
## theta: 0.8612
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -177.6063 -19.7350 4.6851 19.5105 252.8743
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) -57.834415 28.898935 -2.0013 0.04674 *
## value 0.109781 0.010493 10.4627 < 2e-16 ***
## capital 0.308113 0.017180 17.9339 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 2381400
## Residual Sum of Squares: 548900
## R-Squared: 0.7695
## Adj. R-Squared: 0.76716
## F-statistic: 328.837 on 2 and 197 DF, p-value: < 2.22e-16Или в первых разностях:
grun_fd <- plm(inv ~ value + capital, data = Grunfeld, model = "fd")
summary(grun_fd)## Oneway (individual) effect First-Difference Model
##
## Call:
## plm(formula = inv ~ value + capital, data = Grunfeld, model = "fd")
##
## Balanced Panel: n = 10, T = 20, N = 200
## Observations used in estimation: 190
##
## Residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -202.05 -15.23 -1.76 -1.39 7.95 199.27
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## value 0.0890628 0.0082341 10.816 < 2.2e-16 ***
## capital 0.2786940 0.0471564 5.910 1.58e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 584410
## Residual Sum of Squares: 345940
## R-Squared: 0.40876
## Adj. R-Squared: 0.40561
## F-statistic: 129.597 on 1 and 188 DF, p-value: < 2.22e-16Табличку со сравнением всех трёх моделей в студию:
model_names <- c('RE-модель', 'FE-модель', 'FD-модель')
htmlreg(list(grun_re, grun_fe, grun_fd),
custom.model.names = model_names,
star.symbol = '\\*')| RE-модель | FE-модель | FD-модель | ||
|---|---|---|---|---|
| (Intercept) | -57.83* | |||
| (28.90) | ||||
| value | 0.11*** | 0.11*** | 0.09*** | |
| (0.01) | (0.01) | (0.01) | ||
| capital | 0.31*** | 0.31*** | 0.28*** | |
| (0.02) | (0.02) | (0.05) | ||
| R2 | 0.77 | 0.77 | 0.41 | |
| Adj. R2 | 0.77 | 0.75 | 0.41 | |
| Num. obs. | 200 | 200 | 190 | |
| ***p < 0.001, **p < 0.01, *p < 0.05 | ||||
Тест Хаусмана на асимптотическое равенство RE и FE оценок:
phtest(grun_re, grun_fe)##
## Hausman Test
##
## data: inv ~ value + capital
## chisq = 2.3304, df = 2, p-value = 0.3119
## alternative hypothesis: one model is inconsistentНулевая гипотеза о том, что обе модели дают одинаковые оценки коэффициентов, не отвергается. Поэтому мы считаем, что обе модели состоятельны :)
Ура-Ура-Ура :)