mvrm {BNSP} | R Documentation |
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.
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, ...)
formula |
a formula defining the responses and the covariates in the mean and variance models e.g. |
data |
a data frame. |
sweeps |
total number of posterior samples, including those discarded in burn-in period
(see argument |
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. |
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.
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. |
Georgios Papageorgiou gpapageo@gmail.com
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.
# 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)