bmonomvn {monomvn} | R Documentation |
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 (currently only lasso is supported), where standard regressions fail, this function can handle an (almost) arbitrary amount of missing data.
bmonomvn(y, pre = TRUE, p = 0.9, B = 100, T = 200, thin = 10, method = c("default", "rjlasso", "rjlsr", "lasso"), capm = method!="lasso", start = NULL, r = 2, delta = 0.1, rao.s2 = TRUE, verb = 1, trace = FALSE)
y |
data matrix were each row is interpreted as a
random sample from a MVN distribution with missing
values indicated by NA |
pre |
logical indicating whether pre-processing of the
y is to be performed. This sorts the columns so that the
number of NA s is non-decreasing with the column index |
p |
when performing regressions, p is the proportion of the
number of columns to rows in the design matrix before an
alternative regression (lasso) is performed as if
least-squares regression “failed”. Least-squares regression is
known to fail when the number of columns equals the number of rows,
hence a default of p = 0.9 close to 1 .
Alternatively, setting p = 0
forces lasso to be used for every regression.
Intermediate settings of p allow the user to control when
least-squares regressions stop and the lasso ones start |
B |
number of Burn-In MCMC sampling rounds, during which samples are discarded |
T |
total number of MCMC sampling rounds to take place after burn-in, during which samples are saved |
thin |
level of thinning in the MCMC, i.e., the number of MCMC rounds must be collected before a sample is saved |
method |
indicates the lasso regression specification
to be used. The "default" uses the Trevor & Park formulation
until p >= n at which point a Reversible Jump (RJ)
is turned on in order to keep the number of non-zero regression
coefficients beta at most n (when capm =
TRUE . The "rjlasso" method uses RJ at all times that
the lasso is used, whereas "rjlsr" uses RJ without the
lasso prior. Finally, the "lasso" method never uses RJ.
See blasso and note below for more details |
capm |
when TRUE this argument indicates that the
number of components of
beta should not exceed n, the number of
response variables in a particular regression. This argument
is ignored when using method = "lasso" |
start |
a list depicting starting values for the parameters
that are use to initaliza the Markov chain. Usually this will be
a "monomvn" -class object depsicting maximum likelihood
estimates output from the monomvn function, but
the relevant fields are the mean vector $mu , covariance
matrix $S , monotone ordering $o (for sanity checking
with intput y ), component vector $ncomp and
penalty parameter vector $lambda . See note below |
r |
alpha (shape) parameter to the gamma distribution prior for the lasso parameter lambda |
delta |
beta (rate) parameter to the gamma distribution prior for the lasso parameter lambda |
rao.s2 |
indicates whether to use Rao-Blackwellized samples for
s^2 should be used (default TRUE ), see
the details section of blasso for more information |
verb |
verbosity level; currently only verb = 0 and
verb = 1 are supported |
trace |
if TRUE then samples from the regressions are
saved to files in the CWD |
If pre = TRUE
then bmonomvn
first re-arranges 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 over a Multivariate
Normal mean vector and covariance matrix are obtained by obtaining
samples from the posterior distribution of Bayesian regression models.
The method by which these samples from the regression posterior(s) are
used to obtain samples from the mean vector and covariance matrix
is outlined in the monomvn
documentation, detailing a
similarly structured maximum likelihood approach. See also 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 regressions – possibly augmented with Reversible
Jump (RJ) for model selection – are used instead.
See the Park & Casella reference below, and the blasso
documentation. The p
argument
can be used to turn on Lasso regressions at other times.
One difference between the Bayesian and MLE approach is that the MLE approach treats the complete (fully observed) columns jointly, without performing regressions. In contrast, the Bayesian only treats the first complete column without regressions. The remaining complete columns are processed via regression.
bmonomvn
returns an object of class "monomvn"
, which is a
list containing 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 y |
S |
estimated covariance matrix with rows and columns
corresponding to the columns of y |
mu.var |
estimated variance of mean vector with columns
corresponding to the columns of y |
S.var |
estimated variance of the individual components of the
covariance matrix with columns corresponding to the columns
of y |
na |
when pre = TRUE this is a vector containing number of
NA entries in each column of y |
o |
when pre = TRUE this is a vector containing the
index of each column in the sorting of the columns of y
obtained by o <- order(na) |
method |
method of regression used on each column, or
"bcomplete" indicating that no regression was used |
lambda2 |
records the mean lambda^2 value
found in the trace of the Bayesian Lasso regressions. This value
will be zero if the corresponding column corresponds to a complete
case or a ordinary least squares regression (these would be
NA entries maximum likelihood estimates are sought
via monomvn |
ncomp |
records the mean number of components
(columns of the design matrix) used in the regression model for
each column of y . If input RJ = FALSE then this simply
corresponds to the monotone ordering (these would correspond to
the NA entries when maximum likelihood estimates are sought
via monomvn . When RJ = TRUE
the monotone ordering is an upper bound (on each entry) |
trace |
if input trace = TRUE then this field contains
traces of the samples of mu in the field $mu and
of S in the field $S , and of all regression parameters
for each of the m = length(mu) columns in the field
$reg |
B |
from inputs: number of Burn-In MCMC sampling rounds, during which samples are discarded |
T |
from inputs: total number of MCMC sampling rounds to take place after burn-in, 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 |
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 and RJ method. Therefore, should this
case be detected and one of the "rjlsr"
or "lasso"
methods is specified, a warning will be printed
It is important than any starting values provided in the
start
input list be compatible with the regression model
specified by inputs RJ
and method
. Any
incompatibilites 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 bobby@statslab.cam.ac.uk
Robert B. Gramacy and Joo Hee Lee (2007).
On estimating covariances between many assets with histories
of highly variable length. Preprint available on arXiv:0710.5837:
http://arxiv.org/abs/0710.5837
Roderick J.A. Little and Donald B. Rubin (2002). Statistical Analysis with Missing Data, Second Edition. Wilely.
Park, T., Casella, G. (2008).
The Bayesian Lasso, (unpublished)
http://www.stat.ufl.edu/~casella/Papers/bayeslasso.pdf
http://www.statslab.cam.ac.uk/~bobby/monomvn.html
blasso
, monomvn
,
em.norm
in the norm package,
and mlest
in the mvnmle package
## 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 methods ## ## generate N=1000 samples from a random MVN xmuS <- randmvn(1000, 100) ## randomly impose monotone missingness xmiss <- rmono(xmuS$x) ## using least squares only when necessary out.b <- bmonomvn(xmiss) out.b kl.norm(out.b$mu, out.b$S, xmuS$mu, xmuS$S, symm=TRUE) out.mle <- monomvn(xmiss, method="lasso") kl.norm(out.mle$mu, out.mle$S, xmuS$mu, xmuS$S, symm=TRUE) ## using least squares sparingly out.b.s <- bmonomvn(xmiss, p=0.25) kl.norm(out.b.s$mu, out.b.s$S, xmuS$mu, xmuS$S, symm=TRUE) out.mle.s <- monomvn(xmiss, p=0.25, method="lasso") kl.norm(out.mle.s$mu, out.mle.s$S, xmuS$mu, xmuS$S, symm=TRUE) ## using the maximum likelihood solution to initialize ## the Markov chain and avoid burn-in. out.b2.s <- bmonomvn(xmiss, p=0.25, B=0, start=out.mle.s, method="rjlasso") kl.norm(out.b2.s$mu, out.b2.s$S, xmuS$mu, xmuS$S, symm=TRUE)