Package 'BNSP'

Bayesian non- and semi-parametric model fitting

Description

Markov chain Monte Carlo algorithms for non- and semi-parametric models: 1. spike-slab variable selection in multivariate mean/variance regression models with function mvrm, 2. joint mean-covariance models for multivariate longitudinal responses with function lmrm, and 3. Dirichlet process mixture models with function dpmj.

Details

Package:	BNSP
Type:	Package
Version:	2.2.3
Date:	2023-05-25
License:	GPL (>=2)

This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

For details on the GNU General Public License see http://www.gnu.org/copyleft/gpl.html or write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

Author(s)

Georgios Papageorgiou (2014)

Maintainer: Georgios Papageorgiou <[email protected]>

References

Papageorgiou, G. (2020). Bayesian semiparametric modelling of covariance matrices for multivariate longitudinal data. https://arxiv.org/abs/2012.09833

Papageorgiou, G. and Marshall, B.C. (2020). Bayesian semiparametric analysis of multivariate continuous responses, with variable selection. Journal of Computational and Graphical Statistics, 29:4, 896-909, DOI: 10.1080/10618600.2020.1739534

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

Papageorgiou, G. (2018). Bayesian density regression for discrete outcomes. Australian and New Zealand Journal of Statistics, arXiv:1603.09706v3 [stat.ME].

Papageorgiou, G., Richardson, S. and Best, N. (2015). Bayesian nonparametric models for spatially indexed data of mixed type. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 77:973-999.

---

Amitriptyline dataset from Johnson and Wichern

Description

Amitriptyline is a prescription antidepressant. The dataset consists of measurements on 17 patients who had over-dosed on amitriptyline.

Usage

data(ami)

Format

A data frame containing 17 rows and 7 columns. The columns represent

tot

total blood plasma level.

ami

amount of amitriptyline found in the plasma.

gen

gender (1 for female).

amt

amount of the drug taken.

pr

PR wave measurement.

bp

diastolic blood pressure.

qrs

QRS wave measurement.

Source

Johnson, R. A., and Wichern, D. W. (2007), Applied Multivariate Statistical Analysis, Essex: Pearson, page 426.

References

Johnson, R. A., and Wichern, D. W. (2007). Applied Multivariate Statistical Analysis, Essex: Pearson.

---

The Cholesky and modified Cholesky decompositions

Description

Computes the Cholesky factorization and modified Cholesky factorizations of a real symmetric positive-definite square matrix.

Usage

chol(x, mod = TRUE, p = 1, ...)

Arguments

x	A symmetric, positive-definite matrix.
mod	Defaults to TRUE. With this choice, the function returns the modified Cholesky decomposition. When mod = FALSE, the function returns the usual Cholesky decomposition.
p	Relevant only when mod = TRUE. It determines the size of the blocks of the block diagonal matrix.
...	other arguments.

Details

The function computes the modified Cholesky decomposition of a real symmetric positive-definite square matrix $\Sigma$. This is given by

$L \Sigma L^{\top} = D,$

where $L$ is a lower tringular matrix with ones on its main diagonal and D is a block diagonal matrix with block size determined by argument p.

Value

The function returns matrices $L$ and $D$.

Author(s)

Georgios Papageorgiou [email protected]

See Also

The default function from base, chol

Examples

Sigma <- matrix(c(1.21,0.18,0.13,0.41,0.06,0.23,
                  0.18,0.64,0.10,-0.16,0.23,0.07,
                  0.13,0.10,0.36,-0.10,0.03,0.18,
                  0.41,-0.16,-0.10,1.05,-0.29,-0.08,
                  0.06,0.23,0.03,-0.29,1.71,-0.10,
                  0.23,0.07,0.18,-0.08,-0.10,0.36),6,6)
LD <- chol(Sigma)
L <- LD$L
D <- LD$D
round(L,5)
round(D,5)
solve(L) %*% D %*% solve(t(L))
LD <- chol(Sigma, p = 2)
L <- LD$L
D <- LD$D
round(L, 5)
round(D, 5)
solve(L) %*% D %*% solve(t(L))

---

Computes the similarity matrix

Description

Computes the similarity matrix.

Usage

clustering(object, ...)

Arguments

object	an object of class "mvrm", usually a result of a call to mvrm.
...	other arguments.

Details

The function computes the similarity matrix for clustering based on corrrelations or variables.

Value

Similarity matrix.

Author(s)

Georgios Papageorgiou [email protected]

See Also

mvrm

Examples

#see \code{mvrm} example

---

Continues the sampler from where it stopped

Description

Allows the user to continue the sampler from the state it stopped in the previous call to mvrm.

Usage

continue(object, sweeps, burn = 0, thin, discard = FALSE,...)

Arguments

object	An object of class "mvrm", usually a result of a call to mvrm.
sweeps	The number of additional sweeps, maintaining the same thinning interval as specified in the original call to mvrm.
burn	length of burn-in period. Defaults to zero.
thin	thinning parameter. Defaults to the thinning parameter chosen for object.
discard	If set to true, the previous samples are discarded.
...	other arguments.

Details

The function allows the sampler to continue from the state it last stopped.

Value

The function returns an object of class mvrm.

Author(s)

Georgios Papageorgiou [email protected]

See Also

mvrm

Examples

