```{r, include = FALSE, echo=FALSE}
source("../setup.R")
```

## Outline

- What is temporal data?
- What is exploratory temporal data analysis?
- Using temporal objects in R: `tsibble`
- Data wrangling: aggregation, creating temporal components, missing values
- Plotting conventions: connect the dots; aspect ratio, landscape or portrait
- Calendar plots: arranging daily records into a calendar format
- Visual inference for temporal data
- tignostics: cognostics for temporal data
- Interactive graphics for temporal data
- Exploring longitudinal data, with the `brolgar` package


## Philosophy

:::: {.columns}
::: {.column}

> Time series analysis is what you do after all the interesting stuff has been done!

[Heike Hofmann, 2005](https://en.wikipedia.org/wiki/Heike_Hofmann)

<img src="../images/heike-headshot.png" style="width: 400px; border-radius: 50%">

:::
::: {.column}

::: {.fragment}
Time series analysis focuses on modeling the temporal dependence. Data needs to have trend, seasonality, anomalies removed first. 

::: {.info}
Exploratory temporal analysis involves exploring and discovering temporal trend, patterns related to seasons, and anomalies. And possibly also unusual temporal dependence.
:::

:::

:::
::::

## What is temporal data?

:::: {.columns}
::: {.column}

```{r}
#| eval: false
#| echo: false
CO2.ptb <- read.table("https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/merged_in_situ_and_flask/monthly/monthly_merge_co2_ptb.csv", sep=",", skip=69)
colnames(CO2.ptb) <- c("year", "month", "dateE", "date", "co2_ppm", "sa_co2", "fit", "sa_fit", "co2f", "sa_co2f")
CO2.ptb$lat <- (-71.3)
CO2.ptb$lon <- (-156.6)
CO2.ptb$stn <- "ptb"
CO2.ptb$co2_ppm <- replace_na(CO2.ptb$co2_ppm, -99.99)
  
save(CO2.ptb, file=here::here("data/CO2_ptb.rda"))
```

```{r}
#| label: CO2
#| echo: false
#| fig-width: 4
#| fig-height: 7
#| out-width: 50%
load(here::here("data/CO2_ptb.rda"))
CO2.ptb <- CO2.ptb |>
  filter(year > 2015) |>
  filter(co2_ppm > 100) # handle missing values
p1 <- ggplot(CO2.ptb, aes(x=date, y=co2_ppm)) + 
  geom_line(size=2, colour="#D93F00") + xlab("") + ylab("CO2 (ppm)")
p2 <- ggplot(CO2.ptb, aes(x=date, y=co2_ppm)) + 
  geom_smooth(se=FALSE, colour="#D93F00", size=2) + 
  xlab("") + ylab("CO2 (ppm)")
p1 + p2 + plot_layout(ncol=1)
```

:::
::: {.column style="font-size: 80%;"}

- Temporal data has  date/time/ordering index variable, call it [time]{.monash-blue2}. 
- A time variable has special structure:
    - it can have *cyclical* patterns, eg seasonality (summer, winter), an over in cricket
    - the cyclical patterns can be *nested*, eg postcode within state, over within innings
- Measurements are also [NOT independent]{.monash-blue2} - yesterday may influence today.
- It still likely has [non-cyclical patterns]{.monash-blue2}, trends and associations with other variables, eg temperature increasing over time, over is bowled by Elise Perry or Sophie Molineaux

:::
::::

## `tsibble`: R temporal object {background-image="https://tsibble.tidyverts.org/reference/figures/logo.png" background-position="5% 15%" background-size="10%"}

<br><br><br><br>
The `tsibble` package provides a data infrastructure for tidy temporal data with wrangling tools. Adapting the tidy data principles, `tsibble` is a data- and model-oriented object. In `tsibble`:

- [Index]{.monash-blue2} is a variable with inherent ordering from past to present.
- [Key]{.monash-blue2} is a set of variables that define observational units over time.
- Each observation should be [uniquely identified]{.monash-blue2} by index and key.
- Each observational unit should be measured at a [common interval]{.monash-blue2}, if regularly spaced.

## Regular vs irregular

:::: {.columns}
::: {.column}

The Melbourne pedestrian sensor data has a [regular]{.monash-blue2} period. Counts are provided for every hour, at numerous locations.

```{r}
#| label: ped-reg
options(width=55)
pedestrian 
```

:::
::: {.column}

::: {.fragment}
In contrast, the US flights data, below, is [irregular]{.monash-blue2}. 

```{r}
#| label: nycflights
options(width=55)
library(nycflights13)
flights_ts <- flights |>
  mutate(dt = ymd_hm(paste(paste(year, month, day, sep="-"), 
                           paste(hour, minute, sep=":")))) |>
  as_tsibble(index = dt, key = c(origin, dest, carrier, tailnum), regular = FALSE)
flights_ts 
```

:::
:::
::::

::: {.fragment}
[Is pedestrian traffic really regular?]{.monash-orange2}
:::

## Getting started

::: {.info}
Wrangling prior to analysing temporal data includes:

- **aggregate** by temporal unit. 
- **construct temporal units** to study seasonality, such as month, week, day of the week, quarter, ...
- checking and imputing **missings**.
:::

<br>
For the airlines data, you can aggregate by multiple quantities, eg number of arrivals, departures, average hourly arrival delay and departure delays. 

## Aggregating

:::: {.columns}
::: {.column}

The US flights data already has some temporal components created, so aggregating by these is easy. Here is departure delay.

```{r}
#| code-fold: true
flights_mth <- flights_ts |> 
  as_tibble() |>
  group_by(month, origin) |>
  summarise(dep_delay = mean(dep_delay, na.rm=TRUE)) |>
  as_tsibble(key=origin, index=month)
ggplot(flights_mth, aes(x=month, y=dep_delay, colour=origin)) +
  geom_point() +
  geom_smooth(se=F) +
  scale_x_continuous("", breaks = seq(1, 12, 1), 
                     labels=c("J","F","M","A","M","J",
                              "J","A","S","O","N","D")) +
  scale_y_continuous("av dep delay (mins)", limits=c(0, 25)) +
  theme(aspect.ratio = 0.5)
```

:::
::: {.column}

Aggregate by month, but examine arrival delays.

```{r}
#| code-fold: true
flights_mth_arr <- flights_ts |> 
  as_tibble() |>
  group_by(month, origin) |>
  summarise(arr_delay = mean(arr_delay, na.rm=TRUE)) |>
  as_tsibble(key=origin, index=month)
ggplot(flights_mth_arr, aes(x=month, y=arr_delay, colour=origin)) +
  geom_point() +
  geom_smooth(se=F) +
  scale_x_continuous("", breaks = seq(1, 12, 1), 
                     labels=c("J","F","M","A","M","J",
                              "J","A","S","O","N","D")) +
  scale_y_continuous("av arr delay (mins)", limits=c(0, 25)) +
  theme(aspect.ratio = 0.5)
```

:::
::::

## Constructing temporal units

:::: {.columns}
::: {.column style="font-size: 80%;"}

Week day vs weekend would be expected to have different patterns of delay, but this is not provided. 

```{r}
#| code-fold: true
#| fig-height: 8
#| fig-width: 5
#| out-width: 48%
flights_wk <- flights_ts |> 
  as_tibble() |>
  mutate(wday = wday(dt, label=TRUE, week_start = 1)) |>
  group_by(wday, origin) |>
  summarise(dep_delay = mean(dep_delay, na.rm=TRUE)) |>
  mutate(weekend = ifelse(wday %in% c("Sat", "Sun"), "yes", "no")) |>
  as_tsibble(key=origin, index=wday)
ggplot(flights_wk, aes(x=wday, y=dep_delay, fill=weekend)) +
  geom_col() +
  facet_wrap(~origin, ncol=1, scales="free_y") +
  xlab("") +
  ylab("av dep delay (mins)") +
  theme(aspect.ratio = 0.5, legend.position = "none")
```

:::
::: {.column style="font-size: 80%;"}

::: {.fragment}
Be careful of times!

```{r}
#| code-fold: true
flights_airtm <- flights |>
  mutate(dep_min = dep_time %% 100,
         dep_hr = dep_time %/% 100,
         arr_min = arr_time %% 100,
         arr_hr = arr_time %/% 100) |>
  mutate(dep_dt = ymd_hm(paste(paste(year, month, day, sep="-"), 
                           paste(dep_hr, dep_min, sep=":")))) |>
  mutate(arr_dt = ymd_hm(paste(paste(year, month, day, sep="-"), 
                           paste(arr_hr, arr_min, sep=":")))) |>
  mutate(air_time2 = as.numeric(difftime(arr_dt, dep_dt)))

fp <- flights_airtm |> 
  sample_n(3000) |>
  ggplot(aes(x=air_time, y=air_time2, label = paste(origin, dest))) + 
    geom_abline(intercept=0, slope=1) +
    geom_point()
ggplotly(fp, width=500, height=500)
```

Why is this not what we expect?
:::

:::
::::

## Checking and filling missings [(1/4)]{.smallest}

:::: {.columns}
::: {.column}

```{r}
#| label: missings-simple
set.seed(328)
harvest <- tsibble(
  year = c(2010, 2011, 2013, 2011, 
           2012, 2013),
  fruit = rep(c("kiwi", "cherry"), 
              each = 3),
  kilo = sample(1:10, size = 6),
  key = fruit, index = year
)
harvest
```

:::
::: {.column}

```{r}
#| label: missing-gaps
has_gaps(harvest, .full = TRUE) 
```

<br>
[Can you see the gaps in time?]{.monash-orange2}

Both levels of the key have missings.

:::
::::

## Checking and filling missings [(2/4)]{.smallest}

:::: {.columns}
::: {.column}

```{r}
#| label: missings-simple
#| echo: false
```

:::
::: {.column}

```{r}
#| label: count-gaps
count_gaps(harvest,  .full=TRUE)
```

<br>
One missing in each level, although it is a different year.

<br> <br>

Notice how `tsibble` handles this summary so neatly.

:::
::::

## Checking and filling missings [(3/4)]{.smallest}

:::: {.columns}
::: {.column}

```{r}
#| label: missings-simple
#| echo: false
```

:::
::: {.column}

Make the implicit missing values [explicit]{.monash-blue2}.

```{r}
#| label: fill-gaps
harvest <- fill_gaps(harvest, 
                     .full=TRUE) 
harvest 
```

:::
::::

## Checking and filling missings [(4/4)]{.smallest}

:::: {.columns}
::: {.column}

```{r}
#| label: missings-simple
#| echo: false
```

:::
::: {.column}

We have already seen `na_ma()` function, that imputes using a moving average. Alternatively, `na_interpolation()` uses the previous and next values to impute.

```{r}
#| label: impute-gaps
harvest_nomiss <- harvest |> 
  group_by(fruit) |> 
  mutate(kilo = 
    na_interpolation(kilo)) |> 
  ungroup()
harvest_nomiss 
```

:::
::::


## Plotting conventions {.transition-slide .center style="text-align: center;"}

## Conventions

- **lines**: connecting sequential time points reminding the reader that the temporal dependence is important.
- **aspect ratio**: wide or tall? [Cleveland, McGill, McGill (1988) ](https://eagereyes.org/basics/banking-45-degrees) argue the average line slope in a line chart should be 45 degrees, which is called banking to 45 degrees. But this is refuted in Talbot, Gerth, Hanrahan (2012) that the conclusion was based on a flawed study. Nevertheless, aspect ratio is an inescapable skill for designing effective plots. For time series, typically a wide aspect ratio is good. 
- **conventions**: 
    - time on the [horizontal]{.monash-blue2} axis, 
    - [ordering of elements]{.monash-blue2} like week day, month. Most software organises by alphabetical order, so this needs to be controlled.

## Aspect ratio

::: {.panel-tabset}

## 📊

```{r CO2_ratio, fig.width=12, fig.height=7, out.width="80%"}
#| label: CO2-ratio
#| fig-width: 12
#| fig-height: 7
#| out-width: 70%
#| echo: false
load(here::here("data/CO2_ptb.rda"))
CO2.ptb <- CO2.ptb |> 
  filter(year > 1980) |>
  filter(co2_ppm > 100) # handle missing values
p <- ggplot(CO2.ptb, aes(x=date, y=co2_ppm)) + 
  geom_line(size=1) + xlab("") + ylab("CO2 (ppm)")
p1 <- p + theme(aspect.ratio = 1) + ggtitle("1 to 1 (may be useless)")
p3 <- p + theme(aspect.ratio = 2) + ggtitle("tall & skinny:  trend")
p2 <- ggplot(CO2.ptb, aes(x=date, y=co2_ppm)) + 
  annotate("text", x=2000, y=375, label="CO2 at \n Point Barrow,\n Alaska", size=8) + theme_solid()
p4 <- p + 
  scale_x_continuous("", breaks = seq(1980, 2020, 5)) + 
  theme(aspect.ratio = 0.2) + ggtitle("short & wide: seasonality")
grid.arrange(p1, p2, p3, p4, layout_matrix= matrix(c(1,2,3,4,4,4), nrow=2, byrow=T))
```

## learn

- Is the trend linear or non-linear? 
    - Yes, slightly non-linear. We could fit a linear regression model, and examine the residuals to better assess non-linear trend.
- Is there a cyclical pattern?
    - Yes, there is a yearly trend. 
<br>
<br>
<br>

*This type of data is easy to model, and forecast.*

## R

```{r}
#| label: CO2-ratio
#| echo: true
#| eval: false
```

:::

## Calendar plot {.transition-slide .center style="text-align: center;"}

## Case study: NYC flights [(1/2)]{.smallest}

::: {.panel-tabset}

## 📊


```{r}
#| label: calendar
#| fig-width: 10
#| fig-height: 6
#| out-width: 70%
#| echo: false
flights_hourly <- flights |>
  group_by(time_hour, origin) |> 
  summarise(count = n(), 
    dep_delay = mean(dep_delay, 
                     na.rm = TRUE)) |> 
  ungroup() |>
  as_tsibble(index = time_hour, 
             key = origin) |>
    mutate(dep_delay = 
    na_interpolation(dep_delay)) 
calendar_df <- flights_hourly |> 
  filter(origin == "JFK") |>
  mutate(hour = hour(time_hour), 
         date = as.Date(time_hour)) |>
  filter(year(date) < 2014) |>
  frame_calendar(x=hour, y=count, date=date, nrow=2) 
p1 <- calendar_df |>
  ggplot(aes(x = .hour, y = .count, group = date)) +
  geom_line() + theme(axis.line.x = element_blank(),
                      axis.line.y = element_blank()) +
  theme(aspect.ratio=0.5)
prettify(p1, size = 3, label.padding = unit(0.15, "lines"))
```

## About calendars

- Draws the daily data in the layout of a regular calendar
- A wonderful way to get a lot of data into a page
- Easy to examine daily patterns, weekly, monthly patterns

## What do we see?

- The daily pattern at JFK is **very** regular. 
- It is similar for every day of the week, and for every month
- There is a peak in early flights, a drop around lunchtime and then the number of flights pick up again.

Something is fishy here. What is it?

## R

```{r}
#| label: calendar
#| echo: true
#| eval: false
```

:::

## Case study: NYC flights [(2/2)]{.smallest}

::: {.panel-tabset}

## 📊


```{r}
#| label: calendar-delay
#| fig-width: 10
#| fig-height: 6
#| out-width: 70%
#| echo: false
calendar_df <- flights_hourly |> 
  filter(origin == "JFK") |>
  mutate(hour = hour(time_hour), 
         date = as.Date(time_hour)) |>
  filter(year(date) < 2014) |>
  frame_calendar(x=hour, y=dep_delay, date=date, nrow=2) 
p1 <- calendar_df |>
  ggplot(aes(x = .hour, y = .dep_delay, group = date)) +
  geom_line() + theme(axis.line.x = element_blank(),
                      axis.line.y = element_blank()) +
  theme(aspect.ratio=0.5)
prettify(p1, size = 3, label.padding = unit(0.15, "lines"))
```

## What do we see?

Delays are much more interesting to examine

- Most days have few delays
- Jun and July seem to have more delays
- A few days, sporadically in the year, have big delays

[Can you find a reason for one of the days with a big delay?]{.monash-orange2}

::: {style="font-size: 60%;"}
From ChatGPT:
*As of my last update in September 2021, a significant late-season snowstorm did affect parts of the United States in April 2013, but it was more focused on the Midwest rather than the Northeast where JFK Airport (John F. Kennedy International Airport) is located. The storm impacted states like Minnesota, Wisconsin, and South Dakota, among others, and brought heavy snowfall and icy conditions.*

*However, weather conditions can have a cascading effect on flight schedules nationwide, so it's possible that there were some delays at JFK related to this or other weather phenomena.*  
:::

## R

```{r}
#| label: calendar-delay
#| echo: true
#| eval: false
```

:::

## Visual inference {.transition-slide .center style="text-align: center;"}

## Temporal patterns: simulation

:::: {.columns}
::: {.column width=70%}

```{r}
#| label: lm-lineup
#| code-fold: true
#| fig-width: 10
#| fig-height: 5
#| out-width: 100%
p_bourke <- pedestrian |>
  as_tibble() |>
  filter(Sensor == "Bourke Street Mall (North)",
         Date >= ymd("2015-05-03"), Date <= ymd("2015-05-16")) |>
  mutate(date_num = 
    as.numeric(difftime(Date_Time,ymd_hms("2015-05-03 00:00:00"),
       units="hours"))+11) |> # UTC to AEST
  mutate(day = wday(Date, label=TRUE, week_start=1)) |>
  select(date_num, Time, day, Count) |>
  rename(time = date_num, hour=Time, count = Count)
# Fit a linear model with categorical hour variable
p_bourke_lm <- glm(count~day+factor(hour), family="poisson", 
  data=p_bourke)
# Function to simulate from a Poisson
simulate_poisson <- function(model, newdata) {
  lambda_pred <- predict(model, newdata, type = "response")
  rpois(length(lambda_pred), lambda = lambda_pred)
}

set.seed(436)
pos <- sample(1:12)
p_bourke_lineup <- bind_cols(.sample = rep(pos[1], 
  nrow(p_bourke)), p_bourke[,-2])
for (i in 1:11) {
  new <- simulate_poisson(p_bourke_lm, p_bourke)
  x <- tibble(time=p_bourke$time, count=new)
  x <- bind_cols(.sample = rep(pos[i+1], 
         nrow(p_bourke)), x)
  p_bourke_lineup <- bind_rows(p_bourke_lineup, x)
}

ggplot(p_bourke_lineup,
  aes(x=time, y=count)) + 
  geom_line() +
  facet_wrap(~.sample, ncol=4) +
  theme(aspect.ratio=0.5, 
        axis.text = element_blank(),
        axis.title = element_blank())
```

:::
::: {.column width=30%}

1. Decide on a model
2. Simulate from the model to generate nulls

:::
::::

## Association: permutation

:::: {.columns}
::: {.column width=70%}

```{r}
#| label: NYC-lineup
#| code-fold: true
#| fig-width: 10
#| fig-height: 5
#| out-width: 100%
set.seed(514)
ggplot(lineup(null_permute("origin"), true=flights_mth, n=12), 
       aes(x=month, y=dep_delay, colour=origin)) +
  geom_point() +
  geom_smooth(se=F) +
  facet_wrap(~.sample, ncol=4) +
  theme(aspect.ratio = 0.5, 
        legend.position = "none",
        axis.text = element_blank(),
        axis.title = element_blank())
```
:::
::: {.column width=30%}

Break association between variables. Here `origin` is permuted which breaks association with `dep_delay`, while keeping `month` fixed.
<br><br>

*Which plot has the biggest difference between the three groups?*
:::
::::

## Tignostics {.transition-slide .center style="text-align: center;"}

## `feasts`: time series features {background-image="https://feasts.tidyverts.org/reference/figures/logo.png" background-position="5% 15%" background-size="10%"}

:::: {.columns}
::: {.column style="font-size: 80%;"}
<br><br><br><br><br><br>
The `feasts` package provides functions to calculate tignostics for time series.

*Remember scagnostics?*
<br>

Compute [tignostics]{.monash-blue2} for each series, for example,

- trend 
- seasonality
- linearity
- spikiness
- peak
- trough

:::

::: {.column style="font-size: 70%;"}

```{r}
#| label: ts-features
#| fig-width: 4
#| fig-height: 4
#| out-width: 80%
#| code-fold: true
tourism_feat <- tourism |>
  features(Trips, feat_stl)
tourism_feat |>
  ggplot(aes(x = trend_strength, y = seasonal_strength_year)) +
  geom_point()  
```

:::
::::

## Interactivity {.transition-slide .center style="text-align: center;"}


## Interactive exploration with `tsibbletalk`

```{r}
#| label: tsibbletalk1
#| code-fold: true
#| fig-height: 5
#| fig-width: 10
#| out-width: 100%
tourism_shared <- tourism |>
  as_shared_tsibble(spec = (State / Region) * Purpose)

tourism_feat <- tourism_shared |>
  features(Trips, feat_stl)

p1 <- tourism_shared |>
  ggplot(aes(x = Quarter, y = Trips)) +
  geom_line(aes(group = Region), alpha = 0.5) +
  facet_wrap(~ Purpose, scales = "free_y") 
p2 <- tourism_feat |>
  ggplot(aes(x = trend_strength, y = seasonal_strength_year)) +
  geom_point(aes(group = Region))
  
subplot(
    ggplotly(p1, tooltip = "Region", width = 1400, height = 700),
    ggplotly(p2, tooltip = "Region", width = 1200, height = 600),
    nrows = 1, widths=c(0.5, 0.5), heights=1) |>
  highlight(dynamic = FALSE)
  
```

## Wrapping series

:::: {.columns}

::: {.column}

Pedestrian counts at Bourke St Mall, has a daily seasonality.

**DEMO**

```{r}
#| label: tsibbletalk3
#| eval: false
pp <- p_bourke |>
        as_tsibble(index = time) |>
        ggplot(aes(x=time, y=count)) + 
          geom_line() +
          theme(aspect.ratio=0.5)
 
  
ui <- fluidPage(tsibbleWrapUI("tswrap"))
server <- function(input, output, session) {
  tsibbleWrapServer("tswrap", pp, period = "1 day")
}

shinyApp(ui, server)
```


:::
::: {.column}


**Famous data: Lynx**

Annual numbers of lynx trappings for 1821–1934 in Canada. [Almost]{.monash-blue2} 10 year cycle. *Explore periodicity by wrapping series on itself.*

**DEMO**

```{r}
#| label: tsibbletalk4
#| eval: false
lynx_tsb <- as_tsibble(lynx) |>
  rename(count = value)
pl <- ggplot(lynx_tsb, 
  aes(x = index, y = count)) +
  geom_line(size = .2) 

ui <- fluidPage(
  tsibbleWrapUI("tswrap"))
server <- function(input, output, 
                   session) {
  tsibbleWrapServer("tswrap", pl, 
       period = "1 year")
}
shinyApp(ui, server)
```


:::
::::

## Longitudinal data {.transition-slide .center style="text-align: center;"}

## Longitudinal vs time series


::: {.info}
Longitudinal data tracks the same sample of individuals at different points in time. It often has different lengths and different time points for each individual.
:::


```{r}
#| label: ts-visual
#| fig-width: 12
#| fig-height: 4
#| out-width: 100%
#| echo: false
pts <- pedestrian |>
  filter(Sensor == "Southern Cross Station") |>
  filter(between(Date, ymd("2015-07-06"), ymd("2015-07-13"))) |> ggplot() +
  geom_line(aes(x=Date_Time, y=Count)) +
  xlab("") +
  ggtitle("Time series") +
  theme(aspect.ratio=0.5)
plong <- wages |>
  sample_n_keys(size = 10) |>
  ggplot() +
  geom_line(aes(x=xp, y=ln_wages, group=id, colour=factor(id))) +
  xlab("Years") + ylab("Wages (log)") +
  ggtitle("Longitudinal") + 
  theme(aspect.ratio=0.5, legend.position="none") 
 pts + plong          
```

When the time points are the same for each individual, it is usually referred to as [panel data]{.monash-blue2}. When different individuals are measured at each time point, it is called [cross-sectional data]{.monash-blue2}.

## Overall trend

:::: {.columns}
::: {.column width=40%}

::: {style="font-size: 80%;"}

Log(wages) of 888 individuals, measured at various times in their employment [US National Longitudinal Survey of Youth]{.smallest}. 

```{r}
#| label: wages-trend1  
#| code-fold: true
#| fig-width: 6
#| fig-height: 4
#| out-width: 100%
wages |>
  ggplot() +
    geom_line(aes(x = xp, y = ln_wages, group = id), alpha=0.1) +
    geom_smooth(aes(x = xp, y = ln_wages), se=F) +
    xlab("years of experience") +
    ylab("wages (log)") +
  theme(aspect.ratio = 0.6)
```

Wages tend to increase as time in the workforce gets longer, on average.
:::

:::
::: {.column width=60%}

::: {.fragment style="font-size: 70%;"}

```{r}
#| label: wages-trend2  
#| code-fold: true
#| fig-width: 8
#| fig-height: 4
#| out-width: 100%
wages |>
  ggplot() +
    geom_line(aes(x = xp, y = ln_wages, group = id), alpha=0.1) +
    geom_smooth(aes(x = xp, y = ln_wages, 
      group = high_grade, colour = high_grade), se=F) +
    xlab("years of experience") +
    ylab("wages (log)") +
  scale_colour_viridis_c("education") +
  theme(aspect.ratio = 0.6)
```

The higher the education level achieved, the higher overall wage, on average.
:::
:::
::::

## Eating spaghetti


:::: {.columns}

::: {.column}

`brolgar` uses `tsibble` as the data object, and provides:

- sampling individuals
- longnostics for individuals
- diagnostics for statistical models

:::

::: {.column style="font-size: 70%;"}

::: {.panel-tabset}

## Sample 1

```{r}
#| label: sample-n1
#| code-fold: true
set.seed(753)
wages |>
  sample_n_keys(size = 10) |> 
  ggplot(aes(x = xp,
             y = ln_wages,
             group = id,
             colour = as.factor(id))) + 
  geom_line() +
  xlim(c(0,13)) + ylim(c(0, 4.5)) +
  xlab("years of experience") +
  ylab("wages (log)") +
  theme(aspect.ratio = 0.6, legend.position = "none")
```
 
## Sample 2

```{r}
#| label: sample-n2
#| code-fold: true
set.seed(749)
wages |>
  sample_n_keys(size = 10) |> 
  ggplot(aes(x = xp,
             y = ln_wages,
             group = id,
             colour = as.factor(id))) + 
  geom_line() +
  xlim(c(0,13)) + ylim(c(0, 4.5)) +
  xlab("years of experience") +
  ylab("wages (log)") +
  theme(aspect.ratio = 0.6, legend.position = "none")
```

## Sample 3

```{r}
#| label: sample-n3
#| code-fold: true
set.seed(757)
wages |>
  sample_n_keys(size = 10) |> 
  ggplot(aes(x = xp,
             y = ln_wages,
             group = id,
             colour = as.factor(id))) + 
  geom_line() +
  xlim(c(0,13)) + ylim(c(0, 4.5)) +
  xlab("years of experience") +
  ylab("wages (log)") +
  theme(aspect.ratio = 0.6, legend.position = "none")
```

:::

*Few individuals experience wages like the overall trend.*

:::  
::::

## Individual patterns

:::: {.columns}

::: {.column}
*Remember scagnostics?*
<br>

Compute [longnostics]{.monash-blue2} for each subject, for example,

- Slope, intercept from simple linear model
- Variance, standard deviation
- Jumps, differences
:::

::: {.column style="font-size: 70%;"}

::: {.panel-tabset}

## Increasing

```{r}
#| label: increasing
#| code-fold: true
wages_slope <- wages |>   
  add_n_obs() |>
  filter(n_obs > 4) |>
  add_key_slope(ln_wages ~ xp) |> 
  as_tsibble(key = id, index = xp) 
wages_spread <- wages |>
  features(ln_wages, feat_spread) |>
  right_join(wages_slope, by="id")

wages_slope |> 
  filter(.slope_xp > 0.3) |> 
  ggplot(aes(x = xp, 
             y = ln_wages, 
             group = id,
             colour = factor(id))) + 
  geom_line() +
  xlim(c(0, 4.5)) +
  ylim(c(0, 4.5)) +
  xlab("years of experience") +
  ylab("wages (log)") +
  theme(aspect.ratio = 0.6, legend.position = "none")
```

## Decreasing

```{r}
#| label: decreasing
#| code-fold: true
wages_slope |> 
  filter(.slope_xp < (-0.4)) |> 
  ggplot(aes(x = xp, 
             y = ln_wages, 
             group = id,
             colour = factor(id))) + 
  geom_line() +
  xlim(c(0, 4.5)) +
  ylim(c(0, 4.5)) +
  xlab("years of experience") +
  ylab("wages (log)") +
  theme(aspect.ratio = 0.6, legend.position = "none")
```

## Consistent

```{r}
#| label: small-sigma
#| code-fold: true
wages_spread |> 
  filter(sd < 0.1) |> 
  ggplot(aes(x = xp, 
             y = ln_wages, 
             group = id,
             colour = factor(id))) + 
  geom_line() +
  xlim(c(0, 12)) +
  ylim(c(0, 4.5)) +
  xlab("years of experience") +
  ylab("wages (log)") +
  theme(aspect.ratio = 0.6, legend.position = "none")
```

## Volatile

```{r}
#| label: large-sigma
#| code-fold: true
wages_spread |> 
  filter(sd > 0.8) |> 
  ggplot(aes(x = xp, 
             y = ln_wages, 
             group = id,
             colour = factor(id))) + 
  geom_line() +
  xlim(c(0, 12)) +
  ylim(c(0, 4.5)) +
  xlab("years of experience") +
  ylab("wages (log)") +
  theme(aspect.ratio = 0.6, legend.position = "none")
```

:::
:::
::::

## Individual summaries

:::: {.columns}

::: {.column style="font-size: 80%;"}

A different style of five number summary: What does average look like?
What do extremes look like?

Find those individuals who are [representative]{.monash-blue2} of the min, median, maximum, etc of a particular feature, e.g. trend, using `keys_near()`. This reports the individual who is closest to a particular statistic.

`wages_threenum()` returns the three individuals: min, max and closest to the median value, for a particular feature. 

`wages_fivenum()` returns the five individuals: min, max and closest to the median, Q1 and Q3 values, for a particular feature.
:::
::: {.column style="font-size: 70%;"}

```{r}
#| label: five-number
#| code-fold: true
#| fig-width: 8
#| fig-height: 5
wages_fivenum <- wages |>   
  add_n_obs() |>
  filter(n_obs > 6) |>
  key_slope(ln_wages ~ xp) |>
  keys_near(key = id,
            var = .slope_xp,
            funs = l_five_num) |> 
  left_join(wages, by = "id") |>
  as_tsibble(key = id, index = xp) 
  
wages_fivenum |>
  ggplot(aes(x = xp,
             y = ln_wages,
             group = id)) + 
  geom_line() + 
  ylim(c(0, 4.5)) +
  facet_wrap(~stat, ncol=3) +
  xlab("years of experience") +
  ylab("wages (log)") +
  theme(aspect.ratio = 0.6, legend.position = "none")
```

:::
::::

## Assessing model fits

:::: {.columns}

::: {.column style="font-size: 80%;"}

Fit a mixed effect model, education as fixed effect, subject random effect using slope. 

Summary of the fit

```{r}
#| label: model-fit
#| code-fold: true
#| fig-width: 8
#| fig-height: 4
wages_fit_int <- 
  lmer(ln_wages ~ xp + high_grade + 
         (xp |id), data = wages) 
wages_aug <- wages |>
  add_predictions(wages_fit_int, 
                  var = "pred_int") |>
  add_residuals(wages_fit_int, 
                var = "res_int")
  
m1 <- ggplot(wages_aug,
       aes(x = xp,
           y = pred_int,
           group = id)) + 
  geom_line(alpha = 0.2) +
  xlab("years of experience") +
  ylab("wages (log)") +
  theme(aspect.ratio = 0.6)
  
m2 <- ggplot(wages_aug,
       aes(x = pred_int,
           y = res_int,
           group = id)) + 
  geom_point(alpha = 0.5) +
  xlab("fitted values") + ylab("residuals")  

m1 + m2 + plot_layout(ncol=2) 
    
```

:::
::: {.column style="font-size: 80%;"}

Diagnosing the fit: each individual model

```{r}
#| label: model-diag
#| code-fold: true
#| fig-width: 8
#| fig-height: 7
#| out-width: 80%
wages_aug |> add_n_obs() |> filter(n_obs > 4) |>
  sample_n_keys(size = 12) |>
  ggplot() + 
  geom_line(aes(x = xp, y = pred_int, group = id, 
             colour = factor(id))) + 
  geom_point(aes(x = xp, y = ln_wages, 
                 colour = factor(id))) + 
  facet_wrap(~id, ncol=3)  +
  xlab("Years of experience") + ylab("Log wages") +
  theme(aspect.ratio = 0.6, legend.position = "none")
  
```

:::
::::


## Resources

- Wang, Cook, Hyndman (2019) [A New Tidy Data Structure to Support Exploration and Modeling of Temporal Data](https://www.tandfonline.com/doi/full/10.1080/10618600.2019.1695624)
- Wang, Cook, Hyndman, O'Hara-Wild (2019) [tsibble](https://tsibble.tidyverts.org)
- O'Hara-Wild, Hyndman, Wang (2020). [fabletools: Core Tools for Packages in the 'fable' Framework](https://CRAN.R-project.org/package=fabletools)
- O'Hara-Wild, Hyndman, Wang (2024). [feasts: Feature Extraction and Statistics for Time Series](https://feasts.tidyverts.org)
- Tierney, Cook, Prvan (2020) [Browse Over Longitudinal Data Graphically and Analytically in R](https://github.com/njtierney/brolgar)
