Матрунич Консалтинг

Анализ данных, визуализация, маркетинговые исследования, язык R

Tag: R (Page 1 of 2)

Визуальная диагностика простой линейной регрессии с ggvis

Продолжаю тему с однофакторым регрессионным анализом в R. Теперь мне интересно применить инструменты визуализации, доступные в R-библиотеке ggvis, для диагностики регрессионной модели. Здесь же мы применяем перенаправление с помощью функции %>%.

Напомню, что я использую данные, сопровожающие книгу An Introduction to Statistical Learning with Applications in R. Подготовка данных и построение модели было произведено в предыдущей публикации. Так что будем считать, что мы продолжаем творить в пространстве той публикации.

Нелинейность взаимосвязи

data.frame(predict = predict(autolm),
           resids  = residuals(autolm)) %>%
  ggvis(~predict, ~resids) %>% 
  layer_points() %>%
  add_axis("x", title = "Предсказанные значения") %>%
  add_axis("y", title = "Остатки") %>%
  layer_smooths(stroke := 'red')
Точечный график со сглаживанием в ggvis

Точечный график со сглаживанием в ggvis

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

Корреляция остатков

data.frame(id = seq_len(dim(Auto)[1]),
           resids  = residuals(autolm)) %>%
  ggvis(~id, ~resids) %>% 
  layer_paths() %>%
  add_axis("x", title = "Порядковый номер наблюдения") %>%
  add_axis("y", title = "Остатки")
Визуальная оценка взаимосвязи остатков

Визуальная оценка взаимосвязи остатков

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

Непостоянная дисперсия остатков

data.frame(predict = predict(autolm),
           rstuds  = rstudent(autolm)) %>%
  ggvis(~predict, ~rstuds) %>% 
  layer_points() %>%
  add_axis("x", title = "Предсказанные значения") %>%
  add_axis("y", title = "Стьюдентизированные остатки")
Стьюдентизированные остатки и предсказанные значения: точечный график в ggvis

Стьюдентизированные остатки и предсказанные значения: точечный график в ggvis

Видно, что с увеличением величины предсказанных значений происходит увеличение разброса остатков: мы наблюдаем изогнутую воронкообразную форму. Это говорит о наличии гетероскедастичности. Отсутствие гетероскедастичности — условие для корректного расчёта стандартной ошибки, доверительных интервалов и тестирования гипотез. В нашем случае можно преобразовать зависимую переменную с помощью какой-нибудь “выгнутой” функции — например, логарифм или корень квадратный.

autolm1 <- update(autolm, log(.) ~ .)
data.frame(predict = predict(autolm1),
           rstuds  = rstudent(autolm1)) %>%
  ggvis(~predict, ~rstuds) %>% 
  layer_points() %>%
  add_axis("x", title = "Предсказанные значения (log)") %>%
  add_axis("y", title = "Стьюдентизированные остатки")
Точечный график стьюдентизированных остатков и натурального логарифма предсказанных значений

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

Выбросы

Auto %>%
  mutate(rstud = ifelse(abs(rstudent(autolm)) > 3, 'Больше 3', 'Не больше 3'),
         label = row.names(Auto),
         label = ifelse(rstud == 'Больше 3', label, '')) %>% 
  ggvis(~kpl, ~ horsepower, fill = ~rstud) %>%
  layer_points() %>%
  add_axis("x", title = "Километров на литр") %>%
  add_axis("y", title = "Лошадиных сил") %>%
  add_legend("fill", title = "Стьюд.остатки") %>%
  layer_model_predictions(stroke := "red", model = 'lm') %>%
  layer_text(text := ~label, dy := -5, dx := 5)
Визуальное определение предикатов-выбросов

Визуальное определение предикатов-выбросов

Только у двух наблюдений абсолютная величина стьюдентизированных остатков превышает 3.

Точки с избыточным влиянием на модель (High Leverage Points)

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

data.frame(leverage = hatvalues(autolm),
           rstuds  = rstudent(autolm),
           label = row.names(Auto)) %>%
  mutate(label = ifelse(leverage > .025, label, '')) %>% 
  ggvis(~leverage, ~rstuds) %>% 
  layer_points() %>%
  add_axis("x", title = "Leverage") %>%
  add_axis("y", title = "Стьюдентизированные остатки") %>%
  layer_text(text := ~label, dy := -5, dx := 5)
Точечный график с оценкой влияния наблюдений на модель

Точечный график с оценкой влияния наблюдений на модель

Обогащённая графика с диагностикой модели

data.frame(predict = predict(autolm),
           resids  = residuals(autolm),
           rstuds  = rstudent(autolm),
           leverage = hatvalues(autolm)) %>%
  ggvis(~predict, ~resids, fill = ~abs(rstuds), size = ~leverage) %>% 
  layer_points() %>%
  add_axis("x", title = "Предсказанные значения") %>%
  add_axis("y", title = "Остатки")  %>%
  add_legend("size", title = "Влияние на модель") %>% 
  add_legend("fill", title = "Модуль стьюд.ост.", 
             properties = legend_props(legend = list(y = 100)))
График, обобщающий оценку модели: выбросы, влияние отдельных значений, нелинейность взаимосвязи

График, обобщающий оценку модели: выбросы, влияние отдельных значений, нелинейность взаимосвязи

Чем больше точка, тем сильнее её влияние на модель. Чем ярче точка, тем больше у этого наблюдения ошибка предсказания. Закраску точки можно и убрать, поскольку величина остатков коррелирует с величиной стьюдентизированных остатков.

Отчётливо видно, что наблюдения с наибольшим влиянием встречаются в том секторе, где наблюдения встречаются реже.

Однофакторный линейный регрессионный анализ в R

На основе книги An Introduction to Statistical Learning with Applications in R авторов Gareth James, Daniela Witten, Trevor Hastie и Robert Tibshirani. Бесплатная копия книги в PDF-формате, или заказ бумажной копии книги на Amazon.

library(ISLR)
library(ggvis)
library(dplyr)

Подготовка данных

Работаем с массивом Auto из пакета ISLR.

# Загружаем массив Auto в пространство пользователя
data(Auto)

