Difference between revisions of "FirstPrinciples.R"
m |
(rv) |
||
Line 236: | Line 236: | ||
print(t.test(exprs(eset)[i,(1:n)+preps],exprs(eset)[i,(1:n)] , paired=FALSE, var.equal=TRUE)$statistic) | print(t.test(exprs(eset)[i,(1:n)+preps],exprs(eset)[i,(1:n)] , paired=FALSE, var.equal=TRUE)$statistic) | ||
} | } | ||
+ | |||
+ | # ----------------------- Subsetting out sigma component example ------------------------ # | ||
+ | |||
+ | MA <- new("MAList") | ||
+ | MA$M <- matrix(rnorm(100),nc=4) | ||
+ | |||
+ | fit <- lmFit(MA) | ||
+ | fit <- eBayes(fit) | ||
+ | topTable(fit, number=1) | ||
+ | |||
+ | # Select first gene | ||
+ | as.numeric(rownames(topTable(fit))[1]) | ||
+ | fit$sigma[as.numeric(rownames(topTable(fit))[1])] <- NA | ||
+ | |||
+ | topTable(fit, number=1) | ||
+ | |||
+ | # First gene now different | ||
+ | fit <- eBayes(fit) | ||
+ | topTable(fit, number=1) | ||
+ | |||
+ | # Check | ||
+ | topTable(fit,number=25) |
Revision as of 22:02, 10 December 2006
library(limma) options(digits=2) packageDescription("limma", field="Version")
- ========================================================================
- 0) cDNA example analysing M matrix in MAList
- ========================================================================
set.seed(2) n <- 4 preps <- 4
- ========================================================================
- 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="") MA$M
design <- rep(1, preps) fit <- lmFit(MA, design=design, weights=NULL) fit$coeff fit <- eBayes(fit) topTable(fit, adjust="none")
- Checking M values
rowMeans(MA$M)
- Ordinary t-statistics use unequal variance
ordinary.t <- fit$coef / fit$stdev.unscaled / fit$sigma ordinary.t
- Checking t-statistics (not paired because we are using M values)
for(i in 1:n) { # var.equal T/F is the same
print(t.test(MA$M[i,], paired=FALSE, var.equal=TRUE)$statistic)
}
- Sigma component - sd(x)
fit$sigma apply(MA$M, 1, sd)
- Stderr component - se(x_bar)
drop(fit$stdev.unscaled) * fit$sigma apply(MA$M, 1, function(x){sd(x)/sqrt(length(x))})
- ========================================================================
- 2) Underlying lm.series fits a linear regression (weighted)
- ========================================================================
x <- as.matrix(design) Y <- t(MA$M)
- lm.series, no weights usage lm.fit(x,Y)
lmfit <- lm.fit(x,Y) lmfit$coef
- Same as rowmeans if no weights vector provided
rowMeans(MA$M)
- lmFit coefficients
c(fit$coef)
- gene 1, lm.fit(x,y)
lmfit <- lm.fit(x,Y) lmfit$coef
- first principles
solve(t(x)%*%x)%*%t(x)%*%Y
- ========================================================================
- 3) Adding weights to genes
- ========================================================================
RG$weights <- matrix(1, n,preps) RG$weights[row(RG$weights)<col(RG$weights)] <- 0 RG$weights
MA <- MA.RG(RG) design <- rep(1, preps) fit <- lmFit(MA, design=design, weights=MA$weights) fit$coeff fit <- eBayes(fit) topTable(fit, adjust="none")
- Sigma component - sd(x)
MA$M[MA$weights==0] <- NA fit$sigma apply(MA$M, 1, sd, na.rm=TRUE)
- Stderr component - se(x_bar)
drop(fit$stdev.unscaled) * fit$sigma apply(MA$M, 1, function(x,...){sd(x,...)/sqrt(length(x[!is.na(x)]))}, na.rm=TRUE)
- ========================================================================
- 3) First principles using lm.wfit
- ========================================================================
- fit on first gene only
- y <- MA$M[1,]
W <- MA$weights
- weighted loop over all four genes lm.wfit(x,y,w) ( remember Y is t(MA$M) )
for(i in 1:preps) {
print(lm.wfit(x, Y[,i], W[i,])$coef)
} fit$coef
- Equivalent subsetting/averaging out 0 weight information
for(i in 1:preps) {
print(Y[,i] %*% W[i,] / sum(W[i,]))
}
- first principles on gene 2
i <- 2 w <- t(MA$weights[i,]) y <- MA$M[i,] solve(t(x)%*%t(w)%*%w%*%x)%*%t(x)%*%t(w)%*%w%*%y
- Modifying the first weight
w[1] <- 0.5 print(lm.wfit(x, y, w)$coef)
- verification by multiplying weights by y
sum(w/sum(w) * y)
- ========================================================================
- Adding blocks, duplicateCorrelation calls randomizedBlockFit (statmod)
- corfit <- duplicateCorrelation(MA, ndups=1, block=c(1,1,2,2))
- fit <- lmFit(MA, block=c(1,1,2,2), correlation=corfit$consensus)
- ========================================================================
- ========================================================================
- 3) Analysis using limma
- ========================================================================
block <- rep(1:2, each=2) corfit <- duplicateCorrelation(MA, ndups=1, block=block) fit <- lmFit(MA, block=block, correlation=corfit$consensus,
weights=RG$weights)
design <- rep(1, preps) fit$coeff fit <- eBayes(fit) topTable(fit)
fit$coef
i <- 4 ndups <- 1 narrays <- preps A <- factor(rep(1:narrays, rep(ndups, narrays))) block <- rep(1:2, each=2) y <- MA$M[i,] X <- as.matrix(design)
- under independence Z <- model.matrix(~0+factor(1:preps))
Z <- model.matrix(~0+block)
w <- MA$weights[i,] s <- randomizedBlockFit(y, X, Z, only.varcomp = TRUE,
maxit = 20)$varcomp
s2 <- randomizedBlockFit(y, X, Z, w, only.varcomp = TRUE,
maxit = 20)$varcomp
- gls.series by first principles
ub <- unique(block) nblocks <- length(ub) Z <- matrix(block, narrays, nblocks) == matrix(ub, narrays,
nblocks, byrow = TRUE)
V <- Z %*% (corfit$cor * t(Z)) diag(V) <- 1 wrs <- 1/sqrt(drop(MA$weights[i, ]))
cholV <- chol(V) y2 <- backsolve(cholV, y, transpose = TRUE) X2 <- backsolve(cholV, X, transpose = TRUE)
out <- lm.fit(X2, y2)
- Verification
fit$coef[4] out$coef
- ========================================================================
- Section from lm.series which calls lm.fit, and lm.wfit
- ========================================================================
- y <- as.vector(M[i, ])
- obs <- is.finite(y)
- if (sum(obs) > 0) {
- X <- design[obs, , drop = FALSE]
- y <- y[obs]
- if (is.null(weights))
- out <- lm.fit(X, y)
- else {
- w <- as.vector(weights[i, obs])
- out <- lm.wfit(X, y, w)
- }
- ========================================================================
fit <- lmFit(MA) names(fit)
- ebayes adds df's, ted) (moderated) lods, F, and pvalues (usually adjusted)
fit2 <- ebayes(fit) names(fit2)
- ----------------------- Affymetrix approach for t-tests ------------------------ #
library(affy) eset <- new("exprSet") eset@exprs <- cbind(log2(RG$R), log2(RG$G))
design <- cbind(rep(1, preps*2), rep(0:1, each=preps)) fit <- lmFit(eset, design=design, weights=NULL) fit$coeff fit <- eBayes(fit) topTable(fit, coef=2, adjust="none")
fit$coef # First coef is M of group 1, second coef is the mean differences between groups for(i in 1:n) {
print(t.test(eset@exprs[i,1:4], eset@exprs[i,5:8])$est)
}
ordinary.t <- fit$coef / fit$stdev.unscaled / fit$sigma ordinary.t
- Checking t-statistics (not paired because we are using M values)
for(i in 1:n) { # var.equal T/F is the same
print(t.test(exprs(eset)[i,(1:n)+preps],exprs(eset)[i,(1:n)] , paired=FALSE, var.equal=TRUE)$statistic)
}
- ----------------------- Subsetting out sigma component example ------------------------ #
MA <- new("MAList") MA$M <- matrix(rnorm(100),nc=4)
fit <- lmFit(MA) fit <- eBayes(fit) topTable(fit, number=1)
- Select first gene
as.numeric(rownames(topTable(fit))[1]) fit$sigma[as.numeric(rownames(topTable(fit))[1])] <- NA
topTable(fit, number=1)
- First gene now different
fit <- eBayes(fit) topTable(fit, number=1)
- Check
topTable(fit,number=25)