---
title: "Dynamic structural equation models"
author: "James T. Thorson"
output: rmarkdown::html_vignette
#output: rmarkdown::pdf_document
vignette: >
  %\VignetteIndexEntry{Dynamic structural equation models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
  %\VignetteDepends{dsem}
---

```{r, include = FALSE}
has_dsem = requireNamespace("dsem", quietly = TRUE)
EVAL <- has_dsem
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = EVAL,
  purl = EVAL
)
# Install locally
#  devtools::install_local( R'(C:\Users\James.Thorson\Desktop\Git\tinyVAST)', force=TRUE )
# Build and PDF
#  setwd(R'(C:\Users\James.Thorson\Desktop\Git\tinyVAST)'); devtools::build_rmd("vignettes/dsem.Rmd"); rmarkdown::render( "vignettes/dsem.Rmd", rmarkdown::pdf_document())
```

```{r setup, echo=TRUE, warning=FALSE, message=FALSE}
library(tinyVAST)
set.seed(101)
options("tinyVAST.verbose" = FALSE)
```

`tinyVAST` includes features to fit a dynamic structural equation model.  We here show this using a bivariate vector autoregressive model for wolf and moose abundance on Isle Royale.

```{r, eval=TRUE, echo=TRUE, message=FALSE, fig.width=6, fig.height=6, warning=FALSE}
data(isle_royale, package="dsem")

# Convert to long-form
data = expand.grid( "time"=isle_royale[,1], "var"=colnames(isle_royale[,2:3]) )
data$logn = unlist(log(isle_royale[2:3]))

# Define cross-lagged DSEM
dsem = "
  # Link, lag, param_name
  wolves -> wolves, 1, arW
  moose -> wolves, 1, MtoW
  wolves -> moose, 1, WtoM
  moose -> moose, 1, arM
  #wolves -> moose, 0, corr
  wolves <-> moose, 0, corr
"

# fit model
mytiny = tinyVAST( spacetime_term = dsem,
                 data = data,
                 times = isle_royale[,1],
                 variables = colnames(isle_royale[,2:3]),
                 formula = logn ~ 0 + var )
mytiny

# Deviance explained relative to both intercepts
# Note that a process-error-only estimate with have devexpl -> 1
deviance_explained( mytiny, 
                    null_formula = logn ~ 0 + var )

# See summary
knitr::kable( summary(mytiny,"spacetime_term"), digits=3 )
```

And we can specifically inspect the estimated interaction matrix:

```{r, eval=TRUE, echo=FALSE, message=FALSE, fig.width=6, fig.height=6}
B = matrix( mytiny$internal$parlist$beta_z[1:4], nrow=2,
            dimnames=list(colnames(isle_royale[,2:3]),colnames(isle_royale[,2:3])) )
knitr::kable( B, digits=3)
```

We can then compare this with package `dsem`

```{r, eval=requireNamespace("dsem"), echo=TRUE, message=FALSE, fig.width=6, fig.height=6, warning=FALSE}
library(dsem)

# Keep in wide-form
dsem_data = ts( log(isle_royale[,2:3]), start=1959)
family = c("normal","normal")

# fit without delta0
# SEs aren't available because measurement errors goes to zero
mydsem = dsem::dsem( sem = dsem,
             tsdata = dsem_data,
             control = dsem_control(getsd=FALSE),
             family = family )
mydsem

# See summary
knitr::kable( summary(mydsem), digits=3 )
```

where we again inspect the estimated interaction matrix:

```{r, eval=TRUE, echo=FALSE, message=FALSE, fig.width=6, fig.height=6}
B = matrix( mydsem$obj$env$parList()$beta_z[1:4], nrow=2,
            dimnames=list(colnames(isle_royale[,2:3]),colnames(isle_royale[,2:3])) )
knitr::kable( B, digits=3)
```