Markov Chain Monte Carlo Small Area Estimation

Description

Fit multi-level models with possibly correlated random effects using MCMC.

Details

Functions to fit multi-level models with Gaussian, binomial, multinomial, negative binomial or Poisson likelihoods using MCMC. Models with a linear predictor consisting of various possibly correlated random effects are supported, allowing flexible modelling of temporal, spatial or other kinds of dependence structures. For Gaussian models the variance can be modelled too. By modelling variances at the unit level the marginal distribution can be changed to a Student-t or Laplace distribution, which may account better for outliers. The package has been developed with applications to small area estimation in official statistics in mind. The posterior samples for the model parameters can be passed to a prediction function to generate samples from the posterior predictive distribution for user-defined quantities such as finite population domain means. For model assessment, posterior predictive checks and DIC/WAIC criteria can easily be computed.

Solve Ax=b by preconditioned conjugate gradients

Description

Solve Ax=b by preconditioned conjugate gradients

Usage

CG(
  b,
  env,
  x = 0 * b,
  max.it = length(b),
  e = 1e+06 * length(b),
  verbose = FALSE,
  ...
)

Arguments

Value

The (approximated) solution to Ax=b.

References

M.R. Hestenes and E. Stiefel (1952). Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards 49(6), 409-436.

Set options for the conjugate gradient (CG) sampler

Description

Set options for the conjugate gradient (CG) sampler

Usage

CG_control(
  max.it = NULL,
  stop.criterion = NULL,
  preconditioner = c("GMRF", "GMRF2", "GMRF3", "identity"),
  scale = 1,
  chol.control = chol_control(),
  verbose = FALSE
)

Arguments

Value

A list of options used by the conjugate gradients algorithm.

Set up a GMRF structure for a generic model component

Description

This function is used to specify a (non-default) GMRF structure to pass to argument strucA of function gen.

Usage

GMRF_structure(
  type = c("default", "bym2", "leroux", "Leroux"),
  scale.precision = (type == "bym2"),
  prior = NULL,
  control = NULL
)

Arguments

Value

An environment defining the desired GMRF structure, for use by other package functions.

References

B. Leroux, X. Lei and N. Breslow (1999). Estimation of Disease Rates in Small Areas: A New Mixed Model for Spatial Dependence. In M. Halloran and D. Berry (Eds.), Statistical Models in Epidemiology, the Environment and Clinical Trials, 135-178.

A. Riebler, S.H. Sorbye, D. Simpson and H. Rue (2016). An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical methods in medical research, 25(4), 1145-1165.

Compute MCMC diagnostic measures

Description

R_hat computes Gelman-Rubin convergence diagnostics based on the MCMC output in a model component, and n_eff computes the effective sample sizes, .i.e. estimates for the number of independent samples from the posterior distribution.

Usage

R_hat(dc)

n_eff(dc, useFFT = TRUE, lag.max, cl = NULL)

Arguments

Value

In case of R_hat the split R-hat convergence diagnostic for each component of the vector parameter, and in case of n_eff the effective number of independent samples for each component of the vector parameter.

References

A. Gelman and D. B. Rubin (1992). Inference from Iterative Simulation Using Multiple Sequences. Statistical Science 7, 457-511.

A. Gelman, J.B. Carlin, H.S. Stern, D.B. Dunson, A. Vehtari and D.B. Rubin (2013). Bayesian Data Analysis, 3rd edition. Chapman & Hall/CRC.

Examples

ex <- mcmcsae_example()
sampler <- create_sampler(ex$model, data=ex$dat)
sim <- MCMCsim(sampler, burnin=100, n.iter=300, thin=2, n.chain=4, store.all=TRUE)
n_eff(sim$beta)
n_eff(sim$v_sigma)
n_eff(sim$v_rho)
R_hat(sim$beta)
R_hat(sim$llh_)
R_hat(sim$v_sigma)

Convert a draws component object to another format

Description

Use to_mcmc to convert a draws component to class mcmc.list, allowing one to use MCMC diagnostic functions provided by package coda. Use as.array to convert to an array of dimension (draws, chains, parameters). The array format is supported by some packages for analysis or visualisation of MCMC simulation results, e.g. bayesplot. Use as.matrix to convert to a matrix, concatenating the chains. Finally, use to_draws_array to convert either a draws component or (a subset of components of) an mcdraws object to a draws_array object as defined in package posterior.

Usage

to_mcmc(x)

to_draws_array(x, components = NULL)

## S3 method for class 'dc'
as.array(x, ...)

## S3 method for class 'dc'
as.matrix(x, colnames = TRUE, ...)

Arguments

Value

The draws component(s) coerced to an mcmc.list object, a draws_array object, an array, or a matrix.

Examples

data(iris)
sampler <- create_sampler(Sepal.Length ~ reg(~ Petal.Length + Species, name="beta"), data=iris)
sim <- MCMCsim(sampler, burnin=100, n.chain=2, n.iter=200)
summary(sim)
if (require("coda", quietly=TRUE)) {
  mcbeta <- to_mcmc(sim$beta)
  geweke.diag(mcbeta)
}
if (require("posterior", quietly=TRUE)) {
  mcbeta <- to_draws_array(sim$beta)
  mcbeta
  draws <- to_draws_array(sim)
  str(draws)
}
str(as.array(sim$beta))
str(as.matrix(sim$beta))

# generate some example data
n <- 250
dat <- data.frame(x=runif(n), f=as.factor(sample(1:5, n, replace=TRUE)))
gd <- generate_data(~ reg(~ x + f, prior=pr_normal(precision=1), name="beta"), data=dat)
dat$y <- gd$y
sampler <- create_sampler(y ~ reg(~ x + f, name="beta"), data=dat)
sim <- MCMCsim(sampler, burnin=100, n.chain=2, n.iter=300)
str(sim$beta)
str(as.array(sim$beta))
bayesplot::mcmc_hist(as.array(sim$beta))
bayesplot::mcmc_dens_overlay(as.array(sim$beta))
# fake data simulation check:
bayesplot::mcmc_recover_intervals(as.array(sim$beta), gd$pars$beta)
bayesplot::mcmc_recover_hist(as.array(sim$beta), gd$pars$beta)

ex <- mcmcsae_example()
plot(ex$dat$fT, ex$dat$y)
sampler <- create_sampler(ex$model, data=ex$dat)
sim <- MCMCsim(sampler, burnin=100, n.chain=2, n.iter=200, store.all=TRUE)
str(sim$beta)
str(as.matrix(sim$beta))
# fake data simulation check:
bayesplot::mcmc_recover_intervals(as.matrix(sim$beta), ex$pars$beta)
bayesplot::mcmc_recover_intervals(as.matrix(sim$u), ex$pars$u)

Run a Markov Chain Monte Carlo simulation

Description

Given a sampler object this function runs a MCMC simulation and stores the posterior draws. A sampler object for a wide class of multilevel models can be created using create_sampler, but users can also define their own sampler functions, see below. MCMCsim allows to choose the parameters for which simulation results must be stored. It is possible to define derived quantities that will also be stored. To save memory, it is also possible to only store Monte Carlo means/standard errors for some large vector parameters, say. Another way to use less memory is to save the simulation results of large vector parameters to file. For parameters specified in plot.trace trace plots or pair plots of multiple parameters are displayed during the simulation.

Usage

MCMCsim(
  sampler,
  from.prior = FALSE,
  n.iter = 1000L,
  n.chain = 3L,
  thin = 1L,
  burnin = if (from.prior) 0L else 250L,
  start = NULL,
  store,
  store.all = FALSE,
  pred = NULL,
  store.mean,
  store.sds = FALSE,
  to.file = NULL,
  filename = "MCdraws_",
  write.single.prec = FALSE,
  verbose = TRUE,
  n.progress = n.iter%/%10L,
  trace.convergence = NULL,
  stop.on.convergence = FALSE,
  convergence.bound = 1.05,
  plot.trace = NULL,
  add.to.plot = TRUE,
  plot.type = "l",
  n.cores = 1L,
  cl = NULL,
  seed = NULL,
  export = NULL
)

Arguments

Details

A sampler object is an environment containing data and functions to use for sampling. The following elements of the sampler object are used by MCMCsim:

start

function to generate starting values.

draw

function to draw samples, typically from a full conditional posterior distribution.

rprior

function to draw from a prior distribution.

coef.names

list of vectors of parameter coefficient names, for vector parameters.

MHpars

vector of names of parameters that are sampled using a Metropolis-Hastings (MH) sampler; acceptance rates are kept for these parameters.

adapt

function of acceptance rates of MHpars to adapt MH-kernel, called every 100 iterations during the burn-in period.

Value

An object of class mcdraws containing posterior draws as well as some meta information.

Examples

# 1. create a sampler function
sampler <- new.env()
sampler$draw <- function(p) list(x=rnorm(1L), y=runif(1L))
# 2. do the simulation
sim <- MCMCsim(sampler, store=c("x", "y"))
str(sim)
summary(sim)

# example that requires start values or a start function
sampler$draw <- function(p) list(x=rnorm(1L), y=p$x * runif(1L))
sampler$start <- function(p) list(x=rnorm(1L), y=runif(1L))
sim <- MCMCsim(sampler, store=c("x", "y"))
summary(sim)
plot(sim, c("x", "y"))

