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

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 <gpapageo@gmail.com>

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

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 gpapageo@gmail.com

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

Details

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

Value

Similarity matrix.

Author(s)

Georgios Papageorgiou gpapageo@gmail.com

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

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 gpapageo@gmail.com

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

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 \simBeta (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.

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 hth 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:

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

Author(s)

Georgios Papageorgiou gpapageo@gmail.com

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

Details

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

Value

Posterior distributions of elements of correlation matrices.

Author(s)

Georgios Papageorgiou gpapageo@gmail.com

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

.

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 ith 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 gpapageo@gmail.com

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

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:

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

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

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

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)

Convert posterior samples from function mvrm into an object of class ‘mcmc’

Description

Reads in files where the posterior samples were written and creates an object of class ‘mcmc’ so that functions like summary and plot from package coda can be used

Usage

mvrm2mcmc(mvrmObj, labels)

Arguments

Value

An object of class ‘mcmc’ that holds the samples from the posterior of the selected parameter.

Author(s)

Georgios Papageorgiou gpapageo@gmail.com

See Also

mvrm

Examples

#see \code{mvrm} example

Creates plots of terms in the mean and/or variance models

Description

This function plots estimated terms that appear in the mean and variance models.

Usage

## S3 method for class 'mvrm'
plot(x, model, term, response, response2, intercept = TRUE, grid = 30, 
centre = mean, quantiles = c(0.1, 0.9), contour = TRUE, static = TRUE, 
centreEffects = FALSE, plotOptions = list(), nrow, ask = FALSE, 
plotEmptyCluster = FALSE, combine = FALSE, ...) 

Arguments

Details

Use this function to obtain predictions.

Value

Predictions along with credible/pediction intervals

Author(s)

Georgios Papageorgiou gpapageo@gmail.com

See Also

mvrm

Examples

#see \code{mvrm} example

Creates plots of the correlation matrices

Description

This function plots the posterior mean and credible intervals of the elements of correlation matrices.

Usage

plotCorr(x, term = "R", centre = mean, quantiles = c(0.1, 0.9), ...)

Arguments

Details

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

Value

Posterior means and credible intervals of elements of correlation matrices.

Author(s)

Georgios Papageorgiou gpapageo@gmail.com

See Also

mvrm

Examples

#see \code{mvrm} example

Model predictions

Description

Provides predictions and posterior credible/prediction intervals for given feature vectors.

Usage

## S3 method for class 'mvrm'
predict(object, newdata, interval = c("none", "credible", "prediction"), 
                       level = 0.95, ind.preds=FALSE, ...)

Arguments

Details

The function returns predictions of new responses or the means of the responses for given feature vectors. Predictions for new responses or the means of new responses are the same. However, the two differ in the associated level of uncertainty: response predictions are associated with wider (prediction) intervals than mean response predictions. To obtain prediction intervals (for new responses) the function samples from the normal distributions with means and variances as sampled during the MCMC run.

Value

Predictions for given covariate/feature vectors.

Author(s)

Georgios Papageorgiou gpapageo@gmail.com

See Also

mvrm

Examples

#see \code{mvrm} example

Prints an mvrm fit

Description

Provides basic information from an mvrm fit.

Usage

## S3 method for class 'mvrm'
print(x, digits = 5, ...)

Arguments

Details

The function prints information about mvrm fits.

Value

The function provides a matched call, the number of posterior samples obtained and marginal inclusion probabilities of the terms in the mean and variance models.

Author(s)

Georgios Papageorgiou gpapageo@gmail.com

See Also

mvrm

Examples

#see \code{mvrm} example

mgcv constructor s

Description

