a1_bayesdfa.Rmd

---
title: "Overview of the bayesdfa package"
author: "Eric J. Ward, Sean C. Anderson, Mary E. Hunsicker, Mike A. Litzow, Luis A. Damiano, Mark D. Scheuerell, Elizabeth E. Holmes, Nick Tolimieri"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Overview of the bayesdfa package}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

Here we will use the bayesdfa package to fit dynamic factor analysis (DFA) model to simulated time series data. In addition to working through an example of DFA for multivariate time series, we'll apply bayesdfa routines for fitting hidden Markov models (HMMs) to the estimated trends to identify latent regimes. Most of the core functions of the package are included here, including `fit_dfa()` and `find_dfa_trends` for fitting DFA models, `plot_trends()`, `plot_fitted()` and `plot_loadings()` for plotting estimates, `find_swans()` for flagging extremes, `fit_regimes()` and `find_regimes()` for fitting HMM models, and `plot_regimes()` for plotting HMM output. 
   
```{r set-knitr-options, cache=FALSE, echo=FALSE, warning=FALSE, message=FALSE}
library("knitr")
opts_chunk$set(message = FALSE, fig.width = 5.5)
```

Let's load the necessary packages:

```{r, message=FALSE, warning=FALSE}
library(bayesdfa)
library(ggplot2)
library(dplyr)
library(rstan)
chains = 1
iter = 10
```

## Introduction to the DFA model

We adopt the same notation used in the MARSS package for dynamic factor analysis models. The DFA model consists of two models, one describing the latent process or states, and an observation or data model linking the process model to observations. Slight variations on this model are described below, but the process model for the basic DFA model is written as a multivariate random walk,

$$\textbf{x}_{t+1}= \textbf{x}_{t} + \textbf{w}_{t}$$
where the matrix $\textbf{x}$ is dimensioned as the number of years $N$ by the number of latent trends $K$. The process error is assumed to be multivariate normal, $\textbf{w}_{t} \sim MVN(0,\textbf{Q})$ where $\textbf{Q}$ is generally assumed to be a $K$-by-$K$ identity matrix. 

The observation or data model linking $\textbf{x}_{t}$ to observed data $\textbf{y}_{t}$ is 

$$\textbf{y}_{t} = \textbf{Zx}_{t} + \textbf{B}\textbf{d}_{t} + \textbf{e}_{t}$$
The matrix $Z$ represents a matrix of estimated loadings, dimensioned as number of time series $P$ by number of latent trends $K$. Optional covariates $\textbf{d}_{t}$ are included in the observation model with estimated coefficients $B$. The residual errors $\textbf{e}_{t}$ are assumed to be normally distributed, e.g. $\textbf{e}_{t} \sim MVN(0, \textbf{R})$. There are a number of choices for $\textbf{R}$ -- these can be a diagonal matrix with equal or unequal elements, or an unconstrained covariance matrix.  

## DFA model with no extreme events  

First, let's simulate some data. We will use the built-in function `sim_dfa()`, but normally you would start with your own data. We will simulate 20 data points from 4 time series, generated from 2 latent processes. For this first dataset, the data won't include extremes, and loadings will be randomly assigned (default).

```{r simulate-data}
set.seed(1)
sim_dat <- sim_dfa(
  num_trends = 2,
  num_years = 20,
  num_ts = 4
)
```

```{r simulate-data-plot, echo=FALSE, fig.align='center', fig.cap="Simulated data, from a model with 2 latent trends and no extremes.\\label{fig:simulate-data-plot}"}
df = data.frame("time" = rep(1:20,4),"y"=c(t(sim_dat$y_sim)), 
                "ts"=as.factor(sort(rep(1:4,20))))
ggplot(df, aes(time,y,group=ts,col=ts)) + geom_line() + theme_bw() + 
  xlab("Time") + ylab("Observed data")
```

