В прошлый раз (в ноябре 2014-го; мне очень стыдно, что я так затянул с продолжением!) я рассказывал о базовых возможностях языка R. Несмотря на наличие всех привычных управляющих конструкций, таких как циклы и условные блоки, классический подход к обработке данных на основе итерации далеко не лучшее решение, поскольку циклы в R необыкновенно медлительны. Поэтому сейчас я расскажу, как на самом деле нужно работать с данными, чтобы процесс вычислений не заставлял тебя выпивать слишком много чашек кофе в ожидании результата. Кроме того, некоторое время я посвящу рассказу о том, как пользоваться современными средствами визуализации данных в R. Потому что удобство представления результатов обработки данных на практике не менее важно, чем сами результаты. Начнем с простого.

Векторные операции

Как мы помним, базовым типом в R является вовсе не число, а вектор, и основные арифметические операции действуют на векторы поэлементно:

> x <- 1:6; y <- 11:17
> x + y
[1] 12 14 16 18 20 22 18
> x > 2
[1] FALSE FALSE  TRUE  TRUE  TRUE  TRUE
> x * y
[1] 11 24 39 56 75 96 17
> x / y
[1] 0.09090909 0.16666667 0.23076923 0.28571429 0.33333333 0.37500000
[7] 0.05882353

Тут все довольно просто, однако вполне логично задаться вопросом: что же будет, если длина векторов не совпадет? Если мы, скажем, напишем k <- 2, то будет ли x * k соответствовать умножению вектора на число в математическом смысле? Короткий ответ — да. В более общем случае, когда длина векторов не совпадает, меньший вектор просто продолжается повторением:

> z <- c(1, 0.5)
> x * z
[1] 1 1 3 2 5 3

Примерно так же обстоят дела и с матрицами.

> x <- matrix(1:4, 2, 2); y <- matrix(rep(2,4), 2, 2)
> x * y
     [,1] [,2]
[1,]    2    6
[2,]    4    8
> x / y
     [,1] [,2]
[1,]  0.5  1.5
[2,]  1.0  2.0

При этом «нормальное», а не поразрядное умножение матриц будет выглядеть так:

> x %*% y
     [,1] [,2]
[1,]    8    8
[2,]   12   12

Все это, конечно, очень хорошо, однако что же делать, когда нам нужно применять свои собственные функции к элементам векторов или матриц, то есть как это можно сделать без цикла? Подход, который используется в R для решения данной проблемы, очень схож с тем, к чему мы привыкли в функциональных языках, — все это напоминает функцию map в Python или Haskell.

Полезная функция lapply и ее друзья

Первая функция в этом семействе — это lapply. Она позволяет применять заданную функцию к каждому элементу списка или вектора. Причем результатом будет именно список независимо от типа аргумента. Простейший пример с применением лямбда-функций:

> q <- lapply(c(1,2,4), function(x) x^2)
> q
[[1]]
[1] 1

[[2]]
[1] 4

[[3]]
[1] 16

Если функция, которую нужно применить к списку или вектору, требует более одного аргумента, то эти аргументы можно передать через lapply.

> q <- lapply(c(1,2,4), function(x, y) x^2 + y, 3)

Со списком функция работает аналогичным образом:

> x <- list(a=rnorm(10), b=1:10)
> lapply(x, mean)

Здесь функция rnorm задает нормальное распределение (в данном случае десять нормально распределенных чисел в диапазоне от 0 до 1), а mean вычисляет среднее значение. Функция sapply полностью аналогична функции lapply за исключением того, что она пытается упростить результат. К примеру, если каждый элемент списка длины 1, то вместо списка вернется вектор:

> sapply(c(1,2,4), function(x) x^2)
[1]  1  4 16

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

> x <- list(1:4, 5:8)
> sapply(x, function(x) x^2)
     [,1] [,2]
[1,]    1   25
[2,]    4   36
[3,]    9   49
[4,]   16   64

