---
title: "Основы программирования в R"
subtitle: "Загрузка данных и их описание"
author: "Алла Тамбовцева, НИУ ВШЭ"
output:
  pdf_document:
    latex_engine: xelatex
    toc: true
mainfont: CMU Serif
header-includes:
- \usepackage[russian]{babel}
- \usepackage{hyperref}
- \hypersetup{colorlinks = true, urlcolor = blue, linkcolor=magenta}
---

## Загрузка CSV-файла и кодировка

Загрузим данные из файла `firtree.csv`, в котором хранятся результаты вымышленного опроса посетителей ёлочного базара, и сохраним их в переменную `tree`.  

**Показатели в файле:**

* `gender` – пол респондента;

* `ftype` – тип хвойного дерева, которое оценивал респондент;

* `height` – высота хвойного дерева в сантиметрах;

* `expenses` – сумма (в рублях), которую респондент готов отдать за хвойное дерево;

* `score` – балл, на который респондент оценил вид хвойного дерева (1 – очень плохо, 5 – отлично);

* `wish` – ответ на вопрос «Хотели бы, чтобы Вам подарили такое хвойное дерево?» (да, нет).

Так как в файле есть текст на кириллице и создавался он на Mac OS/Linux, необходимо добавить аргумент `encoding` со значением кодировки `"UTF-8"`, иначе на Windows вместо русских букв в тексте будут крокозябры:

```{r, eval=FALSE}
tree <- read.csv("D:/Downloads/firtree.csv", 
                 encoding = "UTF-8")
```

```{r, inclide=FALSE}
tree <- read.csv("/Users/allat/Desktop/firtree.csv")
setwd("/Users/allat/Desktop/")
```

Похожая проблема может возникнуть, если наоборот, файл с кириллицей создавался на Windows, а загружается в R на Mac/Linux. Тогда нужно будет выставить кодировку `"Windows-1251"`. 

Если добавление кодировки в `encoding` не решает проблему (текст не отображается в читаемом виде), нужно перед работой с файлом запустить следующую строчку кода:

```{r, eval=FALSE}
Sys.setlocale("LC_CTYPE", "ru_RU.UTF-8")
```

Этот код сохранит настройки языка и кодировки, и файлы на русском языке будут благополучно открываться.

Какие ещё проблемы могут возникнуть при работе с CSV-файлом? Во-первых, в качестве разделителя столбцов вместо привычной запятой может использоваться точка с запятой. Если этот факт мы не учтём, R выдаст ошибку и файл с данными не загрузит. Справиться с этой проблемой поможет опция `sep` (от *separator*). Загрузим файл `test2.csv`, используя ссылку на странице курса:

```{r}
test2 <- read.csv("https://raw.githubusercontent.com/allatambov/PyDat-0919/master/lectures-seminars/7-pandas/test2.csv", sep = ";")
```

Во-вторых, в качестве десятичного разделителя в файле может использоваться запятая, вместо принятой в R точки. В таком случае дробные числа будут считываться R как некоторый текст. Если присмотреться, именно это мы увидим в датафрейме `test2`:

```{r}
test2
```
Столбцы `B` и `C` имеют тип *character*, хотя задумывались они явно как числовые. Исправим это, добавим ещё один аргумент — `dec` (десятичный разделитель):

```{r}
test2 <- read.csv("https://raw.githubusercontent.com/allatambov/PyDat-0919/master/lectures-seminars/7-pandas/test2.csv", sep = ";", dec = ",")
test2
```
Теперь всё в порядке. 

## Загрузка файла Excel

Для загрузки файла Excel (`xls` или `xlsx`), нам понадобится библиотека `readxl`. Установим её:

```{r, eval=FALSE}
install.packages("readxl")
```

Обратимся к библиотеке, чтобы R видел, откуда брать функции для загрузки файлов Excel (аналог импортирования библиотеки через конструкцию `import` в Python):

```{r, message=FALSE, warning=FALSE}
library(readxl)
```

Теперь вызовем функция `read.excel()` и загрузим тестовый файл `test1.xlsx`:

```{r}
test1 <- read_excel("/Users/allat/Desktop/test1.xlsx")
test1
```

