Fit a MARSS Model via Maximum-Likelihood Estimation

This is the main function for fitting multivariate autoregressive state-space (MARSS) models with linear constraints. Scroll down to the bottom to see some short examples. To open a guide to show you how to get started quickly, type RShowDoc("Quick_Start",package="MARSS"). To open the MARSS User Guide from the command line, type RShowDoc("UserGuide",package="MARSS"). To get an overview of the package and all its main functions and how to get output (parameter estimates, fitted values, residuals, Kalmin filter or smoother output, or plots), go to MARSS-package. If MARSS() is throwing errors or warnings that you don't understand, try the Troubleshooting section of the user guide or type MARSSinfo() at the command line.

The default MARSS model form is "marxss", which is Multivariate Auto-Regressive(1) eXogenous inputs State-Space model: $$\mathbf{x}_{t} = \mathbf{B}_t \mathbf{x}_{t-1} + \mathbf{u}_t + \mathbf{C}_t \mathbf{c}_t + \mathbf{G}_t \mathbf{w}_t, \textrm{ where } \mathbf{W}_t \sim \textrm{MVN}(0,\mathbf{Q}_t)$$ $$\mathbf{y}_t = \mathbf{Z}_t \mathbf{x}_t + \mathbf{a}_t + \mathbf{D}_t \mathbf{d}_t + \mathbf{H}_t \mathbf{v}_t, \textrm{ where } \mathbf{V}_t \sim \textrm{MVN}(0,\mathbf{R}_t)$$ $$\mathbf{X}_1 \sim \textrm{MVN}(\mathbf{x0}, \mathbf{V0}) \textrm{ or } \mathbf{X}_0 \sim \textrm{MVN}(\mathbf{x0}, \mathbf{V0}) $$ The parameters are everything except $\mathbf{x}$, $\mathbf{y}$, $\mathbf{v}$, $\mathbf{w}$, $\mathbf{c}$ and $\mathbf{d}$. $\mathbf{y}$ are data (missing values allowed). $\mathbf{c}$ and $\mathbf{d}$ are inputs (no missing values allowed). All parameters (except $\mathbf{x0}$ and $\mathbf{V0}$) can be time-varying but by default, all are time-constant (and the MARSS equation is generally written without the $t$ subscripts on the parameter matrices). All parameters can be zero, including the variance matrices.

The parameter matrices can have fixed values and linear constraints. This is an example of a 3x3 matrix with linear constraints. All matrix elements can be written as a linear function of $a$, $b$, and $c$: $$\left[\begin{array}{c c c} a+2b & 1 & a\\ 1+3a+b & 0 & b \\ 0 & -2 & c\end{array}\right]$$

Values such as $a b$ or $a^2$ or $log(a)$ are not linear constraints.

Usage

MARSS(y, 
    model = NULL, 
    inits = NULL, 
    miss.value = as.numeric(NA), 
    method = c("kem", "BFGS", "TMB", "BFGS_TMB", "nlminb_TMB"), 
    form = c("marxss", "dfa", "marss"), 
    fit = TRUE, 
    silent = FALSE, 
    control = NULL, 
    fun.kf = c("MARSSkfas", "MARSSkfss"), 
    ...)

Arguments

The default settings for the optional arguments are set in MARSSsettings.R and are given below in the details section. For form specific defaults see the form help file (e.g. MARSS.marxss or MARSS.dfa).

y

A n x T matrix of n time series over T time steps. Only y is required for the function. A ts object (univariate or multivariate) can be used and this will be converted to a matrix with time in the columns.

inits

A list with the same form as the list outputted by coef(fit) that specifies initial values for the parameters. See also MARSS.marxss.

model

Model specification using a list of parameter matrix text shortcuts or matrices. See Details and MARSS.marxss for the default form. Or better yet open the Quick Start Guide RShowDoc("Quick_Start",package="MARSS").

miss.value

Deprecated. Denote missing values by NAs in your data.

method

Estimation method. MARSS provides an EM algorithm (method="kem") (see MARSSkem) and the BFGS algorithm (method="BFGS") (see MARSSoptim).

