Формат данных для временных рядов, который избавит вас от головной боли

Анализ временных рядов в R часто сопровождается приличной "головной болью", что связано с необходимостью представления данных в виде объектов таких классов, как ts, zoo или xts. Все эти форматы противоречат принципам организации и хранения т.н. "опрятных данных" ("tidy data"; Wickham 2014) и, как следствие, затрудняют анализ и моделирование с помощью широко используемых сегодня инструментов из группы tidyverse. Для решения этой проблемы группа исследователей под руководством проф. Роба Хиндмана (Rob Hyndman) разработала новый формат для представления и хранения временных рядов, реализованный в пакете tsibble. Ниже описаны основные особенности работы с этим пакетом.

Инсталляция пакета tsibble

Пакет tsibble можно установить обычным образом из хранилища CRAN:

install.packages("tsibble")

Данные, используемые в примерах

В приведенных ниже примерах использованы данные по стоимости 22 различных криптовалют на момент закрытия торгов в период с 1 января 2018 г. по 6 декабря 2019 г. Эти данные были собраны с сайта CoinMarketCap. Скрипт для сбора данных, а также весь код для воспроизведения примеров, можно найти в репозитории ranalytics/tsibble_format. Данные имеют простую структуру:

require(readr)

(dat <- read_csv("../data/cypto_coins.csv"))

# A tibble: 15,510 x 3
       y ds         coin   
   <dbl> <date>     <chr>  
 1 7547  2019-12-06 bitcoin
 2 7448. 2019-12-05 bitcoin
 3 7252. 2019-12-04 bitcoin
 4 7320. 2019-12-03 bitcoin
 5 7322. 2019-12-02 bitcoin
 6 7424. 2019-12-01 bitcoin
 7 7570. 2019-11-30 bitcoin
 8 7761. 2019-11-29 bitcoin
 9 7463. 2019-11-28 bitcoin
10 7532. 2019-11-27 bitcoin
# ... with 15,500 more rows

Столбец y - это стоимость ($) криптовалюты coin, отмеченная в день ds.

Графически эти данные выглядят так:

require(ggplot2)

dat %>% 
  ggplot(., aes(ds, y, group = coin)) +
  geom_line() +
  scale_y_log10() +
  theme_bw()

Рис. 1

Формат данных tsibble

Как отмечено выше, пакет tsibble предназначен для создания объектов с данными временных рядов в соответствии с принципами "опрятных данных", а именно:

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

Для создания объектов класса tsibble служит функция as_tsibble(), которая имеет следующие аргументы:

• x - объект с данными, подлежащий преобразованию в объект класса tsibble (это может быть, например, числовой вектор, матрица, таблица с данными (data.frame или tibble) и др.).
• index - переменная с временными отметками (указывается без кавычек). Допускается использование временных отметок на шкале от наносекунд до года.
• key - одна или несколько группирующих переменных, уникально определяющих каждый хранящийся в таблице временной ряд. Названия переменных указываются без кавычек и объединяются с помощью функции конкатенации c(). По умолчанию этот аргумент равен NULL, т.е. предполагается, что группирующих переменных в таблице нет.
• regular - логический аргумент, указывает на регулярность учета хранящихся в таблице наблюдений. Значение TRUE (принято по умолчанию) предполагает, что учет выполнялся с одинаковым интервалом (например, каждую минуту, час, день, и т.п.).
• validate - логический аргумент (по умолчанию равен TRUE), позволяющий выполнить проверку уникальности каждого наблюдения по сочетанию значений переменных index и key. Если вы уверены, что каждое наблюдение уникально, то эту проверку можно отключить (FALSE), что приведен к более быстрому выполнению команды as_tsibble() в случае с большими наборами данных.
• .drop - логический аргумент, позволяющий исключить из таблицы "пустые" временные ряды, т.е. такие сочетания значений группирующих переменных, для которых значения x отсутствуют.

Применим функцию as_tibble() к данным по стоимости криптовалют:

require(tsibble)

(dat_ts <- as_tsibble(dat, key = coin, index = ds))

# A tsibble: 15,510 x 3 [1D]
# Key:       coin [22]
       y ds         coin 
   <dbl> <date>     <chr>
 1  74.5 2018-01-01 augur
 2  79.5 2018-01-02 augur
 3  77.5 2018-01-03 augur
 4  73.7 2018-01-04 augur
 5  72.8 2018-01-05 augur
 6  77.5 2018-01-06 augur
 7  80.2 2018-01-07 augur
 8  98.1 2018-01-08 augur
 9  91.9 2018-01-09 augur
10 103.  2018-01-10 augur
# ... with 15,500 more rows

Если внимательно присмотреться, можно увидеть некоторые внешние отличия полученной таблицы dat_ts от исходной dat:

