```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning=FALSE, 
                      message=FALSE, results='hide',
                      fig.align = 'center')
library(ggplot2)
library(rstan)
library(atsar)
library(broom)
library(datasets)
library(MARSS)
library(dplyr)
```

## Overview of today's material

* Quick overview of Stan 
* Manipulating and plotting Stan output
* Examples of time series models

## Review of models we've used so far

### Models
* Regression
* ARMA models
* State Space Models
* Dynamic Factor Analysis
* Dynamic Linear Models
* MARSS models (multivariate time series models)

## Why Bayesian?

* Complex hierarchical models  
    + Non-linear models
    + Hierarchical or shared parameters
    + Non-normal data
    
* Inference: what’s the probability that the data are less than some threshold?   

* No bootstrapping!
     + We get credible intervals for parameters and states simultaneously

## Using STAN

### What is Stan?
* Powerful, cross-platform and cross-language (R, Julia, Matlab, etc) that allows users to write custom code that can be called directly from R

* Estimation can be fully or approximate Bayesian inference, or maximum likelihood optimization (BFGS)

* Useful links:
    + [Stan homepage](https://mc-stan.org/)
    + [Stan manual](https://mc-stan.org/users/documentation/)
    + [rstan](https://cran.r-project.org/web/packages/rstan/index.html)

## Options for using Stan in this class

* Write your own code (based on examples in the manual, etc)

* Use an existing package  

* Use our bundled code to get started with simple models (we'll start here)

## Existing packages: rstanarm and bmrs

* Both packages very flexible, and allow same syntax as basic lm/glm or lmer models, e.g.

```{r echo=TRUE, eval=FALSE}
rstan::stan_lm
rstan::stan_glm
rstan::stan_glmer
```

* Vignettes
[brms](http://paul-buerkner.github.io/brms/)
[rstanarm](http://mc-stan.org/rstanarm/articles/rstanarm.html)

## Potential limitations

`brms` includes autocorrelated errors, non-normal data, non-linear smooths, etc. 

* But doesn't allow for AR processes on parameters / latent variables
* Also doesn't include multivariate data

## Advantages

`brms` offers notation that should be very familiar to run many classes of models, 

```{r, echo=TRUE, eval=FALSE}
brms::brm(y ~ x * z + (1|group), data=d)
brms::brm(y01 ~ x * z + (1|group), data=d, family = binomial("logit"))
brms::brm(bf(y ~ s(x)), data=d)
```


## To install code for this class

* We'll need to install these packages to run Stan,  
```{r, eval=FALSE, echo=TRUE}
install.packages("rstan", repos = "https://cloud.r-project.org")
install.packages("devtools", repos = "https://cloud.r-project.org")
```

* And then we can install our custom package for the class with bundled Stan time series models

```{r, eval=FALSE, echo=TRUE}
devtools::install_github(repo="nwfsc-timeseries/atsar")
install.packages("atsar")
```

## Working with Stan output

* We'll start with a simple example, using a classic dataset of flow on the Nile River

```{r}
library(datasets)
data(Nile) #load the data 
plot(Nile,ylab="Flow volume",xlab="Year") 
```

## Working with Stan output

Fitting a regression model to this data doesn't make a lot of sense, but it will introduce us to the basic functionality of our wrapper functions. We'll start with a function named `fit_stan`

```{r, echo=TRUE, eval=TRUE, cache=TRUE}
library(atsar)
lm_intercept = fit_stan(y = as.numeric(Nile), 
                        x = rep(1, length(Nile)), 
                        model_name = "regression")
```

## Working with Stan output

The output of the fitted model can be examined a number of ways. Starting with the simple summaries,

```{r, echo=TRUE}
lm_intercept
```

This model ran 3 MCMC chains (the default) with warmups of 500, followed by 500 more iterations. These latter were stored (so 1500 parameter samples total)

## Working with Stan output

But we're probably more interested in the values for individual parameters. We can pull these out in a few different ways, and plot them with either base graphics or ggplot. As a first option for getting values we can manipulate, we can use the `extract` function,

```{r echo=TRUE}
pars = rstan::extract(lm_intercept)
```

## Plotting with Stan output
Then we can do all kinds of things with this output, like making a histogram
```{r echo=TRUE, fig.height=3, fig.width=4}
hist(pars$beta, 40, col="grey", xlab="Intercept", main="")
```
or calculate summary statistics
```{r echo=TRUE}
quantile(pars$beta, c(0.025,0.5,0.975))
```

## Plotting with Stan output

The object `lm_intercept` is a `stanfit` object, which means there's a lot of other plotting functionality from the rstan package we can use

* First, let's make some traceplots

```{r, echo=TRUE, fig.height=3}
rstan::traceplot(lm_intercept, 
                 pars = c("beta[1]","sigma"), nrow=2, ncol=1)
```

## Plotting with Stan output

* Second, we can examine the correlation between parameters with a pairs plot,

```{r, echo=TRUE, fig.width=4, fig.height=4}
rstan::stan_scat(lm_intercept, pars = c("beta[1]","sigma"))
```

## Plotting with Stan output

* Third, maybe we want to make some density plots or histograms

```{r, echo=TRUE, fig.width=4, fig.height=1.5}
rstan::stan_dens(lm_intercept, pars = c("beta[1]","sigma"))
```

```{r, echo=TRUE, fig.width=4, fig.height=1.5}
rstan::stan_hist(lm_intercept, pars = c("beta[1]","sigma"))
```

## Plotting with Stan output

* Fourth, we can look at implementing some of the plots included in [bayesplot](http://mc-stan.org/bayesplot/)

```{r echo=TRUE, fig.width=7, fig.height=3}
library(bayesplot)
mcmc_areas(as.matrix(lm_intercept),
           pars = c("beta[1]","sigma"),
           prob = 0.8)
```

## Plotting with Stan output

* Another way to show the same uncertainties is with dotplots and credible intervals,

```{r echo=TRUE, fig.width=7, fig.height=3}
rstan::stan_plot(lm_intercept, pars=c("beta[1]","sigma"))
```

## Plotting with Stan output

* We can also look at the autocorrelation of each parameter, which is another useful diagnostic for convergence (we want very low levels of autocorrelation)
```{r echo=TRUE, fig.width=7, fig.height=3}
rstan::stan_ac(lm_intercept, pars=c("beta[1]","sigma"))
```

## Plotting with Stan output
These plots only the tip of the iceberg for plotting. For more great examples of the kinds of plots avaialable, see these vignettes:

* [Examples on Stan](https://mc-stan.org/users/interfaces/bayesplot)

* [Jonah Gabry's introduction to bayesplot](https://cran.r-project.org/web/packages/bayesplot/vignettes/plotting-mcmc-draws.html)

* [Matthew Kay's introduction to bayesplot and tidybayes](https://cran.r-project.org/web/packages/tidybayes/vignettes/tidybayes.html)

## Tidy summaries from Stan output
Using the `broom` package, we can also extract some tidy summaries of the output 
```{r, echo=TRUE, results='markup'}
coef = broom::tidy(lm_intercept)
head(coef)
```

## Tidy summaries from Stan output
These tidy summaries can then be fed into ggplot for example
```{r, echo=TRUE, fig.height=3, fig.width=3}
coef = broom::tidy(lm_intercept)
ggplot(coef[grep("pred",coef$term),], aes(x = 1:100,y=estimate)) + 
  geom_point() + ylab("Estimate +/- SE")+ xlab("")+
  geom_errorbar(aes(ymin=estimate-std.error, ymax=estimate+std.error))
```

## Preserving chain order
For models with multiple chains, we might want to preserve the chain ID to look at individual chain diagnostics. Remember that

* Each chain is independent
* `extract` defaults to merging samples from all chains together, e.g.
```{r echo=TRUE, eval=FALSE}
extract(object, pars, permuted = TRUE)
```
* But summaries can be generated for each combination of parameters-chains by setting
```{r echo=TRUE, eval=FALSE}
extract(object, pars, permuted = FALSE)
```

## More time series models: random walk

This model should be familiar, 

$$E\left[ { Y }_{ t } \right] =E\left[ { Y }_{ t-1 } \right] +{ e }_{ t-1 }$$

* We’ll fit model to temperature data
* Note that the use of the argument `model_name` and `est_drift`. By not estimating drift, we asssume the process is stationary with respect to the mean
```{r echo=TRUE, cache=TRUE}
data(airquality)
Temp = airquality$Temp # air temperature
rw = fit_stan(y = Temp, est_drift = FALSE, model_name = "rw") 
```

## More time series models: random walk

Did the model converge? 

* One quick check is to look at the value of R-hat for each parameter (generally should be small, < 1.05 or smaller)
```{r echo=TRUE, results='markup'}
rw_summary <- summary(rw, pars = c("sigma"), 
                      probs = c(0.1, 0.9))$summary
print(rw_summary)
```

## More time series models: univariate state space models

State equation: 
$${ x }_{ t }={ \phi x }_{ t-1 }+{ \varepsilon  }_{ t-1 }$$
where ${ \varepsilon  }_{ t-1 } \sim Normal(0, q)$

Observation equation:
$${ Y }_{ t } \sim Normal(x_{t}, r)$$

* Let’s compare models with and without the AR parameter $\phi$ in the process model



## More time series models: univariate state space models

We can first run the model with $\phi$,
```{r cache=TRUE, echo=TRUE}
ss_ar = fit_stan(y = Temp, est_drift=FALSE, model_name = "ss_ar")
```

then without,

```{r cache=TRUE, echo=TRUE}
ss_rw = fit_stan(y = Temp, est_drift=FALSE, model_name = "ss_rw")
```

## More time series models: univariate state space models

Estimates from the AR(1) state space model
```{r}
pars_ar = broom::tidy(ss_ar)
pars_ar = pars_ar[grep("pred",pars_ar$term),]
pars_ar$low = pars_ar$estimate - 1.96*pars_ar$std.error
pars_ar$hi = pars_ar$estimate + 1.96*pars_ar$std.error
pars_ar$x = seq(1,nrow(pars_ar))

df = data.frame("x"=pars_ar$x, "y" = Temp, stringsAsFactors=FALSE)
ggplot(pars_ar, aes(x, estimate)) + 
geom_ribbon(aes(ymin=low,ymax=hi),fill="grey70") + 
xlab("Time") + ylab("Estimate (AR)") + 
geom_point(data=df, aes(x,y), col="blue", alpha=0.6, size=2)
```

## More time series models: univariate state space models

Estimates from the RW state space model
```{r}
pars_rw = broom::tidy(ss_rw)
pars_rw = pars_rw[grep("pred",pars_rw$term),]
pars_rw$low = pars_rw$estimate - 1.96*pars_rw$std.error
pars_rw$hi = pars_rw$estimate + 1.96*pars_rw$std.error
pars_rw$x = seq(1,nrow(pars_rw))

df = data.frame("x"=pars_rw$x, "y" = Temp, stringsAsFactors=FALSE)
ggplot(pars_rw, aes(x, estimate)) + 
geom_ribbon(aes(ymin=low,ymax=hi),fill="grey70") + 
xlab("Time") + ylab("Estimate (AR)") + 
geom_point(data=df, aes(x,y), col="blue", alpha=0.6, size=2)
```

## More time series models: univariate state space models

Estimates from both models (note the difference in credible interval widths)
```{r}
pars_rw$model = "RW"
pars_ar$model = "AR"
pars = rbind(pars_rw,pars_ar)
df = data.frame("x"=pars_rw$x, "y" = Temp, stringsAsFactors=FALSE)
ggplot(pars, aes(x, estimate, fill=model, group=model)) + 
geom_ribbon(aes(ymin=low,ymax=hi,fill=model, alpha=0.3)) + 
xlab("Time") + ylab("Estimate") + 
geom_point(data=df, aes(x,y,group=NA,fill=NA), col="blue", alpha=0.6, size=2)
```

## More time series models: univariate state space models

We might be also interested in looking at posterior estimates for these models.

* What is the posterior distribution of $\phi$?

```{r echo=TRUE, fig.width=3, fig.height=1.5}
rstan::stan_dens(ss_ar, c("phi"))
```

* This shows with a value near 1, the behavior of the model should be very similar to the random walk (from the predictions, it is!). The multi-modal distribution probably is an indicator of convergence issues

## More time series models: univariate state space models

We might also be interested in derived quantities. 

* For example, what's the probability of temperature exceeding 90 degrees? 
```{r, echo=TRUE, results='markup'}
pars = extract(ss_ar)
p = length(which(pars$pred > 90)) / length(pars$pred)
print(round(p,3))
```

## Dynamic Factor Analysis (DFA)
As an example of Bayesian DFA, we'll load the WA plankton dataset. There's 3 versions of the data, and we'll use the one that's been transformed. As a reminder,   

    + 0s replaced with NAs
    + data have been z-scored
    + we only use data from 1980-1989 for simplicity  
```{r echo=TRUE}
dat = lakeWAplanktonTrans
plankdat = dat[dat[,"Year"]>=1980 & dat[,"Year"]<1990,]
phytoplankton = c("Cryptomonas", "Diatoms", "Greens",
                   "Unicells", "Other.algae")
dat.spp.1980 = plankdat[,phytoplankton]

```

## Dynamic Factor Analysis (DFA)

```{r}
reshape2::melt(dat.spp.1980) %>% 
  rename(Time = Var1, species = Var2) %>% 
ggplot(aes(Time, value)) + 
  geom_line() + facet_wrap(~species) + 
  xlab("Time") + ylab("Ln (abundance)")
```

## Running the model

For starters, we'll try a 3-trend model, 

```{r echo=TRUE, eval=FALSE}
mod_3 = fit_dfa(y = t(dat.spp.1980), num_trends=3)
```

## Trends need to be rotated (like MARSS)

Why? Identifiability

* Try re-ordering time series, and running DFA on each new dataset
* Results will be sensitive to time series order, BUT not after rotation
* Like MARSS we’ll use varimax rotation

Use function we’ve written, rotate_trends
```{r echo=FALSE, eval=FALSE}
rot = rotate_trends(mod_3)
```

## Attributes of rotated object

`Z_rot`, rotated Z matrix for each MCMC draw  
`trends`, rotated trends for each MCMC draw  
`Z_rot_mean`, mean Z across draws  
`trends_mean`, mean trends across draws  
`trends_lower`, lower 2.5% on trends  
`trends_upper`, upper 97.5% on trends  

## Other variance structures

By default, the structure is `diagonal and equal`
```{r echo=FALSE, eval=FALSE}
mod_3 = fit_dfa(y = t(dat.spp.1980), num_trends=3)
```

Diagonal and unequal or shared variances can also be specified using `varIndx` argument. The diagonal and unequal structure would be called with 
```{r echo=TRUE, eval=FALSE}
mod_3 = fit_dfa(y = t(dat.spp.1980), num_trends=3, varIndx = rep(1,5))
```

## Uncertainty intervals on states

We often just present the trend estimates for DFA – but not uncertainty

Let’s look at effect of missing data on DFA for the harbor seal dataset

```{r echo=TRUE}
data("harborSealWA")
```

## Model
Assume single trend for the population

```{r}
matplot(harborSealWA[,-1], type="l", lwd=3, ylab="Ln (abundance)")
```


## Extracting the predicted trend

We'll extract predictions from the best model, 
```{r echo=TRUE, cache=TRUE, fig.height=2, fig.align="center"}
pars = extract(fit_dfa(y = t(harborSealWA[,-1]), num_trends = 1))
df = data.frame("time"=1:nrow(harborSealWA),
                "pred"=apply(pars$x[,1,], 2, mean),
                "low"=apply(pars$x[,1,], 2, quantile,0.025),
                "hi"=apply(pars$x[,1,], 2, quantile,0.975))
ggplot(df, aes(time, pred)) + 
  geom_ribbon(aes(ymin=low,ymax=hi), fill="grey30", alpha=0.5)+
  geom_line()
```

## Fitting a DLM with time varying intercept

For comparison to MARSS, we'll use Mark's example of logit-transformed survival from the Columbia River. We can think about setting the DLM up in the slope or the intercept. For this first example, we'll do the latter. 

```{r fig.height=1.5}
ggplot(SalmonSurvCUI, aes(year, logit.s)) + geom_line() + geom_point()
```

```{r cache=TRUE, echo=TRUE}
mod = fit_stan(y = SalmonSurvCUI$logit.s, 
               model_name="dlm-intercept")
```

## Fitting a DLM with time varying intercept
```{r echo=TRUE, fig.height=2}
df = data.frame("year"=1:42,
"pred"=apply(extract(mod, "intercept", permuted=FALSE), 3, mean),
"low"=apply(extract(mod, "intercept", permuted=FALSE), 3, quantile,0.025),
"hi"=apply(extract(mod, "intercept", permuted=FALSE), 3, quantile,0.975))
ggplot(df, aes(year,pred)) + 
  geom_ribbon(aes(ymin=low,ymax=hi),fill="grey30",alpha=0.5) + 
  geom_line() + ylab("Intercept") + xlab("Time")
```



## Fitting a DLM with single intercept and time-varying slope

```{r cache=TRUE, echo=TRUE}
mod = fit_stan(y = SalmonSurvCUI$logit.s, 
x = SalmonSurvCUI$CUI.apr, model_name="dlm-slope")
```

```{r, fig.height=2}
df = data.frame("year"=1:42,
"pred"=apply(extract(mod, "beta", permuted=FALSE), 3, mean),
"low"=apply(extract(mod, "beta", permuted=FALSE), 3, quantile,0.025),
"hi"=apply(extract(mod, "beta", permuted=FALSE), 3, quantile,0.975))
ggplot(df, aes(year,pred)) + 
  geom_ribbon(aes(ymin=low,ymax=hi),fill="grey30",alpha=0.5) + 
  geom_line() + ylab("Beta[1]") + xlab("Time")
```


## Fitting a DLM time-varying intercept and time-varying slope

* Use model.matrix() to specify x
```{r echo=TRUE, eval=FALSE}
lmmod = lm(SalmonSurvCUI$logit.s ~ SalmonSurvCUI$CUI.apr)
x = model.matrix(lmmod)
```

```{r echo=TRUE, eval=FALSE}
lmmod = lm(SalmonSurvCUI$logit.s ~ SalmonSurvCUI$CUI.apr)
mod = fit_stan(y = SalmonSurvCUI$logit.s, 
               x = model.matrix(lmmod), 
               model_name="dlm"
```

## Summary

* Bayesian implementation of time series models in Stan can do everything that MARSS can do and more!

* Very flexible language, great developer community 

* Widely used by students in SAFS / UW / QERM / etc

* Please come to us with questions or modeling issues!


