class: center, middle, inverse, title-slide .title[ # Временные ряды ] .subtitle[ ## Визуализация и анализ географических данных на языке R ] .author[ ### Тимофей Самсонов ] .institute[ ### МГУ имени Ломоносова, Географический факультет ] .date[ ### 2023-10-24 ] --- ## Используемые пакеты: .pull-left[ __tidyverse__ ```r # tidyverse library(tidyverse) library(readxl) library(googlesheets4) ``` __Данные__ ```r library(tsibbledata) library(nasapower) ``` ] .pull-right[ __Временные данные__ ```r # пакеты для временных данных library(gganimate) library(zoo) library(trend) library(tsibble) library(feasts) library(fable) ``` ] --- ## Основные аспекты работы с временными рядами 1. Типы данных времени 2. Структуры данных для временных рядов 3. Преобразования времени и временных рядов 4. Декомпозиция временных рядов 5. Анализ трендов и автокорреляции 6. Построение моделей и прогнозирование 7. Визуализация и анимация --- ## Создание и преобразование дат и времени Текущее время и дата: ```r (date = Sys.Date()) ## [1] "2023-10-24" (time = Sys.time()) ## [1] "2023-10-24 11:31:54 MSK" ``` Полученные объекты имеют типы `Date` и `POSIXct`: ```r class(date) ## [1] "Date" class(time) ## [1] "POSIXct" "POSIXt" ``` --- ## Система счета времени Время `POSIXct` хранится в целочисленном формате и отсчитывается начиная с начала эпохи UNIX — `\(1\)` января `\(1970\)` года по гринвичскому (_UTC_) времени: ```r as.integer(date) ## [1] 19654 difftime(date, as.Date('1970-01-01')) ## Time difference of 19654 days as.integer(time) ## [1] 1698136314 difftime(time, as.POSIXct('1970-01-01 00:00:00', tz = 'UTC'), units = 'secs') ## Time difference of 1698136315 secs ``` --- ## Работа с данными в lubridate __Создание__ дат возможно на основе целочисленных и строковых значений: ```r ymd(20150515) ## [1] "2015-05-15" dmy('15052015') ## [1] "2015-05-15" ``` Для создания отметки времени необходимо сформировать строку, которая может быть интерпретирована должным образом. Часовой пояс указывается в `tz`: ```r ymd_hms('2015-05-15 22:15:34') # по умолчанию Гринвичское время ## [1] "2015-05-15 22:15:34 UTC" ymd_hms('2015-05-15 22:15:34', tz = "Europe/Moscow") ## [1] "2015-05-15 22:15:34 MSK" ``` --- ## Работа с данными в lubridate .pull-left[ __Извлечение__ компоненты даты/времени: ```r year(time) ## [1] 2023 month(time) ## [1] 10 week(time) ## [1] 43 day(time) ## [1] 24 hour(time) ## [1] 11 second(time) ## [1] 54.60705 ``` > Недели отсчитываются от начала года, а не месяца. ] .pull-right[ __Замена__ компонент даты/времени: ```r year(time) = 2015 month(time) = 01 time ## [1] "2015-01-24 11:31:54 MSK" ``` __Округление__ даты/времени. Например, получить первый день в текущем году можно так: ```r floor_date(Sys.Date(), unit = 'year') ## [1] "2023-01-01" ``` ] --- ## Работа с данными в lubridate __Периоды__ (_periods_) — это промежутки дат с точностью до _дня_. С помощью периода можно, например, к дате прибавить 1 год, 4 месяца, 3 недели и 2 дня: ```r date ## [1] "2023-10-24" date + years(1) + months(4) + weeks(3) + days(2) ## [1] "2025-03-19" ``` __Длительности__ (_durations_) — это промежутки времени с точностью до _секунды_. Работают они в целом аналогично периодам: ```r dweeks(1) ## [1] "604800s (~1 weeks)" time ## [1] "2015-01-24 11:31:54 MSK" time + dweeks(1) ## [1] "2015-01-31 11:31:54 MSK" time + weeks(1) ## [1] "2015-01-31 11:31:54 MSK" ``` --- ## Работа с данными в lubridate __Интервалы__ — это отрезки между двумя датами. Интервал можно преобразовывать в периоды и длительности: ```r int = lubridate::interval(Sys.time(), time) int ## [1] 2023-10-24 11:31:54 MSK--2015-01-24 11:31:54 MSK as.period(int, 'days') ## [1] "-3195d 0H 0M -0.139196157455444S" as.duration(int) ## [1] "-276048000.139196s (~-8.75 years)" ``` --- ## Восстановление пропусков В качестве источника данных будем использовать данные об уровне воды на гидропосте Паялка (р. Умба, Мурманская область) с 1932 по 2014 год. Отсутствие информации в файле данных закодировано числом `9999`: ```r src = read_delim('../r-geo-course/data/in_Umba.txt', delim = ' ', col_names = c('day', 'month', 'year', 'level'), na = '9999') head(src) ## # A tibble: 6 × 4 ## day month year level ## <dbl> <dbl> <dbl> <dbl> ## 1 1 1 1932 49.4 ## 2 2 1 1932 NA ## 3 3 1 1932 NA ## 4 4 1 1932 NA ## 5 5 1 1932 NA ## 6 6 1 1932 NA ``` --- ## Восстановление пропусков Сформируем даты на основе первых трёх столбцов и проверим, все ли из них корректны. Если компоненты даты некорректны, то функция `yms()` вернет `NA`: .left-60[ ```r tab = src |> mutate(Date = ymd(paste(year, month, day))) tab |> filter(is.na(Date)) ## # A tibble: 1 × 5 ## day month year level Date ## <dbl> <dbl> <dbl> <dbl> <date> ## 1 29 2 1941 27.4 NA ``` ] .right-40[ <br> Проверка показала, что преобразование в дату оказалось невозможно только для одной строки. В этой строке оператором была введена несуществующая дата — 29 февраля 1941 года. Для дальнейшей работы необходимо избавиться от некорректных дат. ] --- ## Создание временного фрейма данных (tsibble) __tsibble__ — это временной фрейм данных, у которого есть индекс в виде временного столбца: ```r tsbl = tab |> filter(!is.na(Date)) |> distinct(Date, .keep_all = T) |> as_tsibble(index = Date) tsbl ## # A tsibble: 30,314 x 5 [1D] ## day month year level Date ## <dbl> <dbl> <dbl> <dbl> <date> ## 1 1 1 1932 49.4 1932-01-01 ## 2 2 1 1932 NA 1932-01-02 ## 3 3 1 1932 NA 1932-01-03 ## 4 4 1 1932 NA 1932-01-04 ## 5 5 1 1932 NA 1932-01-05 ## 6 6 1 1932 NA 1932-01-06 ## 7 7 1 1932 NA 1932-01-07 ## 8 8 1 1932 NA 1932-01-08 ## 9 9 1 1932 NA 1932-01-09 ## 10 10 1 1932 47 1932-01-10 ## # ℹ 30,304 more rows ``` --- ## Восстановление пропусков Для того чтобы выяснить, есть ли пропуски, можно использовать `scan_gaps()` (индивидуальные пропуски) или `count_gaps()` (периоды пропусков) ```r scan_gaps(tsbl) ## # A tsibble: 2 x 1 [1D] ## Date ## <date> ## 1 1941-12-31 ## 2 1942-01-01 count_gaps(tsbl) ## # A tibble: 1 × 3 ## .from .to .n ## <date> <date> <int> ## 1 1941-12-31 1942-01-01 2 ``` --- ## Восстановление пропусков Для устранения пропусков дат используем `fill_gaps()`: ```r tsbl = fill_gaps(tsbl) scan_gaps(tsbl) ## # A tsibble: 0 x 1 [?] ## # ℹ 1 variable: Date <date> ``` --- ## Заполнение пропусков данных Создание таблицы пропусков: ```r find_gaps = function(df, variable) { df |> as_tibble() |> mutate(type = if_else(is.na(df[variable]), 'gap', 'data'), num = with(rle(type), rep(seq_along(lengths), lengths))) |> group_by(num) |> summarise(start_date = min(Date), end_date = max(Date), duration = end_date - start_date + 1, type = first(type)) } timerep = find_gaps(tab, 'level') ``` --- ## Заполнение пропусков
--- ## Автокорреляция Критическую длину пропуска, заполняемого интерполяцией, можно связать с пороговым значением __автокорреляции__ — коэффициента корреляции исходного ряда данных и его копии, полученной со сдвигом `\(\tau\)`. .pull-left[ ```r acorr = tsbl |> ACF(level) acorr ## # A tsibble: 41 x 2 [1D] ## lag acf ## <cf_lag> <dbl> ## 1 1D 0.994 ## 2 2D 0.982 ## 3 3D 0.964 ## 4 4D 0.943 ## 5 5D 0.919 ## 6 6D 0.893 ## 7 7D 0.866 ## 8 8D 0.839 ## 9 9D 0.812 ## 10 10D 0.786 ## # ℹ 31 more rows ``` ] .pull-right[ ```r autoplot(acorr) ``` ![](08_TemporalStats_files/figure-html/unnamed-chunk-22-1.png)<!-- --> ] --- ## Анализ АКФ За допустимую величину автокорреляции при восстановлении рядов данных принято брать некоторое значение. Например, гидрологи часто используют `\(0.7\)`. Найдем индекс первого элемента менее данной величины, используя функцию `detect_index()` из пакета __purrr__: ```r max_dur = purrr::detect_index(acorr$acf, ~ .x < 0.7) max_dur ## [1] 14 ``` Полученное значение говорит нам о том, что при заданном допуске допустимо интерполировать значения в пропусках данных короче, чем 14 дней. --- ## Интерполяция по времени Для начала визуализируем участок в окрестности одного из пробелов ```r autoplot(tsbl, level, size = 1) + scale_x_date(limits = c(ymd(19320101), ymd(19360101))) ``` ![](08_TemporalStats_files/figure-html/unnamed-chunk-24-1.png)<!-- --> --- ## Интерполяция по времени (линейная) Для выполнения обычной линейной интерполяции между пропускаи воспользуемся функцией `na.approx()` из пакета __zoo__ и округлим полученные значения до одного знака после запятой (что соответствует точности исходных данных): ```r (tsbl_int = tsbl |> mutate(level_interp = zoo::na.approx(level, maxgap = max_dur) |> round(1))) ## # A tsibble: 30,316 x 6 [1D] ## day month year level Date level_interp ## <dbl> <dbl> <dbl> <dbl> <date> <dbl> ## 1 1 1 1932 49.4 1932-01-01 49.4 ## 2 2 1 1932 NA 1932-01-02 49.1 ## 3 3 1 1932 NA 1932-01-03 48.9 ## 4 4 1 1932 NA 1932-01-04 48.6 ## 5 5 1 1932 NA 1932-01-05 48.3 ## 6 6 1 1932 NA 1932-01-06 48.1 ## 7 7 1 1932 NA 1932-01-07 47.8 ## 8 8 1 1932 NA 1932-01-08 47.5 ## 9 9 1 1932 NA 1932-01-09 47.3 ## 10 10 1 1932 47 1932-01-10 47 ## # ℹ 30,306 more rows ``` --- ## Интерполяция по времени (линейная) Визуализируем результат: ```r autoplot(tsbl_int, level_interp, color = 'red', size = 1) + geom_line(aes(Date, level), tsbl, size = 1) + scale_x_date(limits = c(ymd(19320101), ymd(19360101))) ``` ![](08_TemporalStats_files/figure-html/unnamed-chunk-26-1.png)<!-- --> --- ## Интерполяция по времени (с авторегрессией) Ту же задачу можно решить с применением модели [ARIMA](https://ru.wikipedia.org/wiki/ARIMA), которая подыскивает интерполятор с учетом авторегрессии и сглаживающего среднего. ```r tsbl_int2 = tsbl |> model(mdl = ARIMA(level ~ 0 + pdq(0,1,0) + PDQ(0,0,0))) |> interpolate(tsbl) ``` --- ## Интерполяция по времени (с авторегрессией) Однако в простейшем случае здесь есть проблемы с немонотонными участками (начало 1935 года). ```r autoplot(tsbl_int2, level, color = 'magenta', size = 1) + geom_line(aes(Date, level), tsbl, size = 1) + scale_x_date(limits = c(ymd(19320101), ymd(19360101))) ``` ![](08_TemporalStats_files/figure-html/unnamed-chunk-28-1.png)<!-- --> --- ## Интерполяция с учетом предикторов Используем данные по температуре и осадкам и посчитаем по ним скользящее среднее за 10 и 5 дней соответственно: ```r meteo = read_delim('../r-geo-course/data/in_Umba_meteo.txt', delim = ' ', col_names = c('day', 'month', 'year', 'temp', 'prec'), na = '9999') |> mutate(Date = ymd(paste(year, month, day)), tempsm = zoo::rollmean(temp, 10, fill = NA, align = 'center'), precsum = zoo::rollsum(prec, 5, fill = NA, align = 'center')) |> select(Date, temp, tempsm, prec, precsum) tsbl_all = left_join(tsbl_int, meteo, by = 'Date') tsbl_all ## # A tsibble: 30,316 x 10 [1D] ## day month year level Date level_interp temp tempsm prec precsum ## <dbl> <dbl> <dbl> <dbl> <date> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 1 1 1932 49.4 1932-01-01 49.4 -14.9 NA 0.184 NA ## 2 2 1 1932 NA 1932-01-02 49.1 -15.1 NA 0.578 NA ## 3 3 1 1932 NA 1932-01-03 48.9 -10.4 NA 0.942 3.97 ## 4 4 1 1932 NA 1932-01-04 48.6 -7.58 NA 1.44 4.07 ## 5 5 1 1932 NA 1932-01-05 48.3 -8.46 -11.3 0.826 4.42 ## 6 6 1 1932 NA 1932-01-06 48.1 -13.0 -10.7 0.288 4.64 ## 7 7 1 1932 NA 1932-01-07 47.8 -13.3 -9.72 0.925 3.95 ## 8 8 1 1932 NA 1932-01-08 47.5 -10.4 -8.83 1.16 3.79 ## 9 9 1 1932 NA 1932-01-09 47.3 -10.4 -8.33 0.75 4.53 ## 10 10 1 1932 47 1932-01-10 47 -9.16 -7.62 0.662 4.49 ## # ℹ 30,306 more rows ``` --- ## Интерполяция с учетом предикторов Выполним интерполяцию методом ARIMA с внешними предикторами: ```r tsbl_int_met = tsbl_all |> model(mdl = ARIMA(level_interp ~ tempsm + precsum)) |> interpolate(tsbl_all) ``` --- ## Интерполяция с учетом предикторов Визуализируем результат: ```r autoplot(tsbl_int_met, level_interp, color = 'green', size = 1) + geom_line(aes(Date, level_interp), tsbl_all, size = 1) + scale_x_date(limits = c(ymd(19320101), ymd(19360101))) ``` ![](08_TemporalStats_files/figure-html/unnamed-chunk-31-1.png)<!-- --> --- ## Пересчет на другую временную сетку Загрузим в качестве примера данные по метеостанции __NETATMO__: ```r tab = read_csv('../r-geo-course/data/70_ee_50_00_8e_1a.csv') ```
--- ## Пересчет на другую временную сетку Определим расчетные интервалы, кратные 30 минутам: ```r tab = tab |> select(id, temperature, time_temperature) tmin = ceiling_date(min(tab$time_temperature), unit = '30 minutes') tmax = floor_date(max(tab$time_temperature), unit = '30 minutes') ``` --- ## Пересчет на другую временную сетку Проинтерполуруем данные на единую регулярную временную сетку через `\(30\)` минут, используя функцию [`approx()`](https://stat.ethz.ch/R-manual/R-devel/library/stats/html/approxfun.html) из базового __R__: ```r time_interp = tibble( datetime = seq(tmin, tmax, by = '30 min'), temp = round(approx(tab$time_temperature, tab$temperature, xout = datetime)$y, 1) ) ``` --- ## Пересчет на другую временную сетку
--- ## Обеспеченность данными Надежность каждого интерполированного значения можно оценить, задав максимальную продолжительность интерполируемого интервала. Предположим, что она не должна превышать `\(60\)` минут: ```r idx = sapply(time_interp$datetime, \(x) match(TRUE, tab$time_temperature >= x)) time_interp = time_interp |> mutate(before_time = datetime - tab$time_temperature[idx-1], after_time = tab$time_temperature[idx] - datetime, valid = after_time - before_time <= minutes(60)) cat('Valid are ', round(100 * sum(time_interp$valid) / nrow(time_interp), 1), '% of interpolated values', sep = '') ## Valid are 86.3% of interpolated values ``` --- ## Обеспеченность данными
--- ## Прогнозирование по времени Рассмотрим задачу прогнозирования количества [потребляемой электроэнергии](https://otexts.com/fpp3/forecasting.html#example-forecasting-electricity-demand) в штате Виктория (Австралия) в зависимости от температуры воздуха. .code-small[ ```r data(vic_elec, package = 'tsibbledata') vic_elec_daily = vic_elec |> filter(year(Time) == 2014) |> index_by(Date = date(Time)) |> summarise( Demand = sum(Demand) / 1e3, Temperature = max(Temperature), Holiday = any(Holiday) ) |> mutate(Day_Type = case_when( Holiday ~ "Holiday", wday(Date) %in% 2:6 ~ "Weekday", TRUE ~ "Weekend" )) ``` ] --- ## Прогнозирование по времени Визуализируем исходные данные: ```r vic_elec_daily |> pivot_longer(c(Demand, Temperature)) |> ggplot(aes(x = Date, y = value)) + geom_line() + facet_grid(name ~ ., scales = "free_y") + ylab("") ``` ![](08_TemporalStats_files/figure-html/unnamed-chunk-40-1.png)<!-- --> --- ## Прогнозирование по времени Выполним построение модели и сформируем прогнозные данные по температуре: ```r fit = vic_elec_daily |> model(ARIMA(Demand ~ Temperature + I(Temperature^2) + (Day_Type == "Weekday"))) vic_elec_future = new_data(vic_elec_daily, 14) |> mutate( Temperature = 26, Holiday = c(TRUE, rep(FALSE, 13)), Day_Type = case_when( Holiday ~ "Holiday", wday(Date) %in% 2:6 ~ "Weekday", TRUE ~ "Weekend" ) ) ``` --- ## Прогнозирование по времени Выполним прогнозирование: ```r forecast(fit, vic_elec_future) |> autoplot(vic_elec_daily) + labs(title="Daily electricity demand: Victoria", y="GW") + scale_x_date(limits = c(ymd(20141001), ymd(20150115))) ``` ![](08_TemporalStats_files/figure-html/unnamed-chunk-42-1.png)<!-- --> --- ## Аппроксимация по времени В качестве примера возьмем данные межгодичных изменений характеристик стока реки Мезень на посту Малонисогорская: ```r (tab = read_csv('../r-geo-course/data/Mezen.csv')) ## # A tibble: 75 × 57 ## year_number Year1 Year2 datestart datepolend Qy Qmax datemax Qygr ## <dbl> <dbl> <dbl> <date> <date> <dbl> <dbl> <date> <dbl> ## 1 1 1938 1939 2000-03-23 2000-06-08 98.1 747 2000-04-03 19.9 ## 2 2 1939 1940 2000-03-16 2000-07-11 66.9 487 2000-04-28 14.8 ## 3 3 1940 1941 2000-03-09 2000-07-12 97.4 995 2000-04-28 18.9 ## 4 4 1941 1942 2000-03-24 2000-07-21 214. 2030 2000-05-01 31.4 ## 5 5 1942 1943 2000-04-02 2000-07-31 242. 2790 2000-05-08 26.9 ## 6 6 1943 1944 2000-03-24 2000-06-15 71.3 451 2000-05-02 26.5 ## 7 7 1944 1945 2000-02-26 2000-06-30 89.0 530 2000-04-28 27.0 ## 8 8 1945 1946 2000-03-25 2000-07-13 128. 1220 2000-04-27 33.9 ## 9 9 1946 1947 2000-03-15 2000-06-30 179. 2150 2000-04-28 25.8 ## 10 10 1947 1948 2000-03-06 2000-07-15 153. 1720 2000-04-16 31.3 ## # ℹ 65 more rows ## # ℹ 48 more variables: Qmmsummer <dbl>, monmmsummer <date>, Qmmwin <dbl>, ## # nommwin <date>, Q30s <dbl>, date30s1 <date>, date30s2 <date>, Q30w <dbl>, ## # date30w1 <date>, date30w2 <date>, Q10s <dbl>, date10s1 <date>, ## # date10s2 <date>, Q10w <dbl>, date10w1 <date>, date10w2 <date>, Q5s <dbl>, ## # date5s1 <date>, date5s2 <date>, Q5w <dbl>, date5w1 <date>, date5w2 <date>, ## # Wy <dbl>, Wgr <dbl>, Wpol1 <dbl>, Wpol2 <dbl>, Wpol3 <dbl>, Wpavs1 <dbl>, … ``` --- ## Аппроксимация по времени Построим график межгодичной изменчивости объема грунтового стока (переменная `Wgr` в `\(км^3\)`). Для того чтобы нанести линию глобального тренда, добавим грамматику `geom_smooth(method = 'lm')`: .pull-left[ .code-small[ ```r ggplot( tab, mapping = aes(Year1, Wgr) ) + geom_line() + geom_area(alpha = 0.5) + geom_smooth(method = 'lm') + labs( title = 'Объем грунтового стока', x = 'Год', y = expression(км^3) ) ``` ] ] .pull-right[ ![](08_TemporalStats_files/figure-html/wgr-out-1.png)<!-- --> ] --- ## Аппроксимация по времени Добавим линию локальной регрессии `geom_smooth(method = 'loess')`: .pull-left[ .code-small[ ```r ggplot( tab, mapping = aes(Year1, Wgr) ) + geom_line() + geom_area(alpha = 0.5) + geom_smooth(method = 'lm') + geom_smooth( method = 'loess', color = 'orangered', span = 0.2 ) + labs( title = 'Объем грунтового стока', x = 'Год', y = expression(км^3) ) ``` ] ] .pull-right[ ![](08_TemporalStats_files/figure-html/wgr2-out-1.png)<!-- --> ] --- ## Локальная регрессия Нанесем несколько кривых локальной регрессии, используя разный охват данных: .pull-left[ .code-small[ ```r ggplot( tab, mapping = aes(Year1, Wgr) ) + geom_line() + geom_area(alpha = 0.5) + geom_smooth(color = 'yellow', span = 0.1) + geom_smooth(color = 'orangered', span = 0.2) + geom_smooth(color = 'darkred', span = 0.4) + labs( title = 'Объем грунтового стока', x = 'Год', y = expression(км^3) ) ``` ] ] .pull-right[ ![](08_TemporalStats_files/figure-html/loess-out-1.png)<!-- --> ] --- ## Статистики При анализе временных рядов часто используются: - Тест Манна-Кендалла на значимость линейного тренда - Тест Петтитт на точку перелома - Оценка величины тренда по методу Тейла-Сена Эти тесты доступны в пакете __trend__: .code-small[ ```r (mk = mk.test(tab$Wgr)) ## ## Mann-Kendall trend test ## ## data: tab$Wgr ## z = 7.8129, n = 75, p-value = 5.589e-15 ## alternative hypothesis: true S is not equal to 0 ## sample estimates: ## S varS tau ## 1.709000e+03 4.779167e+04 6.158559e-01 ``` ] --- ## Статистики .code-small[ ```r (pt = pettitt.test(tab$Wgr)) ## ## Pettitt's test for single change-point detection ## ## data: tab$Wgr ## U* = 1324, p-value = 4.131e-11 ## alternative hypothesis: two.sided ## sample estimates: ## probable change point at time K ## 38 (ts = sens.slope(tab$Wgr)) ## ## Sen's slope ## ## data: tab$Wgr ## z = 7.8129, n = 75, p-value = 5.589e-15 ## alternative hypothesis: true z is not equal to 0 ## 95 percent confidence interval: ## 0.01378903 0.02202750 ## sample estimates: ## Sen's slope ## 0.01810404 ``` ] --- ## Статистики .pull-left[ .code-small[ ```r ggplot(tab, aes(Year1, Wgr)) + geom_line() + geom_area(alpha = 0.5) + geom_smooth(method = 'lm') + geom_vline( xintercept = tab$Year1[pt$estimate], color = "red", size = 0.5 ) + annotate( "text", label = tab$Year1[pt$estimate], x = tab$Year1[pt$estimate] + 4, y = max(tab$Wgr), size = 4, colour = "red" ) + labs( title = 'Объем грунтового стока', x = 'Год', y = expression(км^3) ) ``` ] ] .pull-right[ Для наглядности нанесем полученную точку перелома на график: ![](08_TemporalStats_files/figure-html/pt-out-1.png)<!-- --> ] --- ## Анимация Возможности анимаций в __gganimate__ реализуются посредством добавления новых грамматик к построенному графику __ggplot2__. К числу этих грамматик относятся: - `transition_*()` — распределение данных по времени; - `view_*()` — поведение осей координат во времени; - `shadow_*()` — отображение данных, не относящихся к текущему временному срезу; - `enter_*()/exit_*()` — характер появления/исчезновения данных в процессе анимации; - `ease_aes()` — порядок смягчения (интерполяции) графических переменных в моменты перехода. --- ## Анимация В качестве первого примера используем данные реанализа NASA POWER суточного осреднения, выгрузив информацию по точкам в трех городах (Мурманск, Москва, Краснодар) за 2018 год: ```r cities = list(Мурманск = c(33, 69), Москва = c(38, 56), Краснодар = c(39, 45)) tab = purrr::imap(cities, function(coords, city){ get_power( community = "AG", lonlat = coords, pars = c("RH2M", "T2M", "PRECTOT"), dates = c("2018-01-01", "2018-12-31"), temporal_average = "DAILY" ) |> mutate(CITY = city, MONTH = month(YYYYMMDD)) }) |> bind_rows() ``` --- ## Переход по времени ```r ggplot(tab, aes(CITY, T2M)) + geom_boxplot() + labs(title = "Температура воздуха", * subtitle = 'Месяц: {round(frame_time)}') + xlab('Город') + ylab('Т, °С') + * transition_time(MONTH) ``` <img src="img/boxplot_anim.gif" width="45%" /> --- ## Переход по времени .pull-left[ ```r ggplot(tab, aes(T2M, fill = CITY)) + * geom_density(alpha = 0.5) + labs( title = "Температура воздуха", subtitle = 'Месяц: {round(frame_time)}', fill = 'Город' ) + xlab('Т, °С') + ylab('Плотность распределения') + transition_time(MONTH) ``` ] .pull-right[ Плотность распределения: ![](img/density_anim.gif)<!-- --> ] --- ## Переход по времени Данные Gapminder по соотношению продолжительности жизни и ВВП: .code-small[ ```r countries = read_excel('../r-geo-course/data/gapminder.xlsx', 2) |> select(Country = name, Region = eight_regions) |> mutate(Country = factor(Country, levels = Country[order(.$Region)])) # ВВП на душу населения gdpdf_tidy = read_sheet('1cxtzRRN6ldjSGoDzFHkB8vqPavq1iOTMElGewQnmHgg') |> pivot_longer(cols = `1764`:`2018`, names_to = 'year', values_to = 'gdp') |> rename(Country = 1) # численность населения popdf_tidy = read_sheet('1IbDM8z5XicMIXgr93FPwjgwoTTKMuyLfzU6cQrGZzH8') |> pivot_longer(cols = `1800`:`2015`, names_to = 'year', values_to = 'pop') |> rename(Country = 1) ``` ] --- ## Переход по времени Данные Gapminder по соотношению продолжительности жизни и ВВП: .code-small[ ```r # продолжительность жизни lifedf_tidy = read_sheet('1H3nzTwbn8z4lJ5gJ_WfDgCeGEXK3PVGcNjQ_U5og8eo') |> pivot_longer(cols = `1800`:`2016`, names_to = 'year', values_to = 'lifexp') |> rename(Country = 1) tab = gdpdf_tidy |> inner_join(lifedf_tidy) |> inner_join(popdf_tidy) |> inner_join(countries) |> mutate(year = as.integer(year)) |> drop_na() ``` ] --- ## Переход по времени .pull-left[ .code-small[ ```r options(scipen = 999) ggplot( tab, aes( gdp, lifexp, size = pop, color = Region ) ) + geom_point(alpha = 0.5) + scale_x_log10() + labs(title = 'Year: {round(frame_time)}') + theme_bw() + transition_time(year) ``` ] ] .pull-right[ ![](img/scatter_anim.gif)<!-- --> ] --- # Переход по состояниям - В ряде случаев вместо перехода по времени целесообразно использовать __переход по состояниям__. В частности, такой подход оказывается удобен, когда сопоставляются данные за аналогичные временные срезы разных периодов. - Подобную стратегию можно использовать для визуализации изменений внутригодичного распределения величины между годами. - Чтобы обеспечить сопоставимость аналогичных дат за разные года, необходимо сохранить у них только месяц и день. - Поскольку такого формата даты не существует, в качестве "трюка" можно просто заменить года всех дат на `\(2000\)` и записать результат в новое поле. --- # Переход по состояниям Необходимые преобразования реализуются следующим образом: ```r flt_tab = tsbl_int |> filter(!is.na(year)) |> group_by(year) |> filter(!all(is.na(level_interp))) |> ungroup() |> mutate(yDate = Date) # фиктивное поле, в котором аналогичные даты за разные года совпадают year(flt_tab$yDate) <- 2000 ``` --- # Переход по состояниям ```r anim = ggplot(flt_tab, mapping = aes(x = yDate, y = level_interp)) + geom_ribbon(aes(ymin = 0, ymax = level_interp), alpha = 0.5) + geom_line() + scale_x_date(date_breaks = "1 month", date_labels = "%b") + labs(title = "Расход воды на гидропосту Паялка", subtitle = 'Год: {closest_state}') + xlab('Дата') + ylab('куб.м/с') + theme(text = element_text(size = 18, family = 'Open Sans')) + * transition_states(year, state_length = 0) + * view_follow(fixed_y = TRUE) animate(anim, fps = 20, # число кадров в секунду nframes = 10 * length(unique(flt_tab$year)), # общее число кадров width = 800, height = 600) ``` --- # Переход по состояниям <img src="img/hydrograph_anim.gif" width="65%" /> > Текущий срез при выполнении анимации по состояниям доступен в переменной окружения `closest_state`. --- <!-- --- --> <!-- ## Локальная регрессия --> <!-- .red[Метод __локальной регрессии__ изначально был разработан для построения кривых регрессии в случае когда зависимость между переменными ведет себя сложным образом и не может быть описана в терминах традиционной линейной и нелинейной регрессии — глобальных методов.] --> <!-- .pull-left[ --> <!-- <br> --> <!-- .small[.blue[Область значений независимой переменной `\(X\)` можно покрыть конечным числом отрезков, для каждого из которых далее находят регрессию традиционным методом — как правило, линейную или квадратичную. Данный метод получил название _LOWESS (Locally weighted scatterplot smoothing)_. В дальнейшем эта аббревиатура была редуцирована до _LOESS_.]] --> <!-- ] --> <!-- .pull-right[ --> <!-- ```{r wgr2-out, ref.label='wgr2', echo=F, fig.height=5.5} --> <!-- ``` --> <!-- ] --> <!-- --- --> <!-- ## Локальная регрессия --> <!-- В классической постановке метод __LOESS__ реализуется следующим образом: --> <!-- - Пусть дано `\(n\)` точек исходных данных с координатами `\(x\)` (независимая переменная) и `\(y\)` (зависимая). --> <!-- - Задается число `\(0 < \alpha \leq 1\)`, которое обозначает долю от общего количества точек `\(n,\)` выбираемую в окрестности каждой точки для построения регрессии. --> <!-- - В абсолютном исчислении количество ближайших точек будет равно `\(r = [\alpha n]\)`, где `\([\cdot]\)` — округление до ближайшего целого. --> <!-- Тогда вес, который будет иметь каждая `\(k\)`-я точка исходных данных в уравнении регрессии для `\(i\)`-й точки исходных данных будет определяться по формуле: --> <!-- `$$w_k (x_i) = W\big((x_k - x_i)h_i^{-1}\big),$$` --> <!-- где `\(h_i\)` — расстояние до `\(r\)`-го по близости соседа точки `\(x_i\)`, а `\(W\)` — весовая функция: --> <!-- --- --> <!-- ## Локальная регрессия --> <!-- `$$w_k (x_i) = W\big((x_k - x_i)h_i^{-1}\big)$$` --> <!-- 1. `\(W(x) > 0\)` если `\(|x| < 1\)`; --> <!-- 2. `\(W(-x) = W(x)\)`; --> <!-- 3. `\(W(x)\)` невозрастающая функция для `\(x \geq 0\)`; --> <!-- 4. `\(W(x) = 0\)` если `\(|x| \geq 1\)`. --> <!-- Одним из стандартных вариантов весовой функции является _"трикубическая"_ функция, определяемая как: --> <!-- $$ --> <!-- W(x) = \begin{cases} --> <!-- (1 - |x|^3)^3, & \text{если } |x| < 1, \\ --> <!-- 0, & \text{если } |x| \geq 1. --> <!-- \end{cases} --> <!-- $$ --> <!-- Более близкие к `\(x_i\)` точки оказывают большее влияние на коэффициенты регрессии. Помимо этого за пределами расстояния `\(h_i\)` веса всех точек исходных данных будут обнуляться. --> <!-- --- --> <!-- ## Локальная регрессия --> <!-- Сглаженная оценка `\(\hat{y}_i\)` в точке `\(x_i\)` получается в виде полинома степени `\(d\)`: --> <!-- `$$\hat{y}_i = \sum_{j=0}^d \hat{\beta}_j (x_i) x_i^j,$$` --> <!-- где коэффициенты `\(\hat{\beta}_j\)` находятся методом наименьших квадратов путем минимизации ошибки: --> <!-- `$$\sum_{k=1}^n w_k (x_i) (y_k - \beta_0 - \beta_1 x_k - ... - \beta_d x_k^d)^2$$` --> <!-- Процедура поиска коэффициентов регрессии повторяется для каждой из `\(n\)` точек. --> <!-- --- --> <!-- ## Локальная регрессия --> <!-- В методе __LOESS__ используются степени регрессии `\(d = 0, 1, 2\)`. Кубические и более высокие степени полиномов на практике не применяются. --> <!-- > При степени равной 0 метод носит название _сглаживающего среднего_. --> <!-- Для построения линии локальной регрессии используйте функцию `geom_smooth()` без параметра `method` или с явным указанием параметра `method = 'loess'`: --> <!-- Поведением локальной регрессии можно управлять, задавая параметры --> <!-- - `n` (количество ближайших точек `\(r\)`), --> <!-- - `span` (доля ближайших точек `\(\alpha\)`) --> <!-- - `formula` (формула аппроксимируемой зависимости). --> <!-- По умолчанию используется регрессия первой степени (`formula = y ~ x`), значения `n = 80` и `span = 0.75`. -->