---
title: "Dynamic Factor Analysis"
output: 
  rmarkdown::html_vignette:
    toc: yes
vignette: >
  %\VignetteIndexEntry{Dynamic Factor Analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

\newcommand{\QQ}{\mathbf{Q}}
\newcommand{\RR}{\mathbf{R}}
\newcommand{\xx}{\mathbf{x}}

```{r setup, include = FALSE}
# knitr options
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
run.comparisons <- TRUE
set.seed(1234)
```

The `MARSS()` function allows you to fit DFAs with the same form as `MARSS(x, form="dfa")`. This has a diagonal $\QQ$ with 1 on the diagonal and a stochastic $\xx_1$ with mean 0 and variance of 5 (diagonal variance-covariance matrix). There are only 3 options allowed for $\RR$:
diagonal and equal, diagonal and unequal, and unconstrained.

## Example data

```{r}
library(MARSS)
library(dplyr)
data(lakeWAplankton, package = "MARSS")
phytoplankton <- c("Cryptomonas", "Diatoms", "Greens", "Unicells", "Other.algae")
dat <- as.data.frame(lakeWAplanktonTrans) |>
  subset(Year >= 1980 & Year <= 1989) |>
  subset(select=phytoplankton) |>
  t() |>
  MARSS::zscore()
```

## Fit models without covariates

```{r}
mod.list <- list(R='unconstrained', m=1, tinitx=1)
```

Fit with MARSS with EM or optim and BFGS.
```{r}
m1 <- MARSS(dat, model=mod.list, form='dfa', z.score=FALSE, silent = TRUE)
m2 <- MARSS(dat, model=mod.list, form='dfa', z.score=FALSE, silent = TRUE, method="BFGS")
```

Fit with TMB.
```{r}
library(marssTMB)
m3 <- dfaTMB(dat, model=list(m=1, R='unconstrained'))
m4 <- MARSS(dat, model=mod.list, form='dfa', method="TMB", silent=TRUE)
m5 <- MARSS(dat, model=mod.list, method="BFGS_TMB", form='dfa', silent=TRUE)
```

### Log likelihoods

```{r echo=FALSE}
df <- data.frame(name=c("MARSS-EM", "MARSS-BFGS", "dfaTMB-nlminb", "MARSS_tmb-nlminb", "MARSS_tmb-optim-BFGS"),
                 logLik=c(m1$logLik, m2$logLik, m3$logLik, m4$logLik, m5$logLik))
knitr::kable(df)
```

### Compare parameter estimates

```{r echo=FALSE}
library(tidyr)
library(ggplot2)
# BFGS was swapped in sign so just mult by -1
adj <- 1
if(sign(coefficients(m1)$Z[1]) != sign(coefficients(m2)$Z[1])) adj <- -1
pars <- data.frame(
  name = c(paste0("R", rownames(coefficients(m1)$R)), paste0("Z", rownames(coefficients(m1)$Z))),
  EM = c(coefficients(m1)$R, coefficients(m1)$Z),
  BFGS = c(coefficients(m2)$R, adj*coefficients(m2)$Z),
  TMB_nlminb = c(as.vector(m3$Estimates$R[lower.tri(m3$Estimates$R,diag=TRUE)]), as.vector(m3$Estimates$Z)),
  TMB_BFGS = c(coefficients(m4)$R, coefficients(m4)$Z))
pars <- pars %>% tidyr::pivot_longer(!name, names_to = "model")

dodge <- position_dodge(width=0.5)
ggplot(pars, aes(x=name, y=value, col=model)) +
  geom_point(position=dodge) +
  ggtitle("same estimates")
```

## Add example with covariates

For `form="dfa"`, pass in covariates with `covariates=xyz`. If using the default model form (not dfa), then pass in covariates with `model$d` or `model$c`.

### Fit model

```{r}
# use a simpler R
mod.list2 <- list(m=1, R='diagonal and unequal', tinitx=1)
# add a temperature covariate
temp <- as.data.frame(lakeWAplanktonTrans) |>
    subset(Year >= 1980 & Year <= 1989) |>
    subset(select=Temp)
covar <- t(temp) |> zscore()
m.cov1 <- MARSS(dat, model=mod.list2, form="dfa", covariates=covar, silent = TRUE, z.score = FALSE, method="TMB")
```

Add a 2nd covariate

```{r}
TP <- as.data.frame(lakeWAplanktonTrans) |>
    subset(Year >= 1980 & Year <= 1989) |>
    subset(select=TP)
covar <- rbind(covar, t(TP)) |> zscore()
m.cov2 <- MARSS(dat, model=mod.list2, form="dfa", covariates=covar, silent = TRUE, z.score=FALSE, method="TMB")
```

### Parameter estimates

Model with one covariate.
```{r}
tidy(m.cov1)
```

Model with two covariates.
```{r}
tidy(m.cov2)
```


### Plot estimates

The estimated factors for one covariate model.
```{r}
library(ggplot2)
autoplot(m.cov1, plot.type="xtT")
```

The estimated factors for two covariate model.
```{r}
library(ggplot2)
autoplot(m.cov2, plot.type="xtT")
```