Next, we'll fit a 1-trend, 2-trend, and 3-trend DFA model to the simulated time series using the `fit_dfa()` function. Starting with the 1-trend model, we'll estimate the posterior distributions of the trends and loadings. Note that this example uses 1 MCMC chain and 50 iterations --- for real examples, you'll want to use more (say 4 chains, 5000 iterations).

```{r fit-1-trend, message=FALSE, warning=FALSE, results='hide'}
f1 <- fit_dfa(
  y = sim_dat$y_sim, num_trends = 1, scale="zscore",
  iter = iter, chains = chains, thin = 1
)
```

Convergence of DFA models can be evaluated with our `is_converged()` function. This function takes a fitted object, and specified `threshold` argument representing the maximum allowed Rhat value (default = 1.05). The convergence test isn't that useful for a model with such a short number of iterations, but is called with

```{r}
is_converged(f1, threshold = 1.05)
```

This function evaluates Rhat values for all parameters and log likelihood values - so be sure to check what's not converging if the model is not passing this test. 

Before we extract the trends from the model, we need to rotate the loadings matrix and trends. By default we use the varimax rotation, implemented in the `rotate_trends()` function. An optional argument is the `conf_level` argument, which calculates the specified confidence (credible) interval of the estimates (by default, this is set to 0.95).

```{r rot-1-trend, warning=FALSE, message=FALSE, results='hide'}
r <- rotate_trends(f1)
```

The rotated object has several quantities of interest, including the mean values of the trends "trends_mean" and loadings "Z_rot_mean",

```{r}
names(r)
```

We can then plot the trends and intervals, with 
```{r plot-1-trend, fig.align='center', fig.cap="Estimated trend and 95% CI for a 1-trend DFA model applied to simulated data.\\label{fig:simulate-data-plot}"}
plot_trends(r) + theme_bw()
```

We can also plot the estimated loadings (we'll show that plot for the more complicated 2-trend model below because it's not as interesting for the 1-trend model) and the fitted values. To plot the fitted values from the 1-trend model, we'll use the `plot_fitted()` function (predicted values can also be returned without a plot, with the `predicted()`) function.

The trends and intervals are plotted, faceting by time series, with 
```{r plot-1-fitted-example, fig.align='center', fig.cap="Model predicted values from the 1-trend DFA model applied to simulated data.\\label{fig:fitted-example}"}
plot_fitted(f1) + theme_bw()
```

Moving to a more complex model, we'll fit the 2-trend and 3-trend models. All other arguments stay the same as before,

```{r fit-models, warning=FALSE, results='hide', message=FALSE}
f2 <- fit_dfa(
  y = sim_dat$y_sim, num_trends = 2, scale="zscore",
  iter = iter, chains = chains, thin = 1
)
r2 <- rotate_trends(f2)

f3 <- fit_dfa(
  y = sim_dat$y_sim, num_trends = 3, scale="zscore",
  iter = chains, chains = chains, thin = 1
)
r3 <- rotate_trends(f3)
```

The fits from the 2-trend model look considerably better than that from the 1-trend model,
```{r plot-2-fitted-example, fig.align='center', fig.cap="Model predicted values from the 2-trend DFA model applied to simulated data.\\label{fig:fitted-example2}"}
plot_fitted(f2) + theme_bw()
```

The loadings from the 1-trend model aren't as interesting because for a 1-trend model the loadings are a 1-dimensional vector. For the 2 trend model, there's a separate loading of each time series on each trend,

```{r}
round(r2$Z_rot_mean, 2)
```

These loadings can also be plotted with the `plot_loadings()` function. This shows the distribution of the densities as violin plots, with color proportional to being different from 0.

```{r plot-loadings, fig.align='center', fig.cap="Estimated loadings from the 2-trend DFA model.\\label{fig:plot-loadings}"}
plot_loadings(r2) + theme_bw()
```