# Структура массива
str(Auto)
## 'data.frame':    392 obs. of  9 variables:
##  $ mpg         : num  18 15 18 16 17 15 14 14 14 15 ...
##  $ cylinders   : num  8 8 8 8 8 8 8 8 8 8 ...
##  $ displacement: num  307 350 318 304 302 429 454 440 455 390 ...
##  $ horsepower  : num  130 165 150 150 140 198 220 215 225 190 ...
##  $ weight      : num  3504 3693 3436 3433 3449 ...
##  $ acceleration: num  12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
##  $ year        : num  70 70 70 70 70 70 70 70 70 70 ...
##  $ origin      : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ name        : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...
# Первые строки массива
head(Auto)
##   mpg cylinders displacement horsepower weight acceleration year origin
## 1  18         8          307        130   3504         12.0   70      1
## 2  15         8          350        165   3693         11.5   70      1
## 3  18         8          318        150   3436         11.0   70      1
## 4  16         8          304        150   3433         12.0   70      1
## 5  17         8          302        140   3449         10.5   70      1
## 6  15         8          429        198   4341         10.0   70      1
##                        name
## 1 chevrolet chevelle malibu
## 2         buick skylark 320
## 3        plymouth satellite
## 4             amc rebel sst
## 5               ford torino
## 6          ford galaxie 500
# Конвертируем миль (1.609 км) на галлон (3.79 л) в км на литр
Auto$kpl <- Auto$mpg * 1.609 / 3.79

Построение модели

Построим однофакторную регрессионную модель, чтобы выявлить, как мощность двигателя влияет на пробег на километр топлива, и ответим на следующие вопросы:

  1. Существует ли взаимосвязь между мощностью и пробегом.
  2. Насколько сильна взаимосвязь между мощностью и пробегом.
  3. Взаимосвязь между мощностью и пробегом прямая или обратная.
  4. Каково прогнозное значение пробега на километр при мощности 98 лошадиных сил. Каковы 95% доверительный интервал и интервал предсказания.
# Задаем модель
autolm <- lm(kpl ~ horsepower, data = Auto)

# Результаты анализа
summary(autolm)
## 
## Call:
## lm(formula = kpl ~ horsepower, data = Auto)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.761 -1.384 -0.146  1.173  7.185 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 16.95430    0.30461    55.7   <2e-16 ***
## horsepower  -0.06701    0.00274   -24.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.08 on 390 degrees of freedom
## Multiple R-squared:  0.606,  Adjusted R-squared:  0.605 
## F-statistic:  600 on 1 and 390 DF,  p-value: <2e-16

Коэффициент horsepower составляет -0.067, т.е. при увеличении мощности автомобиля на одну лошадиную силу пробег снижается на 0.07 км на литр топлива. Наличие этой взаимосвязи подтверждается p-значением t-статистики: t равен -24.5, а вероятность получить такое значение t при независимых переменных (т.е. p-значение) стремится к нулю (<2e-16).

Прогнозирование

# Прогнозируем пробег при мощности 98 л.с.
# Доверительный интервал
confint <- predict(autolm,
                   newdata = data.frame(horsepower = 98),
                   interval = 'confidence')
confint
##     fit   lwr  upr
## 1 10.39 10.18 10.6

Если мы возьмём несколько разных автомобилей с мощностью двигателя 98 л.с., то средний пробег этих машин на литр топлива с вероятностью 95% попадёт в интервал от 10.2 до 10.6 км.

# Интервал предсказания
predint <- predict(autolm,
                   newdata = data.frame(horsepower = 98),
                   interval = 'prediction')
predint
##     fit   lwr   upr
## 1 10.39 6.287 14.49

Если мы возьмём какой-то конкретный автомобиль с мощностью двигателя 98 л.с., то его пробег на литр топлива с вероятностью 95% попадёт в интервал от 6.3 до 14.5 км.

Auto %>% 
  ggvis(~horsepower, ~kpl) %>%
  layer_points() %>%
  layer_model_predictions(model = 'lm', se = T)

Точечный график и линия тренда

Диагностика

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

  1. Нелинейная взаимосвязь между независимыми и зависимой переменными.
  2. Корреляция остатков между собой.
  3. Непостоянная дисперсия остатков (гетероскедастичность).
  4. Наличие выбросов.
  5. Избыточное влияние отдельных наблюдений.
  6. Коллинеарность предикатов.
par(mfrow = c(2, 3))
plot(autolm, which = 1:6, caption = 
       list("Прогноз и остатки", "Квантили", "Разброс-Положение",
            "Расстояние Кука", "Влияние и остатки", 
            expression("Влияние и Кук " * h[ii] / (1 - h[ii]))))
Визуальная диагностика регрессионной модели

Визуальная диагностика регрессионной модели

Точечный график предсказанных значений и остатков позволяет обнаружить закономерности распределения последних. График квантилей помогает определить нормальность остатков. График “Разброс-положение” даёт возможность оценить монотонность распределения остатков (гетероскедастичность). С помощью графика “Остатки-влияние” можно определить точки с большим влиянием на модель.

Что такое грамматика графики ggplot2 (ч.3)

Вы открыли продолжение статьи об основах ggplot2. Читайте начало в

  1. Что такое грамматика графики ggplot2 (ч.1)
  2. ggplot2: сопоставления, визуальные средства, статистические преобразования и слои

Подписи осей графика в ggplot2

Снова вернёмся к начальному графику из первой части:

fig <- ggplot(data = cars, aes(x = speed, y = dist)) + geom_point()
fig

ggplot(data = cars, aes(x = speed, y = dist)) + geom_point()

Обратите внимание, здесь я сохранил объект класса ggplot в переменной fig. То есть вся информация о нашем графике теперь хранится в этой переменной. Во второй строке я вызываю эту переменную, и здесь, собственно, происходит отрисовка.

Ggplot2 автоматически взял названия для осей из имён переменных, которые мы использовали для построения графика. Авторы ggplot2 предлагают сразу несколько вариантов того, как можно изменить названия осей. Можно использовать функции ylab(), xlab():

fig + xlab("Скорость, миль в час") +
  ylab("Тормозной путь, футы")
ggplot(data = cars, aes(x = speed, y = dist)) + geom_point() + xlab("Скорость, миль в час") + ylab("Тормозной путь, футы")

ggplot(data = cars, aes(x = speed, y = dist)) + geom_point() + xlab("Скорость, миль в час") + ylab("Тормозной путь, футы")

Аналогично можно использовать функцию labs():

fig + labs(x = "Скорость, миль в час",
           y = "Тормозной путь, футы")

Здесь мы получим аналогичный график, поэтому не привожу его.

Следующий вариант будет удобен, если помимо названия вы хотите сделать в оси дополнительные изменения: например, изменить расположение основных и вспомогательных линий сетки, задать собственный диапазон значений оси. В нашем случае функция для определения этих параметров называется scale_x_continuous(). К чему такое сложное название функции? Давайте вспомним про то, что пакет ggplot2 — это не очередное R-расширение для создания диаграмм, а целая концепция графической грамматики. И с этой точки зрения, подпись горизонтальной оси — это не просто текстовая строка, расположенная под абсциссой, а элемент инструкции (guide), которая призвана помочь человеку, смотрящему на график, правильно понять, что за информация передаётся с помощью горизонтальной оси. Поэтому логично поместить все параметры одного элемента графика в одну функцию. А функции labs() и ylab()/xlab() — это лишь удобные ярлыки, созданные для уменьшения нагрузки на суставы ваших пальцев и клавиатуру.

