mcmc.ssp: Function for sampling from the posterior under multivariate...

In maryclare/mvshrink:

Description

Function for sampling values of beta and possibly Sigma from the posterior distribution under multivariate normal, structured normal-gamma, structured power/bridge and structured product normal priors

Usage

samples <- mcmc.ssp(X = X, y = y, Sigma = diag(dim(X)[3]), sigma.sq = 1, prior = "sno", num.samp = 100)

Arguments

X	n by t by r array of covariate values, y[i] = vec(X_i)'vec(B) + z[i], i = 1,...,n
y	n by 1 response vector or m by p response matrix (if matrix, require mp=n)
Sigma	either a fixed r by r positive semidefinite covariance matrix (the provided Σ can always be used by the multivariate normal and SPN priors, but for the SNG and SPB priors, the provided value of Σ may not be achievable and the closest positive definite matrix will be used) or NULL, in which case values of Σ are resampled under either an inverse-Wishart prior for Σ^{-1}, or Σ is assumed to be diagonal with inverse-Wishart distributed elements or Σ is proportional to an identity matrix with inverse-gamma scale, the "str" parameter determines the prior used for Σ, prior hyperparameters given by the "pr.V.inv" and "pr.df" parameters
sigma.sq	positive constant, error variance
prior	string equal to "sno" for multivariate normal prior, "sng" for normal-gamma prior, "spb" for power/bridge prior, "spn" for normal product prior
c	positive constant, required when "sng" prior is used, default is NULL
q	positive constant, required when "spb" prior is used, default is NULL
Psi	positive semi-definite matrix, used when prior="spn" and Σ is provided and fixed
num.samp	integer; total number of posterior samples to return, defaults to 100
burn.in	integer; total number of posterior samples to discard, defaults to 1
thin	integer; number of posterior samples to be discarded between those returned, defaults to 0
print.iter	logical; if TRUE, the function prints a count of the number of samples from the posterior, defaults to FALSE
str	string equal to "uns", "het" or "con" which determines the prior used for Sigma. "con" forces Σ to be proportional to the identity matrix, "het" forces Σ to be diagonal, "str" allows Σ to have arbitrary structure
pr.V.inv	prior scale matrix for inverse-Wishart prior, set by default to I_r, if a diagonal Σ is used provide a diagonal matrix with the inverse-gamma rate parameters for each diagonal element of Σ, if Σ \propto I_r is used, provide the rate parameter for the inverse-gamma prior on the scale divided by r
pr.df	prior degrees of freedom for inverse-Wishart prior, set by default to p + 2, if a diagonal Σ is used provide the shape parameter for the inverse-gamma priors on the diagonal elements of Σ, if Σ \propto I_r is used provide the shape parameter for the inverse-gamma prior on the scale divided by r
pr.shape	prior shape parameter for inverse-gamma prior on the noise variance, set to 3/2 by default
pr.rate	prior rate parameter for inverse-gamma prior on the noise variance, set to 1/2 by default
W	n by l matrix of unpenalized covariates to include in regression
reg	string equal to "linear" for linear regression, "logit" for logistic regression, "nb" for negative binomial regression
pr.ssqi.V.inv	If y has matrix structure, this is the scale matrix for an inverse wishart prior on column covariance matrix
pr.ssqi.df	If y has matrix structure, this is the scale matrix for an inverse wishart on column covariance matrix
str.ssqi	Structure to be used for covariance matrix of columns of Y
joint.beta	logical; if TRUE, the function updates all elements of matrix of penalized regression coefficients jointly, otherwise the function updates the matrix of penalized regression coefficients one row at a time, set to TRUE by default
nb.r	positive scalar; overdispersion parameter for negative binomial regression, defaults to NULL in which case a gamma prior is assumed for r and values of r are resampled
pr.nb.r.shape	positive scalar; shape parameter for gamma prior on r, negative binomial overdisperson parameter. defaults to 1/2
pr.nb.r.rate	positive scalar; rate parameter for gamma prior on r, negative binomial overdispersion parameter. defaults to 1/2
eps	positive scalar or length p vector; variance for random-walk metropolis-hastings for updating log(r), defaults to 0.5