Finally, we might be interested in comparing some measure of model selection across these models to identify whether the data support the 1-trend, 2-trend, or 3-trend models. The Leave One Out Information Criterion can be calculated with the `loo()` function, for example the LOOIC for the 1-trend model can be accessed with

```{r loo, warning=FALSE, message=FALSE}
loo1 <- loo(f1)
loo1$estimates
```

where `r loo1$estimates["looic","Estimate"]` is the estimate and `r loo1$estimates["looic","SE"]` is the standard error. 

As an alternative to fitting each model individually as we did above, we also developed the `find_dfa_trends()` to automate fitting a larger number of models. In addition to evaluating different trends, this function allows the user to optionally evaluate models with normal and Student-t process errors, and alternative variance structures (observation variance of time series being equal, or not). For example, to fit models with 1:5 trends, both Student-t and normal errors, and equal and unequal variances, the call would be

```{r eval=FALSE}
m <- find_dfa_trends(
  y = s$y_sim, iter = iter,
  kmin = 1, kmax = 5, chains = chains, compare_normal = TRUE,
  variance = c("equal", "unequal")
)
```

## DFA model with extreme events 

In this example, we'll simulate data with an extreme anomaly. The biggest difference between this model and the conventional model is that in the DFA process model, 

$$\textbf{x}_{t+1}= \textbf{x}_{t} + \textbf{w}_{t}$$

instead of $\textbf{w}_{t}$ being normally distributed, we assume $\textbf{w}_{t}$ is Student-t distributed. With multiple trends, this becomes a multivariate Student-t, 

$$\textbf{w}_{t} \sim MVT(\nu, 0, \textbf{Q})$$
The parameter $\nu$ controls how much the tails of this distribution deviate from the normal, with smaller values ($\nu$ closer to 2) resulting in more extreme anomalies, and large values ($\nu$ closer to 30) resulting in behavior similar to a normal distribution. 

As before, this will be 20 data points from 4 time series, generated from 2 latent processes. The `sim_dfa()` function's arguments `extreme_value` and `extreme_loc` allow the user to specify the magnitude of the extreme (as an additive term in the random walk), and the location of the extreme (defaults to the midpoint of the time series). Here we'll include an extreme value of 6,

```{r simulate-data2}
set.seed(1)
sim_dat <- sim_dfa(
  num_trends = 2,
  num_years = 20, num_ts = 4,
  extreme_value = 6
)
```

Plotting the data shows the anomaly occurring between time step 9 and 10,
```{r simulate-data-plot2, fig.align='center', fig.cap="Simulated data, from a model with 2 latent trends and an extreme in the midpoint of the time series.\\label{fig:simulate-data-plot2}"}
matplot(t(sim_dat$y_sim), type = "l", ylab = "Response", xlab = "Time")
```

Though the plot is a little more clear if we standardize the time series first,

```{r simulate-data-plot3, fig.align='center', fig.cap="Simulated data (z-scored), from a model with 2 latent trends and an extreme in the midpoint of the time series.\\label{fig:simulate-data-plot3}"}
matplot(scale(t(sim_dat$y_sim)), type = "l", ylab = "Response", xlab = "Time")
```

Instead of fitting a model with normal process deviations, we may be interested in fitting the model with Student-t deviations. We can turn on the estimation of `nu` with the `estimate_nu` argument. [Alternatively, nu can also be fixed a priori by setting the argument `nu_fixed`]. Here's the code for a 2-trend model with Student-t deviations,

```{r fit-2-trend-extreme, message=FALSE, warning=FALSE, results='hide'}
t2 <- fit_dfa(
  y = sim_dat$y_sim, num_trends = 2, scale="zscore",
  iter = iter, chains = chains, thin = 1, estimate_nu = TRUE
)
```

Again we have to rotate the trends before plotting,
```{r fit-extreme-dfa, fig.align='center', fig.cap="Estimated trends, from a model with 2 latent trends Student-t deviations.\\label{fig:fit-extreme-dfa}"}
r <- rotate_trends(t2)
plot_trends(r) + theme_bw()
```

