LmFit-design.R

Code snipits and programs written in R, S or S-PLUS

1. ========================================================================
2. This script examines the effect of the design matrix sigma estimates
3. ========================================================================

library(limma) options(digits=2) packageDescription("limma", field="Version")

set.seed(2) n <- 50 preps <- 10

RG <- new("RGList") RG$R <- matrix(rnorm(n*preps,8,2), nc=preps) RG$G <- matrix(rnorm(n*preps,8,2), nc=preps) RG$printer <- structure(list(ngrid.c=1, ngrid.r=1, nspot.c=2, nspot.r=2),

                        class="printlayout")

MA <- MA.RG(RG) colnames(MA$M) <- paste("array", 1:preps, sep="") rownames(MA$M) <- paste("gene", 1:n, sep="")

1. ========================================================================
2. 1) Analysis, two groups common reference
3. ========================================================================

design <- cbind(1, rep(c(0,1), each=5)) colnames(design) <- c("A","B-A") design

fit <- lmFit(MA, design=design, weights=NULL) fit <- eBayes(fit) topTable(fit, coef="B-A", adjust="none")

1. Sigma component - sd(x). The design matrix parametrisation is fitting a
2. more complex model, decreasing the rank by 1 which gives slightly
3. modified estimates of sigma.

1. The code in lm.series which calculates fit$sigma is;
2. fit$sigma <- sqrt(colMeans(fit$effects[(fit$rank + 1):narrays,
3. , drop = FALSE]^2))

A <- apply(MA$M[,1:5], 1, sd) B <- apply(MA$M[,6:10], 1, sd) sqrt((A^2+B^2)/2) fit$sigma

apply(MA$M , 1, sd, na.rm=TRUE) # Naive sd calculation different, (design parametrisation)

1. Stderr component - se(x_bar)

A<- apply(MA$M[,1:5], 1, function(x,...){sd(x,...)/sqrt(length(x[!is.na(x)]))}, na.rm=TRUE) B<- apply(MA$M[,6:10], 1, function(x,...){sd(x,...)/sqrt(length(x[!is.na(x)]))}, na.rm=TRUE) sqrt(A^2+B^2) drop(fit$stdev.unscaled[,2]) * fit$sigma

1. ========================================================================
2. 2) Analysis, two groups common reference (alternate parametrization)
3. ========================================================================

design <- cbind(rep(c(1,0), each=5), rep(c(0,1), each=5)) colnames(design) <- c("A","B") design

cont.matrix <- makeContrasts("B-A" = B-A, levels=design)

fit2 <- lmFit(MA, design=design, weights=NULL) fit2 <- contrasts.fit(fit2, cont.matrix) fit2 <- eBayes(fit2) topTable(fit2, adjust="none")

fit$sigma # identical fit2$sigma

1. ========================================================================
2. 3) Analysis, of affy expression type set
3. ========================================================================

library(convert)

exprSet <- as(RG, "exprSet")

design <- cbind(rep(1, 20), rep(c(1,0), 10)) colnames(design) <- c("A","B-A") design

BminusA <- design[,2]==1

fit <- lmFit(exprSet, design=design) fit <- eBayes(fit) topTable(fit, adjust="none", coef="B-A") fit$df.residual

A <- apply(exprs(exprSet)[,BminusA], 1, sd) B <- apply(exprs(exprSet)[,!BminusA], 1, sd) sqrt((A^2+B^2)/2) fit$sigma

1. Stderr component - se(x_bar)

A<- apply(exprs(exprSet)[,BminusA], 1, function(x,...){sd(x,...)/sqrt(length(x[!is.na(x)]))}, na.rm=TRUE) B<- apply(exprs(exprSet)[,!BminusA], 1, function(x,...){sd(x,...)/sqrt(length(x[!is.na(x)]))}, na.rm=TRUE) sqrt(A^2+B^2) drop(fit$stdev.unscaled[,2]) * fit$sigma

1. ========================================================================
2. 4) Analysis, of affy expression type set (alternate parametrization)
3. ========================================================================

design <- cbind(rep(c(0,1), 10), rep(c(1,0), 10)) colnames(design) <- c("A","B") design

cont.matrix <- makeContrasts("B-A" = B-A, levels=design) fit2 <- lmFit(exprSet, design=design) fit2 <- contrasts.fit(fit2, cont.matrix) fit2 <- eBayes(fit2) topTable(fit2, adjust="none") fit2$df.residual fit2$sigma