Source

Based on the material in the working paper "Structured Shrinkage Priors"

Examples

set.seed(1)
library(RColorBrewer)
# Simulate some data
n <- 30; t <- 5; r <- 4
X <- array(rnorm(n*t*r), dim = c(n, t, r))
Sigma <- (1 - 0.5)*diag(r) + 0.5*diag(r)
beta <- c(matrix(rnorm(t*r), nrow = t, ncol = r)%*%chol(Sigma))
y <- t(apply(X, 1, "c"))%*%beta + rnorm(n)

# Fit a few models to this data, first fixing Sigma then not
mvn.fs <- mcmc.ssp(X = X, y = y, Sigma = Sigma, sigma.sq = 1, prior = "sno")
spn.fs <- mcmc.ssp(X = X, y = y, Sigma = Sigma, sigma.sq = 1, prior = "spn", Psi = sqrt(abs(Sigma)))
sng.fs <- mcmc.ssp(X = X, y = y, Sigma = Sigma, sigma.sq = 1, prior = "sng", c = 1)
spb.fs <- mcmc.ssp(X = X, y = y, Sigma = Sigma, sigma.sq = 1, prior = "spb", q = 1)

# Take a look at some trace plots, compare very small sample posterior means
par(mfrow = c(2, 2))
plot(mvn.fs$betas[, 1], xlab = "", ylab = expression(beta[1]))
plot(spn.fs$betas[, 1], xlab = "", ylab = "")
plot(sng.fs$betas[, 1], xlab = "Sample", ylab = expression(beta[1]))
plot(spb.fs$betas[, 1], xlab = "Sample", ylab = "")

par(mfrow = c(1, 3))
plot(colMeans(mvn.fs$betas), colMeans(spn.fs$betas), xlab = "MVN", ylab = "SPN")
abline(a = 0, b = 1, lty = 2)
plot(colMeans(mvn.fs$betas), colMeans(sng.fs$betas), xlab = "MVN", ylab = "SNG")
abline(a = 0, b = 1, lty = 2)
plot(colMeans(mvn.fs$betas), colMeans(spb.fs$betas), xlab = "MVN", ylab = "SPB")
abline(a = 0, b = 1, lty = 2)

mvn.vs <- mcmc.ssp(X = X, y = y, Sigma = NULL, sigma.sq = NULL, prior = "sno", str = "uns")
spn.vs <- mcmc.ssp(X = X, y = y, Sigma = NULL, sigma.sq = NULL, prior = "spn", str = "uns")
sng.vs <- mcmc.ssp(X = X, y = y, Sigma = NULL, sigma.sq = NULL, prior = "sng", str = "uns", c = 1)
spb.vs <- mcmc.ssp(X = X, y = y, Sigma = NULL, sigma.sq = NULL, prior = "spb", str = "uns", q = 1)

# Plot covariance matrix posterior means
cols <- brewer.pal(11, "RdBu")
cols <- cols[length(cols):1]
breaks <- seq(-1, 1, length.out = 12)
par(mfrow = c(2, 2))
image(matrix(rowMeans(apply(mvn.vs$Sigmas, 1, function(x) {cov2cor(x)})), nrow = r, ncol = r)[r:1, r:1], breaks = breaks, col = cols, main = "MVN")
image(matrix(rowMeans(apply(spn.vs$Sigmas, 1, function(x) {cov2cor(x)})), nrow = r, ncol = r)[r:1, r:1], breaks = breaks, col = cols, main = "SPN")
image(matrix(rowMeans(apply(sng.vs$Sigmas, 1, function(x) {cov2cor(x)})), nrow = r, ncol = r)[r:1, r:1], breaks = breaks, col = cols, main = "SNG")
image(matrix(rowMeans(apply(spb.vs$Sigmas, 1, function(x) {cov2cor(x)})), nrow = r, ncol = r)[r:1, r:1], breaks = breaks, col = cols, main = "SPB")

maryclare/mvshrink documentation built on May 25, 2019, 9:34 p.m.