class: inverse background-image: linear-gradient(to right, rgba(150, 150, 150, .1), rgba(150, 150, 150, .4)), url("resources/hourglass.jpg") background-size: cover <style type="text/css"> /* custom.css */ .left-code { color: #777; width: 40%; height: 92%; float: left; } .right-plot { width: 58%; float: right; padding-left: 1%; } </style> .title[fable] .sticker-float[] ## Forecasting with multiple seasonality .bottom[ ### Mitchell O'Hara-Wild (<svg style="height:0.8em;top:.04em;position:relative;fill:#1da1f2;" viewBox="0 0 512 512"><path d="M459.37 151.716c.325 4.548.325 9.097.325 13.645 0 138.72-105.583 298.558-298.558 298.558-59.452 0-114.68-17.219-161.137-47.106 8.447.974 16.568 1.299 25.34 1.299 49.055 0 94.213-16.568 130.274-44.832-46.132-.975-84.792-31.188-98.112-72.772 6.498.974 12.995 1.624 19.818 1.624 9.421 0 18.843-1.3 27.614-3.573-48.081-9.747-84.143-51.98-84.143-102.985v-1.299c13.969 7.797 30.214 12.67 47.431 13.319-28.264-18.843-46.781-51.005-46.781-87.391 0-19.492 5.197-37.36 14.294-52.954 51.655 63.675 129.3 105.258 216.365 109.807-1.624-7.797-2.599-15.918-2.599-24.04 0-57.828 46.782-104.934 104.934-104.934 30.213 0 57.502 12.67 76.67 33.137 23.715-4.548 46.456-13.32 66.599-25.34-7.798 24.366-24.366 44.833-46.132 57.827 21.117-2.273 41.584-8.122 60.426-16.243-14.292 20.791-32.161 39.308-52.628 54.253z"/></svg>[@mitchoharawild](https://twitter.com/mitchoharawild)) ### 25 November 2020 ### Slides @ [slides.mitchelloharawild.com/nhs2020](https://slides.mitchelloharawild.com/nhs2020) ] --- class: inverse background-image: linear-gradient(to right, rgba(150, 150, 150, .1), rgba(150, 150, 150, .4)), url("resources/timezone.jpg") background-size: cover .title[Hello!] --- class: center background-image: linear-gradient(to right, rgba(255, 255, 255, .1), rgba(150, 150, 150, .1)), url("resources/season.jpg") background-size: cover ## Understanding multiple seasonality --- background-image: linear-gradient(to right, rgba(255, 255, 255, .1), rgba(150, 150, 150, .1)), url("resources/season.jpg") background-size: cover .center[ ## What is seasonality? ] -- > A consistent pattern which repeats over a fixed period of time. -- .box-12.opaque[ <img src="figure/ped-seasonality-1.svg" style="display: block; margin: auto;" /> ] --- background-image: linear-gradient(to right, rgba(255, 255, 255, .1), rgba(150, 150, 150, .1)), url("resources/season.jpg") background-size: cover .center[ ## What is *multiple* seasonality? ] -- > Two or more seasonal patterns. -- .box-12.opaque[ <img src="figure/ped-multi-seasonality-1.svg" style="display: block; margin: auto;" /> ] --- class: inverse background-image: linear-gradient(to right, rgba(50, 50, 50, .5), rgba(50, 50, 50, .4)), url("resources/tree.jpg") background-size: cover .center[ ## How can seasonal patterns be modelled? ] <hr> -- <br><br> .center[ ### Seasonal dummy variables ### Lagged terms ### Fourier terms ### Exogenous regressors ] --- ## Modelling seasonal patterns with... ## Seasonal dummy variables .pull-left[ Produces distinct intercepts for each discrete season. ✅ Simplest to interpret and produce. ❌ Requires many terms for longer periods. ❌ Poorly handles complex seasonality. ] .pull-right[
$$d_{1,t}$$
$$d_{2,t}$$
$$d_{3,t}$$
$$d_{4,t}$$
$$d_{5,t}$$
$$d_{6,t}$$
Monday
1
0
0
0
0
0
Tuesday
0
1
0
0
0
0
Wednesday
0
0
1
0
0
0
Thursday
0
0
0
1
0
0
Friday
0
0
0
0
1
0
Saturday
0
0
0
0
0
1
Sunday
0
0
0
0
0
0
Monday
1
0
0
0
0
0
...
...
...
...
...
...
...
] --- ## Modelling seasonal patterns with... ## Lagged terms .pull-left[ *'Today will be similar to yesterday.'* ✅ Very difficult to beat accuracy of. ❌ Limited forecast horizon. ❌ Cannot handle non-integer seasonality. ] .pull-right[ <img src="figure/lag-plot-1.svg" style="display: block; margin: auto;" /> ] --- ## Modelling seasonal patterns with... ## Fourier terms `$$x_{1,t} = \sin\left(\textstyle\frac{2\pi t}{m}\right), x_{2,t} = \cos\left(\textstyle\frac{2\pi t}{m}\right),$$` `$$x_{3,t} = \sin\left(\textstyle\frac{4\pi t}{m}\right), x_{4,t} = \cos\left(\textstyle\frac{4\pi t}{m}\right), \ldots$$` <img src="figure/fourier-plot-1.svg" style="display: block; margin: auto;" /> --- ## Modelling seasonal patterns with... ## Fourier terms .pull-left[ Captures seasonal distribution with a continuous periodic functions. ✅ Handles non-integer seasonality (days in year, weeks in month). ✅ Flexibility of seasonal pattern controllable by number of harmonics. ❌ More complicated to interpret. ] .pull-right[ <img src="figure/fourier-plot-2-1.svg" style="display: block; margin: auto;" /> ] --- ## Modelling seasonal patterns with... ## Exogenous regressors .pull-left[ Capture seasonality by providing underlying factors like temperature, and events. ✅ Better handling of non-seasonal anomalies. ❌ Data may be hard to find. ❌ Future values must be forecasted. ] .pull-right[ <img src="figure/exog-plot-1.svg" style="display: block; margin: auto;" /> ] --- class: inverse background-image: linear-gradient(to right, rgba(50, 50, 50, .5), rgba(50, 50, 50, .4)), url("resources/tree.jpg") background-size: cover .center[ ## Which seasonal method should I use? ] <hr> -- <br> .center[ ### For *multiple* seasonality, *multiple* methods may be needed. <br> ### Use methods that are available and work best for your data! ] --- ## Models for forecasting multiple seasonalities Capability for including seasonal terms. <table class="table"> <thead> <tr> <th>Model</th> <th>Dummies</th> <th>Lags</th> <th>Fourier</th> <th>Regressors</th> </tr> </thead> <tbody> <tr> <td>MLR</td> <td class="success">Yes</td> <td class="success">Yes</td> <td class="success">Yes</td> <td class="success">Yes</td> </tr> <tr> <td>GAM</td> <td class="success">Yes</td> <td class="success">Yes</td> <td class="success">Yes</td> <td class="success">Yes</td> </tr> <tr> <td>Prophet</td> <td class="danger">No</td> <td class="success">Yes</td> <td class="success">Yes</td> <td class="success">Yes</td> </tr> <tr> <td>DSHW</td> <td class="success">Yes</td> <td class="danger">No</td> <td class="danger">No</td> <td class="danger">No</td> </tr> <tr> <td>BATS/TBATS</td> <td class="success">Yes</td> <td class="warning">Partial</td> <td class="success">Yes</td> <td class="danger">No</td> </tr> <tr> <td>FASSTER</td> <td class="success">Yes</td> <td class="success">Yes</td> <td class="success">Yes</td> <td class="success">Yes</td> </tr> <tr> <td>STL + ???</td> <td class="warning">???</td> <td class="warning">???</td> <td class="warning">???</td> <td class="warning">???</td> </tr> </tbody> </table> .footnote[(This *non-exhaustive* list includes the most common statistical methods only.)] --- ## Models for forecasting multiple seasonalities <br> <table class="table"> <thead> <tr> <th>Model</th> <th>Flexible</th> <th>Speed</th> <th>Accuracy</th> <th>Decompose Components</th> <th>Evolving Terms</th> </tr> </thead> <tbody> <tr> <td>Regression<br>(MLR, GAM, Prophet)</td> <td></td> <td></td> <td></td> <td></td> <td></td> </tr> <tr> <td>State Space<br>(DSHW, TBATS, FASSTER)</td> <td></td> <td></td> <td></td> <td></td> <td></td> </tr> <tr> <td>Decomposition<br>(STL + ???)</td> <td></td> <td> </td> <td></td> <td></td> <td></td> </tr> </tbody> </table> --- ## Models for forecasting multiple seasonalities <br> <table class="table"> <thead> <tr> <th>Model</th> <th>Flexible</th> <th>Speed</th> <th>Accuracy</th> <th>Decompose Components</th> <th>Evolving Terms</th> </tr> </thead> <tbody> <tr> <td>Regression<br>(MLR, GAM, Prophet)</td> <td class="warning">Depends</td> <td class="warning">Depends</td> <td class="warning">Depends</td> <td class="warning">Depends</td> <td class="warning">Depends</td> </tr> <tr> <td>State Space<br>(DSHW, TBATS, FASSTER)</td> <td class="warning">Depends</td> <td class="warning">Depends</td> <td class="warning">Depends</td> <td class="warning">Depends</td> <td class="warning">Depends</td> </tr> <tr> <td>Decomposition<br>(STL + ???)</td> <td class="warning">Depends</td> <td class="warning">Depends </td> <td class="warning">Depends</td> <td class="warning">Depends</td> <td class="warning">Depends</td> </tr> </tbody> </table> --- ## Models for forecasting multiple seasonalities <br> <table class="table"> <thead> <tr> <th>Model</th> <th>Flexible</th> <th>Speed</th> <th>Accuracy</th> <th>Decompose Components</th> <th>Evolving Terms</th> </tr> </thead> <tbody> <tr> <td>Regression<br>(MLR, GAM, Prophet)</td> <td class="success">Yes</td> <td class="success">Fast</td> <td class="warning">Okay</td> <td class="warning">Depends</td> <td class="danger">No</td> </tr> <tr> <td>State Space<br>(DSHW, TBATS, FASSTER)</td> <td class="danger">No</td> <td class="danger">Slow</td> <td class="warning">Okay</td> <td class="success">Yes</td> <td class="success">Yes</td> </tr> <tr> <td>Decomposition<br>(STL + ???)</td> <td class="warning">Depends</td> <td class="warning">Moderate</td> <td class="warning">Okay</td> <td class="success">Yes</td> <td class="warning">Depends</td> </tr> </tbody> </table> .footnote[(This table is a *massive* generalisation, everything depends on data and model complexity.)] --- class: center ## Forecasting multiple seasonality in R .sticker[] -- .animated.fadeIn[ .sticker[] .sticker[] .sticker[] ## [tidyverts.org](http://www.tidyverts.org) ] --- class: top
--- class: inverse, center background-image: linear-gradient(to right, rgba(50, 50, 50, .5), rgba(50, 50, 50, .5)), url("resources/disk.jpg") background-size: cover .title[Tidy] <hr> <br> ## Data preparation --- class: top # Hospital admissions ```r fs::dir_tree("~/github/webinar-msc-nhs/") ``` ``` #> ~/github/webinar-msc-nhs/ #> ├── README.md #> ├── data #> │ ├── ae_uk.csv #> │ ├── holiday.csv #> │ └── temp.csv #> ├── rscript.R #> └── webinar-msc-nhs.Rproj ``` .footnote[*https://github.com/bahmanrostamitabar/webinar-msc-nhs/*] --- class: top # Admissions event data ```r admissions <- read_csv( file = "~/github/webinar-msc-nhs/data/ae_uk.csv", col_types = cols(arrival_time=col_datetime(format = "%d/%m/%Y %H:%M")) ) %>% select(-ID) # ID is a redundant row identifier print(admissions) ``` ``` #> # A tibble: 773,779 x 3 #> gender type_injury arrival_time #> <chr> <chr> <dttm> #> 1 male minor 2010-01-01 00:35:00 #> 2 male minor 2010-01-01 00:57:00 #> 3 male minor 2010-01-01 00:42:00 #> 4 male minor 2010-01-01 00:02:00 #> 5 male minor 2010-01-01 00:24:00 #> 6 female major 2010-01-01 00:48:00 #> 7 female minor 2010-01-01 00:10:00 #> 8 male minor 2010-01-01 00:39:00 #> 9 female minor 2010-01-01 00:56:00 #> 10 female minor 2010-01-01 00:05:00 #> # … with 773,769 more rows ``` --- class: top # Holiday data ```r holidays <- read_csv( file = "~/github/webinar-msc-nhs/data/holiday.csv", col_types = cols(date=col_date(format = "%d/%m/%Y")) ) %>% filter(!is.na(date)) # Some extra empty rows exist that need removing print(holidays) ``` ``` #> # A tibble: 2,282 x 7 #> date public_holiday_type public_holiday school_holiday school_holiday_type celebratory_day celebratory_day_ty… #> <date> <chr> <dbl> <dbl> <chr> <dbl> <chr> #> 1 2010-01-01 New Years Day 1 0 0 0 0 #> 2 2010-01-02 0 0 1 Winter School Holid… 0 0 #> 3 2010-01-03 0 0 1 Winter School Holid… 0 0 #> 4 2010-01-04 0 0 0 0 0 0 #> 5 2010-01-05 0 0 0 0 0 0 #> 6 2010-01-06 0 0 0 0 0 0 #> 7 2010-01-07 0 0 0 0 0 0 #> 8 2010-01-08 0 0 0 0 0 0 #> 9 2010-01-09 0 0 0 0 0 0 #> 10 2010-01-10 0 0 0 0 0 0 #> # … with 2,272 more rows ``` --- class: top # Temperature data ```r temperatures <- read_csv( file = "~/github/webinar-msc-nhs/data/temp.csv", col_types = cols(date=col_date(format = "%d/%m/%Y")) ) print(temperatures) ``` ``` #> # A tibble: 2,282 x 3 #> date actual_temp forecast_temp #> <date> <dbl> <dbl> #> 1 2010-01-01 4 4.3 #> 2 2010-01-02 4.8 5 #> 3 2010-01-03 6 6.4 #> 4 2010-01-04 3 1.9 #> 5 2010-01-05 4.1 3.8 #> 6 2010-01-06 3.5 2.2 #> 7 2010-01-07 1.6 1 #> 8 2010-01-08 3 2.7 #> 9 2010-01-09 1.5 1 #> 10 2010-01-10 1.4 1.3 #> # … with 2,272 more rows ``` --- class: top .sticker-float[] # Tidy temporal data structures
--- class: top .sticker-float[] # Tidy temperature data Identify the index, key and measurement variables: ```r temperatures <- as_tsibble(temperatures, * index = ???, key = ???) ``` ``` #> # A tibble: 2,282 x 3 #> date actual_temp forecast_temp #> <date> <dbl> <dbl> #> 1 2010-01-01 4 4.3 #> 2 2010-01-02 4.8 5 #> 3 2010-01-03 6 6.4 #> 4 2010-01-04 3 1.9 #> 5 2010-01-05 4.1 3.8 #> 6 2010-01-06 3.5 2.2 #> 7 2010-01-07 1.6 1 #> 8 2010-01-08 3 2.7 #> 9 2010-01-09 1.5 1 #> 10 2010-01-10 1.4 1.3 #> # … with 2,272 more rows ``` --- class: top .sticker-float[] # Tidy temperature data Identify the index, key and measurement variables: ```r temperatures <- as_tsibble(temperatures, * index = date, key = NULL) ``` ``` #> # A tsibble: 2,282 x 3 [1D] #> date actual_temp forecast_temp #> <date> <dbl> <dbl> #> 1 2010-01-01 4 4.3 #> 2 2010-01-02 4.8 5 #> 3 2010-01-03 6 6.4 #> 4 2010-01-04 3 1.9 #> 5 2010-01-05 4.1 3.8 #> 6 2010-01-06 3.5 2.2 #> 7 2010-01-07 1.6 1 #> 8 2010-01-08 3 2.7 #> 9 2010-01-09 1.5 1 #> 10 2010-01-10 1.4 1.3 #> # … with 2,272 more rows ``` --- class: top .sticker-float[] # Tidy holiday data Identify the index, key and measurement variables: ```r holidays <- as_tsibble(holidays, * index = ???, key = ???) ``` ``` #> # A tibble: 2,282 x 7 #> date public_holiday_type public_holiday school_holiday school_holiday_type celebratory_day celebratory_day_ty… #> <date> <chr> <dbl> <dbl> <chr> <dbl> <chr> #> 1 2010-01-01 New Years Day 1 0 0 0 0 #> 2 2010-01-02 0 0 1 Winter School Holid… 0 0 #> 3 2010-01-03 0 0 1 Winter School Holid… 0 0 #> 4 2010-01-04 0 0 0 0 0 0 #> 5 2010-01-05 0 0 0 0 0 0 #> 6 2010-01-06 0 0 0 0 0 0 #> 7 2010-01-07 0 0 0 0 0 0 #> 8 2010-01-08 0 0 0 0 0 0 #> 9 2010-01-09 0 0 0 0 0 0 #> 10 2010-01-10 0 0 0 0 0 0 #> # … with 2,272 more rows ``` --- class: top .sticker-float[] # Tidy holiday data Identify the index, key and measurement variables: ```r holidays <- as_tsibble(holidays, * index = date, key = NULL) ``` ``` #> # A tsibble: 2,282 x 7 [1D] #> date public_holiday_type public_holiday school_holiday school_holiday_type celebratory_day celebratory_day_ty… #> <date> <chr> <dbl> <dbl> <chr> <dbl> <chr> #> 1 2010-01-01 New Years Day 1 0 0 0 0 #> 2 2010-01-02 0 0 1 Winter School Holid… 0 0 #> 3 2010-01-03 0 0 1 Winter School Holid… 0 0 #> 4 2010-01-04 0 0 0 0 0 0 #> 5 2010-01-05 0 0 0 0 0 0 #> 6 2010-01-06 0 0 0 0 0 0 #> 7 2010-01-07 0 0 0 0 0 0 #> 8 2010-01-08 0 0 0 0 0 0 #> 9 2010-01-09 0 0 0 0 0 0 #> 10 2010-01-10 0 0 0 0 0 0 #> # … with 2,272 more rows ``` --- class: top .sticker-float[] # Tidy admissions data Identify the index, key and measurement variables: ```r admissions <- as_tsibble(admissions, * index = ???, key = ???) ``` ``` #> # A tibble: 773,779 x 3 #> gender type_injury arrival_time #> <chr> <chr> <dttm> #> 1 male minor 2010-01-01 00:35:00 #> 2 male minor 2010-01-01 00:57:00 #> 3 male minor 2010-01-01 00:42:00 #> 4 male minor 2010-01-01 00:02:00 #> 5 male minor 2010-01-01 00:24:00 #> 6 female major 2010-01-01 00:48:00 #> 7 female minor 2010-01-01 00:10:00 #> 8 male minor 2010-01-01 00:39:00 #> 9 female minor 2010-01-01 00:56:00 #> 10 female minor 2010-01-01 00:05:00 #> # … with 773,769 more rows ``` --- class: top .sticker-float[] # Tidy admissions data Identify the index, key and measurement variables: ```r admissions <- as_tsibble(admissions, * index = arrival_time, key = c(gender, type_injury)) ``` ``` #> Error: A valid tsibble must have distinct rows identified by key and index. #> ℹ Please use `duplicates()` to check the duplicated rows. ``` -- ```r duplicates(admissions, index = arrival_time, key = c(gender, type_injury)) ``` ``` #> # A tibble: 82,578 x 3 #> gender type_injury arrival_time #> <chr> <chr> <dttm> #> 1 female minor 2010-01-01 01:37:00 #> 2 female minor 2010-01-01 01:37:00 #> 3 female minor 2010-01-01 10:07:00 #> 4 female minor 2010-01-01 10:07:00 #> 5 female minor 2010-01-01 11:15:00 #> 6 female minor 2010-01-01 11:15:00 #> 7 male minor 2010-01-01 12:01:00 #> 8 female minor 2010-01-01 12:24:00 #> 9 female minor 2010-01-01 12:24:00 #> 10 male minor 2010-01-01 12:01:00 #> # … with 82,568 more rows ``` --- class: top .sticker-float[] # Tidy admissions data Identify the index, key and measurement variables: ```r admissions <- admissions %>% * count(gender, type_injury, arrival_time, name = "arrivals") %>% as_tsibble(index = arrival_time, key = c(gender, type_injury)) ``` ``` #> # A tsibble: 731,906 x 4 [1m] <UTC> #> # Key: gender, type_injury [4] #> gender type_injury arrival_time arrivals #> <chr> <chr> <dttm> <int> #> 1 female major 2010-01-01 00:48:00 1 #> 2 female major 2010-01-01 04:19:00 1 #> 3 female major 2010-01-01 05:57:00 1 #> 4 female major 2010-01-01 12:17:00 1 #> 5 female major 2010-01-01 12:42:00 1 #> 6 female major 2010-01-01 13:34:00 1 #> 7 female major 2010-01-01 16:59:00 1 #> 8 female major 2010-01-02 23:56:00 1 #> 9 female major 2010-01-03 00:24:00 1 #> 10 female major 2010-01-03 13:14:00 1 #> # … with 731,896 more rows ``` --- class: top .sticker-float[] # Bringing the data together Can be a non-trivial task! -- ```r interval(admissions) ``` ``` #> <interval[1]> #> [1] 1m ``` ```r interval(holidays) ``` ``` #> <interval[1]> #> [1] 1D ``` ```r interval(temperatures) ``` ``` #> <interval[1]> #> [1] 1D ``` --- class: top .sticker-float[] # Choosing appropriate temporal intervals What temporal granularity should `admissions` have? -- A balancing act between: * The temporal detail needed for the forecasts * Enough detail/signal in the data for a good model -- This choice has a big impact on the seasonal structures that will exist in the data! --- class: top .sticker-float[] # Choosing appropriate temporal intervals Option A: Keep it as is (event level) ```r admissions %>% filter(yearmonth(arrival_time) == yearmonth("2010 Jan")) %>% * summarise(arrivals = sum(arrivals)) %>% * fill_gaps(arrivals = 0) %>% autoplot() # # ``` <img src="figure/temporal-a-1.svg" style="display: block; margin: auto;" /> --- class: top .sticker-float[] # Choosing appropriate temporal intervals Option B: 30 minutes ```r admissions %>% filter(yearmonth(arrival_time) == yearmonth("2010 Jan")) %>% * index_by(time = floor_date(arrival_time, "30 minutes")) %>% summarise(arrivals = sum(arrivals)) %>% fill_gaps(arrivals = 0) %>% autoplot() ``` <img src="figure/temporal-b-1.svg" style="display: block; margin: auto;" /> --- class: top .sticker-float[] # Choosing appropriate temporal intervals Option C: 1 hour ```r admissions %>% filter(yearmonth(arrival_time) == yearmonth("2010 Jan")) %>% * index_by(time = floor_date(arrival_time, "1 hour")) %>% summarise(arrivals = sum(arrivals)) %>% fill_gaps(arrivals = 0) %>% autoplot() ``` <img src="figure/temporal-c-1.svg" style="display: block; margin: auto;" /> --- class: top .sticker-float[] # Choosing appropriate temporal intervals Option D: 1 day ```r admissions %>% filter(year(arrival_time) == 2010) %>% * index_by(time = as.Date(floor_date(arrival_time, "1 day"))) %>% summarise(arrivals = sum(arrivals)) %>% fill_gaps(arrivals = 0) %>% autoplot() ``` <img src="figure/temporal-d-1.svg" style="display: block; margin: auto;" /> --- class: top .sticker-float[] # Choosing appropriate temporal intervals Option E: 1 month ```r admissions %>% * index_by(time = yearmonth(floor_date(arrival_time, "1 day"))) %>% summarise(arrivals = sum(arrivals)) %>% fill_gaps(arrivals = 0) %>% autoplot() # ``` <img src="figure/temporal-e-1.svg" style="display: block; margin: auto;" /> --- class: top .sticker-float[] # Choosing appropriate temporal intervals So what's best? > * Has the temporal detail needed for the forecasts > * Enough detail/signal in the data for a good model -- .center[*It depends!*] -- .center[ # 🙄 ] -- .center[For a multiple seasonality webinar, probably B or C. Let's go with hourly data (C) for as it has better signal.] --- class: top .sticker-float[] # Choosing appropriate temporal intervals ```r admissions <- admissions %>% * index_by(time = floor_date(arrival_time, "1 hour")) %>% summarise(arrivals = sum(arrivals)) %>% ungroup() %>% fill_gaps(arrivals = 0) ``` -- ``` #> # A tsibble: 54,768 x 2 [1h] <UTC> #> time arrivals #> <dttm> <dbl> #> 1 2010-01-01 00:00:00 13 #> 2 2010-01-01 01:00:00 15 #> 3 2010-01-01 02:00:00 12 #> 4 2010-01-01 03:00:00 15 #> 5 2010-01-01 04:00:00 11 #> 6 2010-01-01 05:00:00 7 #> 7 2010-01-01 06:00:00 6 #> 8 2010-01-01 07:00:00 4 #> 9 2010-01-01 08:00:00 6 #> 10 2010-01-01 09:00:00 10 #> # … with 54,758 more rows ``` --- class: middle background-image: linear-gradient(to right, rgba(50, 50, 50, .5), rgba(50, 50, 50, .5)), url("resources/many_clocks.jpg") background-size: cover .box-12[ ## Aside Advanced forecasters may even choose multiple granularities! <br> Forecasts at many levels can be combined using **temporal reconciliation** to get more accurate results. This process ensures the hourly forecasts add to the daily forecasts, and so on. <br> <svg style="height:0.8em;top:.04em;position:relative;" viewBox="0 0 576 512"><path d="M528 0H48C21.5 0 0 21.5 0 48v320c0 26.5 21.5 48 48 48h192l-16 48h-72c-13.3 0-24 10.7-24 24s10.7 24 24 24h272c13.3 0 24-10.7 24-24s-10.7-24-24-24h-72l-16-48h192c26.5 0 48-21.5 48-48V48c0-26.5-21.5-48-48-48zm-16 352H64V64h448v288z"/></svg> ISF 2020 slides: [slides.mitchelloharawild.com/isf2020](https://slides.mitchelloharawild.com/isf2020) <svg style="height:0.8em;top:.04em;position:relative;" viewBox="0 0 576 512"><path d="M336.2 64H47.8C21.4 64 0 85.4 0 111.8v288.4C0 426.6 21.4 448 47.8 448h288.4c26.4 0 47.8-21.4 47.8-47.8V111.8c0-26.4-21.4-47.8-47.8-47.8zm189.4 37.7L416 177.3v157.4l109.6 75.5c21.2 14.6 50.4-.3 50.4-25.8V127.5c0-25.4-29.1-40.4-50.4-25.8z"/></svg> Video recording: [youtu.be/6D7rNHZ5E-Q](https://youtu.be/6D7rNHZ5E-Q?t=1120) ] --- class: top .sticker-float[] # Bringing the data together ```r interval(admissions) ``` ``` #> <interval[1]> #> [1] 1h ``` ```r interval(holidays) ``` ``` #> <interval[1]> #> [1] 1D ``` ```r interval(temperatures) ``` ``` #> <interval[1]> #> [1] 1D ``` -- It's okay for `holidays` to be daily as the holiday affects the entire day. Ideally, `temperatures` would be hourly as it varies throughout the day. --- class: top .sticker-float[] # Bringing the data together ```r admissions %>% # Compute the common variable mutate(date = as.Date(time)) ``` ``` #> # A tsibble: 54,768 x 3 [1h] <UTC> #> time arrivals date #> <dttm> <dbl> <date> #> 1 2010-01-01 00:00:00 13 2010-01-01 #> 2 2010-01-01 01:00:00 15 2010-01-01 #> 3 2010-01-01 02:00:00 12 2010-01-01 #> 4 2010-01-01 03:00:00 15 2010-01-01 #> 5 2010-01-01 04:00:00 11 2010-01-01 #> 6 2010-01-01 05:00:00 7 2010-01-01 #> 7 2010-01-01 06:00:00 6 2010-01-01 #> 8 2010-01-01 07:00:00 4 2010-01-01 #> 9 2010-01-01 08:00:00 6 2010-01-01 #> 10 2010-01-01 09:00:00 10 2010-01-01 #> # … with 54,758 more rows ``` --- class: top .sticker-float[] # Bringing the data together ```r admissions %>% # Compute the common variable mutate(date = as.Date(time)) %>% # Join in holiday data (keeping only 'public_holiday' for simplicity) left_join(holidays %>% select(public_holiday), by = "date") ``` ``` #> # A tsibble: 54,768 x 4 [1h] <UTC> #> time arrivals date public_holiday #> <dttm> <dbl> <date> <dbl> #> 1 2010-01-01 00:00:00 13 2010-01-01 1 #> 2 2010-01-01 01:00:00 15 2010-01-01 1 #> 3 2010-01-01 02:00:00 12 2010-01-01 1 #> 4 2010-01-01 03:00:00 15 2010-01-01 1 #> 5 2010-01-01 04:00:00 11 2010-01-01 1 #> 6 2010-01-01 05:00:00 7 2010-01-01 1 #> 7 2010-01-01 06:00:00 6 2010-01-01 1 #> 8 2010-01-01 07:00:00 4 2010-01-01 1 #> 9 2010-01-01 08:00:00 6 2010-01-01 1 #> 10 2010-01-01 09:00:00 10 2010-01-01 1 #> # … with 54,758 more rows ``` --- class: top .sticker-float[] # Bringing the data together ```r hospital <- admissions %>% # Compute the common variable mutate(date = as.Date(time)) %>% # Join in holiday data (keeping only 'public_holiday' for simplicity) left_join(holidays %>% select(public_holiday), by = "date") %>% # Join in temperature data (keeping only 'actual_temp' for simplicity) left_join(temperatures %>% select(actual_temp), by = "date") ``` ``` #> # A tsibble: 54,768 x 5 [1h] <UTC> #> time arrivals date public_holiday actual_temp #> <dttm> <dbl> <date> <dbl> <dbl> #> 1 2010-01-01 00:00:00 13 2010-01-01 1 4 #> 2 2010-01-01 01:00:00 15 2010-01-01 1 4 #> 3 2010-01-01 02:00:00 12 2010-01-01 1 4 #> 4 2010-01-01 03:00:00 15 2010-01-01 1 4 #> 5 2010-01-01 04:00:00 11 2010-01-01 1 4 #> 6 2010-01-01 05:00:00 7 2010-01-01 1 4 #> 7 2010-01-01 06:00:00 6 2010-01-01 1 4 #> 8 2010-01-01 07:00:00 4 2010-01-01 1 4 #> 9 2010-01-01 08:00:00 6 2010-01-01 1 4 #> 10 2010-01-01 09:00:00 10 2010-01-01 1 4 #> # … with 54,758 more rows ``` --- class: inverse, center background-image: linear-gradient(to right, rgba(50, 50, 50, .5), rgba(50, 50, 50, .5)), url("resources/graph.jpg") background-size: cover .title[Visualise] <hr> <br> ## Data exploration --- # Visualising multiple seasonality Plotting hourly data observed over six years can be tricky! ```r hospital %>% autoplot(arrivals) ``` <img src="figure/admissions-plot-1.png" style="display: block; margin: auto;" /> --- # Visualising multiple seasonality One possibility is to view a portion of the data (but you might miss things!) ```r library(dygraphs) tsbox::ts_xts(hospital) %>% dygraph() %>% dyRangeSelector(dateWindow = c("2010-01-01", "2010-02-01")) ```
--- .sticker-float[] # Visualising multiple seasonality Seeing the big picture may require aggregating the data. ```r library(feasts) hospital %>% index_by(date = yearmonth(time)) %>% summarise(arrivals = sum(arrivals)) %>% gg_subseries(arrivals, period = "year") ``` <img src="figure/subseries-month-1.svg" style="display: block; margin: auto;" /> --- .sticker-float[] # Visualising multiple seasonality ```r hospital %>% index_by(date = as.Date(time)) %>% summarise(arrivals = sum(arrivals)) %>% gg_subseries(arrivals, period = "week") ``` <img src="figure/subseries-week-1.png" style="display: block; margin: auto;" /> --- .sticker-float[] # Visualising multiple seasonality ```r hospital %>% gg_subseries(arrivals, period = "day") ``` <img src="figure/subseries-day-1.png" style="display: block; margin: auto;" /> --- .sticker-float[] # Visualising multiple seasonality ```r hospital %>% ggplot(aes(x = hour(time), y = arrivals)) + geom_jitter(alpha = 0.01) + geom_smooth() + facet_grid(cols = vars(wday(time, label = TRUE, week_start = 1))) ``` <img src="figure/ggplot-day-1.png" style="display: block; margin: auto;" /> --- .sticker-float[] # Visualising multiple seasonality ```r *library(sugrrants) p <- hospital %>% mutate(date = as.Date(time)) %>% filter(year(time) == 2012) %>% * frame_calendar(x = hour(time), y = arrivals, date = date) %>% ggplot(aes(x = `.hour(time)`, y = .arrivals, group = date)) + geom_line() prettify(p) ``` --- .sticker-float[] # Visualising multiple seasonality <img src="figure/cal-plot-1.svg" style="display: block; margin: auto;" /> --- class: inverse, center background-image: linear-gradient(to right, rgba(50, 50, 50, .7), rgba(50, 50, 50, .7)), url("resources/books.jpg") background-size: cover .title[Making forecasts] <hr> <br> ## Data modelling .left[ .footnote[ ⚠️ The specific models/forecasts shown are in no way recommended or good. They simply demonstrate how multiple seasonalities can be modelled.] ] --- ## Models for forecasting multiple seasonalities <br> <table class="table"> <thead> <tr> <th>Model</th> <th>Flexible</th> <th>Speed</th> <th>Accuracy</th> <th>Decompose Components</th> <th>Evolving Terms</th> </tr> </thead> <tbody> <tr> <td>Regression<br>(MLR, GAM, <b>Prophet</b>)</td> <td class="success">Yes</td> <td class="success">Fast</td> <td class="warning">Okay</td> <td class="warning">Depends</td> <td class="danger">No</td> </tr> <tr> <td>State Space<br>(DSHW, TBATS, <b>FASSTER</b>)</td> <td class="danger">No</td> <td class="danger">Slow</td> <td class="warning">Okay</td> <td class="success">Yes</td> <td class="success">Yes</td> </tr> <tr> <td>Decomposition<br>(<b>STL + ???</b>)</td> <td class="warning">Depends</td> <td class="warning">Moderate</td> <td class="warning">Okay</td> <td class="success">Yes</td> <td class="warning">Depends</td> </tr> </tbody> </table> .footnote[(This table is a *massive* generalisation, everything depends on data and model complexity.)] --- ## Regression modelling <br> <table class="table"> <thead> <tr> <th>Model</th> <th>Flexible</th> <th>Speed</th> <th>Accuracy</th> <th>Decompose Components</th> <th>Evolving Terms</th> </tr> </thead> <tbody> <tr> <td>Regression<br>(MLR, GAM, <b>Prophet</b>)</td> <td class="success">Yes</td> <td class="success">Fast</td> <td class="warning">Okay</td> <td class="warning">Depends</td> <td class="danger">No</td> </tr> </tbody> </table> Uses time series features / regressors with known future values to produce forecasts. ✅ Generally fast, as model estimation is time independent. ❌ Model behaviour does not update over time for changing dynamics. Highly recommended read on GAM forecasting: https://petolau.github.io/Analyzing-double-seasonal-time-series-with-GAM-in-R/ --- ## Prophet Facebook's forecasting model based on a GAM which uses regressors for growth, seasonality, and holidays. $$ y(t) = g(t) + s_1(t) + s_2(t) + h(t) + \varepsilon_t $$ where `\(g(t)\)` is the growth function, `\(s(t)\)` are fourier terms, and `\(h(t)\)` are holiday regressors. Available directly in the [`{prophet}`](https://facebook.github.io/prophet/) package, or a fable interface from [`{fable.prophet}`](http://pkg.mitchelloharawild.com/fable.prophet/). --- ## Specifying a prophet model Using default settings ```r library(fable.prophet) hospital %>% model(prophet(arrivals)) %>% forecast(h = "2 weeks") %>% autoplot(tail(hospital, 24*7*4)) ``` <img src="figure/prophet-1.svg" style="display: block; margin: auto;" /> --- ## Specifying a prophet model Customising the `season()` options. ```r library(fable.prophet) hospital %>% model(prophet(arrivals ~ season("day", 7, type = "mult") + season("week", 4))) %>% forecast(h = "2 weeks") %>% autoplot(tail(hospital, 24*7*4)) ``` <img src="figure/prophet-custom-1.svg" style="display: block; margin: auto;" /> --- ## State space modelling <br> <table class="table"> <thead> <tr> <th>Model</th> <th>Flexible</th> <th>Speed</th> <th>Accuracy</th> <th>Decompose Components</th> <th>Evolving Terms</th> </tr> </thead> <tbody> <tr> <td>State Space<br>(DSHW, TBATS, <b>FASSTER</b>)</td> <td class="danger">No</td> <td class="danger">Slow</td> <td class="warning">Okay</td> <td class="success">Yes</td> <td class="success">Yes</td> </tr> </tbody> </table> Iterates 'state equations' over time to produce forecasts. ✅ Model behaviour updates for changes in the data. ❌ Relatively slow, as model estimation is done sequentially over time. --- ## FASSTER A model which provides flexible regression-style model specification with evolving terms of state space models via a DLM. `\begin{align*} y_t &= F_t\theta_t + v_t, & v_t &\sim \mathcal{N}(0,V)\\ \theta_t &= G\theta_{t-1} + w_t, & w_t&\sim \mathcal{N}(0,W) \end{align*}` The underlying states are given by `\(\theta\)`, which can be thought of as components linearly combined by `\(F_t\)` to produce the response `\(y_t\)`. The `\(G\)` matrix defines the behaviour of each state (seasonality, trend, regression, ...). Available with direct fable compatibility in the [`{fasster}`](http://fasster.tidyverts.org/) package. --- ## Specifying seasonality with FASSTER FASSTER allows all aforementioned seasonal terms to be specified.
Seasonal term
Special
Seasonal dummies
y ~ season(p)
Lag terms
y ~ lag(y, p)
Fourier terms
y ~ fourier(p, K)
Regressors
y ~ x
--- ## Specifying several seasonalities with FASSTER ```r library(fasster) ped_train <- pedestrian %>% filter(Sensor == "Southern Cross Station", yearmonth(Date) == yearmonth("2015 July")) ped_train %>% model(fasster(Count ~ trend(1) + fourier("week", 3) + (wday(Date_Time) %in% c(1,7)) %S% fourier("day", 10))) %>% forecast(h = "2 weeks") %>% autoplot(ped_train) ``` <img src="figure/fasster-simple-1.svg" style="display: block; margin: auto;" /> --- ## Decomposition modelling <br> <table class="table"> <thead> <tr> <th>Model</th> <th>Flexible</th> <th>Speed</th> <th>Accuracy</th> <th>Decompose Components</th> <th>Evolving Terms</th> </tr> </thead> <tbody> <tr> <td>Decomposition<br>(<b>STL + ???</b>)</td> <td class="warning">Depends</td> <td class="warning">Moderate</td> <td class="warning">Okay</td> <td class="success">Yes</td> <td class="warning">Depends</td> </tr> </tbody> </table> Decomposes time series into components, allowing seasonality to be forecasted separately with simpler methods. ✅ Allows complex time series to be forecasted with simpler models. ❌ Highly dependent on the ability to decompose the time series. --- ## STL + ??? Uses repeated STL decompositions to decompose multiple seasonalities into separate time series. Then other models (`???`) can be used to forecast the simpler components. `$$y_{t} = S_{\text{year},t} + S_{\text{week},t} + S_{\text{day},t} + T_{t} + R_t$$` Available with direct fable compatibility in the [`{feasts}`](http://feasts.tidyverts.org/) package. --- ## Decomposing a time series Using an `STL` decomposition with default settings: ```r library(feasts) hospital_dcmp <- hospital %>% model(STL(arrivals)) components(hospital_dcmp) %>% autoplot() ``` <img src="figure/dcmp-stl-1.png" style="display: block; margin: auto;" /> --- ## Decomposing a time series Using an `STL` decomposition with default settings: ```r components(hospital_dcmp) %>% tail(24*7*5) %>% autoplot() ``` <img src="figure/dcmp-stl-plot-1.svg" style="display: block; margin: auto;" /> --- ## Decomposing a time series Refining the loess windows for the sub-daily data frequency: ```r hospital_dcmp <- hospital %>% model(STL(arrivals ~ trend(window = 24*7*5) + season("year", window = Inf) + season("week", window = 24*7*5) + season("day", window = 24*7))) components(hospital_dcmp) %>% autoplot() ``` <img src="figure/dcmp-custom-1.png" style="display: block; margin: auto;" /> --- ## Decomposing a time series Refining the loess windows for the sub-daily data frequency: ```r components(hospital_dcmp) %>% tail(24*7*5) %>% autoplot() ``` <img src="figure/dcmp-custom-plot-1.svg" style="display: block; margin: auto;" /> --- ## Modelling decomposed components As seasonality changes slowly, it is easily modelled with naive lags. ```r components(hospital_dcmp)[-1] %>% model(SNAIVE(season_week ~ lag("week"))) %>% forecast(h = "2 weeks") %>% autoplot(tail(components(hospital_dcmp)[-1], 24*7*4)) ``` <img src="figure/season-week-fc-1.svg" style="display: block; margin: auto;" /> --- ## Modelling decomposed components As seasonality changes slowly, it is easily modelled with naive lags. ```r components(hospital_dcmp)[-1] %>% model(SNAIVE(season_day ~ lag("day"))) %>% forecast(h = "2 weeks") %>% autoplot(tail(components(hospital_dcmp)[-1], 24*7*4)) ``` <img src="figure/season-day-fc-1.svg" style="display: block; margin: auto;" /> --- ## Modelling decomposed components The seasonally adjusted data (trend + remainder) is typically forecasted together. ```r components(hospital_dcmp)[-1] %>% model(ARIMA(season_adjust ~ 0 + pdq(3,0,3) + PDQ(1,0,0))) %>% forecast(h = "3 days") %>% autoplot(tail(components(hospital_dcmp)[-1], 24*4)) ``` <img src="figure/seasadjust-arima-1.svg" style="display: block; margin: auto;" /> --- ## Combining individual forecasts Use `decomposition_model()` to specify all component models: ```r dcmp_spec <- decomposition_model( dcmp = STL(arrivals ~ trend(window = 24*7*5) + season("year", window = Inf) + season("week", window = 24*7*5) + season("day", window = 24*7)), ARIMA(season_adjust ~ 0 + pdq(3,0,3) + PDQ(1,0,0)), SNAIVE(season_week ~ lag("week")), SNAIVE(season_day ~ lag("day")) ) ``` -- Notice that I haven't provided a model for `season_year`. A `SNAIVE()` model will be used by default for seasonal components. --- ## Combining individual forecasts Then use it as a regular model, the recombination is done automatically. ```r hospital %>% model(dcmp_spec) %>% forecast(h = "2 weeks") %>% autoplot(tail(hospital, 24*7*4)) ``` <img src="figure/dcmp-model-fc-1.svg" style="display: block; margin: auto;" /> --- class: inverse, top background-image: linear-gradient(to right, rgba(150, 150, 150, .1), rgba(150, 150, 150, .4)), url("resources/hourglass.jpg") background-size: cover .sticker-float[] .title[Thanks! <svg style="height:0.8em;top:.04em;position:relative;fill:white;" viewBox="0 0 576 512"><path d="M416 192c0-88.4-93.1-160-208-160S0 103.6 0 192c0 34.3 14.1 65.9 38 92-13.4 30.2-35.5 54.2-35.8 54.5-2.2 2.3-2.8 5.7-1.5 8.7S4.8 352 8 352c36.6 0 66.9-12.3 88.7-25 32.2 15.7 70.3 25 111.3 25 114.9 0 208-71.6 208-160zm122 220c23.9-26 38-57.7 38-92 0-66.9-53.5-124.2-129.3-148.1.9 6.6 1.3 13.3 1.3 20.1 0 105.9-107.7 192-240 192-10.8 0-21.3-.8-31.7-1.9C207.8 439.6 281.8 480 368 480c41 0 79.1-9.2 111.3-25 21.8 12.7 52.1 25 88.7 25 3.2 0 6.1-1.9 7.3-4.8 1.3-2.9.7-6.3-1.5-8.7-.3-.3-22.4-24.2-35.8-54.5z"/></svg>] .larger[ <svg style="height:0.8em;top:.04em;position:relative;fill:white;" viewBox="0 0 496 512"><path d="M336.5 160C322 70.7 287.8 8 248 8s-74 62.7-88.5 152h177zM152 256c0 22.2 1.2 43.5 3.3 64h185.3c2.1-20.5 3.3-41.8 3.3-64s-1.2-43.5-3.3-64H155.3c-2.1 20.5-3.3 41.8-3.3 64zm324.7-96c-28.6-67.9-86.5-120.4-158-141.6 24.4 33.8 41.2 84.7 50 141.6h108zM177.2 18.4C105.8 39.6 47.8 92.1 19.3 160h108c8.7-56.9 25.5-107.8 49.9-141.6zM487.4 192H372.7c2.1 21 3.3 42.5 3.3 64s-1.2 43-3.3 64h114.6c5.5-20.5 8.6-41.8 8.6-64s-3.1-43.5-8.5-64zM120 256c0-21.5 1.2-43 3.3-64H8.6C3.2 212.5 0 233.8 0 256s3.2 43.5 8.6 64h114.6c-2-21-3.2-42.5-3.2-64zm39.5 96c14.5 89.3 48.7 152 88.5 152s74-62.7 88.5-152h-177zm159.3 141.6c71.4-21.2 129.4-73.7 158-141.6h-108c-8.8 56.9-25.6 107.8-50 141.6zM19.3 352c28.6 67.9 86.5 120.4 158 141.6-24.4-33.8-41.2-84.7-50-141.6h-108z"/></svg> Learn more: [fable.tidyverts.org](https://fable.tidyverts.org/) <svg style="height:0.8em;top:.04em;position:relative;fill:white;" viewBox="0 0 512 512"><path d="M496 384H64V80c0-8.84-7.16-16-16-16H16C7.16 64 0 71.16 0 80v336c0 17.67 14.33 32 32 32h464c8.84 0 16-7.16 16-16v-32c0-8.84-7.16-16-16-16zM464 96H345.94c-21.38 0-32.09 25.85-16.97 40.97l32.4 32.4L288 242.75l-73.37-73.37c-12.5-12.5-32.76-12.5-45.25 0l-68.69 68.69c-6.25 6.25-6.25 16.38 0 22.63l22.62 22.62c6.25 6.25 16.38 6.25 22.63 0L192 237.25l73.37 73.37c12.5 12.5 32.76 12.5 45.25 0l96-96 32.4 32.4c15.12 15.12 40.97 4.41 40.97-16.97V112c.01-8.84-7.15-16-15.99-16z"/></svg> Keep updated: [tidyverts.org](http://www.tidyverts.org) <svg style="height:0.8em;top:.04em;position:relative;fill:white;" viewBox="0 0 576 512"><path d="M528 0H48C21.5 0 0 21.5 0 48v320c0 26.5 21.5 48 48 48h192l-16 48h-72c-13.3 0-24 10.7-24 24s10.7 24 24 24h272c13.3 0 24-10.7 24-24s-10.7-24-24-24h-72l-16-48h192c26.5 0 48-21.5 48-48V48c0-26.5-21.5-48-48-48zm-16 352H64V64h448v288z"/></svg> Review slides: [slides.mitchelloharawild.com/nhs2020](https://slides.mitchelloharawild.com/nhs2020) <svg style="height:0.8em;top:.04em;position:relative;fill:white;" viewBox="0 0 512 512"><path d="M459.37 151.716c.325 4.548.325 9.097.325 13.645 0 138.72-105.583 298.558-298.558 298.558-59.452 0-114.68-17.219-161.137-47.106 8.447.974 16.568 1.299 25.34 1.299 49.055 0 94.213-16.568 130.274-44.832-46.132-.975-84.792-31.188-98.112-72.772 6.498.974 12.995 1.624 19.818 1.624 9.421 0 18.843-1.3 27.614-3.573-48.081-9.747-84.143-51.98-84.143-102.985v-1.299c13.969 7.797 30.214 12.67 47.431 13.319-28.264-18.843-46.781-51.005-46.781-87.391 0-19.492 5.197-37.36 14.294-52.954 51.655 63.675 129.3 105.258 216.365 109.807-1.624-7.797-2.599-15.918-2.599-24.04 0-57.828 46.782-104.934 104.934-104.934 30.213 0 57.502 12.67 76.67 33.137 23.715-4.548 46.456-13.32 66.599-25.34-7.798 24.366-24.366 44.833-46.132 57.827 21.117-2.273 41.584-8.122 60.426-16.243-14.292 20.791-32.161 39.308-52.628 54.253z"/></svg> Say hello: [@mitchoharawild](twitter.com/mitchoharawild/) <br> .bottom[This work is licensed as <svg style="height:0.8em;top:.04em;position:relative;fill:white;" viewBox="0 0 496 512"><path d="M245.83 214.87l-33.22 17.28c-9.43-19.58-25.24-19.93-27.46-19.93-22.13 0-33.22 14.61-33.22 43.84 0 23.57 9.21 43.84 33.22 43.84 14.47 0 24.65-7.09 30.57-21.26l30.55 15.5c-6.17 11.51-25.69 38.98-65.1 38.98-22.6 0-73.96-10.32-73.96-77.05 0-58.69 43-77.06 72.63-77.06 30.72-.01 52.7 11.95 65.99 35.86zm143.05 0l-32.78 17.28c-9.5-19.77-25.72-19.93-27.9-19.93-22.14 0-33.22 14.61-33.22 43.84 0 23.55 9.23 43.84 33.22 43.84 14.45 0 24.65-7.09 30.54-21.26l31 15.5c-2.1 3.75-21.39 38.98-65.09 38.98-22.69 0-73.96-9.87-73.96-77.05 0-58.67 42.97-77.06 72.63-77.06 30.71-.01 52.58 11.95 65.56 35.86zM247.56 8.05C104.74 8.05 0 123.11 0 256.05c0 138.49 113.6 248 247.56 248 129.93 0 248.44-100.87 248.44-248 0-137.87-106.62-248-248.44-248zm.87 450.81c-112.54 0-203.7-93.04-203.7-202.81 0-105.42 85.43-203.27 203.72-203.27 112.53 0 202.82 89.46 202.82 203.26-.01 121.69-99.68 202.82-202.84 202.82z"/></svg> BY-NC 4.0.] ]