#see \code{mvrm} example

---

Dirichlet process mixtures of joint models

Description

Fits Dirichlet process mixtures of joint response-covariate models, where the covariates are of mixed type while the discrete responses are represented utilizing continuous latent variables. See ‘Details’ section for a full model description and Papageorgiou (2018) for all technical details.

Usage

dpmj(formula, Fcdf, data, offset, sampler = "truncated", Xpred, offsetPred,
     StorageDir, ncomp, sweeps, burn, thin = 1, seed, H, Hdf, d, D,
     Alpha.xi, Beta.xi, Alpha.alpha, Beta.alpha, Trunc.alpha, ...)

Arguments

formula	a formula defining the response and the covariates e.g. y ~ x.
Fcdf	a description of the kernel of the response variable. Currently five options are supported: 1. "poisson", 2. "negative binomial", 3. "generalized poisson", 4. "binomial" and 5. "beta binomial". The first three kernels are used for count data analysis, where the third kernel allows for both over- and under-dispersion relative to the Poisson distribution. The last two kernels are used for binomial data analysis. See ‘Details’ section for some of the kernel details.
data	an optional data frame, list or environment (or object coercible by ‘as.data.frame’ to a data frame) containing the variables in the model. If not found in ‘data’, the variables are taken from ‘environment(formula)’.
offset	this can be used to specify an a priori known component to be included in the model. This should be ‘NULL’ or a numeric vector of length equal to the sample size. One ‘offset’ term can be included in the formula, and if more are required, their sum should be used.
sampler	the MCMC algorithm to be utilized. The two options are sampler = "slice" which implements a slice sampler (Walker, 2007; Papaspiliopoulos, 2008) and sampler = "truncated" which proceeds by truncating the countable mixture at ncomp components (see argument ncomp).
Xpred	an optional design matrix the rows of which include the values of the covariates $x$ for which the conditional distribution of $Y|x,D$ (where $D$ denotes the data) is calculated. These are treated as ‘new’ covariates i.e. they do not contribute to the likelihood. The matrix shouldn't include a column of 1's. NA's can be included to obtain averaged effects.
offsetPred	the offset term associated with the new covariates Xpred. It is of dimension one i.e. the same offset term is used for all rows of Xpred. If Fcdf is one of "poisson" or "negative binomial" or "generalized poisson", then offsetPred is the Poisson offset term. If Fcdf is one of "binomial" or "beta binomial", then offsetPred is the number of Binomial trials. If offsetPred is missing, it is taken to be the mean of offset, rounded to the nearest integer.
StorageDir	a directory to store files with the posterior samples of models parameters and other quantities of interest. If a directory is not provided, files are created in the current directory and removed when the sampler completes.
ncomp	number of mixture components. It defines where the countable mixture of densities [in (1) below] is truncated. Even if sampler="slice" is chosen, ncomp needs to be specified as it is used in the initialization process.
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.
H	optional scale matrix of the Wishart-like prior assigned to the restricted covariance matrices $\Sigma_h^*$. See ‘Details’ section.
Hdf	optional degrees of freedom of the prior Wishart-like prior assigned to the restricted covariance matrices $\Sigma_h^*$. See ‘Details’ section.
d	optional prior mean of the mean vector $\mu_h$. See ‘Details’ section.
D	optional prior covariance matrix of the mean vector $\mu_h$. See ‘Details’ section.
Alpha.xi	an optional parameter that depends on the specified Fcdf argument. If Fcdf = "poisson", this argument is parameter $\alpha_{\xi}$ of the prior of the Poisson rate: $\xi \sim$ Gamma($\alpha_{\xi},\beta_{\xi}$). If Fcdf = "negative binomial", this argument is a two-dimensional vector that includes parameters $\alpha_{1\xi}$ and $\alpha_{2\xi}$ of the priors: $\xi_1 \sim$ Gamma($\alpha_{1\xi},\beta_{1\xi}$) and $\xi_2 \sim$ Gamma($\alpha_{2\xi},\beta_{2\xi}$), where $\xi_1$ and $\xi_2$ are the two parameters of the Negative Binomial pmf. If Fcdf = "generalized poisson", this argument is a two-dimensional vector that includes parameters $\alpha_{1\xi}$ and $\alpha_{2\xi}$ of the priors: $\xi_1 \sim$ Gamma($\alpha_{1\xi},\beta_{1\xi}$) and $\xi_2 \sim$ N($\alpha_{2\xi},\beta_{2\xi})I[\xi_2 \in R_{\xi_2}]$, where $\xi_1$ and $\xi_2$ are the two parameters of the Generalized Poisson pmf. Parameter $\xi_2$ is restricted in the range $R_{\xi_2} = (0.05,\infty)$ as it is a dispersion parameter. If Fcdf = "binomial", this argument is parameter $\alpha_{\xi}$ of the prior of the Binomial probability: $\xi \sim$ Beta($\alpha_{\xi},\beta_{\xi}$). If Fcdf = "beta binomial", this argument is a two-dimensional vector that includes parameters $\alpha_{1\xi}$ and $\alpha_{2\xi}$ of the priors: $\xi_1 \sim$ Gamma($\alpha_{1\xi},\beta_{1\xi}$) and $\xi_2 \sim$ Gamma($\alpha_{2\xi},\beta_{2\xi}$), where $\xi_1$ and $\xi_2$ are the two parameters of the Beta Binomial pmf. See ‘Details’ section.
Beta.xi	an optional parameter that depends on the specified family. If Fcdf = "poisson", this argument is parameter $\beta_{\xi}$ of the prior of the Poisson rate: $\xi \sim$ Gamma($\alpha_{\xi},\beta_{\xi}$). If Fcdf = "negative binomial", this argument is a two-dimensional vector that includes parameters $\beta_{1\xi}$ and $\beta_{2\xi}$ of the priors: $\xi_1 \sim$ Gamma($\alpha_{1\xi},\beta_{1\xi}$) and $\xi_2 \sim$ Gamma($\alpha_{2\xi},\beta_{2\xi}$), where $\xi_1$ and $\xi_2$ are the two parameters of the Negative Binomial pmf. If Fcdf = "generalized poisson", this argument is a two-dimensional vector that includes parameters $\beta_{1\xi}$ and $\beta_{2\xi}$ of the priors: $\xi_1 \sim$ Gamma($\alpha_{1\xi},\beta_{1\xi}$) and $\xi_2 \sim$ Normal($\alpha_{2\xi},\beta_{2\xi})I[\xi_2 \in R_{\xi_2}]$, where $\xi_1$ and $\xi_2$ are the two parameters of the Generalized Poisson pmf. Parameter $\xi_2$ is restricted in the range $R_{\xi_2} = (0.05,\infty)$ as it is a dispersion parameter. Note that $\beta_{2\xi}$ is a standard deviation. If Fcdf = "binomial", this argument is parameter $\beta_{\xi}$ of the prior of the Binomial probability: $\xi \sim$ Beta($\alpha_{\xi},\beta_{\xi}$). If Fcdf = "beta binomial", this argument is a two-dimensional vector that includes parameters $\beta_{1\xi}$ and $\beta_{2\xi}$ of the priors: $\xi_1 \sim$ Gamma($\alpha_{1\xi},\beta_{1\xi}$) and $\xi_2 \sim$ Gamma($\alpha_{2\xi},\beta_{2\xi}$), where $\xi_1$ and $\xi_2$ are the two parameters of the Beta Binomial pmf. See ‘Details’ section.
Alpha.alpha	optional shape parameter $\alpha_{\alpha}$ of the Gamma prior assigned to the concentration parameter $\alpha$. See ‘Details’ section.
Beta.alpha	optional rate parameter $\beta_{\alpha}$ of the Gamma prior assigned to concentration parameter $\alpha$. See ‘Details’ section.
Trunc.alpha	optional truncation point $c_{\alpha}$ of the Gamma prior assigned to concentration parameter $\alpha$. See ‘Details’ section.
...	Other options that will be ignored.

