library(tinyVAST)
set.seed(101)
options("tinyVAST.verbose" = 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
#> 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.5426672 secs
#> 
#> Family: 
#> $obs
#> 
#> Family: gaussian 
#> Link function: identity 
#> 
#> 
#> 
#> 
#> sdreport(.) result
#>               Estimate   Std. Error
#> alpha_j     3.32526205 2.483494e-01
#> alpha_j     6.44165376 2.116033e-01
#> beta_z      0.89304267 8.420630e-02
#> beta_z      0.01420903 1.279150e-01
#> beta_z     -0.13239998 3.454996e-02
#> beta_z      0.86147617 7.107382e-02
#> beta_z     -0.01285884 5.063467e-02
#> beta_z      0.37727135 3.504315e-02
#> beta_z      0.17062264 1.584742e-02
#> log_sigma -12.67123025 2.093639e+04
#> Maximum gradient component: 6.355499e-06 
#> 
#> 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.89304267 0.08420630 10.6054137
#> 2     1 wolves  moose         2  <NA>   1  0.01420903 0.12791499  0.1110818
#> 3     1  moose wolves         3  <NA>   1 -0.13239998 0.03454996 -3.8321308
#> 4     1  moose  moose         4  <NA>   1  0.86147617 0.07107382 12.1208650
#> 5     2  moose wolves         5  <NA>   0 -0.01285884 0.05063467 -0.2539534
#> 6     2 wolves wolves         6  <NA>   0  0.37727135 0.03504315 10.7659086
#> 7     2  moose  moose         7  <NA>   0  0.17062264 0.01584742 10.7665891
#>        p_value
#> 1 2.812165e-26
#> 2 9.115515e-01
#> 3 1.270382e-04
#> 4 8.189021e-34
#> 5 7.995316e-01
#> 6 4.986720e-27
#> 7 4.950011e-27
#> 
#> Fixed terms: 
#>           Estimate Std_Error  z_value       p_value
#> varwolves 3.325262 0.2483494 13.38945  6.969496e-41
#> varmoose  6.441654 0.2116033 30.44212 1.523171e-203
#> 
#> Sanity check:

# 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 )

library(dsem)

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

# fit without delta0
# Fitting with family = "fixed" seems stable on CRAN checks
mydsem = dsem::dsem( 
  sem = dsem,
  tsdata = dsem_data,
  control = dsem_control(
    getsd = FALSE
  ), 
)
#>   Coefficient_name Number_of_coefficients   Type
#> 1           beta_z                      7  Fixed
#> 2             mu_j                      2 Random
mydsem
#> $par
#>       beta_z       beta_z       beta_z       beta_z       beta_z       beta_z 
#>  0.895834720  0.007358847 -0.124879928  0.874884847 -0.014394603  0.378522244 
#>       beta_z 
#>  0.172997994 
#> 
#> $objective
#> [1] 7.739638
#> attr(,"logarithm")
#> [1] TRUE
#> 
#> $convergence
#> [1] 0
#> 
#> $iterations
#> [1] 41
#> 
#> $evaluations
#> function gradient 
#>       56       42 
#> 
#> $message
#> [1] "relative convergence (4)"

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