## Выбор рабочей папки

Рабочая папка — папка, из которой R запускается по умолчанию. На практике это означает, что R по умолчанию видит только те файлы, которые находятся в рабочей папке, поэтому при обращении к ним не нужно прописывать полный путь, достаточно указать только название. Узнать, какая папка является рабочей, можно с помощью функции `getwd()`:

```{r}
# wd - working directory
getwd()
```

Чтобы изменить рабочую папку, нужно вызвать функцию `setwd()` и указать путь к новой папке:

```{r}
setwd("/Users/allat/Desktop/")
```

Теперь файл `firtree.csv`, лежащий на рабочем столе (`Desktop`), можно загрузить довольно просто:

```{r}
tree <- read.csv("firtree.csv")
```

Иногда, если очень не хочется возиться с папками и путями, можно обратиться к функции `file.choose()`, она откроет обычное окно для выбора файла:

```{r, eval=FALSE}
tree <- read.csv(file.choose())
```

## Просмотр данных

Вернёмся к файлу с данными опроса на ёлочном базаре и датафрейму `tree`. Посмотрим на датафрейм — функция `View()` открывает датафрейм в отдельной вкладке:

```{r}
View(tree)
```

**Внимание:** первая буква у `View()` заглавная!

Теперь запросим первые строки датафрейма:

```{r}
head(tree)
```

По умолчанию функция `head()` выдает первые 6 строк, но это можно изменить:

```{r}
# первые 8 строк
head(tree, 8)
```

Аналогичным образом можно вывести последние строки датафрейма:

```{r}
View(tail(tree))
```

Здесь функцию `tail()` мы заключили в `View()`, чтобы строки выводились в удобном формате — не в консоль, а в отдельном окне.

## Техническое описание данных

Для начала запросим размерность датафрейма: число строк и число столбцов.

```{r}
# 1200 строк и 7 столбцов
dim(tree)
```

Функция `dim()` возвращает вектор из двух элементов, причем на первом месте всегда идёт число строк, на втором — число столбцов. Если нам нужно только число строк или только число столбцов, можно выбрать нужный элемент по индексу, а можно поступить проще — воспользоваться готовыми функциями.

Функция `ncol()` возвращает число столбцов, а функция `nrow()` — число строк. 

```{r}
ncol(tree)
```

```{r}
nrow(tree)
```

Если мы хотим получить техническое описание датафрейма — сколько в нём строк и столбцов, какого типа эти столбцы, можно воспользоваться функцией `str()`. Эта функция (`str` от  *structure*) возвращает структуру любого объекта, не только датафрейма, поэтому, если не совсем ясно, какой объект выдала какая-нибудь функция из неизвестной библиотеки, можно смело её использовать. 

Посмотрим на структуру датафрейма `tree`:

```{r}
str(tree)
```

Со столбцами `X`, `height`, `score` и `expenses` всё понятно, это обычные целочисленные столбцы типа *integer*. С остальными столбцами интереснее — они имеют тип `factor`. `Levels` здесь — это уникальные значения в векторе.

Тип `factor` используется в тех случаях, когда нечисловые, качественные, значения кодируются числами. Другими словами, когда числа «ненастоящие», когда с ними нельзя работать как с числами в математике. Например, если вместо значений `"female"` и `"male"` в столбце `gender` мы будем ставить 0 и 1, мы всё равно не сможем говорить, что 1 здесь больше 0, это какие-то наши условные обозначения, результат договоренности. Или, например, если мы будем кодировать любимый цвет респондента числами от 1 до 4 (красный, жёлтый, зелёный, синий), мы не сможем сравнивать эти числа и утверждать, что 4 в два раза больше 2, потому что это то же самое, что сравнивать слова «жёлтый» и «синий». Считать среднее значение по такому набору чисел тоже неправильно, даже если технически мы можем все числа сложить и поделить на их количество, потому что результат будет неинтерпретируемым. Ведь непонятно, что такое средний цвет, равный, к примеру, 2.5.

