Description Usage Arguments Details Value Note Author(s) References See Also Examples
Bayesian estimation via sampling from the posterior distribution of the of the mean and covariance matrix of multivariate normal (MVN) distributed data with a monotone missingness pattern, via Gibbs Sampling. Through the use of parsimonious/shrinkage regressions (lasso/NG & ridge), where standard regressions fail, this function can handle an (almost) arbitrary amount of missing data
1 2 3 4 5 
y 
data 
pre 
logical indicating whether preprocessing of the

p 
when performing regressions, 
B 
number of BurnIn MCMC sampling rounds, during which samples are discarded 
T 
total number of MCMC sampling rounds to take place after burnin, during which samples are saved 
thin 
multiplicative thinning in the MCMC. Each Bayesian
(lasso) regression will discard 
economy 
indicates whether memory should be economized at
the expense of speed. When 
method 
indicates the Bayesian parsimonious regression
specification to be used, choosing between the lasso (default)
of Park \& Casella, the NG extension, the horseshoe,
a ridge regression special case, and leastsquares.
The 
RJ 
indicates the Reversible Jump strategy to be employed.
The default argument of 
capm 
when 
start 
a list depicting starting values for the parameters
that are use to initialize the Markov chain. Usually this will be
a 
mprior 
prior on the number of nonzero regression coefficients
(and therefore covariates) 
rd 