Для работы с матрицами удобно использовать функцию apply:

> x <- matrix(rnorm(50), 5, 10)
> apply(x, 2, mean)
> apply(x, 1, sum)

Здесь для начала мы создаем матрицу из пяти строк и десяти столбцов, потом сначала считаем среднее по столбцам, а затем сумму в строках. Для полноты картины следует отметить, что задачи вычисления среднего и суммы по строкам настолько часто встречаются, что в R для этого предусмотрены специальные функции rowSums, rowMeans, colSums и colMeans.
Также функцию apply можно применять для многомерных массивов:

> arr <- array(rnorm(2 * 2 * 10), c(2, 2, 10))
> apply(arr, c(1,2), mean)

Последний вызов можно заменить на более удобный для чтения вариант:

> rowMeans(arr, dim = 2)

Перейдем к функции mapply, представляющей собой многомерный аналог lapply. Начнем с простого примера, который можно найти прямо в стандартной документации к R:

> mapply(rep, 1:4, 4:1)
[[1]]
[1] 1 1 1 1

[[2]]
[1] 2 2 2

[[3]]
[1] 3 3

[[4]]
[1] 4

Как можно видеть, здесь происходит применение функции rep к набору параметров, которые генерируются из двух последовательностей. Сама функция rep просто повторяет первый аргумент то число раз, которое указано в качестве второго аргумента. Таким образом, предыдущий код просто эквивалентен следующему:

> list(rep(1,4), rep(2,3), rep(3,2), rep(4,1))

Иногда бывает необходимо применить функцию к какой-то части массива. Это можно сделать с помощью функции tapply. Давай рассмотрим следующий пример:

> x <- c(rnorm(10, 1), runif(10), rnorm(10,2))
> f <- gl(3,10)
> tapply(x,f,mean)

Сначала мы создаем вектор, части которого формируются из случайных величин с различным распределением, далее мы генерируем вектор из факторов, который представляет собой не что иное, как десять единиц, потом десять двоек и столько же троек. Затем вычисляем среднее по соответствующим группам. Функция tapply по умолчанию пытается упростить результат. Эту опцию можно выключить, указав в качестве параметра simplify=FALSE.

> tapply(x, f, range, simplify=FALSE)

Когда говорят о функциях apply, обычно также говорят о функции split, которая разбивает вектор на части, аналогично tapply. Так, если мы вызовем split(x, f) то получим список из трех векторов. Таким образом, пара lapply/split работает так же, как и tapply со значением simplify, равным FALSE:

> lapply(split(x, f), mean)

Функция split полезна и за пределами работы с векторами: ее также можно использовать и для работы с фреймами данных. Рассмотрим следующий пример (я позаимствовал его из курса R Programming на Coursera):

> library(datasets)
> head(airquality)
    Ozone Solar.R Wind Temp Month Day
1    41     190  7.4   67     5   1
2    36     118  8.0   72     5   2
3    12     149 12.6   74     5   3
4    18     313 11.5   62     5   4
5    NA      NA 14.3   56     5   5
6    28      NA 14.9   66     5   6

> s <- split(airquality, airquality$Month)
> lapply(s, function(x) colMeans(x[, c("Ozone", "Solar.R", "Wind")]))

Здесь мы работаем с набором данных, который содержит информацию о состоянии воздуха (содержание озона, солнечная радиация, ветер, температура в градусах Фаренгейта, месяц и день). Мы можем легко сделать отчет о среднемесячных показателях, используя split и lapply, как показано в коде. Использование sapply, однако, даст нам результат в более удобном виде:

> sapply(s, function(x) colMeans(x[, c("Ozone", "Solar.R", "Wind")]))
               5         6          7        8        9
Ozone         NA        NA         NA       NA       NA
Solar.R       NA 190.16667 216.483871       NA 167.4333
Wind    11.62258  10.26667   8.941935 8.793548  10.1800

