dsdive.gibbs.obs.cov: Gibbs sampler for parameters of a model for dives across...
In jmhewitt/dsdive: Discrete Space Models for Dive Trajectories of Animals

This function differs from dsdive.gibbs.obs in that this function allows the model to be estimated with dive-specific covariates.

dsdive.gibbs.obs.cov(
  dsobs.list,
  covs,
  t.stages.list,
  beta.init,
  alpha.init,
  verbose = FALSE,
  maxit,
  checkpoint.fn,
  checkpoint.interval = 3600,
  beta1.prior,
  beta2.prior,
  alpha1.prior,
  alpha2.prior,
  alpha3.prior,
  tstep,
  depth.bins,
  T1.prior.params,
  T2.prior.params,
  max.width,
  max.width.offset,
  t0.prior.params,
  tf.prior.params,
  offsets,
  offsets.tf,
  cl,
  pi.formula,
  lambda.formula,
  warmup = Inf,
  delta = 1e-10,
  optim.maxit = 1000,
  adaptive = FALSE,
  adaptation.frequency = 10,
  gapprox = TRUE
)

`dsobs.list`	list of `dsobs` objects, which describe the observation times and depths of a collection of dives
`covs`	matrix of covariates associated with `dsobs.list`. Each row ov `covs` should contain all covariates for a single dive.
`t.stages.list`	list of initial stage transition times for dives observed in `dsobs.list`
`beta.init`	Initial values for directional preference model parameters. See `dsdive.tx.params` for more details. Should be a list, in which each component contains the initial value for coefficients for logit(pi^(s).)
`alpha.init`	Initial values for diving rate model parameters. See `dsdive.tx.params` for more details. Should be a list, in which each component contains the initial value for coefficients for log(lambda^(s)).
`verbose`	`TRUE` to output sampler status while running
`maxit`	number of Gibbs iterations to run
`checkpoint.fn`	User-defined function to run during a checkpoint step; gives the user an opportunity to save partial output from the sampler
`checkpoint.interval`	Number of seconds between calls to `checkpoint.fn`
`beta1.prior`	List containing `mean` and `sd` parameters for independent Normal prior distributions on the coefficients for the dive-stage model parameter π^{(1)}. The length of the `mean` and `sd` vectors should be equal to the number of coefficients for this model component.
`beta2.prior`	List containing `mean` and `sd` parameters for independent Normal prior distributions on the coefficients for the dive-stage model parameter π^{(3)}. The length of the `mean` and `sd` vectors should be equal to the number of coefficients for this model component.
`alpha1.prior`	List containing `mean` and `sd` parameters for independent Normal prior distributions on the coefficients for the dive-stage model parameter λ^{(1)}. The length of the `mean` and `sd` vectors should be equal to the number of coefficients for this model component.
`alpha2.prior`	List containing `mean` and `sd` parameters for independent Normal prior distributions on the coefficients for the dive-stage model parameter λ^{(2)}. The length of the `mean` and `sd` vectors should be equal to the number of coefficients for this model component.
`alpha3.prior`	List containing `mean` and `sd` parameters for independent Normal prior distributions on the coefficients for the dive-stage model parameter λ^{(3)}. The length of the `mean` and `sd` vectors should be equal to the number of coefficients for this model component.
`tstep`	Time between observations in `dsobs.list`
`depth.bins`	n x 2 Matrix that defines the depth bins. The first column defines the depth at the center of each depth bin, and the second column defines the half-width of each bin.
`T1.prior.params`	`shape` and `rate` parameters for Gamma prior on the descent-stage duration.
`T2.prior.params`	`shape` and `rate` parameters for Gamma prior on the bottom-stage duration.
`max.width`	The stage transition times are updated with a piecewise proposal distribution. `max.width` controls the maximum width of the intervals for the proposal distribution. This is a tuning parameter that controls the numerical stability of the proposal distribution, which is sampled via inverse CDF techniques.
`max.width.offset`	The t0 and tf offsets are updated with a piecewise proposal distribution. `max.width.offset` controls the maximum width of the intervals for the proposal distribution. This is a tuning parameter that controls the numerical stability of the proposal distribution, which is sampled via inverse CDF techniques.
`t0.prior.params`	`shape1` and `shape2` parameters for the scaled and shifted Beta prior distribution for the t0 offset.
`tf.prior.params`	`shape1` and `shape2` parameters for the scaled and shifted Beta prior distribution for the tf offset.
`offsets`	vector with initial values for t0 offsets.
`offsets.tf`	vector with initial values for tf offsets.
`cl`	Shared-memory cluster to be used to distribute some computations. The cluster must be fully initialized before running this function. The cluster requires random seeds and `Rdsm`-initialization.
`pi.formula`	List of formula objects, each of which defines the linear model that combines covariates in the `covs` matrix into dive-specific diving preferences logit(pi^(s)).
`lambda.formula`	List of formula objects, each of which defines the linear model that combines covariates in the `covs` matrix into dive-specific diving rates log(lambda^(s)).
`warmup`	number of iterations during which the proposal distributions will be updated at each step
`delta`	If `delta>0`, then the observation matrix and raw components computed will be for a transition matrix whose generator is perturbed to allow much faster computation.
`optim.maxit`	maximum number of steps to take during numerical optimization to compute Gaussian approximation to full conditional posteriors used to propose model parameters
`adaptive`	`TRUE` to use adaptive Random walk Metropolis-Hastings samplers to update model parameters instead of Gaussian approximations to the full conditional posteriors.
`adaptation.frequency`	Random walk proposals will only be updated at intervals of this step count
`gapprox`	If `NULL`, then exact-likelihood Gaussian approximations to the full conditional posteriors will be used to update model parameters. Otherwise, `gapprox` should be a list that defines the interpolation grids used to construct approximate likelihoods, for building Gaussian approximations. If `gapprox==TRUE`, then default approximation settings will be used