# example using create_sampler; first generate some data
n <- 100
dat <- data.frame(x=runif(n), f=as.factor(sample(1:4, n, replace=TRUE)))
gd <- generate_data(~ reg(~ x + f, prior=pr_normal(precision=1), name="beta"), data=dat)
dat$y <- gd$y
sampler <- create_sampler(y ~ x + f, data=dat)
sim <- MCMCsim(sampler, burnin=100, n.iter=400, n.chain=2)
(summary(sim))
gd$pars

S4 methods for products of matrix objects

Description

Several methods for products of matrix objects. Here a matrix object can be an ordinary (dense) matrix or a (sparse) Matrix of class ddiMatrix-class, tabMatrix, dgCMatrix-class, or dsCMatrix-class. The return types are restricted to matrix or numeric in case of dense objects, and dgCMatrix, dsCMatrix, ddiMatrix or tabMatrix in case of sparse objects.

Usage

## S4 method for signature 'ddiMatrix,matrix'
x %*% y

## S4 method for signature 'dgCMatrix,matrix'
x %*% y

## S4 method for signature 'dsCMatrix,matrix'
x %*% y

## S4 method for signature 'tabMatrix,matrix'
x %*% y

## S4 method for signature 'matrix,tabMatrix'
x %*% y

## S4 method for signature 'tabMatrix,numLike'
x %*% y

## S4 method for signature 'ddiMatrix,dgCMatrix'
x %*% y

## S4 method for signature 'ddiMatrix,tabMatrix'
x %*% y

## S4 method for signature 'CsparseMatrix,tabMatrix'
x %*% y

## S4 method for signature 'tabMatrix,CsparseMatrix'
x %*% y

## S4 method for signature 'matrix,dgCMatrix'
tcrossprod(x, y)

## S4 method for signature 'matrix,tabMatrix'
tcrossprod(x, y)

## S4 method for signature 'matrix,ddiMatrix'
tcrossprod(x, y)

## S4 method for signature 'tabMatrix,missing'
crossprod(x, y)

## S4 method for signature 'tabMatrix,matrix'
crossprod(x, y)

## S4 method for signature 'tabMatrix,dgCMatrix'
crossprod(x, y)

## S4 method for signature 'tabMatrix,tabMatrix'
crossprod(x, y)

## S4 method for signature 'dgCMatrix,matrix'
crossprod(x, y)

Arguments

Value

A matrix object. In case one of the arguments is a regular (dense) matrix the result is a matrix as well.

Simulation based calibration

Description

Simulation based calibration

Usage

SBC_test(
  ...,
  pars,
  n.draw = 25L,
  n.sim = 20L * n.draw,
  burnin = 25L,
  thin = 2L,
  show.progress = TRUE,
  verbose = TRUE,
  n.cores = 1L,
  cl = NULL,
  seed = NULL,
  export = NULL
)

Arguments

Value

A matrix with ranks.

References

M. Modrak, A.H. Moon, S. Kim, P. Buerkner, N. Huurre, K. Faltejskova, A. Gelman and A. Vehtari (2023). Simulation-based calibration checking for Bayesian computation: The choice of test quantities shapes sensitivity. Bayesian Analysis, 1(1), 1-28.

Examples

## Not run: 
# this example may take a long time
n <- 10L
dat <- data.frame(x=runif(n))
ranks <- SBC_test(~ reg(~ 1 + x, prior=pr_normal(mean=c(0.25, 1), precision=1), name="beta"),
  family=f_gaussian(var.prior=pr_invchisq(df=1, scale=list(df=1, scale=1))),
  data=dat,
  pars=list(mu="beta[1]", beta_x="beta[2]", sigma="sigma_"),
  n.draw=9L, n.sim=10L*20L, thin=2L, burnin=20L
)
ranks

## End(Not run)

Functions for specifying the method and corresponding options for sampling from a possibly truncated and degenerate multivariate normal distribution

Description

These functions are intended for use in the method argument of create_TMVN_sampler.

Usage

m_direct(use.cholV = NULL)

m_Gibbs(slice = FALSE, eps = sqrt(.Machine$double.eps), diagnostic = FALSE)

m_HMC(Tsim = pi/2, max.events = .Machine$integer.max, diagnostic = FALSE)

m_HMCZigZag(
  Tsim = 1,
  scale = 1,
  prec.eq = NULL,
  diagnostic = FALSE,
  max.events = .Machine$integer.max,
  adapt = FALSE
)

m_softTMVN(
  sharpness = 100,
  useV = FALSE,
  CG = NULL,
  PG.approx = TRUE,
  PG.approx.m = -2L
)

Arguments

Value

A method object, for internal use only.

Compute autocovariance or autocorrelation function via Wiener-Khinchin theorem using Fast Fourier Transform

Description

Compute autocovariance or autocorrelation function via Wiener-Khinchin theorem using Fast Fourier Transform

Usage

ac_fft(x, demean = TRUE)

Arguments

Value

A matrix of the same size as x with autocovariances at all lags from 1 to the number of rows.

Return Metropolis-Hastings acceptance rates

Description

Return Metropolis-Hastings acceptance rates

Usage

acceptance_rates(obj, aggregate.chains = FALSE)

Arguments

Value

A list of acceptance rates.

Examples

ex <- mcmcsae_example()
# specify a model that requires MH sampling (in this case for a modelled
#   degrees of freedom parameter in the variance part of the model)
sampler <- create_sampler(ex$model, data=ex$dat,
  family = f_gaussian(var.model = ~vfac(factor="fA", prior=pr_invchisq(df="modeled")))
)
sim <- MCMCsim(sampler, burnin=100, n.iter=300, thin=2, n.chain=4, store.all=TRUE)
(summary(sim))
acceptance_rates(sim)

Utility function to construct a sparse aggregation matrix from a factor

Description

Utility function to construct a sparse aggregation matrix from a factor

Usage

aggrMatrix(fac, w = 1, mean = FALSE, facnames = FALSE)

Arguments

Value

A sparse aggregation matrix of class tabMatrix.

Examples

n <- 1000
f <- sample(1:100, n, replace=TRUE)
x <- runif(n)
M <- aggrMatrix(f)
all.equal(crossprod_mv(M, x), as.vector(tapply(x, f, sum)))

S4 method for generic 'anyNA' and signature 'tabMatrix'

Description

S4 method for generic 'anyNA' and signature 'tabMatrix'

Usage

## S4 method for signature 'tabMatrix'
anyNA(x, recursive = FALSE)

Arguments

Value

Whether the tabMatrix object contains missings or not.

Create a model component object for a BART (Bayesian Additive Regression Trees) component in the linear predictor

Description

This function is intended to be used on the right hand side of the formula argument to create_sampler or generate_data. It creates a BART term in the model's linear predictor. To use this model component one needs to have R package dbarts installed.

Usage

brt(
  formula,
  X = NULL,
  n.trees = 75L,
  name = "",
  debug = FALSE,
  keepTrees = FALSE,
  ...
)

Arguments

Value

An object with precomputed quantities and functions for sampling from prior or conditional posterior distributions for this model component. Intended for internal use by other package functions.

References

H.A. Chipman, E.I. Georgea and R.E. McCulloch (2010). BART: Bayesian additive regression trees. The Annals of Applied Statistics 4(1), 266-298.

J.H. Friedman (1991). Multivariate adaptive regression splines. The Annals of Statistics 19, 1-67.

Examples

# generate data, based on an example in Friedman (1991)
gendat <- function(n=150L, p=10L, sigma=1) {
  x <- matrix(runif(n * p), n, p)
  mu <- 10*sin(pi*x[, 1] * x[, 2]) + 20*(x[, 3] - 0.5)^2 + 10*x[, 4] + 5*x[, 5]
  y <- mu + sigma * rnorm(n)
  data.frame(x=x, mu=mu, y=y)
}

train <- gendat()
test <- gendat(n=25)

# keep trees for later prediction based on new data
sampler <- create_sampler(
  y ~ brt(~ . - y, name="bart", keepTrees=TRUE),
  family = f_gaussian(var.prior=pr_invchisq(df=3, scale=var(train$y))),
  data = train
)
# increase burnin and n.iter below to improve MCMC convergence
sim <- MCMCsim(sampler, n.chain=2, burnin=100, n.iter=200, thin=2,
  store.all=TRUE, verbose=FALSE)
(summ <- summary(sim))
plot(train$mu, summ$bart[, "Mean"]); abline(0, 1)
# NB prediction is currently slow

pred <- predict(sim, newdata=test,
  iters=sample(seq_len(n_draws(sim)), 50),
  show.progress=FALSE
)
(summpred <- summary(pred))
plot(test$mu, summpred[, "Mean"]); abline(0, 1)

Set options for Cholesky decomposition

Description

These options are only effective in case the matrix to be decomposed is sparse, i.p. of class dsCMatrix-class.

Usage

chol_control(perm = NULL, super = NA, ordering = 0L, inplace = TRUE)

Arguments

Value

A list with specified options used for Cholesky decomposition.

Combine multiple mcdraws objects into a single one by combining their chains

Description

This function can be used to combine the results of parallel simulations.

Usage

combine_chains(...)

Arguments

Value

A combined object of class mcdraws where the number of stored chains equals the sum of the numbers of chains in the input objects.

Combine multiple mcdraws objects into a single one by combining their draws

Description

This function is used to combine the results of parallel posterior predictive simulations.

Usage

combine_iters(...)

Arguments

Value