Как видно, некоторые значения величин не определены (и для этого используется зарезервированное значение NA). Это означает, что какие-то (хотя бы одно) значения в колонках Ozone и Solar.R были также не определены. В этом смысле функция colMeans ведет себя совершенно корректно: если есть какие-то неопределенные значения, то и среднее, таким образом, не определено. Проблему можно решить, принудив функцию не учитывать значения NA с помощью параметра na.rm=TRUE:

> sapply(s, function(x) colMeans(x[, c("Ozone", "Solar.R", "Wind")], na.rm=TRUE))
                5         6          7          8         9
Ozone    23.61538  29.44444  59.115385  59.961538  31.44828
Solar.R 181.29630 190.16667 216.483871 171.857143 167.43333
Wind     11.62258  10.26667   8.941935   8.793548  10.18000

Зачем нужно такое количество функций для решения очень схожих между собой задач? Думаю, такой вопрос задаст каждый второй человек, который все это прочитал. Все эти функции на самом деле пытаются решить проблему обработки векторных данных без использования циклов. Но одно дело — добиться высокой скорости обработки данных и совершенно другое — получить хотя бы часть той гибкости и контроля, которую обеспечивают такие управляющие конструкции, как циклы и условные операторы.

Визуализация данных

Система R необыкновенно богата на средства визуализации данных. И тут передо мной стоит непростой выбор — о чем вообще говорить, если область так велика. Если в случае программирования есть какой-то базовый набор функций, без которого ничего не сделать, то в визуализации огромное количество различных задач и каждая из них (как правило) может быть решена несколькими способами, каждый из которых имеет свои плюсы и минусы. Более того, всегда есть множество опций и пакетов, позволяющих решать эти задачи различным образом.
Про стандартные средства визуализации в R написано очень много, поэтому здесь мне бы хотелось рассказать о чем-то более интересном. В последние годы все более популярным становится пакет ggplot2, вот про него и поговорим.

Для того чтобы начать работать с ggplot2, нужно установить библиотеку с помощью команды install.package("ggplot2"). Далее подключаем ее для использования:

> library("ggplot2")
> head(diamonds)
  carat       cut color clarity depth table price    x    y    z
1  0.23     Ideal     E     SI2  61.5    55   326 3.95 3.98 2.43
2  0.21   Premium     E     SI1  59.8    61   326 3.89 3.84 2.31
3  0.23      Good     E     VS1  56.9    65   327 4.05 4.07 2.31
4  0.29   Premium     I     VS2  62.4    58   334 4.20 4.23 2.63
5  0.31      Good     J     SI2  63.3    58   335 4.34 4.35 2.75
6  0.24 Very Good     J    VVS2  62.8    57   336 3.94 3.96 2.48
> head(mtcars)
                   mpg cyl disp  hp drat    wt  qsec vs am gear carb
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

Данные diamonds и mtcars являются частью пакета ggplot2, и именно с ними мы будем сейчас работать. С первым все понятно — это данные о бриллиантах (чистота, цвет, стоимость и прочее), а второй сет — это данные дорожных тестов (количество миль на галлон, количество цилиндров...) автомобилей 1973–1974 годов выпуска из американского журнала Motor Trends. Более подробную информацию о данных (к примеру, размерность) можно получить, набрав ?diamonds или ?mtcars.

Для визуализации в пакете предусмотрено много функций, из которых для нас сейчас будет наиболее важна qplot. Функция ggplot предоставляет существенно больше контроля над процессом. Все, что можно сделать с помощью qplot, также можно сделать и с помощью ggplot. Рассмотрим это на простом примере:

> qplot(clarity, data=diamonds, fill=cut, geom="bar")

Того же самого эффекта можно достичь и функцией ggplot:

> ggplot(diamonds, aes(clarity, fill=cut)) + geom_bar()

