[Robast-commits] r276 - in pkg/RobLoxBioC: R inst/scripts man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Apr 8 16:41:53 CEST 2009
Author: stamats
Date: 2009-04-08 16:41:53 +0200 (Wed, 08 Apr 2009)
New Revision: 276
Modified:
pkg/RobLoxBioC/R/KolmogorovMinDist.R
pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R
pkg/RobLoxBioC/man/KolmogorovMinDist.Rd
Log:
changed return value of KolmogorovMinDist, now also returns the corresponding sample sizes.
Modified: pkg/RobLoxBioC/R/KolmogorovMinDist.R
===================================================================
--- pkg/RobLoxBioC/R/KolmogorovMinDist.R 2009-04-08 10:51:15 UTC (rev 275)
+++ pkg/RobLoxBioC/R/KolmogorovMinDist.R 2009-04-08 14:41:53 UTC (rev 276)
@@ -15,19 +15,19 @@
Med <- rowMedians(x, na.rm = TRUE)
Mad <- pmax(rowMedians(abs(x-Med), na.rm = TRUE)/qnorm(0.75), mad0)
startPars <- cbind(Med, Mad)
- n <- nrow(x)
- res <- numeric(n)
- for(i in seq_len(n)){
+ M <- nrow(x)
+ n <- ncol(x)
+ res <- matrix(NA, nrow = M, ncol = 2)
+ for(i in seq_len(M)){
temp <- try(optim(par = startPars[i,], fn = ksdist, x = x[i,],
method = "L-BFGS-B", lower = c(-Inf, 1e-15),
- upper = c(Inf, Inf))$value)
- if(inherits(temp, "try-error")){
- res[i] <- NA
- }else{
- res[i] <- temp
+ upper = c(Inf, Inf))$value, silent = TRUE)
+ if(!inherits(temp, "try-error")){
+ res[i,1] <- temp
+ res[i,2] <- n
}
}
- res
+ list("dist" = res[,1], "n" = res[,2])
})
setMethod("KolmogorovMinDist", signature(x = "AffyBatch",
@@ -49,6 +49,7 @@
intensData <- intensity(x)
kmd <- numeric(m*n)
+ ns <- numeric(m*n)
if(pmcorrect){
if(verbose) cat("Calculating PM/MM ...")
div <- function(INDEX, x){
@@ -75,12 +76,16 @@
ind <- which(NROW == k)
temp <- matrix(do.call(rbind, res[ind]), nrow = k)
ind1 <- as.vector(sapply(seq_len(n)-1, function(x, ind, m){ ind + x*m }, ind = ind, m = m))
- kmd[ind1] <- KolmogorovMinDist(log2(t(temp)), D = D)
+ temp.res <- KolmogorovMinDist(log2(t(temp)), D = D)
+ kmd[ind1] <- temp.res[["dist"]]
+ ns[ind1] <- temp.res[["n"]]
}
if(verbose) cat(" done.\n")
kmd.mat <- matrix(kmd, nrow = m)
+ ns.mat <- matrix(ns, nrow = m)
dimnames(kmd.mat) <- list(ids, sampleNames(x))
- kmd.mat
+ dimnames(ns.mat) <- list(ids, sampleNames(x))
+ list("dist" = kmd.mat, "n" = ns.mat)
})
setMethod("KolmogorovMinDist", signature(x = "BeadLevelList",
@@ -138,23 +143,27 @@
pr <- pr[!nasinf]
if (imagesPerArray == 1) {
G.kmd <- matrix(0, nrow = noprobes, ncol = len)
- colnames(G.kmd) <- arraynms
+ G.ns <- matrix(0, nrow = noprobes, ncol = len)
+ colnames(G.kmd) <- colnames(G.ns) <- arraynms
if (BLData at arrayInfo$channels == "two" && !is.null(BLData[[arraynms[1]]]$R) && whatelse == "R")
- R.kmd <- G.kmd
- else R.kmd <- NULL
+ R.kmd <- R.ns <- G.kmd
+ else R.kmd <- R.ns <- NULL
}
else if (imagesPerArray == 2) {
G.kmd <- matrix(0, nrow = noprobes, ncol = (len/2))
- colnames(G.kmd) <- arraynms[seq(1, len, by = 2)]
+ G.ns <- matrix(0, nrow = noprobes, ncol = (len/2))
+ colnames(G.kmd) <- colnames(G.ns) <- arraynms[seq(1, len, by = 2)]
if (BLData at arrayInfo$channels == "two" && !is.null(BLData[[arraynms[1]]]$R) && whatelse == "R")
- R.kmd <- G.kmd
- else R.kmd <- NULL
+ R.kmd <- R.ns <- G.kmd
+ else R.kmd <- R.ns <- NULL
}
i <- j <- 1
while (j <= len) {
probeIDs <- as.integer(pr)
start <- 0
- G.kmd[,i] <- kmdBeadLevel(x = finten, D = D, probeIDs = probeIDs, probes = probes)
+ G.res <- kmdBeadLevel(x = finten, D = D, probeIDs = probeIDs, probes = probes)
+ G.kmd[,i] <- G.res[["dist"]]
+ G.ns[,i] <- G.res[["n"]]
if (BLData at arrayInfo$channels == "two" && !is.null(BLData[[arraynms[i]]]$R) && whatelse == "R") {
if (imagesPerArray == 1) {
finten <- getArrayData(BLData, what = whatelse, log = log, array = arraynms[i])[sel]
@@ -170,7 +179,9 @@
binten <- rep(0, length(finten))
}
start <- 0
- R.kmd[,i] <- kmdBeadLevel(x = finten, D = D, probeIDs = probeIDs, probes = probes)
+ R.res[,i] <- kmdBeadLevel(x = finten, D = D, probeIDs = probeIDs, probes = probes)
+ R.kmd[,i] <- R.res[["dist"]]
+ R.ns[,i] <- R.res[["n"]]
}
j <- j + imagesPerArray
i <- i + 1
@@ -197,12 +208,13 @@
}
}
if (whatelse == "R") {
- rownames(G.kmd) <- rownames(R.kmd) <- probes
- res <- list(G = G.kmd, R = R.kmd)
+ rownames(G.kmd) <- rownames(G.ns) <- rownames(R.kmd) <- rownames(R.ns) <- probes
+ res <- list(G = list("dist" = G.kmd, "n" = G.ns),
+ R = list("dist" = R.kmd, "n" = R.ns))
}
else {
- rownames(G.kmd) <- probes
- res <- G.kmd
+ rownames(G.kmd) <- rownames(G.ns) <- probes
+ res <- list("dist" = G.kmd, "n" = G.ns)
}
return(res)
})
@@ -215,6 +227,7 @@
probes1 <- comIDs
len1 <- length(probes1)
kmd1 <- numeric(len1)
+ ns1 <- numeric(len1)
for(i in seq(along = noBeads.uni)){
index <- noBeads == noBeads.uni[i]
IDs <- probes1[index]
@@ -222,14 +235,19 @@
kmd1[index] <- 0.5
}else{
temp <- matrix(x[probeIDs %in% IDs], ncol = noBeads.uni[i], byrow = TRUE)
- kmd1[index] <- KolmogorovMinDist(temp, D = D)
+ res <- KolmogorovMinDist(temp, D = D)
+ kmd1[index] <- res[["dist"]]
+ ns1[index] <- res[["n"]]
}
}
len <- length(probes)
kmd <- numeric(len)
+ ns <- numeric(len)
nas <- !(probes %in% comIDs)
kmd[nas] <- NA
kmd[!nas] <- kmd1
+ ns[nas] <- NA
+ ns[!nas] <- ns1
- return(kmd)
+ list("dist" = kmd, "n" = ns)
}
Modified: pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R
===================================================================
--- pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R 2009-04-08 10:51:15 UTC (rev 275)
+++ pkg/RobLoxBioC/inst/scripts/AffymetrixExample.R 2009-04-08 14:41:53 UTC (rev 276)
@@ -75,10 +75,18 @@
## takes more than 130 min on Intel P9500 (64bit Linux, 4 GByte RAM)
#system.time(minKD.hgu133a <- KolmogorovMinDist(spikein.hgu133a, Norm()))
-## load the results from r-forge ...
-load(file = "https://r-forge.r-project.org/minKD_hgu95a.RData")
-(minKD.hgu133a)
+## load the results from R-forge ...
+con <- url("http://robast.r-forge.r-project.org/data/minKD_hgu95a.RData")
+load(file = con)
+close(con)
+con <- url("http://robast.r-forge.r-project.org/data/minKD_hgu133a.RData")
+load(file = con)
+close(con)
+boxplot(as.data.frame(minKD.hgu95a$dist), main = "HGU95a")
+boxplot(as.data.frame(minKD.hgu133a$dist), main = "HGU133a")
+
+
###########################################################
## assessments for MAS 5.0 and RMA including dilution data from package affycomp
###########################################################
Modified: pkg/RobLoxBioC/man/KolmogorovMinDist.Rd
===================================================================
--- pkg/RobLoxBioC/man/KolmogorovMinDist.Rd 2009-04-08 10:51:15 UTC (rev 275)
+++ pkg/RobLoxBioC/man/KolmogorovMinDist.Rd 2009-04-08 14:41:53 UTC (rev 276)
@@ -49,7 +49,9 @@
So far, only the minimum distance to the set of normal distributions can be
computed.
}
-\value{ Numeric vector with minimum Kolmogorov distances. }
+\value{ List with components \code{dist} containing a numeric vector
+or matrix with minimum Kolmogorov distances and \code{n} a numeric vector
+or matrix with the corresponding sample sizes. }
\references{
Huber, P.J. (1981) \emph{Robust Statistics}. New York: Wiley.
@@ -68,15 +70,21 @@
data(SpikeIn)
probes <- log2(pm(SpikeIn))
(res <- KolmogorovMinDist(probes, Norm()))
-boxplot(res)
+boxplot(res$dist)
\dontrun{
## "Not run" just because of computation time
require(affydata)
data(Dilution)
res <- KolmogorovMinDist(Dilution[,1], Norm())
-summary(res)
-boxplot(res)
+summary(res$dist)
+boxplot(res$dist)
+plot(res$n, res$dist, pch = 20, main = "Kolmogorov distance vs. sample size",
+ xlab = "sample size", ylab = "Kolmogorov distance",
+ ylim = c(0, max(res$dist)))
+uni.n <- min(res$n):max(res$n)
+lines(uni.n, 1/(2*uni.n), col = "orange", lwd = 2)
+legend("topright", legend = "minimal possible distance", fill = "orange")
}
## using Illumina-Data
@@ -85,8 +93,15 @@
data(BLData)
res <- KolmogorovMinDist(BLData, Norm(), arrays = 1)
res1 <- KolmogorovMinDist(BLData, log = TRUE, Norm(), arrays = 1)
-summary(cbind(res, res1))
-boxplot(list(res, res1), names = c("raw", "log-raw"))
+summary(cbind(res$dist, res1$dist))
+boxplot(list(res$dist, res1$dist), names = c("raw", "log-raw"))
+sort(unique(res1$n))
+plot(res1$n, res1$dist, pch = 20, main = "Kolmogorov distance vs. sample size",
+ xlab = "sample size", ylab = "Kolmogorov distance",
+ ylim = c(0, max(res1$dist)), xlim = c(min(res1$n), 56))
+uni.n <- min(res1$n):56
+lines(uni.n, 1/(2*uni.n), col = "orange", lwd = 2)
+legend("topright", legend = "minimal possible distance", fill = "orange")
}
}
\concept{normal location and scale}
More information about the Robast-commits
mailing list