DoY_Trends_Model.Rmd

---
title: "Model of Long-term Trends in Day-Of-Year Temperatures"
output: html_document
editor_options: 
  chunk_output_type: console
---

R version 4.1.0 (2021-05-18) -- "Camp Pontanezen"

## Get the required packages
```{r}
library(tidyverse)
library(lubridate)
```

Daily high temperature data.
```{r}
TLH.df <- read.csv(file = 'Data/TLH_Daily1940-2020.csv') |>
  mutate(Date = as.Date(DATE)) |>
  mutate(Year = year(Date), 
         month = month(Date, label = TRUE, abbr = TRUE),
         doy = yday(Date),
         MaxTemp = TMAX,
         MinTemp = TMIN) |>
  filter(Year <= 2020 & Year >= 1990)

TLH.df$MaxTemp[TLH.df$Date == "2005-07-08"] <- 95
TLH.df$MinTemp[TLH.df$Date == "2005-07-08"] <- 72
TLH.df$MinTemp[TLH.df$Date == "2002-10-08"] <- 67 # Quincy

TLH.df <-
  TLH.df |>
  select(MaxTemp, Year, doy, Date) |>
  mutate(MaxTempS = scale(MaxTemp),
         YearS = scale(Year),
         doyS = scale(doy))
```

Trends in the data
```{r}
library(broom)

Trends <- 
  TLH.df |> 
    group_by(doy) |>
    do(tidy(lm(MaxTemp ~ Year, data = .))) |>
    filter(term == "Year") |>
    mutate(Date = as.Date(doy - 1, origin = "2020-01-01"),
           estimate = estimate,
           std.error = std.error,
           ymax = estimate + std.error,
           ymin = estimate - std.error)
```

## Linear mixed effects models

```{r}
library(lme4)

fit1 <- lmer(MaxTemp ~  ( Year  | doy), 
             data = TLH.df)

library(blme)

fit1.0 <- blmer(MaxTemp ~ Year + (Year | doy), 
                data = TLH.df)

library(brms)

fit1.1 <-
  brm(data = TLH.df,
      formula = MaxTempS ~ YearS + (0 + YearS | doyS),
      prior = prior(normal(0, 10), class = b),
      iter = 2000, chains = 2, cores = 4,
      seed = 31415)

ps <- fit1.1 |>
  posterior::as_draws_df()
```


## Structural time-series model
```{r}
library(bsts)
```

From: https://www.unofficialgoogledatascience.com/2017/07/fitting-bayesian-structural-time-series.html

Test model
```{r}
data(iclaims)     # bring the initial.claims data into scope

ss <- AddLocalLinearTrend(list(), initial.claims$iclaimsNSA)
ss <- AddSeasonal(ss, initial.claims$iclaimsNSA, nseasons = 52)
model1 <- bsts(initial.claims$iclaimsNSA,
               state.specification = ss,
               niter = 1000)
```

Models
```{r}
x <- data.frame(MaxTemp = TLH.df$MaxTemp, Year = TLH.df$Year, doy = TLH.df$doy)
TLH.z <- zoo(x, order.by = TLH.df$Date)

y <- TLH.z$MaxTemp
ss <- list()
ss <- bsts::AddLocalLinearTrend(ss, y)
ss <- bsts::AddSeasonal(ss, y, nseasons = 31, season.duration = 366)
model1 <- bsts(y,
               state.specification = ss,
               niter = 200)

#plot(model1)
#plot(model1, "components")

ss <- list()
ss <- bsts::AddLocalLinearTrend(ss, TLH.z$MaxTemp)
model2 <- bsts(TLH.z$MaxTemp,
               state.specification = ss,
               niter = 200)

ss <- list()
ss <- bsts::AddSeasonal(ss, TLH.z$MaxTemp, nseasons = 31, season.duration = 366)
model6 <- bsts(MaxTemp ~ Year,
               state.specification = ss,
               niter = 200,
               data = TLH.z)

ss <- list()
ss <- bsts::AddAutoAr(ss, TLH.z$MaxTemp, lags = 1)
model3 <- bsts(MaxTemp ~ Year * doy,
               state.specification = ss,
               niter = 200,
               data = TLH.z)

ss <- list()
ss <- bsts::AddAutoAr(ss, TLH.z$MaxTemp, lags = 14)
ss <- bsts::AddSeasonal(ss, TLH.z$MaxTemp, nseasons = 31, season.duration = 366)
model4 <- bsts(MaxTemp ~ Year,
               state.specification = ss,
               niter = 200,
               data = TLH.z)

CompareBstsModels(list("Model 7" = model7,
                       "Model 3" = model3),
                  colors = c("black", "green"))
```

