Сегодня обсудим интересную базу данных – базу с экспериментальными данными по шоколадным тортикам (к политологии вернемся на семинаре).
Описание базы данных. Детали по эксперименту.
Загрузим базу данных по ссылке:
df <- read.csv("https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/lme4/cake.csv")Посмотрим на нее:
View(df)Для вывода описательных статистик в R есть специальная функция summary():
summary(df)## X replicate recipe temperature angle
## Min. : 1.00 Min. : 1 A:90 Min. :175 Min. :18.00
## 1st Qu.: 68.25 1st Qu.: 4 B:90 1st Qu.:185 1st Qu.:26.00
## Median :135.50 Median : 8 C:90 Median :200 Median :31.00
## Mean :135.50 Mean : 8 Mean :200 Mean :32.12
## 3rd Qu.:202.75 3rd Qu.:12 3rd Qu.:215 3rd Qu.:36.75
## Max. :270.00 Max. :15 Max. :225 Max. :63.00
## temp
## Min. :175
## 1st Qu.:185
## Median :200
## Mean :200
## 3rd Qu.:215
## Max. :225Для количественных переменных эта функция выдает минимальное и максимальное значение, среднее арифметическое, медиану, нижний (1st Qu.) и верхний (3rd Qu.) квартиль. Нижний квартиль – значение, которое 25% значений в выборке не превышают, а верхний квартиль – значение, которое 75% значений в выборке не превышают. Для качественных переменных (текстовые, факторные), R будет выводить количество значений по каждой группе (уровню).
В данном случае по выдаче R мы можем определить следующее. Всего в базе данных у нас 270 наблюдений (переменная X здесь служит id наблюдений, а ее максимальное значение 270), значит, в рамках исследования было приготовлено 270 шоколадных тортов. Минимальная температура, при которой выпекали торты, равна 175 градусам, максимальная - 225. Средняя температура, при которой выпекали торты, равна 200. Медианное значение температуры в данном случае совпадает со средним значением – в половине случаев температура при выпечке не превышала 200 градусов. Нижний квартиль равен 185 градусам – в 25% случаев торты выпекались при температуре не выше 185 градусов, верхний квартиль равен 215 – 75% случаев температура не превышала 215 градусов (или в 25% случаев превышала!).
Необязательно выводить описательные статистики для всех переменных в базе данных, можно вывести описание одной переменной:
summary(df$temperature)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 175 185 200 200 215 225Или проделать то же самое для нескольких переменных:
summary(df[5:6])## angle temp
## Min. :18.00 Min. :175
## 1st Qu.:26.00 1st Qu.:185
## Median :31.00 Median :200
## Mean :32.12 Mean :200
## 3rd Qu.:36.75 3rd Qu.:215
## Max. :63.00 Max. :225Или так:
summary(df[ , c("angle", "temperature")])## angle temperature
## Min. :18.00 Min. :175
## 1st Qu.:26.00 1st Qu.:185
## Median :31.00 Median :200
## Mean :32.12 Mean :200
## 3rd Qu.:36.75 3rd Qu.:215
## Max. :63.00 Max. :225Точно так же необязательно выводить все статистики сразу. Можно запрашивать по отдельности:
min(df$angle) # минимум## [1] 18max(df$angle) # максимум## [1] 63mean(df$angle) # среднее## [1] 32.12222median(df$angle) # медиана## [1] 31В данном случае все показатели считаются без проблем, потому что в базе данных все строки полностью заполнены. В прошлый раз было хуже – вместо среднего значения R выдавал NA, так как в переменной были пропущенные значения. Напомню, чтобы решить эту проблему, мы прописывали дополнительный аргумент na.rm = TRUE, который говорил R не учитывать пропущенные значения при расчете статистик (из самой базы значения при этом не выкидываются!).
mean(df$angle, na.rm = TRUE) ## [1] 32.12222Это не единственные описательные статистики, которые можно вывести. Часто нас интересует не только среднее (или медианное) значение, а разброс значений относительно этого среднего. Для этого можем посчитать выборочную дисперсию или стандартное отклонение.
var(df$temperature) # дисперсия## [1] 292.7509sd(df$temperature) # стандартное отклонение## [1] 17.10997Но сами по себе эти значения не очень информативны – по ним сложно понять, насколько однородны наши данные (сильно ли они разбросаны относительно среднего значения). Для того, чтобы оценить степень однородности наших данных, нашей выборки, можно воспользоваться таким показателем как коэффициент вариации. Коэффициент вариации считается несложно: стандартное отклонение нужно поделить на среднее значение. Обычно значение коэффициента вариации, взятое по модулю, лежит в пределах от 0 до 1, но иногда, если данные очень разнородны (стандартное отклонение большое), оно может быть больше 1.
Вопрос: а когда коэффициент вариации нужно использовать с осторожностью?
Часто коэффициент вариации выражают в процентах. Давайте напишем код, который будет считать коэффициент вариации в процентах для переменной angle. Ваш ход код, маэстро!
# код
sd(df$angle) / mean(df$angle) * 100## [1] 25.56657Конечно, перечисленными примерами не ограничивается описание количественных данных. Есть еще процентили, децили и прочее… Чтобы не запутаться, обсуждать их не будем, почитать про их вычисление в R и интерпретацию можно, например, здесь.
Пока мы обсудили только описательные статистики для количественных переменных. А как быть с качественными? Какую информацию по ним можно получить? Число наблюдений, соответствующих каждому значению (классу):
table(df$recipe)##
## A B C
## 90 90 90В прошлый раз мы говорили о том, что для показателей, измеренных в качественной шкале, вычислять среднее или медиану бессмысленно, нужно смотреть на моду. Искать специальную функцию не нужно, достаточно помнить, что мода – это значение, которое встречается в выборке чаще всего (да, мода может быть не одна).
Кстати, для количественных переменных функция table() тоже хорошо работает:
table(df$temperature)##
## 175 185 195 205 215 225
## 45 45 45 45 45 45Прежде чем перейти к более красивым вещам, предлагаю выполнить маленькое задание.
Задание
temperature. Однородны ли значения температуры, при которой выпекали торты?recipe.temperature. Проинтерпретируйте.Решение
mean(df$angle) #1## [1] 32.12222sd(df$temperature) / mean(df$temperature) * 100 # 2## [1] 8.554983table(df$recipe) # 3##
## A B C
## 90 90 90quantile(df$temperature, prob = c(0.1, 0.2, 0.3)) # 4## 10% 20% 30%
## 175 185 185Обратимся к библиотеке:
library(dplyr)Выведем стандартную таблицу, которую часто можно увидеть в исследованиях в разделе, посвященном описанию данных. В этой таблице будет указано число наблюдений, среднее значение и стандартное отклонение.
df %>% summarise(observations = n(),
avg = mean(temperature),
sd = sd(temperature)) ## observations avg sd
## 1 270 200 17.10997Добавим в таблицу коэффициент вариации:
df %>% summarise(observations = n(),
avg = mean(temperature),
sd = sd(temperature)) %>% mutate(cv = sd / avg) * 100## observations avg sd cv
## 1 27000 20000 1710.997 8.554983Не очень важная, но полезная деталь: в качестве названий можно использовать и “сложные названия” – включающие скобки и прочие символы. Но тогда название столбца нужно вводить в кавычках:
df %>% summarise(observations = n(),
avg = mean(temperature),
sd = sd(temperature)) %>% mutate("cv (in %)" = sd / avg) * 100## observations avg sd cv (in %)
## 1 27000 20000 1710.997 8.554983Можем вывести описательные статистики для групп наблюдений. Сгруппируем наблюдения по рецепту, по которому готовился торт.
df %>% group_by(recipe) %>%
summarise(observations = n(),
avg = mean(temperature),
sd = sd(temperature))## # A tibble: 3 x 4
## recipe observations avg sd
## <fctr> <int> <dbl> <dbl>
## 1 A 90 200 17.17393
## 2 B 90 200 17.17393
## 3 C 90 200 17.17393R позволяет выгружать не только таблицы с результатами регрессии, но и почти любые таблицы, например, таблицу с описательными статистиками. Единственная проблема: просто так в Word выгрузить ничего не получится, R дружелюбен только по отношению к LaTeX :) Но эту проблему можно решить.
Для начала загрузим библиотеку stargazer. Она используется для выгрузки таблиц “во внешний мир”.
install.packages("stargazer")library(stargazer)##
## Please cite as:## Hlavac, Marek (2015). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2. http://CRAN.R-project.org/package=stargazerДля того, чтобы вывести таблицу с описательными статистиками, нужно просто набрать stargazer() и указать базу данных:
stargazer(df) ##
## % Table created by stargazer v.5.2 by Marek Hlavac, Harvard University. E-mail: hlavac at fas.harvard.edu
## % Date and time: Сб, окт 14, 2017 - 04:37:22
## \begin{table}[!htbp] \centering
## \caption{}
## \label{}
## \begin{tabular}{@{\extracolsep{5pt}}lccccc}
## \\[-1.8ex]\hline
## \hline \\[-1.8ex]
## Statistic & \multicolumn{1}{c}{N} & \multicolumn{1}{c}{Mean} & \multicolumn{1}{c}{St. Dev.} & \multicolumn{1}{c}{Min} & \multicolumn{1}{c}{Max} \\
## \hline \\[-1.8ex]
## X & 270 & 135.500 & 78.086 & 1 & 270 \\
## replicate & 270 & 8.000 & 4.329 & 1 & 15 \\
## temperature & 270 & 200.000 & 17.110 & 175 & 225 \\
## angle & 270 & 32.122 & 8.213 & 18 & 63 \\
## temp & 270 & 200.000 & 17.110 & 175 & 225 \\
## \hline \\[-1.8ex]
## \end{tabular}
## \end{table}Для тех, кто умеет работать в LaTeX, все просто – R выдал теховский код, который можно просто скопировать в tex-файл. Обратите внимание: в начале кода в качестве комментария всегда указывается, есть ли необходимость догружать (указывать в преамбуле) специальные пакеты или нет. В данном случае не нужно.
Тем, кто в LaTeX не работает, легче не стало. Но выход есть. Целых два.
Выход 1. Stargazer умеет не только выводить код в консоль, но и сохранять результат в отдельный файл. Этот файл может быть теховским файлом (.tex), текстовым (.txt) и html-файлом (.htmlили .htm). Word умеет открывать файлы с расширением .htm. И Libre Office тоже.
##
## <table style="text-align:center"><tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Statistic</td><td>N</td><td>Mean</td><td>St. Dev.</td><td>Min</td><td>Max</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">X</td><td>270</td><td>135.500</td><td>78.086</td><td>1</td><td>270</td></tr>
## <tr><td style="text-align:left">replicate</td><td>270</td><td>8.000</td><td>4.329</td><td>1</td><td>15</td></tr>
## <tr><td style="text-align:left">temperature</td><td>270</td><td>200.000</td><td>17.110</td><td>175</td><td>225</td></tr>
## <tr><td style="text-align:left">angle</td><td>270</td><td>32.122</td><td>8.213</td><td>18</td><td>63</td></tr>
## <tr><td style="text-align:left">temp</td><td>270</td><td>200.000</td><td>17.110</td><td>175</td><td>225</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr></table>Даже если редактор документов не хочет открывать файл или делает это странно, не нужно расстраиваться. Html-файлы, как и страницы в Интернете, можно просто окрыть в браузере. И нужную таблицу отскринить.
Выход 2. Воспользоваться онлайн-редактором для LaTeX, скопировать в него код из консоли и скомпилировать симпатичный pdf-файл. Можно зарегистрироваться на сайте ShareLatex, это бесплатно. Затем зайти и нажать Создать проект – Новый проект. Когда мы дадим имя проекту, откроется поле для работы. В левой части – код, в правой части – pdf-файл. Нужно добавить в преамбулу (часть перед \begin{document}) строку \usepackage[english, russian]{babel}, чтобы в документе нормально отображалась кириллица. Затем мы можем убрать строку для заголовка (\maketitle) и раздела (\section) и вставить туда код из консоли:
stargazer(df)##
## % Table created by stargazer v.5.2 by Marek Hlavac, Harvard University. E-mail: hlavac at fas.harvard.edu
## % Date and time: Сб, окт 14, 2017 - 04:37:22
## \begin{table}[!htbp] \centering
## \caption{}
## \label{}
## \begin{tabular}{@{\extracolsep{5pt}}lccccc}
## \\[-1.8ex]\hline
## \hline \\[-1.8ex]
## Statistic & \multicolumn{1}{c}{N} & \multicolumn{1}{c}{Mean} & \multicolumn{1}{c}{St. Dev.} & \multicolumn{1}{c}{Min} & \multicolumn{1}{c}{Max} \\
## \hline \\[-1.8ex]
## X & 270 & 135.500 & 78.086 & 1 & 270 \\
## replicate & 270 & 8.000 & 4.329 & 1 & 15 \\
## temperature & 270 & 200.000 & 17.110 & 175 & 225 \\
## angle & 270 & 32.122 & 8.213 & 18 & 63 \\
## temp & 270 & 200.000 & 17.110 & 175 & 225 \\
## \hline \\[-1.8ex]
## \end{tabular}
## \end{table}Нажать Компилировать.
Теперь для простоты вернемся к html-файлам и будем приводить таблицу в порядок. Для начала дадим ей заголовок.
stargazer(df, type = "html",
title = "Summary statistics",
out = "cake_summary.htm") ##
## <table style="text-align:center"><caption><strong>Summary statistics</strong></caption>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Statistic</td><td>N</td><td>Mean</td><td>St. Dev.</td><td>Min</td><td>Max</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">X</td><td>270</td><td>135.500</td><td>78.086</td><td>1</td><td>270</td></tr>
## <tr><td style="text-align:left">replicate</td><td>270</td><td>8.000</td><td>4.329</td><td>1</td><td>15</td></tr>
## <tr><td style="text-align:left">temperature</td><td>270</td><td>200.000</td><td>17.110</td><td>175</td><td>225</td></tr>
## <tr><td style="text-align:left">angle</td><td>270</td><td>32.122</td><td>8.213</td><td>18</td><td>63</td></tr>
## <tr><td style="text-align:left">temp</td><td>270</td><td>200.000</td><td>17.110</td><td>175</td><td>225</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr></table>Теперь уберем лишние переменные:
stargazer(subset(df, select = c("angle", "temperature", "replicate")),
type = "html",
title = "Summary statistics",
out = "cake_summary.htm") ##
## <table style="text-align:center"><caption><strong>Summary statistics</strong></caption>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Statistic</td><td>N</td><td>Mean</td><td>St. Dev.</td><td>Min</td><td>Max</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">angle</td><td>270</td><td>32.122</td><td>8.213</td><td>18</td><td>63</td></tr>
## <tr><td style="text-align:left">temperature</td><td>270</td><td>200.000</td><td>17.110</td><td>175</td><td>225</td></tr>
## <tr><td style="text-align:left">replicate</td><td>270</td><td>8.000</td><td>4.329</td><td>1</td><td>15</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr></table>Хотим округлить все значения до второго знака после запятой:
stargazer(subset(df, select = c("angle", "temperature", "replicate")),
type = "html",
out = "cake_summary.htm",
title = "Summary statistics",
digits = 2)##
## <table style="text-align:center"><caption><strong>Summary statistics</strong></caption>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Statistic</td><td>N</td><td>Mean</td><td>St. Dev.</td><td>Min</td><td>Max</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">angle</td><td>270</td><td>32.12</td><td>8.21</td><td>18</td><td>63</td></tr>
## <tr><td style="text-align:left">temperature</td><td>270</td><td>200.00</td><td>17.11</td><td>175</td><td>225</td></tr>
## <tr><td style="text-align:left">replicate</td><td>270</td><td>8.00</td><td>4.33</td><td>1</td><td>15</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr></table>Хотим добавить комментарий (текст после таблицы):
stargazer(subset(df, select = c("angle", "temperature", "replicate")),
type = "html",
out = "cake_summary.htm",
title = "Summary statistics",
digits = 2, notes = "Notes:")##
## <table style="text-align:center"><caption><strong>Summary statistics</strong></caption>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Statistic</td><td>N</td><td>Mean</td><td>St. Dev.</td><td>Min</td><td>Max</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">angle</td><td>270</td><td>32.12</td><td>8.21</td><td>18</td><td>63</td></tr>
## <tr><td style="text-align:left">temperature</td><td>270</td><td>200.00</td><td>17.11</td><td>175</td><td>225</td></tr>
## <tr><td style="text-align:left">replicate</td><td>270</td><td>8.00</td><td>4.33</td><td>1</td><td>15</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td colspan="6" style="text-align:left">Notes:</td></tr>
## </table>