Итак, пересилим себя и зададим названия осей с помощью функций scale_*:

fig + scale_x_continuous(name = "Скорость, миль в час") +
  scale_y_continuous(name = "Тормозной путь, футы")

Расширим наше использование scale_* и «поиграем» с характеристиками осей.

fig + scale_x_continuous(name = "Скорость, миль в час",
                         breaks = c(7, 10, 20),
                         limits = c(6, 22)) +
  scale_y_continuous(name = "Тормозной путь, футы",
                     limits = c(5, 50),
                     breaks = NULL)
Точечный график ggplot2 с модификацией осей через scale_continuous()

Точечный график ggplot2 с модификацией осей через scale_continuous()

Здесь по горизонатальной оси мы «обрезали» края графика, выходящие за пределы менее 6 и более 22 миль в час, и указали, чтобы линии сетки проходили только по точкам 7, 10 и 20 миль в час. Аналогичная операция проведена для вертикальной оси, только целиком убраны горизонтальные линии сетки и соответствующие подписи.

Работа с легендами графика

Добавим к графику цвет.

fig2 <- ggplot(data = cars, aes(speed, dist, color = speed)) + geom_point()
fig2
Точечный график ggplot2 с окраской точек

Точечный график ggplot2 с окраской точек

Обратите внимание, я не стал в функции aes() указывать названия аргументов x и y. Так можно делать, так как, если вы не указываете в aes() имена первых двух аргументов, то aes() «думает», что первым аргументом является положение по горизонтали, а вторым — по вертикали.

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

fig2 + scale_color_continuous("Скорость") #Снова не указываю имя аргумента
Точечный крафик ggplot2 с исправленным названием цветовой легенды

Точечный крафик ggplot2 с исправленным названием цветовой легенды

Название графика

Название графика можно добавить одной из двух функций: ggtitle() или labs().

fig2 + scale_color_continuous("Скорость") +
  ggtitle("Влияние скорости автомобиля\nна длину тормозного пути") +
  xlab("Скорость, миль в час") +
  ylab("Тормозной путь, футов")
График ggplot2 с названием, именованными осями и легендой

График ggplot2 с названием, именованными осями и легендой

Дополнительную информацию о настройке внешнего вида графика можно почерпнуть на англоязычном сайте одного из разработчиков ggplot2 Уинстона Чанга Cookbook for R

Вы прочитали продолжение статьи. Начало читайте в:

  1. Что такое грамматика графики ggplot2 (ч.1)
  2. ggplot2: сопоставления, визуальные средства, статистические преобразования и слои

Что такое грамматика графики ggplot2 (ч.2)

Продолжение. Читайте начало в Что такое грамматика графики ggplot2 (ч.1)

Читайте продолжение в Названия элементов графика ggplot2: название графика, осей, легенд

Визуальные средства и геометрические объекты

Визуальные средства для отрисовки геометрических объектов типа «точка» (point) — это расположение по горизонтали (x) и вертикали (y), цвет границы (colour) и цвет заливки (fill), размер (size), форма (shape), прозрачность (alpha). Обязательными для создания точечного графика являются только первые два параметра: положение по вертикали и горизонтали. В примере с автомобилями мы привязали скорость к положению по горизонтали, а тормозной путь — к положению по вертикали. Оставшиеся визуальные средства можно было бы использовать для передачи на графике другой информации, если бы она содержалась в исходном массиве данных cars. Или мы можем продублировать информацию, передаваемую через положение, другими визуальными средствами. «Нагрузим» цвет точки скоростью.

ggplot(data = cars, aes(x = speed, y = dist, color = speed)) + geom_point()
ggplot(data = cars, aes(x = speed, y = dist, color = speed)) + geom_point()

ggplot(data = cars, aes(x = speed, y = dist, color = speed)) + geom_point()

Алгоритм определил, что speed — это числовая (интервальная) переменная и, по умолчанию для этого типа, предложил палитру в виде градации от тёмно-синего до голубого. Чем выше скорость, тем светлее точка. Также, поскольку наш график «обзавёлся» дополнительной информацией, передаваемой через цвет, то алгоритм сам добавил в картинку легенду со значениями цвета. В наших силах также выбрать другую цветовую палитру, дать другое название легенде и разместить легенду в другом месте графика, но сейчас не будем на это отвлекаться.

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

ggplot(data = cars, aes(x = speed, 
                        y = dist, 
                        color = speed, 
                        size = dist)) + geom_point()
ggplot(data = cars, aes(x = speed, y = dist, color = speed, size = dist)) + geom_point()

ggplot(data = cars, aes(x = speed, y = dist, color = speed, size = dist)) + geom_point()

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

Статистические преобразования

Вернёмся к самому первому графику из первой части:

ggplot(data = cars, aes(x = speed, y = dist)) + geom_point()

ggplot(data = cars, aes(x = speed, y = dist)) + geom_point()

В этом примере схема построения графика на первый взгляд выглядит так:

  1. Берётся массив данных cars.
  2. Переменная speed сопоставляется положению по горизонтали (абсциссе), переменная dist — положению по вертикали (ординате).
  3. На графике отображаются точки с соответствующими координатами.

Между вторым и третим пунктом остался незамеченным важный элемент: статистическое преобразование данных. Подход ggplot2 подразумевает, что для отрисовки графического объекта исходные данные нужно преобразовать. Например, для перед построением гистограммы алгоритм соответствующего статистического преобразования должен посчитать, сколько наблюдений попадает в каждый интервал, а затем передать данные о границах каждого интервала и о количестве элементов в алгоритм построения геометрического объекта. В нашем случае, когда рисуются точки, данные никак преобразовывать не нужно, поэтому на этапе преобразования используется stat_identity() — тождественное преобразование, то есть без изменений.

Взглянем на распределение длин тормозных путей и построим график, требующий более сложных вычислений: «Ящик с усами». Он покажет нам информацию о медиане, квартилях и наличии выбросов. Поскольку график такого типа удобно использовать для сравнения групп, то мы разделим наш массив на две группы: в первой автомобили, у которых скорость 15 миль в час и меньше, а во второй — у кого больше.

