LmFit-Weights.R

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

1. ========================================================================
2. This script examines the effect of weights and NA's on sigma estimates
3. ========================================================================

library(limma) options(digits=3) 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, no weights, no NA's
3. ========================================================================

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

1. Sigma component - sd(x)

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

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) Analysis, adding some weights to MAList
3. ========================================================================

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

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) Analysis, adding NA's only to MAList
3. ========================================================================

MA <- MA.RG(RG) MA$M[row(MA$M)<col(MA$M)] <- NA MA$M

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. 4) Added weights and NA's simultaneously at random for MAList
3. ========================================================================

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)

1. ========================================================================
2. 4) Added weights of 0.1 for an MAList
3. ========================================================================

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

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="")

MA$weights <- matrix(0.01, nr=n, nc=preps) MA$weights[,1] <- 1

design <- rep(1, preps) fit <- lmFit(MA, design, weights=MA$weights)

1. Weighting the coefficients

fit$coef[1]

x <- MA$M[1, ,drop=FALSE] weights <- MA$weights[1,] sum(x * weights)/ sum(weights)

1. Full residual df given

fit$df.residual

1. Error calculation

lmfit <- lm.fit(as.matrix(design), t(x)) sqrt(sum(lmfit$res^2)/lmfit$df.residual) # unweighted linear model

lmwfit <- lm.wfit(as.matrix(design), t(x), weights) fit$sigma sqrt(sum(weights * lmwfit$res^2)/lmwfit$df.residual) # weighted model

fit <- eBayes(fit) topTable(fit)