Details

Function dpmj returns samples from the posterior distributions of the parameters of the model:

$f(y_i,x_i) = \sum_{h=1}^{\infty} \pi_h f(y_i,x_i|\theta_h), \hspace{200pt} (1)$

where $y_i$ is a univariate discrete response, $x_i$ is a $p$-dimensional vector of mixed type covariates, and $\pi_h, h \geq 1,$ are obtained according to Sethuraman's (1994) stick-breaking construction: $\pi_1 = v_1$, and for $l \geq 2, \pi_l = v_l \prod_{j=1}^{l-1} (1-v_j)$, where $v_k$ are iid samples $v_k \sim$Beta $(1,\alpha), k \geq 1.$

Let $Z$ denote a discrete variable (response or covariate). It is represented as discretized version of a continuous latent variable $Z^*$. Observed discrete $Z$ and continuous latent variable $Z^*$ are connected by:

$z = q \iff c_{q-1} < z^* < c_{q}, q=0,1,2,\dots,$

where the cut-points are obtained as: $c_{-1} = -\infty$, while for $q \geq 0$, $c_{q} = c_{q}(\lambda) = \Phi^{-1}\{F(q;\lambda)\}.$ Here $\Phi(.)$ is the cumulative distribution function (cdf) of a standard normal variable and $F()$ denotes an appropriate cdf. Further, latent variables are assumed to independently follow a $N(0,1)$ distribution, where the mean and variance are restricted to be zero and one as they are non-identifiable by the data. Choices for $F()$ are described next.

For counts, three options are supported. First, $F(.;\lambda_i)$ can be specified as the cdf of a Poisson$(H_i \xi_h)$ variable. Here $\lambda_i=(\xi_h,H_i)^T, \xi_h$ denotes the Poisson rate associated with cluster $h$, and $H_i$ the offset term associated with sampling unit $i$. Second, $F(.;\lambda_i)$ can be specified as the negative binomial cdf, where $\lambda_i= (\xi_{1h},\xi_{2h},H_i)^T$. This option allows for overdispersion within each cluster relative to the Poisson distribution. Third, $F(.;\lambda_i)$ can be specified as the Generalized Poisson cdf, where, again, $\lambda_i=(\xi_{1h},\xi_{2h},H_i)^T$. This option allows for both over- and under-dispersion within each cluster.

