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 ###