Однако вызов qplot выглядит проще. На рис. 1 можно увидеть, как строится зависимость количества бриллиантов с различным качеством огранки (cut) от чистоты (clarity).

Рис. 1. Зависимость количества бриллиантов разного качества огранки от чистоты
Рис. 1. Зависимость количества бриллиантов разного качества огранки от чистоты

Теперь построим зависимость пробега на единицу топлива автомобилей от их массы. Полученная точечная диаграмма (или диаграмма рассеивания scatter plot) представлена
на рис. 2.

> qplot(wt, mpg, data=mtcars)

Можно также еще добавить цветовое отображение показателя времени разгона на четверть мили (qsec):

> qplot(wt, mpg, data=mtcars, color=qsec)

Рис. 2. Диаграмма зависимости экономичности автомобиля от массы
Рис. 2. Диаграмма зависимости экономичности автомобиля от массы

При визуализации также можно преобразовывать данные:

> qplot(log(wt), mpg - 10, data=mtcars)

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

> qplot(wt, mpg, data=mtcars, color=factor(cyl))

Рис. 3. Диаграмма зависимости экономичности автомобиля от массы с отображением информации о количестве цилиндров
Рис. 3. Диаграмма зависимости экономичности автомобиля от массы с отображением информации о количестве цилиндров

Также можно менять размер точек, используя, к примеру, size=3. Если ты собираешься печатать графики на черно-белом принтере, то лучше не использовать цвета, а вместо этого менять форму маркера в зависимости от фактора. Это можно сделать, заменив color=factor(cyl) на shape=factor(cyl).
Тип графика задается с помощь параметра geom, и в случае точечных диаграмм значение этого параметра равно "points".

Теперь пускай мы хотим просто построить гистограмму по количеству автомобилей с соответствующим значением цилиндров:

> qplot(factor(cyl), data=mtcars, geom="bar")
> qplot(factor(cyl), data=mtcars, geom="bar", color=factor(cyl))
> qplot(factor(cyl), data=mtcars, geom="bar", fill=factor(cyl))

Первый вызов просто рисует три гистограммы для различных значений цилиндров. Надо сказать, что первая попытка придать цвет гистограмме не приведет к ожидаемому результату — черные столбики так и будут черными, только получат цветной контур. А вот последний вызов qplot сделает красивую гистограмму, как показано на рис. 4.

Тут следует внести ясность. Дело в том, что текущий построенный нами объект не является гистограммой в строгом смысле слова. Обычно под гистограммой понимают аналогичное отображение для непрерывных данных. В английском языке bar chart (это то, что мы только что сделали) и histogram — это два различных понятия (см. соответствующие статьи в Википедии). Здесь я, с некоторой тяжестью на душе, буду использовать слово «гистограмма» для обоих понятий, полагая, что сама природа данных говорит за себя.

Рис. 4. Распределение автомобилей по числу цилиндров в нашем наборе данных
Рис. 4. Распределение автомобилей по числу цилиндров в нашем наборе данных

Если вернуться к рис. 1, то ggplot2 предоставляет несколько полезных опций в позиционирование графиков (по умолчанию используется значение position="stack"):

> qplot(clarity, data=diamonds, geom="bar", fill=cut, position="dodge")
> qplot(clarity, data=diamonds, geom="bar", fill=cut, position="fill")
> qplot(clarity, data=diamonds, geom="bar", fill=cut, position="identity")

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

Рис. 5. Сравнительный анализ качества огранки бриллиантов с заданной чистотой
Рис. 5. Сравнительный анализ качества огранки бриллиантов с заданной чистотой

Рис. 6. Доля бриллиантов с различным качеством огранки в общем объеме бриллиантов с заданной чистотой
Рис. 6. Доля бриллиантов с различным качеством огранки в общем объеме бриллиантов с заданной чистотой

Теперь рассмотрим пример настоящей гистограммы:

> qplot(carat, data=diamonds, geom="histogram", bandwidth=0.1)
> qplot(carat, data=diamonds, geom="histogram", bandwidth=0.05)

