[Robast-commits] r1037 - in pkg/RobLoxBioC: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jul 23 22:31:11 CEST 2018
Author: ruckdeschel
Date: 2018-07-23 22:31:10 +0200 (Mon, 23 Jul 2018)
New Revision: 1037
Modified:
pkg/RobLoxBioC/DESCRIPTION
pkg/RobLoxBioC/NAMESPACE
pkg/RobLoxBioC/R/AffySimStudyFunction.R
pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R
pkg/RobLoxBioC/R/robloxbiocBeadLevelData.R
pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd
pkg/RobLoxBioC/man/robloxbioc.Rd
Log:
[RobLoxBioC] merged branch 1.1 into trunk
Modified: pkg/RobLoxBioC/DESCRIPTION
===================================================================
--- pkg/RobLoxBioC/DESCRIPTION 2018-07-23 20:29:47 UTC (rev 1036)
+++ pkg/RobLoxBioC/DESCRIPTION 2018-07-23 20:31:10 UTC (rev 1037)
@@ -1,19 +1,18 @@
-Package: RobLox
+Package: RobLoxBioC
Version: 1.1.0
Date: 2018-07-17
-Title: Optimally Robust Influence Curves and Estimators for Location and Scale
+Title: Infinitesimally Robust Estimators for Preprocessing -Omics Data
Description: Functions for the determination of optimally robust influence curves and
- estimators in case of normal location and/or scale.
-Depends: R(>= 2.14.0), stats, distrMod(>= 2.5.2), RobAStBase(>= 0.9)
-Imports: methods, lattice, RColorBrewer, Biobase, RandVar(>= 0.9.2), distr(>= 2.5.2)
-Suggests: MASS
-Authors at R: c(person("Matthias", "Kohl", role=c("cre", "cph"),
- email="Matthias.Kohl at stamats.de"), person("Peter", "Ruckdeschel", role=c("aut",
- "cph")))
+ estimators for preprocessing omics data, in particular gene expression data.
+Depends: R(>= 2.14.0), methods, distr(>= 2.5.2), affy
+Imports: Biobase, BiocGenerics, beadarray, RobLox(>= 0.9.2), distrMod(>= 2.5.2), lattice,
+ RColorBrewer, AnnotationDbi
+Suggests: affydata, hgu95av2cdf, beadarrayExampleData, illuminaHumanv3.db
+Authors at R: person("Matthias", "Kohl", role=c("aut", "cre", "cph"), email="Matthias.Kohl at stamats.de")
ByteCompile: yes
License: LGPL-3
-Encoding: latin1
URL: http://robast.r-forge.r-project.org/
+Encoding: latin1
LastChangedDate: {$LastChangedDate$}
LastChangedRevision: {$LastChangedRevision$}
VCS/SVNRevision: 940
Modified: pkg/RobLoxBioC/NAMESPACE
===================================================================
--- pkg/RobLoxBioC/NAMESPACE 2018-07-23 20:29:47 UTC (rev 1036)
+++ pkg/RobLoxBioC/NAMESPACE 2018-07-23 20:31:10 UTC (rev 1037)
@@ -1,48 +1,23 @@
-importFrom("Biobase", "rowMedians")
-importFrom("lattice","stripplot")
+import("methods", "distr")
+importMethodsFrom("distrMod","estimate")
importFrom("RColorBrewer","brewer.pal")
-importFrom("distr", "p", "q", "r", "q.l", "Reals")
-importClassesFrom("distr","Reals")
-importFrom("RandVar", "EuclRandVarList", "RealRandVariable")
-importClassesFrom("RandVar", "EuclRandVarList", "RealRandVariable")
-importFrom("distrMod", "NormLocationFamily", "NormLocationScaleFamily", "symmetricBias", "NormType")
-importClassesFrom("distrMod", "L2LocationFamily", "NormLocationFamily", "NormLocationScaleFamily")
-importMethodsFrom("distrMod", "distribution", "main", "addInfo<-", "Infos", "Infos<-",
- "normtype", "biastype","normtype<-","biastype<-", "estimate")
-importFrom("RobAStBase","generateIC","ContNeighborhood", "IC")
+importFrom("lattice","stripplot")
+importMethodsFrom("BiocGenerics","append","annotation")
+importFrom("Biobase","assayDataNew")
+importFrom("AnnotationDbi", "revmap")
+importClassesFrom("Biobase","ExpressionSet","AssayData","AnnotatedDataFrame")
+importMethodsFrom("Biobase", "rowMedians", "featureNames", "sampleNames",
+ "phenoData", "experimentData", "exprs", "exprs<-", "assayData","assayData<-",
+ "featureData<-", "phenoData<-")
+importFrom("affy","bg.correct.mas","getCdfInfo", "tukey.biweight")
+importMethodsFrom("affy","intensity")
+importMethodsFrom("AnnotationDbi", "mappedkeys")
+import("beadarray")
+importFrom("RobLox", "rowRoblox", "finiteSampleCorrection")
importFrom("grDevices", "dev.new")
-importFrom("graphics", "abline", "boxplot", "curve", "layout",
- "legend", "par", "plot")
-importFrom("methods", "is", "new")
-importFrom("stats", "approx", "complete.cases", "dnorm", "integrate",
- "mad", "median", "na.omit", "optim", "optimize", "pnorm",
- "qnorm", "rbinom", "rnorm", "uniroot")
-importClassesFrom("RobAStBase", "HampelWeight", "ALEstimate", "kStepEstimate", "ContNeighborhood")
-importMethodsFrom("RobAStBase", "clip", "cent", "stand", "weight", "clip<-", "cent<-", "stand<-", "weight<-",
- "getweight", "makeIC", "CallL2Fam", "CallL2Fam<-", "neighborRadius", "modifyIC", "addInfo<-")
-
-export(finiteSampleCorrection,
- rlsOptIC.AL,
- rlsOptIC.M,
- rlsOptIC.BM,
- rlsOptIC.Hu1,
- rlsOptIC.Hu2a,
- rlsOptIC.Hu2,
- rlsOptIC.Hu3,
- rlsOptIC.HuMad,
- rlsOptIC.Ha3,
- rlsOptIC.Ha4,
- rlsOptIC.HaMad,
- rlsOptIC.An1,
- rlsOptIC.An2,
- rlsOptIC.AnMad,
- rlsOptIC.Tu1,
- rlsOptIC.Tu2,
- rlsOptIC.TuMad,
- rlsOptIC.MM2,
- rlOptIC,
- rsOptIC,
- roblox,
- rowRoblox,
- colRoblox,
- showdown)
+importFrom("graphics", "abline", "curve", "layout", "legend", "par")
+importFrom("stats", "optim", "pnorm", "qnorm", "rbinom", "rnorm",
+ "uniroot")
+importFrom("utils", "packageDescription")
+exportMethods("robloxbioc", "KolmogorovMinDist")
+export("AffySimStudy", "IlluminaSimStudy")
Modified: pkg/RobLoxBioC/R/AffySimStudyFunction.R
===================================================================
--- pkg/RobLoxBioC/R/AffySimStudyFunction.R 2018-07-23 20:29:47 UTC (rev 1036)
+++ pkg/RobLoxBioC/R/AffySimStudyFunction.R 2018-07-23 20:31:10 UTC (rev 1037)
@@ -1,120 +1,120 @@
-###############################################################################
-## Function to perform simulation study comparing Tukey's biweight with
-## rmx estimators
-###############################################################################
-AffySimStudy <- function(n, M, eps, seed = 123, eps.lower = 0, eps.upper = 0.05,
- steps = 3L, fsCor = TRUE, contD,
- plot1 = FALSE, plot2 = FALSE, plot3 = FALSE){
- stopifnot(n >= 3)
- stopifnot(eps >= 0, eps <= 0.5)
- if(plot1){
- from <- min(-6, q(contD)(1e-15))
- to <- max(6, q(contD)(1-1e-15))
- curve(pnorm, from = from, to = to, lwd = 2, n = 201,
- main = "Comparison: ideal vs. real", ylab = "cdf")
- fun <- function(x) (1-eps)*pnorm(x) + eps*p(contD)(x)
- curve(fun, from = from, to = to, add = TRUE, col = "orange",
- lwd = 2, n = 201, ylab = "cdf")
- legend("topleft", legend = c("ideal", "real"),
- fill = c("black", "orange"))
- }
-
- set.seed(seed)
- r <- rbinom(n*M, prob = eps, size = 1)
- Mid <- rnorm(n*M)
- Mcont <- r(contD)(n*M)
- Mre <- matrix((1-r)*Mid + r*Mcont, ncol = n)
- ind <- rowSums(matrix(r, ncol = n)) >= n/2
- while(any(ind)){
- M1 <- sum(ind)
- cat("Samples to re-simulate:\t", M1, "\n")
- r <- rbinom(n*M1, prob = eps, size = 1)
- Mid <- rnorm(n*M1)
- Mcont <- r(contD)(n*M1)
- Mre[ind,] <- (1-r)*Mid + r*Mcont
- ind[ind] <- rowSums(matrix(r, ncol = n)) >= n/2
- }
- rm(Mid, Mcont, r, ind)
-
-
- if(plot2){
- ind <- if(M <= 20) 1:M else sample(1:M, 20)
- if(plot1) dev.new()
- M1 <- min(M, 20)
- print(
- stripplot(rep(1:M1, each = n) ~ as.vector(Mre[ind,]),
- ylab = "samples", xlab = "x", pch = 20,
- main = ifelse(M <= 20, "Samples", "20 randomly chosen samples"))
- )
- }
-
- ## ML-estimator: mean and sd
- Mean <- rowMeans(Mre)
- Sd <- sqrt(rowMeans((Mre-Mean)^2))
- ## Median and MAD
- Median <- rowMedians(Mre)
- Mad <- rowMedians(abs(Mre - Median))/qnorm(0.75)
- ## Tukey 1-step + MAD
- Tukey <- apply(Mre, 1, function(x) tukey.biweight(x))
- Tukey <- cbind(Tukey, Mad)
-
- ## Radius-minimax estimator
- RadMinmax <- estimate(rowRoblox(Mre, eps.lower = eps.lower,
- eps.upper = eps.upper, k = steps,
- fsCor = fsCor))
-
- if(plot3){
- Ergebnis1 <- list(Mean, Median, Tukey[,1], RadMinmax[,1])
- Ergebnis2 <- list(Sd, Mad, RadMinmax[,2])
- myCol <- brewer.pal(4, "Dark2")
- if(plot1 || plot2) dev.new()
- layout(matrix(c(1, 1, 1, 1, 3, 2, 2, 2, 2, 3), ncol = 2))
- boxplot(Ergebnis1, col = myCol, pch = 20, main = "Location")
- abline(h = 0)
- boxplot(Ergebnis2, col = myCol[c(1,2,4)], pch = 20, main = "Scale")
- abline(h = 1)
- op <- par(mar = rep(2, 4))
- plot(c(0,1), c(1, 0), type = "n", axes = FALSE)
- legend("center", c("ML", "Med/MAD", "biweight", "rmx"),
- fill = myCol, ncol = 4, cex = 1.5)
- on.exit(par(op))
- }
-
- ## ML-estimator
- MSE1.1 <- n*mean(Mean^2)
- ## Median + MAD
- MSE2.1 <- n*mean(Median^2)
- ## Tukey
- MSE3.1 <- n*mean(Tukey[,1]^2)
- ## Radius-minimax
- MSE4.1 <- n*mean(RadMinmax[,1]^2)
- empMSE <- data.frame(ML = MSE1.1, Med = MSE2.1, Tukey = MSE3.1, "rmx" = MSE4.1)
- rownames(empMSE) <- "n x empMSE (loc)"
- relMSE <- empMSE[1,]/empMSE[1,4]
- empMSE <- rbind(empMSE, relMSE)
- rownames(empMSE)[2] <- "relMSE (loc)"
-
- ## ML-estimator
- MSE1.2 <- n*mean((Sd-1)^2)
- ## Median + MAD
- MSE2.2 <- n*mean((Mad-1)^2)
- ## Tukey
- MSE3.2 <- MSE2.2
- ## Radius-minimax
- MSE4.2 <- n*mean((RadMinmax[,2]-1)^2)
- empMSE <- rbind(empMSE, c(MSE1.2, MSE2.2, MSE3.2, MSE4.2))
- rownames(empMSE)[3] <- "n x empMSE (scale)"
- relMSE <- empMSE[3,]/empMSE[3,4]
- empMSE <- rbind(empMSE, relMSE)
- rownames(empMSE)[4] <- "relMSE (scale)"
- empMSE <- rbind(empMSE, c(MSE1.1 + MSE1.2, MSE2.1 + MSE2.2, MSE3.1 + MSE3.2,
- MSE4.1 + MSE4.2))
- rownames(empMSE)[5] <- "n x empMSE (loc + scale)"
- relMSE <- empMSE[5,]/empMSE[5,4]
- empMSE <- rbind(empMSE, relMSE)
- rownames(empMSE)[6] <- "relMSE (loc + scale)"
-
- empMSE
-}
-
-
+###############################################################################
+## Function to perform simulation study comparing Tukey's biweight with
+## rmx estimators
+###############################################################################
+AffySimStudy <- function(n, M, eps, seed = 123, eps.lower = 0, eps.upper = 0.05,
+ steps = 3L, fsCor = TRUE, contD,
+ plot1 = FALSE, plot2 = FALSE, plot3 = FALSE){
+ stopifnot(n >= 3)
+ stopifnot(eps >= 0, eps <= 0.5)
+ if(plot1){
+ from <- min(-6, q.l(contD)(1e-15))
+ to <- max(6, q.l(contD)(1-1e-15))
+ curve(pnorm, from = from, to = to, lwd = 2, n = 201,
+ main = "Comparison: ideal vs. real", ylab = "cdf")
+ fun <- function(x) (1-eps)*pnorm(x) + eps*p(contD)(x)
+ curve(fun, from = from, to = to, add = TRUE, col = "orange",
+ lwd = 2, n = 201, ylab = "cdf")
+ legend("topleft", legend = c("ideal", "real"),
+ fill = c("black", "orange"))
+ }
+
+ set.seed(seed)
+ r <- rbinom(n*M, prob = eps, size = 1)
+ Mid <- rnorm(n*M)
+ Mcont <- r(contD)(n*M)
+ Mre <- matrix((1-r)*Mid + r*Mcont, ncol = n)
+ ind <- rowSums(matrix(r, ncol = n)) >= n/2
+ while(any(ind)){
+ M1 <- sum(ind)
+ cat("Samples to re-simulate:\t", M1, "\n")
+ r <- rbinom(n*M1, prob = eps, size = 1)
+ Mid <- rnorm(n*M1)
+ Mcont <- r(contD)(n*M1)
+ Mre[ind,] <- (1-r)*Mid + r*Mcont
+ ind[ind] <- rowSums(matrix(r, ncol = n)) >= n/2
+ }
+ rm(Mid, Mcont, r, ind)
+
+
+ if(plot2){
+ ind <- if(M <= 20) 1:M else sample(1:M, 20)
+ if(plot1) dev.new()
+ M1 <- min(M, 20)
+ print(
+ stripplot(rep(1:M1, each = n) ~ as.vector(Mre[ind,]),
+ ylab = "samples", xlab = "x", pch = 20,
+ main = ifelse(M <= 20, "Samples", "20 randomly chosen samples"))
+ )
+ }
+
+ ## ML-estimator: mean and sd
+ Mean <- rowMeans(Mre)
+ Sd <- sqrt(rowMeans((Mre-Mean)^2))
+ ## Median and MAD
+ Median <- rowMedians(Mre)
+ Mad <- rowMedians(abs(Mre - Median))/qnorm(0.75)
+ ## Tukey 1-step + MAD
+ Tukey <- apply(Mre, 1, function(x) tukey.biweight(x))
+ Tukey <- cbind(Tukey, Mad)
+
+ ## Radius-minimax estimator
+ RadMinmax <- estimate(rowRoblox(Mre, eps.lower = eps.lower,
+ eps.upper = eps.upper, k = steps,
+ fsCor = fsCor))
+
+ if(plot3){
+ Ergebnis1 <- list(Mean, Median, Tukey[,1], RadMinmax[,1])
+ Ergebnis2 <- list(Sd, Mad, RadMinmax[,2])
+ myCol <- brewer.pal(4, "Dark2")
+ if(plot1 || plot2) dev.new()
+ layout(matrix(c(1, 1, 1, 1, 3, 2, 2, 2, 2, 3), ncol = 2))
+ boxplot(Ergebnis1, col = myCol, pch = 20, main = "Location")
+ abline(h = 0)
+ boxplot(Ergebnis2, col = myCol[c(1,2,4)], pch = 20, main = "Scale")
+ abline(h = 1)
+ op <- par(mar = rep(2, 4))
+ plot(c(0,1), c(1, 0), type = "n", axes = FALSE)
+ legend("center", c("ML", "Med/MAD", "biweight", "rmx"),
+ fill = myCol, ncol = 4, cex = 1.5)
+ on.exit(par(op))
+ }
+
+ ## ML-estimator
+ MSE1.1 <- n*mean(Mean^2)
+ ## Median + MAD
+ MSE2.1 <- n*mean(Median^2)
+ ## Tukey
+ MSE3.1 <- n*mean(Tukey[,1]^2)
+ ## Radius-minimax
+ MSE4.1 <- n*mean(RadMinmax[,1]^2)
+ empMSE <- data.frame(ML = MSE1.1, Med = MSE2.1, Tukey = MSE3.1, "rmx" = MSE4.1)
+ rownames(empMSE) <- "n x empMSE (loc)"
+ relMSE <- empMSE[1,]/empMSE[1,4]
+ empMSE <- rbind(empMSE, relMSE)
+ rownames(empMSE)[2] <- "relMSE (loc)"
+
+ ## ML-estimator
+ MSE1.2 <- n*mean((Sd-1)^2)
+ ## Median + MAD
+ MSE2.2 <- n*mean((Mad-1)^2)
+ ## Tukey
+ MSE3.2 <- MSE2.2
+ ## Radius-minimax
+ MSE4.2 <- n*mean((RadMinmax[,2]-1)^2)
+ empMSE <- rbind(empMSE, c(MSE1.2, MSE2.2, MSE3.2, MSE4.2))
+ rownames(empMSE)[3] <- "n x empMSE (scale)"
+ relMSE <- empMSE[3,]/empMSE[3,4]
+ empMSE <- rbind(empMSE, relMSE)
+ rownames(empMSE)[4] <- "relMSE (scale)"
+ empMSE <- rbind(empMSE, c(MSE1.1 + MSE1.2, MSE2.1 + MSE2.2, MSE3.1 + MSE3.2,
+ MSE4.1 + MSE4.2))
+ rownames(empMSE)[5] <- "n x empMSE (loc + scale)"
+ relMSE <- empMSE[5,]/empMSE[5,4]
+ empMSE <- rbind(empMSE, relMSE)
+ rownames(empMSE)[6] <- "relMSE (loc + scale)"
+
+ empMSE
+}
+
+
Modified: pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R
===================================================================
--- pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R 2018-07-23 20:29:47 UTC (rev 1036)
+++ pkg/RobLoxBioC/R/IlluminaSimStudyFunction.R 2018-07-23 20:31:10 UTC (rev 1037)
@@ -1,128 +1,128 @@
-###############################################################################
-## Function to perform simulation study comparing Illumina's default method
-## with rmx estimators
-###############################################################################
-IlluminaSimStudy <- function(n, M, eps, seed = 123,
- eps.lower = 0, eps.upper = 0.05,
- steps = 3L, fsCor = TRUE, contD,
- plot1 = FALSE, plot2 = FALSE, plot3 = FALSE){
- stopifnot(n >= 3)
- stopifnot(eps >= 0, eps <= 0.5)
- if(plot1){
- from <- min(-6, q(contD)(1e-15))
- to <- max(6, q(contD)(1-1e-15))
- curve(pnorm, from = from, to = to, lwd = 2, n = 201,
- main = "Comparison: ideal vs. real", ylab = "cdf")
- fun <- function(x) (1-eps)*pnorm(x) + eps*p(contD)(x)
- curve(fun, from = from, to = to, add = TRUE, col = "orange",
- lwd = 2, n = 201, ylab = "cdf")
- legend("topleft", legend = c("ideal", "real"),
- fill = c("black", "orange"))
- }
-
- set.seed(seed)
- r <- rbinom(n*M, prob = eps, size = 1)
- Mid <- rnorm(n*M)
- Mcont <- r(contD)(n*M)
- Mre <- matrix((1-r)*Mid + r*Mcont, ncol = n)
- ind <- rowSums(matrix(r, ncol = n)) >= n/2
- while(any(ind)){
- M1 <- sum(ind)
- cat("Samples to re-simulate:\t", M1, "\n")
- r <- rbinom(n*M1, prob = eps, size = 1)
- Mid <- rnorm(n*M1)
- Mcont <- r(contD)(n*M1)
- Mre[ind,] <- (1-r)*Mid + r*Mcont
- ind[ind] <- rowSums(matrix(r, ncol = n)) >= n/2
- }
- rm(Mid, Mcont, r, ind)
-
-
- if(plot2){
- ind <- if(M <= 20) 1:M else sample(1:M, 20)
- if(plot1) dev.new()
- M1 <- min(M, 20)
- print(
- stripplot(rep(1:M1, each = n) ~ as.vector(Mre[ind,]),
- ylab = "samples", xlab = "x", pch = 20,
- main = ifelse(M <= 20, "Samples", "20 randomly chosen samples"))
- )
- }
-
- ## ML-estimator: mean and sd
- Mean <- rowMeans(Mre)
- Sd <- sqrt(rowMeans((Mre-Mean)^2))
- ## Median and MAD
- Median <- rowMedians(Mre)
- Mad <- rowMedians(abs(Mre - Median))/qnorm(0.75)
- ## Illumina method
- ind <- (Mre < (Median - 3*Mad)) | (Mre > (Median + 3*Mad))
- x.ill <- Mre
- x.ill[ind] <- NA
- Illum.mean <- rowMeans(x.ill, na.rm = TRUE)
- n.ill <- rowSums(!is.na(x.ill))
- n.ill[n.ill < 1] <- NA
- Illum.sd <- sqrt(rowSums((x.ill - Illum.mean)^2, na.rm = TRUE)/(n.ill-1))
- Illum <- cbind(Illum.mean, Illum.sd)
-
- ## Radius-minimax estimator
- RadMinmax <- estimate(rowRoblox(Mre, eps.lower = eps.lower,
- eps.upper = eps.upper, k = steps,
- fsCor = fsCor))
-
- if(plot3){
- Ergebnis1 <- list(Mean, Median, Illum[,1], RadMinmax[,1])
- Ergebnis2 <- list(Sd, Mad, Illum[,2], RadMinmax[,2])
- myCol <- brewer.pal(4, "Dark2")
- if(plot1 || plot2) dev.new()
- layout(matrix(c(1, 1, 1, 1, 3, 2, 2, 2, 2, 3), ncol = 2))
- boxplot(Ergebnis1, col = myCol, pch = 20, main = "Location")
- abline(h = 0)
- boxplot(Ergebnis2, col = myCol, pch = 20, main = "Scale")
- abline(h = 1)
- op <- par(mar = rep(2, 4))
- plot(c(0,1), c(1, 0), type = "n", axes = FALSE)
- legend("center", c("ML", "Med/MAD", "Illumina", "rmx"),
- fill = myCol, ncol = 4, cex = 1.5)
- op$cin <- op$cra <- op$csi <- op$cxy <- op$din <- NULL
- on.exit(par(op))
- }
-
- ## ML-estimator
- MSE1.1 <- n*mean(Mean^2)
- ## Median + MAD
- MSE2.1 <- n*mean(Median^2)
- ## Illumina's default method
- MSE3.1 <- n*mean(Illum[,1]^2)
- ## Radius-minimax
- MSE4.1 <- n*mean(RadMinmax[,1]^2)
- empMSE <- data.frame(ML = MSE1.1, Med = MSE2.1, Illumina = MSE3.1, "rmx" = MSE4.1)
- rownames(empMSE) <- "n x empMSE (loc)"
- relMSE <- empMSE[1,]/empMSE[1,4]
- empMSE <- rbind(empMSE, relMSE)
- rownames(empMSE)[2] <- "relMSE (loc)"
-
- ## ML-estimator
- MSE1.2 <- n*mean((Sd-1)^2)
- ## Median + MAD
- MSE2.2 <- n*mean((Mad-1)^2)
- ## Illumina's default method
- MSE3.2 <- n*mean((Illum[,2]-1)^2)
- ## Radius-minimax
- MSE4.2 <- n*mean((RadMinmax[,2]-1)^2)
- empMSE <- rbind(empMSE, c(MSE1.2, MSE2.2, MSE3.2, MSE4.2))
- rownames(empMSE)[3] <- "n x empMSE (scale)"
- relMSE <- empMSE[3,]/empMSE[3,4]
- empMSE <- rbind(empMSE, relMSE)
- rownames(empMSE)[4] <- "relMSE (scale)"
- empMSE <- rbind(empMSE, c(MSE1.1 + MSE1.2, MSE2.1 + MSE2.2, MSE3.1 + MSE3.2,
- MSE4.1 + MSE4.2))
- rownames(empMSE)[5] <- "n x empMSE (loc + scale)"
- relMSE <- empMSE[5,]/empMSE[5,4]
- empMSE <- rbind(empMSE, relMSE)
- rownames(empMSE)[6] <- "relMSE (loc + scale)"
-
- empMSE
-}
-
-
+###############################################################################
+## Function to perform simulation study comparing Illumina's default method
+## with rmx estimators
+###############################################################################
+IlluminaSimStudy <- function(n, M, eps, seed = 123,
+ eps.lower = 0, eps.upper = 0.05,
+ steps = 3L, fsCor = TRUE, contD,
+ plot1 = FALSE, plot2 = FALSE, plot3 = FALSE){
+ stopifnot(n >= 3)
+ stopifnot(eps >= 0, eps <= 0.5)
+ if(plot1){
+ from <- min(-6, q.l(contD)(1e-15))
+ to <- max(6, q.l(contD)(1-1e-15))
+ curve(pnorm, from = from, to = to, lwd = 2, n = 201,
+ main = "Comparison: ideal vs. real", ylab = "cdf")
+ fun <- function(x) (1-eps)*pnorm(x) + eps*p(contD)(x)
+ curve(fun, from = from, to = to, add = TRUE, col = "orange",
+ lwd = 2, n = 201, ylab = "cdf")
+ legend("topleft", legend = c("ideal", "real"),
+ fill = c("black", "orange"))
+ }
+
+ set.seed(seed)
+ r <- rbinom(n*M, prob = eps, size = 1)
+ Mid <- rnorm(n*M)
+ Mcont <- r(contD)(n*M)
+ Mre <- matrix((1-r)*Mid + r*Mcont, ncol = n)
+ ind <- rowSums(matrix(r, ncol = n)) >= n/2
+ while(any(ind)){
+ M1 <- sum(ind)
+ cat("Samples to re-simulate:\t", M1, "\n")
+ r <- rbinom(n*M1, prob = eps, size = 1)
+ Mid <- rnorm(n*M1)
+ Mcont <- r(contD)(n*M1)
+ Mre[ind,] <- (1-r)*Mid + r*Mcont
+ ind[ind] <- rowSums(matrix(r, ncol = n)) >= n/2
+ }
+ rm(Mid, Mcont, r, ind)
+
+
+ if(plot2){
+ ind <- if(M <= 20) 1:M else sample(1:M, 20)
+ if(plot1) dev.new()
+ M1 <- min(M, 20)
+ print(
+ stripplot(rep(1:M1, each = n) ~ as.vector(Mre[ind,]),
+ ylab = "samples", xlab = "x", pch = 20,
+ main = ifelse(M <= 20, "Samples", "20 randomly chosen samples"))
+ )
+ }
+
+ ## ML-estimator: mean and sd
+ Mean <- rowMeans(Mre)
+ Sd <- sqrt(rowMeans((Mre-Mean)^2))
+ ## Median and MAD
+ Median <- rowMedians(Mre)
+ Mad <- rowMedians(abs(Mre - Median))/qnorm(0.75)
+ ## Illumina method
+ ind <- (Mre < (Median - 3*Mad)) | (Mre > (Median + 3*Mad))
+ x.ill <- Mre
+ x.ill[ind] <- NA
+ Illum.mean <- rowMeans(x.ill, na.rm = TRUE)
+ n.ill <- rowSums(!is.na(x.ill))
+ n.ill[n.ill < 1] <- NA
+ Illum.sd <- sqrt(rowSums((x.ill - Illum.mean)^2, na.rm = TRUE)/(n.ill-1))
+ Illum <- cbind(Illum.mean, Illum.sd)
+
+ ## Radius-minimax estimator
+ RadMinmax <- estimate(rowRoblox(Mre, eps.lower = eps.lower,
+ eps.upper = eps.upper, k = steps,
+ fsCor = fsCor))
+
+ if(plot3){
+ Ergebnis1 <- list(Mean, Median, Illum[,1], RadMinmax[,1])
+ Ergebnis2 <- list(Sd, Mad, Illum[,2], RadMinmax[,2])
+ myCol <- brewer.pal(4, "Dark2")
+ if(plot1 || plot2) dev.new()
+ layout(matrix(c(1, 1, 1, 1, 3, 2, 2, 2, 2, 3), ncol = 2))
+ boxplot(Ergebnis1, col = myCol, pch = 20, main = "Location")
+ abline(h = 0)
+ boxplot(Ergebnis2, col = myCol, pch = 20, main = "Scale")
+ abline(h = 1)
+ op <- par(mar = rep(2, 4))
+ plot(c(0,1), c(1, 0), type = "n", axes = FALSE)
+ legend("center", c("ML", "Med/MAD", "Illumina", "rmx"),
+ fill = myCol, ncol = 4, cex = 1.5)
+ op$cin <- op$cra <- op$csi <- op$cxy <- op$din <- NULL
+ on.exit(par(op))
+ }
+
+ ## ML-estimator
+ MSE1.1 <- n*mean(Mean^2)
+ ## Median + MAD
+ MSE2.1 <- n*mean(Median^2)
+ ## Illumina's default method
+ MSE3.1 <- n*mean(Illum[,1]^2)
+ ## Radius-minimax
+ MSE4.1 <- n*mean(RadMinmax[,1]^2)
+ empMSE <- data.frame(ML = MSE1.1, Med = MSE2.1, Illumina = MSE3.1, "rmx" = MSE4.1)
+ rownames(empMSE) <- "n x empMSE (loc)"
+ relMSE <- empMSE[1,]/empMSE[1,4]
+ empMSE <- rbind(empMSE, relMSE)
+ rownames(empMSE)[2] <- "relMSE (loc)"
+
+ ## ML-estimator
+ MSE1.2 <- n*mean((Sd-1)^2)
+ ## Median + MAD
+ MSE2.2 <- n*mean((Mad-1)^2)
+ ## Illumina's default method
+ MSE3.2 <- n*mean((Illum[,2]-1)^2)
+ ## Radius-minimax
+ MSE4.2 <- n*mean((RadMinmax[,2]-1)^2)
+ empMSE <- rbind(empMSE, c(MSE1.2, MSE2.2, MSE3.2, MSE4.2))
+ rownames(empMSE)[3] <- "n x empMSE (scale)"
+ relMSE <- empMSE[3,]/empMSE[3,4]
+ empMSE <- rbind(empMSE, relMSE)
+ rownames(empMSE)[4] <- "relMSE (scale)"
+ empMSE <- rbind(empMSE, c(MSE1.1 + MSE1.2, MSE2.1 + MSE2.2, MSE3.1 + MSE3.2,
+ MSE4.1 + MSE4.2))
+ rownames(empMSE)[5] <- "n x empMSE (loc + scale)"
+ relMSE <- empMSE[5,]/empMSE[5,4]
+ empMSE <- rbind(empMSE, relMSE)
+ rownames(empMSE)[6] <- "relMSE (loc + scale)"
+
+ empMSE
+}
+
+
Modified: pkg/RobLoxBioC/R/robloxbiocBeadLevelData.R
===================================================================
--- pkg/RobLoxBioC/R/robloxbiocBeadLevelData.R 2018-07-23 20:29:47 UTC (rev 1036)
+++ pkg/RobLoxBioC/R/robloxbiocBeadLevelData.R 2018-07-23 20:31:10 UTC (rev 1037)
@@ -181,7 +181,7 @@
rownames(eMat) <- rownames(varMat) <- rownames(nObs) <- as.character(IlluminaIDs)
status <- rep("Unknown", length(probeIDs))
annoPkg <- paste("illumina", annoName, ".db", sep = "")
- annoVers <- packageDescription(annoPkg, field = "Version")
+ annoVers <- packageDescription(annoPkg, fields = "Version")
message(paste("Annotating control probes using package ", annoPkg, " Version:", annoVers, "\n", sep = ""))
mapEnv <- as.name(paste("illumina", annoName, "REPORTERGROUPNAME", sep = ""))
t <- try(eval(mapEnv), silent = TRUE)
Modified: pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd
===================================================================
--- pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd 2018-07-23 20:29:47 UTC (rev 1036)
+++ pkg/RobLoxBioC/man/0RobLoxBioC-package.Rd 2018-07-23 20:31:10 UTC (rev 1037)
@@ -12,15 +12,15 @@
\details{
\tabular{ll}{
Package: \tab RobLoxBioC \cr
-Version: \tab 0.9 \cr
-Date: \tab 2013-09-12 \cr
-Depends: \tab R(>= 2.14.0), methods, Biobase, affy, beadarray, distr, RobLox,
-lattice, RColorBrewer \cr
-LazyLoad: \tab yes \cr
+Version: \tab 1.1.0 \cr
+Date: \tab 2018-07-17 \cr
+Depends:\tab R(>= 2.14.0), methods, distr(>= 2.5.2), affy \cr
+Imports:\tab Biobase, BiocGenerics, beadarray, RobLox(>= 0.9.2), distrMod(>= 2.5.2), lattice, RColorBrewer \cr
+Suggests:\tab affydata, hgu95av2cdf, beadarrayExampleData, illuminaHumanv3.db \cr
ByteCompile: \tab yes \cr
License: \tab LGPL-3 \cr
URL: \tab http://robast.r-forge.r-project.org/\cr
-SVNRevision: \tab 696 \cr
+VCS/SVNRevision: \tab 940 \cr
Encoding: \tab latin1 \cr
}
}
@@ -43,7 +43,7 @@
Rieder, H., Kohl, M. and Ruckdeschel, P. (2008) The Costs of not Knowing
the Radius. \emph{Statistical Methods and Applications} \bold{17}(1) 13-40.
- Extended version: \url{http://www.stamats.de/RRlong.pdf}
+ Extended version: \url{http://r-kurs.de/RRlong.pdf}
}
\seealso{
\code{\link[RobLox]{roblox}}, \code{\link[RobLox]{rowRoblox}}
Modified: pkg/RobLoxBioC/man/robloxbioc.Rd
===================================================================
--- pkg/RobLoxBioC/man/robloxbioc.Rd 2018-07-23 20:29:47 UTC (rev 1036)
+++ pkg/RobLoxBioC/man/robloxbioc.Rd 2018-07-23 20:31:10 UTC (rev 1037)
@@ -125,7 +125,7 @@
Rieder, H., Kohl, M. and Ruckdeschel, P. (2008) The Costs of not Knowing
the Radius. \emph{Statistical Methods and Applications} \bold{17}(1) 13-40.
- Extended version: \url{http://www.stamats.de/RRlong.pdf}
+ Extended version: \url{http://r-kurs.de/RRlong.pdf}
}
\author{Matthias Kohl \email{Matthias.Kohl at stamats.de},
@@ -147,13 +147,13 @@
robloxbioc(X, eps = 0.05)
robloxbioc(X, eps = 0.05, steps = 5)
-## Don't run to reduce check time on CRAN
-\dontrun{
+\donttest{
+## \donttest to reduce check time
## the function is designed for large scale problems
X <- matrix(rnorm(50000*20, mean = 1), nrow = 50000)
system.time(robloxbioc(X))
-## using Affymetrix-Data
+## using Affymetrix data
## confer example to generateExprVal.method.mas
## A more worked out example can be found in the scripts folder
## of the package.
@@ -167,13 +167,15 @@
points(concentrations, rl[,1], pch = 20, col="orange", type="b")
legend("topleft", c("MAS", "roblox"), pch = c(1, 20))
-require(affydata)
+## Affymetrix dilution data
+library(affydata)
data(Dilution)
eset <- robloxbioc(Dilution)
## Affymetrix scale normalization
eset1 <- robloxbioc(Dilution, normalize = TRUE)
-## "Not run" just because of computation time
-require(beadarrayExampleData)
+
+## Illumina bead level data
+library(beadarrayExampleData)
data(exampleBLData)
res <- robloxbioc(exampleBLData, eps.upper = 0.5)
res
More information about the Robast-commits
mailing list