Browse Projects

This program (VBA implemented in Excel) animates the path of a tagged animal in a VR2 hydrohpone network.

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.

A series of modules for estimation of PVA parameters from time series data. Uses kalman filters, REML, and slope methods.

Miscellaneous matlab plotting functions that come in handy.

This program allows one to walk through the steps required to conduct a population viability analysis, or PVA, using a population time time series. The model outputs probabilities of extinction as a function of time steps into the future, and importantly, gives confidence intervals for these probabilities.

This tool has two major advantages over traditional PVA techniques:

1) It uses a state-space Kalman filter that allows for both process and non-process error.

So what's the big deal? --Functionally it filters the data, and allows a more accurate fit for population parameters of interest.
2) It uses a Bayesian sampling-importance-resampling algorithm to fully address uncertainty in the parameter estimates given the data.
So what's the big deal? -- Rather than developing a single function that describes the probability of population extinction through time, we can use the uncertainty in parameter estimates to develop 'probabilities of probabilities', or, the uncertainty surrounding the probability of extinction through time.

These are function for generating various standard types of stochastic population processes: random walks, Ornstein-Uhlenbeck, discrete Gompertz, etc. Also some random number generators needed by these are here. These matlab files need the Matlat Statistics Toolbox.

This is some matlab labs for teaching basic concepts about stochastic population trajectories and estimation for those trajectories.