For Binomial data, two options are supported. First, $F(.;\lambda_i)$ may be taken to be the cdf of a Binomial$(H_i,\xi_h)$ variable, where $\xi_h$ denotes the success probability of cluster $h$ and $H_i$ the number of trials associated with sampling unit $i$. Second, $F(.;\lambda_i)$ may be specified to be the beta-binomial cdf, where $\lambda=(\xi_{1h},\xi_{2h},H_i)^T$.

The special case of Binomial data is treated as

$Z = 0 \iff z^* < 0, z^* \sim N(\mu_z^{*},1).$

Details on all kernels are provided in the two tables below. The first table provides the probability mass functions and the mean in the presence of an offset term (which may be taken to be one). The column ‘Sample’ indicates for which parameters the routine provides posterior samples. The second table provides information on the assumed priors along with the default values of the parameters of the prior distributions and it also indicates the function arguments that allow the user to alter these.

Kernel	PMF	Offset	Mean	Sample
Poisson	$\exp(-H\xi) (H\xi)^y /y!$	$H$	$H \xi$	$\xi$
Negative Binomial	$\frac{\Gamma(y+\xi_1)}{\Gamma(\xi_1)\Gamma(y+1)}(\frac{\xi_2}{H+\xi_2})^{\xi_1}(\frac{H}{H+\xi_2})^{y}$	$H$	$H \xi_1/\xi_2$	$\xi_1, \xi_2$
Generalized Poisson	$\xi_1 \{\xi_1+(\xi_2-1)y\}^{y-1} \xi_2^{-y} \times$	$H$	$H\xi_1$	$\xi_1,\xi_2$
	$~~ \exp\{-[\xi_1+(\xi_2-1)y]/\xi_2\}/y!$
Binomial	${N \choose y} \xi^y (1-\xi)^{N-y}$	$N$	$N \xi$	$\xi$
Beta Binomial	${N \choose y} \frac{{Beta}{(y+\xi_1,N-y+\xi_2)}}{{Beta}{(\xi_1,\xi_2)}}$	$N$	$N \xi_1/(\xi_1+\xi_2)$	$\xi_1,\xi_2$

Kernel	Priors	Default Values
Poisson	$\xi \sim$ Gamma$(\alpha_{\xi},\beta_{\xi})$	Alpha.xi = 1.0, Beta.xi = 0.1
Negative Binomial	$\xi_i \sim$ Gamma$(\alpha_{\xi_i},\beta_{\xi_i}), i=1,2$	Alpha.xi = c(1.0,1.0), Beta.xi = c(0.1,0.1)
Generalized Poisson	$\xi_1 \sim$ Gamma$(\alpha_{\xi_1},\beta_{\xi_1})$
	$\xi_2 \sim$ N$(\alpha_{\xi_2},\beta_{\xi_2})I[\xi_2 > 0.05]$	Alpha.xi = c(1.0,1.0), Beta.xi = c(0.1,1.0)
	where $\beta_{\xi_2}$ denotes st.dev.
Binomial	$\xi \sim$ Beta$(\alpha_{\xi},\beta_{\xi})$	Alpha.xi = 1.0, Beta.xi = 1.0
Beta Binomial	$\xi_i \sim$ Gamma$(\alpha_{\xi_i},\beta_{\xi_i}), i=1,2$	Alpha.xi = c(1.0,1.0), Beta.xi = c(0.1,0.1)

Let $z_i = (y_i,x_{i}^T)^T$ denote the joint vector of observed continuous and discrete variables and $z_i^*$ the corresponding vector of continuous observed and latent variables. With $\theta_h$ denoting model parameters associated with the $h$th cluster, the joint density $f(z_{i}|\theta_h)$ takes the form

$f(z_i|\theta_h) = \int_{R(y)} \int_{R(x_{d})} N_{q}(z_i^*;\mu^*_h,\Sigma^*_h) dx_{d}^{*} dy^{*},$

where

$\begin{array}{ll} \mu^*_h = \left( \begin{array}{l} 0 \\ \mu_h \\ \end{array} \right), & \Sigma^*_h=\left[ \begin{array}{ll} C_h & \nu_h^T \\ \nu_h & \Sigma_h \\ \end{array} \right] \end{array},$

where $C_h$ is the covariance matrix of the latent continuous variables and it has diagonal elements equal to one i.e. it is a correlation matrix.

In addition to the priors defined in the table above, we specify the following:

1. The restricted covariance matrix $\Sigma^*_h$ is assigned a prior distribution that is based on the Wishart distribution with degrees of freedom set by default to dimension of matrix plus two and diagonal scale matrix, with the sub-matrix that corresponds to discrete variables taken to be the identity matrix and with sub-matrix that corresponds to continuous variables having entries equal to 1/8 of the square of the observed data range. Default values can be changed using arguments H and Hdf.

2. The prior on $\mu_h$, the non-zero part of $\mu_h^*$, is taken to be multivariate normal $\mu_h \sim N(d,D)$. The mean $d$ is taken to be equal to the center of the dataset. The covariance matrix $D$ is taken to be diagonal. Its elements that correspond to continuous variables are set equal to 1/8 of the square of the observed data range while the elements that correspond to binary variables are set equal to 5. Arguments Mu.mu and Sigma.mu allow the user to change the default values.