A combined object of class mcdraws where the number of stored draws equals the sum of the numbers of draws in the input objects.

Compute a list of design matrices for all terms in a model formula, or based on a sampler environment

Description

If sampler is provided instead of formula, the design matrices are based on the model used to create the sampler environment. In that case, if data is NULL, the design matrices stored in sampler are returned, otherwise the design matrices are computed for the provided data based on the sampler's model. The output is a list of dense or sparse design matrices for the model components with respect to data.

Usage

computeDesignMatrix(formula = NULL, data = NULL, labels = TRUE)

Arguments

Value

A list of design matrices.

Examples

n <- 1000
dat <- data.frame(
  x = rnorm(n),
  f = factor(sample(1:50, n, replace=TRUE))
)
str(computeDesignMatrix(~ x, dat)[[1]])
model <- ~ reg(~x, name="beta") + gen(~x, factor=~f, name="v")
X <- computeDesignMatrix(model, dat)
str(X)

Correlation factor structures in generic model components

Description

Element 'factor' of a model component created using function gen is a formula composed of several possible terms described below. It is used to derive a (typically sparse) precision matrix for a set of coefficients, and possibly a matrix representing a set of linear constraints to be imposed on the coefficient vector.

iid(f)

Independent effects corresponding to the levels of factor f.

RW1(f, circular=FALSE, w=NULL)

First-order random walk over the levels of factor f. The random walk can be made circular and different (fixed) weights can be attached to the innovations. If specified, w must be a positive numeric vector of length one less than the number of factor levels. For example, if the levels correspond to different times, it would often be reasonable to choose w proportional to the reciprocal time differences. For equidistant times there is generally no need to specify w.

RW2(f)

Second-order random walk.

AR1(f, phi, w=NULL, control=NULL)

First-order autoregressive correlation structure among the levels of f. Argument phi can be a single numerical value of the autoregressive parameter, or an appropriate prior specification if phi should be inferred. If not supplied, a uniform prior on (-1, 1] is assumed. For irregularly spaced AR(1) processes weights can be specified, in the same way as for RW1.

season(f, period)

Dummy seasonal with period period.

spatial(f, graph, snap, queen)

CAR spatial correlation. Argument graph can either be an object of (S4) class SpatialPolygonsDataFrame or an object of (S3) class sf. The latter can be obtained, e.g., by reading in a shape file using function st_read. Arguments snap and queen are passed to poly2nb, which computes a neighbours list. Alternatively, a neighbours list object of class nb can be passed directly as argument graph.

splines(f, knots, degree)

P-splines, i.e. penalised B-splines structure over the domain of a quantitative variable f. Arguments knots and degree are passed to splineDesign. If knots is a single value it is interpreted as the number of knots, otherwise as a vector of knot positions. By default 40 equally spaced knots are used, and a degree of 3.

custom(f, D=NULL, Q=NULL, R=NULL, derive.constraints=NULL)

Either a custom precision or incidence matrix associated with factor f can be passed to argument Q or D. Optionally a constraint matrix can be supplied as R, or constraints can be derived from the null space of the precision matrix by setting derive.constraints=TRUE.

Usage

iid(name)

RW1(name, circular = FALSE, w = NULL)

RW2(name)

AR1(name, phi, w = NULL, control = NULL)

season(name, period)

spatial(
  name,
  graph = NULL,
  snap = sqrt(.Machine$double.eps),
  queen = TRUE,
  poly.df = NULL,
  derive.constraints = FALSE
)

splines(name, knots, degree)

custom(name, D = NULL, Q = NULL, R = NULL, derive.constraints = NULL)

Arguments

References

B. Allevius (2018). On the precision matrix of an irregularly sampled AR(1) process. arXiv:1801.03791v2.

H. Rue and L. Held (2005). Gaussian Markov Random Fields. Chapman & Hall/CRC.

Examples

# example of CAR spatial random effects
if (requireNamespace("sf")) {
  # 1. load a shape file of counties in North Carolina
  nc <- sf::st_read(system.file("shape/nc.shp", package="sf"))
  # 2. generate some data according to a model with a few regression
  # effects, as well as spatial random effects
  gd <- generate_data(
    ~ reg(~ AREA + BIR74, prior=pr_normal(precision=1), name="beta") +
      gen(factor = ~ spatial(NAME, graph=nc), name="vs"),
    family=f_gaussian(var.prior = pr_invchisq(df=10, scale=1)),
    data = nc
  )
  # add the generated target variable and the spatial random effects to the
  # spatial dataframe object
  nc$y <- gd$y
  nc$vs_true <- gd$pars$vs
  # 3. fit a model to the generated data, and see to what extent the
  #    parameters used to generate the data, gd$pars, are reproduced
  sampler <- create_sampler(
    y ~ reg(~ AREA + BIR74, prior=pr_normal(precision=1), name="beta") +
    gen(factor = ~ spatial(NAME, graph=nc), name="vs"),
    data=nc
  )
  # increase burnin and n.iter below to improve MCMC convergence
  sim <- MCMCsim(sampler, store.all=TRUE, burnin=100, n.iter=200, n.chain=2, verbose=FALSE)
  (summ <- summary(sim))
  nc$vs <- summ$vs[, "Mean"]
  plot(nc[c("vs_true", "vs")])
  plot(gd$pars$vs, summ$vs[, "Mean"]); abline(0, 1, col="red")
}

Set up a sampler object for sampling from a possibly truncated and degenerate multivariate normal distribution

Description

This function sets up an object for multivariate normal sampling based on a specified precision matrix. Linear equality and inequality restrictions are supported. For sampling under inequality restrictions four algorithms are available. The default in that case is an exact Hamiltonian Monte Carlo algorithm (Pakman and Paninski, 2014). A related algorithm is the zig-zag Hamiltonian Monte Carlo method (Nishimura et al., 2021) in which momentum is sampled from a Laplace instead of normal distribution. Alternatively, a Gibbs sampling algorithm can be used (Rodriguez-Yam et al., 2004). The fourth option is a data augmentation method that samples from a smooth approximation to the truncated multivariate normal distribution (Souris et al., 2018).

Usage

create_TMVN_sampler(
  Q,
  mu = NULL,
  Xy = NULL,
  update.Q = FALSE,
  update.mu = update.Q,
  name = "x",
  coef.names = NULL,
  constraints = NULL,
  check.constraints = FALSE,
  method = NULL,
  reduce = NULL,
  chol.control = chol_control(),
  debug = FALSE
)

Arguments

Details

The componentwise Gibbs sampler uses univariate truncated normal samplers as described in Botev and L'Ecuyer (2016). These samplers are implemented in R package TruncatedNormal, but here translated to C++ for an additional speed-up.

Value

An environment for sampling from a possibly degenerate and truncated multivariate normal distribution.

Author(s)

Harm Jan Boonstra, with help from Grzegorz Baltissen

References

Z.I. Botev and P. L'Ecuyer (2016). Simulation from the Normal Distribution Truncated to an Interval in the Tail. in VALUETOOLS.

Y. Cong, B. Chen and M. Zhou (2017). Fast simulation of hyperplane-truncated multivariate normal distributions. Bayesian Analysis 12(4), 1017-1037.

Y. Li and S.K. Ghosh (2015). Efficient sampling methods for truncated multivariate normal and student-t distributions subject to linear inequality constraints. Journal of Statistical Theory and Practice 9(4), 712-732.

A. Nishimura, Z. Zhang and M.A. Suchard (2024). Zigzag path connects two Monte Carlo samplers: Hamiltonian counterpart to a piecewise deterministic Markov process. Journal of the American Statistical Association, 1-13.

A. Pakman and L. Paninski (2014). Exact Hamiltonian Monte Carlo for truncated multivariate gaussians. Journal of Computational and Graphical Statistics 23(2), 518-542.

G. Rodriguez-Yam, R.A. Davis and L.L. Scharf (2004). Efficient Gibbs sampling of truncated multivariate normal with application to constrained linear regression. Unpublished manuscript.

H. Rue and L. Held (2005). Gaussian Markov Random Fields. Chapman & Hall/CRC.

A. Souris, A. Bhattacharya and P. Debdeep (2019). The Soft Multivariate Truncated Normal Distribution with Applications to Bayesian Constrained Estimation. arXiv:1807.09155v2.

K.A. Valeriano, C.E. Galarza and L.A. Matos (2023). Moments and random number generation for the truncated elliptical family of distributions. Statistics and Computing 33(1), 1-20.

Examples

S <- cbind(diag(2), c(-1, 1), c(1.1, -1))  # inequality matrix
# S'x >= 0 represents the wedge x1 <= x2 <= 1.1 x1
C <- set_constraints(S=S)  # create a constraints object
# example taken from Pakman and Paninski (2014)
# 1. exact Hamiltonian Monte Carlo (Pakman and Paninski, 2014)
sampler <- create_TMVN_sampler(Q=diag(2), mu=c(4, 4), constraints=C, method="HMC")
sim <- MCMCsim(sampler, n.iter=600, verbose=FALSE)
summary(sim)
plot(as.matrix(sim$x), pch=".")
# 2. exact Hamiltonian Monte Carlo with Laplace momentum (Nishimura et al., 2021)
sampler <- create_TMVN_sampler(Q=diag(2), mu=c(4, 4), constraints=C, method="HMCZigZag")
sim <- MCMCsim(sampler, n.iter=600, verbose=FALSE)
summary(sim)
plot(as.matrix(sim$x), pch=".")
# 3. Gibbs sampling approach (Rodriguez-Yam et al., 2004)
sampler <- create_TMVN_sampler(Q=diag(2), mu=c(4, 4), constraints=C, method="Gibbs")
sim <- MCMCsim(sampler, burnin=500, n.iter=2000, verbose=FALSE)
summary(sim)
plot(as.matrix(sim$x), pch=".")
# 4. soft TMVN approximation (Souris et al., 2018)
sampler <- create_TMVN_sampler(Q=diag(2), mu=c(4, 4), constraints=C, method="softTMVN")
sim <- MCMCsim(sampler, n.iter=600, verbose=FALSE)
summary(sim)
plot(as.matrix(sim$x), pch=".")

