ETC5521: Diving Deeply into Data Exploration
Exploring data having a space and time context Part I
Philosophy
Time series analysis is what you do after all the interesting stuff has been done!
Time series analysis focuses on modeling the temporal dependence. Data needs to have trend, seasonality, anomalies removed first.
Exploratory temporal analysis involves exploring and discovering temporal trend, patterns related to seasons, and anomalies. And possibly also unusual temporal dependence.
Aggregating
The US flights data already has some temporal components created, so aggregating by these is easy. Here is departure delay.
Code
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)
Aggregate by month, but examine arrival delays.
Code
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
Week day vs weekend would be expected to have different patterns of delay, but this is not provided.
Code
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")
Be careful of times!
Code
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?
Aspect ratio
- 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.
- Yes, there is a yearly trend.
This type of data is easy to model, and forecast.
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))Case study: NYC flights (1/2)
- 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
- 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?
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"))Case study: NYC flights (2/2)
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?
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.
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"))Temporal patterns: simulation
Code
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())
- Decide on a model
- Simulate from the model to generate nulls
Association: permutation
Code
Break association between variables. Here origin is permuted which breaks association with dep_delay, while keeping month fixed.
Which plot has the biggest difference between the three groups?
feasts: time series features
The feasts package provides functions to calculate tignostics for time series.
Remember scagnostics?
Compute tignostics for each series, for example,
- trend
- seasonality
- linearity
- spikiness
- peak
- trough
Longitudinal vs time series
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.
When the time points are the same for each individual, it is usually referred to as panel data. When different individuals are measured at each time point, it is called cross-sectional data.
Overall trend
Log(wages) of 888 individuals, measured at various times in their employment US National Longitudinal Survey of Youth.
Code
Wages tend to increase as time in the workforce gets longer, on average.
Code
The higher the education level achieved, the higher overall wage, on average.
Eating spaghetti
brolgar uses tsibble as the data object, and provides:
- sampling individuals
- longnostics for individuals
- diagnostics for statistical models
Few individuals experience wages like the overall trend.
Individual patterns
Remember scagnostics?
Compute longnostics for each subject, for example,
- Slope, intercept from simple linear model
- Variance, standard deviation
- Jumps, differences
Code
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")
Individual summaries
A different style of five number summary: What does average look like? What do extremes look like?
Find those individuals who are representative 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.
Code
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
Fit a mixed effect model, education as fixed effect, subject random effect using slope.
Summary of the fit
Code
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) 
Diagnosing the fit: each individual model
Code
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")