• R теперь "знает", что таблица dat_ts является объектом класса tsibble с 15510 наблюдениями, учтенными с дневным интервалом (см. комментарий "# A tsibble: 15,510 x 3 [1D]");
• у таблицы dat_ts есть группирующая переменная coin с 22 уровнями (см. "# Key: coin [22]");
• все наблюдения каждого временного ряда упорядочены по возрастанию значений индексирующей переменной ds.

Однако это всего лишь внешние различия. Гораздо важнее то, что к данным из таблицы dat_ts теперь можно применять методы анализа, специфичные для временных рядов, и делать это обычным для инструментов tidyverse образом. Рассмотрим некоторые примеры.

Агрегирование по календарному периоду

Обычной задачей при работе с временными рядам является расчет агрегированных показателей за определенный календарный период (например, средние значения за неделю, месяц, квартал и т.п.). В пакете tsibble для таких вычислений используется связка из двух функций: index_by() + summarise(). Первая из этих функций входит в состав пакета tsibble и аналогична group_by() из пакета dplyr. Как следует из ее названия, index_by() группирует данные по заданному пользователем календарному периоду. Для указания необходимого периода используются такие функции из пакета tsibble, как yearweek() (неделя),  yearmonth() (месяц) и yearquarter() (квартал). Кроме того, можно также использовать базовую as.Date() и многие функции из пакета lubridate. Вторая функция из указанной выше связки - summarise() - входит в состав пакета dplyr и служит для вычисления соответствующих агрегированных величин. В качестве примера рассчитаем среднемесячную стоимость каждой криптовалюты, а также общее число наблюдений, учтенных в каждом месяце:

require(dplyr)

monthly <- dat_ts %>% 
  group_by_key() %>% 
  index_by(year_month = ~ yearmonth(.)) %>% 
  summarise(
    avg_y = mean(y),
    n = n()
  )

monthly

# A tsibble: 528 x 4 [1M]
# Key:       coin [22]
   coin  year_month avg_y     n
   <chr>      <mth> <dbl> <int>
 1 augur   2018 Jan  84.1    31
 2 augur   2018 Feb  51.1    28
 3 augur   2018 Mar  35.9    31
 4 augur   2018 Apr  32.5    30
 5 augur   2018 May  45.8    31
 6 augur   2018 Jun  34.5    30
 7 augur   2018 Jul  31.9    31
 8 augur   2018 Aug  21.7    31
 9 augur   2018 Sep  14.7    30
10 augur   2018 Oct  13.1    31
# ... with 518 more rows

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

На приведенном ниже рисунке показан результат вычисления среднемесячной стоимости анализируемых крипотовалют:

monthly %>% 
  ggplot(., aes(year_month, avg_y, group = coin)) + 
  geom_line() +
  scale_y_log10() +
  theme_bw()

Рис. 2

Агрегирование по скользящему окну

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

• slide_ - выполняют вычисления в пределах перекрывающихся окон размером .size (размер "нахлеста" можно контролировать с помощью аргумента .step);
• tile_ - вычисления ведутся в пределах следующих друг за другом, но неперекрывающихся окон размером .size;
• stretch_ - исходная ширина окна .init постепенно увеличивается на .step временных отметок.

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

Приведенный ниже анимированный рисунок поможет понять суть этих вычислений (заимствован из официальной документации по пакету tsibble):

Рис. 3

В качестве примера вычислим скользящую среднюю стоимость каждой криптовалюты в пределах окна шириной в 7 дней:

dat_ts %>% 
  group_by_key() %>% 
  mutate(moving_avg_y = slide_dbl(y, ~mean(.), .size = 7))

# A tsibble: 15,510 x 4 [1D]
# Key:       coin [22]
# Groups:    coin [22]
       y ds         coin  moving_avg_y
   <dbl> <date>     <chr>        <dbl>
 1  74.5 2018-01-01 augur         NA  
 2  79.5 2018-01-02 augur         NA  
 3  77.5 2018-01-03 augur         NA  
 4  73.7 2018-01-04 augur         NA  
 5  72.8 2018-01-05 augur         NA  
 6  77.5 2018-01-06 augur         NA  
 7  80.2 2018-01-07 augur         76.5
 8  98.1 2018-01-08 augur         79.9
 9  91.9 2018-01-09 augur         81.7
10 103.  2018-01-10 augur         85.3
# ... with 15,500 more rows

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

Обработка пропущенных наблюдений

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

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

has_gaps(dat_ts)

# A tibble: 22 x 2
   coin      .gaps
   <chr>     <lgl>
 1 augur     FALSE
 2 bitcoin   FALSE
 3 cardano   FALSE
 4 chainlink FALSE
 5 dash      FALSE
 6 decred    FALSE
 7 dogecoin  FALSE
 8 eos       FALSE
 9 ethereum  FALSE
10 iota      FALSE
11 litecoin  FALSE
12 maker     FALSE
13 monero    FALSE
14 nano      FALSE
15 neo       FALSE
16 qtum      FALSE
17 stellar   FALSE
18 tether    FALSE
19 tezos     FALSE
20 tron      FALSE
21 xrp       FALSE
22 zcash     FALSE

