ANLY 505 - NYCFlightsDataModel
Week 3
Shaonan Hu
2019-11-25
Summarize Data
daily <- flights %>%
mutate(date = make_date(year, month, day)) %>%
group_by(date) %>%
summarize(n = n())
ggplot(daily, aes(date, n)) +
geom_line()
Investigate Daily-Weekly Pattern
daily <- daily %>%
mutate(wday = wday(date, label = TRUE))
ggplot(daily, aes(wday,n)) +
geom_boxplot()
mod = lm(n ~ wday, data = daily, na.action = na.warn)
grid <- daily %>%
data_grid(wday) %>%
add_predictions(mod, "n")
ggplot(daily, aes(wday, n)) +
geom_boxplot() +
geom_point(data = grid, color = "orange", size = 4)
Investigate residuals
daily <- daily %>%
add_residuals(mod)
daily %>%
ggplot(aes(date, resid)) +
geom_ref_line(h = 0) +
geom_line()
ggplot(daily, aes(date, resid, color = wday)) +
geom_ref_line(h = 0, colour = "red") +
geom_line()
daily %>%
filter(resid < -100)## # A tibble: 11 x 4
## date n wday resid
## <date> <int> <ord> <dbl>
## 1 2013-01-01 842 周二 -109.
## 2 2013-01-20 786 周日 -105.
## 3 2013-05-26 729 周日 -162.
## 4 2013-07-04 737 周四 -229.
## 5 2013-07-05 822 周五 -145.
## 6 2013-09-01 718 周日 -173.
## 7 2013-11-28 634 周四 -332.
## 8 2013-11-29 661 周五 -306.
## 9 2013-12-24 761 周二 -190.
## 10 2013-12-25 719 周三 -244.
## 11 2013-12-31 776 周二 -175.daily %>%
ggplot(aes(date, resid)) +
geom_ref_line(h = 0, colour = "red", size = 1) +
geom_line(color = "grey50") +
geom_smooth(se = FALSE, span = 0.20)## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Seasonal Saturday effect
daily %>%
filter(wday == "Sat") %>%
ggplot(aes(date, n)) +
geom_point()+
geom_line() +
scale_x_date(
NULL,
date_breaks = "1 month",
date_labels = "%b"
)
Add Seasonal Variable
term <- function(date) {
cut(date,
breaks = ymd(20130101, 20130605, 20130825, 20140101),
labels = c("spring", "summer", "fall")
)
}
daily <- daily %>%
mutate(term = term(date))
daily %>%
filter(wday == "Sat") %>%
ggplot(aes(date, n, color = term)) +
geom_point(alpha = 1/3)+
geom_line() +
scale_x_date(
NULL,
date_breaks = "1 month",
date_labels = "%b"
)
daily %>%
ggplot(aes(wday, n, color = term)) +
geom_boxplot()
mod1 <- lm(n ~ wday, data = daily, na.action = na.warn)
mod2 <- lm(n ~ wday * term, data = daily, na.action = na.warn)
daily %>%
gather_residuals(without_term = mod1, with_term = mod2) %>%
ggplot(aes(date, resid, color = model)) +
geom_line(alpha = 0.75)
grid <- daily %>%
data_grid(wday, term) %>%
add_predictions(mod2, "n")
ggplot(daily, aes(wday, n)) +
geom_boxplot() +
geom_point(data = grid, color = "red") +
facet_wrap(~ term)
Better model for outliers (Robust regression)
mod3 <- MASS::rlm(n ~ wday * term, data = daily, na.action = na.warn)
daily %>%
add_residuals(mod3, "resid") %>%
ggplot(aes(date, resid)) +
geom_hline(yintercept = 0, size = 2, color = "red") +
geom_line()
Computed Variables
# If you are creating variables it might be a good idea to bundle the creation of the variables up into a function
compute_vars <- function(data) {
data %>%
mutate(term = term(date),
wday = wday(date, label = TRUE)
)
}
# Another option would be to put the transformations directly in the model formula:
wday2 <- function(x) wday(x, label = TRUE)
mod3 <- lm(n ~ wday2(date) * term(date), data = daily, na.action = na.warn)Time of Year: An Alternative Approach
# We could use a more flexible model to capture the pattern of school term in the data
library(splines)
mod <- MASS::rlm(n ~ wday * ns(date, 5), data = daily, na.action = na.warn)
daily %>%
data_grid(wday, date = seq_range(date, n = 13)) %>%
add_predictions(mod) %>%
ggplot(aes(date, pred, color = wday)) +
geom_line() +
geom_point()
# We see a strong pattern in the numbers of Sat flights. This is reassuring, because we also saw that pattern in the raw data. It's a good sign when you get the same signal from different approaches.Question #1
Why are there fewer than expected flights on January 20, May 26 and September 1? (Hint: they all have the same explanation.) How would these days generalize into another year?
# Use this chunk to answer question 1Jan 21 is Martin Luther King Jr. Day, May 27 is Memorial day, and Sep 2 is labot day. Since these holiday are specific days based on week, they can’t be generalised to other years because those dates depend on weekdays.
Question #2
What do the three days with high positive residuals represent? How would these days generalize to another year?
# Use this chunk to answer question 2
daily %>%
top_n(3, resid)## # A tibble: 3 x 5
## date n wday resid term
## <date> <int> <ord> <dbl> <fct>
## 1 2013-11-30 857 周六 112. fall
## 2 2013-12-01 987 周日 95.5 fall
## 3 2013-12-28 814 周六 69.4 fallThese are the Saturdays close to Thanksgiving Day in the United States , Christmas Day or Hanukkah, and New Year’s Day .Traveling around popular holidays is in high demand. 2014 would witness the same case.
Question #3
Create a new variable that splits the “wday” variable into terms, but only for Saturdays, i.e., it should have Thurs, Fri, but Sat-summer, Sat-spring, Sat-fall. How does this model compare with the model with every combination of “wday” and “term”?
# Use this chunk to answer question 3
daily <- daily %>%
mutate(term = term(date)) %>%
mutate(term2 = ifelse(wday == 'Sat',paste0(wday,"-",term),as.character(term) ))
mod3 <- lm(n~ wday * term2, data = daily)
daily %>%
gather_residuals(mod3,mod2)%>%
arrange(date)%>%
ggplot(aes(date,resid,color = model))+
geom_line(alpha = 0.75)
Question #4
Create a new “wday” variable that combines the day of week, term(for Saturdays), and public holidays. What do the residuals of the model look like?
# Use this chunk to answer question 4
daily_holidays <-
daily %>%
mutate(holidays = case_when(date %in% ymd(c(20130101,
20130121,20130527,20130704, 20130902, 20131028, 20131111,20131128,20131225)) ~ "holiday",TRUE ~ "None")) %>%
unite(term, term2, holidays)
mod4 <- lm(n ~ term, data = daily_holidays)
daily_holidays %>%
add_residuals(mod4) %>%
ggplot(aes(date, resid)) +
geom_line()
Question #5
What happens if you fit a day-of-week effect that varies by month (i.e.m n ~ wday*month)? Why is this not very helpful?
# Use this chunk to answer question 5
daily <- flights %>%
mutate(date = make_date(year, month, day)) %>%
group_by(date,month) %>%
summarise(n = n())
daily <- daily %>%
mutate(wday = wday(date, label = TRUE)) %>%
mutate(term = term(date))
mod2 <- lm(n ~ wday * term, data = daily)
mod5 <- lm(n~ wday * month, data = daily)
daily %>%
gather_residuals(mod5,mod2)%>%
arrange(date)%>%
ggplot(aes(date,resid,color = model))+
geom_line(alpha = 0.75)
Question #6
What would you expect the model n ~ wday + ns(date,5) to look like? Knowing what you know about the data, why would you expect it not to be particularly effective?
# Use this chunk to answer question 6
mod5 <- lm(n ~ wday + splines::ns(date, 5), data = daily)
daily %>%
gather_residuals(mod2,mod5)%>%
arrange(date)%>%
ggplot(aes(date,resid,color = model))+
geom_line(alpha = 0.75)
Question #7
We hypothesized that people leaving on Sundays are more likely to be business travelers who need to be somewhere on Monday. Explore the hypothesis by seeing how if breaks down based on distance and time: if it’s true, you’d expect to see more Sunday evening flights to places that are far away.
# Use this chunk to answer question 7
week_relevel <- function(x) {
fct_relevel(x, "Sun", after = 7) # after 7 means night flights.
}
daily %>%
mutate(wday = week_relevel(wday)) %>%
ggplot(aes(wday, n)) + ggtitle("Sunday Night Flights") +
geom_boxplot()## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun
## Warning: Unknown levels in `f`: Sun