3. The concentration parameter $\alpha$ is assigned a Gamma$(\alpha_{\alpha},\beta_{\alpha})$ prior over the range $(c_{\alpha},\infty)$, that is, $f(\alpha) \propto \alpha^{\alpha_{\alpha}-1} \exp\{-\alpha \beta_{\alpha}\} I[\alpha > c_{\alpha}]$, where $I[.]$ is the indicator function. The default values are $\alpha_{\alpha}=2.0, \beta_{\alpha}=5.0$, and $c_{\alpha}=0.25$. Users can alter the default using using arguments Alpha.alpha, Beta.alpha and Turnc.alpha.

Value

Function dpmj returns the following:

call	the matched call.
seed	the seed that was used (in case replication of the results is needed).
meanReg	if Xpred is specified, the function returns the posterior mean of the conditional expectation of the response $y$ given each new covariate $x$.
medianReg	if Xpred is specified, the function returns the posterior mean of the conditional 50% quantile of the response $y$ given each new covariate $x$.
q1Reg	if Xpred is specified, the function returns the posterior mean of the conditional 25% quantile of the response $y$ given each new covariate $x$.
q3Reg	if Xpred is specified, the function returns the posterior mean of the conditional 75% quantile of the response $y$ given each new covariate $x$.
modeReg	if Xpred is specified, the function returns the posterior mean of the conditional mode of the response $y$ given each new covariate $x$.
denReg	if Xpred is specified, the function returns the posterior mean conditional density of the response $y$ given each new covariate $x$. Results are presented in a matrix the rows of which correspond to the different $x$s.
denVar	if Xpred is specified, the function returns the posterior variance of the conditional density of the response $y$ given each new covariate $x$. Results are presented in a matrix the rows of which correspond to the different $x$s.

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

alpha.txt	this file contains samples from the posterior of the concentration parameters $\alpha$. The file is arranged in (sweeps-burn)/thin lines and one column, each line including one posterior sample.
compAlloc.txt	this file contains the allocations to clusters obtained during posterior sampling. It consists of (sweeps-burn)/thin lines, that represent the posterior samples, and $n$ columns, that represent the sampling units. Clusters are represented by integers ranging from 0 to ncomp-1.
MeanReg.txt	this file contains the conditional means of the response $y$ given covariates $x$ obtained during posterior sampling. The rows represent the (sweeps-burn)/thin posterior samples. The columns represent the various covariate values $x$ for which the means are obtained.
MedianReg.txt	this file contains the 50% conditional quantile of the response $y$ given covariates $x$ obtained during posterior sampling. The rows represent the (sweeps-burn)/thin posterior samples. The columns represent the various covariate values $x$ for which the medians are obtained.
muh.txt	this file contains samples from the posteriors of the $p$-dimensional mean vectors $\mu_h, h=1,2,\dots$,ncomp. The file is arranged in ((sweeps-burn)/thin)*ncomp lines and $p$ columns. In more detail, sweeps create ncomp lines representing samples $\mu_h^{(sw)}, h=1,\dots,$ncomp, where superscript $sw$ represents a particular sweep. The elements of $\mu_h^{(sw)}$ are written in the columns of the file.
nmembers.txt	this file contains (sweeps-burn)/thin lines and ncomp columns, where the lines represent posterior samples while the columns represent the components or clusters. The entries represent the number of sampling units allocated to each component.
Q05Reg.txt	this file contains the 5% conditional quantile of the response $y$ given covariates $x$ obtained during posterior sampling. The rows represent the (sweeps-burn)/thin posterior samples. The columns represent the various covariate values $x$ for which the quantiles are obtained.
Q10Reg.txt	as above, for the 10% conditional quantile.
Q15Reg.txt	as above, for the 15% conditional quantile.
Q20Reg.txt	as above, for the 20% conditional quantile.
Q25Reg.txt	as above, for the 25% conditional quantile.
Q75Reg.txt	as above, for the 75% conditional quantile.
Q80Reg.txt	as above, for the 80% conditional quantile.
Q85Reg.txt	as above, for the 85% conditional quantile.
Q90Reg.txt	as above, for the 90% conditional quantile.
Q95Reg.txt	as above, for the 95% conditional quantile.
Sigmah.txt	this file contains samples from the posteriors of the $q \times q$ restricted covariance matrices $\Sigma_h^*, h=1,2,\dots,$ncomp. The file is arranged in ((sweeps-burn)/thin)*ncomp lines and $q^2$ columns. In more detail, sweeps create ncomp lines representing samples $\Sigma_h^{(sw)}, h=1,\dots,$ncomp, where superscript $sw$ represents a particular sweep. The elements of $\Sigma_h^{(sw)}$ are written in the columns of the file.
xih.txt	this file contains samples from the posteriors of parameters $\xi_h$, $h=1,2,\dots,$ncomp. The file is arranged in ((sweeps-burn)/thin)*ncomp lines and one or two columns, depending on the number of parameters in the selected Fcdf. Sweeps write in the file ncomp lines representing samples $\xi_h^{(sw)}, h=1,\dots,$ncomp, where superscript $sw$ represents a particular sweep.
Updated.txt	this file contains (sweeps-burn)/thin lines with the number of components updated at each iteration of the sampler (relevant for slice sampling).

Author(s)

Georgios Papageorgiou [email protected]

References

Consul, P. C. & Famoye, G. C. (1992). Generalized Poisson regression model. Communications in Statistics - Theory and Methods, 1992, 89-109.