Как видим, в нашем наборе данных пропущенных наблюдений нет - все значения столбца .gaps в полученной таблице с результатами проверки равны FALSE. Посмотрим, что получится, если мы случайным образом удалим несколько строк из таблицы dat_ts:

set.seed(42)

dat_ts_na <- dat %>% 
  sample_n(., 15000) %>% 
  as_tsibble(., key = coin, index = ds)

has_gaps(dat_ts_na)

# A tibble: 22 x 2
   coin      .gaps
   <chr>     <lgl>
 1 augur     TRUE 
 2 bitcoin   TRUE 
 3 cardano   TRUE 
 4 chainlink TRUE 
 5 dash      TRUE 
 6 decred    TRUE 
 7 dogecoin  TRUE 
 8 eos       TRUE 
 9 ethereum  TRUE 
10 iota      TRUE 
11 litecoin  TRUE 
12 maker     TRUE 
13 monero    TRUE 
14 nano      TRUE 
15 neo       TRUE 
16 qtum      TRUE 
17 stellar   TRUE 
18 tether    TRUE 
19 tezos     TRUE 
20 tron      TRUE 
21 xrp       TRUE 
22 zcash     TRUE

Теперь пропущенные наблюдения есть в каждом временном ряду. С помощью функции scan_gaps() можно выяснить, какие именно наблюдения отсутствуют в каком из рядов:

scan_gaps(dat_ts_na)

# A tsibble: 509 x 2 [1D]
# Key:       coin [22]
   coin  ds        
   <chr> <date>    
 1 augur 2018-01-06
 2 augur 2018-02-19
 3 augur 2018-03-29
 4 augur 2018-04-06
 5 augur 2018-04-23
 6 augur 2018-05-27
 7 augur 2018-06-03
 8 augur 2018-06-12
 9 augur 2018-06-27
10 augur 2018-09-15
# ... with 499 more rows

Функция count_gaps() дает более развернутый отчет - в частности, она определяет длину каждого "пробела" в каждом временном ряду:

(gaps <- count_gaps(dat_ts_na))

# A tibble: 496 x 4
   coin  .from      .to           .n
   <chr> <date>     <date>     <int>
 1 augur 2018-01-06 2018-01-06     1
 2 augur 2018-02-19 2018-02-19     1
 3 augur 2018-03-29 2018-03-29     1
 4 augur 2018-04-06 2018-04-06     1
 5 augur 2018-04-23 2018-04-23     1
 6 augur 2018-05-27 2018-05-27     1
 7 augur 2018-06-03 2018-06-03     1
 8 augur 2018-06-12 2018-06-12     1
 9 augur 2018-06-27 2018-06-27     1
10 augur 2018-09-15 2018-09-15     1
# ... with 486 more rows

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

gaps %>% 
  ggplot(., aes(x = coin)) +
  geom_linerange(aes(ymin = .from, ymax = .to)) +
  geom_point(aes(y = .from)) +
  geom_point(aes(y = .to)) +
  coord_flip() +
  theme_bw()

Рис. 4

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

dat_ts_na_filled <- dat_ts_na %>% 
  fill_gaps()

dat_ts_na_filled

# A tsibble: 15,510 x 3 [1D]
# Key:       coin [22]
       y ds         coin 
   <dbl> <date>     <chr>
 1  74.5 2018-01-01 augur
 2  79.5 2018-01-02 augur
 3  77.5 2018-01-03 augur
 4  73.7 2018-01-04 augur
 5  72.8 2018-01-05 augur
 6  NA   2018-01-06 augur
 7  80.2 2018-01-07 augur
 8  98.1 2018-01-08 augur
 9  91.9 2018-01-09 augur
10 103.  2018-01-10 augur
# ... with 15,500 more rows

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

dat_ts_na_filled_median <- dat_ts_na %>% 
  group_by_key() %>% 
  fill_gaps(y = median(y, na.rm = TRUE))

dat_ts_na_filled_median

# A tsibble: 15,510 x 3 [1D]
# Key:       coin [22]
# Groups:    coin [22]
       y ds         coin 
   <dbl> <date>     <chr>
 1  74.5 2018-01-01 augur
 2  79.5 2018-01-02 augur
 3  77.5 2018-01-03 augur
 4  73.7 2018-01-04 augur
 5  72.8 2018-01-05 augur
 6  15.3 2018-01-06 augur
 7  80.2 2018-01-07 augur
 8  98.1 2018-01-08 augur
 9  91.9 2018-01-09 augur
10 103.  2018-01-10 augur
# ... with 15,500 more rows

Заключение

В этом сообщении были рассмотрены основные особенности работы с пакетом tsibble, в котором реализован новый формат данных для временных рядов, а также ряд функций по подготовке таких данных для последующего анализа с использованием мощного арсенала инструментов tidyverse. Если вы еще используете "классические" для R форматы ts, zoo или xts, я бы порекомендовал перейти на tsibble.