mda.norm: Monotone data augmentation for incomplete multivariate normal data

Description

Monotone data augmentation under the usual noninformative prior, as described in Chapter 6 of Schafer (1996). This function simulates one or more iterations of a single Markov chain. One iteration consists of a random imputation of the missing data given the observed data and the current parameter value (I-step), followed by a draw from the posterior distribution of the parameter given the observed data and the imputed data (P-step). The I-step imputes only enough data to complete a monotone pattern, which typically makes this algorithm converge more quickly than da.norm, particularly when the observed data are nearly monotone. The order of the variables in the original data matrix determines the monotone pattern to be completed. For fast convergence, it helps to order the variables according to their rates of missingness, with the most observed (least missing) variable on the left and the least observed variable on the right.

Usage

mda.norm(s, theta, steps=1, showits=FALSE)

Arguments

summary list of an incomplete normal data matrix produced by the function prelim.norm.

theta

starting value of the parameter. This is a parameter vector in packed storage, such as one created by the function makeparam.norm. One obvious choice for a starting value is an ML estimate or posterior mode produced by em.norm.

steps

number of monotone data augmentation iterations to be simulated.

showits

if TRUE, reports the iterations so the user can monitor the progress of the algorithm.

Value

Returns a parameter vector, the result of the last P-step. If the value of steps was large enough to guarantee approximate stationarity, then this parameter can be regarded as a proper draw from the observed-data posterior, independent of start.

WARNING

Before this function may be used, the random number generator seed must be initialized with rngseed at least once in the current S session.

References

Chapter 6 of Schafer (1996).

Examples

Run this code

# NOT RUN {
data(mdata)
s <- prelim.norm(mdata)
thetahat <- em.norm(s)   #find the MLE for a starting value
rngseed(1234567)   #set random number generator seed
theta <- mda.norm(s,thetahat,steps=20,showits=TRUE)  # take 20 steps
getparam.norm(s,theta) # look at result
# }

Run the code above in your browser using DataCamp Workspace