get.IDR: Posterior class propbabilites and local, and adjusted IDR.

Description

Functions for computing posterior cluster probabilities (get.prob) in the general GMCM as well as local and adjusted irreproducibility discovery rates (get.IDR) in the special GMCM.

Usage

get.IDR(x, par, threshold = 0.05, ...)
get.idr(x, theta, ...)
get.prob(x, theta, ...)

Arguments

A matrix of observations where rows corresponds to features and colums to studies.

par

A vector of length 4 where par[1] is mixture proportion of the irreproducible component, par[2] is the mean value, par[3] is the standard deviation, and par[4] is the correlation of the reproducible comp

threshold

The thresholding level of the IDR rate.

theta

A list of parameters for the full model as described in rtheta.

...

Arguments passed to qgmm.marginal.

Value

get.IDR returns a list of length 5 with elements:
idrA vector of the local idr values. I.e. the posterior probability that x[i, ] belongs to the irreprodicible component.
IDRA vector of the adjusted IDR values.
lThe number of reproducible features at the specified threshold.
thresholdThe IDR threshold at which features are deemed reproducible.
KhatA vector signifying whether the corresponding feature is reproducible or not.
get.prob returns a matrix where entry (i,j) is the posterior probability that the observation x[i, ] belongs to cluster j.

References

Li, Q., Brown, J. B. J. B., Huang, H., & Bickel, P. J. (2011). Measuring reproducibility of high-throughput experiments. The Annals of Applied Statistics, 5(3), 1752-1779. doi:10.1214/11-AOAS466

Tewari, A., Giering, M., & Raghunathan, A. (2011). Parametric Characterization of Multimodal Distributions with Non-gaussian Modes. IEEE 11th International Conference on Data Mining Workshops, 2011, 286-292. doi:10.1109/ICDMW.2011.135

Examples

Run this code

set.seed(1123)

# True parameters
true.par <- c(0.9, 2, 0.7, 0.6)

# Simulation of data from the GMCM model
data <-  SimulateGMCMData(n = 1000, par = true.par, d = 2)

# Initial parameters
init.par <- c(0.5, 1, 0.5, 0.9)

# Nelder-Mead optimization
nm.par   <- fit.meta.GMCM(data$u, init.par = init.par, method = "NM")

# Get IDR values
res <- get.IDR(data$u, nm.par, threshold = 0.05)

# Plot results
plot(data$u, col = res$Khat, pch = c(3,16)[data$K])

Run the code above in your browser using DataLab