Set up a a function for direct sampling from a constrained multivariate normal distribution

Description

Set up a a function for direct sampling from a constrained multivariate normal distribution

Usage

create_block_cMVN_sampler(
  mbs,
  X,
  Q,
  R = NULL,
  r = NULL,
  sampler,
  name = "x",
  chol.control
)

Arguments

Value

An environment with precomputed quantities and functions for sampling from a multivariate normal distribution subject to equality constraints.

Set up a function for direct sampling from a constrained multivariate normal distribution

Description

Set up a function for direct sampling from a constrained multivariate normal distribution

Usage

create_cMVN_sampler(
  D = NULL,
  Q = NULL,
  update.Q = FALSE,
  R = NULL,
  r = NULL,
  eps1 = sqrt(.Machine$double.eps),
  eps2 = sqrt(.Machine$double.eps),
  chol.control = chol_control(perm = TRUE)
)

Arguments

Value

An environment with precomputed quantities and a method 'draw' for sampling from a multivariate normal distribution subject to equality constraints.

Create a sampler object

Description

This function sets up a sampler object, based on the specification of a model. The object contains functions to draw a set of model parameters from their prior and conditional posterior distributions, and to generate starting values for the MCMC simulation. The functions share a common environment containing precomputed quantities such as design matrices based on the model and the data. The sampler object is the main input for the MCMC simulation function MCMCsim.

Usage

create_sampler(
  formula,
  data = NULL,
  family = "gaussian",
  ny = NULL,
  ry = NULL,
  r.mod = NULL,
  sigma.fixed = NULL,
  sigma.mod = NULL,
  Q0 = NULL,
  formula.V = NULL,
  logJacobian = NULL,
  linpred = NULL,
  compute.weights = FALSE,
  block = NULL,
  prior.only = FALSE,
  control = sampler_control()
)

Arguments

Details

The right hand side of the formula argument to create_sampler can be used to specify additive model components. Currently four model components are supported: reg(...) for regression or 'fixed' effects, gen(...) for generic random effects, mec(...) for measurement in covariates effects, and brt(...) for a Bayesian additive regression trees component. Note that an offset can be added separately, in the usual way using offset(...).

For gaussian models, formula.V can be used to specify the variance structure of the model. Currently two specialised variance model components are supported, vreg(...) for regression effects predicting the log-variance and vfac(...) for modelled variance factors.

Value

A sampler object, which is the main input for the MCMC simulation function MCMCsim. The sampler object is an environment with precomputed quantities and functions. The main functions are rprior, which returns a sample from the prior distributions, draw, which returns a sample from the full conditional posterior distributions, and start, which returns a list with starting values for the Gibbs sampler. If prior.only is TRUE, functions draw and start are not created.

References

J.H. Albert and S. Chib (1993). Bayesian analysis of binary and polychotomous response data. Journal of the American statistical Association 88(422), 669-679.

D. Bates, M. Maechler, B. Bolker and S.C. Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software 67(1), 1-48.

S.W. Linderman, M.J. Johnson and R.P. Adams (2015). Dependent multinomial models made easy: Stick-breaking with the Polya-Gamma augmentation. Advances in Neural Information Processing Systems, 3456-3464.

P.A. Parker, S.H. Holan and R. Janicki (2024). Conjugate Modeling Approaches for Small Area Estimation with Heteroscedastic Structure. Journal of Survey Statistics and Methodology 12(4), 1061-1080.

N. Polson, J.G. Scott and J. Windle (2013). Bayesian Inference for Logistic Models Using Polya-Gamma Latent Variables. Journal of the American Statistical Association 108(504), 1339-1349.

H. Rue and L. Held (2005). Gaussian Markov Random Fields. Chapman & Hall/CRC.

Examples

# first generate some data
n <- 200
x <- rnorm(n)
y <- 0.5 + 2*x + 0.3*rnorm(n)
# create a sampler for a simple linear regression model
sampler <- create_sampler(y ~ x)
sim <- MCMCsim(sampler)
(summary(sim))

y <- rbinom(n, 1, 1 / (1 + exp(-(0.5 + 2*x))))
# create a sampler for a binary logistic regression model
sampler <- create_sampler(y ~ x, family="binomial")
sim <- MCMCsim(sampler)
(summary(sim))

Specify a binomial sampling distribution

Description

This function can be used in the family argument of create_sampler or generate_data to specify a binomial sampling distribution. This includes the special case of binary (Bernoulli) data.

Usage

f_binomial(link = c("logit", "probit"), n.trial = NULL)

Arguments

Value

A family object.

Specify a Gamma sampling distribution

Description

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

Usage

f_gamma(
  link = "log",
  shape.vec = ~1,
  shape.prior = pr_gamma(0.1, 0.1),
  control = set_MH(type = "RWLN", scale = 0.2, adaptive = TRUE)
)

Arguments

Value

A family object.

References

J.W. Miller (2019). Fast and Accurate Approximation of the Full Conditional for Gamma Shape Parameters. Journal of Computational and Graphical Statistics 28(2), 476-480.

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
)

Arguments

Value

A family object.

Specify a Gaussian-Gamma sampling distribution

Description

This function can be used in the family argument of create_sampler or generate_data to specify a Gaussian-Gamma sampling distribution, i.e., a Gaussian sampling distribution whose variances are observed subject to error according to a Gamma distribution.

Usage

f_gaussian_gamma(link = "identity", var.model, ...)

Arguments

Value

A family object.

Specify a multinomial sampling distribution

Description

This function can be used in the family argument of create_sampler or generate_data to specify a multinomial sampling distribution. This includes the special case of categorical (multinoulli) data.

Usage

f_multinomial(link = "logit", n.trial = NULL, K = NULL)

Arguments

Value

A family object.

Specify a negative binomial sampling distribution

Description

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

Usage

f_negbinomial(
  link = "log",
  shape.vec = ~1,
  inv.shape.prior = pr_invchisq(df = 1, scale = 1),
  control = negbin_control()
)

Arguments

Details

The negative binomial distribution with shape r and probability p has density

p(y|r, p) = {\Gamma(y + r)\over y!\Gamma(r)}(1-p)^r p^y

with mean \mu = E(y|r,p) = {rp\over 1-p} and variance V(y|r,p) = \mu(1 + \mu/r). The second term of the variance can be interpreted as overdispersion with respect to a Poisson distribution, which would correspond to the limit r \rightarrow \infty. So the reciprocal shape 1/r is an overdispersion parameter, which typically is inferred. It is assigned a default prior, which may be changed through argument inv.shape.prior.

The only supported link function is "log". Strictly speaking the relation between mean \mu and linear predictor \eta is

\log\mu = \log r + \log{p\over 1-p} = \log r + \eta

This way the likelihood function has the same form as that of logistic binomial regression, so that a Polya-Gamma data augmentation sampling algorithm can be employed. Note that the fact that the linear predictor \eta does not include \log r effectively changes the interpretation of its intercept.

Value

A family object.

References

N. Polson, J.G. Scott and J. Windle (2013). Bayesian Inference for Logistic Models Using Polya-Gamma Latent Variables. Journal of the American Statistical Association 108(504), 1339-1349.

M. Zhou and L. Carin (2015). Negative Binomial Process Count and Mixture Modeling. IEEE Transactions on Pattern Analysis and Machine Intelligence 37(2), 307-320.

Specify a Poisson sampling distribution

Description

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

Usage

f_poisson(link = "log", control = poisson_control())

Arguments

Value

A family object.

Create a model component object for a generic random effects component in the linear predictor

Description

This function is intended to be used on the right hand side of the formula argument to create_sampler or generate_data.

Usage

gen(
  formula = ~1,
  factor = NULL,
  remove.redundant = FALSE,
  drop.empty.levels = FALSE,
  X = NULL,
  var = NULL,
  prior = NULL,
  Q0 = NULL,
  PX = NULL,
  priorA = NULL,
  strucA = GMRF_structure(),
  GMRFconstr = NULL,
  constraints0 = NULL,
  constraintsA = NULL,
  formula.gl = NULL,
  a = 1000,
  name = "",
  sparse = NULL,
  control = gen_control(),
  debug = FALSE
)

Arguments

Value

An object with precomputed quantities and functions for sampling from prior or conditional posterior distributions for this model component. Intended for internal use by other package functions.

References

J. Besag and C. Kooperberg (1995). On Conditional and Intrinsic Autoregression. Biometrika 82(4), 733-746.

C.M. Carvalho, N.G. Polson and J.G. Scott (2010). The horseshoe estimator for sparse signals. Biometrika 97(2), 465-480.

