NMixMCMC {mixAK} | R Documentation |
This function runs MCMC for a model in which unknown density is specified as a normal mixture with either known or unknown number of components. With a prespecified number of components, MCMC is implemented through Gibbs sampling (see Diebolt and Robert, 1994) and dimension of the data can be arbitrary. With unknown number of components, currently only univariate case is implemented using the reversible jump MCMC (Richardson and Green, 1997).
Further, the data are allowed to be censored in which case additional Gibbs step is used within the MCMC algorithm
NMixMCMC(y0, y1, censor, scale, prior, init, init2, RJMCMC, nMCMC=c(burn=10, keep=10, thin=1, info=10), PED, keep.chains=TRUE, onlyInit=FALSE, dens.zero=1e-300) ## S3 method for class 'NMixMCMC': print(x, dic, ...) ## S3 method for class 'NMixMCMClist': print(x, ped, dic, ...)
y0 |
numeric vector of length n or n x p matrix with observed data. It contains exactly observed, right-censored, left-censored data and lower limits for interval-censored data. |
y1 |
numeric vector of length n or n x p
matrix with upper limits for interval-censored data. Elements
corresponding to exactly observed, right-censored or left-censored
data are ignored and can be filled arbitrarily (by
NA 's) as well.
It does not have to be supplied if there are no interval-censored data. |
censor |
numeric vector of length n or n x p
matrix with censoring indicators. The following values indicate:
If it is not supplied then it is assumed that all values are exactly observed. |
scale |
a list specifying how to scale the data before running
MCMC. It should have two components:
If there is no censoring, and argument scale is missing
then the data are scaled to have zero mean and unit variances, i.e.,
scale(y0) is used for MCMC. In the case there is censoring
and scale is missing,
scale$shift is taken to be a sample mean of init$y and
scale$scale are sample standard deviations of columns of init$y .
If you do not wish to scale the data before running MCMC, specify scale=list(shift=0, scale=1) .
|
prior |
a list with the parameters of the prior distribution. It
should have the following components (for some of them,
the program can assign default values and the user does not have to
specify them if he/she wishes to use the defaults):
|
init |
a list with the initial values for the MCMC. All initials
can be determined by the program if they are not specified. The list
may have the following components:
|
init2 |
a list with the initial values for the second chain
needed to estimate the penalized expected deviance of Plummer
(2008). The list init2 has the same structure as the list
init . All initials in init2 can be determined by the
program (differently than the values in init ) if they are not
specified.
Ignored when PED is FALSE .
|
RJMCMC |
a list with the parameters needed to run reversible jump
MCMC for mixtures with varying number of components. It does not
have to be specified if the number of components is fixed. Most of
the parameters can be determined by the program if they are not
specified. The list may have the following components:
|
nMCMC |
numeric vector of length 4 giving parameters of the MCMC
simulation. Its components may be named (ordering is then unimportant) as:
|
PED |
a logical value which indicates whether the penalized
expected deviance (see Plummer, 2008 for more details)
is to be computed (which requires two parallel
chains). If not specified, PED is set to TRUE
for models with fixed number of components and is set to
FALSE for models with numbers of components estimated using RJ-MCMC.
|
keep.chains |
logical. If FALSE , only summary statistics
are returned in the resulting object. This might be useful in the
model searching step to save some memory. |
onlyInit |
logical. If TRUE then the function only
determines parameters of the prior distribution, initial values,
values of scale and
parameters for the reversible jump MCMC. |
dens.zero |
a small value used instead of zero when computing deviance related quantities. |
x |
an object of class NMixMCMC or NMixMCMClist to
be printed. |
dic |
logical which indicates whether DIC should be printed. By default, DIC is printed only for models with a fixed number of mixture components. |
ped |
logical which indicates whether PED should be printed. By default, PED is printed only for models with a fixed number of mixture components. |
... |
additional arguments passed to the default print method. |
See accompanying paper (Komárek, 2009). In the rest of the helpfile, the same notation is used as in the paper, namely, n denotes the number of observations, p is dimension of the data, K is the number of mixture components, w[1],...,w[K] are mixture weights, mu[1],...,mu[K] are mixture means, Sigma_1,...,Sigma_K are mixture variance-covariance matrices, Q_1,...,Q_K are their inverses.
For the data y[1],...,y[n] the following g(y) density is assumed
g(y) = |S|^(-1) * sum[j=1]^K w[j] phi(S^(-1)*(y - m) | mu[j], Sigma[j]),
where phi(. | mu, Sigma) denotes a density of the (multivariate) normal distribution with mean mu and a~variance-covariance matrix Σ. Finally, S is a pre-specified diagonal scale matrix and m is a pre-specified shift vector. Sometimes, by setting m to sample means of components of y and diagonal of S to sample standard deviations of y (considerable) improvement of the MCMC algorithm is achieved.
An object of class NMixMCMC
or class NMixMCMClist
.
Object of class NMixMCMC
is returned if PED
is
FALSE
. Object of class NMixMCMClist
is returned if
PED
is TRUE
.
Objects of class NMixMCMC
have the following components:
iter |
index of the last iteration performed. |
nMCMC |
used value of the argument nMCMC . |
dim |
dimension p of the distribution of data |
prior |
a list containing the used value of the argument prior . |
init |
a list containing the used value of the argument init . |
RJMCMC |
a list containing the used value of the argument RJMCMC . |
scale |
a list containing the used value of the argument scale . |
state |
a list having the components labeled
y , K , w , mu , Li , Q , Sigma ,
gammaInv , r containing the last sampled values of
generic parameters. |
freqK |
frequency table of K based on the sampled chain. |
propK |
posterior distribution of K based on the sampled chain. |
DIC |
a data.frame having columns labeled
DIC , pD , D.bar , D.in.bar containing
values used to compute deviance information criterion
(DIC). Currently only DIC[3] of Celeux et al. (2006) is
implemented.
|
moves |
a data.frame which summarizes the acceptance
probabilities of different move types of the sampler. |
K |
numeric vector with a chain for K (number of mixture components). |
w |
numeric vector or matrix with a chain for w (mixture weights). It is a matrix with K columns when K is fixed. Otherwise, it is a vector with weights put sequentially after each other. |
mu |
numeric vector or matrix with a chain for mu (mixture means). It is a matrix with p*K columns when K is fixed. Otherwise, it is a vector with means put sequentially after each other. |
Q |
numeric vector or matrix with a chain for lower triangles of Q (mixture inverse variances). It is a matrix with (p*(p+1)2)*K columns when K is fixed. Otherwise, it is a vector with lower triangles of Q matrices put sequentially after each other. |
Sigma |
numeric vector or matrix with a chain for lower triangles of Sigma (mixture variances). It is a matrix with (p*(p+1)2)*K columns when K is fixed. Otherwise, it is a vector with lower triangles of Sigma matrices put sequentially after each other. |
Li |
numeric vector or matrix with a chain for lower triangles of Cholesky decompositions of Q matrices. It is a matrix with (p*(p+1)2)*K columns when K is fixed. Otherwise, it is a vector with lower triangles put sequentially after each other. |
gammaInv |
matrix with p columns with a chain for inverses of the hyperparameter gamma. |
order |
numeric vector or matrix with order indeces of mixture components. It is a matrix with K columns when K is fixed. Otherwise it is a vector with orders put sequentially after each other. |
rank |
numeric vector or matrix with rank indeces of mixture components. It is a matrix with K columns when K is fixed. Otherwise it is a vector with ranks put sequentially after each other. |
mixture |
data.frame with columns labeled
y.Mean.* , y.SD.* , y.Corr.*.* ,
z.Mean.* , z.SD.* , z.Corr.*.* containing the
chains for the means, standard deviations and correlations of the
distribution of the original (y ) and scaled (z ) data
based on a normal mixture at each iteration.
|
deviance |
data.frame with columns labeles
LogL0 , LogL1 , dev.complete , dev.observed
containing the chains of quantities needed to compute DIC.
|
pm.y |
a data.frame with p columns with posterior
means for (latent) values of observed data (useful when there is
censoring). |
pm.z |
a data.frame with p columns with posterior
means for (latent) values of scaled observed data (useful when there is censoring). |
pm.indDev |
a data.frame with columns labeled
LogL0 , LogL1 , dev.complete ,
dev.observed , pred.dens containing posterior means of
individual contributions to the deviance.
|
pred.dens |
a numeric vector with the predictive density of the
data based on the MCMC sample evaluated at data points.
Note that when there is censoring, this is not exactly the predictive density as it is computed as the average of densities at each iteration evaluated at sampled values of latent observations at iterations. |
poster.comp.prob1 |
a matrix which is present in the output object
if the number of mixture components in the distribution of random
effects is fixed and equal to K. In that case,
poster.comp.prob1 is a matrix with K columns and n
rows with estimated posterior component probabilities
– posterior means of the components of the underlying 0/1
allocation vector.
These can be used for possible clustering of the subjects. |
poster.comp.prob2 |
a matrix which is present in the output object
if the number of mixture components in the distribution of random
effects is fixed and equal to K. In that case,
poster.comp.prob2 is a matrix with K columns and n
rows with estimated posterior component probabilities
– posterior mean over model parameters.
These can also be used for possible clustering of the subjects |
summ.y.Mean |
Posterior summary statistics based on chains stored
in y.Mean.* columns of the data.frame mixture . |
summ.y.SDCorr |
Posterior summary statistics based on chains
stored in y.SD.* and y.Corr.*.* columns of the
data.frame mixture . |
summ.z.Mean |
Posterior summary statistics based on chains stored
in z.Mean.* columns of the data.frame mixture . |
summ.z.SDCorr |
Posterior summary statistics based on chains
stored in z.SD.* and z.Corr.*.* columns of the data.frame mixture . |
poster.mean.w |
a numeric vector with posterior means of mixture weights after re-labeling. It is computed only if K is fixed and even then I am not convinced that these are useful posterior summary statistics. In any case, they should be used with care. |
poster.mean.mu |
a matrix with posterior means of mixture means after re-labeling. It is computed only if K is fixed and even then I am not convinced that these are useful posterior summary statistics. In any case, they should be used with care. |
poster.mean.Q |
a list with posterior means of mixture inverse variances after re-labeling. It is computed only if K is fixed and even then I am not convinced that these are useful posterior summary statistics. In any case, they should be used with care. |
poster.mean.Sigma |
a list with posterior means of mixture variances after re-labeling. It is computed only if K is fixed and even then I am not convinced that these are useful posterior summary statistics. In any case, they should be used with care. |
poster.mean.Li |
a list with posterior means of Cholesky decompositions of mixture inverse variances after re-labeling. It is computed only if K is fixed and even then I am not convinced that these are useful posterior summary statistics. In any case, they should be used with care. |
Objects of class NMixMCMClist
is the list having two components
of class NMixMCMC
representing two parallel chains and
additionally the following components:
D.expect
=
estimated expected deviance, where the estimate is based on two
parallel chains; popt
= estimated penalty, where the
estimate is based on simple MCMC average based on two parallel
chains; PED
= estimated penalized expected deviance
= D.expect
+ popt
; wpopt
=
estimated penalty, where the estimate is based on weighted MCMC average
(through importance sampling) based on two parallel chains;
wPED
= estimated penalized expected deviance =
D.expect
+ wpopt
.dens.zero
) and was not taken into account when
computing D.expect
.popt
.
wpopt
.Arnošt Komárek arnost.komarek[AT]mff.cuni.cz
Celeux, G., Forbes, F., Robert, C. P., and Titterington, D. M. (2006). Deviance information criteria for missing data models. Bayesian Analysis, 1, 651-674.
Cappé, Robert and Rydén (2003). Reversible jump, birth-and-death and more general continuous time Markov chain Monte Carlo samplers. Journal of the Royal Statistical Society, Series B, 65, 679-700.
Diebolt, J. and Robert, C. P. (1994). Estimation of finite mixture distributions through Bayesian sampling. Journal of the Royal Statistical Society, Series B, 56, 363–375.
Komárek, A. (2009). A new R package for Bayesian estimation of multivariate normal mixtures allowing for selection of the number of components and interval-censored data. Computational Statistics and Data Analysis, 53, 3932–3947.
Plummer, M. (2008). Penalized loss functions for Bayesian model comparison. Biostatistics, 9, 523-539.
Richardson, S. and Green, P. J. (1997). On Bayesian analysis of mixtures with unknown number of components (with Discussion). Journal of the Royal Statistical Society, Series B, 59, 731-792.
Spiegelhalter, D. J.,Best, N. G., Carlin, B. P., and van der Linde, A. (2002). Bayesian measures of model complexity and fit (with Discussion). Journal of the Royal Statistical Society, Series B, 64, 583-639.
NMixPredDensMarg
, NMixPredDensJoint2
.
## See also additional material available in ## YOUR_R_DIR/library/mixAK/doc/ ## or YOUR_R_DIR/site-library/mixAK/doc/ ## - files Galaxy.pdf, Faithful.pdf, Tandmob.pdf ## ============================================== ## Simple analysis of Anderson's iris data ## ============================================== library("colorspace") data(iris, package="datasets") summary(iris) VARS <- names(iris)[1:4] #COLS <- rainbow_hcl(3, start = 60, end = 240) COLS <- c("red", "darkblue", "darkgreen") names(COLS) <- levels(iris[, "Species"]) ### Prior distribution and the length of MCMC Prior <- list(priorK = "fixed", Kmax = 3) nMCMC <- c(burn=5000, keep=10000, thin=5, info=1000) ### Run MCMC set.seed(20091230) fit <- NMixMCMC(y0 = iris[, VARS], prior = Prior, nMCMC = nMCMC) ### Basic posterior summary print(fit) ### Univariate marginal posterior predictive densities ### based on chain #1 pdens1 <- NMixPredDensMarg(fit[[1]], lgrid=150) plot(pdens1) plot(pdens1, main=VARS, xlab=VARS) ### Bivariate (for each pair of margins) predictive densities ### based on chain #1 pdens2a <- NMixPredDensJoint2(fit[[1]]) plot(pdens2a) plot(pdens2a, xylab=VARS) plot(pdens2a, xylab=VARS, contour=TRUE) ### Determine the grid to compute bivariate densities grid <- list(Sepal.Length=seq(3.5, 8.5, length=75), Sepal.Width=seq(1.8, 4.5, length=75), Petal.Length=seq(0, 7, length=75), Petal.Width=seq(-0.2, 3, length=75)) pdens2b <- NMixPredDensJoint2(fit[[1]], grid=grid) plot(pdens2b, xylab=VARS) ### Plot with contours ICOL <- rev(heat_hcl(20, c=c(80, 30), l=c(30, 90), power=c(1/5, 2))) oldPar <- par(mfrow=c(2, 3), bty="n") for (i in 1:3){ for (j in (i+1):4){ NAME <- paste(i, "-", j, sep="") MAIN <- paste(VARS[i], "x", VARS[j]) image(pdens2b$x[[i]], pdens2b$x[[j]], pdens2b$dens[[NAME]], col=ICOL, xlab=VARS[i], ylab=VARS[j], main=MAIN) contour(pdens2b$x[[i]], pdens2b$x[[j]], pdens2b$dens[[NAME]], add=TRUE, col="brown4") } } ### Plot with data for (i in 1:3){ for (j in (i+1):4){ NAME <- paste(i, "-", j, sep="") MAIN <- paste(VARS[i], "x", VARS[j]) image(pdens2b$x[[i]], pdens2b$x[[j]], pdens2b$dens[[NAME]], col=ICOL, xlab=VARS[i], ylab=VARS[j], main=MAIN) for (spec in levels(iris[, "Species"])){ Data <- subset(iris, Species==spec) points(Data[,i], Data[,j], pch=16, col=COLS[spec]) } } } ### Set the graphical parameters back to their original values par(oldPar) ### Clustering based on posterior summary statistics of component allocations ### or on the posterior distribution of component allocations ### (these are two equivalent estimators of probabilities of belonging ### to each mixture components for each observation) p1 <- fit[[1]]$poster.comp.prob1 p2 <- fit[[1]]$poster.comp.prob2 ### Clustering based on posterior summary statistics of mixture weight, means, variances p3 <- NMixPlugDA(fit[[1]], iris[, VARS]) p3 <- p3[, paste("prob", 1:3, sep="")] ### Observations from "setosa" species (all would be allocated in component 1) apply(p1[1:50,], 2, quantile, prob=seq(0, 1, by=0.1)) apply(p2[1:50,], 2, quantile, prob=seq(0, 1, by=0.1)) apply(p3[1:50,], 2, quantile, prob=seq(0, 1, by=0.1)) ### Observations from "versicolor" species (almost all would be allocated in component 2) apply(p1[51:100,], 2, quantile, prob=seq(0, 1, by=0.1)) apply(p2[51:100,], 2, quantile, prob=seq(0, 1, by=0.1)) apply(p3[51:100,], 2, quantile, prob=seq(0, 1, by=0.1)) ### Observations from "virginica" species (all would be allocated in component 3) apply(p1[101:150,], 2, quantile, prob=seq(0, 1, by=0.1)) apply(p2[101:150,], 2, quantile, prob=seq(0, 1, by=0.1)) apply(p3[101:150,], 2, quantile, prob=seq(0, 1, by=0.1))