mvrm {BNSP}R Documentation

Bayesian semiparametric analysis of multivariate continuous responses, with variable selection

Description

Implements an MCMC algorithm for posterior sampling based on a semiparametric model for continuous multivariate responses and additive models for the mean and variance functions. The model utilizes spike-slab priors for variable selection and regularization. See ‘Details’ section for a full description of the model.

Usage

mvrm(formula, data = list(), sweeps, burn = 0, thin = 1, seed, StorageDir,
c.betaPrior = "IG(0.5, 0.5 * n * p)", pi.muPrior = "Beta(1, 1)", 
c.alphaPrior = "IG(1.1, 1.1)", sigmaPrior = "HN(2)", pi.sigmaPrior = "Beta(1, 1)",
mu.RPrior = "N(0, 1)", sigma.RPrior = "HN(1)", corr.Model = c("common", nClust = 1),
DP.concPrior = "Gamma(5, 2)", tuneAlpha, tuneSigma2, tuneCb, tuneCa, tuneR, 
tuneSigma2R, tau, FT = 1, ...) 

Arguments

formula

a formula defining the responses and the covariates in the mean and variance models e.g. y1 | y2 ~ x | z or for smooth effects y1 | y2 ~ sm(x) | sm(z). The package uses the extended formula notation, where the responses are defined on the left of ~ and the mean and variance models on the right.

data

a data frame.

sweeps

total number of posterior samples, including those discarded in burn-in period (see argument burn) and those discarded by the thinning process (see argument thin).

burn

length of burn-in period.

thin

thinning parameter.

seed

optional seed for the random generator.

StorageDir

a required directory to store files with the posterior samples of models parameters.

c.betaPrior

The inverse Gamma prior of c_{β}. The default is "IG(0.5,0.5*n*p)", that is, an inverse Gamma with parameters 1/2 and np/2, where n is the number of sampling units and p is the length of the response vector.

pi.muPrior

The Beta prior of π_{μ}. The default is "Beta(1,1)". It can be of dimension 1, of dimension K (the number of effects that enter the mean model), or of dimension pK

c.alphaPrior

The inverse Gamma prior of c_{α}. The default is "IG(1.1,1.1)". Half-normal priors for √{c_{α}} are also available, declared using "HN(a)", where "a" is a positive number. It can be of dimension 1 or p (the length of the multivariate response).

sigmaPrior

The prior of σ. The default is "HN(2)", a half-normal prior for σ with variance equal to two, σ \sim N(0,2) I[σ>0]. Inverse Gamma priors for σ^2 are also available, declared using "IG(a,b)". It can be of dimension 1 or p (the length of the multivariate response).

pi.sigmaPrior

The Beta prior of π_{σ}. The default is "Beta(1,1)". It can be of dimension 1, of dimension Q (the number of effects that enter the variance model), or of dimension pQ

mu.RPrior

The normal prior for μ_{R}. The default is the standard normal distribution.

sigma.RPrior

The half normal prior for σ_{R}. The default is the half normal distribution with variance one.

corr.Model

Specifies the model for the correlation matrix R. The three choices supported are "common", that specifies a common correlations model, "groupC", that specifies a grouped correlations model, and "groupV", that specifies a grouped variables model. When the model chosen is either "groupC" or "groupV", the upper limit on the number of clusters can also be specified, using corr.Model = c("groupC", nClust = d) or corr.Model = c("groupV", nClust = p). If the number of clusters is left unspecified, for the "groupV" model, it is taken to be p, the number of responses. For the "groupC" model, it is taken to be d = p(p-1)/2, the number of free elements in the correlation matrix.

DP.concPrior

The Gamma prior for the Dirichlet process concentration parameter.

tuneAlpha

Starting value of the tuning parameter for sampling regression coefficients of the variance model α. Defaults at 5.

tuneSigma2

Starting value of the tuning parameter for sampling variances σ^2_j. Defaults at 1.

tuneCb

Starting value of the tuning parameter for sampling c_{β}. Defaults at 10.

tuneCa

Starting value of the tuning parameter for sampling c_{α}. Defaults at 1.

tuneR

Starting value of the tuning parameter for sampling correlation matrices. Defaults at 100(p+2).