ggplot(data = cars, aes(x = speed > 15, y = dist)) + geom_boxplot()
ggplot(data = cars, aes(x = speed > 15, y = dist)) + geom_boxplot()

ggplot(data = cars, aes(x = speed > 15, y = dist)) + geom_boxplot()

Более полная схема построения графика следующая:

  1. Берётся массив данных cars.
  2. Результат выполнения выражения cars$speed > 15 сопоставляется абсциссе, переменная dist — ординате. При этом абсцисса (x) из числовой переменной (мили в час) превратилась в категориальную (логическую: больше 15 миль в час или нет), поэтому все наблюдения разбились на эти две категории.
  3. Данные передаются статистическому преобразованию stat_boxplot(), которое из предыдущего шага знает, что cars$speed > 15 нужно использовать для деления массива на группы, а dist — для расчёта статистик, необходимых для отрисовки «ящиков с усами».
  4. Результаты работы stat_boxplot(), а это для каждой группы медиана, квартили, границы нижнего и верхнего «усов», передаются алгоритму отрисовки.

Нам не понадобилось явно задавать stat_boxplot() в коде, т.к. это статистическое преобразование вызывается по умолчанию при отрисовке «ящиков с усами». Явно указывать статистику имело бы смысл в такой ситуации, когда наш массив уже содержит все необходимые данные для geom_boxplot() (медиану, квартили и пр.), и мы бы захотели изменить способ статистического преобразования: geom_boxplot(stat="identity")

Слои графика

Каждый раз, когда вы в R добавляете к своему ggplot2-объекту новые geom_*() или stat_*(), на вашем графике появляется новый слой. Добавим к нашим «ящикам» слой с точками.

ggplot(data = cars, aes(x = speed > 15, y = dist)) +
  geom_boxplot() + 
  geom_point()
ggplot(data = cars, aes(x = speed > 15, y = dist)) + geom_boxplot() + geom_point()

ggplot(data = cars, aes(x = speed > 15, y = dist)) + geom_boxplot() + geom_point()

Здесь сначала по старой схеме были отрисованы «ящики», а затем на основе тех же самых данных, тех же самых сопоставлений добавлены точки.

Читайте первую часть статьи в Что такое грамматика графики ggplot2 (ч.1)

Читайте продолжение статьи в Названия элементов графика ggplot2: название графика, осей, легенд

Оптимизация расходов на услуги сотовой связи (черновик)

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

Управленческая проблема

Я задался вопросом, схожим с тем, который мучил интервьюеров, когда решил сменить мой текущий тарифный план «Дело чести» от «Теле2» на что-нибудь свеженькое. Принимать решение «на основе интуиции» не хочется, когда под боком есть статистические данные, и можно принять data-driven decision («решение на основе данных»). Подобная задача может стать перед организацией, которая хочет оптимизировать расходы на связь подбором более выгодных тарифов.

Я решил создать пробный алгоритм на языке R для выбора наиболее выгодного тарифа на основе имеющейся истории использования услуг сотовой связи и рассмотреть перспективы создания сервиса, предлагающего людям и компаниям услугу оптимизации расходов на связь.

Программа исследования

Проблема и цель исследования

Минимум
Проблема заключается в том, что клиенту неизвестно, какой из предлагаемых тарифов позволит понести наименьшие затраты. Цель — определить такой тариф. Поскольку проект пилотный, а, к тому же, я не собираюсь менять своего оператора, то остановимся на тарифах только от Теле2.
Максимум
Непонятно, какие технические трудности могут появиться при реализации идеи по подбору тарифа в виде вэб-сервиса. Цель — в процессе работы над пилотным проектом рассмотреть возможности его переноса на бизнес-рельсы.

Объект, предмет и выборочная совокупность

Объектом исследования становится использование мной услуг сотовой связи. Предмет — различия в среднемесячной стоимости услуг при их расчёте по разным тарифам Теле2.

В качестве выборки мы возьмём мою историю звонков/сообщений и трафика за апрель 2014 года.

Схема моделирования и анализа

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

Чтобы определить выгодный тариф для моего «усреднённого» месяца, нам потребуется перенести выводы эксперимента с выборочной выборки (апрель 2014) на генеральную совокупность. Для этого мы воспользуемся методом Bootstrapping.

Вся обработка данных, расчёты и визуализация произведены на языке программирования R. Исходный код скриптов опубликован в репозитории на github.

Результаты пилотного исследования

Перспективы создания веб-сервиса или предложения услуги по выбору тарифа

Приложение: описание резализации пилотного проекта

Подготовка данных

Подготовка данных обычно занимает самую большую часть времени. Наш проект — не исключение.

Преобразование PDF в Excel

Истончиком данных послужила детализация расходов (стоимость — 15 рублей). Рассматривался вариант использования истории напрямую из телефона, но у детализации от оператора есть преимущества: в документе указываются сведения о стоимости каждой операции, и это позволяет проверить, правильно ли работает наш алгоритм, а также гораздо проще разобраться, находился ли телефон в роуминге, на какой телефон совершён вызов (на псковский Теле2 или на Теле2 в соседнем регионе и пр.). Плюс к этому, из истории телефона не вытащишь посессионные объёмы трафика (по крайней мере, из моего).

В личном кабинете Теле2 не предлагается выбрать формат, в котором клиент получает детализацию расходов. Все отчёты только в формате PDF. Я смог бесплатно преобразовать полученный файл в формат электронной таблицы с помощью онлайн-сервиса. К сожалению, я забыл его название, и в следующий раз придётся искать заново. После конвертации пришлось кое-где добавить ячеек, т.к. где-то дата и время попали в одну ячейку, а где-то в разные. Для однократной операции и небольшой истории затраты времени не слишком ощутимы, но в корпоративных случаях придётся договариваться с сотовыми провайдерами о чём-то более дружественном.

Первые строки массива после импорта в R приведены ниже. Все телефонные номера заменены на случайные.

> head(data)
                  date       plan                 type             number seconds  sum  pay
1  2014-04-01 08:42:48 Дело чести      Internet-трафик Трафик: 13915 байт      NA 0.63 0.10
2  2014-04-01 08:43:54 Дело чести      Internet-трафик Трафик: 53514 байт      NA 0.63 0.10
3  2014-04-01 09:00:05         БВ         Входящее SMS       +12468335597       0 0.00 0.00
4  2014-04-01 09:26:15 Дело чести Исходящий на Мегафон       +79211162647      14 0.35 0.35
5  2014-04-01 10:29:36         БВ         Входящее SMS               7031       0 0.00 0.00

Перенос выводов о выгодности от конкетного месяца на «усреднённый» месяц