Здесь параметр bandwidth как раз показывает, какой ширины полоса в гистограмме. Гистограмма показывает, сколько данных приходится на какой диапазон. Результаты представлены на рис. 7 и 8.

Рис. 7. Гистограмма распределения бриллиантов по массе для значения bandwidth=0.1
Рис. 7. Гистограмма распределения бриллиантов по массе для значения bandwidth=0.1

Рис. 8. Гистограмма распределения бриллиантов по массе для значения bandwidth=0.05
Рис. 8. Гистограмма распределения бриллиантов по массе для значения bandwidth=0.05

Иногда, когда нам нужно построить модель (линейную или, скажем, полиномиальную), мы можем сделать это прямо в qplot и увидеть результат. К примеру, мы можем построить график зависимости mpg от массы wt прямо поверх точечной диаграммы:

> qplot(wt, mpg, data=mtcars, geom=c("point", "smooth"))

По умолчанию в качестве модели будет использоваться локальная полиномиальная регрессия (method="loess"). Результат работы будет выглядеть, как показано на рис. 9, где темно-серая полоса — это стандартная ошибка. Она отображается по умолчанию, ее отображение можно выключить, написав se=FALSE.

Рис. 9. График, полученный из точечной диаграммы с помощью локальной полиномиальной регрессии. На графике также отображается стандартная ошибка
Рис. 9. График, полученный из точечной диаграммы с помощью локальной полиномиальной регрессии. На графике также отображается стандартная ошибка

Если же мы хотим попробовать натянуть линейную модель на эти данные, то это можно сделать, просто указав method=lm (рис. 10).

Рис. 10. Линейная регрессия на основе данных из предыдущей диаграммы
Рис. 10. Линейная регрессия на основе данных из предыдущей диаграммы

И напоследок, конечно же, нужно показать, как строить круговые диаграммы:

> t <- ggplot(mtcars, aes(x=factor(1), fill=factor(cyl))) + geom_bar(width=1)
> t + coord_polar(theta="y")

Здесь мы воспользуемся более гибкой функцией ggplot. Это работает так: сначала мы строим график, отображающий доли автомобилей с разным количеством цилиндров в общей массе (рис. 11), затем переводим график в полярные координаты (рис. 12).

Рис. 11. Промежуточный вариант из которого мы планируем получить круговую диаграмму
Рис. 11. Промежуточный вариант из которого мы планируем получить круговую диаграмму

Рис. 12. Финальный вариант круговой диаграммы
Рис. 12. Финальный вариант круговой диаграммы

Вместо заключения

Вот мы и освоились с использованием R. Что дальше? Понятно, что здесь даны самые базовые возможности ggplot2 и рассмотрены вопросы, связанные с векторизацией. Есть несколько хороших книг по R, которые стоит упомянуть, и к ним, вне всякого сомнения, стоит обращаться чаще, чем к услугам корпорации очень навязчивого добра. Во-первых, это книга Нормана Матлоффа (Norman Matloff) The Art of R Programming. Если же у тебя уже есть опыт в программировании на R, то тебе пригодится The R Inferno, написанная П. Бернсом (Patrick Burns). Классическая книга Software for Data Analysis Джона Чамберса (John Chambers) также вполне уместна.

Если говорить о визуализации в R, то есть хорошая книга R Graphics Cookbook В. Чанга (Winston Chang). Примеры для ggplot2 в этой статье были взяты из Tutorial: ggplot2. До встречи в следующей статье «Анализ данных и машинное обучение в R»!

1 комментарий

  1. om52

    20.02.2016 at 19:01

    Виталий!
    полезные статьи!
    требуем продолжения!)

Оставить мнение

Check Also

Google как средство взлома. Разбираем актуальные рецепты Google Dork Queries

Тесты на проникновение обычно требуют набора специальных утилит, но одна из них доступна к…