# Dynamic structural equation models

#### James T. Thorson

    library(tinyVAST)
    set.seed(101)
    options("tinyVAST.verbose" = FALSE)

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

    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
    #> Call:
    #> tinyVAST(formula = logn ~ 0 + var, data = data, spacetime_term = dsem,
    #>     times = isle_royale[, 1], variables = colnames(isle_royale[,
    #>         2:3]))
    #>
    #> Run time:
    #> Time difference of 0.349545 secs
    #>
    #> Family:
    #> $obs
    #>
    #> Family: gaussian
    #> Link function: identity
    #>
    #>
    #>
    #>
    #> sdreport(.) result
    #>               Estimate   Std. Error
    #> alpha_j     3.32526197 2.483494e-01
    #> alpha_j     6.44165343 2.116035e-01
    #> beta_z      0.89304266 8.420632e-02
    #> beta_z      0.01420908 1.279150e-01
    #> beta_z     -0.13239996 3.454997e-02
    #> beta_z      0.86147627 7.107386e-02
    #> beta_z     -0.01285890 5.063467e-02
    #> beta_z      0.37727131 3.504314e-02
    #> beta_z      0.17062266 1.584742e-02
    #> log_sigma -12.62369789 2.032840e+04
    #> Maximum gradient component: 3.603845e-05
    #>
    #> Proportion conditional deviance explained:
    #> [1] 1
    #>
    #> spacetime_term:
    #>   heads     to   from parameter start lag    Estimate  Std_Error    z_value
    #> 1     1 wolves wolves         1  <NA>   1  0.89304266 0.08420632 10.6054118
    #> 2     1 wolves  moose         2  <NA>   1  0.01420908 0.12791496  0.1110823
    #> 3     1  moose wolves         3  <NA>   1 -0.13239996 0.03454997 -3.8321293
    #> 4     1  moose  moose         4  <NA>   1  0.86147627 0.07107386 12.1208601
    #> 5     2  moose wolves         5  <NA>   0 -0.01285890 0.05063467 -0.2539545
    #> 6     2 wolves wolves         6  <NA>   0  0.37727131 0.03504314 10.7659099
    #> 7     2  moose  moose         7  <NA>   0  0.17062266 0.01584742 10.7665881
    #>        p_value
    #> 1 2.812223e-26
    #> 2 9.115511e-01
    #> 3 1.270389e-04
    #> 4 8.189516e-34
    #> 5 7.995307e-01
    #> 6 4.986648e-27
    #> 7 4.950067e-27
    #>
    #> Fixed terms:
    #>           Estimate Std_Error  z_value       p_value
    #> varwolves 3.325262 0.2483494 13.38945  6.969636e-41
    #> varmoose  6.441653 0.2116035 30.44210 1.524035e-203

    # 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 )
    #> [1] 1

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

heads | to | from | parameter | start | lag | Estimate | Std_Error | z_value | p_value
--- | --- | --- | --- | --- | --- | --- | --- | --- | ---
1 | wolves | wolves | 1 | NA | 1 | 0.893 | 0.084 | 10.605 | 0.000
1 | wolves | moose | 2 | NA | 1 | 0.014 | 0.128 | 0.111 | 0.912
1 | moose | wolves | 3 | NA | 1 | -0.132 | 0.035 | -3.832 | 0.000
1 | moose | moose | 4 | NA | 1 | 0.861 | 0.071 | 12.121 | 0.000
2 | moose | wolves | 5 | NA | 0 | -0.013 | 0.051 | -0.254 | 0.800
2 | wolves | wolves | 6 | NA | 0 | 0.377 | 0.035 | 10.766 | 0.000
2 | moose | moose | 7 | NA | 0 | 0.171 | 0.016 | 10.767 | 0.000

And we can specifically inspect the estimated interaction matrix:

| wolves | moose
--- | --- | ---
wolves | 0.893 | -0.132
moose | 0.014 | 0.861

We can then compare this with package `dsem`

    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 )
    #>   Coefficient_name Number_of_coefficients   Type
    #> 1           beta_z                      7  Fixed
    #> 2        lnsigma_j                      2  Fixed
    #> 3             x_tj                    122 Random
    mydsem
    #> $par
    #>       beta_z       beta_z       beta_z       beta_z       beta_z       beta_z
    #>   0.84199231   0.06629794  -0.05542667   1.02617142   0.37035485   0.38054542
    #>       beta_z    lnsigma_j    lnsigma_j
    #>   0.82581950 -11.71750317 -12.72180788
    #>
    #> $objective
    #> [1] 102.5012
    #> attr(,"logarithm")
    #> [1] TRUE
    #>
    #> $convergence
    #> [1] 0
    #>
    #> $iterations
    #> [1] 71
    #>
    #> $evaluations
    #> function gradient
    #>      104       72
    #>
    #> $message
    #> [1] "relative convergence (4)"

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

path | lag | name | start | parameter | first | second | direction | Estimate
--- | --- | --- | --- | --- | --- | --- | --- | ---
wolves -> wolves | 1 | arW | NA | 1 | wolves | wolves | 1 | 0.842
moose -> wolves | 1 | MtoW | NA | 2 | moose | wolves | 1 | 0.066
wolves -> moose | 1 | WtoM | NA | 3 | wolves | moose | 1 | -0.055
moose -> moose | 1 | arM | NA | 4 | moose | moose | 1 | 1.026
wolves <-> moose | 0 | corr | NA | 5 | wolves | moose | 2 | 0.370
wolves <-> wolves | 0 | V[wolves] | NA | 6 | wolves | wolves | 2 | 0.381
moose <-> moose | 0 | V[moose] | NA | 7 | moose | moose | 2 | 0.826

where we again inspect the estimated interaction matrix:

| wolves | moose
--- | --- | ---
wolves | 0.842 | -0.055
moose | 0.066 | 1.026

Runtime for this vignette: 1.28 secs

# Works cited

Thorson, James T., Alexander G. Andrews III, Timothy E. Essington, and Scott I. Large. 2024. “Dynamic Structural Equation Models Synthesize Ecosystem Dynamics Constrained by Ecological Mechanisms.” Methods in Ecology and Evolution 15 (4): 744–55. https://doi.org/10.1111/2041-210X.14289.