form

The equation form used in the MARSS() call. The default is "marxss". See MARSS.marxss or MARSS.dfa.

fit

TRUE/FALSE Whether to fit the model to the data. If FALSE, a marssMLE object with only the model is returned.

silent

Setting to TRUE(1) suppresses printing of full error messages, warnings, progress bars and convergence information. Setting to FALSE(0) produces error output. Setting silent=2 will produce more verbose error messages and progress information.

fun.kf

What Kalman filter function to use. MARSS has two: MARSSkfas() which is based on the Kalman filter in the KFAS package based on Koopman and Durbin and MARSSkfss() which is a native R implementation of the Kalman filter and smoother in Shumway and Stoffer. The KFAS filter is much faster. MARSSkfas() modifies the input and output in order to output the lag-one covariance smoother needed for the EM algorithm (per page 321 in Shumway and Stoffer (2000).

control

Estimation options for the maximization algorithm. The typically used control options for method="kem" are below but see marssMLE for the full list of control options. Note many of these are not allowed if method="BFGS"; see MARSSoptim for the allowed control options for this method.

minit: The minimum number of iterations to do in the maximization routine (if needed by method). If method="kem", this is an easy way to up the iterations and see how your estimates are converging. (positive integer)
maxit: Maximum number of iterations to be used in the maximization routine (if needed by method) (positive integer).
min.iter.conv.test: Minimum iterations to run before testing convergence via the slope of the log parameter versus log iterations.
conv.test.deltaT=9: Number of iterations to use for the testing convergence via the slope of the log parameter versus log iterations.
conv.test.slope.tol: The slope of the log parameter versus log iteration to use as the cut-off for convergence. The default is 0.5 which is a bit high. For final analyses, this should be set lower. If you want to only use abstol as your convergence test, then to something very large, for example conv.test.slope.tol=1000. Type MARSSinfo(11) to see some comments on when you might want to do this.
abstol: The logLik.(iter-1)-logLik.(iter) convergence tolerance for the maximization routine. To meet convergence both the abstol and slope tests must be passed.
allow.degen: Whether to try setting $\mathbf{Q}$ or $\mathbf{R}$ elements to zero if they appear to be going to zero.
trace: An integer specifying the level of information recorded and error-checking run during the algorithms. trace=0, specifies basic error-checking and brief error-messages; trace>0 will print full error messages. In addition if trace>0, the Kalman filter output will be added to the outputted marssMLE object. Additional information recorded depends on the method of maximization. For the EM algorithm, a record of each parameter estimate for each EM iteration will be added. See optim for trace output details for the BFGS method. trace=-1 will turn off most internal error-checking and most error messages. The internal error checks are time expensive so this can speed up model fitting. This is particularly useful for bootstrapping and simulation studies. It is also useful if you get an error saying that MARSS() stops in MARSSkfss() due to a chol() call. MARSSkfss() uses matrix inversions and for some models these are unstable (high condition value). MARSSkfss() is used for error-checks and does not need to be called normally.
safe: Setting safe=TRUE runs the Kalman smoother after each parameter update rather than running the smoother only once after updated all parameters. The latter is faster but is not a strictly correct EM algorithm. In most cases, safe=FALSE (default) will not change the fits. If this setting does cause problems, you will know because you will see an error regarding the log-likelihood dropping and it will direct you to set safe=TRUE.

...

Optional arguments passed to function specified by form.

Details

The model argument specifies the structure of your model. There is a one-to-one correspondence between how you would write your model in matrix form on the whiteboard and how you specify the model for MARSS(). Many different types of multivariate time-series models can be converted to the MARSS form. See the User Guide and Quick Start Guide for examples.

The MARSS package has two forms for standard users: marxss and dfa.

MARSS.marxss: This is the default form. This is a MARSS model with (optional) inputs $\mathbf{c}_t$ or $\mathbf{d}_t$. Most users will want this help page.
MARSS.dfa: This is a model form to allow easier specification of models for Dynamic Factor Analysis. The $\mathbf{Z}$ parameters has a specific form and the $\mathbf{Q}$ is set at i.i.d (diagonal) with variance of 1.

Those looking to modify or understand the base code, should look at MARSS.marss and MARSS.vectorized. These describe the forms used by the base functions. The EM algorithm uses the MARSS model written in vectorized form. This form is what allows linear constraints.

The likelihood surface for MARSS models can be multimodal or with strong ridges. It is recommended that for final analyses the estimates are checked by using a Monte Carlo initial conditions search; see the chapter on initial conditions searches in the User Guide. This requires more computation time, but reduces the chance of the algorithm terminating at a local maximum and not reaching the true MLEs. Also it is wise to check the EM results against the BFGS results (if possible) if there are strong ridges in the likelihood. Such ridges seems to slow down the EM algorithm considerably and can cause the algorithm to report convergence far from the maximum-likelihood values. EM steps up the likelihood and the convergence test is based on the rate of change of the log-likelihood in each step. Once on a strong ridge, the steps can slow dramatically. You can force the algorithm to keep working by setting minit. BFGS seems less hindered by the ridges but can be prodigiously slow for some multivariate problems. BFGS tends to work better if you give it good initial conditions (see Examples below for how to do this).

If you are working with models with time-varying parameters, it is important to notice the time-index for the parameters in the process equation (the $\mathbf{x}$ equation). In some formulations (e.g. in KFAS), the process equation is $\mathbf{x}_t=\mathbf{B}_{t-1}\mathbf{x}_{t-1}+\mathbf{w}_{t-1}$ so $\mathbf{B}_{t-1}$ goes with $\mathbf{x}_t$ not $\mathbf{B}_t$. Thus one needs to be careful to line up the time indices when passing in time-varying parameters to MARSS(). See the User Guide for examples.

Value

An object of class marssMLE. The structure of this object is discussed below, but if you want to know how to get specific output (like residuals, coefficients, smoothed states, confidence intervals, etc), see print.marssMLE(), tidy.marssMLE(), MARSSresiduals() and plot.marssMLE().

The outputted marssMLE object has the following components:

model: MARSS model specification. It is a marssMODEL object in the form specified by the user in the MARSS() call. This is used by print functions so that the user sees the expected form.
marss: The marssMODEL object in marss form. This form is needed for all the internal algorithms, thus is a required part of a marssMLE object.
call: All the information passed in in the MARSS() call.
start: List with specifying initial values that were used for each parameter matrix.
control: A list of estimation options, as specified by arguments control.
method: Estimation method.

If fit=TRUE, the following are also added to the marssMLE object. If fit=FALSE, a marssMLE object ready for fitting via the specified method is returned.

par: A list of estimated parameter values in marss form. Use print(), tidy() or coef() for outputing the model estimates in the MARSS() call (e.g. the default "marxss" form).
states: The expected value of $\mathbf{X}$ conditioned on all the data, i.e. smoothed states.
states.se: The standard errors of the expected value of $\mathbf{X}$.
ytT: The expected value of $\mathbf{Y}$ conditioned on all the data. Note this is just $y$ for those $y$ that are not missing.
ytT.se: The standard errors of the expected value of $\mathbf{Y}$. Note this is 0 for any non-missing $y$.
numIter: Number of iterations required for convergence.
convergence: Convergence status. 0 means converged successfully, 3 means all parameters were fixed (so model did not need to be fit) and -1 means call was made with fit=FALSE and parameters were not fixed (thus no $par element and Kalman filter/smoother cannot be run). Anything else is a warning or error. 2 means the marssMLE object has an error; the object is returned so you can debug it. The other numbers are errors during fitting. The error code depends on the fitting method. See MARSSkem and MARSSoptim.
logLik: Log-likelihood.
AIC: Akaike's Information Criterion.
AICc: Sample size corrected AIC.

If control$trace is set to 1 or greater, the following are also added to the marssMLE object.

kf: A list containing Kalman filter/smoother output from MARSSkf(). This is not normally added to a marssMLE object since it is verbose, but can be added using MARSSkf().
Ey: A list containing output from MARSShatyt. This isn't normally added to a marssMLE object since it is verbose, but can be computed using MARSShatyt().

References

The MARSS User Guide: Holmes, E. E., E. J. Ward, and M. D. Scheuerell (2012) Analysis of multivariate time-series using the MARSS package. NOAA Fisheries, Northwest Fisheries Science Center, 2725 Montlake Blvd E., Seattle, WA 98112 Type RShowDoc("UserGuide",package="MARSS") to open a copy.

Holmes, E. E. (2012). Derivation of the EM algorithm for constrained and unconstrained multivariate autoregressive state-space (MARSS) models. Technical Report. arXiv:1302.3919 [stat.ME]

Holmes, E. E., E. J. Ward and K. Wills. (2012) MARSS: Multivariate autoregressive state-space models for analyzing time-series data. R Journal 4: 11-19.

Author

Eli Holmes, Eric Ward and Kellie Wills, NOAA, Seattle, USA.

Examples

dat <- t(harborSealWA)
dat <- dat[2:4, ] # remove the year row
# fit a model with 1 hidden state and 3 observation time series
kemfit <- MARSS(dat, model = list(
  Z = matrix(1, 3, 1),
  R = "diagonal and equal"
))
#> Success! abstol and log-log tests passed at 20 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 20 iterations. 
#> Log-likelihood: 15.70271 
#> AIC: -19.40541   AICc: -17.57932   
#>  
#>         Estimate
#> A.SJI    0.79691
#> A.EBays  0.27579
#> R.diag   0.02236
#> U.U      0.06075
#> Q.Q      0.00909
#> x0.x0    6.10491
#> Initial states (x0) defined at t=0
#> 
#> Standard errors have not been calculated. 
#> Use MARSSparamCIs to compute CIs and bias estimates.
#> 
kemfit$model # This gives a description of the model
#> 
#> Model form is marxss. Model Structure is
#> m: 1 state process(es) named X1
#> n: 3 observation time series named SJF SJI EBays
#> 
#> Z :  design matrix with rows: X1 X1 X1 (3 x 1) 
#> A :  see summary() (3 x 1) 
#> R :  diagonal and equal (3 x 3) 
#> B :  fixed and all one (1 x 1) 
#> U :  scalar (1 x 1) 
#> Q :  scalar (1 x 1) 
#> x0 :  scalar (1 x 1) 
#> V0 :  fixed and zero (1 x 1) 
#> D :  fixed and zero (3 x 1) 
#> C :  fixed and zero (1 x 1) 
#> d :  fixed and zero (1 x 1) 
#> c :  fixed and zero (1 x 1) 
#> G :  fixed and all one (1 x 1) 
#> H :  identity (3 x 3) 
#> L :  fixed and all one (1 x 1) 
print(kemfit$model) # same as kemfit$model
#> 
#> Model form is marxss. Model Structure is
#> m: 1 state process(es) named X1
#> n: 3 observation time series named SJF SJI EBays
#> 
#> Z :  design matrix with rows: X1 X1 X1 (3 x 1) 
#> A :  see summary() (3 x 1) 
#> R :  diagonal and equal (3 x 3) 
#> B :  fixed and all one (1 x 1) 
#> U :  scalar (1 x 1) 
#> Q :  scalar (1 x 1) 
#> x0 :  scalar (1 x 1) 
#> V0 :  fixed and zero (1 x 1) 
#> D :  fixed and zero (3 x 1) 
#> C :  fixed and zero (1 x 1) 
#> d :  fixed and zero (1 x 1) 
#> c :  fixed and zero (1 x 1) 
#> G :  fixed and all one (1 x 1) 
#> H :  identity (3 x 3) 
#> L :  fixed and all one (1 x 1) 
summary(kemfit$model) # This shows the model structure
#> Model Structure is
#> m: 1 state process(es)
#> n: 3 observation time series
#> $Z
#>       X1
#> SJF    1
#> SJI    1
#> EBays  1
#> 
#> $A
#>       Y      
#> SJF   0      
#> SJI   "SJI"  
#> EBays "EBays"
#> 
#> $R
#>       SJF    SJI    EBays 
#> SJF   "diag" 0      0     
#> SJI   0      "diag" 0     
#> EBays 0      0      "diag"
#> 
#> $B
#>    X1
#> X1  1
#> 
#> $U
#>    X  
#> X1 "U"
#> 
#> $Q
#>    X1 
#> X1 "Q"
#> 
#> $x0
#>    X   
#> X1 "x0"
#> 
#> $V0
#>    X1
#> X1  0
#> 
#> $D
#>      [,1]
#> [1,]    0
#> [2,]    0
#> [3,]    0
#> 
#> $C
#>      [,1]
#> [1,]    0
#> 
#> $d
#>      [,1]
#> [1,]    0
#> 
#> $c
#>      [,1]
#> [1,]    0
#> 
#> $G
#>      [,1]
#> [1,]    1
#> 
#> $H
#>      [,1] [,2] [,3]
#> [1,]    1    0    0
#> [2,]    0    1    0
#> [3,]    0    0    1
#> 
#> $L
#>      [,1]
#> [1,]    1
#> 
#> $tinitx
#> [1] 0
#> 

# add CIs to a marssMLE object
# default uses an estimated Hessian matrix
kem.with.hess.CIs <- MARSSparamCIs(kemfit)
kem.with.hess.CIs
#> 
#> MARSS fit is
#> Estimation method: kem 
#> Convergence test: conv.test.slope.tol = 0.5, abstol = 0.001
#> Estimation converged in 20 iterations. 
#> Log-likelihood: 15.70271 
#> AIC: -19.40541   AICc: -17.57932   
#>  
#>          ML.Est Std.Err   low.CI  up.CI
#> A.SJI   0.79691 0.04984  0.69922 0.8946
#> A.EBays 0.27579 0.05070  0.17642 0.3752
#> R.diag  0.02236 0.00512  0.01233 0.0324
#> U.U     0.06075 0.02145  0.01871 0.1028
#> Q.Q     0.00909 0.00569 -0.00206 0.0202
#> x0.x0   6.10491 0.13330  5.84365 6.3662
#> Initial states (x0) defined at t=0
#> 
#> CIs calculated at alpha = 0.05 via method=hessian 
#> 

# fit a model with 3 hidden states (default)
kemfit <- MARSS(dat, silent = TRUE) # suppress printing
kemfit
#> 
#> MARSS fit is
#> Estimation method: kem 
#> Convergence test: conv.test.slope.tol = 0.5, abstol = 0.001
#> Estimation converged in 44 iterations. 
#> Log-likelihood: 17.84491 
#> AIC: -15.68982   AICc: -10.45173   
#>  
#>                     Estimate
#> R.diag               0.00582
#> U.X.SJF              0.06833
#> U.X.SJI              0.07084
#> U.X.EBays            0.04221
#> Q.(X.SJF,X.SJF)      0.04150
#> Q.(X.SJI,X.SJI)      0.01271
#> Q.(X.EBays,X.EBays)  0.00807
#> x0.X.SJF             5.97602
#> x0.X.SJI             6.70656
#> x0.X.EBays           6.63306
#> Initial states (x0) defined at t=0
#> 
#> Standard errors have not been calculated. 
#> Use MARSSparamCIs to compute CIs and bias estimates.
#> 

# Fit the above model with BFGS using a short EM fit as initial conditions
kemfit <- MARSS(dat, control=list(minit=5, maxit=5))
#> Warning! Reached maxit before parameters converged. Maxit was 5.
#>  abstol not reached and no log-log test info since maxit less than min.iter.conv.test.
#> 
#> MARSS fit is
#> Estimation method: kem 
#> Convergence test: conv.test.slope.tol = 0.5, abstol = 0.001
#> WARNING: Abstol convergence only no info on log-log convergnece.
#>  No log-log convergence info because maxit (=5) < min.iter.conv.test (=15).
#>  The likelihood and params might not be at the ML values.
#>  Try setting control$maxit higher.
#> Log-likelihood: 15.95084 
#> 
#>  
#>                     Estimate
#> R.diag                0.0136
#> U.X.SJF               0.0683
#> U.X.SJI               0.0716
#> U.X.EBays             0.0419
#> Q.(X.SJF,X.SJF)       0.0221
#> Q.(X.SJI,X.SJI)       0.0145
#> Q.(X.EBays,X.EBays)   0.0121
#> x0.X.SJF              6.0171
#> x0.X.SJI              6.7335
#> x0.X.EBays            6.6519
#> Initial states (x0) defined at t=0
#> 
#> Standard errors have not been calculated. 
#> Use MARSSparamCIs to compute CIs and bias estimates.
#> 
#> Convergence warnings
#>  No convergence testing performed.
#> 
bffit <- MARSS(dat, method="BFGS", inits=kemfit)
#> Success! Converged in 17 iterations.
#> Function MARSSkfas used for likelihood calculation.
#> 
#> MARSS fit is
#> Estimation method: BFGS 
#> Estimation converged in 17 iterations. 
#> Log-likelihood: 17.852 
#> AIC: -15.704   AICc: -10.4659   
#>  
#>                     Estimate
#> R.diag               0.00539
#> U.X.SJF              0.06832
#> U.X.SJI              0.07069
#> U.X.EBays            0.04229
#> Q.(X.SJF,X.SJF)      0.04309
#> Q.(X.SJI,X.SJI)      0.01312
#> Q.(X.EBays,X.EBays)  0.00833
#> x0.X.SJF             5.97474
#> x0.X.SJI             6.70366
#> x0.X.EBays           6.62806
#> Initial states (x0) defined at t=0
#> 
#> Standard errors have not been calculated. 
#> Use MARSSparamCIs to compute CIs and bias estimates.
#> 

# fit a model with 3 correlated hidden states
# with one variance and one  covariance
# maxit set low to speed up example, but more iters are needed for convergence
kemfit <- MARSS(dat, model = list(Q = "equalvarcov"), control = list(maxit = 50))
#> Warning! Reached maxit before parameters converged. Maxit was 50.
#>  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  50  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: 28.72123 
#> AIC: -39.44246   AICc: -35.25641   
#>  
#>            Estimate
#> R.diag       0.0109
#> U.X.SJF      0.0669
#> U.X.SJI      0.0759
#> U.X.EBays    0.0386
#> Q.diag       0.0114
#> Q.offdiag    0.0113
#> x0.X.SJF     5.9930
#> x0.X.SJI     6.6645
#> x0.X.EBays   6.6372
#> 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  logLik  parameter value has not converged.
#>  Type MARSSinfo("convergence") for more info on this warning.
#> 
# use Q="unconstrained" to allow different variances and covariances

# fit a model with 3 independent hidden states
# where each observation time series is independent
# the hidden trajectories 2-3 share their U parameter
kemfit <- MARSS(dat, model = list(U = matrix(c("N", "S", "S"), 3, 1)))
#> Success! abstol and log-log tests passed at 45 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 45 iterations. 
#> Log-likelihood: 17.48233 
#> AIC: -16.96465   AICc: -12.77861   
#>  
#>                     Estimate
#> R.diag               0.00542
#> U.N                  0.06832
#> U.S                  0.05330
#> Q.(X.SJF,X.SJF)      0.04279
#> Q.(X.SJI,X.SJI)      0.01382
#> Q.(X.EBays,X.EBays)  0.00862
#> x0.X.SJF             5.97495
#> x0.X.SJI             6.72656
#> x0.X.EBays           6.61015
#> Initial states (x0) defined at t=0
#> 
#> Standard errors have not been calculated. 
#> Use MARSSparamCIs to compute CIs and bias estimates.
#> 

# same model, but with fixed independent observation errors
# and the 3rd x processes are forced to have a U=0
# Notice how a list matrix is used to combine fixed and estimated elements
# all parameters can be specified in this way using list matrices
kemfit <- MARSS(dat, model = list(U = matrix(list("N", "N", 0), 3, 1), R = diag(0.01, 3)))
#> Success! algorithm run for 15 iterations. abstol and log-log tests passed.
#> 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
#> Algorithm ran 15 (=minit) iterations and convergence was reached. 
#> Log-likelihood: 15.47987 
#> AIC: -16.95974   AICc: -14.47085   
#>  
#>                     Estimate
#> U.N                   0.0709
#> Q.(X.SJF,X.SJF)       0.0296
#> Q.(X.SJI,X.SJI)       0.0107
#> Q.(X.EBays,X.EBays)   0.0116
#> x0.X.SJF              5.9871
#> x0.X.SJI              6.7310
#> x0.X.EBays            6.7139
#> Initial states (x0) defined at t=0
#> 
#> Standard errors have not been calculated. 
#> Use MARSSparamCIs to compute CIs and bias estimates.
#> 

# fit a model with 2 hidden states (north and south)
# where observation time series 1-2 are north and 3 is south
# Make the hidden state process independent with same process var
# Make the observation errors different but independent
# Make the growth parameters (U) the same
# Create a Z matrix as a design matrix that assigns the "N" state to the first 2 rows of dat
# and the "S" state to the 3rd row of data
Z <- matrix(c(1, 1, 0, 0, 0, 1), 3, 2)
# You can use factor is a shortcut making the above design matrix for Z
# Z <- factor(c("N","N","S"))
# name the state vectors
colnames(Z) <- c("N", "S")
kemfit <- MARSS(dat, model = list(
  Z = Z,
  Q = "diagonal and equal", R = "diagonal and unequal", U = "equal"
))
#> Success! abstol and log-log tests passed at 32 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 32 iterations. 
#> Log-likelihood: 24.7002 
#> AIC: -33.4004   AICc: -30.12767   
#>  
#>                 Estimate
#> A.SJI            0.79810
#> R.(SJF,SJF)      0.02473
#> R.(SJI,SJI)      0.00581
#> R.(EBays,EBays)  0.00305
#> U.1              0.05641
#> Q.diag           0.01160
#> x0.N             5.94175
#> x0.S             6.58600
#> Initial states (x0) defined at t=0
#> 
#> Standard errors have not been calculated. 
#> Use MARSSparamCIs to compute CIs and bias estimates.
#> 

# print the model followed by the marssMLE object
kemfit$model
#> 
#> Model form is marxss. Model Structure is
#> m: 2 state process(es) named N S
#> n: 3 observation time series named SJF SJI EBays
#> 
#> Z :  design matrix with rows: N S N S (3 x 2) 
#> A :  see summary() (3 x 1) 
#> R :  diagonal and unequal (3 x 3) 
#> B :  identity (2 x 2) 
#> U :  all equal (2 x 1) 
#> Q :  diagonal and equal (2 x 2) 
#> x0 :  unconstrained (2 x 1) 
#> V0 :  fixed and zero (2 x 2) 
#> D :  fixed and zero (3 x 1) 
#> C :  fixed and zero (2 x 1) 
#> d :  fixed and zero (1 x 1) 
#> c :  fixed and zero (1 x 1) 
#> G :  identity (2 x 2) 
#> H :  identity (3 x 3) 
#> L :  identity (2 x 2) 

if (FALSE) {
# simulate some new data from our fitted model
sim.data <- MARSSsimulate(kemfit, nsim = 10, tSteps = 10)

# Compute bootstrap AIC for the model; this takes a long, long time
kemfit.with.AICb <- MARSSaic(kemfit, output = "AICbp")
kemfit.with.AICb
}

if (FALSE) {
# Many more short examples can be found in the
# Quick Examples chapter in the User Guide
RShowDoc("UserGuide", package = "MARSS")

# You can find the R scripts from the chapters by
# going to the index page
RShowDoc("index", package = "MARSS")
}