mda.norm.RdMonotone 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.
mda.norm(s, theta, steps=1, showits=FALSE)summary list of an incomplete normal data matrix produced by the
function prelim.norm.
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.
number of monotone data augmentation iterations to be simulated.
if TRUE, reports the iterations so the user can monitor the progress
of the algorithm.
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.
Before this function may be used, the random number generator seed
must be initialized with rngseed at least once in the current S
session.
Chapter 6 of Schafer (1996).
rngseed, em.norm, prelim.norm, and getparam.norm.
data(mdata)
s <- prelim.norm(mdata)
thetahat <- em.norm(s) #find the MLE for a starting value
#> Iterations of EM:
#> 1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...
rngseed(1234567) #set random number generator seed
theta <- mda.norm(s,thetahat,steps=20,showits=TRUE) # take 20 steps
#> Steps of Monotone Data Augmentation:
#> 1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...18...19...20...
getparam.norm(s,theta) # look at result
#> $mu
#> [1] 42.241085 36.478323 13.558893 54054.686293 2.383784
#>
#> $sigma
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1.065414e+02 56.499139 12.558988 135664.856 3.634417
#> [2,] 5.649914e+01 37.038244 8.996055 90189.082 1.494470
#> [3,] 1.255899e+01 8.996055 26.496029 61020.445 -2.598673
#> [4,] 1.356649e+05 90189.081775 61020.444688 653204090.852 -4411.142222
#> [5,] 3.634417e+00 1.494470 -2.598673 -4411.142 1.350514
#>