Суть метода «бутстреп» заключается в том, что из всей моей апрельской истории (в ней 520 телефонных событий) случайным (вероятностным) образом составляется множество (несколько тысяч) выборок объёмом 520 элементов. Каждый элемент может попасть в выборку несколько раз или не попасть вообще. Т.е. мы получаем массив из нескольких тысяч выборок, которые почти одинаковы, но немного (а иногда сильно) отличаются друг от друга. Другими словами, мы создали несколько тысяч альтернатив моему «телефонному» апрелю 2014 года. Поскольку апрель, по моим ощущениям, в плоскости использования телефона не имел радикальных отличий от других месяцев (в апреле, как и в другие месяцы, я не звонил на спутниковые телефоны, не болтал часами в международном роуминге, не отправлял SMS в какие-нибудь голосования), я могу с достаточной доли уверенности сказать, что многие из имеющихся сгенерированных выборок могли также быть получены и из других месяцев. Посчитав по каждой из нескольких тысяч выборок стоимость услуг для каждого из тарифов, мы можем определить среднюю цену моего месяца по тарифу и доверительные интервалы. Это даст нам возможность сравнить тарифы, понять, есть ли между ними значимые отличия и выбрать самый выгодный.

Если описание алгоритма Bootstrap показалось вам достаточно сложным, не переживайте: он реализован в R, и нам останется лишь передать данные и параметры подсчёта программе, а затем получить результаты.

Сравнение речей президентов России, Белоруссии и Украины с помощью облака слов

Сравнение поздравительных речей Путина, Лукашенко и Турчинова на День Победы

Сравним речи, которые произнесли руководители России, Белоруссии и Украины 9 мая 2014 года.
Сделаем это с помощью визуализации текстов в виде облака слов.

Облако слов со сравнением частот слов в речах Путина, Лукашенко и Турчинова

Облако слов со сравнением частот слов в речах Путина, Лукашенко и Турчинова

Как интерпретировать график? Размер шрифта, которым изображено слово, показывает относительную частоту, с которой данное слово встречалось в речи. Цвет показывает, в речи какого из президентов данное слово занимает наибольшую долю. Например, слово «победа» (и его формы) встречалось у Лукашенко 5 раз, у Путина – 3, у Турчинова – 5. Понятно, что слово должно попасть или к Лукашенко, или к Турчинову. Но у Лукашенко слово «победа» занимает 0,89% от массива учтённых слов, а у Турчинова – 2,8%, соответственно оно попадает в украинскую часть графика.

У Лукашенко много слов с маленьким размеров шрифта – это значит, что он говорил больше, чем остальные спикеры, размер словаря его выступления больше, и на каждое слово у белорусского президента приходится меньшая относительная частота. Действительно, у Лукашенко зафиксировано 562 слова, у Путина – 229, у Турчинова – 177.

Чем схожи выступления спикеров? Это можно проверить с помощью облака сходства. Здесь приведены слова, которые встречались одновременно у всех трёх выступающих. Чем больше размер шрифта, тем выше средняя доля встречаемости слова.

Облако сходства текстов выступлений Лукашенко, Путина и Турчинова на День Победы

Облако сходства текстов выступлений Лукашенко, Путина и Турчинова на День Победы

При подготовке из текстов выступлений были удалены стоп-слова (предлоги и пр.), а оставшиеся слова были приведены к исходной словоформе с помощью программы MyStem от компании Яндекс. Поскольку определение первичной словоформы происходит в автоматическом режиме без учёта контекста, некоторые слова преобразованы некорректно. Например, слово «дуг» в речи Путина в оригинале звучало как «Дуге». Обработка данных произведена в среде статистической обработки данных R. Более подробно о создании графиков в виде облака слов можно ознакомиться в публикации Обработка естественно-языковых текстов в R: облако слов.

Источники информации

Обработка естественно-языковых текстов в R: облако слов

Обработка естественно-языковых текстов (Natural language processing) – это одна из областей, в которых применяется R.

В этой публикации вы познакомитесь с базовыми инструментами анализа, основанного на данных о частоте встречаемости слов. В частности, мы рассмотрим функции из расширений tm и wordcloud: подготовим текстовые документы для частотного анализа и сделаем на их основе облако слов.

Материалом для работы у нас послужит выборка новостных сообщений о работе Псковского городского молодёжного центра, опубликованных в течение 2013 года в различных Интернет-СМИ. Поиск публикаций и сбор отобранных текстов произвела Екатерина Банщикова в рамках работы над своей дипломным проектом. Всего из более чем сотни новостных сообщений случайным образом было выбрано 20 текстов.

Текстовый корпус

Источники формирования корпуса

Первым шагом при обработке текстов является подготовка простого корпуса текстов. В библиотеке tm поддерживается несколько источников формирования корпуса. Например, на моём компьютере после установки этой библиотеки мне стали доступны функции по созданию корпуса из массива данных (DataFrameSource), из каталога файлов, каждый из которых содержит отдельный документ (DirSource), из вектора (VectorSource), из текстов, доступных по URI-ссылке (URISource) и из набора новостных сообщений агентства Reuters, входящих в состав расширения tm (ReutersSource). Помимо типов, которые входят в библиотеку tm, пользователь может оперативно установить расширения из CRAN, позволяющие импортировать тексты из веб-источников в форматах XML, JSON, HTML (tm.plugin.webmining), из почтовых файлов eml и mbox (tm.plugin.mail), из файлов французской программы текстового анализа Alceste (tm.plugin.alceste), из файлов, полученных от новостных агрегаторов Europresse (tm.plugin.europresse), LexisNexis (tm.plugin.lexisnexis), Dow Jones Factiva (tm.plugin.factiva).

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

Чтение источника

Документы могут храниться в источнике в разных форматах. Библиотека tm поддерживает чтение документов в форматах простого текста (readPlain), PDF (readPDF), Microsoft Word (readDOC), XML (readXML), табличной структуры (readTabular – для источников DataFrameSource и VectorSource) и несколько Reuters-форматов. Наш вариант – простой текст.

library(tm)
# Файлы хранятся в каталоге articles
dirpath <- "articles"
articles <- Corpus(DirSource(dirpath),
                   readerControl = list(reader = readPlain,
                                      language = "ru",
                                      load = T))

Предполётная обработка документов корпуса

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

Библиотека tm содержит инструменты по удалению лишних пробелов, приведению всех букв к строчному виду, удалению цифр, знаков пунктуации и стоп-слов. Приятно то, что присутствует в библиотеке и список стоп-слов для русского языка (проверить его содержание можно с помощью команды stopwords("russian")).

Команда tm_map позволяет применить заданную функцию к каждому документу в корпусе.