Papageorgiou, G. (2018). Bayesian density regression for discrete outcomes. arXiv:1603.09706v3 [stat.ME].

Papaspiliopoulos, O. (2008). A note on posterior sampling from Dirichlet mixture models. Technical report, University of Warwick.

Sethuraman, J. (1994). A constructive definition of Dirichlet priors. Statistica Sinica, 4, 639-650.

Walker, S. G. (2007). Sampling the Dirichlet mixture model with slices. Communications in Statistics Simulation and Computation, 36(1), 45-54.

Examples

#Bayesian nonparametric joint model with binomial response Y and one predictor X
data(simD)
pred<-seq(with(simD,min(X))+0.1,with(simD,max(X))-0.1,length.out=30)
npred<-length(pred)
# fit1 and fit2 define the same model but with different numbers of
# components and posterior samples
fit1 <- dpmj(cbind(Y,(E-Y))~X, Fcdf="binomial", data=simD, ncomp=10, sweeps=20,
             burn=10, sampler="truncated", Xpred=pred, offsetPred=30)
fit2 <- dpmj(cbind(Y,(E-Y))~X, Fcdf="binomial", data=simD, ncomp=50, sweeps=5000,
               burn=1000, sampler="truncated", Xpred=pred, offsetPred=30)
plot(with(simD,X),with(simD,Y)/with(simD,E))
lines(pred,fit2$medianReg/30,col=3,lwd=2)
# with discrete covariate
simD<-data.frame(simD,Xd=sample(c(0,1),300,replace=TRUE))
pred<-c(0,1)
fit3 <- dpmj(cbind(Y,(E-Y))~Xd, Fcdf="binomial", data=simD, ncomp=10, sweeps=20,
             burn=10, sampler="truncated", Xpred=pred, offsetPred=30)

---

Creates plots of correlation matrices

Description

This function plots the posterior distribution of the elements of correlation matrices.

Usage

histCorr(x, term = "R", plotOptions = list(),...)

Arguments

x	an object of class ‘mvrm’, as generated by function mvrm.
term	Admits two possible values: "R" to plot samples from the posterior of the correlation matrix $R$, and "muR" to plot samples from the posterior of the means $\mu_R$.
plotOptions	ggplot type options.
...	other arguments.

Details

Use this function to visualize the elements of a correlation matrix.

Value

Posterior distributions of elements of correlation matrices.

Author(s)

Georgios Papageorgiou [email protected]

See Also

mvrm

Examples

#see \code{mvrm} example

---

Bayesian semiparametric modelling of covariance matrices for multivariate longitudinal data

Description

Implements an MCMC algorithm for posterior sampling based on a semiparametric model for continuous longitudinal multivariate responses. The overall model consists of 5 regression submodels and it utilizes spike-slab priors for variable selection and function regularization. See ‘Details’ section for a full description of the model.

Usage

lmrm(formula, data = list(), centre=TRUE, id, time,
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)", pi.phiPrior = "Beta(1,1)", c.psiPrior = "HN(2)",
sigmaPrior = "HN(2)", pi.sigmaPrior = "Beta(1,1)", 
corr.Model = c("common", nClust = 1), DP.concPrior = "Gamma(5,2)",
c.etaPrior = "IG(0.5,0.5*samp)", pi.nuPrior = "Beta(1,1)", 
pi.fiPrior = "Beta(1,1)", c.omegaPrior = "IG(1.1,1.1)", sigmaCorPrior = "HN(2)",
tuneCalpha, tuneSigma2, tuneCbeta, tuneAlpha, tuneSigma2R, tuneR, tuneCpsi, 
tuneCbCor, tuneOmega, tuneComega, tau, FT = 1,...)

Arguments