And the loadings,

```{r plot-extreme-loadings, fig.align='center', fig.cap="Estimated loadings, from a model with 2 latent trends Student-t deviations.\\label{fig:plot-extreme-loadings}"}
plot_loadings(r) + theme_bw()
```

One way to look for extremes is using the `find_swans()` function, which evaluates the probability of observing a deviation in the estimated trend (or data) greater than what is expected from a normal distribution. This function takes a `threshold` argument, which specifies the cutoff. For example, to find extremes greater than 1 in 1000 under a normal distribution, the function call is 

```{r eval=FALSE}
find_swans(r, plot = FALSE, threshold = 1 / 1000)
```
Setting plot to TRUE also creates a time series plot that flags these values.

We can also look at the estimated `nu` parameter, which shows some support for using the Student-t distribution (values greater than ~ 30 lead to similar behavior as a normal distribution), 

```{r summarize-nu}
summary(rstan::extract(t2$model, "nu")[[1]])
```

## Fitting DFA models with non-Gaussian families

We've implemented a number of alternative families for cases when the response variable might be non-normally distributed. These alternative families may be specified with the family argument as a text string in the `fit_dfa` function, e.g.

```{r echo=TRUE, eval=FALSE}
f <- fit_dfa(..., family = "poisson")
```

The currently supported families can be specified as any of the following -- the link functions are
currently hard-coded, and included in the table below.

```{r echo=FALSE, results='asis'}
m = cbind(c("gaussian","lognormal","gamma","binomial","poisson","nbinom2"), 
                 c("identity","log","log","logit","log","log"))
colnames(m) = c("Family","link")
knitr::kable(m)
```

## Alternative loadings for DFA models

By default, the loadings matrix in a DFA is constrained by zeros. For example, a 3-trend model applied to 5 time series would have a loadings matrix that was constrained as 

```{r, echo=FALSE, results='asis'}
Z = matrix("",5,3)
for(i in 1:5) {
  for(j in 1:3) {
    Z[i,j] = paste0("z[",i,",",j,"]")
  }
}
Z[1,2:3] = "0"
Z[2,3] = "0"

colnames(Z) = c("Trend 1","Trend 2","Trend 3")
knitr::kable(Z)
```

As an alternative, we may wish to fit a model where each time series arises as a mixture of the trends. In this case, the loadings matrix would be 

```{r, echo=FALSE, results='asis'}
Z = matrix("",5,3)
for(i in 1:5) {
  for(j in 1:3) {
    Z[i,j] = paste0("z[",i,",",j,"]")
  }
}
colnames(Z) = c("Trend 1", "Trend 2", "Trend 3")
knitr::kable(Z)
```

And the added constraint is that each row sums to 1, e.g.

$$\sum_{j=1:3} Z_{1,j} = 1$$
## Including autoregressive (AR) or moving-average (MA) components on trends

For some models, it may be appropriate to include autoregressive or moving average components to model the latent trends. We've implemented 1st - order components on each, though by default these are not included. 

To include the AR(1) component on the trend, you can specify

```{r eval = FALSE}
fit <- fit_dfa(..., estimate_trend_ar = TRUE)
```

This results in a model where trend $i$ is modeled as 

$$x_{i,t+1} = \phi_{i}*x_{i,t} + \delta_{i,t}$$
Each trend is allowed to have a unique AR(1) parameter, $\phi_{i}$.

In conventional DFA models, the process deviations are assumed to be independent, e.g. $\delta_{i,t} ~ Normal(0, r)$. By including a MA(1) component on the trends, these terms may be modeled as 

$$delta_{i,t+1} \sim Normal( \theta_{i}*delta_{i,t}, q_{i})$$

where $\theta_{i}$ is the trend-specific MA parameter, $q$ is the process variance, and usually constrained to be not estimated and fixed at 1. 

## Applying Hidden Markov Models to identify latent regimes