```{r}
ss <- list()
ss <- bsts::AddAutoAr(ss, TLH.z$MaxTemp, lags = 14)
model7 <- bsts(MaxTemp ~ doy, 
               state.specification = ss,
               niter = 200,
               data = TLH.z)

ss <- list()
ss <- bsts::AddLocalLevel(ss, TLH.z$MaxTemp)
ss <- bsts::AddDynamicRegression(ss, TLH.z$MaxTemp ~ TLH.z$Year)
#ss <- bsts::AddAutoAr(ss, TLH.z$MaxTemp, lags = 14)

model7 <- bsts(TLH.z$MaxTemp, 
               state.specification = ss,
               niter = 200,
               data = TLH.z)


```



## Multilevel model

From: https://paul-buerkner.github.io/brms/articles/brms_distreg.html

Data
```{r}
library(brms)

dat_smooth <- mgcv::gamSim(eg = 6, n = 200, scale = 2, verbose = FALSE)
head(dat_smooth)
```

Model
```{r}
fit_smooth <- brm(
  bf(y ~ s(x1) + s(x2) + (1|fac), sigma ~ s(x0) + (1|fac)),
  data = dat_smooth, family = gaussian(),
  chains = 2, control = list(adapt_delta = 0.95)
)
```

Summaries and plots
```{r}
summary(fit_smooth)
plot(conditional_effects(fit_smooth), points = TRUE, ask = FALSE)
```

Posterior samples of the parameters.
```{r}
ps <- fit_smooth |>
  posterior::as_draws_df()
```

Posterior samples of the outcome.
```{r}
pp <- fit_smooth |>
  brms::posterior_predict() |>
  as.data.frame()
```

See `PredictionModel.Rmd`.

Repeat with daily high temperature data. Problem: The random effect is not seasonal. The first day knows nothing about the 366th day. Should use a structural time series model?

Prepare data to match `dat_smooth`.
```{r}
TLH.df <- read.csv(file = 'Data/TLH_Daily1940-2020.csv') |>
  mutate(Date = as.Date(DATE)) |>
  mutate(Year = year(Date), 
         month = month(Date, label = TRUE, abbr = TRUE),
         doy = yday(Date),
         MaxTemp = TMAX,
         MinTemp = TMIN) |>
  filter(Year <= 2020 & Year >= 1981)

TLH.df$MaxTemp[TLH.df$Date == "2005-07-08"] <- 95
TLH.df$MinTemp[TLH.df$Date == "2005-07-08"] <- 72
TLH.df$MinTemp[TLH.df$Date == "2002-10-08"] <- 67 # Quincy
```

Data
```{r}
max_temp <- 
  TLH.df |>
  dplyr::select(MaxTemp, Year, doy) |>
  mutate(doyF = as.factor(doy))
```

Model
```{r}
fit3 <- brm(
  bf(MaxTemp ~ Year + (1 + Year|doy), sigma ~ (1|doy)),
  data = max_temp, family = gaussian(),
  chains = 2, control = list(adapt_delta = 0.95)
)

plot(conditional_effects(fit3), points = FALSE, ask = FALSE)

max_temp <-
  max_temp |>
  mutate(MaxTempS = scale(MaxTemp))

fit4 <- brm(
  bf(MaxTempS ~ YearS * doyF + (YearS|doyF)),
  data = max_temp, family = gaussian(),
  chains = 2, control = list(adapt_delta = 0.95)
)

plot(conditional_effects(fit4), points = FALSE, ask = FALSE)
```

Compare with lmer models
```{r}
library(lme4)
```


First make a tibble composed of tibbles with `group_by()` followed by `nest()`.
```{r}
by_doy <- TLH.df |>
  group_by(doy) |>
  nest()

head(by_doy)

by_doy$data[[1]]
```

Next fit a linear model to the first doy. Could not get the `update()` function to work. Crashes the session!
```{r}
fit0.1 <- lm(formula = MaxTemp ~ Year, 
           data = by_doy$data[[1]])

library(brms)

fit1.1 <-
  brm(data = by_doy$data[[1]],
      formula = MaxTemp ~ Year,
      prior = prior(normal(0, 10), class = b),
      iter = 4000, chains = 4, cores = 4,
      seed = 31415,
      file = "fits/fit1.1")

fit1.2 <-
  update(fit1.1, 
         newdata = by_doy$data[[2]],
         control = list(adapt_delta = .9),
         file = "fits/fit1.2")

fit1.2 <-
  brm(data = by_doy$data[[2]],
      formula = MaxTemp ~ Year,
      prior = prior(normal(0, 10), class = b),
      iter = 4000, chains = 4, cores = 4,
      seed = 31415)
```