# Текстовые файлы были сохранены в кодировке Windows. Конвертируем.
articles <- tm_map(articles, iconv, 'cp1251', 'UTF-8')
articles <- tm_map(articles, stripWhitespace)
articles <- tm_map(articles, tolower)
articles <- tm_map(articles, removeNumbers)
articles <- tm_map(articles, removeWords, stopwords("russian"))
articles <- tm_map(articles, removePunctuation)

Облако слов: первый заход

library(wordcloud)
wordcloud(articles, random.order=F, max.words=80, 
          colors=brewer.pal(6,"Oranges"))

Облако слов без стемматизации

Облако слов без стемматизации

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

Стемматизация текстов

Стемматизация по-русски, или плач о национальной свободе

Любую околонаучно-исследовательскую работу крайне желательно проводить с помощью таких компьютерных программ, чей исходный код открыт и доступен для проверки публикой. Иначе, если вы, например, производите расчёты в Excel или SPSS, выводы и воспроизводимость вашего исследования можно поставить под сомнение, т.к. алгоритмы, по которым эти расчёты производились, недоступны ни вам, ни сообществу, а только разработчику программы.

Подозреваю, что алгоритмы стемматизации, в частности, для русского языка - это вещь крайне наукоёмкая, требующая для её создания больших умственных трудозатрат. В России, где отношение к нарушению прав на интеллектуальную собственность достаточно лояльное, а уровень понимания того, что такое лицензионный договор, открытый исходный код и свободные лицензии, у нашего научного сообщества низкий. Поэтому неудивительно, что эффективной программы для стемминга русскоязычных текстов нет, а для создания, например, «Национального корпуса русского языка» с участием РАН использовалась программа с закрытым исходным кодом MyStem от Яндекса. Это при том, что на разработку корпуса выделялось несколько грантов из государственного бюджета. Выглядит это грустно-позорно. Благо, обнаружился Открытый корпус русского языка, все тексты и связанное с ним программное обеспечение распространяются под свободными лицензиями - сразу жить стало приятнее.

Что же есть доступного в R для стемминга русского языка? Библиотека SnowballC реализует стеммер Портера. Например, на наши экземпляры дублирующих слов выше мы получем следующее:

wordStem(c('молодёжный', "молодёжного", "молодёжном",
           "года", "году", "участники", "участников"), language='ru')

[1] "молодёжн" "молодёжн" "молодёжн" "год" "год" "участник" "участник"

Вроде бы ничего, но не каждый поймёт, что такое «молодёжн». Ещё один пример стеммера с открытым исходным кодом - это отечественная разработка Андрея Коваленко Stemka. Написан на C++, может быть использован в качестве библиотеки, так что энтузиасты могут сделать библиотеку для R. На аналогичную строку русских слов Stemka выдаст следующий результат:

молодёжн|ый
молодёжн|ог|о
молодёжн|ом|
год|а
год|у
участник|ам|
участник|и

Удалив в каждой строке вертикальную черту и символы справа от неё мы получим результат, аналогичный портеровскому. Кстати, слово «кровать» он превращает в «крова», как и Портер.

Существует открытый проект АОТ, согласно описаниям предлагающий «COM-объекты, которые можно скомпилировать под Visual Studio и использовать под Delphi». Судя по активности в их svn-репозитории, проект жив. Спасибо разработчикам, что помимо доступных COM-объектов, для компиляции которых нужно купить Visual Studio, Delphi и операционную систему Windows, они предоставили демо-версию своих продуктов онлайн. Мы можем воспользоваться примером морфологического анализа http://aot.ru/demo/morph.html. Данный анализатор успешно справляется и со словом «кровать», и «молодёжного» он превращает в «молодёжный». Для немассированной обработки текстов в R можно отправлять автоматические запросы через демо-версию и «вытаскивать» нужное слово из получаемых в ответ страниц с помощью XPath. Но здесь нужно подумать, насколько такие действия сочетаются с условиями использования демо-версии.

Исходя из вышесказанного, придётся грусто-позорно пользоваться программой с закрытым исходным кодом.

Стемматизация с помощью MyStem от Яндекса

Стеммер MyStem разработан в 1998 году сооснователем Яндекса Ильёй Сегаловичем. Программа доступна для различных операционных систем, исходный код закрыт. И хотя сам сайт Яндекса уверяет, что программа доступна только для некоммерческого использования, внимательно ознакомившись с лицензионным договором, в пункте 3.2 вы обнаружите, что «Программа может использоваться в коммерческих целях для разработки/создания каких-либо сервисов или программ, включаться и использоваться по прямому функциональному назначению в составе таких сервисов или программ, а также использоваться иным образом в процессе оказания услуг/выполнения работ». Исключения составляют те случаи, когда пользователь хочет с помощью MyStem заниматься спамом, раскруткой сайтов или создать свой поисковый движок. Спасибо Яндексу за это!

Правда, помимо закрытого кода у MyStem есть ещё одна ложка дёгтя: в пункте 3.11 соглашения читаем, что «при установке на персональный компьютер каждой копии Программы присваивается индивидуальный номер, который автоматически сообщается Правообладателю». Понятно, что Яндекс хочет знать, сколько копий его программы используется (ведь часть пользователей MyStem точно приносит экономический вред Яндексу, когда пользуются им для раскрутки сайтов), но всё равно неприятно.

Для использования mystem в R создадим одноименную функцию. Для работы требуется, чтобы операционная система знала, где лежит исполняемый файл mystem (в Windows, видимо, mystem.exe).

mystem <- function(doc) {
  library(stringr)
  sdoc <- system('mystem -nl -e utf-8 ', intern=T, input=doc)
  # При получении нескольких вариантов mystem разделяет их
  # вертикальной чертой. Удалим черту и варианты.
  sdoc <- str_replace(sdoc, '\\|.*$', '')
  # Если mystem сомневается в результате, он добавляет знак вопроса. Удаляем.
  sdoc <- str_replace(sdoc, '\\?', '')
  sdoc <- paste(sdoc, collapse=" ")
  attributes(sdoc) <- attributes(doc)
  sdoc
}

Теперь мы обладаем инструментом для стемминга и можем возвращаться в облака.

Финальный аккорд

Дообработка текстов корпуса

Добавляем к обработке текстов стемматизацию. Также логично удалить стоп-слова уже после приведения слов к исходной форме.