data('dive.sim')
attach(dive.sim)
attach(dive.sim$params)
library(parallel)
library(Rdsm)

cl = makeCluster(2, 'SOCK')
clusterEvalQ(cl, library(dsdive))
clusterEvalQ(cl, library(Rdsm))

mgrinit(cl)

t.stages = sim$times[c(FALSE,diff(sim$stages)==1)]
  
tstep = diff(sim.obs$times[1:2])

obstx.mat = lapply(1:3, function(s) {
  dsdive.obstx.matrix(depth.bins = depth.bins, beta = beta, 
                      lambda = lambda, s0 = s, tstep = tstep, 
                      include.raw = TRUE, delta = 1e-10)
})

lambda.priors = list(
  list(mu = 0, sd = 3),
  list(mu = 0, sd = 3),
  list(mu = 0, sd = 3)
)

beta.priors = list(
  list(mu = 0, sd = 3),
  list(mu = 0, sd = 3)
)

T1.prior.params = c(25, .04)
T2.prior.params = c(56, .06)

dsobs.list = list(sim.obs, sim.obs)
t.stages.list = list(t.stages, t.stages)

beta.init = list(
  c(qlogis(.5), 0),
  c(0, 1)
)

alpha.init = c(beta.init, list(c(-1,-1)))

covs = data.frame(x1 = c(.5, 1), x2 = c(0, .3))

pi.formula = ~x1
lambda.formula = ~x1:x2

fit = dsdive.gibbs.obs.cov(
  dsobs.list = dsobs.list, t.stages.list = t.stages.list, 
  beta.init = beta.init, alpha.init = alpha.init, verbose = TRUE, maxit = 1, 
  beta1.prior = beta.priors[[1]], beta2.prior = beta.priors[[2]], 
  alpha1.prior = lambda.priors[[1]], alpha2.prior = lambda.priors[[2]], 
  alpha3.prior = lambda.priors[[3]], tstep = tstep, depth.bins = depth.bins, 
  T1.prior.params = T1.prior.params, T2.prior.params = T2.prior.params, 
  max.width = 100, max.width.offset = 30, t0.prior.params = c(1,1), 
  tf.prior.params = c(1,1), offsets = 0, offsets.tf = 0, warmup = 1, 
  covs = covs, pi.formula = pi.formula,  lambda.formula = lambda.formula, 
  cl = cl, optim.maxit = 1, delta = 1e-10)


detach(dive.sim$params)
detach(dive.sim)

stopCluster(cl)