In stephenslab/mvebnm: Ultimate Deconvolution for Multivariate Normal Means

knitr::opts_chunk$set(comment = "#",collapse = TRUE,results = "hold",
                      fig.align = "center")

Here we illustrate the features of the Ultimate Deconvolution (UD) package.

library(udr)

We simulate 4,000 data points from a UD model with 4 mixture components: $$ x_i \sim w_1 N(0,V + U_1) + \cdots + w_4 N(0,V + U_4). $$

set.seed(1)
n <- 4000
V <- rbind(c(0.8,0.2),
           c(0.2,1.5))
U <- list(none   = rbind(c(0,0),
                         c(0,0)),
          shared = rbind(c(1.0,0.9),
                         c(0.9,1.0)),
          only1  = rbind(c(1,0),
                         c(0,0)),
          only2  = rbind(c(0,0),
                         c(0,1)))
w <- c(0.8,0.1,0.075,0.025)
X <- simulate_ud_data(n,w,U,V)

This is the most basic usage, using the default settings:

Note: in this toy simulation example, we set model parameter $V$ to $V$ as we know the simulation truth. In practice, it's likely one needs to estimate $V$ and use $\hat V$ for $V$.

fit0 <- ud_init(X, V = V)
fit1 <- ud_fit(fit0)

The default model is a mixture with $K = 10$ mixture components: 2 scaled scaled and 8 unconstrained prior covariance matrices $U_k$:

summary(fit1)

The function ud_init is a flexible interface for defining a variety of UD models. For example,

fit0 <- ud_init(X,U_scaled = U,n_rank1 = 1,n_unconstrained = 1,V = V)
fit2 <- ud_fit(fit0)
summary(fit2)

The ud_fit interface also allows for flexibly re-fitting a model; that is, a ud_fit output can be used in another call to ud_fit. There is also a more advanced model fitting interface for more fine-grained control; see help(ud_fit_advanced) for details.

Notice that the data do not need to be passed to ud_fit because they are already embedded in the ud_fit output. However, one can also pass a new data set X to ud_fit.

The control argument to ud_fit controls the optimization of the model parameters. For example, the unconstrained prior covariance matrices are by default optimized using the "truncated eigenvalue decomposition" (TED) updates, and we can replace these with the Extreme Deconvolution (ED) updates:

fit3 <- ud_fit(fit0,control = list(unconstrained.update = "ed"))

Notice that the likelihood is slightly higher with the TED updates, reflecting the fact that the TED updates have greater freedom to optimize the covariance matrices.

These are the current model fitting defaults:

unlist(ud_fit_control_default())

The most intensive computations are implemented in R (control$version = R) and C++ (control$version = Rcpp); both versions should return the same, or nearly the same output, but the C++ version can sometimes be much faster.

In some cases each sample might have different measurement error, and udr can handle such cases by specifying V as a list instead of a matrix:

fit0 <- ud_init(X,V = rep(list(V),n))
fit4 <- ud_fit(fit0)

The TED updates cannot be applied to this setting, so the ED updates must be used instead.

Session info

This is the version of R and the packages that were used to generate these results:

sessionInfo()

stephenslab/mvebnm documentation built on June 4, 2024, 6:37 a.m.