Provides interface between mgcv::s and BNSP. s(...) calls mgcv::smoothCon(mgcv::s(...),...

Usage

s(..., data, knots = NULL, absorb.cons = FALSE, scale.penalty = TRUE, 
n = nrow(data), dataX = NULL, null.space.penalty = FALSE, sparse.cons = 0, 
diagonal.penalty = FALSE, apply.by = TRUE, modCon = 0, k = -1, fx = FALSE, 
bs = "tp", m = NA, by = NA, xt = NULL, id = NULL, sp = NULL, pc = NULL)

Arguments

Details

The most relevant arguments for BNSP users are the list of variables ..., knots, absorb.cons, bs, and by.

Value

A design matrix that specifies a smooth term in a model.

Author(s)

Georgios Papageorgiou gpapageo@gmail.com

Simulated dataset

Description

Just a simulated dataset to illustrate the DO mixture model. The success probability and the covariate have a non-linear relationship.

Usage

data(simD)

Format

A data frame with 300 independent observations. Three numerical vectors contain information on

Y

number of successes.

E

number of trials.

X

explanatory variable.

Simulated dataset

Description

A simulated dataset to illustrate the multivariate longitudinal model. It consists of a bivariate vector of responses observed over 6 time points.

Usage

data(simD2)

Format

A data frame that includes observations on 50 sampling units. The data frame has 300 rows for the 50 sampling units observed over 6 time points. It has 4 columns

Y1

first response.

Y2

second response.

time

the time of observation.

id

unique sampling unit identifier.

Sinusoid terms in mvrm formulae

Description

Function used to define sinusoidal curves in the mean formula of function mvrm. The function is used internally to construct design matrices.

Usage

sinusoid(..., harmonics = 1, amplitude = 1, period = 0, periodRange = NULL, 
breaks = NULL, knots = NULL)

Arguments

Details

Use this function within calls to function mvrm to specify sinusoidal waves in the mean function of a regression model.

Consider the sinusoidal curve

y_t = \beta_0 + A(t) \sin(2\pi t/\omega+\varphi) + \epsilon_t,

where y_t is the response at time t, \beta_0 is an intercept term, A(t) is a time-varying amplitude, \varphi \in [0,2\pi] is the phase shift parameter, \omega is the period taken to be known, and \epsilon_t is the error term.

The period \omega is defined by the argument period.

The time-varying amplitude is represented using A(t) = \sum_{j=1}^{K} \beta_{Aj} \phi_{Aj}(t), where K, the number of knots, is defined by argument amplitude. If amplitude = 1, then the amplitude is taken to be fixed: A(t)=A.

Further, \sin(2\pi t/\omega+\varphi) is represented utilizing \sin(2\pi t/\omega+\varphi) = \sum_{k=1}^{L} a_k \sin(2k\pi t/\omega) + b_k \cos(2k\pi t/\omega), where L, the number of harmonics, is defined by argument harmonics.

Value

Specifies the design matrices of an mvrm call

Author(s)

Georgios Papageorgiou gpapageo@gmail.com

See Also

mvrm

Examples

# Simulate and fit a sinusoidal curve
# First load releveant packages
require(BNSP)
require(ggplot2)
require(gridExtra)
require(Formula)
# Simulate the data 
mu <- function(u) {cos(0.5 * u) * sin(2 * pi * u + 1)}
set.seed(1)
n <- 100
u <- sort(runif(n, min = 0, max = 2*pi))
y <- rnorm(n, mu(u), 0.1)
data <- data.frame(y, u)
# Define the model and call function \code{mvrm} that perfomes posterior sampling for the given 
# dataset and defined model 
model <- y ~ sinusoid(u, harmonics = 2, amplitude = 20, period = 1)
## Not run: 
m1 <- mvrm(formula = model, data = data, sweeps = 10000, burn = 5000, thin = 2, seed = 1, 
           StorageDir = getwd())
# Plot
x1 <- seq(min(u), max(u), length.out = 100)
plotOptionsM <- list(geom_line(aes_string(x = x1, y = mu(x1)), col = 2, alpha = 0.5, lty = 2),
                     geom_point(data = data, aes(x = u, y = y)))
plot(x = m1, term = 1, plotOptions = plotOptionsM, intercept = TRUE, 
     quantiles = c(0.005, 0.995), grid = 100, combine = 1)

## End(Not run)

Smooth terms in mvrm formulae

Description

Function used to define smooth effects in the mean and variance formulae of function mvrm. The function is used internally to construct the design matrices.

Usage

sm(..., k = 10, knots = NULL, bs = "rd")

Arguments

Details

Use this function within calls to function mvrm to specify smooth terms in the mean and/or variance function of the regression model.

Univariate radial basis functions with q basis functions or q-1 knots are defined by

\mathcal{B}_1 = \left\{\phi_{1}(u)=u , \phi_{2}(u)=||u-\xi_{1}||^2 \log\left(||u-\xi_{1}||^2\right), \dots, \phi_{q}(u)=||u-\xi_{q-1}||^2 \log\left(||u-\xi_{q-1}||^2\right)\right\},

where ||u|| denotes the Euclidean norm of u and \xi_1,\dots,\xi_{q-1} are the knots that are chosen as the quantiles of the observed values of explanatory variable u, with \xi_1=\min(u_i), \xi_{q-1}=\max(u_i) and the remaining knots chosen as equally spaced quantiles between \xi_1 and \xi_{q-1}.

Thin plate splines are defined by

\mathcal{B}_2 = \left\{\phi_{1}(u)=u , \phi_{2}(u)=(u-\xi_{1})_{+}, \dots, \phi_{q}(u)=(u-\xi_{q})_{+}\right\},

where (a)_+ = \max(a,0).

Radial basis functions for bivariate smooths are defined by

\mathcal{B}_3 = \left\{u_1,u_2,\phi_{3}(u)=||u-\xi_{1}||^2 \log\left(||u-\xi_{1}||^2\right), \dots, \phi_{q}(u)=||u-\xi_{q-1}||^2 \log\left(||u-\xi_{q-1}||^2\right)\right\}.

Value

Specifies the design matrices of an mvrm call

Author(s)

Georgios Papageorgiou gpapageo@gmail.com

See Also

mvrm

Examples

#see \code{mvrm} example

Summary of an mvrm fit

Description

Provides basic information from an mvrm fit.

Usage

## S3 method for class 'mvrm'
summary(object, nModels = 5, digits = 5, printTuning = FALSE, ...)

Arguments

Details

Use this function to summarize mvrm fits.

Value

The functions provides a detailed description of the specified model and priors. In addition, the function provides information about the Markov chain ran (length, burn-in, thinning) and the folder where the files with posterior samples are stored. Lastly, the function provides the mean posterior and null deviance and the mean/variance models visited most often during posterior sampling.

Author(s)

Georgios Papageorgiou gpapageo@gmail.com

See Also

mvrm

Examples

#see \code{mvrm} example

mgcv constructor te

Description

Provides interface between mgcv::te and BNSP. te(...) calls mgcv::smoothCon(mgcv::te(...),...

Usage

te(..., data, knots = NULL, absorb.cons = FALSE, scale.penalty = TRUE, 
n = nrow(data), dataX = NULL, null.space.penalty = FALSE, sparse.cons = 0, 
diagonal.penalty = FALSE, apply.by = TRUE, modCon = 0, k = NA, bs = "cr", 
m = NA, d = NA, by = NA, fx = FALSE, np = TRUE, xt = NULL, id = NULL, 
sp = NULL, pc = NULL)

Arguments

Details

The most relevant arguments for BNSP users are the list of variables ..., knots, absorb.cons, bs, and by.

Value

A design matrix that specifies a smooth term in a model.

Author(s)

Georgios Papageorgiou gpapageo@gmail.com

mgcv constructor ti

Description

Provides interface between mgcv::ti and BNSP. ti(...) calls mgcv::smoothCon(mgcv::ti(...),...

Usage

ti(..., data, knots = NULL, absorb.cons = FALSE, scale.penalty = TRUE, 
n = nrow(data), dataX = NULL, null.space.penalty = FALSE, sparse.cons = 0, 
diagonal.penalty = FALSE, apply.by = TRUE, modCon = 0, k = NA, bs = "cr", 
m = NA, d = NA, by = NA, fx = FALSE, np = TRUE, xt = NULL, id = NULL, 
sp = NULL, mc = NULL, pc = NULL)

Arguments

Details

The most relevant arguments for BNSP users are the list of variables ..., knots, absorb.cons, bs, and by.

Value

A design matrix that specifies a smooth term in a model.

Author(s)

Georgios Papageorgiou gpapageo@gmail.com