LmFit-Weights.R

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

1. ========================================================================
2. 0) cDNA example analysing M matrix in MAList
3. ========================================================================

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

1. ========================================================================
2. 1) Generating a RGList for transformation into MAList
3. ========================================================================

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

                        class="printlayout")

1. Ordering genes same ass topTable for convenience
2. RG <- RG[c(4,2,1,3),]

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

design <- rep(c(1,-1), preps/2)

1. design <- cbind(c(1,0,1,0,1,0,1,0,1,0),c(0,-1,0,-1,0,-1,0,-1,0,-1))
2. colnames(design) <- c("NHC1", "NHC1db")
3. cont.matrix <- makeContrasts("NHC1 vs NHC1db" = NHC1 - NHC1db, levels=design)

fit <- lmFit(MA, design=design, weights=NULL)

1. fit2 <- contrasts.fit(fit, cont.matrix)

fit <- eBayes(fit) topTable(fit, adjust="none")

if( !is.null(dim(design)) ) {

 design <- drop(design %*% c(1,1))

}

1. Sigma component - sd(x)

fit$sigma apply(MA$M * outer(rep(1, n), design), 1, sd) apply(MA$M, 1, sd) plot(1:10)

1. Stderr component - se(x_bar)

drop(fit$stdev.unscaled) * fit$sigma apply(MA$M * outer(rep(1, n), design), 1, function(x){sd(x)/sqrt(length(x))})

1. ========================================================================
2. 2) Adding some weights for MAList
3. ========================================================================

RG$weights <- matrix(1, n,preps) RG$weights[row(RG$weights)<col(RG$weights)] <- 0 RG$weights

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

design <- rep(c(1,-1), preps/2) fit <- lmFit(MA, design=design, weights=MA$weights) fit <- eBayes(fit) topTable(fit, adjust="none") fit$df.residual

1. Zero weights should be identical to NA's

MA$M[MA$weights==0] <- NA designMat <- outer(rep(1, n), design)

1. Sigma component - sd(x)

fit$sigma apply(MA$M * designMat, 1, sd, na.rm=TRUE)

1. Stderr component - se(x_bar)

drop(fit$stdev.unscaled) * fit$sigma apply(MA$M * designMat, 1, function(x,...){sd(x,...)/sqrt(length(x[!is.na(x)]))}, na.rm=TRUE)

1. ========================================================================
2. 3) Added some NA's for MAList (several lines above)
3. ========================================================================

design <- rep(c(1,-1), preps/2) fit <- lmFit(MA, design=design, weights=MA$weights) fit <- eBayes(fit) topTable(fit, adjust="none") fit$df.residual

1. Sigma component - sd(x)

fit$sigma apply(MA$M * designMat, 1, sd, na.rm=TRUE)

1. Stderr component - se(x_bar)

drop(fit$stdev.unscaled) * fit$sigma apply(MA$M * designMat, 1, function(x,...){sd(x,...)/sqrt(length(x[!is.na(x)]))}, na.rm=TRUE)

1. ========================================================================
2. 3) Added weights and NA's simultaneously at random for MAList
3. ========================================================================

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)

ncoords <- sample(n,round(n/2)) pcoords <- sample(preps,round(n/2), replace=TRUE)

MA$M[cbind(ncoords, pcoords)] <- NA

ncoords <- sample(n,round(n/2)) pcoords <- sample(preps,round(n/2), replace=TRUE)

MA$weights <- matrix(1, nc=preps, nr=n) MA$weights[cbind(ncoords, pcoords)] <- 0

design <- rep(c(1,-1), preps/2) fit <- lmFit(MA, design=design, weights=MA$weights) fit$coeff fit <- eBayes(fit) topTable(fit, adjust="none") fit$df.residual

1. Zero weights should be identical to NA's

MA$M[MA$weights==0] <- NA

1. Sigma component - sd(x)

fit$sigma apply(MA$M * designMat, 1, sd, na.rm=TRUE)

1. Stderr component - se(x_bar)

drop(fit$stdev.unscaled) * fit$sigma apply(MA$M * designMat, 1, function(x,...){sd(x,...)/sqrt(length(x[!is.na(x)]))}, na.rm=TRUE)