Finally, we might be interested in evaluating evidence for the estimated trends changing between multiple regimes. We've implemented HMMs in the functions `fit_regimes()` and `find_regimes()`. `fit_regimes` fits a HMM with a pre-specified number of latent states, while `find_regimes()` loops over a range of models so that support can be evaluated with LOOIC.

We'll illustrate an example of the `fit_regimes()` function. Note that in the current version of the package, this has to be applied to each trend separately. Also, uncertainty estimates from the DFA trends may also be included as data (instead of estimated) -- this may be particularly useful for datasets with lots of missing values in portions of the time series (where the uncertainty balloons up).

Applying the 2-regime model to the first trend,
```{r eval = FALSE, warning=FALSE, message=FALSE, results='hide'}
reg_mod <- fit_regimes(
  y = r$trends_mean[1, ],
  sds = (r$trends_upper - r$trends_mean)[1, ] / 1.96,
  n_regimes = 2,
  iter = 50, chains = 1
)
```

In addition to getting diagnostics and quantities like LOOIC out of the `reg_mod` object, we can plot the estimated states. 

```{r plot-regime, eval=FALSE, fig.align='center', fig.cap="Estimated regimes, from a HMM model applied to the first trend of a 2-trend DFA model with Student-t deviations.\\label{fig:plot-regime}"}
plot_regime_model(reg_mod)
```

In this case, there might be some support for the 2-regime model. Sometimes (but not in this example) label switching makes the plots a little challenging to interpret, and the labels need to be reversed. We've included the argument `flip_regimes` to help with this,

```{r plot-regime-flipped, eval=FALSE, fig.align='center', fig.cap="Estimated regimes (after flipping), from a HMM model applied to the first trend of a 2-trend DFA model with Student-t deviations.\\label{fig:plot-regime-flipped}"}
plot_regime_model(reg_mod, flip_regimes = TRUE)
```

## DFA model with weights

We allow weights to be used in DFA models in two ways. In the first form, inverse variance weighting is used to adjust observations based on some standard error associated with each observation. Specifically, the weights are included by modifying each variance to be $$\sigma^2 / w_i$$. As a concrete example, we'll simulate a dataset, add some examples of standard errors on the survey indices, and then perform the DFA.

Our simulated standard errors are the same for all surveys -- except time series 2, which is much more precise.  

```{r}
set.seed(1)
sim_dat <- sim_dfa(
  num_trends = 2,
  num_years = 20,
  num_ts = 4
)

df <- data.frame(obs = c(sim_dat$y_sim), time = sort(rep(1:20,4)),
                 ts = rep(1:4,20))
df$se <- runif(nrow(df), 0.6, 0.8)
df$se[which(df$ts == 2)] = 0.2
```

Next we can generate the weights (this is redundant, and "se" could be used instead in the function call below). Because the weights are used as an offset, $$\sigma^2 / w_i$$, we don't want to use the SE alone as a weight but make them inversely related to the SE. As a quick note, the scale of these may affect estimation and some additional normalization may be needed (rather than standard errors, it may be more helpful to think about the sample size each data point represents).

```{r}
df$weights <- (1 / df$se)^2
```

And fit the model with the weights argument

```{r}
f2 <- fit_dfa(
  y = df, num_trends = 2, scale="zscore",
  iter = 500, chains = 1, thin = 1,
  inv_var_weights = "weights", data_shape = "long"
)
```

As a second type of weighting, we also have implemented weights in the same form used in other widely used packages (glmmTMB, sdmTMB, brms, etc). In this case, weights are used as multipliers on the log-likelihood of each observation. To specify these kinds of weights, we use the `likelihood_weights` argument instead, where observations with higher weights contribute more to the total log likelihood.

```{r eval=FALSE}
f2 <- fit_dfa(
  y = df, num_trends = 2, scale="zscore",
  iter = 500, chains = 1, thin = 1,
  likelihood_weights = "weights", data_shape = "long"
)
```