mark.scheuerell

This set of Matlab scripts conducts a first order MAR(1) model to estimate interactive effects, and covariate effects, on a time series of data from a community.

It uses a Gibbs sampler to estimate parameters, and currently is set up with diffuse priors on all parameters for the model. It is pretty basic at this point, but it works.

Note on the Gamma distribution:

Going back and forth between Matlab/R/WinBUGS can be confusing, because of the different parameterizations of the gamma pdf. Here's the Matlab/BUGS forms and the R equivalents:

Matlab: X ~ g(a,b) E[X] = ab

R equivalent: X ~ g(shape=a,scale=b)

BUGS: X ~ g(a,b) E[X] = a/b

R equivalent: X ~ g(shape=a,rate=b) OR

X ~ g(shape=a,scale=1/b)

THE FOLLOWING FUNCTION CONDUCTS MULTIPLE LINEAR REGRESSION PARAMETER ESTIMATION VIA MCMC USING GIBBS SAMPLING.



The point of this code is to show how to do it with a conjugate distribution such that the method can easily be plugged into LAMBDA.



The function assumes that:

1)the B coefficients are normally distributed ~N(prior_mean,prior_s2)and thus have a normal distribution prior. Large prior variances produce suitably "diffuse" priors

2) the model variance is gamma distributed Ga(alpha, beta), with an associated prior. Small values of alpha and beta (<.001) produce a suitably "diffuse" prior. Note that Matlab's parameterization of their gamma function is a bit wonky compared to the standard. Basically, Matlab uses Ga(alpha,1/beta)-- ie. inverse of the standard scale parameter.

The R code contains 2 versions of Bayesian linear regression. The first (univariate.all) uses the multivariate normal distribution to sample regression coefficients; the second (univariate.vat) does variable at a time sampling. Both R versions are nearly identical, allow a flexible range of priors for B, and allow the prior variance to be specified either as the shape-scale (e.g. gamma, scaled inv chi sq) or in terms of the mean-variance. Output includes saved parameter draws, which can be fed into CODA, and the calculation of the log(Bayes Factor) to be used for model selection.

This is the DEVELOPMENT site for the MARSS.  For the current MARSS  release go to CRAN or download straight from the R GUI using "Install Packages" menu.

MARSS fits mulitvariate autoregressive state-space (MARSS) models with Gaussian errors to multivariate time series data.  A MARSS model is:

x(t) = B x(t-1) + u + v(t), v(t)~MVN(0,Q)

y(t) = Z x(t) + a + w(t), w(t)~MVN(0,R)

Project news (Jan 27, 2012): MARSS 2.8 uploaded to CRAN.  2.3+ implements the algorithm for a fully unconstrained MARSS model with fixed and shared elements in all parameters.  See changes.pdf for the fixes since 2.5.  MARSS 2.9 is in the works which will allow more typical specification of covariates. I'm also working on MARSS 3.0, which implements the general EM algorithm with linear constraints on the parameter matrices. MARSS 2.8 is implementing a constrained version of the more general algorithm in the EMDerivation paper.  The change for 3.0 requires changing the 'wrappers' so user can specify the linear constraints.  The actual algorithm code (MARSSkem.r) is unchanged.

Developers: Eli Holmes, Eric Ward and Kellie Wills

Current known issues:

• lap-p models not working with method=kem since kemcheck is blocking.  Need to review EM algorithm per constraint that B subblock corresponding to diag(Q)=0 must be diagonal.    Use method="BFGS" as a work around. 
• moving average models not tested.
• The likelihood for the covariate "trick" where you want covariate to affect process only and do that by setting R=zero is wrong.  There shouldn't be a likelihood of the covariate added in, but will be since Q=1.I added a note to manual for 2.8 saying you need to subtract off that extra LL.
  
• For cross-platform compatibility with MacOS, I think all of the source files need a hard return / empty line as the last line.