formula	a formula defining the responses and the covariates in the 5 regression models e.g. y1 | y2 ~ x | w | z | t | t or for smooth effects y1 | y2 ~ sm(x) | sm(w) | sm(z) | sm(t) | sm(t). The package uses the extended formula notation, where the responses are defined on the left of ~ and the models on the right.
data	a data frame.
centre	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.
id	identifiers of the individuals or other sampling units that are observed over time.
time	a vector input that specifies the time of observation
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_{\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.
pi.muPrior	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 $p K$
c.alphaPrior	The inverse Gamma prior of $c_{\alpha}^2$. The default is "IG(1.1,1.1)". Half-normal priors for $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).
pi.phiPrior	The Beta prior of $\pi_{\phi}$. The default is "Beta(1,1)". It can be of dimension $1$, of dimension $B$ (the number of effects that enter the dependence model), or of dimension $p^2 B$
c.psiPrior	The prior of $c_{\psi}^2$. The default is "HN(2)", a half-normal prior for $c_{\psi}$ with variance equal to two, $c_{\psi} \sim N(0,2) I[c_{\psi} > 0]$. Inverse Gamma priors for $c_{\psi}^2$ are also available, declared using "IG(a,b)". It can be of dimension $1$ or $p^2$ (the number of dependence models).
sigmaPrior	The prior of $\sigma_k^2, k=1,\dots,p$. The default is "HN(2)", a half-normal prior for $\sigma_k$ with variance equal to two, $\sigma_k \sim N(0,2) I[\sigma>0]$. Inverse Gamma priors for $\sigma_k^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 $\pi_{\sigma}$. The default is "Beta(1,1)". It can be of dimension $1$, of dimension $L$ (the number of effects that enter the variance model), or of dimension $pL$
corr.Model	Specifies the model for the correlation matrices $R_t$. 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 matrices.
DP.concPrior	The Gamma prior for the Dirichlet process concentration parameter.
c.etaPrior	The inverse Gamma prior of $c_{\eta}$. The default is "IG(0.5,0.5*samp)", that is, an inverse Gamma with parameters $1/2$ and $samp/2$, where $samp$ is the number of correlations observed over time, that is $samp=M*d$ where $M$ is the number of unique observation time points and $d$ is the number of non-redundant elements of $R$.
pi.nuPrior	The Beta prior of $\pi_{\nu}$. The default is "Beta(1,1)". It can be of dimension $1$.
pi.fiPrior	The Beta prior of $\pi_{\varphi}$. The default is "Beta(1,1)". It can be of dimension $1$.
c.omegaPrior	The prior of $c_{\omega}^2$. The default is "HN(2)", a half-normal prior for $c_{\omega}$ with variance equal to two, $c_{\omega} \sim N(0,2) I[c_{\omega} > 0]$. Inverse Gamma priors for $c_{\omega}^2$ are also available, declared using "IG(a,b)". It can be of dimension $1$.
sigmaCorPrior	The prior of $\sigma^2$. The default is "HN(2)", a half-normal prior for $\sigma^2$ 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$.
tuneCalpha	Starting value of the tuning parameter for sampling $c_{\alpha k}, k=1,\dots,p$. Defaults at a vector of $p$ ones. It could be of dimension $p$.
tuneSigma2	Starting value of the tuning parameter for sampling $\sigma^2_k, k=1,\dots,p$. Defaults at a vector of $p$ ones. It could be of dimension $p$.
tuneCbeta	Starting value of the tuning parameter for sampling $c_{\beta}$. Defaults at 100.
tuneAlpha	Starting value of the tuning parameter for sampling regression coefficients of the variance models $\alpha_k, k=1,\dots,p$. Defaults at a vector of 5s. It could be of dimension $L p$
tuneSigma2R	Starting value of the tuning parameter for sampling $\sigma^2$. Defaults at 1.
tuneR	Starting value of the tuning parameter for sampling correlation matrices. Defaults at $40*(p+2)^3$. Can be of dimension $1$ or $M$ is the number of unique observation time points.
tuneCpsi	Starting value of the tuning parameter for sampling variances $c_{\psi}^2$. Defaults at 5. Can be of dimension $1$ or $p^2$

.

tuneCbCor	Starting value of the tuning parameter for sampling $c_{\eta}^2$. Defaults at 10.
tuneOmega	Starting value of the tuning parameter for sampling regression coefficients of the variance models $\omega$. Defaults at 5.
tuneComega	Starting value of the tuning parameter for sampling $c_{\omega}$. 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 lmrm returns samples from the posterior distributions of the parameters of a regression model with normally distributed multivariate longitudinal responses. To describe the model, let $Y_{ij} = (Y_{ij1},\dots,Y_{ijp})^{\top}$ denote the vector of $p$ responses observed on individual $i, i=1,\dots,n,$ at time point $t_{ij}, j=1,\dots,n_i$. The observational time points $t_{ij}$ are allowed to be unequally spaced. Further, let $u_{ij}$ denote the covariate vector that is observed along with $Y_{ij}$ and that may include time, other time-dependent covariates and time-independent ones. In addition, let $Y_{i}=(Y_{i1}^{\top},\dots,Y_{in_i}^{\top})^{\top}$ denote the $i$th response vector. With $\mu_i=E(Y_{i})$ and $\Sigma_i =$ cov$(Y_i)$, the model assumes multivariate normality, $Y_{i} \sim N(\mu_i, \Sigma_i), i=1,2,\dots,n$. The means $\mu_i$ and covariance matrices $\Sigma_i$ are modelled semiparametrically in terms of covariates. For the means one can specify semiparametric models,

$\mu_{ijk} = \beta_{k0} + \sum_{l=1}^{K_1} u_{ijl} \beta_{kl} + \sum_{l=K_1+1}^{K} f_{\mu,k,l}(u_{ijl}).$

This is the first of the 5 regression submodels.

To model the covariance matrix, first consider the modified Cholesky decomposition, $L_i \Sigma_i L_i^{\top} = D_i,$ where $L_i$ is a unit block lower triangular matrix and $D_i$ is a block diagonal matrix,

$\begin{array}{cc} L_i = \left[ \begin{array}{cccc} I & 0 & \dots & 0 \\ -\Phi_{i21} & I & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ -\Phi_{in_i1} & -\Phi_{in_i1} & \dots & I \\ \end{array} \right], & D_i = \left[ \begin{array}{cccc} D_{1} & 0 & \dots & 0 \\ 0 & D_{2} & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & D_{n_i} \\ \end{array} \right], \end{array}$

For modelling $D_{ij}, i=1,\dots,n, j=1,\dots,n_i$ in terms of covariates, first we separate the variances and the correlations $D_{ij} = S_{ij}^{1/2} R_{ij} S_{ij}^{1/2}$. It is easy to model matrix $S_{ij}$ in terms of covariates as the only requirement on its diagonal elements is that they are nonnegative,