tuneSigma2R

Starting value of the tuning parameter for sampling σ_{R}^2. Defaults at 1.

tau

The tau of the shadow prior. Defaults at 0.01.

FT

Binary indicator. If set equal to 1, the Fisher's z transform of the correlations is modelled, otherwise if set equal to 0, the untransformed correlations are modelled.

...

Other options that will be ignored.

Details

Function mvrm returns samples from the posterior distributions of the parameters of a regression model with normally distributed multivariate responses and mean and variance functions modeled in terms of covariates. For instance, in the presence of two responses (y_1, y_2) and two covariates in the mean model (u_1, u_2) and two in the variance model (w_1, w_2), we may choose to fit

μ_u = β_0 + β_1 u_1 + f_{μ}(u_2),

\log(σ^2_W) = α_0 + α_1 w_1 + f_{σ}(w_2),

parametrically modelling the effects of u_1 and w_1 and non-parametrically modelling the effects of u_2 and w_2. Smooth functions, such as f_{μ} and f_{σ}, are represented by basis function expansion,

f_{μ}(u_2) = ∑_{j} β_{j} φ_{j}(u_2),

f_{σ}(w_2) = ∑_{j} α_{j} φ_{j}(w_2),

where φ are the basis functions and β and α are regression coefficients.

The variance model can equivalently be expressed as

σ^2_W = \exp(α_0) \exp(α_1 w_1 + f_{σ}(w_2)) = σ^2 \exp(α_1 w_1 + f_{σ}(w_2)),

where σ^2 = \exp(α_0). This is the parameterization that we adopt in this implementation.

Positive prior probability that the regression coefficients in the mean model are exactly zero is achieved by defining binary variables γ that take value γ=1 if the associated coefficient β \neq 0 and γ = 0 if β = 0. Indicators δ that take value δ=1 if the associated coefficient α \neq 0 and δ = 0 if α = 0 for the variance function are defined analogously. We note that all coefficients in the mean and variance functions are subject to selection except the intercepts, β_0 and α_0.

Prior specification:

For the vector of non-zero regression coefficients β_{γ} we specify a g-prior

β_{γ} | c_{β}, σ^2, γ, α, δ \sim N(0,c_{β} σ^2 (\tilde{X}_{γ}^{\top} \tilde{X}_{γ} )^{-1}).

where \tilde{X} is a scaled version of design matrix X of the mean model.

For the vector of non-zero regression coefficients α_{δ} we specify a normal prior

α_{δ} | c_{α}, δ \sim N(0,c_{α} I).

Independent priors are specified for the indicators variables γ and δ as P(γ = 1 | π_{μ}) = π_{μ} and P(δ = 1 | π_{σ}) = π_{σ}. Further, Beta priors are specified for π_{μ} and π_{σ}

π_{μ} \sim Beta(c_{μ},d_{μ}), π_{σ} \sim Beta(c_{σ},d_{σ}).

We note that blocks of regression coefficients associated with distinct covariate effects have their own probability of selection (π_{μ} or π_{σ}) and this probability has its own prior distribution.

Further, we specify inverse Gamma priors for c_{β} and c_{α}

c_{β} \sim IG(a_{β},b_{β}), c_{α} \sim IG(a_{α},b_{α})

For σ^2 we consider inverse Gamma and half-normal priors

σ^2 \sim IG(a_{σ},b_{σ}), |σ| \sim N(0,φ^2_{σ}).

Lastly, for the elements of the correlation matrix, we specify normal distributions with mean μ_R and variance σ^2_R, with the priors on these two parameters being normal and half-normal, respectively. This is the common correlations model. Further, the grouped correlations model can be specified. It considers a mixture of normal distributions for the means μ_R. The grouped correlations model can also be specified. It clusters the variables instead of the correlations.

Value

Function mvrm returns the following:

call

the matched call.

formula

model formula.

seed

the seed that was used (in case replication of the results is needed).

data

the dataset

X

the mean model design matrix.

Z

the variance model design matrix.

LG

the length of the vector of indicators γ.

LD

the length of the vector of indicators δ.

mcpar

