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]