articles <- tm_map(articles, iconv, 'cp1251', 'UTF-8')
articles <- tm_map(articles, stripWhitespace)
articles <- tm_map(articles, tolower)
articles <- tm_map(articles, removeNumbers)
articles <- tm_map(articles, removePunctuation)
articles <- tm_map(articles, mystem)
articles <- tm_map(articles, removeWords, c(stopwords("russian"),
                                            "это", "также",
                                            "быть", "мочь",
                                            "май", "апрель", "март",
                                            "ноябрь", "который",
                                            "псковский", "молодежный",
                                            "псков", "центр", 
                                            "городской"))

Рисуем облако слов

wordcloud(articles, random.order=F, max.words=80, 
          colors=brewer.pal(6,"Oranges"))
Облако слов после стемматизации. Слова с низкой частотой еле видны

Облако слов после стемматизации. Слова с низкой частотой еле видны

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

В частотах весь фокус: в предыдущей версии графика частоты одного и того же слова «размазывавались» по разным словоформам, а теперь частоты скопились в базовой словоформе. И, естественно, слова-чемпионы добрали гораздо больше, чем слова-аутсайдеры. В итоге, разрыв в частотах между лидерами и отстающими увеличился, а чтобы этот разрыв показать, алгоритм увеличил разницу в размере шрифта и интенсивности цвета. Что делать?

Как правильно использовать цвета в статистике или немного грамматики графики

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

Облако слов без использования цвета. Все слова хорошо читаются

Облако слов без использования цвета. Все слова хорошо читаются

В таком случае все слова хорошо читаются. А если заказчик всё-таки хочет цветную картинку? Надо подумать...

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

articles <- tm_map(articles, removeWords, c(stopwords("russian"),
                                            "это", "также",
                                            "быть", "мочь", "октябрь",
                                            "май", "апрель", "март",
                                            "ноябрь", "который",
                                            "псковский", "молодежный",
                                            "псков", "центр", 
                                            "городской"))
wordcloud(articles, random.order=F, max.words=80)
Облако слов из новостей о Псковском городском молодёжном центре

Облако слов из новостей о Псковском городском молодёжном центре

Опубликованы слайды выступления “Язык R: анализ и визуализация данных”

Что такое грамматика графики ggplot2 (ч.1)

Начало. Продолжение в

Что такое ggplot2 и ggvis

ggplot2 – это расширение языка R, предназначенное для визуализации данных. Для создания графики ggplot2 использует систему абстрактных понятий: массив данных, визуальные средства, геометрические объекты, сопоставление переменных из массива визуальным средствам, статистическое преобразование переменных, системы координат и пр. Пакет ggplot2 является одним из самых популярных среди R-пользователей: в январе-мае 2013 года его скачали не менее 82 тысяч раз, что ставит его на третье место среди всех R-расширений (если учесть, что второе место занимает пакет digest, который используется для вычисления криптографических хеш-значений объектов в R, т.е. вряд ли он часто применяется пользователями напрямую, а гораздо чаще является зависимостью для других расширений, в т.ч. и для ggplot2, то ggplot2 переезжает на второе место). Функции ggplot2 берут на себя решение многих второстепенных вопросов (нужна ли легенда, где её разместить, какие границы выбрать для осей и пр.) и позволяют пользователю сконцентрироваться на самом важном. Для эффективного использования функционала ggplot2 пользователю необходимо понять принципы, заложенные в этом расширении.

Необходимо добавить несколько слов про ggvis: это совсем свежий R-пакет, реализующий идеи грамматики графики в интерактивном веб-ключе. Если ggplot2 вы будете использовать для создания статичных изображений, то ggvis предназначен для производства визуализаций аналогичного плана, но они будут интерактивными. В качестве бэкэнда визуализации ggvis использует vega, которая, в свою очередь, покоится на плечах D3.js, а для организации взаимодействия с пользователем ggvis использует R-расширение shiny.

Краткая предыстория создания ggplot2

ggplot2 начинался с того, что его будущий создатель Хэдли Викхэм прочитал книгу Лэланда ВилкинсонаThe Grammar of Graphics“. Среди прочего Вилкинсон – это тот товарищ, который в статистическом пакете SPSS поменял предыдущую систему визуализации на новую, основанную на оригинальном языке программирования графики GPL (не путать с названием лицензионного договора GNU General Public License!).

И в 2005 году Викхэм создаёт своё первое расширение для R – он реализует идеи, представленные в книге Вилкинсона, в R-коде пакета ggplot. Это был его первый пакет для R, и, как позже говорит сам Хэдли, написан он был не самым лучшим образом (с точки зрения позднего Хэдли 🙂 ). Сегодня Хэдли – автор уже не одного десятка разных пакетов, включая лидера из упомянутого выше рейтинга. В 2007 году Хэдли выпустил радикально переделанный ggplot и, чтобы не ломать скрипты и пакеты, написанные другими пользователями под синтаксис первой версии, он меняет название на ggplot2. В 2009 году Хэдли публикует книгу ggplot2: Elegant Graphics for Data Analysis. В 2012 году выходит очередная версия пакета под номером 0.9, одно из значительных улучшений в которой – это переделанная внутренняя структура, которая упрощает сторонним пользователям процесс подключения к разработке ggplot2. В начале 2014 года в git-репозитории ggplot2 насчитывается 21 разработчик. Хэдли – молодец, он смог грамотно перевести ggplot2 из формата личного проекта на рельсы разработки открытым сообществом. А это, я подозреваю, не просто! Нужно было наступить на горло своей песне, сократить время, затрачиваемое на чистую разработку и добавление функций, а сконцентрироваться на такой “скучной ерунде”, как написание документации, комментирование кода и т.д.

25 февраля 2014 года в почтовой группе, посвящённой ggplot2, Хэдли Викхэм официально объявил, что отныне работа над ggplot2 будет вестись в режиме поддержки: сами разработчики больше не будут добавлять в расширение новые функции, будут только вылавливать ошибки. Добавлением новых функций будут заниматься сами пользователи, которые могут дорабатывать ggplot2 и предлагать включать свои наработки в официальную версию. Соответственно, чтобы подчеркнуть переход ggplot2 на новый этап, версия расширения сменится на 1.0.0, а Хэдли ещё напишет популярным языком инструкцию о том, как сообщать об ошибках и предлагать свои наработки (чувствуете, опять работа с сообществом!).

Основы графической грамматики

В ggplot2 график является результатом взаимодействия ряда элементов:

  • Массив данных (data)
  • Схема соответствия переменных из массива визуальным средствам (aesthetic)
  • Геометрический объект (geom)
  • Статистическое преобразование (stat)
  • Координатная система (coord)
  • Ориентиры (guide)
  • Панели (facet)
  • Художественное оформление (theme)

Например, пользователь сообщает компьютеру, что хочет использовать массив данных про автомобили, переменная “скорость” будет выражена через положение по горизонтали, переменная “тормозной путь” через положение по вертикали, всё это нужно нарисовать с помощью геометрических объектов типа “точка”.

