На 3-м курсе факультета экономики вышки есть Самая Важная дисциплина. Это — Эконометрика. А этот текст — моя рабочая тетрадь R, в которой я готовлю семинары.
Страница курса, https://github.com/bdemeshev/em301/wiki — там размещены все наборы данных и другие материалы по курсу.
Если есть вопросы — пишите на boris.demeshev@gmail.com.
Борис Демешев
Linux/MacOs
1. Узнать текущую рабочую директорию:
pwd
cd target_folder
Не забывайте, что в командной строке работает
ls
mkdir the_new_folder
Rscript file_to_run.R
Windows:
1. Узнать текущую рабочую директорию:
cd
cd target_folder
Не забывайте, что в командной строке работает
dir
mkdir the_new_folder
C:\path_to_R\Rscript file_to_run.R
По умолчанию в Windows для запуска программы нужно указывать её полный путь. Чтобы этого не делать каждый раз, а писать коротко Rscript file_to_run.R, нужно внести путь к папке с R в переменную $PATH. Для этого…getwd()
## [1] "/home/yuliana/boris_essex/em301"
setwd("/home/boris/")
## Error: невозможно сменить рабочий каталог
Не забывайтие, что в R-studio работает
3. Посмотреть список файлов в текущей директории:
dir()
## [1] "datasets" "em301_ht1.txt" "em301_ht1.txt~"
## [4] "em301_r.html" "em301_r.md" "em301_r.Rmd"
## [7] "em301_r.Rmd~" "figure" "hse_notes"
## [10] "README.md" "README.md~" "the_new_folder"
## [13] "working_plans.txt" "working_plans.txt~"
Можно указать пролистываемую папку, dir('folder')
4. Создать папку
dir.create("the_new_folder")
## Warning: 'the_new_folder' уже существует
savehistory("log_file.txt")
## Error: нет истории, чтобы сохранить
source("script_to_run.txt")
## Warning: не могу открыть файл 'script_to_run.txt': Нет такого файла или
## каталога
## Error: не могу открыть соединение
quit() или коротко q()На дворе 21 век, тем не менее:
Справка по конкретной команде
`?`(t.test)
Поиск нужной команды по ключевому слову
`?`(`?`(anova))
А дальше имеет смысл искать в интернете. Имеет смысл гуглить по фразе “something I need in R”. Если поиск не увенчался успехом, то тогда имеет смысл задать вопрос на одном из форумов:
Самый простой тип — количественные переменные.
x = c(1.5, 2.4, 3.6)
y = c("Аня", "Маша", "Ира", "Света")
z = c(TRUE, FALSE, TRUE)
w = c("Да", "Да", "Не знаю", "Нет", "Не знаю")
Более интересный тип, качественные данные. Текстовые можно легко преобразовать в качественные:
w2 = as.factor(w)
Качественные данные занимают в памяти комьютера меньше места и для них есть множество специальных функций. Частный случай, качественные порядковые переменные
stud = ordered(c("Школьник", "Бакалавр", "Бакалавр",
"Школьник", "Магистр", "Школьник", "Магистр",
"Бакалавр"), levels = c("Школьник", "Бакалавр",
"Магистр"))
stud
## [1] Школьник Бакалавр Бакалавр Школьник Магистр Школьник Магистр Бакалавр
## Levels: Школьник < Бакалавр < Магистр
Переменные могут быть организованы:
В векторы:
В матрицы:
Маленькая техническая подробность. На самом деле R хранит любую матрицу в виде вектора с дополнительной пометкой, атрибутом, о размере матрицы.
Узнаём размер матрицы:
dim(mat)
## Error: объект 'mat' не найден
Можно использовать и более длинную, но зато универсальную команду, подходяющую для всех атрибутов, attr(mat,'dim').
В списки:
Несколько переменных может быть организовано в таблицу данных dataframe
Переменную из таблицы данных можно вытащить с помощью df['x'] или, короче, df$x.
Помимо переменных у таблицы с данными может быть куча атрибутов. Атрибут может нести любую вспомогательную информацию. Самым распространённый атрибут, имена переменных, attr(df,'names'). В силу его распространённости для него есть короткая команда names(df).
Настоятельно рекомендую добавлять атрибут var_labels с более подробным описанием
переменных, подобно variable labels в Стате
attr(df, "var_labels") = c("Прожиточный минимум в 3042 году согласно исследованиям 2012 года",
"Просто какая-то переменная")
Чтобы узнать, что за объект перед нами, используем, str()
str(stud)
## Ord.factor w/ 3 levels "Школьник"<"Бакалавр"<..: 1 2 2 1 3 1 3 2
Загрузка файла STATA
library(foreign)
df = read.dta("myfile.dta")
## Error: unable to open file: 'Нет такого файла или каталога'
В качестве атрибутов сохраняются названия и описания переменных:
attr(df, "names")
## NULL
attr(df, "var.labels")
## NULL
Вместо attr(df,'names') можно писать коротко names(df). Атрибут var.labels используется реже, поэтому для него более короткой команды нет. (хм, а вдруг есть?)
Чтобы посмотреть, что ещё оказалось в таблице с данными, можно использовать функцию str(df).
Можно рисовать графики встроенными средствами R, но для более качественных графиков мы будем использовать пакет ggplot2
library(ggplot2)
Самые часто употребимые типы
plot(cars)
Линейная регрессия оценивается командой:
res <- lm(y ~ x + z, data = df)
## Error: аргумент 'data' неправильного типа
Отчёт о результатах
summary(res)
## Error: объект 'res' не найден
Из списка результатов для дальнейшей работы можно раздобыть коэффициенты, вектор остатков, вектор оценённых зависимых переменных и многое другое:
coef(res)
## Error: объект 'res' не найден
residuals(res)
## Error: объект 'res' не найден
fitted.values(res)
## Error: объект 'res' не найден
Посмотреть, что ещё есть в списке легко:
names(res)
## Error: объект 'res' не найден
Можно также использовать str(res).
Оценка ковариационной матрицы считается дополнительно,
vcov(res)
## Error: объект 'res' не найден
Для установки пакетов нужно соединение с интернетом. Установим, например, пакет vcd для мозаичных графиков
install.packages("vcd")
## Installing package(s) into
## '/home/yuliana/R/i686-redhat-linux-gnu-library/2.15' (as 'lib' is
## unspecified)
## Error: trying to use CRAN without setting a mirror
После установки пакета на компьютер, его можно загрузить
library(vcd)
## Loading required package: MASS
## Loading required package: grid
## Loading required package: colorspace
И использовать…
В реальности редко приходится вручную максимизировать функцию правдоподобия. Для популярных моделей в R уже есть готовые функции. Для замысловатых экзотических моделей скорее всего имеет смысл воспользоваться MCMC, так как ML будет работать плохо. Однако в учебных целях это просто обязательно!
Хорошие мелочи:
mlog <- function(params, data) {
# это функция вычисляет функцию правдоподобия со знаком минус
logl = tra - ta - ta
return(-logl)
}
Можно применять метод MCMC используя только сам R. Но это не самый продуктивный способ, для этого существенно более быстрые специальные пакеты. Мы будем использовать JAGS. Это отдельный пакет, в котором есть свой язык описания моделей, и его, в принципе, можно использовать отдельно от R. Для вызова JAGS из R, нужно установить дополнительный пакет, rjags. Делается это командой install.packages('rjags') из R.
library(rjags)
## Loading required package: coda
## Loading required package: lattice
## linking to JAGS 3.2.0
## module basemod loaded
## module bugs loaded
Байесовский подход за 15 минут.
Мы предполагаем некое априорное распределение параметров \( \theta \).
Мы предполагаем некую модель, как устроены наши данные.
Если модель проста, то применив формулу условной вероятности, можно получить апостериорное распределение параметров в явном виде.
Если модель сложна, то в явном виде апостериорное распределение получить не удаётся, но зато с помощью алгоритмов MCMC можно получить сгенерить выборку случайных величин из апостериорного распределения.
Пример 1. (одномерный с явными функциями)
Изначально я верю, что …
Получаем апостериорное распределение аналитически.
Посмотрим, как его получит JAGS
MCMC строит последовательность случайных величин, распределение которых стремится к искомому апостериорному распределению. … Поэтому для описания апостериорного распределения мы игнорируем начало последовательности случайных величин. burn in
Как JAGS это сделал?
Пример 2.
Происшествия на английских угольных шахтах
модель
Эта модель уже настолько сложна, что апостериорное распределение параметров в явном виде не выписывается!
Происшествия на английских угольных шахтах с JAGS:
Что делал JAGS за кадром…
Из кода видно, что нам не приходится выводить вручную условные распределения. На самом деле мы всего лишь описываем модель в дерева (направленного ациклического графа).
(картинка)
Из описания модели в виде графа можно получить условную функцию плотности параметра \( \theta \) по теореме
\[
f(\theta|G-\theta) \sim f(\theta|parents(\theta))\cdot \prod_{a
\in children(\theta)}f(a|parents(a))
\]
Схематично JAGS поступает так:
Что делать когда объемы данных - большие и памяти не хватает…