co2 <- read_csv("../data/CO2-monthly.csv") |>
mutate(date = dmy(date))
co2_long <- co2 |>
pivot_longer(ljo:spo, names_to = "stn", values_to = "co2") |>
filter(stn %in% c("ptb", "ljo", "spo", "mlf")) |>
mutate(stn = factor(stn, levels = c("ptb", "ljo", "mlf", "spo"))) |>
filter(!is.na(co2)) |>
filter(co2 < 430)
co2_p1 <- ggplot(co2_long, aes(x=date, y=co2, colour=stn)) +
geom_point() +
facet_wrap(~stn, ncol=2) +
theme(legend.position="none") +
ylim(c(330, 430)) +
theme(aspect.ratio = 0.5)
co2_fit <- co2_long |>
mutate(time = as.numeric(date)-7304) |>
nest(dat = -stn) |>
mutate(mod_trend = map(dat, ~lm(co2~time+I(time^2), data=.x)),
mod_season = map(dat, ~lm(co2~time+I(time^2)+month, data=.x)),
coef_trend = map(mod_trend, tidy),
coef_seasonal = map(mod_season, tidy),
fit_trend = map(mod_trend, augment),
fit_season = map(mod_season, augment)) |>
select(stn, starts_with('fit')) |>
pivot_longer(
cols = starts_with('fit')
) |>
unnest(value)
co2_p2 <- co2_fit |>
filter(name == "fit_trend") |>
ggplot(aes(x=time, y=.fitted, colour=stn)) +
geom_point() +
ylim(c(330, 430))
co2_p3 <- co2_fit |>
filter(name == "fit_season") |>
mutate(year = time %/% 356 + 1990) |>
group_by(stn, year) |>
mutate(.fitted_seas = .fitted - mean(.fitted, na.rm=TRUE)) |>
ggplot(aes(x=time, y=.fitted_seas, colour=stn)) +
geom_line() +
facet_wrap(~stn, ncol=2) +
theme(aspect.ratio=0.5)