Nothing
## ----knit-int, echo=FALSE, include=FALSE---------------------------------
library("knitr")
set.seed(5828993)
opts_knit$set(self.contained = TRUE)
## ----load-packages, include=TRUE-----------------------------------------
#install.packages("GMCM") # Uncomment to install the GMCM package
library("GMCM")
## ----load-data, include=TRUE, echo=TRUE----------------------------------
x <- get(data("u133VsExon"))
head(x, n = 4)
## ------------------------------------------------------------------------
x <- x[1:5000, ]
## ----preprocess-data-----------------------------------------------------
u <- Uhat(1 - x)
## ------------------------------------------------------------------------
par(mfrow = c(1,2)) # Visualizing P-values and the ranked P-values
plot(x, cex = 0.5, pch = 4, col = "tomato", main = "P-values",
xlab = "P-value (U133)", ylab = "P-value (Exon)")
plot(u, cex = 0.5, pch = 4, col = "tomato", main = "Ranked P-values",
xlab = "rank(1-P) (U133)", ylab = "rank(1-P) (Exon)")
## ----show-initial-params, include=TRUE, echo=TRUE------------------------
init_par <- c(pie1 = 0.6, mu = 1, sigma = 1, rho = 0.2)
## ----fit_model, error=TRUE-----------------------------------------------
par <- fit.meta.GMCM(u = u,
init.par = init_par,
method = "NM",
max.ite = 1000,
verbose = FALSE)
print(par)
## ----compute_probs-------------------------------------------------------
meta_IDR_thres <- 0.05
out <- get.IDR(x, par = par,
threshold = meta_IDR_thres) # Compute IDR
str(out)
out <- get.IDR(u, par = par, threshold = meta_IDR_thres)
below <- out[["IDR"]] < meta_IDR_thres
out$l <- sum(below)
out$Khat <- ifelse(below, 2, 1)
## ----classes_table-------------------------------------------------------
table(out$Khat)
## ----plot_results--------------------------------------------------------
plot(x, col = out$Khat, asp = 1) # Plot of raw values
plot(u, col = out$Khat, asp = 1) # Plot of copula values
z <- GMCM:::qgmm.marginal(u, theta = meta2full(par, d = ncol(u))) # Estimate latent process
plot(z, col = out$Khat, asp = 1) # Plot of estimated latent process
## ----plot_theta----------------------------------------------------------
plot(meta2full(par, d = ncol(u)))
## ----session-info--------------------------------------------------------
sessionInfo()
## ----citation, echo=FALSE, results='asis'--------------------------------
cites <- lapply(c("GMCM", "knitr", "rmarkdown"), citation)
fmt_cites <- unlist(lapply(cites, format, style = "text"))
cat(paste0(" ", seq_along(fmt_cites), ". ", fmt_cites, "\n"))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.