Вы спросите, а почему не заданы статистическое преобразование, координатная система, ориентиры и прочее. Все явно не указанные пользователем элементы графика берутся из значений по умолчанию. Например, если пользователь в качестве графического объекта указал тип “точки”, то по умолчанию статистических преобразований производиться не будет. А если он укажет тип “столбик”, то наблюдения в исходной переменной будут сгруппированы, а результатом применения статистики станет количество наблюдений в каждой группе.

На языке ggplot2 (да-да, я не оговорился: язык ggplot2 – это отдельный язык, а не R, но об этом позже) пример с точками будет выглядеть так:

ggplot(data = cars, aes(x = speed, y = dist)) + geom_point()

ggplot(data = cars, aes(x = speed, y = dist)) + geom_point()
Аналогичный результат мы получим, если укажем не геометрию (geom_point), а соответствующую ей по умолчанию статистику (stat_identity).

ggplot(data = cars, aes(x = speed, y = dist)) + stat_identity()

Пример со столбиками:

ggplot(data = cars, aes(x = speed)) + geom_bar()

ggplot(data = cars, aes(x = speed)) + geom_bar()
Умолчания распространяются и на другие элементы графика. Координатная система по умолчанию – декартова (coord_cartesian). Можем заменить на полярные координаты (coord_polar):

ggplot(data = cars, aes(x = speed)) + geom_bar() + coord_polar()

ggplot(data = cars, aes(x = speed)) + geom_bar() + coord_polar()
А что за дурацкий серый фон у графиков – часто звучит такой вопрос. Это тоже умолчание, но теперь уже в оформлении (theme_grey). Предполагается, что пользователь поместит график в публикацию, т.е. окружит его текстом. В обычной публикации цвет фона (бумаги) – белый, чернила – чёрные. “Средний” цвет страницы – серый. Если график сделать на белом фоне, то он будет резко выделяться на странице и вносить дисгармонию. А серый фон предназначен для уравновешивания цветового баланса страницы.

# Чёрно-белое оформление
ggplot(data = cars, aes(x = speed, y = dist)) + stat_identity() + theme_bw()

ggplot(data = cars, aes(x = speed, y = dist)) + stat_identity() + theme_bw()


Как нарисовать картограмму с помощью ggplot2 и ggmap


Продолжение в

Выборка частей массивов data.frame с помощью [квадратных скобок] в R

Функция [ (квадратная скобка) в R служит для отбора отдельных элементов или наборов элементов из векторов, списков (list), матриц и массивов (data.frame). Разберём, как это делать в data.frame.

Возьмём для примера массив “amis” из пакета “boot”.

library(boot)
data(amis)

В массиве содержатся результаты эксперимента по оценке влияния предупредительных дорожных знаков на скорость автомобилей. Массив состоит из четырёх колонок. speed – скорость автомобиля, period – время замера скорости: до установки знака / сразу после установки / через время после установки, warning – тип участка, на котором проводился замер: тестовый или контрольный, pair – порядковый номер пары участков, на которой проведён замер. В массиве 8437 строк, т.е. замеров скоростей.

head(amis)
##   speed period warning pair
## 1    26      1       1    1
## 2    26      1       1    1
## 3    26      1       1    1
## 4    26      1       1    1
## 5    27      1       1    1
## 6    28      1       1    1

Отбор по условию

Отберём те строки, которые относятся к контрольной группе (warning==2).

x <- amis[amis$warning == 2, ]
head(x)
##     speed period warning pair
## 201    21      1       2    1
## 202    25      1       2    1
## 203    25      1       2    1
## 204    26      1       2    1
## 205    26      1       2    1
## 206    26      1       2    1

В данном случае мы применили функцию квадратной скобки к массиву amis. Выражение, указанное между открывающей квадратной скобкой и запятой, задаёт требования для отбора строк: такие строки, в которых переменная warning равна двум. Между запятой и закрывающей квадратной скобкой указывается, какие колонки мы хотим отобрать. Пустота в нашем случае - это указание на отбор всех колонок.

Добавим к предыдущему примеру отбор только колонки со скоростью.

x <- amis[amis$warning == 2, "speed"]
head(x)
## [1] 21 25 25 26 26 26

Наш двухмерный массив превратился в одномерный вектор. В фукнцию “[” заложена такая логика, что если после отбора измерение становится лишним, то оно удаляется. Иногда это полезно, иногда нет. Чтобы оставить два измерения, необходимо добавить аргумент drop=FALSE.

x <- amis[amis$warning == 2, "speed", drop = FALSE]
head(x)
##     speed
## 201    21
## 202    25
## 203    25
## 204    26
## 205    26
## 206    26

Отбор по номерам строк и колонок

Снова отберём колонку со скоростью, но указываем колонку не по имени, а по её порядковому номеру (она в массиве первая).

x <- amis[amis$warning == 2, 1, drop = FALSE]
head(x)
##     speed
## 201    21
## 202    25
## 203    25
## 204    26
## 205    26
## 206    26

Номера можно указывать и при отборе строк. Отберём строки 7, 9, 16 и 18 из колонок №1, 2.

x <- amis[c(7, 9, 16, 18), 1:2]
head(x)
##    speed period
## 7     28      1
## 9     28      1
## 16    29      1
## 18    30      1

Объединение условий для отбора

Критерии для отбора можно комбинировать через операторы. Отберём строки из тестовой группы (warnning=1) и из временного периода до установки знака (period=1). Объединение происходит через логический оператор & (логическое И).

x <- amis[amis$warning == 1 & amis$period == 1, ]
head(x)
##   speed period warning pair
## 1    26      1       1    1
## 2    26      1       1    1
## 3    26      1       1    1
## 4    26      1       1    1
## 5    27      1       1    1
## 6    28      1       1    1

Отберём наблюдения из контрольной группы (warning=2), относящиеся к временным периодам сразу после установки знака (period=2) и через время после установки знака (period-3). Возможно несколько вариантов.

x <- amis[amis$warning == 2 & (amis$period == 2 | amis$period == 3), ]
x1 <- amis[amis$warning == 2 & amis$period > 1, ]
x2 <- amis[amis$warning == 2 & amis$period %in% 2:3, ]

В x мы использовали оператор “логическое ИЛИ” и скобки, чтобы установить нужный нам порядок вычисления условий.

В x2 мы использовали функцию %in%, которая проверяет, равен ли первый аргумент (amis$period) какому-либо из значений второго арумента (2 или 3).

Page 1 of 2

Powered by WordPress & Theme by Anders Norén