Increasing rlm iterations for convergence in limma's normalizeRobustSpline
0
1
Entering edit mode
Paul Boutros ▴ 340
@paul-boutros-371
Last seen 9.5 years ago
Hi again, I've come across some cDNA array data that doesn't converge when I use limma's robust-spline normalization. I get this from: > library(limma); > setwd("c:/paul/dev/rat-mouse/limma/"); > targets <- readTargets(); > design <- cbind(c(1,1,1,1,-1,-1,-1,-1), c(0,0,0,0,1,1,1,1)); > colnames(design) <- c("Species", "Dye"); > read.series(targets$FileName,suffix=NULL,skip=29,sep="\t"); [1] "8 slides read" > layout <- list(ngrid.r=12, ngrid.c=4, nspot.r=25, nspot.c=26); > RG <- kooperberg(names=targets$FileName,layout=layout); > RG$genes <- readGAL(); > MA <- normalizeWithinArrays(RG, layout, method="robustspline"); Warning message: rlm failed to converge in 20 steps in: rlm.default(x, y, weights = w, method = method, maxit = maxit) To get around this I modified normalizeRobustSpline and normalizeWithinArrays to accept a maxit parameter to pass to rlm (the code is below). This allows: > MA <- normalizeWithinArrays(RG, layout, method="robustspline", maxit=20); Warning message: rlm failed to converge in 20 steps in: rlm.default(x, y, weights = w, method = method, maxit = maxit) > MA <- normalizeWithinArrays(RG, layout, method="robustspline", maxit=21); Warning message: rlm failed to converge in 21 steps in: rlm.default(x, y, weights = w, method = method, maxit = maxit) > MA <- normalizeWithinArrays(RG, layout, method="robustspline", maxit=22); > MA <- normalizeWithinArrays(RG, layout, method="robustspline", maxit=23); Which indicated that 22 iterations were necessary for convergence. To test that this didn't break anything else I also: > MA.old <- normalizeWithinArrays(RG, layout, method="printtiploess"); > # changed functions here: removed for email-display # > MA.new <- normalizeWithinArrays(RG, layout, method="printtiploess", maxit=25); > write.table(MA.old$M, 'old.txt'); > write.table(MA.new$M, 'new.txt'); And then did a file-diff on the files: C:\Program Files\R\rw1091>comp Name of first file to compare: old.txt Name of second file to compare: new.txt Option: Comparing old.txt and new.txt... Files compare OK The code is below -- hopefully it's useful to other people as well. Now, my question is: should the fact that my data took more than 20 iterations to converge be sounding alarm-bells in my head? Paul ### BEGIN CODE ### normalizeWithinArrays <- function (object, layout = object$printer, method = "printtiploess", weights = object$weights, span = 0.3, iterations = 4, controlspots = NULL, df = 5, robust = "M", maxit=20) { if (!is(object, "MAList")) object <- MA.RG(object) choices <- c("none", "median", "loess", "printtiploess", "composite", "robustspline") method <- match.arg(method, choices) if (method == "none") return(object) narrays <- ncol(object$M) if (method == "median") { for (j in 1:narrays) object$M[, j] <- object$M[, j] - median(object$M[, j], na.rm = TRUE) return(object) } switch(method, loess = { for (j in 1:narrays) { y <- object$M[, j] x <- object$A[, j] w <- weights[, j] object$M[, j] <- loessFit(y, x, w, span = span, iterations = iterations)$residuals } }, printtiploess = { if (is.null(layout)) stop("Layout argument not specified") ngr <- layout$ngrid.r ngc <- layout$ngrid.c nspots <- layout$nspot.r * layout$nspot.c for (j in 1:narrays) { spots <- 1:nspots for (gridr in 1:ngr) for (gridc in 1:ngc) { y <- object$M[spots, j] x <- object$A[spots, j] w <- weights[spots, j] object$M[spots, j] <- loessFit(y, x, w, span = span, iterations = iterations)$residuals spots <- spots + nspots } } }, composite = { if (is.null(layout)) stop("Layout argument not specified") if (is.null(controlspots)) stop("controlspots argument not specified") ntips <- layout$ngrid.r * layout$ngrid.c nspots <- layout$nspot.r * layout$nspot.c for (j in 1:narrays) { y <- object$M[, j] x <- object$A[, j] w <- weights[, j] f <- is.finite(y) & is.finite(x) & is.finite(w) y[!f] <- NA fit <- loess(y ~ x, weights = w, span = span, subset = controlspots, na.action = na.exclude, degree = 0, surface = "direct", family = "symmetric", trace.hat = "approximate", iterations = iterations) alpha <- global <- y global[f] <- predict(fit, newdata = x[f]) alpha[f] <- (rank(x[f]) - 1)/sum(f) spots <- 1:nspots for (tip in 1:ntips) { y <- object$M[spots, j] x <- object$A[spots, j] w <- weights[spots, j] local <- loessFit(y, x, w, span = span, iterations = iterations)$fitted object$M[spots, j] <- object$M[spots, j] - alpha[spots] * global[spots] - (1 - alpha[spots]) * local spots <- spots + nspots } } }, robustspline = { if (is.null(layout)) stop("Layout argument not specified") for (j in 1:narrays) object$M[, j] <- normalizeRobustSpline(object$M[, j], object$A[, j], layout, df = df, method = robust, maxit=maxit) }) object } normalizeRobustSpline <- function (M, A, layout, df = 5, method = "M", maxit=20) { require(MASS) require(splines) ngrids <- layout$ngrid.r * layout$ngrid.c nspots <- layout$nspot.r * layout$nspot.c O <- is.finite(M) & is.finite(A) X <- matrix(NA, ngrids * nspots, df) X[O, ] <- ns(A[O], df = df, intercept = TRUE) x <- X[O, , drop = FALSE] y <- M[O] w <- rep(1, length(y)) s <- summary(rlm(x, y, weights = w, method = method, maxit=maxit), method = "XtWX", correlation = FALSE) beta0 <- s$coefficients[, 1] covbeta0 <- s$cov * s$stddev^2 beta <- array(1, c(ngrids, 1)) %*% array(beta0, c(1, df)) covbeta <- array(0, c(ngrids, df, df)) spots <- 1:nspots for (i in 1:ngrids) { o <- O[spots] y <- M[spots][o] if (length(y)) { x <- X[spots, ][o, , drop = FALSE] r <- qr(x)$rank if (r < df) x <- x[, 1:r, drop = FALSE] w <- rep(1, length(y)) s <- summary(rlm(x, y, weights = w, method = method, maxit=maxit), method = "XtWX", correlation = FALSE) beta[i, 1:r] <- s$coefficients[, 1] covbeta[i, 1:r, 1:r] <- s$cov * s$stddev^2 } spots <- spots + nspots } res.cov <- cov(beta) - apply(covbeta, c(2, 3), mean) Sigma0 <- covbeta0 * max(0, sum(diag(res.cov))/sum(diag(covbeta0))) spots <- 1:nspots for (i in 1:ngrids) { beta[i, ] <- beta0 + Sigma0 %*% solve(Sigma0 + covbeta[i, , ], beta[i, ] - beta0) o <- O[spots] x <- X[spots, ][o, ] M[spots][o] <- M[spots][o] - x %*% beta[i, ] M[spots][!o] <- NA spots <- spots + nspots } M } ### END CODE ###
Normalization Normalization • 2.7k views
ADD COMMENT

Login before adding your answer.

Traffic: 618 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6