theta 
the rate parameter ( 
rao.s2 
indicates whether to RaoBlackwellized samples for
s^2 should be used (default 
QP 
if non 
verb 
verbosity level; currently only 
trace 
if 
If pre = TRUE
then bmonomvn
first rearranges the columns
of y
into nondecreasing order with respect to the number of
missing (NA
) entries. Then (at least) the first column should
be completely observed.
Samples from the posterior distribution of the MVN mean vector and
covariance matrix are obtained sampling
from the posterior distribution of Bayesian regression models.
The methodology for converting these to samples from the mean vector
and covariance matrix is outlined in the monomvn
documentation, detailing a similarly structured maximum likelihood
approach. Also see the references below.
Whenever the regression model is ill–posed (i.e., when there are
more covariates than responses, or a
“big p
small n
” problem) then
Bayesian lasso or ridge regressions – possibly augmented with Reversible
Jump (RJ) for model selection – are used instead.
See the Park \& Casella reference below, and the blasso
documentation. To guarantee each regression is well posed the
combination setting of method="lsr"
and RJ="none"
is not allowed.
As in monomvn
the p
argument can be used to
turn on lasso or ridge regressions (possibly with RJ) at other times.
The exception is the "factor"
method which always involves
an OLS regression on (a subset of) the first p
columns of y
.
Samples from a function of samples of mu
and Sigma
can be obtained by specifying a Quadratic program via the
argument QP
. The idea is to allow for the calculation of
the distribution of minimum variance and mean–variance portfolios,
although the interface is quite general. See default.QP
for more details, as default.QP(ncol(y))
is used
when the argument QP = TRUE
is given. When a positive integer
is given, then the first QP
columns of y
are treated
as factors by using
default.QP(ncol(y)  QP)
instead. The result is that the corresponding components of (samples of)
mu
and rows/cols of S
are not factored into the
specification of the resulting Quadratic Program
bmonomvn
returns an object of class "monomvn"
,
which is a list
containing the inputs above and a
subset of the components below.
call 
a copy of the function call as used 
mu 
estimated mean vector with columns corresponding to the
columns of 
S 
estimated covariance matrix with rows and columns
corresponding to the columns of 
mu.var 
estimated variance of the mean vector with columns
corresponding to the columns of 
mu.cov 
estimated covariance matrix of the mean vector
with columns corresponding to the columns of 
S.var 
estimated variance of the individual components of the
covariance matrix with columns and rows corresponding to the columns
of 
mu.map 
estimated maximum a' posteriori (MAP) of the
mean vector with columns corresponding to the columns of 
S.map 
estimated MAP of the individual
components of the covariance matrix with columns and rows
corresponding to the columns of 
S.nz 
posterior probability that the individual entries of the covariance matrix are non–zero 
Si.nz 
posterior probability that the individual entries of the inverse of the covariance matrix are non–zero 
nu 
when 
lpost.map 
log posterior probability of the MAP estimate 
which.map 
gives the time index of the sample corresponding to the MAP estimate 
llik 
a trace of the log likelihood of the data 
llik.norm 
a trace of the log likelihood
under the Normal errors model when sampling under the
Studentt model; i.e., it is not present unless 
na 
when 
o 
when 
method 
method of regression used on each column, or

thin 
the (actual) number of thinning rounds used for the
regression ( 
lambda2 
records the mean lambda^2 value
found in the trace of the Bayesian Lasso regressions. Zerovalues
result when the column corresponds to a complete
case or an ordinary least squares regression (these would be

ncomp 
records the mean number of components
(columns of the design matrix) used in the regression model for
each column of 
trace 
if input 
R 
gives a 
B 
from inputs: number of BurnIn MCMC sampling rounds, during which samples are discarded 
T 
from inputs: total number of MCMC sampling rounds to take place after burnin, during which samples are saved 
r 
from inputs: alpha (shape) parameter to the gamma distribution prior for the lasso parameter lambda 
delta 
from inputs: beta (rate) parameter to the gamma distribution prior for the lasso parameter lambda 
QP 
if a valid (non– 
W 
when input 
Whenever the bmonomvn
algorithm requires a regression
where p >= n
, i.e., if any of the columns in the y
matrix have fewer non–NA
elements than the number of
columns with more non–NA
elements, then it is helpful
to employ both lasso/ridge and the RJ method.
It is important that any starting values provided in the
start
be compatible with the regression model
specified by inputs RJ
and method
. Any
incompatibilities will result with a warning that
(alternative) default action was taken and may result in
an undesired (possibly inferior) model being fit
Robert B. Gramacy rbg@vt.edu
R.B. Gramacy and E. Pantaleo (2010). Shrinkage regression for multivariate inference with missing data, and an application to portfolio balancing. Bayesian Analysis. 5(1), 237262. Preprint available on arXiv:0710.5837 http://arxiv.org/abs/0907.2135
Roderick J.A. Little and Donald B. Rubin (2002). Statistical Analysis with Missing Data, Second Edition. Wilely.
http://bobby.gramacy.com/r_packages/monomvn
blasso
, monomvn
,
default.QP
, em.norm
in the now defunct
norm
and mvnmle
packages, and returns
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83  ## standard usage, duplicating the results in
## Little and Rubin, section 7.4.3
data(cement.miss)
out < bmonomvn(cement.miss)
out
out$mu
out$S
##
## A bigger example, comparing the various
## parsimonious methods
##
## generate N=100 samples from a 10d random MVN
xmuS < randmvn(100, 20)
## randomly impose monotone missingness
xmiss < rmono(xmuS$x)
## using least squares only when necessary,
obl < bmonomvn(xmiss)
obl
## look at the posterior variability
par(mfrow=c(1,2))
plot(obl)
plot(obl, "S")
## compare to maximum likelihood
Ellik.norm(obl$mu, obl$S, xmuS$mu, xmuS$S)
oml < monomvn(xmiss, method="lasso")
Ellik.norm(oml$mu, oml$S, xmuS$mu, xmuS$S)
##
## a minvariance portfolio allocation example
##
## get the returns data, and use 20 random cols
data(returns)
train < returns[,sample(1:ncol(returns), 20)]
## missingness pattern requires DA; also gather
## samples from the solution to a QP
obl.da < bmonomvn(train, p=0, QP=TRUE)
## plot the QP weights distribution
plot(obl.da, "QP", xaxis="index")
## get ML solution: will warn about monotone violations
suppressWarnings(oml.da < monomvn(train, method="lasso"))
## add mean and MLE comparison, requires the
## quadprog library for the solve.QP function
add.pe.QP(obl.da, oml.da)
## now consider adding in the market as a factor
data(market)
mtrain < cbind(market, train)
## fit the model using only factor regressions
obl.daf < bmonomvn(mtrain, method="factor", p=1, QP=1)
plot(obl.daf, "QP", xaxis="index", main="using only factors")
suppressWarnings(oml.daf < monomvn(mtrain, method="factor"))
add.pe.QP(obl.daf, oml.daf)
##
## a Bayes/MLE comparison using least squares sparingly
##
## fit Bayesian and classical lasso
p < 0.25
obls < bmonomvn(xmiss, p=p)
Ellik.norm(obls$mu, obls$S, xmuS$mu, xmuS$S)
omls < monomvn(xmiss, p=p, method="lasso")
Ellik.norm(omls$mu, omls$S, xmuS$mu, xmuS$S)
## compare to ridge regression
obrs < bmonomvn(xmiss, p=p, method="ridge")
Ellik.norm(obrs$mu, obrs$S, xmuS$mu, xmuS$S)
omrs < monomvn(xmiss, p=p, method="ridge")
Ellik.norm(omrs$mu, omrs$S, xmuS$mu, xmuS$S)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.