L. Fahrmeir, T. Kneib and S. Lang (2004). Penalized Structured Additive Regression for Space-Time Data: a Bayesian Perspective. Statistica Sinica 14, 731-761.

A. Gelman (2006). Prior distributions for variance parameters in hierarchical models. Bayesian Analysis 1(3), 515-533.

A. Gelman, D.A. Van Dyk, Z. Huang and W.J. Boscardin (2008). Using Redundant Parameterizations to Fit Hierarchical Models. Journal of Computational and Graphical Statistics 17(1), 95-122.

T. Park and G. Casella (2008). The Bayesian Lasso. Journal of the American Statistical Association 103(482), 681-686.

H. Rue and L. Held (2005). Gaussian Markov Random Fields. Chapman & Hall/CRC.

Set computational options for the sampling algorithms used for a 'gen' model component

Description

Set computational options for the sampling algorithms used for a 'gen' model component

Usage

gen_control(MHprop = c("GiG", "LNRW"))

Arguments

Value

A list of computational options regarding a 'gen' model component.

Generate a data vector according to a model

Description

This function generates draws from the prior predictive distribution. Parameter values are drawn from their priors, and consequently data is generated from the sampling distribution given these parameter values.

Usage

generate_data(
  formula,
  data = NULL,
  family = "gaussian",
  ny = NULL,
  ry = NULL,
  r.mod = NULL,
  sigma.fixed = NULL,
  sigma.mod = NULL,
  Q0 = NULL,
  formula.V = NULL,
  linpred = NULL
)

Arguments

Value

A list with a generated data vector and a list of prior means of the parameters. The parameters are drawn from their priors.

Examples

n <- 250
dat <- data.frame(
  x = rnorm(n),
  g = factor(sample(1:10, n, replace=TRUE)),
  ny = 10
)
gd <- generate_data(
  ~ reg(~ 1 + x, prior=pr_normal(precision=10, mean=c(0, 1)), name="beta") +
    gen(factor = ~ g, name="v"),
  family=f_binomial(n.trial = ~ ny), data=dat
)
gd
plot(dat$x, gd$y)

Extract a list of parameter values for a single draw

Description

Extract a list of parameter values for a single draw

Usage

get_draw(obj, iter, chain)

Arguments

Value

A list with all parameter values of draw iter from chain chain.

Examples

ex <- mcmcsae_example(n=50)
sampler <- create_sampler(ex$model, data=ex$dat)
sim <- MCMCsim(sampler, burnin=100, n.iter=300, thin=2, n.chain=4, store.all=TRUE)
get_draw(sim, iter=20, chain=3)

Create a model object for group-level regression effects within a generic random effects component.

Description

This function is intended to be used to specify the formula.gl argument to the gen model component specification function. Group-level predictors and hierarchical centring are not used by default, and they currently cannot be used in a model component that is sampled together with another model component in the same Gibbs block.

Usage

glreg(
  formula = NULL,
  remove.redundant = FALSE,
  prior = NULL,
  Q0 = NULL,
  data = NULL,
  name = ""
)

Arguments

Value

An object with precomputed quantities for sampling from prior or conditional posterior distributions for this model component. Only intended for internal use by other package functions.

Get and set the variable labels of a draws component object for a vector-valued parameter

Description

Get and set the variable labels of a draws component object for a vector-valued parameter

Usage

## S3 method for class 'dc'
labels(object, ...)

labels(object) <- value

Arguments

Value

The extractor function returns the variable labels.

Examples

ex <- mcmcsae_example()
sampler <- create_sampler(ex$model, data=ex$dat)
sim <- MCMCsim(sampler, burnin=50, n.iter=100, n.chain=1, store.all=TRUE)
labels(sim$beta)
labels(sim$v)
labels(sim$beta) <- c("a", "b")
labels(sim$beta)

Fast matrix-vector multiplications

Description

Functions for matrix-vector multiplies like %*% and crossprod, but often faster for the matrix types supported. The return value is always a numeric vector.

Usage

M %m*v% v

crossprod_mv(M, v)

Arguments

Value

For %m*v% the vector Mv and for crossprod_mv the vector M'v where M' denotes the transpose of M.

Examples

M <- matrix(rnorm(10*10), 10, 10)
x <- rnorm(10)
M %m*v% x
crossprod_mv(M, x)
M <- Matrix::rsparsematrix(100, 100, nnz=100)
x <- rnorm(100)
M %m*v% x
crossprod_mv(M, x)

Maximise the log-likelihood or log-posterior as defined by a sampler closure

Description

Maximise the log-likelihood or log-posterior as defined by a sampler closure

Usage

maximize_log_lh_p(
  sampler,
  type = c("llh", "lpost"),
  method = "BFGS",
  control = list(fnscale = -1),
  ...
)

Arguments

Value

A list of parameter values that, provided the optimisation was successful, maximize the (log-)likelihood or (log-)posterior.

Examples

n <- 1000
dat <- data.frame(
  x = rnorm(n),
  f = factor(sample(1:50, n, replace=TRUE))
)
df <- generate_data(
  ~ reg(~x, name="beta", prior=pr_normal(precision=1)) + gen(~x, factor=~f, name="v"),
  family=f_gaussian(var.prior=pr_fixed(value=1)), data=dat
)
dat$y <- df$y
sampler <- create_sampler(y ~ x + gen(~x, factor=~f, name="v"), data=dat)
opt <- maximize_log_lh_p(sampler)
str(opt)
plot(df$par$v, opt$par$v); abline(0, 1, col="red")

Create a model component object for an offset, i.e. fixed, non-parametrised term in the linear predictor

Description

This function is intended to be used on the right hand side of the formula argument to create_sampler or generate_data.

Usage

mc_offset(formula, value = NULL, name = "")

Arguments

Value

An model component object with data and methods needed for dealing with an offset term in model estimation, and prior and posterior prediction. Intended for internal use by other package functions.

Generate artificial data according to an additive spatio-temporal model

Description

This function is used to generate data for several examples.

Usage

mcmcsae_example(n = 100L, family = "gaussian")

Arguments

Value

A list containing the generated dataset, the values of the model parameters, and the model specification as a formula.

Examples

ex <- mcmcsae_example()
str(ex)

Create a model component object for a regression (fixed effects) component in the linear predictor with measurement errors in quantitative covariates

Description

This function is intended to be used on the right hand side of the formula argument to create_sampler or generate_data. It creates an additive regression term in the model's linear predictor. Covariates are assumed to be measured subject to normally distributed errors with zero mean and variance specified using the formula or V arguments. Note that this means that formula should only contain quantitative variables, and no intercept. By default, the prior for the regression coefficients is improper uniform. A proper normal prior can be set up using function pr_normal, and passed to argument prior. It should be noted that pr_normal expects a precision matrix as input for its second argument, and that the prior variance (matrix) is taken to be the inverse of this precision matrix, where in case the model's family is "gaussian" this matrix is additionally multiplied by the residual scalar variance parameter sigma_^2.

Usage

mec(
  formula = ~1,
  sparse = NULL,
  X = NULL,
  V = NULL,
  prior = NULL,
  Q0 = NULL,
  b0 = NULL,
  constraints = NULL,
  name = "",
  debug = FALSE
)

Arguments

Value

An object with precomputed quantities and functions for sampling from prior or conditional posterior distributions for this model component. Intended for internal use by other package functions.

References

L.M. Ybarra and S.L. Lohr (2008). Small area estimation when auxiliary information is measured with error. Biometrika 95(4), 919-931.

S. Arima, G.S. Datta and B. Liseo (2015). Bayesian estimators for small area models when auxiliary information is measured with error. Scandinavian Journal of Statistics 42(2), 518-529.

Examples

# example of Ybarra and Lohr (2008)
m <- 50
X <- rnorm(m, mean=5, sd=3)  # true covariate values
v <- rnorm(m, sd=2)
theta <- 1 + 3*X + v  # true values
psi <- rgamma(m, shape=4.5, scale=2)
e <- rnorm(m, sd=sqrt(psi))  # sampling error
y <- theta + e  # direct estimates
C <- c(rep(3, 10), rep(0, 40))  # measurement error for first 10 values
W <- X + rnorm(m, sd=sqrt(C))  # covariate subject to measurement error

# fit Ybarra-Lohr model
sampler <- create_sampler(
  y ~ 1 + mec(~ 0 + W, V=C, name="ME") + gen(factor=~local_),
  family=f_gaussian(var.prior=pr_fixed(1), var.vec=~psi),
  linpred="fitted"
)
sim <- MCMCsim(sampler, n.iter=800, n.chain=2, store.all=TRUE, verbose=FALSE)
(summ <- summary(sim))
plot(X, W, xlab="true X", ylab="inferred X")
points(X, summ$ME_X[, "Mean"], col="green")
abline(0, 1, col="red")
legend("topleft", legend=c("prior mean", "posterior mean"), col=c("black", "green"), pch=c(1,1))

Compute DIC, WAIC and leave-one-out cross-validation model measures

Description

Compute the Deviance Information Criterion (DIC) or Watanabe-Akaike Information Criterion (WAIC) from an object of class mcdraws output by MCMCsim. Method waic.mcdraws computes WAIC using package loo. Method loo.mcdraws also depends on package loo to compute a Pareto-smoothed importance sampling (PSIS) approximation to leave-one-out cross-validation.

Usage

compute_DIC(x, use.pV = FALSE)

