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, distribution = "normal", data = list(), centre = TRUE, 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)", c.psiPrior = "HN(1)", phiPrior = "HN(2)",
pi.omegaPrior = "Beta(1, 1)", mu.RPrior = "N(0, 1)", sigma.RPrior = "HN(1)",
corr.Model = c("common", nClust = 1), DP.concPrior = "Gamma(5, 2)",
breaksPrior = "SBeta(1, 2)", tuneCbeta, tuneCalpha, tuneAlpha, tuneSigma2,
tuneCpsi, tunePsi, tunePhi, tuneR, tuneSigma2R, tuneHar, tuneBreaks, tunePeriod, tau,
FT = 1, compDeviance = FALSE, ...)Function mvrm returns the following:
the matched call.
model formula.
the seed that was used (in case replication of the results is needed).
the dataset
the mean model design matrix.
the variance model design matrix.
the length of the vector of indicators \(\gamma\).
the length of the vector of indicators \(\delta\).
the MCMC parameters: length of burn in period, total number of samples, thinning period.
total number of posterior samples
the storage directory
Further, function mvrm creates files where the posterior samples are written.
These files are (with all file names preceded by `BNSP.'):
contains samples from the posterior of vector \(\alpha\). Rows represent posterior samples and columns represent the regression coefficient, and they are in the same order as the columns of design matrix Z.
contains samples from the posterior of vector \(\beta\). Rows represent posterior samples and columns represent the regression coefficients, and they are in the same order as the columns of design matrix X.
contains samples from the posterior of the vector of the indicators \(\gamma\). Rows represent posterior samples and columns represent the indicator variables, and they are in the same order as the columns of design matrix X.
contains samples from the posterior of the vector of the indicators \(\delta\). Rows represent posterior samples and columns represent the indicator variables, and they are in the same order as the columns of design matrix Z.
contains samples from the posterior of the error variance \(\sigma^2\) of each response.
contains samples from the posterior of \(c_{\beta}\).
contains samples from the posterior of \(c_{\alpha}\).
contains samples from the posterior of the correlation matrix \(R\).
contains samples from the posterior of \(\theta\) of the shadow prior (probably not needed).
contains samples from the posterior of \(\mu_R\).
contains samples from the posterior of \(\sigma^2_{R}\).
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
The cluster at which the correlations were allocated, \(\lambda_{kl}\). These are integers from zero to the specified number of clusters minus one.
The numbers of correlations assigned to each cluster.
Contains samples from the posterior of the Dirichlet process concentration parameter.
In addition to the above, for models that cluster the variables we have
The cluster at which the variables were allocated, \(\lambda_{k}\). These are integers from zero to the specified number of clusters minus one.
The numbers of variables assigned to each cluster.
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.
The distribution for the response variables. Currently two options are supported: "normal" and "t".
a data frame.
Binary indicator. If set equal to TRUE, the design matrices are centred, to have column mean equl to zero, otherwise, if set to FALSE, the columns are not centred.
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).
length of burn-in period.
thinning parameter.
optional seed for the random generator.
a required directory to store files with the posterior samples of models parameters.
The inverse Gamma prior of \(c_{\beta}\). 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.
The Beta prior of \(\pi_{\mu}\). 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\)
The prior of \(c_{\alpha}\). The default is "IG(1.1,1.1)". Half-normal priors for \(\sqrt{c_{\alpha}}\) 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).
The prior of \(\sigma\). The default is "HN(2)", a half-normal prior for \(\sigma\) with variance equal to two, \(\sigma \sim N(0,2) I[\sigma>0]\). Inverse Gamma priors for \(\sigma^2\) are also available, declared using "IG(a,b)". It can be of dimension \(1\) or \(p\) (the length of the multivariate response).
The Beta prior of \(\pi_{\sigma}\). 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\)
The prior of \(c_{\psi}\). The default is half-normal for \(\sqrt{c_{\psi}}\), declared using "HN(a)", where "a" is a positive number. The default value for "a" is one. Inverse-gamma priors are also available and they can be declared using "IG(a,b)", where "a" and "b" are positive constants. The prior can be of dimension \(1\) or \(p\) (the length of the multivariate response).
The prior of \(varphi^2\). The default is half-normal for \(varphi\), declared using "HN(a)", where "a" is a positive number. The default value for "a" is two. Inverse-gamma priors are also available and they can be declared using "IG(a,b)", where "a" and "b" are positive constants. The prior can be of dimension \(1\) or \(p\) (the length of the multivariate response).
The Beta prior of \(\pi_{\omega}\). The default is "Beta(1,1)". It can be of dimension \(1\), of dimension \(B\) (the number of effects that enter the shape parameter model), or of dimension \(pB\)
The normal prior for \(\mu_{R}\). The default is the standard normal distribution.
The half normal prior for \(\sigma_{R}\). The default is the half normal distribution with variance one.
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.
The Gamma prior for the Dirichlet process concentration parameter.
The prior for the shifts associated with the growth break points. The shift is taken to have a scaled Beta prior with support the (0, p) interval, where p is the period of the sin curve. The default SBeta(1, 2) is a scaled Beta(1, 2) distribution supported in the (0, p) interval. The shifts are in increasing order.
Starting value of the tuning parameter for sampling \(c_{\beta}\). Defaults at 20.
Starting value of the tuning parameter for sampling \(c_{\alpha}\). Defaults at 1.
Starting value of the tuning parameter for sampling regression coefficients of the variance model \(\alpha\). Defaults at 5.
Starting value of the tuning parameter for sampling variances \(\sigma^2_j\). Defaults at 1.
Starting value of the tuning parameter for sampling \(c_{\psi}\). Defaults at 1.
Starting value of the tuning parameter for sampling \(\psi\). Defaults at 5.
Starting value of the tuning parameter for sampling \(\varphi\). Defaults at 0.5.
Starting value of the tuning parameter for sampling correlation matrices. Defaults at \(40*(p+2)^3\).
Starting value of the tuning parameter for sampling \(\sigma_{R}^2\). Defaults at 1.
Starting value of the tuning parameter for sampling the regression coefficients of the harmonics. Defaults at 100.
Starting value of the tuning parameter for sampling the shift parameters associated with growth breaks. Defaults at 0.01 times the period of the sin wave.
Starting value of the tuning parameter for sampling the period parameter of the sin curve. Defaults to 0.01.
The \(tau\) of the shadow prior. Defaults at 0.01.
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.
Binary indicator. If set equal to 1, the deviance is computed.
Other options that will be ignored.
Georgios Papageorgiou gpapageo@gmail.com
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
$$
\mu_u = \beta_0 + \beta_1 u_1 + f_{\mu}(u_2),$$
$$
\log(\sigma^2_W) = \alpha_0 + \alpha_1 w_1 + f_{\sigma}(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_{\mu}\)
and \(f_{\sigma}\), are represented by basis function expansion,
$$
f_{\mu}(u_2) = \sum_{j} \beta_{j} \phi_{j}(u_2),
$$
$$
f_{\sigma}(w_2) = \sum_{j} \alpha_{j} \phi_{j}(w_2),
$$
where \(\phi\) are the basis functions and \(\beta\) and \(\alpha\) are regression coefficients.
The variance model can equivalently be expressed as $$ \sigma^2_W = \exp(\alpha_0) \exp(\alpha_1 w_1 + f_{\sigma}(w_2)) = \sigma^2 \exp(\alpha_1 w_1 + f_{\sigma}(w_2)),$$ where \(\sigma^2 = \exp(\alpha_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 \(\gamma\) that take value \(\gamma=1\) if the associated coefficient \(\beta \neq 0\) and \(\gamma = 0\) if \(\beta = 0\). Indicators \(\delta\) that take value \(\delta=1\) if the associated coefficient \(\alpha \neq 0\) and \(\delta = 0\) if \(\alpha = 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, \(\beta_0\) and \(\alpha_0\).
Prior specification:
For the vector of non-zero regression coefficients \(\beta_{\gamma}\) we specify a g-prior $$ \beta_{\gamma} | c_{\beta}, \sigma^2, \gamma, \alpha, \delta \sim N(0,c_{\beta} \sigma^2 (\tilde{X}_{\gamma}^{\top} \tilde{X}_{\gamma} )^{-1}). $$ where \(\tilde{X}\) is a scaled version of design matrix \(X\) of the mean model.
For the vector of non-zero regression coefficients \(\alpha_{\delta}\) we specify a normal prior $$ \alpha_{\delta} | c_{\alpha}, \delta \sim N(0,c_{\alpha} I). $$
Independent priors are specified for the indicators variables \(\gamma\) and \(\delta\) as \(P(\gamma = 1 | \pi_{\mu}) = \pi_{\mu}\) and \(P(\delta = 1 | \pi_{\sigma}) = \pi_{\sigma}\). Further, Beta priors are specified for \(\pi_{\mu}\) and \(\pi_{\sigma}\) $$ \pi_{\mu} \sim Beta(c_{\mu},d_{\mu}), \pi_{\sigma} \sim Beta(c_{\sigma},d_{\sigma}). $$ We note that blocks of regression coefficients associated with distinct covariate effects have their own probability of selection (\(\pi_{\mu}\) or \(\pi_{\sigma}\)) and this probability has its own prior distribution.
Further, we specify inverse Gamma priors for \(c_{\beta}\) and \(c_{\alpha}\) $$ c_{\beta} \sim IG(a_{\beta},b_{\beta}), c_{\alpha} \sim IG(a_{\alpha},b_{\alpha})$$
For \(\sigma^2\) we consider inverse Gamma and half-normal priors $$ \sigma^2 \sim IG(a_{\sigma},b_{\sigma}), |\sigma| \sim N(0,\phi^2_{\sigma}). $$
Lastly, for the elements of the correlation matrix, we specify normal distributions with mean \(\mu_R\) and variance \(\sigma^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 \(\mu_R\). The grouped correlations model can also be specified. It clusters the variables instead of the correlations.
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()
if (FALSE) 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)
Run the code above in your browser using DataLab