$\log \sigma^2_{ijk} = \alpha_{k0} + \sum_{l=1}^{L_1} w_{ijl} \alpha_{kl} + \sum_{l=L_1+1}^{L} f_{\sigma,k,l}(w_{ijl})$

This is the second of the 5 regression submodels.

For $\phi_{ijklm}$, the $(l,m)$ element of $\Phi_{ijk}, l,m=1,\dots,p$, one can specify semiparametric models

$\phi_{ijklm} = \psi_{lm0} + \sum_{b=1}^{B_1} v_{ijkb} \psi_{lmb} + \sum_{b=B_1+1}^{B} f_{\phi,l,m,b}(v_{ijkb})$

This is the third of the 5 regression submodels.

The elements of the correlations matrices $R_{ij}$ are modelled in terms of covariate time only, hence they are denoted by $R_t$. Subject to the positive definiteness constraint, the elements of $R_t$ are modelled using a normal distribution with location and scale parameters, $\mu_{ct}$ and $\sigma^2_{ct}$, modelled as

$\mu_{ct} = \eta_0 + f_{\mu}(t),$

$\log \sigma^2_{ct} = \omega_0 + f_{\sigma}(t),$

and these are the last 2 of the 5 submodels.

Value

Function lmrm returns samples from the posteriors of the model parameters.

Author(s)

Georgios Papageorgiou [email protected]

References

Papageorgiou, G. (2020). Bayesian semiparametric modelling of covariance matrices for multivariate longitudinal data. arXiv:2012.09833.

Examples

# Fit a joint mean-covariance model on the simulated dataset simD2
	require(ggplot2)
	data(simD2)
	model <- Y1 | Y2 ~ time | sm(time) | sm(lag) | sm(time) | 1 
	# the above defines the responses and the regression models on the left and 
	# right of "~", respectively 
	# the first model, for the mean, is a linear function of time, this is sufficient as 
	# the 2 responses have constant mean.
	# the second model, for the variances, is a smooth function of time
	# the third model, for the dependence structure, is a smooth function of lag, 
	# that lmrm figures out and it does not need to be computed by the user
	# the fourth model, for location of the correlations, is a smooth function of time
	# the fifth model, for scale of the correlations, is just an intercept model
	## Not run: 
	m1 <- lmrm(formula = model, corr.Model = c("common", nClust = 1), data = simD2,
		       id = id, time = time, sweeps = 2500, burn = 500, thin = 2, 
		       StorageDir = getwd(), seed = 1)
	plot(m1)
	
## End(Not run)

---

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, 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, ...)

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.
distribution	The distribution for the response variables. Currently two options are supported: "normal" and "t".
data	a data frame.
centre	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.
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_{\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.
pi.muPrior	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$
c.alphaPrior	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).
sigmaPrior	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).
pi.sigmaPrior	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$
c.psiPrior	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).
phiPrior	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).
pi.omegaPrior	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$
mu.RPrior	The normal prior for $\mu_{R}$. The default is the standard normal distribution.
sigma.RPrior	The half normal prior for $\sigma_{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.
breaksPrior	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.
tuneCbeta	Starting value of the tuning parameter for sampling $c_{\beta}$. Defaults at 20.
tuneCalpha	Starting value of the tuning parameter for sampling $c_{\alpha}$. Defaults at 1.
tuneAlpha	Starting value of the tuning parameter for sampling regression coefficients of the variance model $\alpha$. Defaults at 5.
tuneSigma2	Starting value of the tuning parameter for sampling variances $\sigma^2_j$. Defaults at 1.
tuneCpsi	Starting value of the tuning parameter for sampling $c_{\psi}$. Defaults at 1.
tunePsi	Starting value of the tuning parameter for sampling $\psi$. Defaults at 5.
tunePhi	Starting value of the tuning parameter for sampling $\varphi$. Defaults at 0.5.
tuneR	Starting value of the tuning parameter for sampling correlation matrices. Defaults at $40*(p+2)^3$.
tuneSigma2R	Starting value of the tuning parameter for sampling $\sigma_{R}^2$. Defaults at 1.
tuneHar	Starting value of the tuning parameter for sampling the regression coefficients of the harmonics. Defaults at 100.
tuneBreaks	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.
tunePeriod	Starting value of the tuning parameter for sampling the period parameter of the sin curve. Defaults to 0.01.
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.
compDeviance	Binary indicator. If set equal to 1, the deviance is computed.
...	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

$\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.

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 $\gamma$.
LD	the length of the vector of indicators $\delta$.
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 $\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.
beta.txt	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.
gamma.txt	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.
delta.txt	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.
sigma2.txt	contains samples from the posterior of the error variance $\sigma^2$ of each response.
cbeta.txt	contains samples from the posterior of $c_{\beta}$.
calpha.txt	contains samples from the posterior of $c_{\alpha}$.
R.txt	contains samples from the posterior of the correlation matrix $R$.
theta.txt	contains samples from the posterior of $\theta$ of the shadow prior (probably not needed).
muR.txt	contains samples from the posterior of $\mu_R$.
sigma2R.txt	contains samples from the posterior of $\sigma^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, $\lambda_{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, $\lambda_{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 [email protected]

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)

---