mcmcsae (version 0.8.1)

f_gaussian: Specify a Gaussian sampling distribution

Description

This function can be used in the family argument of create_sampler or generate_data to specify a Gaussian sampling distribution.

Usage

f_gaussian(
  link = "identity",
  var.prior = pr_invchisq(df = 0, scale = 1),
  var.vec = ~1,
  prec.mat = NULL,
  var.model = NULL,
  logJacobian = NULL,
  control = gaussian_control()
)

Value

A family object.

Arguments

link

the name of a link function. Currently the only allowed link function for a Gaussian distribution is "identity".

var.prior

prior for the variance parameter of a Gaussian sampling distribution. This can be specified by a call to one of the prior specification functions pr_invchisq, pr_exp, pr_gig or pr_fixed for inverse chi-squared, exponential, generalised inverse gaussian or degenerate prior distribution, respectively. The default is an improper prior pr_invchisq(df=0, scale=1). A half-t prior on the standard deviation can be specified using pr_invchisq with a chi-squared distributed scale parameter.

var.vec

a formula to specify unequal variances, i.e. heteroscedasticity. The default corresponds to equal variances.

prec.mat

a possibly non-diagonal positive-definite symmetric matrix interpreted as the precision matrix, i.e. inverse of the covariance matrix. If this argument is specified var.vec is ignored.

var.model

a formula specifying the terms of a variance model in the case of a Gaussian likelihood. Several types of terms are supported: a regression term for the log-variance specified with vreg(...), and a term vfac(...) for multiplicative modelled factors at a certain level specified by a factor variable. By using unit-level inverse-chi-squared factors the marginal sampling distribution becomes a Student-t distribution, and by using unit-level exponential factors it becomes a Laplace or double exponential distribution. In addition, reg and gen can be used to specify regression or random effect terms. In that case the prior distribution of the coefficients is not exactly normal, but instead Multivariate Log inverse Gamma (MLiG), see also pr_MLiG.

logJacobian

if the data are transformed the logarithm of the Jacobian can be supplied so that it is incorporated in all log-likelihood computations. This can be useful for comparing information criteria for different transformations. It should be supplied as a vector of the same size as the response variable. For example, when a log-transformation is used on response vector y, the vector -log(y) should be supplied.

control

a list with computational options. These options can be specified using function gaussian_control.

Examples

if (FALSE) {
n <- 4000
m <- 25
dat <- data.frame(
  x = rnorm(n),
  g = factor(sample(1:m, n, replace=TRUE), levels=1:m)
)
v <- rnorm(m, sd=0.6)
dat$y <- rnorm(n, mean = with(dat, 1 - 0.5*x + v[g]), sd=0.4)

sampler <- create_sampler(
  y ~ x + (1|g), data=dat
)
sim <- MCMCsim(sampler, store.all=TRUE)
compute_DIC(sim)
summary(sim)

# more explicit specification, allowing non-default names, priors etc.
sampler <- create_sampler(
  y ~ reg(~ x, name="beta") + gen(~1, factor = ~ g, name="v"),
  data=dat,
  family = f_gaussian(var.prior = pr_fixed(value = 0.4^2))
)
sim <- MCMCsim(sampler, store.all=TRUE)
compute_DIC(sim)
summary(sim)
bayesplot::mcmc_recover_intervals(as.array(sim$beta), c(1, -0.5))
bayesplot::mcmc_recover_hist(as.array(sim$v_sigma), 0.6)
bayesplot::mcmc_recover_scatter(as.array(sim$v), v)

# heteroscedastic data
dat$y <- rnorm(n,
  mean = with(dat, 1 - 0.5*x + v[g]),
  sd = with(dat, exp(0.5*(0.2 + 0.5*x)))
)

sampler <- create_sampler(
  y ~ reg(~ x, name="beta") + gen(~1, factor = ~ g, name="v"),
  data=dat,
  family = f_gaussian(
    var.prior = pr_fixed(value = 1),
    var.model = ~ 1 + x
  )
)
sim <- MCMCsim(sampler, store.all=TRUE)
compute_DIC(sim)
compute_WAIC(sim)
summary(sim)
bayesplot::mcmc_recover_intervals(as.array(sim$beta), c(1, -0.5))
bayesplot::mcmc_recover_intervals(as.array(sim$vreg1), c(0.2, 0.5))
bayesplot::mcmc_recover_hist(as.array(sim$v_sigma), 0.6)
bayesplot::mcmc_recover_scatter(as.array(sim$v), v)
}