Before you begin with this code, make sure you download the latest release of OpenBUGS(now 3.0.3) and it is installed to the default directory, "C:\Program Files\OpenBUGS". We've also found that it helps to have a cleaned up registry + defregmented hard drive before using OpenBUGS. The version of R that you're using should also be noted - R 2.5.1 and 2.7.0 work great with OpenBUGS, but there was a glitch in R 2.6.0 that prevented OpenBUGS from ever running through R. When using BUGS, I think that you can obviously use WinBUGS or OpenBUGS. While both programs should yield the same result, I had much more success with OpenBUGS, particularly in terms of the program not crashing during the burn-in (anyone familiar with 'trap' error messages knows what I'm talking about).

There are 2 project files. The first (writeModel.r) is sourced by the second file (runDataCloning2.r). If you want to modify the bounds on parameters, you'll need to edit writeModel.r - otherwise, everything will be done behind the scenes. I've tried to predict many of the errors that might occur - for example, it doesn't make sense to have all sites belonging to the same group (m = 1) and an unconstrained covariance matrix.

These files won't automate every possible MSSM that you want to fit, but they will do the majority of them. The R script files actually write the BUGS code for you, so you don't have to know any BUGS coding. You should however, be aware of the priors. Uniform priors are used on SDs and growth rates - to change these, look anywhere in the file for (~ dunif(a,b), where a and b are the limits). The code also includes the option of including an interaction matrix (multiple species) or density dependence matrix (multiple populations). The priors for all of these terms are bounded (0,1).

Included are several examples of summarizing parameters, including the median, mean, and using density to estimate the mode. The actual MLEs are going to be the points associated with the best (lowest) density.

These routines allow you to take a matrix of MCMC samples and calculate the Bayes factor based on the harmonic mean algorithm proposed by Gelfand and Dey (1994). Caution: Bayes factors tend to be numerically unstable!

This is project involves estimating catastrophes in the presence of observation and process error. It was later expanded into a paper on fur seals (Ward et al. 2007, CJFAS). The model may be run in WinBugs / OpenBugs, or through R. The example has been made to be generic and includes some dummy data. The priors included have been simplified for demonstration purposes.

These files are a collection of Java functions I've written to calculate PDFs/log likelihoods, and to do a variety of random number generator. You can use the built in Java random number generator, but I chose to link these files to the free RngPack library.

This is a collection of routines for numerical sampling and searching of multivariate surfaces. It contains things like gibbs samplers, metropolis algorithms, mcmc and simulated annealing algorithms.

LAMBDA is a MatLab toolkit designed to do MAR-1 based data analysis on long-term datasets and is based on the methods described in Ives et al. 2003, Ecological Monographs 73:301-330. LAMBDA is designed to allow the user to step through the entire modeling process, from importing the data, to obtaining descriptive statistics of the dataset, to, finally, performing a MAR-1 regression model and obtaining output parameters. A MAR-1 process is a Multivariate, Auto-Regressive first (1st) order process. Essentially, it is a means of estimating interactions between multiple variates from time series data, using matrix algebra. A MAR-1 model is a stochastic, non-mechanistic model that uses time series data on species numbers and co-variates to deduce inter-population interactions and the effects of covariates (e.g., physical variables) on populations.

Where does it come from?

LAMBDA is a product of the Mathematical Biology program at the Northwest Fisheries Science Center in Seattle, WA, and was developed with support by NOAA/NMFS and the National Research Council. It is open-source software released under the GNU GPL license, meaning you are free to use and modify it in (almost) any way you see fit. LAMBDA was developed by Steven Viscido while on a National Research Council postdoctoral associateship with Elizabeth Holmes.

Credits

LAMBDA is based on the techniques outlined in the paper Ives et al. 2003, Ecological Monographs 73:301-330. The code for the actual MAR-1 regression was based on the "MARbasic.m" MatLab script written by Tony Ives (available at the Ecological Archives). The parameter search code was based on an unpublished script written by Tony Ives. All other code was written by Steven Viscido.

Executable version

Download the LAMBDA executable along the installation instructions from the links below.

LAMBDA_MCR_pkg.exe 0.9.2 Warning: This is a 138 MB file!

Installation Instructions

The executable version of LAMBDA does not require MatLab. Its system requirements are Windows XP/2000, 256 MB of RAM, and 150 MB of Hard Drive Space.

If Installation hangs

This is a bug on MatLab's side. To work around it, you will need to install vcredist_x86.exe (32 bit systems) or vcredist_x64.exe (64 bit system) first and then repeat the LAMBDA installation. But read the If_Installation_Hangs.txt readme file if this happens to you. This bug affects about 20% of computers, randomly it would seem.

Matlab Source Code

This is not needed if you are using the executable version. Current source code version is LAMBDA_0.9.2Beta.zip. The source code can be downloaded below. Requirements for running LAMBDA from the source code are

• MatLab version 7.0.1 (R14) w/service Pack 1, or later
• MatLab's Statistics and Optimization toolboxes
• At least 256 MB of RAM on your system
• 5 MB of Hard Drive space for the LAMBDA installation

These are a series of Matlab routines to carry out a hierarchical linear model of fish population trajectories that includes a linear adjustment for a "reserve effect" at sites designated as MPAs. %The model is set up as follows: %FIRST LEVEL------------------------------------------------------------------------------------------------------- %MODEL:***** y~ intercept + abundance_index * (S_g + R_t +R_nt) + error *** %where: %S_g is the slope for a given species in a given geographic region (drawn %from a normally distributed species specific hyperparameter) %R_t is a species specific additive adjustment to the slope that is drawn %from a normally distributed hyper parameter representing the reserve effect on harvested species. Activated by a dummy variable. %R_nt is the species specific additive adjustment to the slope that is drawn from a normally distributed hyper parameter %representing the reserve effect on non-harvested species. Activated by a dummy variable. %both intercept and error terms are nuisance parameters and do not need hyperparameters. I suppose you could put in an %error hyper parameter with an informative prior just to keep the model in check. %SECOND LEVEL (hyper parameters all normally distributed)----------------------------------------------------------- %S_g_hat -- These are slopes for a given species from which S_g's are drawn %R_t_hat -- This is the reserve effect for targeted species from which R_t's across targeted species are drawn %R_nt_hat -- This is the reserve effect for non-targed species from which %R_nt's across non-targeted species are drawn %--------------------------------------------------------------------------------------------------------------------------

These files give code for ecological null model analyses.

This script writes and runs the data cloning procedure described by Lele et al. (2007, Ecology Letters). The routine is in R, and requires the user to install OpenBUGS before running. The number of clones can (and should) be modified. For this simple example, 100 clones takes ~ 2 minutes on my 2.5 year old laptop (1.66Ghz, 1GB RAM). More details and diagnostics can be found in the multivariate cloning project.

This project is a simple example of how to estimate the number of species using the sequential broken stick (lognormal) model. References include MacArthur (1957), MacArthur and Wilson (1967), Sugihara (1980), Solow (1994). This method that a sample of individuals (m) is drawn from a community. The observed species in the sample (Nm) represent a sample of the true number of species (N). N may be estimated using ML or Bayesian techniques by integrating out uncertainty in the array of relative abundances (p).

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)

This code is written to match the MCMC package SIAR in R. Two versions are included, 1 with no residual error, and 1 that models the variation as a mixture of residual + source variance.

This code is written to match the MCMC package SIAR in R. Two versions are included, 1 with no residual error, and 1 that models the variation as a mixture of residual + source variance.

This is just the list of animals used for the Ward et al. paper in Frontiers in Zoology.