In bdemeshev/r_manual_book: R manual for beginners (in russian)
library(knitr)
library(tikzDevice)
activate_tikz <- function() {
# tikz plots options
options(tikzDefaultEngine = "xetex")
# cash font metrics for speed:
# options(tikzMetricsDictionary = "./tikz_metrics")
add_xelatex <- c("\\defaultfontfeatures{Ligatures=TeX, Scale=MatchLowercase}",
"\\setmainfont{Linux Libertine O}",
"\\setmonofont{Linux Libertine O}",
"\\setsansfont{Linux Libertine O}",
"\\newfontfamily{\\cyrillicfonttt}{Linux Libertine O}",
"\\newfontfamily{\\cyrillicfont}{Linux Libertine O}",
"\\newfontfamily{\\cyrillicfontsf}{Linux Libertine O}")
options(tikzXelatexPackages = c(getOption("tikzXelatexPackages"),
add_xelatex))
# does remove warnings:
# it is important to remove fontenc package wich is loaded by default
options(tikzUnicodeMetricPackages = c("\\usetikzlibrary{calc}",
"\\usepackage{fontspec, xunicode}", add_xelatex))
opts_chunk$set(dev = "tikz", dev.args = list(pointsize = 11))
}
activate_png <- function() {
opts_chunk$set(dev = "png", dpi = 300)
}
colFmt <- function(x, color) {
outputFormat <- opts_knit$get("rmarkdown.pandoc.to")
if (outputFormat == "latex") {
result <- paste0("\\textcolor{", color, "}{", x, "}")
} else if (outputFormat %in% c("html", "epub")) {
result <- paste0("<font color='", color, "'>", x, "</font>")
} else {
result <- x
}
return(result)
}
outputFormat <- opts_knit$get("rmarkdown.pandoc.to")
if (outputFormat == "latex") {
# activateTikz()
activate_png()
# другую тему для ggplot2 выставить?
}
library(tidyverse) # рабочие кони
library(sandwich) # оценка Var для гетероскедастичности
library(lmtest) # тест Бройша-Пагана
library(data.table) # манипуляции с данными
library(reshape2) # преобразование длинных таблиц в широкие
library(broom)
Статистика и не только {#statistics_and_more}
Пора вспомнить о том, что R ориентирован на статистику :)
Генерирование случайных величин {#random_variables}
Для решения задач по теории вероятностей или исследования свойств статистических алгоритмов может потребоваться сгененировать случайную выборку из заданного закона распределения.
Генерируем 10 равномерных на отрезке $[4;10.5]$ случайных величин:
runif(10, min = 4, max = 10.5)
Генерируем 10 нормальных $N(2;9)$ случайных величин с математическим ожиданием $2$ и дисперсией $9=3^2$:
rnorm(10, mean = 2, sd = 3)
Например, с помощью симуляций легко оценить математическое ожидание $E(1/X)$, где $X \sim N(2;9)$. Для этого мы вспомним Закон Больших Чисел. Он говорит, что арифметическое среднее по большой выборке стремится по вероятности и почти наверное к математическому ожиданию. Поэтому мы просто сгенерируем большую выборку в миллион наблюдений:
n_obs <- 10^6
x <- rnorm(n_obs, mean = 2, sd = 3)
mean(1/x)
Также легко оценить многие вероятности. Например, оценим вероятность $P(X_1 + X_2 + X_3^2 > 5)$, где величины $X_i$ независимы и одинаково распределены $X_i \sim U[0;2]$:
n_obs <- 10^6
x_1 <- runif(n_obs, min = 0, max = 2)
x_2 <- runif(n_obs, min = 0, max = 2)
x_3 <- runif(n_obs, min = 0, max = 2)
success <- x_1 + x_2 + x_3^2 > 5
sum(success) / n_obs
Здесь вектор success будет содержать значение TRUE там, где условие x_1 + x_2 + x_3^2 > 5 выполнено, и FALSE там, где условие не выполнено. При сложении командой sum() каждое TRUE будет посчитано как единица, а каждое FALSE как ноль. Поэтому sum(success) даст количество раз, когда условие x_1 + x_2 + x_3^2 > 5 выполнено.
С любым распределением [xxx] в R связано четыре функции: r[xxx], d[xxx], p[xxx] и q[xxx]. Для примера возьмём нормальное распределение $N(2;9)$:
- Функция для создания случайной выборки из нормального $N(2; 9)$ распределения —
rnorm:
x <- rnorm(100, mean = 2, sd = 3) # случайная выборка из 100 нормальных N(2; 9) величин
head(x, 10) # первые 10 элементов вектора
- Функция плотности —
dnorm:
x <- seq(from = -6, to = 10, length.out = 100)
y <- dnorm(x, mean = 2, sd = 3)
qplot(x = x, y = y, geom = "line") + xlab("Значения случайной величины $X$") + ylab("Функция плотности")
Для дискретных распределений с буквы d начинается название функции, возвращающей вероятность получить заданное значение. Найдём для случайной величины $W$, имеющей Пуассоновское распределение с параметром $\lambda = 2$ вероятность $P(W = 3)$:
dpois(3, lambda = 2)
- Функция распределения, $F(t)=P(X\leq t)$ —
pnorm:
x <- seq(from = -6, to = 10, length.out = 100)
y <- pnorm(x, mean = 2, sd = 3)
qplot(x = x, y = y, geom = "line") + xlab("Значения случайной величины $X$") + ylab("Функция распределения")
- Квантильная функция или обратная функция распределения, $q(x) = F^{-1}(x)$ —
qnorm:
Найдём перцентили 5\%, 15\% и 90\% для нормального $N(2; 9)$ распределения:
x <- c(0.05, 0.15, 0.9)
qnorm(x, mean = 2, sd = 3)
Иногда бывает полезно получить случайную выборку из заданного вектора без повторений:
sample(1:100, size = 20)
Или с повторениями:
sample(c("Орёл", "Решка"), size = 10, replace = TRUE)
Можно добавить неравные вероятности:
sample(c("Орёл", "Решка"), size = 10, replace = TRUE, prob = c(0.3, 0.7))
Если выполнить команду rnorm(10, mean = 2, sd = 3) на двух разных компьютерах или два раза на одном и том же, то результат будет разный. Не зря же они случайные :) Однако генерирование случайных величин никак не противоречит идее абсолютно точной воспроизводимости исследований. Для того, чтобы получились одинаковые результаты, необходимо синхронизировать генераторы случайных чисел на этих двух компьютерах. Делается это путём задания зерна генератора случайных чисел (seed). Зерно также называют стартовым значением. В качестве зерна подойдёт любое целое число.
И в результате запуска кода
set.seed(42)
rnorm(1, mean = 2, sd = 3)
все компьютеры выведут число $6.112875$.
```{block, type="warning"}
Если код содержит генерерирование случайных чисел, то необходимо задавать зерно генератора случайных чисел!
## Базовые статистические тесты {#statistical_tests}
...
## Множественная регрессия {#classic_regression}
...
Эконометристы любят копаться в остатках :)
[https://www.r-bloggers.com/visualising-residuals/](https://www.r-bloggers.com/visualising-residuals/)
## Квантильная регрессия {#quantile_regression}
Незаслуженно забытой оказывается квантильная регрессия. Коэнкер (ссылка...) утверждает, что развитие эконометрики началось именно с квантильной регрессии. Для оценок квантильной регрессии не существует формул в явном виде, поэтому она проиграла классической регрессии среднего с готовой формулой $\hat\beta = (X'X)^{-1}X'y$. Сейчас компьютер позволяет начихать на отсутствие явных формул :)
...
## Инструментальные переменные {#instrumental_variables}
Каждый исследователь мечтает обнаружить не просто статистическую связь, а причинно-следственную. К сожалению, это не так просто. Записывая количество людей с зонтиками на улице и количество осадков но не зная физическую природу дождя, невозможно определить, вызывают ли люди с зонтиками дождь или наоборот. Для выяснения причинно-следственных связей необходим случайный эксперимент. Например, можно выбрать несколько случайных дней в году и выгнать толпу знакомых с зонтами на улицу, а затем посмотреть, было ли больше осадков в эти дни.
Использование инструментальных переменных не гарантирует нахождение причинно-следственной связи! Однако они могут быть полезны.
```{block, type="warning"}
Для обнаружения причинно-следственной связи необходимо либо использовать данные случайного эксперимента, либо прекрасно разбираться в закономерностях происходящего.
Этот раздел не о том, как обнаружить причинно-следственную связь, а о том, как реализовать метод инструментальных переменных в R и как его проинтепретировать.
Гетероскедастичность {#heteroskedasticity}
Подключаем нужные пакеты:
library("ggplot2") # графики
library("sandwich") # оценка Var для гетероскедастичности
library("lmtest") # тест Бройша-Пагана
library("dplyr") # манипуляции с данными
Условная гетероскедастичность — $Var(\epsilon_i | x_i)\neq const$.
Последствия
- Нарушены предпосылки теоремы Гаусса-Маркова. Получаемые $\hat{\beta}_{ols}$ не являются эффективными.
- Несмещенность оценок $\hat{\beta}_{ols}$ сохраняется
- Нарушены предпосылки для использования $t$-статистик. При использовании стандартных формул $t$-статистики не имеют $t$-распределения. Доверительные интервалы и проверка гипотез по стандартным формулам даёт неверные результаты.
Что делать?
- Если нужны несмещенные оценки — ничего
- Если нужно проверять гипотезы или строить доверительные интервалы, то нужно использовать правильную формулу для $\widehat{Var}(\hat{\beta})$.
- Если нужны эффективные оценки и есть представление о том, какого вида может быть гетероскедастичность то можно взвесить наблюдения. Оценить регрессию:
[
\frac{y_i}{\hat{\sigma}_i}=\beta_1 \frac{1}{\hat{\sigma}_i}+\beta_2 \frac{x_i}{\hat{\sigma}_i}+\frac{\epsilon_i}{\hat{\sigma}_i}
]
- Задать другой вопрос:
- Вместо регрессии $y_i=\beta_1+\beta_2 x_i+\epsilon_i$ попробовать оценить регрессию $\ln(y_i)=\beta_1+\beta_2 \ln(x_i)+\epsilon_i$.
- Оценить квантильную регрессию
Файл с данными по стоимости квартир в Москве доступен по ссылке goo.gl/zL5JQ
filename <- "datasets/flats_moscow.txt"
flats <- read.table(filename, header = TRUE)
ggplot(flats) + geom_point(aes(x = totsp, y = price)) +
labs(x = "Общая площадь квартиры, кв.м.", y = "Цена квартиры, тысяч у.е.",
title = "Стоимость квартир в Москве")
Строим регрессию стоимости на общую площадь.
m1 <- lm(price ~ totsp, data = flats)
summary(m1)
Покажем доверительные интервалы для коэффициентов из регрессии.
confint(m1)
Посмотрим на $\widehat{Var}(\hat{\beta}|X)$.
vcov(m1)
Если не хотим смотреть на всё summary, а только на то, что касается бет с крышкой.
coeftest(m1)
Как было сказано ранее: если в данных присутствует гетероскедастичность, но нужно проверить гипотезы и строить доверительные интервалы, то нужно использовать правильную оценку ковариационной матрицы, она выглядит так:
[
\widehat{Var}_{HC}(\hat{\beta}|X)=(X'X)^{-1}X'\hat{\Omega}X(X'X)^{-1}
]
[
\hat{\Omega}=diag(\hat{\sigma}^2_1,\ldots,\hat{\sigma}^2_n)
]
Есть разные способы состоятельно оценить отдельные дисперсии $\hat{\sigma}^2_i$, подробно можно почитать в описании пакета sandwich. Наиболее популярный на сегодня — это так называемый "HC3":
[
\hat{\sigma}^2_i=\frac{\hat{\epsilon}i^2}{(1-h{ii})^2}
]
Здесь $h_{ii}$ — диагональный элемент матрицы-шляпницы (hat matrix), $\hat{y}=Hy$, $H=X(X'X)^{-1}X'$. Матрицу-шляпницу также называют проектором, т.к. она проецирует вектор $y$ на линейную оболочку регрессоров.
Посмотрим на матрицу $\widehat{Var}_{HC}(\hat{\beta}|X)$ в R.
vcovHC(m1)
Применяя правильную $\widehat{Var}_{HC}(\hat{\beta}|X)$, проверим гипотезы.
coeftest(m1, vcov = vcovHC(m1))
Доверительный интервал надо строить руками. За основу мы возьмем табличку с правильными стандартными ошибками и прибавим и вычтем их от оценок коэффициентов.
conftable <- coeftest(m1, vcov = vcovHC(m1))
ci <- data_frame(estimate = conftable[, 1],
se_hc = conftable[, 2],
left_95 = estimate - qnorm(0.975)*se_hc,
right_95 = estimate + qnorm(0.975)*se_hc)
ci
Ещё есть малоизвестный пакет Greg Макса Гордона и в нём есть готовая функция:
Greg::confint_robust(m1)
Пакет Greg не доступен на официальном репозитории CRAN, ставится с github командой devtools::install_github("gforge/Greg"). Пакет Greg тянет за собой такую кучу всего, что не хочется его ставить ради одной функции. Здесь уместно упражнение: сделайте функцию confint_robust самостоятельно :)
В условиях гетероскедастичности обычная $t$-статистика не имеет ни $t$-распределения при нормальных $\varepsilon_i$, ни даже асимптотически нормального. И такая же проблема возникает с $F$-статистикой для сравнения вложенных моделей. Отныне привычная $F$-статистика
[
F=\frac{(RSS_R-RSS_{UR})/r}{RSS_{UR}/(n-k_{UR})}
]
не имеет ни $F$-распределения при нормальных остатках, ни $\chi^2$-распределения асимптотически. В частности, для проверки незначимости регрессии в целом некорректно использовать $F$-статистику равную $F=\frac{ESS/(k-1)}{RSS/(n-k)}$.
Если для теста Вальда использовать корректную оценку ковариационной матрицы, то его можно применять в условиях гетероскедастичности. Допустим, мы хотим проверить гипотезу $H_0$: $A\beta=b$, состоящую из $q$ уравнений, против альтернативной о том, что хотя бы одно равенство нарушено. Тогда статистика Вальда имеет вид:
[
W = (A\beta - b)'(A\cdot \widehat{Var}_{HC}(\hat\beta)A')^{-1}(A\beta - b)
]
и имеет асимптотически $\chi^2$-распределение с $q$ степенями свободы.
В R реализуется с помощью опции для команды waldtest. Покажем на примере гипотезы о незначимости регрессии в целом:
m0 <- lm(price ~ 1, data = flats) # оцениваем ограниченную модель (регрессия на константу)
waldtest(m0, m1, vcov = vcovHC(m1))
Переходим ко взвешенному МНК. Взвешенный МНК требует знания структуры гетероскедастичности, что редко встречается на практике. Предположим, что в нашем случае $Var(\epsilon_i | totsp_i) = const \cdot totsp_i$.
В R вектор весов weights должен быть обратно пропорционален дисперсиям, т.е. $w_i=1/\sigma^2_i$.
m2 <- lm(price ~ totsp, weights = I(1 / totsp), data = flats)
summary(m2)
Тесты на гетероскедастичность {#hetero_testing}
Графическое обнаружение гетероскедастичности:
- Оценивается исходная регрессия и из неё получаются остатки $\hat{\varepsilon}_i$.
- Если верить в гомоскедастичность, то дисперсия вектора остатков определяется по формуле $Var(\hat{\varepsilon}|X)=\sigma^2 (I-X(X'X)^{-1}X')$. Т.е. дисперсия у остатков разная (!) даже если $Var(\varepsilon_i|X)=\sigma^2$. Поэтому остатки стандартизуют по формуле:
[
s_i = \frac{\varepsilon_i}{\sqrt{\hat{\sigma}^2(1-h_{ii})}}
]
Остатки полученные после такого преобразования называют стандартизированными или иногда стьюдентизированными.
- Можно построить график зависимости величины $s_i^2$ или $|s_i|$ от переменной, предположительно повинной в гетроскедастичности.
В R:
m1.st.resid <- rstandard(m1) # получаем стандартизированные остатки
ggplot(aes(x = totsp, y = abs(m1.st.resid)), data =flats) + geom_point(alpha = 0.2) +
labs(x = "Общая площадь квартиры, кв.м", y = expression(paste("Стандартизированные остатки, ", s[i])),
title = "Графическое обнаружение гетероскедастичности")
Иногда для простоты пропускают второй шаг со стандартизацией остатков и в результате могут ошибочно принять принять разную $Var(\hat{\varepsilon}_i)$ за гетероскедастичность, т.е. за разную $Var(\varepsilon_i)$. Впрочем, если гетероскедастичность сильная, то её будет видно и на нестандартизированных остатках.
Тест Бройша-Пагана, Breusch-Pagan: {#hetero_bp_test}
Есть две версии теста Бройша-Пагана: авторская и современная модификация.
Предпосылки авторской версии:
- нормальность остатков, $\varepsilon_i \sim N(0,\sigma^2_i)$
- $\sigma^2_i=h(\alpha_1+\alpha_2 z_{i2}+\ldots+\alpha_{p}z_{ip})$
- у функции $h$ существуют первая и вторая производные
- тест асимптотический
Суть теста: Используя метод максимального правдоподобия посчитаем LM-статистику. При верной $H_0$ она имеет хи-квадрат распределение с $p-1$ степенью свободы.
Оказывается, что LM-статистику можно получить с помощью вспомогательной регрессии. Авторская процедура:
- Оценивается регрессия $y_i=\beta_1+\beta_2 x_i +\varepsilon_i$
- Переходим к $g_i=\frac{n}{SSR}\hat{\varepsilon_i}$
- Строим регрессию $g_i$ на $\alpha_1+\alpha_2 z_{i2}+\ldots+\alpha_{p}z_{ip}$
- $LM=\frac{ESS}{2}$
Современная модификация выглядит (неизвестный рецензент Коэнкера) так:
- Оценивается регрессия $y_i=\beta_1+\beta_2 x_i +\varepsilon_i$
- Оценивается регрессия $\hat{\varepsilon_i^2}$ на переменные вызывающие гетероскедастичность
- При верной $H_0$ асимптотически:
[
nR^2 \sim \chi^2_{p-1},
]
где $p$ — число оцениваемых коэффициентов во вспомогательной регрессии. По смыслу $(p-1)$ — это количество факторов, от которых потенциально может зависеть дисперсия $Var(\varepsilon_i)$.
Тест Бройша-Пагана в R:
bptest(m1) # версия Коэнкера
bptest(m1, studentize = FALSE) # классика Бройша-Пагана
У современной модификации Бройша-Пагана более слабые предпосылки:
- $H_0$: $Var(\varepsilon_i)=\sigma^2$, нормальность не требуется
- $E(\varepsilon_i^4)=const$
- тест асимптотический
- (? Коэнкер в статье критикует что-то про мощность ?)
Тест Уайта. White test {#white_test}
Предпосылки:
- $H_0$: $Var(\varepsilon_i)=\sigma^2$, нормальность не требуется
- $E(\varepsilon_i^4)=const$
- тест асимптотический
Процедура:
- Оценивается регрессия $y_i=\beta_1+\beta_2 x_i +\varepsilon_i$
- Строим регрессию $\hat{\varepsilon_i^2}$ на исходные регрессоры и их квадраты
- Асимптотически $nR^2$ имеет хи-квадрат распределение
Тест Уайта в R:
bptest(m1, varformula = ~ totsp + I(totsp^2), data = flats)
Тест Уайта является частным случаем современной модификации теста Бройша-Пагана. В тесте Бройша-Пагана во вспомогательной регресии можно брать любые объясняющие переменные. В тесте Уайта берутся исходные регрессоры, их квадраты и попарные произведения. Если в R попросить сделать тест Бройша-Пагана и не указать спецификацию вспомогательной регрессии, то он возьмет только все регрессоры исходной.
Тест Goldfeld-Quandt {#gq_test}
Предпосылки:
- нормальность остатков, $\varepsilon_i \sim N(0,\sigma^2_i)$
- Наблюдения упорядочены по возростанию гетероскедастичности
- тест точный (неасимптотический)
Процедура:
- Упорядочить наблюдения в том порядке, в котором подозревается рост гетероскедастичности
- Выкинуть некий процент наблюдений по середине, чтобы подчеркнуть разницу в дисперсии между начальными и конечными наблюдениями.
- Оценить исходную модель по наблюдениям из начала выборки и по наблюдениям конца выборки
- Получить, соответственно, $RSS_1$ и $RSS_2$
При верной $H_0$
[
\frac{RSS_2}{RSS_1}\sim F_{r-k,r-k},
]
где $r$ — размер подвыборки в начале и в конце.
Тест Голдфельда-Квандта в R:
flats2 <- flats[order(flats$totsp), ] # сменим порядок строк в табличке h
m3 <- lm(price ~ totsp, data = flats2)
gqtest(m3, fraction = 0.2) # проведем GQ тест выкинув посередине 20% наблюдений
Тестирование и борьба с гетероскедастичностью {#hetero_fail}
Чайники часто используют такой ошибочный подход:
Протестировать гетероскедастичность. Если тесты выявили гетероскедастичность, то бороться с ней (использовать устойчивые стандартные ошибки или применять взвешенный мнк). Если тесты не выявили гетероскедастичность, то не бороться с ней.
Этот подход неверен!
Правильно делать так:
- Если есть теоретические основания ожидать гетероскедастичность и наблюдений достаточно, то использовать устойчивые стандартные ошибки вне зависимости от результатов тестирования.
- Если есть теоретические основания ожидать гетероскедастичность и наблюдений мало, то отказаться от девичьих грёз об эффективности оценок и проверке гипотез. Довольствоваться возможной несмещенностью оценок.
Тут обычно спрашивают: "А зачем же тогда тестировать на гетероскедастичность, если это решение ни на что не влияет?". Ответ прост — чтобы узнать, есть ли гетероскедастичность. :) Впрочем, если есть теоретические основания её ожидать, то можно и не тестировать.
Прогнозирование при гетероскедастичности {#hetero_forecasting}
При прогнозировании строят доверительный интервал для среднего значения зависимой переменной и предиктивный интервал для конкретного будущего значения. С доверительным интервалом похоже надо код писать руками:
# ...
С предиктивным интервалом чуть сложнее. Поскольку предиктивный интервал строиться для конкретного наблюдения, то нам надо знать дисперсию этого наблюдения. Идея робастных ошибок состоит в том, чтобы проверять гипотезы не делая точных предположений о структуре гетероскедастичности. Поэтому при использовании робастных стандартных ошибок традиционно предиктивные интервалы не строят.
Хм. А если взять парадигму рандомных регрессоров и оценку безусловной дисперсии?
Если же случилось чудо и структура гетероскедастичности доподлинно известна, то мы спокойно используем взвешенный МНК и указываем вес для новых наблюдений.
# ...
Дополнительно: {#hetero_links}
- Описание пакета
sandwich с очень правильным изложением современного взгляда на гетероскедастичность
- Breusch, Pagan, Simple test for heteroskedasticity
- White, A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity
- Koenker, A note on studentizing a test for heteroscedasticity
Работа с качественными переменными {#factor_variables}
...
Логит и пробит с визуализацией {#logit_probit}
...
Метод главных компонент {#pca}
...
Мультиколлинеарность {#multicollinearity}
...
LASSO {#lasso}
Оказывается, что оценки LASSO для всех значений параметра $\lambda$ можно получить очень быстро! Для этого нам потребуется алгоритм LARS (Least Angle Regresssion):
-
Положим стартовые значения всех $\beta_j$ равными нулю.
-
Найдём регрессор $x_j$ сильнее всего коррелированный с $y$ по модулю
-
Изменяем коэффициент $\beta_j$ в направлении знака корреляции $x_j$ и $y_j$. Попутно считаем остатки $\hat u = y - \hat y$. Прекращаем изменять коэффициент $\beta_j$, как только некоторый регрессор $x_k$ будет настолько же сильно коррелирован с $r$ как $x_j$
-
Увеличиваем пару $(\beta_j, \beta_k)$ в направлении задаваемом линейной регрессией, до тех пор пока у некоторого $x_m$ не будет высокой корреляции с $r$
-
Если ненулевой коэффициент уходит в ноль, то выкидываем его из множества активных регрессоров.
-
Продолжаем до тех пор, пока все регрессоры не вошли в модель.
Метод максимального правдоподобия {#maximum_likelihood}
...
Метод опорных векторов {#svm}
...
Случайный лес {#random_forest}
...
Экспоненциальное сглаживание {#exponential_smoothing}
...
ARMA модели {#arma}
...
GARCH {#garch}
...
VAR и BVAR {#var_bvar}
...
Я не верю в пользу от структурных BVAR, поэтому их здесь нет :)
Панельные данные {#panel_data}
...
Байесовский подход: первые шаги {#bayesian_first_steps}
...
Байесовский подход: STAN {#stan}
...
Карты {#maps}
Где карта, Билли?
Дифференциальные уравнения {#differential_equations}
...
Задачи оптимизации {#optimisation}
...
bdemeshev/r_manual_book documentation built on May 12, 2019, 3:41 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(knitr) library(tikzDevice) activate_tikz <- function() { # tikz plots options options(tikzDefaultEngine = "xetex") # cash font metrics for speed: # options(tikzMetricsDictionary = "./tikz_metrics") add_xelatex <- c("\\defaultfontfeatures{Ligatures=TeX, Scale=MatchLowercase}", "\\setmainfont{Linux Libertine O}", "\\setmonofont{Linux Libertine O}", "\\setsansfont{Linux Libertine O}", "\\newfontfamily{\\cyrillicfonttt}{Linux Libertine O}", "\\newfontfamily{\\cyrillicfont}{Linux Libertine O}", "\\newfontfamily{\\cyrillicfontsf}{Linux Libertine O}") options(tikzXelatexPackages = c(getOption("tikzXelatexPackages"), add_xelatex)) # does remove warnings: # it is important to remove fontenc package wich is loaded by default options(tikzUnicodeMetricPackages = c("\\usetikzlibrary{calc}", "\\usepackage{fontspec, xunicode}", add_xelatex)) opts_chunk$set(dev = "tikz", dev.args = list(pointsize = 11)) } activate_png <- function() { opts_chunk$set(dev = "png", dpi = 300) } colFmt <- function(x, color) { outputFormat <- opts_knit$get("rmarkdown.pandoc.to") if (outputFormat == "latex") { result <- paste0("\\textcolor{", color, "}{", x, "}") } else if (outputFormat %in% c("html", "epub")) { result <- paste0("<font color='", color, "'>", x, "</font>") } else { result <- x } return(result) } outputFormat <- opts_knit$get("rmarkdown.pandoc.to") if (outputFormat == "latex") { # activateTikz() activate_png() # другую тему для ggplot2 выставить? }
library(tidyverse) # рабочие кони library(sandwich) # оценка Var для гетероскедастичности library(lmtest) # тест Бройша-Пагана library(data.table) # манипуляции с данными library(reshape2) # преобразование длинных таблиц в широкие library(broom)
Статистика и не только {#statistics_and_more}
Пора вспомнить о том, что R ориентирован на статистику :)
Генерирование случайных величин {#random_variables}
Для решения задач по теории вероятностей или исследования свойств статистических алгоритмов может потребоваться сгененировать случайную выборку из заданного закона распределения.
Генерируем 10 равномерных на отрезке $[4;10.5]$ случайных величин:
runif(10, min = 4, max = 10.5)
Генерируем 10 нормальных $N(2;9)$ случайных величин с математическим ожиданием $2$ и дисперсией $9=3^2$:
rnorm(10, mean = 2, sd = 3)
Например, с помощью симуляций легко оценить математическое ожидание $E(1/X)$, где $X \sim N(2;9)$. Для этого мы вспомним Закон Больших Чисел. Он говорит, что арифметическое среднее по большой выборке стремится по вероятности и почти наверное к математическому ожиданию. Поэтому мы просто сгенерируем большую выборку в миллион наблюдений:
n_obs <- 10^6 x <- rnorm(n_obs, mean = 2, sd = 3) mean(1/x)
Также легко оценить многие вероятности. Например, оценим вероятность $P(X_1 + X_2 + X_3^2 > 5)$, где величины $X_i$ независимы и одинаково распределены $X_i \sim U[0;2]$:
n_obs <- 10^6 x_1 <- runif(n_obs, min = 0, max = 2) x_2 <- runif(n_obs, min = 0, max = 2) x_3 <- runif(n_obs, min = 0, max = 2) success <- x_1 + x_2 + x_3^2 > 5 sum(success) / n_obs
Здесь вектор success будет содержать значение TRUE там, где условие x_1 + x_2 + x_3^2 > 5 выполнено, и FALSE там, где условие не выполнено. При сложении командой sum() каждое TRUE будет посчитано как единица, а каждое FALSE как ноль. Поэтому sum(success) даст количество раз, когда условие x_1 + x_2 + x_3^2 > 5 выполнено.
С любым распределением [xxx] в R связано четыре функции: r[xxx], d[xxx], p[xxx] и q[xxx]. Для примера возьмём нормальное распределение $N(2;9)$:
- Функция для создания случайной выборки из нормального $N(2; 9)$ распределения —
rnorm:
x <- rnorm(100, mean = 2, sd = 3) # случайная выборка из 100 нормальных N(2; 9) величин head(x, 10) # первые 10 элементов вектора
- Функция плотности —
dnorm:
x <- seq(from = -6, to = 10, length.out = 100) y <- dnorm(x, mean = 2, sd = 3) qplot(x = x, y = y, geom = "line") + xlab("Значения случайной величины $X$") + ylab("Функция плотности")
Для дискретных распределений с буквы d начинается название функции, возвращающей вероятность получить заданное значение. Найдём для случайной величины $W$, имеющей Пуассоновское распределение с параметром $\lambda = 2$ вероятность $P(W = 3)$:
dpois(3, lambda = 2)
- Функция распределения, $F(t)=P(X\leq t)$ —
pnorm:
x <- seq(from = -6, to = 10, length.out = 100) y <- pnorm(x, mean = 2, sd = 3) qplot(x = x, y = y, geom = "line") + xlab("Значения случайной величины $X$") + ylab("Функция распределения")
- Квантильная функция или обратная функция распределения, $q(x) = F^{-1}(x)$ —
qnorm:
Найдём перцентили 5\%, 15\% и 90\% для нормального $N(2; 9)$ распределения:
x <- c(0.05, 0.15, 0.9) qnorm(x, mean = 2, sd = 3)
Иногда бывает полезно получить случайную выборку из заданного вектора без повторений:
sample(1:100, size = 20)
Или с повторениями:
sample(c("Орёл", "Решка"), size = 10, replace = TRUE)
Можно добавить неравные вероятности:
sample(c("Орёл", "Решка"), size = 10, replace = TRUE, prob = c(0.3, 0.7))
Если выполнить команду rnorm(10, mean = 2, sd = 3) на двух разных компьютерах или два раза на одном и том же, то результат будет разный. Не зря же они случайные :) Однако генерирование случайных величин никак не противоречит идее абсолютно точной воспроизводимости исследований. Для того, чтобы получились одинаковые результаты, необходимо синхронизировать генераторы случайных чисел на этих двух компьютерах. Делается это путём задания зерна генератора случайных чисел (seed). Зерно также называют стартовым значением. В качестве зерна подойдёт любое целое число.
И в результате запуска кода
set.seed(42) rnorm(1, mean = 2, sd = 3)
все компьютеры выведут число $6.112875$.
```{block, type="warning"} Если код содержит генерерирование случайных чисел, то необходимо задавать зерно генератора случайных чисел!
## Базовые статистические тесты {#statistical_tests}
...
## Множественная регрессия {#classic_regression}
...
Эконометристы любят копаться в остатках :)
[https://www.r-bloggers.com/visualising-residuals/](https://www.r-bloggers.com/visualising-residuals/)
## Квантильная регрессия {#quantile_regression}
Незаслуженно забытой оказывается квантильная регрессия. Коэнкер (ссылка...) утверждает, что развитие эконометрики началось именно с квантильной регрессии. Для оценок квантильной регрессии не существует формул в явном виде, поэтому она проиграла классической регрессии среднего с готовой формулой $\hat\beta = (X'X)^{-1}X'y$. Сейчас компьютер позволяет начихать на отсутствие явных формул :)
...
## Инструментальные переменные {#instrumental_variables}
Каждый исследователь мечтает обнаружить не просто статистическую связь, а причинно-следственную. К сожалению, это не так просто. Записывая количество людей с зонтиками на улице и количество осадков но не зная физическую природу дождя, невозможно определить, вызывают ли люди с зонтиками дождь или наоборот. Для выяснения причинно-следственных связей необходим случайный эксперимент. Например, можно выбрать несколько случайных дней в году и выгнать толпу знакомых с зонтами на улицу, а затем посмотреть, было ли больше осадков в эти дни.
Использование инструментальных переменных не гарантирует нахождение причинно-следственной связи! Однако они могут быть полезны.
```{block, type="warning"}
Для обнаружения причинно-следственной связи необходимо либо использовать данные случайного эксперимента, либо прекрасно разбираться в закономерностях происходящего.
Этот раздел не о том, как обнаружить причинно-следственную связь, а о том, как реализовать метод инструментальных переменных в R и как его проинтепретировать.
Гетероскедастичность {#heteroskedasticity}
Подключаем нужные пакеты:
library("ggplot2") # графики library("sandwich") # оценка Var для гетероскедастичности library("lmtest") # тест Бройша-Пагана library("dplyr") # манипуляции с данными
Условная гетероскедастичность — $Var(\epsilon_i | x_i)\neq const$.
Последствия
- Нарушены предпосылки теоремы Гаусса-Маркова. Получаемые $\hat{\beta}_{ols}$ не являются эффективными.
- Несмещенность оценок $\hat{\beta}_{ols}$ сохраняется
- Нарушены предпосылки для использования $t$-статистик. При использовании стандартных формул $t$-статистики не имеют $t$-распределения. Доверительные интервалы и проверка гипотез по стандартным формулам даёт неверные результаты.
Что делать?
- Если нужны несмещенные оценки — ничего
- Если нужно проверять гипотезы или строить доверительные интервалы, то нужно использовать правильную формулу для $\widehat{Var}(\hat{\beta})$.
- Если нужны эффективные оценки и есть представление о том, какого вида может быть гетероскедастичность то можно взвесить наблюдения. Оценить регрессию:
[ \frac{y_i}{\hat{\sigma}_i}=\beta_1 \frac{1}{\hat{\sigma}_i}+\beta_2 \frac{x_i}{\hat{\sigma}_i}+\frac{\epsilon_i}{\hat{\sigma}_i} ]
- Задать другой вопрос:
- Вместо регрессии $y_i=\beta_1+\beta_2 x_i+\epsilon_i$ попробовать оценить регрессию $\ln(y_i)=\beta_1+\beta_2 \ln(x_i)+\epsilon_i$.
- Оценить квантильную регрессию
Файл с данными по стоимости квартир в Москве доступен по ссылке goo.gl/zL5JQ
filename <- "datasets/flats_moscow.txt" flats <- read.table(filename, header = TRUE) ggplot(flats) + geom_point(aes(x = totsp, y = price)) + labs(x = "Общая площадь квартиры, кв.м.", y = "Цена квартиры, тысяч у.е.", title = "Стоимость квартир в Москве")
Строим регрессию стоимости на общую площадь.
m1 <- lm(price ~ totsp, data = flats) summary(m1)
Покажем доверительные интервалы для коэффициентов из регрессии.
confint(m1)
Посмотрим на $\widehat{Var}(\hat{\beta}|X)$.
vcov(m1)
Если не хотим смотреть на всё summary, а только на то, что касается бет с крышкой.
coeftest(m1)
Как было сказано ранее: если в данных присутствует гетероскедастичность, но нужно проверить гипотезы и строить доверительные интервалы, то нужно использовать правильную оценку ковариационной матрицы, она выглядит так: [ \widehat{Var}_{HC}(\hat{\beta}|X)=(X'X)^{-1}X'\hat{\Omega}X(X'X)^{-1} ]
[ \hat{\Omega}=diag(\hat{\sigma}^2_1,\ldots,\hat{\sigma}^2_n) ]
Есть разные способы состоятельно оценить отдельные дисперсии $\hat{\sigma}^2_i$, подробно можно почитать в описании пакета sandwich. Наиболее популярный на сегодня — это так называемый "HC3":
[
\hat{\sigma}^2_i=\frac{\hat{\epsilon}i^2}{(1-h{ii})^2}
]
Здесь $h_{ii}$ — диагональный элемент матрицы-шляпницы (hat matrix), $\hat{y}=Hy$, $H=X(X'X)^{-1}X'$. Матрицу-шляпницу также называют проектором, т.к. она проецирует вектор $y$ на линейную оболочку регрессоров.
Посмотрим на матрицу $\widehat{Var}_{HC}(\hat{\beta}|X)$ в R.
vcovHC(m1)
Применяя правильную $\widehat{Var}_{HC}(\hat{\beta}|X)$, проверим гипотезы.
coeftest(m1, vcov = vcovHC(m1))
Доверительный интервал надо строить руками. За основу мы возьмем табличку с правильными стандартными ошибками и прибавим и вычтем их от оценок коэффициентов.
conftable <- coeftest(m1, vcov = vcovHC(m1)) ci <- data_frame(estimate = conftable[, 1], se_hc = conftable[, 2], left_95 = estimate - qnorm(0.975)*se_hc, right_95 = estimate + qnorm(0.975)*se_hc) ci
Ещё есть малоизвестный пакет Greg Макса Гордона и в нём есть готовая функция:
Greg::confint_robust(m1)
Пакет Greg не доступен на официальном репозитории CRAN, ставится с github командой devtools::install_github("gforge/Greg"). Пакет Greg тянет за собой такую кучу всего, что не хочется его ставить ради одной функции. Здесь уместно упражнение: сделайте функцию confint_robust самостоятельно :)
В условиях гетероскедастичности обычная $t$-статистика не имеет ни $t$-распределения при нормальных $\varepsilon_i$, ни даже асимптотически нормального. И такая же проблема возникает с $F$-статистикой для сравнения вложенных моделей. Отныне привычная $F$-статистика
[ F=\frac{(RSS_R-RSS_{UR})/r}{RSS_{UR}/(n-k_{UR})} ]
не имеет ни $F$-распределения при нормальных остатках, ни $\chi^2$-распределения асимптотически. В частности, для проверки незначимости регрессии в целом некорректно использовать $F$-статистику равную $F=\frac{ESS/(k-1)}{RSS/(n-k)}$.
Если для теста Вальда использовать корректную оценку ковариационной матрицы, то его можно применять в условиях гетероскедастичности. Допустим, мы хотим проверить гипотезу $H_0$: $A\beta=b$, состоящую из $q$ уравнений, против альтернативной о том, что хотя бы одно равенство нарушено. Тогда статистика Вальда имеет вид:
[ W = (A\beta - b)'(A\cdot \widehat{Var}_{HC}(\hat\beta)A')^{-1}(A\beta - b) ]
и имеет асимптотически $\chi^2$-распределение с $q$ степенями свободы.
В R реализуется с помощью опции для команды waldtest. Покажем на примере гипотезы о незначимости регрессии в целом:
m0 <- lm(price ~ 1, data = flats) # оцениваем ограниченную модель (регрессия на константу) waldtest(m0, m1, vcov = vcovHC(m1))
Переходим ко взвешенному МНК. Взвешенный МНК требует знания структуры гетероскедастичности, что редко встречается на практике. Предположим, что в нашем случае $Var(\epsilon_i | totsp_i) = const \cdot totsp_i$.
В R вектор весов weights должен быть обратно пропорционален дисперсиям, т.е. $w_i=1/\sigma^2_i$.
m2 <- lm(price ~ totsp, weights = I(1 / totsp), data = flats) summary(m2)
Тесты на гетероскедастичность {#hetero_testing}
Графическое обнаружение гетероскедастичности:
- Оценивается исходная регрессия и из неё получаются остатки $\hat{\varepsilon}_i$.
- Если верить в гомоскедастичность, то дисперсия вектора остатков определяется по формуле $Var(\hat{\varepsilon}|X)=\sigma^2 (I-X(X'X)^{-1}X')$. Т.е. дисперсия у остатков разная (!) даже если $Var(\varepsilon_i|X)=\sigma^2$. Поэтому остатки стандартизуют по формуле: [ s_i = \frac{\varepsilon_i}{\sqrt{\hat{\sigma}^2(1-h_{ii})}} ] Остатки полученные после такого преобразования называют стандартизированными или иногда стьюдентизированными.
- Можно построить график зависимости величины $s_i^2$ или $|s_i|$ от переменной, предположительно повинной в гетроскедастичности.
В R:
m1.st.resid <- rstandard(m1) # получаем стандартизированные остатки ggplot(aes(x = totsp, y = abs(m1.st.resid)), data =flats) + geom_point(alpha = 0.2) + labs(x = "Общая площадь квартиры, кв.м", y = expression(paste("Стандартизированные остатки, ", s[i])), title = "Графическое обнаружение гетероскедастичности")
Иногда для простоты пропускают второй шаг со стандартизацией остатков и в результате могут ошибочно принять принять разную $Var(\hat{\varepsilon}_i)$ за гетероскедастичность, т.е. за разную $Var(\varepsilon_i)$. Впрочем, если гетероскедастичность сильная, то её будет видно и на нестандартизированных остатках.
Тест Бройша-Пагана, Breusch-Pagan: {#hetero_bp_test}
Есть две версии теста Бройша-Пагана: авторская и современная модификация.
Предпосылки авторской версии:
- нормальность остатков, $\varepsilon_i \sim N(0,\sigma^2_i)$
- $\sigma^2_i=h(\alpha_1+\alpha_2 z_{i2}+\ldots+\alpha_{p}z_{ip})$
- у функции $h$ существуют первая и вторая производные
- тест асимптотический
Суть теста: Используя метод максимального правдоподобия посчитаем LM-статистику. При верной $H_0$ она имеет хи-квадрат распределение с $p-1$ степенью свободы.
Оказывается, что LM-статистику можно получить с помощью вспомогательной регрессии. Авторская процедура:
- Оценивается регрессия $y_i=\beta_1+\beta_2 x_i +\varepsilon_i$
- Переходим к $g_i=\frac{n}{SSR}\hat{\varepsilon_i}$
- Строим регрессию $g_i$ на $\alpha_1+\alpha_2 z_{i2}+\ldots+\alpha_{p}z_{ip}$
- $LM=\frac{ESS}{2}$
Современная модификация выглядит (неизвестный рецензент Коэнкера) так:
- Оценивается регрессия $y_i=\beta_1+\beta_2 x_i +\varepsilon_i$
- Оценивается регрессия $\hat{\varepsilon_i^2}$ на переменные вызывающие гетероскедастичность
- При верной $H_0$ асимптотически: [ nR^2 \sim \chi^2_{p-1}, ] где $p$ — число оцениваемых коэффициентов во вспомогательной регрессии. По смыслу $(p-1)$ — это количество факторов, от которых потенциально может зависеть дисперсия $Var(\varepsilon_i)$.
Тест Бройша-Пагана в R:
bptest(m1) # версия Коэнкера bptest(m1, studentize = FALSE) # классика Бройша-Пагана
У современной модификации Бройша-Пагана более слабые предпосылки:
- $H_0$: $Var(\varepsilon_i)=\sigma^2$, нормальность не требуется
- $E(\varepsilon_i^4)=const$
- тест асимптотический
- (? Коэнкер в статье критикует что-то про мощность ?)
Тест Уайта. White test {#white_test}
Предпосылки:
- $H_0$: $Var(\varepsilon_i)=\sigma^2$, нормальность не требуется
- $E(\varepsilon_i^4)=const$
- тест асимптотический
Процедура:
- Оценивается регрессия $y_i=\beta_1+\beta_2 x_i +\varepsilon_i$
- Строим регрессию $\hat{\varepsilon_i^2}$ на исходные регрессоры и их квадраты
- Асимптотически $nR^2$ имеет хи-квадрат распределение
Тест Уайта в R:
bptest(m1, varformula = ~ totsp + I(totsp^2), data = flats)
Тест Уайта является частным случаем современной модификации теста Бройша-Пагана. В тесте Бройша-Пагана во вспомогательной регресии можно брать любые объясняющие переменные. В тесте Уайта берутся исходные регрессоры, их квадраты и попарные произведения. Если в R попросить сделать тест Бройша-Пагана и не указать спецификацию вспомогательной регрессии, то он возьмет только все регрессоры исходной.
Тест Goldfeld-Quandt {#gq_test}
Предпосылки:
- нормальность остатков, $\varepsilon_i \sim N(0,\sigma^2_i)$
- Наблюдения упорядочены по возростанию гетероскедастичности
- тест точный (неасимптотический)
Процедура:
- Упорядочить наблюдения в том порядке, в котором подозревается рост гетероскедастичности
- Выкинуть некий процент наблюдений по середине, чтобы подчеркнуть разницу в дисперсии между начальными и конечными наблюдениями.
- Оценить исходную модель по наблюдениям из начала выборки и по наблюдениям конца выборки
- Получить, соответственно, $RSS_1$ и $RSS_2$
При верной $H_0$ [ \frac{RSS_2}{RSS_1}\sim F_{r-k,r-k}, ] где $r$ — размер подвыборки в начале и в конце.
Тест Голдфельда-Квандта в R:
flats2 <- flats[order(flats$totsp), ] # сменим порядок строк в табличке h m3 <- lm(price ~ totsp, data = flats2) gqtest(m3, fraction = 0.2) # проведем GQ тест выкинув посередине 20% наблюдений
Тестирование и борьба с гетероскедастичностью {#hetero_fail}
Чайники часто используют такой ошибочный подход:
Протестировать гетероскедастичность. Если тесты выявили гетероскедастичность, то бороться с ней (использовать устойчивые стандартные ошибки или применять взвешенный мнк). Если тесты не выявили гетероскедастичность, то не бороться с ней.
Этот подход неверен!
Правильно делать так:
- Если есть теоретические основания ожидать гетероскедастичность и наблюдений достаточно, то использовать устойчивые стандартные ошибки вне зависимости от результатов тестирования.
- Если есть теоретические основания ожидать гетероскедастичность и наблюдений мало, то отказаться от девичьих грёз об эффективности оценок и проверке гипотез. Довольствоваться возможной несмещенностью оценок.
Тут обычно спрашивают: "А зачем же тогда тестировать на гетероскедастичность, если это решение ни на что не влияет?". Ответ прост — чтобы узнать, есть ли гетероскедастичность. :) Впрочем, если есть теоретические основания её ожидать, то можно и не тестировать.
Прогнозирование при гетероскедастичности {#hetero_forecasting}
При прогнозировании строят доверительный интервал для среднего значения зависимой переменной и предиктивный интервал для конкретного будущего значения. С доверительным интервалом похоже надо код писать руками:
# ...
С предиктивным интервалом чуть сложнее. Поскольку предиктивный интервал строиться для конкретного наблюдения, то нам надо знать дисперсию этого наблюдения. Идея робастных ошибок состоит в том, чтобы проверять гипотезы не делая точных предположений о структуре гетероскедастичности. Поэтому при использовании робастных стандартных ошибок традиционно предиктивные интервалы не строят.
Хм. А если взять парадигму рандомных регрессоров и оценку безусловной дисперсии?
Если же случилось чудо и структура гетероскедастичности доподлинно известна, то мы спокойно используем взвешенный МНК и указываем вес для новых наблюдений.
# ...
Дополнительно: {#hetero_links}
- Описание пакета
sandwichс очень правильным изложением современного взгляда на гетероскедастичность - Breusch, Pagan, Simple test for heteroskedasticity
- White, A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroskedasticity
- Koenker, A note on studentizing a test for heteroscedasticity
Работа с качественными переменными {#factor_variables}
...
Логит и пробит с визуализацией {#logit_probit}
...
Метод главных компонент {#pca}
...
Мультиколлинеарность {#multicollinearity}
...
LASSO {#lasso}
Оказывается, что оценки LASSO для всех значений параметра $\lambda$ можно получить очень быстро! Для этого нам потребуется алгоритм LARS (Least Angle Regresssion):
-
Положим стартовые значения всех $\beta_j$ равными нулю.
-
Найдём регрессор $x_j$ сильнее всего коррелированный с $y$ по модулю
-
Изменяем коэффициент $\beta_j$ в направлении знака корреляции $x_j$ и $y_j$. Попутно считаем остатки $\hat u = y - \hat y$. Прекращаем изменять коэффициент $\beta_j$, как только некоторый регрессор $x_k$ будет настолько же сильно коррелирован с $r$ как $x_j$
-
Увеличиваем пару $(\beta_j, \beta_k)$ в направлении задаваемом линейной регрессией, до тех пор пока у некоторого $x_m$ не будет высокой корреляции с $r$
-
Если ненулевой коэффициент уходит в ноль, то выкидываем его из множества активных регрессоров.
-
Продолжаем до тех пор, пока все регрессоры не вошли в модель.
Метод максимального правдоподобия {#maximum_likelihood}
...
Метод опорных векторов {#svm}
...
Случайный лес {#random_forest}
...
Экспоненциальное сглаживание {#exponential_smoothing}
...
ARMA модели {#arma}
...
GARCH {#garch}
...
VAR и BVAR {#var_bvar}
...
Я не верю в пользу от структурных BVAR, поэтому их здесь нет :)
Панельные данные {#panel_data}
...
Байесовский подход: первые шаги {#bayesian_first_steps}
...
Байесовский подход: STAN {#stan}
...
Карты {#maps}
Где карта, Билли?
Дифференциальные уравнения {#differential_equations}
...
Задачи оптимизации {#optimisation}
...
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.