Description Usage Arguments Details Value Author(s) References See Also Examples
Simulates parameters and missing values from a joint posterior distribution under a normal model using Markov chain Monte Carlo.
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 | mcmcNorm(obj, ...)
## Default S3 method:
mcmcNorm(obj, x = NULL, intercept = TRUE,
starting.values, iter = 1000, multicycle = NULL,
seeds = NULL, prior = "uniform",
prior.df = NULL, prior.sscp = NULL, save.all.series = TRUE,
save.worst.series = FALSE, worst.linear.coef = NULL,
impute.every = NULL, ...)
## S3 method for class 'formula'
mcmcNorm(formula, data, starting.values,
iter = 1000, multicycle = NULL, seeds = NULL, prior = "uniform",
prior.df = NULL, prior.sscp = NULL, save.all.series = TRUE,
save.worst.series = FALSE, worst.linear.coef = NULL,
impute.every=NULL, ...)
## S3 method for class 'norm'
mcmcNorm(obj, starting.values = obj$param,
iter = 1000, multicycle = obj$multicycle,
seeds = NULL, prior = obj$prior, prior.df = obj$prior.df,
prior.sscp = obj$prior.sscp,
save.all.series = !(obj$method=="MCMC" & is.null( obj$series.beta )),
save.worst.series = !is.null( obj$worst.linear.coef ),
worst.linear.coef = obj$worst.linear.coef,
impute.every = obj$impute.every, ...)
|
obj |
an object used to select a method. It may be |
x |
a numeric matrix, vector or data frame of covariates to be
used as predictors for |
intercept |
if |
formula |
an object of class |
data |
an optional data frame, list or environment (or object
coercible by |
starting.values |
starting values for the model
parameters. This must be a list with two named components,
|
iter |
number of iterations to be performed. By default, each
iteration consists of one Imputation or I-step followed by
one Posterior or P-step, but this can be changed by
|
multicycle |
number of cycles per iteration, with
|
seeds |
two integers to initialize the random number generator; see DETAILS. |
prior |
should be |
prior.df |
prior degrees of freedom for a ridge
( |
prior.sscp |
prior sums of squares and cross-products (SSCP)
matrix for an inverted Wishart prior ( |
save.all.series |
if |
save.worst.series |
if |
worst.linear.coef |
vector or coefficients that define the worst
linear function of the parameters. Under ordinary circumstances,
these are provided automatically in the result from |
impute.every |
how many iterations to perform between
imputations? If |
... |
values to be passed to the methods. |
There are three different ways to specify the data and model when
calling mcmcNorm
:
by directly supplying as the initial argument a matrix of
numeric response variables y
, along with an optional
matrix of predictor variables x
;
by supplying a model specification
formula
, along with an optional data frame data
; or
by supplying an object of class
"norm"
, which was produced by an earlier call to
emNorm
or mcmcNorm
.
In the first case, the matrix y
is assumed to have a
multivariate normal
linear regression on x
with coefficients beta
and
covariance matrix sigma
, where
dim(beta)=c(ncol(x),ncol(y))
and
dim(sigma)=c(ncol(y),ncol(y))
. Missing values NA
are allowed in y
but not in x
.
In the second case, formula
is a formula for a (typically
multivariate) linear regression model in the manner expected by
lm
. A formula is given as y ~ model
, where
y
is either a single numeric variable or a matrix of numeric
variables bound together with the function cbind
. The
right-hand side of the formula (everything to the right of ~
) is a
linear predictor, a series of terms separated by operators +
,
:
or *
to specify main effects and
interactions. Factors are allowed on the right-hand side and will
enter the model as contrasts among the levels
. The
intercept term 1
is included by default; to remove the
intercept, use -1
.
In the third case, the initial argument to mcmcNorm
is an
object of class
"norm"
returned by a previous call to emNorm
or mcmcNorm
. The value of the parameters
carried in this object (the estimates from the last iteration of
EM or the simulated values from the last iteration of MCMC) will be
used as the starting values.
The matrix y
is assumed to have a multivariate normal
linear regression on x
with coefficients beta
and
covariance matrix sigma
, where
dim(beta)=c(ncol(x),ncol(y))
and
dim(sigma)=c(ncol(y),ncol(y))
.
Starting values for the parameters must be provided. In most cases
these will be the result of a previous call to emNorm
or
mcmcNorm
. If the starting
values are close to the mode (i.e., if they are the result of an EM
run that converged) then the worst linear function of the
parameters will be saved at each iteration. If the starting
values are the result of a previous run of MCMC, then the new
run will be a continuation of the same Markov chain.
If multicycle=
k for some k>1,
then the length of the saved parameter
series will be reduced by a factor of k, and the serial
correlation in the series will also be reduced. This option is
useful in large problems with many parameters and in slowly
converging problems for which many iterations are needed.
norm2
functions use their own internal random number generator which
is seeded by two integers, for example, seeds=c(123,456)
,
which allows results to be reproduced in the future. If
seeds=NULL
then
the function will seed itself with two random
integers from R. Therefore, results can also be made reproducible by
calling set.seed
beforehand and taking seeds=NULL
.
If prior="invwish"
then an inverted Wishart prior distribution
is applied to sigma
with hyperparameters prior.df
(a
scalar) and prior.sscp
(a symmetric, positive definite matrix
of the same dimension as sigma
). Using the device of imaginary
results, we can interpret prior.sscp/prior.df
as a prior guess
for sigma
, and prior.df
as the prior degrees of
freedom on which this guess is based.
The usual noninformative prior for the normal regression model
(prior="jeffreys"
) is equivalent to the inverted
Wishart density with prior.df
equal to 0 and
prior.sscp
equal to a matrix of 0's.
The ridge prior (prior="ridge"
) is a special case of the
inverted Wishart (Schafer, 1997). The prior
guess for sigma
is a diagonal matrix with diagonal elements
estimated by regressing the observed values in each column of
y
on the corresponding rows of x
. When
prior="ridge"
, the user must supply a value for
prior.df
, which
determines how strongly the estimated correlations are smoothed
toward zero.
If the first argument to mcmcNorm
is an object of class
"norm"
, then the parameter values stored in that object will
automatically be used as starting values.
For details of the MCMC algorithm, see the manual distributed
with the NORM package in the subdirectory doc
.
a list whose
class
attribute has been set to "norm"
.
This object may be
passed as the first argument in subsequent calls to emNorm
,
mcmcNorm
, impNorm
,
loglikNorm
or logpostNorm
. The
object also carries the original data and specifies the prior
distribution, so that these do not need to be provided again.
To see a summary of
this object, use the generic function summary
,
which passes the object to summary.norm
.
Components of the list may also be directly accessed and examined by the user. Components which may be of interest include:
iter |
number of MCMC iterations performed. |
param |
a list with elements |
loglik |
a numeric vector of length |
logpost |
a numeric vector of length |
series.worst |
a time-series object (class |
series.beta |
a multivariate time-series object (class
|
series.sigma |
a multivariate time-series object (class
|
imp.list |
a list containing the multiple imputations. Each
component of this list is a data matrix resembling |
miss.patt |
logical matrix with |
miss.patt.freq |
integer vector of length
|
which.patt |
integer vector of length |
Joe Schafer Joseph.L.Schafer@census.gov
Schafer, J.L. (1997) Analysis of Incomplete Multivariate
Data. London: Chapman & Hall/CRC Press.
For more information about this function and other functions in
the NORM package, see User's Guide for norm2
in the library subdirectory doc
.
emNorm
, summary.norm
,
impNorm
,
loglikNorm
,
logpostNorm
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | ## run EM for marijuana data with ridge prior
data(marijuana)
emResult <- emNorm(marijuana, prior="ridge", prior.df=0.5)
## run MCMC for 5,000 iterations starting from the
## posterior mode using the same prior
mcmcResult <- mcmcNorm(emResult, iter=5000)
## summarize and plot worst linear function
summary(mcmcResult)
plot(mcmcResult$series.worst)
acf(mcmcResult$series.worst, lag.max=50)
## generate 25 multiple imputations, taking
## 100 steps between imputations, and look st
## the first imputed dataset
mcmcResult <- mcmcNorm(emResult, iter=2500, impute.every=100)
mcmcResult$imp.list[[1]]
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.