Difference between revisions of "LmFit-Weights.R"
m (rv) |
(added affy and reference design examples) |
||
| Line 135: | Line 135: | ||
drop(fit$stdev.unscaled) * fit$sigma | drop(fit$stdev.unscaled) * fit$sigma | ||
apply(MA$M * designMat, 1, function(x,...){sd(x,...)/sqrt(length(x[!is.na(x)]))}, na.rm=TRUE) | apply(MA$M * designMat, 1, function(x,...){sd(x,...)/sqrt(length(x[!is.na(x)]))}, na.rm=TRUE) | ||
| + | |||
| + | |||
| + | # ======================================================================== | ||
| + | # Two groups common reference style analysis | ||
| + | # ======================================================================== | ||
| + | |||
| + | 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") | ||
| + | # Ordering genes same ass topTable for convenience | ||
| + | #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 <- cbind(1, rep(c(0,1), each=5)) | ||
| + | colnames(design) <- c("A","B-A") | ||
| + | fit <- lmFit(MA, design=design, weights=NULL) | ||
| + | fit <- eBayes(fit) | ||
| + | topTable(fit, coef="B-A", adjust="none") | ||
| + | |||
| + | # Sigma component - sd(x) | ||
| + | fit$sigma | ||
| + | apply(MA$M , 1, sd, na.rm=TRUE) # Not quite the same (lmfit rank different in design) | ||
| + | |||
| + | plot(fit$sigma, apply(MA$M , 1, sd, na.rm=TRUE)) | ||
| + | |||
| + | #apply(MA$M * designMat, 1, sd, na.rm=TRUE) | ||
| + | |||
| + | # Stderr component - se(x_bar) | ||
| + | drop(fit$stdev.unscaled[,2]) * fit$sigma # Check differences here! | ||
| + | apply(MA$M, 1, function(x,...){sd(x,...)/sqrt(length(x[!is.na(x)]))}, na.rm=TRUE) | ||
| + | |||
| + | # ======================================================================== | ||
| + | # Two groups common reference style analysis (alternate parametrisation) | ||
| + | # ======================================================================== | ||
| + | |||
| + | design <- cbind(rep(c(1,0), each=5), rep(c(0,1), each=5)) | ||
| + | colnames(design) <- c("A","B") | ||
| + | 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 | ||
| + | |||
| + | # ======================================================================== | ||
| + | # Affy style analysis (using R, G channels here) | ||
| + | # ======================================================================== | ||
| + | library(convert) | ||
| + | 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") | ||
| + | |||
| + | exprSet <- as(RG, "exprSet") | ||
| + | |||
| + | design <- cbind(rep(1, 20), rep(c(1,0), 10)) | ||
| + | colnames(design) <- c("A","B-A") | ||
| + | |||
| + | fit <- lmFit(exprSet, design=design) | ||
| + | fit <- eBayes(fit) | ||
| + | topTable(fit, adjust="none", coef="B-A") | ||
| + | fit$df.residual | ||
| + | fit$sigma | ||
| + | |||
| + | design <- cbind(rep(c(0,1), 10), rep(c(1,0), 10)) | ||
| + | colnames(design) <- c("A","B") | ||
| + | cont.matrix <- makeContrasts("B-A" = B-A, levels=design) | ||
| + | |||
| + | fit <- lmFit(exprSet, design=design) | ||
| + | fit <- contrasts.fit(fit, cont.matrix) | ||
| + | fit <- eBayes(fit) | ||
| + | topTable(fit, adjust="none") | ||
| + | fit$df.residual | ||
| + | fit$sigma | ||
| + | |||
| + | |||
| + | # Sigma component - sd(x) | ||
| + | fit$sigma | ||
| + | apply(exprs(exprSet), 1, sd) | ||
| + | plot(fit$sigma, apply(exprs(exprSet), 1, sd)) | ||
| + | abline(0,1) | ||
| + | #apply(exprs(exprSet) * outer(rep(1, n), (design[,2]*2)-1), 1, sd) | ||
| + | #apply(exprs(exprSet) * outer(rep(1, n), design), 1, sd) | ||
| + | |||
| + | # 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))}) | ||
Revision as of 22:48, 2 August 2006
library(limma) options(digits=2) packageDescription("limma", field="Version")
- ========================================================================
- 0) cDNA example analysing M matrix in MAList
- ========================================================================
set.seed(2) n <- 50 preps <- 10
- ========================================================================
- 1) Generating a RGList for transformation into MAList
- ========================================================================
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")
- Ordering genes same ass topTable for convenience
- 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) fit <- lmFit(MA, design=design, weights=NULL) fit <- eBayes(fit) topTable(fit, adjust="none")
- Sigma component - sd(x)
fit$sigma apply(MA$M * outer(rep(1, n), design), 1, sd)
- 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))})
- ========================================================================
- 2) Adding some weights for MAList
- ========================================================================
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
- Zero weights should be identical to NA's
MA$M[MA$weights==0] <- NA designMat <- outer(rep(1, n), design)
- Sigma component - sd(x)
fit$sigma apply(MA$M * designMat, 1, sd, na.rm=TRUE)
- 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)
- ========================================================================
- 3) Added some NA's for MAList (several lines above)
- ========================================================================
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
- Sigma component - sd(x)
fit$sigma apply(MA$M * designMat, 1, sd, na.rm=TRUE)
- 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)
- ========================================================================
- 3) Added weights and NA's simultaneously at random for MAList
- ========================================================================
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
- Zero weights should be identical to NA's
MA$M[MA$weights==0] <- NA
- Sigma component - sd(x)
fit$sigma apply(MA$M * designMat, 1, sd, na.rm=TRUE)
- 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)
- ========================================================================
- Two groups common reference style analysis
- ========================================================================
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")
- Ordering genes same ass topTable for convenience
- 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 <- cbind(1, rep(c(0,1), each=5)) colnames(design) <- c("A","B-A") fit <- lmFit(MA, design=design, weights=NULL) fit <- eBayes(fit) topTable(fit, coef="B-A", adjust="none")
- Sigma component - sd(x)
fit$sigma apply(MA$M , 1, sd, na.rm=TRUE) # Not quite the same (lmfit rank different in design)
plot(fit$sigma, apply(MA$M , 1, sd, na.rm=TRUE))
- apply(MA$M * designMat, 1, sd, na.rm=TRUE)
- Stderr component - se(x_bar)
drop(fit$stdev.unscaled[,2]) * fit$sigma # Check differences here! apply(MA$M, 1, function(x,...){sd(x,...)/sqrt(length(x[!is.na(x)]))}, na.rm=TRUE)
- ========================================================================
- Two groups common reference style analysis (alternate parametrisation)
- ========================================================================
design <- cbind(rep(c(1,0), each=5), rep(c(0,1), each=5)) colnames(design) <- c("A","B") 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
- ========================================================================
- Affy style analysis (using R, G channels here)
- ========================================================================
library(convert) 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")
exprSet <- as(RG, "exprSet")
design <- cbind(rep(1, 20), rep(c(1,0), 10)) colnames(design) <- c("A","B-A")
fit <- lmFit(exprSet, design=design) fit <- eBayes(fit) topTable(fit, adjust="none", coef="B-A") fit$df.residual fit$sigma
design <- cbind(rep(c(0,1), 10), rep(c(1,0), 10)) colnames(design) <- c("A","B") cont.matrix <- makeContrasts("B-A" = B-A, levels=design)
fit <- lmFit(exprSet, design=design) fit <- contrasts.fit(fit, cont.matrix) fit <- eBayes(fit) topTable(fit, adjust="none") fit$df.residual fit$sigma
- Sigma component - sd(x)
fit$sigma apply(exprs(exprSet), 1, sd) plot(fit$sigma, apply(exprs(exprSet), 1, sd)) abline(0,1)
- apply(exprs(exprSet) * outer(rep(1, n), (design[,2]*2)-1), 1, sd)
- apply(exprs(exprSet) * outer(rep(1, n), design), 1, sd)
- 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))})