the MCMC parameters: length of burn in period, total number of samples, thinning period.

nSamples

total number of posterior samples

DIR

the storage directory

Further, function mvrm creates files where the posterior samples are written. These files are (with all file names preceded by ‘BNSP.’):

alpha.txt

contains samples from the posterior of vector α. Rows represent posterior samples and columns represent the regression coefficient, and they are in the same order as the columns of design matrix Z.

beta.txt

contains samples from the posterior of vector β. Rows represent posterior samples and columns represent the regression coefficients, and they are in the same order as the columns of design matrix X.

gamma.txt

contains samples from the posterior of the vector of the indicators γ. Rows represent posterior samples and columns represent the indicator variables, and they are in the same order as the columns of design matrix X.

delta.txt

contains samples from the posterior of the vector of the indicators δ. Rows represent posterior samples and columns represent the indicator variables, and they are in the same order as the columns of design matrix Z.

sigma2.txt

contains samples from the posterior of the error variance σ^2 of each response.

cbeta.txt

contains samples from the posterior of c_{β}.

calpha.txt

contains samples from the posterior of c_{α}.

R.txt

contains samples from the posterior of the correlation matrix R.

theta.txt

contains samples from the posterior of θ of the shadow prior (probably not needed).

muR.txt

contains samples from the posterior of μ_R.

sigma2R.txt

contains samples from the posterior of σ^2_{R}.

deviance.txt

contains the deviance, minus twice the log likelihood evaluated at the sampled values of the parameters.

In addition to the above, for models that cluster the correlations we have

compAlloc.txt

The cluster at which the correlations were allocated, λ_{kl}. These are integers from zero to the specified number of clusters minus one.

nmembers.txt

The numbers of correlations assigned to each cluster.

DPconc.txt

Contains samples from the posterior of the Dirichlet process concentration parameter.

In addition to the above, for models that cluster the variables we have

compAllocV.txt

The cluster at which the variables were allocated, λ_{k}. These are integers from zero to the specified number of clusters minus one.

nmembersV.txt

The numbers of variables assigned to each cluster.

Author(s)

Georgios Papageorgiou gpapageo@gmail.com

References

Papageorgiou, G. and Marshall, B.C. (2019). Bayesian semiparametric analysis of multivariate continuous responses, with variable selection. arXiv.

Papageorgiou, G. (2018). BNSP: an R Package for fitting Bayesian semiparametric regression models and variable selection. The R Journal, 10(2):526-548.

Chan, D., Kohn, R., Nott, D., & Kirby, C. (2006). Locally adaptive semiparametric estimation of the mean and variance functions in regression models. Journal of Computational and Graphical Statistics, 15(4), 915-936.

Examples

# Fit a mean/variance regression model on the cps71 dataset from package np. 
#This is a univariate response model
require(np)
require(ggplot2)
data(cps71)
model <- logwage ~ sm(age,k=30,bs="rd") | sm(age,k=30,bs="rd")
DIR<-getwd()
## Not run: m1 <- mvrm(formula=model,data=cps71,sweeps=10000,burn=5000,thin=2, seed=1,StorageDir=DIR)
#Print information and summarize the model
print(m1)
summary(m1)
#Summarize and plot one parameter of interest
alpha<-mvrm2mcmc(m1,"alpha")
summary(alpha)
plot(alpha)
#Obtain a plot of a term in the mean model
wagePlotOptions<-list(geom_point(data=cps71,aes(x=age,y=logwage)))
plot(x=m1,model="mean",term="sm(age)",plotOptions=wagePlotOptions)
plot(m1)
#Obtain predictions for new values of the predictor "age"
predict(m1,data.frame(age=c(21,65)),interval="credible")

# Fit a bivariate mean/variance model on the marks dataset from package ggm
# two responses: marks mechanics and vectors, and one covariate: marks on algebra 
model2 <- mechanics | vectors ~ sm(algebra,k=5) | sm(algebra,k=3)
m2 <- mvrm(formula=model2, data=marks, sweeps = 100000, burn = 50000, 
                       thin = 2, seed = 1, StorageDir = DIR)
plot(m2)

## End(Not run)

[Package BNSP version 2.1.6 Index]