Особого внимания заслуживает столбец `wish`. Помимо очевидных значений `"да"` и `"нет"` здесь есть значение `""`. На самом деле это пустые ячейки, которые считались в R таким образом. Чтобы они нам не мешали, давайте ещё раз загрузим файл, добавив опцию `na.strings = ""`, которая принудит R считать такие ячейки за полноценные пропущенные значения `NA`:

```{r, eval=FALSE}
tree <- read.csv("D:/Downloads/firtree.csv", 
                 encoding = "UTF-8", na.strings = "")
```

```{r, include=FALSE}
tree <- read.csv("/Users/allat/Desktop/firtree.csv", 
                 encoding = "UTF-8", na.strings = "")
```

Посмотрим на структуру обновленного датафрейма:

```{r}
str(tree)
```

Всё исправилось!

## Поиск пропущенных значений

Теперь мы точно знаем, что в некоторых столбцах есть пропущенные значения (NA's). Попробуем их посчитать. Для начала воспользуемся функцией `complete.cases()`, которая вернёт нам вектор из значений `TRUE` и `FALSE`, где `TRUE` означает, что строка в таблице не содержит пропущенные значения (*case* — это строка, то есть одно наблюдение). Выведем первые несколько значений вектора:

```{r}
head(complete.cases(tree))
```

Теперь, чтобы посчитать число полностью заполненных строк, нам достаточно посчитать число `TRUE`. Сделать это очень просто: R воспринимает значения `TRUE` как 1, а `FALSE` — как 0, поэтому можно просто суммировать все значения в векторе выше:

```{r}
sum(complete.cases(tree))
```

Но нам нужен противоположный набор значений, ведь мы хотим посчитать число строк с пропущенными значениями! Поэтому к `complete.cases()` нужно добавить отрицание. Отрицание в программировании обычно задаётся с помощью восклицательного знака. Поставим его перед функцией и получим «перевёрнутый» вектор, где `TRUE` и `FALSE` поменялись местами.

```{r}
sum(!complete.cases(tree))
```

Получается, в датафрейме `tree` у нас есть две строки, в которых есть хотя бы одно пропущенное значение.

**Важно:** в R есть ещё одна функция для поиска пропущенных значений — `is.na()`: 

```{r}
sum(is.na(tree))
```

В нашем случае результаты с `complete.cases()` и `is.na()` совпадают, но так будет не всегда. Функция `complete.cases()` проверяет заполненность строк, а функция `is.na()` — заполненность ячеек. Допустим, у нас есть маленький датафрейм такого вида:

```{r}
test <- cbind.data.frame(a = c(NA, 2, 3), 
                         b = c(NA, NA, 1))
test
```

В нём две строки, содержащие хотя бы один `NA`, но всего пропущенных значений три. Сравним результаты:

```{r}
sum(!complete.cases(test))
```

```{r}
sum(is.na(test))
```

## Содержательное описание данных

Выведем описательные статистики по всему датафрейму `tree` с помощью функции `summary()`:

```{r}
summary(tree)
```

Для количественных показателей функция возвращает минимальное и максимальное значения (`Min.` и `Max`), среднее арифметическое и медиану (`Mean` и `Median`), а также нижний и верхний квантили (`1st Qu.` и `3rd Qu.`). Так, для столбца `height` получаем:

* высота 50% деревьев в данных не превышает значение 157 см;
* высота 25% деревьев в данных не превышает значение 115 см;
* высота 75% деревьев в данных не превышает значение 203.2 см.

Для текстовых показателей функция не возвращает ничего интересного. Чтобы это исправить, нужно считать текстовые столбцы как факторные (`factor`). Факторный тип — особый тип данных в R, к нему можно относиться как к типу, который хранит неколичественные значения, но при этом присваивает им числовые метки. Так, факторный вектор может хранить ответы «да» и «нет», но при этом R будет знать, что значению «да» соответствует число 1, а  значению «нет» — 2.

Для того, чтобы текстовые столбцы считались как факторные, при загрузке файла нужно добавить аргумент `stringsAsfactors=TRUE`:

```{r}
tree <- read.csv("firtree.csv", na.strings = "", stringsAsFactors = TRUE)
summary(tree)
```

Теперь для нечисловых (факторных) столбцов  функция `summary()` показывает частоты — сколько раз то или иное значение встречается в столбце. Количество пропущенных значений тоже учитывается. 

## Паттерны пропущенных значений

Для дальнейшей работы с пропущенными значениями нам понадобится дополнительная библиотека `VIM`. Установим её.

```{r, eval=FALSE}
install.packages("VIM")
```

Обратимся к ней:

```{r, message=FALSE, warning=FALSE}
library(VIM)
```

Выведем графики, которые покажут, в каких переменных пропущенных значений больше всего и как выглядит таблица с пропущенными значениями (паттерны пропущенных значений).

На графике слева показано, с какой частотой встречаются пропущенные значения в той или иной переменной. На графике справа показано, в каких комбинациях эти пропущенные значения встречаются.

```{r}
aggr(tree)
```

Следующий график отвечает за заполненность наблюдений (красным цветом отмечены пропущенные значения, остальное — заполненные значения, чем темнее цвет, тем больше значение). По вертикальной оси — номер строки в датафрейме, id наблюдения.

```{r}
matrixplot(tree)
```

Так как в датафрейме `tree` всего две строчки с пропущенными значениями, и они не рядом, на графике их почти не видно. Но если пропусков много, этот график их покажет, сразу станет видно красные «дыры» на фоне серых и черных полосок. Для примера можем посмотреть на тот же график для `test`:

```{r}
matrixplot(test)
```

Датафрейм маленький, и по графику сразу видно, что ячеек с пропущенными значениями много, если сравнивать с общим числом ячеек в датафрейме.

## Описание качественных данных

Если нас интересует отдельный столбец датафрейма, его можно выбрать через `$`:

```{r}
head(tree$wish) # первые несколько значений
```

Выбрать, а дальше описывать отдельно. Если показатель качественный (текстовый или факторный), для него логично определить уникальные значения:

```{r}
unique(tree$wish)
```

И соответствующие им частоты:

```{r}
table(tree$wish)
```

Потом эту таблицу частот можно поместить внутрь функции `barplot()` и построить столбиковую диаграмму:

```{r, eval = FALSE}
barplot(table(tree$wish))
```

```{r, include = FALSE}
dev.copy(png, "bar.png")
barplot(table(tree$wish))
dev.off()
```
\begin{figure}
\centering
\includegraphics[height=10cm]{bar.png}
\end{figure}

Можем добавить цвета:

```{r, eval = FALSE}
barplot(table(tree$wish), 
        col = c("lavender", "darkviolet"))
```

```{r, include=FALSE}
dev.copy(png, "bar2.png")
barplot(table(tree$wish), 
        col = c("lavender", "darkviolet"))
dev.off()
```

\begin{figure}
\centering
\includegraphics[height=10cm]{bar2.png}
\end{figure}

График далёк от идеального: подписей нет, вертикальная ось коротковата... Но настройкой графиков мы будем заниматься позже, пока просто смотрим, что возможность быстро построить график есть.

## Описание количественных данных

Уже знакомую нам функцию `summary()` мы можем применить и к отдельному столбцу (и к вектору вне датафрейма тоже):

```{r}
summary(tree$expenses)
```

Здесь уже всё знакомо.

Теперь посмотрим на более подробную выдачу R с описательными статистиками. Чтобы это сделать, нам понадобится библиотека `psych`, которая содержит набор функций, часто используемых в психометрических исследованиях. Установим её:

```{r, eval=FALSE}
install.packages("psych")
```

Обратимся к библиотеке через `library()`:

```{r, message=FALSE, warning=FALSE}
library(psych)
```

Теперь запросим описательные статистики для столбца `expenses` с помощью функции `describe()`:

```{r}
describe(tree$expenses)
```

Что есть что?

* `vars`: число описываемых переменных (не путать с `var` для дисперсии);
* `n`: число наблюдений;
* `mean`: среднее арифметическое, выборочное среднее;
* `sd`: стандартное отклонение;
* `median`: медиана;
* `trimmed`: усечённое среднее, среднее по цензурированной выборке (см. ниже);
* `mad`: медианное значение абсолютного отклонения от медианы (нам не понадобится);
* `min`, `max`: минимальное и максимальное значение;
* `range`: размах;
* `skew`: коэффициент асимметрии или скошенности (см.нижк);
* `kurtosis`: коэффициент эксцесса (см. ниже);
* `se`: стандартная ошибка среднего;

Подробнее про некоторые статистики.

*Усечённое среднее, среднее по цензурированной выборке*

Считается так: выборка упорядочивается по возрастанию, из неё убирается 5% наблюдений слева и справа (наименьшие и наибольшие), потом по такой усечённой или цензурированной выборке считается обычное среднее арифметическое.

Наравне с медианой считается более устойчивой оценкой среднего, так как после усечения выборки такой показатель уже несильно зависит от слишком больших или слишном маленьких (нетипичных) значений в выборке. То есть, при наличии нетипичных наблюдений в выборке (выбросов) такое среднее более адекватно отражает реальность, чем обычное среднее арифметическое.

*Коэффициент асимметрии*

Показатель принимает значения примерно от $-3$ до $3$. Значение $0$ соответствует симметричному распределению (например, нормальному, вспомните график плотности, симметричный относительно математического ожидания). Значения меньше $0$ соответствуют распределению, которое скошено влево (длинный хвост «слева»), значения больше $0$ соответствуют распределению, которое скошено вправо (длинный «хвост» справа).

В нашем случае распределение почти симметричное, коэффициент близок к нулю, но при это оно немного скошено вправо, поэтому значение больше $0$.

*Коэффициент эксцесса*

Показатель принимает значения примерно от $-3$ до $3$ и отвечает за выраженность пика распределения. Чем больше значение коэффициента, тем более выраженный пик. Стандартное нормальное распределение имеет коэффициент эксцесса равный $0$. Отрицательные значения коэффициента соответствуют более «плоским» и «гладким» распределениям, у которых пик не такой заметный. Посмотрите на картинку здесь и сравните.

В нашем случае распределение несильно отличается от нормального, поэтому коэффициент близок к нулю.

Библиотека `psych` удобна тем, что она содержит функцию `describeBy()`, которая позволяет выводить описательные статистики по группам. Нет необходимости отфильтровывать нужные строки и сохранять их в отдельные датасеты, можно просто указать группирующую переменную. Например, сравним, сколько на хвойные деревья могут тратить мужчины и женщины:

```{r}
describeBy(tree$expenses, tree$gender)
```

Очень удобно!

Если нас интересует только определённая характеристика столбца, можем воспользоваться базовыми, уже знакомыми нам, функциями.

```{r}
min(tree$expenses)
```

```{r}
max(tree$expenses)
```

```{r}
mean(tree$expenses)
```

```{r}
median(tree$expenses)
```

```{r}
var(tree$expenses)
```

```{r}
sd(tree$expenses)
```

Однако у всех этих функций есть одна особенность — они возвращают `NA`, если в столбце или векторе есть хотя бы одно пропущенное значение. Попробуем посчитать среднее для вектора с `NA`:

```{r}
mean(c(7, 5, NA, 9))
```

Нет ответа, плюс, получили предупреждение о наличие `NA`. Чтобы этого избежать, можно добавить опцию `na.rm = TRUE`, сокращение от  *NA remove*:

```{r}
mean(c(7, 5, NA, 9), na.rm = TRUE)
```

Пропущенные значения не удаляются из самого вектора, но не учитываются при вычислении среднего. То же будет актуально и для других характеристик (минимум, медиана и прочие).

Напоследок построим гистограмму:

```{r}
hist(tree$expenses)
```

Добавим цвет: 

```{r}
hist(tree$expenses, col = "limegreen")
```

Выбор цветов в R богатый, список всех цветов с примерами можно посмотреть [здесь](http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf). При желании можно вводить не название цвета, а его код в формате RGB или HEX. Пример с цветом в формате HEX (*hexadecimal*):

```{r}
hist(tree$expenses, col = "#266136")
```

Про форматы цветов можно посмотреть [здесь](https://www.w3schools.com/Colors/default.asp). 

Настройку графиков и наведение красоты мы обсудим позже.
