Skip to contents

Parameter Variance-Covariance Matrix from the Hessian Matrix

MARSShessian.Rd

Calculates an approximate parameter variance-covariance matrix for the parameters using an inverse of the Hessian of the negative log-likelihood function at the MLEs (the observed Fisher Information matrix). It appends $Hessian, $parMean, $parSigma to the marssMLE object.

Usage

MARSShessian(MLEobj, method=c("Harvey1989", "fdHess", "optim"))

Arguments

MLEobj

An object of class marssMLE. This object must have a $par element containing MLE parameter estimates from e.g. MARSSkem.

method

The method to use for computing the Hessian. Options are Harvey1989 to use the Harvey (1989) recursion, which is an analytical solution, fdHess or optim which are two numerical methods. Although optim can be passed to this function, in the internal functions which call this function, fdHess will be used if a numerical estimate is requested.

Details

See MARSSFisherI for a discussion of the observed Fisher Information matrix and references.

Method fdHess uses fdHess from package nlme to numerically estimate the Hessian matrix (the matrix of partial 2nd derivatives of the negative log-likelihood function at the MLE). Method optim uses optim with hessian=TRUE and list(maxit=0) to ensure that the Hessian is computed at the values in the par element of the MLE object. Method Harvey1989 (the default) uses the recursion in Harvey (1989) to compute the observed Fisher Information of a MARSS model analytically.

Note that the parameter confidence intervals computed with the observed Fisher Information matrix are based on the asymptotic normality of maximum-likelihood estimates under a large-sample approximation.

Value

MARSShessian() attaches

Hessian, parMean and parSigma to the marssMLE object that is passed into the function.

Author

Eli Holmes, NOAA, Seattle, USA.

See also

MARSSFisherI(), MARSSharveyobsFI(), MARSShessian.numerical(), MARSSparamCIs(), marssMLE

References

Harvey, A. C. (1989) Section 3.4.5 (Information matrix) in Forecasting, structural time series models and the Kalman filter. Cambridge University Press, Cambridge, UK.

See also J. E. Cavanaugh and R. H. Shumway (1996) On computing the expected Fisher information matrix for state-space model parameters. Statistics & Probability Letters 26: 347-355. This paper discusses the Harvey (1989) recursion (and proposes an alternative).

Examples

dat <- t(harborSeal)
dat <- dat[c(2, 11), ]
MLEobj <- MARSS(dat)
#> Success! abstol and log-log tests passed at 26 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 26 iterations. 
#> Log-likelihood: 11.74016 
#> AIC: -9.480311   AICc: -6.3692   
#>  
#>                                           Estimate
#> R.diag                                      0.0115
#> U.X.CoastalEstuaries                        0.0613
#> U.X.OR.NorthCoast                           0.0510
#> Q.(X.CoastalEstuaries,X.CoastalEstuaries)   0.0147
#> Q.(X.OR.NorthCoast,X.OR.NorthCoast)         0.0122
#> x0.X.CoastalEstuaries                       7.3823
#> x0.X.OR.NorthCoast                          6.2707
#> Initial states (x0) defined at t=0
#> 
#> Standard errors have not been calculated. 
#> Use MARSSparamCIs to compute CIs and bias estimates.
#> 
MLEobj.hessian <- MARSShessian(MLEobj)

# show the approx Hessian
MLEobj.hessian$Hessian
#>                                                  R.diag U.X.CoastalEstuaries
#> R.diag                                    52914.2357125           -154.19136
#> U.X.CoastalEstuaries                       -154.1913603           1662.00142
#> U.X.OR.NorthCoast                           169.2513722              0.00000
#> Q.(X.CoastalEstuaries,X.CoastalEstuaries) 12164.4113628            120.60426
#> Q.(X.OR.NorthCoast,X.OR.NorthCoast)       13737.8120129              0.00000
#> x0.X.CoastalEstuaries                         1.9507559             68.10553
#> x0.X.OR.NorthCoast                           -0.7477657              0.00000
#>                                           U.X.OR.NorthCoast
#> R.diag                                             169.2514
#> U.X.CoastalEstuaries                                 0.0000
#> U.X.OR.NorthCoast                                 2332.1767
#> Q.(X.CoastalEstuaries,X.CoastalEstuaries)            0.0000
#> Q.(X.OR.NorthCoast,X.OR.NorthCoast)               -159.5791
#> x0.X.CoastalEstuaries                                0.0000
#> x0.X.OR.NorthCoast                                  82.0963
#>                                           Q.(X.CoastalEstuaries,X.CoastalEstuaries)
#> R.diag                                                                 12164.411363
#> U.X.CoastalEstuaries                                                     120.604263
#> U.X.OR.NorthCoast                                                          0.000000
#> Q.(X.CoastalEstuaries,X.CoastalEstuaries)                              16821.992564
#> Q.(X.OR.NorthCoast,X.OR.NorthCoast)                                        0.000000
#> x0.X.CoastalEstuaries                                                     -1.525828
#> x0.X.OR.NorthCoast                                                         0.000000
#>                                           Q.(X.OR.NorthCoast,X.OR.NorthCoast)
#> R.diag                                                          13737.8120129
#> U.X.CoastalEstuaries                                                0.0000000
#> U.X.OR.NorthCoast                                                -159.5790963
#> Q.(X.CoastalEstuaries,X.CoastalEstuaries)                           0.0000000
#> Q.(X.OR.NorthCoast,X.OR.NorthCoast)                             23236.9701088
#> x0.X.CoastalEstuaries                                               0.0000000
#> x0.X.OR.NorthCoast                                                  0.7050328
#>                                           x0.X.CoastalEstuaries
#> R.diag                                                 1.950756
#> U.X.CoastalEstuaries                                  68.105533
#> U.X.OR.NorthCoast                                      0.000000
#> Q.(X.CoastalEstuaries,X.CoastalEstuaries)             -1.525828
#> Q.(X.OR.NorthCoast,X.OR.NorthCoast)                    0.000000
#> x0.X.CoastalEstuaries                                 44.922134
#> x0.X.OR.NorthCoast                                     0.000000
#>                                           x0.X.OR.NorthCoast
#> R.diag                                            -0.7477657
#> U.X.CoastalEstuaries                               0.0000000
#> U.X.OR.NorthCoast                                 82.0962990
#> Q.(X.CoastalEstuaries,X.CoastalEstuaries)          0.0000000
#> Q.(X.OR.NorthCoast,X.OR.NorthCoast)                0.7050328
#> x0.X.CoastalEstuaries                              0.0000000
#> x0.X.OR.NorthCoast                                21.5655793

# generate a parameter sample using the Hessian
# this uses the rmvnorm function in the mvtnorm package
hess.params <- mvtnorm::rmvnorm(1,
  mean = MLEobj.hessian$parMean,
  sigma = MLEobj.hessian$parSigma
)