compute_WAIC(
  x,
  diagnostic = FALSE,
  batch.size = NULL,
  show.progress = TRUE,
  cl = NULL,
  n.cores = 1L
)

## S3 method for class 'mcdraws'
waic(x, by.unit = FALSE, ...)

## S3 method for class 'mcdraws'
loo(x, by.unit = FALSE, r_eff = FALSE, n.cores = 1L, ...)

Arguments

Value

For compute_DIC a vector with the deviance information criterion and effective number of model parameters. For compute_WAIC a vector with the WAIC model selection criterion and WAIC effective number of model parameters. Method waic returns an object of class waic, loo, see the documentation for waic in package loo. Method loo returns an object of class psis_loo, see loo.

References

D. Spiegelhalter, N. Best, B. Carlin and A. van der Linde (2002). Bayesian Measures of Model Complexity and Fit. Journal of the Royal Statistical Society B 64 (4), 583-639.

S. Watanabe (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning 11, 3571-3594.

A. Gelman, J. Hwang and A. Vehtari (2014). Understanding predictive information criteria for Bayesian models. Statistics and Computing 24, 997-1016.

A. Vehtari, D. Simpson, A. Gelman, Y. Yao and J. Gabry (2024). Pareto smoothed importance sampling. arXiv:1507.02646v9.

A. Vehtari, A. Gelman and J. Gabry (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing 27, 1413-1432.

P.-C. Buerkner, J. Gabry and A. Vehtari (2021). Efficient leave-one-out cross-validation for Bayesian non-factorized normal and Student-t models. Computational Statistics 36, 1243-1261.

Examples

ex <- mcmcsae_example(n=100)
sampler <- create_sampler(ex$model, data=ex$dat)
sim <- MCMCsim(sampler, burnin=100, n.iter=300, n.chain=4, store.all=TRUE)
compute_DIC(sim)
compute_WAIC(sim)
if (require(loo)) {
  waic(sim)
  loo(sim, r_eff=TRUE)
}

Compute possibly sparse model matrix

Description

Compute possibly sparse model matrix

Usage

model_matrix(
  formula,
  data = NULL,
  contrasts.arg = NULL,
  drop.unused.levels = FALSE,
  sparse = NULL,
  drop0 = TRUE,
  catsep = "",
  by = NULL,
  tabM = FALSE
)

Arguments

Value

Design matrix X, either an ordinary matrix or a sparse dgCMatrix.

Get the number of chains, samples per chain or the number of variables in a simulation object

Description

Get the number of chains, samples per chain or the number of variables in a simulation object

Usage

n_chains(obj)

n_draws(obj)

n_vars(dc)

Arguments

Value

The number of chains or retained samples per chain or the number of variables.

Examples

ex <- mcmcsae_example(n=50)
sampler <- create_sampler(ex$model, data=ex$dat)
sim <- MCMCsim(sampler, burnin=100, n.iter=300, thin=2, n.chain=5, store.all=TRUE)
n_chains(sim); n_chains(sim$beta)
n_draws(sim); n_draws(sim$beta)
n_vars(sim$beta); n_vars(sim$sigma_); n_vars(sim$llh_); n_vars(sim$v)
plot(sim, "beta")
n_chains(subset(sim$beta, chains=1:2))
n_draws(subset(sim$beta, draws=sample(1:n_draws(sim), 100)))
n_vars(subset(sim$u, vars=1:2))

Set computational options for the sampling algorithms

Description

Set computational options for the sampling algorithms

Usage

negbin_control(CRT.approx.m = 20L)

Arguments

Value

A list with computational options for the sampling algorithm.

Get the parameter names from an mcdraws object

Description

Get the parameter names from an mcdraws object

Usage

par_names(obj)

Arguments

Value

The names of the parameters whose MCMC simulations are stored in obj.

Examples

data(iris)
sampler <- create_sampler(Sepal.Length ~ 
    reg(~ Petal.Length + Species, name="beta"), data=iris)
sim <- MCMCsim(sampler, burnin=100, n.iter=400)
(summary(sim))
par_names(sim)

Trace, density and autocorrelation plots for (parameters of a) draws component (dc) object

Description

Trace, density and autocorrelation plots for (parameters of a) draws component (dc) object

Usage

## S3 method for class 'dc'
plot(x, nrows, ncols, ask = FALSE, ...)

Arguments

Examples

ex <- mcmcsae_example(n=50)
sampler <- create_sampler(ex$model, data=ex$dat)
sim <- MCMCsim(sampler, store.all=TRUE)
plot(sim$u)

Trace, density and autocorrelation plots

Description

Trace, density and autocorrelation plots for selected components of an mcdraws object.

Usage

## S3 method for class 'mcdraws'
plot(x, vnames, nrows, ncols, ask = FALSE, ...)

Arguments

Examples

ex <- mcmcsae_example(n=50)
sampler <- create_sampler(ex$model, data=ex$dat)
sim <- MCMCsim(sampler, store.all=TRUE)
plot(sim, c("beta", "u", "u_sigma", "v_sigma"), ask=TRUE)

Plot a set of model coefficients or predictions with uncertainty intervals based on summaries of simulation results or other objects.

Description

This function plots estimates with error bars. Multiple sets of estimates can be compared. The error bars can either be based on standard errors or on explicitly specified lower and upper bounds. The function is adapted from function plot.sae in package hbsae, which in turn was adapted from function coefplot.default from package arm.

Usage

plot_coef(
  ...,
  n.se = 1,
  est.names,
  sort.by = NULL,
  decreasing = FALSE,
  index = NULL,
  maxrows = 50L,
  maxcols = 6L,
  offset = 0.1,
  cex.var = 0.8,
  mar = c(0.1, 2.1, 5.1, 0.1)
)

Arguments

Examples

# create artificial data
set.seed(21)
n <- 100
dat <- data.frame(
  x=runif(n),
  f=factor(sample(1:20, n, replace=TRUE))
)
model <- ~ reg(~ x, prior=pr_normal(precision=1), name="beta") + gen(factor=~f, name="v")
gd <- generate_data(model, data=dat)
dat$y <- gd$y
# fit a base model
model0 <- y ~ reg(~ 1, name="beta") + gen(factor=~f, name="v")
sampler <- create_sampler(model0, data=dat)
sim <- MCMCsim(sampler, store.all=TRUE)
(summ0 <- summary(sim))
# fit 'true' model
model <- y ~ reg(~ x, name="beta") + gen(factor=~f, name="v")
sampler <- create_sampler(model, data=dat)
sim <- MCMCsim(sampler, store.all=TRUE)
(summ <- summary(sim))
# compare random effect estimates against true parameter values
plot_coef(summ0$v, summ$v, list(est=gd$pars$v), n.se=2, offset=0.2,
  maxrows=10, est.names=c("base model", "true model", "true"))

Set computational options for the sampling algorithms

Description

Set computational options for the sampling algorithms

Usage

poisson_control(nb.shape = 100)

Arguments

Value

A list with computational options for the sampling algorithm.

Get means or standard deviations of parameters from the MCMC output in an mcdraws object

Description

Get means or standard deviations of parameters from the MCMC output in an mcdraws object

Usage

get_means(obj, vnames = NULL)

get_sds(obj, vnames = NULL)

Arguments

Value

A list with simulation means or standard deviations.

Examples

ex <- mcmcsae_example(n=50)
sampler <- create_sampler(ex$model, data=ex$dat)
sim <- MCMCsim(sampler, burnin=100, n.iter=300, thin=2, n.chain=4)
get_means(sim)
get_means(sim, "e_")
sim <- MCMCsim(sampler, burnin=100, n.iter=300, thin=2, n.chain=4,
  store.mean=c("beta", "u"), store.sds=TRUE)
summary(sim, "beta")
get_means(sim, "beta")
get_sds(sim, "beta")
get_means(sim, "u")
get_sds(sim, "u")

Create an object representing a Multivariate Log inverse Gamma (MLiG) prior distribution

Description

Create an object representing a Multivariate Log inverse Gamma (MLiG) prior distribution

Usage

pr_MLiG(mean = 0, precision = 0, labels = NULL, a = 1000)

Arguments

Value

An environment representing the specified prior, for internal use.

References

J.R. Bradley, S.H. Holan and C.K. Wikle (2018). Computationally efficient multivariate spatio-temporal models for high-dimensional count-valued data (with discussion). Bayesian Analysis 13(1), 253-310.

Create an object representing beta prior distributions

Description

Create an object representing beta prior distributions

Usage

pr_beta(a = 1, b = 1)

Arguments

Value

An environment representing the specified prior, for internal use.

Create an object representing exponential prior distributions

Description

Create an object representing exponential prior distributions

Usage

pr_exp(scale = 1)

Arguments

Value

An environment representing the specified prior, for internal use.

Create an object representing a degenerate prior fixing a parameter (vector) to a fixed value

Description

Create an object representing a degenerate prior fixing a parameter (vector) to a fixed value

Usage

pr_fixed(value = 1)

Arguments

Value

An environment representing the specified prior, for internal use.

Create an object representing gamma prior distributions

Description

Create an object representing gamma prior distributions

Usage

pr_gamma(shape = 1, rate = 1)

Arguments

Value

An environment representing the specified prior, for internal use.

Create an object representing Generalised Inverse Gaussian (GIG) prior distributions

Description

Create an object representing Generalised Inverse Gaussian (GIG) prior distributions

Usage

pr_gig(a, b, p)

Arguments

Value

An environment representing the specified prior, for internal use.

Create an object representing inverse chi-squared priors with possibly modelled degrees of freedom and scale parameters

Description

Create an object representing inverse chi-squared priors with possibly modelled degrees of freedom and scale parameters

Usage

pr_invchisq(df = 1, scale = 1)

Arguments

Value

An environment representing the specified prior, for internal use.

Create an object representing an inverse Wishart prior, possibly with modelled scale matrix

Description

Create an object representing an inverse Wishart prior, possibly with modelled scale matrix

Usage

pr_invwishart(df = NULL, scale = NULL)

Arguments

Value

An environment representing the specified prior, for internal use.

References

A. Huang and M.P. Wand (2013). Simple marginally noninformative prior distributions for covariance matrices. Bayesian Analysis 8, 439-452.

Create an object representing a possibly multivariate normal prior distribution

Description

Create an object representing a possibly multivariate normal prior distribution

Usage

pr_normal(mean = 0, precision = 0, labels = NULL)

Arguments

Value

An environment representing the specified prior, for internal use.

Create an object representing truncated normal prior distributions

Description

Create an object representing truncated normal prior distributions

Usage

pr_truncnormal(mean = 0, precision = 1, lower = 0, upper = Inf)

Arguments

Value

An environment representing the specified prior, for internal use.

Create an object representing uniform prior distributions

Description

Create an object representing uniform prior distributions

Usage

pr_unif(min = 0, max = 1)

Arguments

Value

An environment representing the specified prior, for internal use.

Generate draws from the predictive distribution

Description

Generate draws from the predictive distribution

Usage

## S3 method for class 'mcdraws'
predict(
  object,
  newdata = NULL,
  X. = if (is.null(newdata)) "in-sample" else NULL,
  type = c("data", "link", "response", "data_cat"),
  weights = NULL,
  fun. = identity,
  labels = NULL,
  ppcheck = FALSE,
  iters = NULL,
  to.file = FALSE,
  filename,
  write.single.prec = FALSE,
  show.progress = TRUE,
  verbose = TRUE,
  n.cores = 1L,
  cl = NULL,
  seed = NULL,
  export = NULL,
  ...
)

Arguments

Value

An object of class dc, containing draws from the posterior (or prior) predictive distribution. If ppcheck=TRUE posterior predictive p-values are returned as an additional attribute. In case to.file=TRUE the file name used is returned.

Examples

n <- 250
dat <- data.frame(x=runif(n))
dat$y <- 1 + dat$x + rnorm(n)
sampler <- create_sampler(y ~ x, data=dat)
sim <- MCMCsim(sampler)
summary(sim)
# in-sample prediction
pred <- predict(sim, ppcheck=TRUE)
hist(attr(pred, "ppp"))
# out-of-sample prediction
pred <- predict(sim, newdata=data.frame(x=seq(0, 1, by=0.1)))
summary(pred)

Display a summary of a dc object

Description

Display a summary of a dc object

Usage

## S3 method for class 'dc_summary'
print(
  x,
  digits = 3L,
  max.lines = 1000L,
  tail = FALSE,
  sort = NULL,
  max.label.length = NULL,
  ...
)

Arguments

Examples

ex <- mcmcsae_example()
sampler <- create_sampler(ex$model, data=ex$dat)
sim <- MCMCsim(sampler, store.all=TRUE)
print(summary(sim$u), sort="n_eff")

Print a summary of MCMC simulation results

Description

Display a summary of an mcdraws object, as output by MCMCsim.

Usage

## S3 method for class 'mcdraws_summary'
print(x, digits = 3L, max.lines = 10L, tail = FALSE, sort = NULL, ...)

Arguments

Examples

ex <- mcmcsae_example()
sampler <- create_sampler(ex$model, data=ex$dat)
sim <- MCMCsim(sampler, store.all=TRUE)
print(summary(sim), sort="n_eff")

Read MCMC draws from a file

Description

Read draws written to file by MCMCsim used with argument to.file.

Usage

read_draws(name, filename = paste0("MCdraws_", name, ".dat"))

Arguments

Value

An object of class dc containing MCMC draws for a (vector) parameter.

Examples

## Not run: 
# NB this example creates a file "MCdraws_e_.dat" in the working directory
n <- 100
dat <- data.frame(x=runif(n), f=as.factor(sample(1:5, n, replace=TRUE)))
gd <- generate_data(~ reg(~ x + f, prior=pr_normal(precision=1), name="beta"), data=dat)
dat$y <- gd$y
sampler <- create_sampler(y ~ reg(~ x + f, name="beta"), data=dat)
# run the MCMC simulation and write draws of residuals to file:
sim <- MCMCsim(sampler, n.iter=500, to.file="e_")
summary(sim)
mcres <- read_draws("e_")
summary(mcres)

## End(Not run)

Create a model component object for a regression (fixed effects) component in the linear predictor

Description

This function is intended to be used on the right hand side of the formula argument to create_sampler or generate_data. It creates an additive regression term in the model's linear predictor. By default, the prior for the regression coefficients is improper uniform. A proper normal prior can be set up using function pr_normal, and passed to argument prior. It should be noted that pr_normal expects a precision matrix as input for its second argument, and that the prior variance (matrix) is taken to be the inverse of this precision matrix, where in case the model's family is "gaussian" this matrix is additionally multiplied by the residual scalar variance parameter sigma_^2.

Usage

reg(
  formula = ~1,
  remove.redundant = FALSE,
  sparse = NULL,
  X = NULL,
  prior = NULL,
  Q0 = NULL,
  b0 = NULL,
  constraints = NULL,
  name = "",
  debug = FALSE
)

Arguments

Value

An object with precomputed quantities and functions for sampling from prior or conditional posterior distributions for this model component. Intended for internal use by other package functions.

Examples

data(iris)
# default: flat priors on regression coefficients
sampler <- create_sampler(Sepal.Length ~
    reg(~ Petal.Length + Species, name="beta"),
  data=iris
)
sim <- MCMCsim(sampler, burnin=100, n.iter=400)
summary(sim)
# (weakly) informative normal priors on regression coefficients
sampler <- create_sampler(Sepal.Length ~
    reg(~ Petal.Length + Species, prior=pr_normal(precision=1e-2), name="beta"),
  data=iris
)
sim <- MCMCsim(sampler, burnin=100, n.iter=400)
summary(sim)
# binary regression
sampler <- create_sampler(Species == "setosa" ~
    reg(~ Sepal.Length, prior=pr_normal(precision=0.1), name="beta"),
  family="binomial", data=iris)
sim <- MCMCsim(sampler, burnin=100, n.iter=400)
summary(sim)
pred <- predict(sim)
str(pred)
# example with equality constrained regression effects
n <- 500
df <- data.frame(x=runif(n))
df$y <- rnorm(n, 1 + 2*df$x)
R <- matrix(1, 2, 1)
r <- 3
C <- set_constraints(R=R, r=r)
sampler <- create_sampler(y ~ reg(~ 1 + x, constraints=C, name="beta"), data=df)
sim <- MCMCsim(sampler)
summary(sim)
plot(sim, "beta")
summary(transform_dc(sim$beta, fun=function(x) crossprod_mv(R, x) - r))

Extract draws of fitted values or residuals from an mcdraws object

Description

For a model created with create_sampler and estimated using MCMCsim, these functions return the posterior draws of fitted values or residuals. In the current implementation the fitted values correspond to the linear predictor and the residuals are computed as the data vector minus the fitted values, regardless of the model's distribution family. For large datasets the returned object can become very large. One may therefore select a subset of draws or chains or use mean.only=TRUE to return a vector of posterior means only.

Usage

## S3 method for class 'mcdraws'
fitted(
  object,
  mean.only = FALSE,
  units = NULL,
  chains = seq_len(n_chains(object)),
  draws = seq_len(n_draws(object)),
  matrix = FALSE,
  type = c("link", "response"),
  ...
)

## S3 method for class 'mcdraws'
residuals(
  object,
  mean.only = FALSE,
  units = NULL,
  chains = seq_len(n_chains(object)),
  draws = seq_len(n_draws(object)),
  matrix = FALSE,
  ...
)

Arguments

Value

Either a draws component object or a matrix with draws of fitted values or residuals. The residuals are always on the response scale, whereas fitted values can be on the scale of the linear predictor or the response depending on type. If mean.only=TRUE, a vector of posterior means.

Examples

ex <- mcmcsae_example(n=50)
sampler <- create_sampler(ex$model, data=ex$dat)
sim <- MCMCsim(sampler, burnin=100, n.iter=300, thin=2, store.all=TRUE)
fitted(sim, mean.only=TRUE)
summary(fitted(sim))
residuals(sim, mean.only=TRUE)
summary(residuals(sim))
bayesplot::mcmc_intervals(as.matrix(subset(residuals(sim), vars=1:20)))

Set computational options for the sampling algorithms

Description

Set computational options for the sampling algorithms

Usage

sampler_control(
  add.outer.R = TRUE,
  add.eps.I = FALSE,
  eps = sqrt(.Machine$double.eps),
  recompute.e = TRUE,
  cMVN.sampler = FALSE,
  CG = NULL,
  block = TRUE,
  block.V = TRUE,
  auto.order.block = TRUE,
  chol.control = chol_control(),
  max.size.cps.template = 100,
  PG.approx = TRUE,
  PG.approx.m = -2L
)

Arguments

Value

A list with specified computational options used by various sampling functions.

References

D. Bates, M. Maechler, B. Bolker and S.C. Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software 67(1), 1-48.

Y. Chen, T.A. Davis, W.W. Hager and S. Rajamanickam (2008). Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate. ACM Transactions on Mathematical Software 35(3), 1-14.

Set options for Metropolis-Hastings sampling

Description

Set options for Metropolis-Hastings sampling

Usage

set_MH(type = "RWTN", scale = 0.025, adaptive = NULL, ...)

Arguments

Value

An environment with variables and methods for Metropolis-Hastings sampling, for use by other package functions.

Set up a system of linear equality and/or inequality constraints

Description

Two-sided inequalities specified by S2, l2, u2 are currently transformed into the one-sided form S'x >= s, combined with any directly specified constraints of this form. Some basic consistency checks are carried out, notably regarding the dimensions of the inputs.

Usage

set_constraints(
  R = NULL,
  r = NULL,
  S = NULL,
  s = NULL,
  S2 = NULL,
  l2 = NULL,
  u2 = NULL,
  scale = FALSE
)

Arguments

Value

An environment with constraint matrices and vectors and a method to check whether a numeric vector satisfies all constraints. Returns NULL in case of no constraints.

Set up conjugate gradient sampler

Description

Set up conjugate gradient sampler

Usage

setup_CG_sampler(mbs, X, sampler, control = CG_control())

Arguments

Value

An environment with precomputed quantities and functions for multiplication by A and by an (inverse) preconditioning matrix.

References

A. Nishimura and M.A. Suchard (2022). Prior-preconditioned conjugate gradient method for accelerated Gibbs sampling in "large n, large p" Bayesian sparse regression. Journal of the American Statistical Association, 1-14.

Set up a cluster for parallel computing

Description

The cluster is set up for a number of workers by loading the mcmcsae package and setting up independent RNG streams.

Usage

setup_cluster(n.cores = NULL, seed = NULL, export = NULL)

Arguments

Value

An object representing the cluster.

Compute a Monte Carlo estimate of the marginal variances of a (I)GMRF

Description

Estimate marginal variances of a (I)GMRF prior defined in terms of a sparse precision matrix and possibly a set of equality constraints. The marginal variances might be used to rescale the precision matrix such that a default prior for a corresponding variance component is more appropriate.

Usage

sim_marg_var(
  D,
  Q = NULL,
  R = NULL,
  r = NULL,
  eps1 = 1e-09,
  eps2 = 1e-09,
  nSim = 100L
)

Arguments

Value

A vector of Monte Carlo estimates of the marginal variances.

References

S.H. Sorbye and H. Rue (2014). Scaling intrinsic Gaussian Markov random field priors in spatial modelling. Spatial Statistics, 8, 39-51.

Stop a cluster

Description

Stop a cluster set up by setup_cluster.

Usage

stop_cluster(cl)

Arguments

Value

NULL.

Select a subset of chains, samples and parameters from a draws component (dc) object

Description

Select a subset of chains, samples and parameters from a draws component (dc) object

Usage

## S3 method for class 'dc'
subset(x, chains = NULL, draws = NULL, vars = NULL, ...)

Arguments

Value

The selected part of the draws component as an object of class dc.

Examples

n <- 300
dat <- data.frame(x=runif(n), f=as.factor(sample(1:7, n, replace=TRUE)))
gd <- generate_data(~ reg(~ x + f, prior=pr_normal(precision=1), name="beta"), data=dat)
dat$y <- gd$y
sampler <- create_sampler(y ~ reg(~ x + f, name="beta"), data=dat)
sim <- MCMCsim(sampler)
(summary(sim$beta))
(summary(subset(sim$beta, chains=1)))
(summary(subset(sim$beta, chains=1, draws=sample(1:n_draws(sim), 100))))
(summary(subset(sim$beta, vars=1:2)))

Summarise a draws component (dc) object

Description

Summarise a draws component (dc) object

Usage

## S3 method for class 'dc'
summary(
  object,
  probs = c(0.05, 0.5, 0.95),
  na.rm = FALSE,
  time = NULL,
  abbr = FALSE,
  batch.size = 100L,
  ...
)

Arguments

Value

A matrix with summaries of class dc_summary.

Examples

ex <- mcmcsae_example()
sampler <- create_sampler(ex$model, data=ex$dat)
sim <- MCMCsim(sampler, store.all=TRUE)
summary(sim$u)

Summarise an mcdraws object

Description

Summarise an mcdraws object

Usage

## S3 method for class 'mcdraws'
summary(
  object,
  vnames = NULL,
  probs = c(0.05, 0.5, 0.95),
  na.rm = FALSE,
  efficiency = FALSE,
  abbr = FALSE,
  batch.size = 100L,
  ...
)

Arguments

Value

A list of class mcdraws_summary summarizing object.

Examples

ex <- mcmcsae_example()
sampler <- create_sampler(ex$model, data=ex$dat)
sim <- MCMCsim(sampler, store.all=TRUE)
summary(sim)
par_names(sim)
summary(sim, c("beta", "v_sigma", "u_sigma"))

S4 method for row and column subsetting a 'tabMatrix'

Description

S4 method for row and column subsetting a 'tabMatrix'

Usage

## S4 method for signature 'tabMatrix,index,missing,logical'
x[i, j, ..., drop = TRUE]

## S4 method for signature 'tabMatrix,index,missing,missing'
x[i, j, ..., drop = TRUE]

## S4 method for signature 'tabMatrix,missing,index,logical'
x[i, j, ..., drop = TRUE]

## S4 method for signature 'tabMatrix,missing,index,missing'
x[i, j, ..., drop = TRUE]

Arguments

Value

The selected rows/columns as a tabMatrix or a vector.

Transform one or more draws component objects into a new one by applying a function

Description

Transform one or more draws component objects into a new one by applying a function

Usage

transform_dc(..., fun, to.matrix = FALSE, labels = NULL)

Arguments

Value

Either a matrix or a draws component object.

Examples

ex <- mcmcsae_example(n=50)
sampler <- create_sampler(ex$model, data=ex$dat)
sim <- MCMCsim(sampler, burnin=100, n.iter=300, thin=2, n.chain=4, store.all=TRUE)
summary(sim$v_sigma)
summary(transform_dc(sim$v_sigma, fun=function(x) x^2))
summary(transform_dc(sim$u, sim$u_sigma, fun=function(x1, x2) abs(x1)/x2))

Create a model component object for a variance factor component in the variance function of a gaussian sampling distribution

Description

This function is intended to be used on the right hand side of the formula.V argument to create_sampler or generate_data.

Usage

vfac(
  factor = "local_",
  prior = pr_invchisq(df = 1, scale = 1),
  name = "",
  debug = FALSE
)

Arguments

Value

An object with precomputed quantities and functions for sampling from prior or conditional posterior distributions for this model component. Intended for internal use by other package functions.

Create a model component object for a regression component in the variance function of a gaussian sampling distribution

Description

This function is intended to be used on the right hand side of the formula.V argument to create_sampler or generate_data.

Usage

vreg(
  formula = NULL,
  remove.redundant = FALSE,
  sparse = NULL,
  X = NULL,
  prior = NULL,
  Q0 = NULL,
  b0 = NULL,
  name = ""
)

Arguments

Value

An object with precomputed quantities and functions for sampling from prior or conditional posterior distributions for this model component. Intended for internal use by other package functions.

References

E. Cepeda and D. Gamerman (2000). Bayesian modeling of variance heterogeneity in normal regression models. Brazilian Journal of Probability and Statistics, 207-221.

T.I. Lin and W.L. Wang (2011). Bayesian inference in joint modelling of location and scale parameters of the t distribution for longitudinal data. Journal of Statistical Planning and Inference 141(4), 1543-1553.

Extract weights from an mcdraws object

Description

Extract weights from an mcdraws object

Usage

## S3 method for class 'mcdraws'
weights(object, ...)

Arguments

Value

A vector with (simulation means of) weights.

Examples

# first create a population data frame
N <- 1000  # population size
pop <- data.frame(x=rnorm(N), area=factor(sample(1:10, N, replace=TRUE)))
pop$y <- 1 + 2*pop$x + seq(-1, to=1, length.out=10)[pop$area] + 0.5*rnorm(N)
pop$sample <- FALSE
pop$sample[sample(seq_len(N), 100)] <- TRUE
# a simple linear regression model:
sampler <- create_sampler(
  y ~ reg(~ x, name="beta"),
  linpred=list(beta=rowsum(model.matrix(~ x, pop), pop$area)), compute.weights=TRUE,
  data=pop[pop$sample, ]
)
sim <- MCMCsim(sampler)
(summary(sim))
str(weights(sim))
crossprod_mv(weights(sim), pop$y[pop$sample])
summary(sim$linpred_)
# a multilevel model:
sampler <- create_sampler(
  y ~ reg(~ x, name="beta") + gen(factor = ~ area, name="v"),
  linpred=list(beta=rowsum(model.matrix(~ x, pop), pop$area), v=diag(10)), compute.weights=TRUE,
  data=pop[pop$sample, ]
)
sim <- MCMCsim(sampler)
(summary(sim))
str(weights(sim))
crossprod_mv(weights(sim), pop$y[pop$sample])
summary(sim$linpred_)