R/var_to_cholvar.R
var_to_cholvar.RdThe free, fixed and par (and start) return the vec of the var-cov matrix M: fixed + free%*%p = vec(M) This function transforms the free, fixed, par and start so that they return the chol of the varcov matrix: vec(chol(M)).
var_to_cholvar(MLEobj)A properly formatted MARSS model as output by MARSS()
A new MARSS::marssMLE object with new fixed, free, par and start
Why does this function solve for the new par rather than just putting the chol values in? Because the user might have scrambled the order of the parameter list. I won't know which par correspond to which elements of the var-cov matrix. The free matrix is what does the sorting/permutation of the estimated parameters into their correct places. Also the matrix might be block-diagonal, partially fixed, etc.
Restrictions
Q and R cannot be time-varying (at the moment)
dat <- t(harborSealWA)[2:4, ]
fit <- MARSS(dat, model=list(Q="unconstrained"))
#> Success! abstol and log-log tests passed at 269 iterations.
#> Alert: conv.test.slope.tol is 0.5.
#> Test with smaller values (<0.1) to ensure convergence.
#>
#> MARSS fit is
#> Estimation method: kem
#> Convergence test: conv.test.slope.tol = 0.5, abstol = 0.001
#> Estimation converged in 269 iterations.
#> Log-likelihood: 30.36813
#> AIC: -34.73625 AICc: -25.40292
#>
#> Estimate
#> R.diag 0.01046
#> U.X.SJF 0.06694
#> U.X.SJI 0.07632
#> U.X.EBays 0.03760
#> Q.(1,1) 0.01283
#> Q.(2,1) 0.01259
#> Q.(3,1) 0.01128
#> Q.(2,2) 0.01238
#> Q.(3,2) 0.01108
#> Q.(3,3) 0.00994
#> x0.X.SJF 5.97015
#> x0.X.SJI 6.65086
#> x0.X.EBays 6.67692
#> Initial states (x0) defined at t=0
#>
#> Standard errors have not been calculated.
#> Use MARSSparamCIs to compute CIs and bias estimates.
#>
fitchol <- var_to_cholvar(fit)
Qchol = coef(fitchol, type="matrix")$Q
Q = coef(fit, type="matrix")$Q
Q
#> X.SJF X.SJI X.EBays
#> X.SJF 0.01283399 0.01259311 0.011283045
#> X.SJI 0.01259311 0.01237582 0.011082500
#> X.EBays 0.01128304 0.01108250 0.009935782
crossprod(Qchol)
#> X.SJF X.SJI X.EBays
#> X.SJF 0.01283399 0.01259311 0.011283045
#> X.SJI 0.01259311 0.01237582 0.011082500
#> X.EBays 0.01128304 0.01108250 0.009935782
Q = coef(fit, type="matrix", what="start")$Q
Qchol = coef(fitchol, type="matrix", what="start")$Q
Q
#> X.SJF X.SJI X.EBays
#> X.SJF 0.05 0.00 0.00
#> X.SJI 0.00 0.05 0.00
#> X.EBays 0.00 0.00 0.05
crossprod(Qchol)
#> X.SJF X.SJI X.EBays
#> X.SJF 0.05 0.00 0.00
#> X.SJI 0.00 0.05 0.00
#> X.EBays 0.00 0.00 0.05
fit <- MARSS(dat, model=list(Q=diag(0.1,3)+0.01), control=list(maxit=15))
#> Warning! Reached maxit before parameters converged. Maxit was 15.
#> neither abstol nor log-log convergence tests were passed.
#>
#> MARSS fit is
#> Estimation method: kem
#> Convergence test: conv.test.slope.tol = 0.5, abstol = 0.001
#> WARNING: maxit reached at 15 iter before convergence.
#> Neither abstol nor log-log convergence test were passed.
#> The likelihood and params are not at the MLE values.
#> Try setting control$maxit higher.
#> Log-likelihood: -2.605166
#> AIC: 19.21033 AICc: 21.69922
#>
#> Estimate
#> R.diag 0.00743
#> U.X.SJF 0.06828
#> U.X.SJI 0.06911
#> U.X.EBays 0.04302
#> x0.X.SJF 5.96983
#> x0.X.SJI 6.68251
#> x0.X.EBays 6.58815
#> Initial states (x0) defined at t=0
#>
#> Standard errors have not been calculated.
#> Use MARSSparamCIs to compute CIs and bias estimates.
#>
#> Convergence warnings
#> Warning: the R.diag parameter value has not converged.
#> Warning: the logLik parameter value has not converged.
#> Type MARSSinfo("convergence") for more info on this warning.
#>
fitchol <- var_to_cholvar(fit)
Qchol = coef(fitchol, type="matrix")$Q
Q = coef(fit, type="matrix")$Q
Q
#> X.SJF X.SJI X.EBays
#> X.SJF 0.11 0.01 0.01
#> X.SJI 0.01 0.11 0.01
#> X.EBays 0.01 0.01 0.11
crossprod(Qchol)
#> X.SJF X.SJI X.EBays
#> X.SJF 0.11 0.01 0.01
#> X.SJI 0.01 0.11 0.01
#> X.EBays 0.01 0.01 0.11
Q = coef(fit, type="matrix", what="start")$Q
Qchol = coef(fitchol, type="matrix", what="start")$Q
Q
#> X.SJF X.SJI X.EBays
#> X.SJF 0.11 0.01 0.01
#> X.SJI 0.01 0.11 0.01
#> X.EBays 0.01 0.01 0.11
crossprod(Qchol)
#> X.SJF X.SJI X.EBays
#> X.SJF 0.11 0.01 0.01
#> X.SJI 0.01 0.11 0.01
#> X.EBays 0.01 0.01 0.11