Difference between revisions of "LmFit-Weights.R"

From Organic Design wiki
(Mods to first 4 examples)
m
 
(3 intermediate revisions by the same user not shown)
Line 1: Line 1:
 +
# {{R}}
 
# ========================================================================
 
# ========================================================================
 
# This script examines the effect of weights and NA's on sigma estimates  
 
# This script examines the effect of weights and NA's on sigma estimates  
Line 4: Line 5:
  
 
library(limma)
 
library(limma)
−
options(digits=2)
+
options(digits=3)
 
packageDescription("limma", field="Version")
 
packageDescription("limma", field="Version")
  
Line 121: Line 122:
 
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
+
# 4)  Added weights of 0.1 for an MAList
 
# ========================================================================
 
# ========================================================================
 +
 +
set.seed(2)
 +
n    <- 50
 +
preps <- 4
  
 
RG <- new("RGList")
 
RG <- new("RGList")
Line 132: Line 136:
 
RG$printer <-  structure(list(ngrid.c=1, ngrid.r=1, nspot.c=2, nspot.r=2),
 
RG$printer <-  structure(list(ngrid.c=1, ngrid.r=1, nspot.c=2, nspot.r=2),
 
                         class="printlayout")
 
                         class="printlayout")
−
# Ordering genes same ass topTable for convenience
 
−
#RG <- RG[c(4,2,1,3),]
 
  
 
MA <- MA.RG(RG)
 
MA <- MA.RG(RG)
Line 139: Line 141:
 
rownames(MA$M) <- paste("gene", 1:n, sep="")
 
rownames(MA$M) <- paste("gene", 1:n, sep="")
  
−
design <- cbind(1, rep(c(0,1), each=5))
+
MA$weights <- matrix(0.01, nr=n, nc=preps)
−
colnames(design) <- c("A","B-A")
+
MA$weights[,1] <- 1
−
fit <- lmFit(MA, design=design, weights=NULL)
 
−
fit <- eBayes(fit)
 
−
topTable(fit, coef="B-A", adjust="none")
 
  
 +
design <- rep(1, preps)
 +
fit <- lmFit(MA, design, weights=MA$weights)
  
−
# Sigma component - sd(x). The design matrix parametrisation is fitting a more complex
+
# Weighting the coefficients
−
# model, decreasing the rank by 1 which gives slightly modified estimates of sigma
+
fit$coef[1]
−
# The code in lm.series which calculates fit$sigma is;
 
  
−
# fit$sigma <- sqrt(colMeans(fit$effects[(fit$rank + 1):narrays, , drop = FALSE]^2))
+
x <- MA$M[1, ,drop=FALSE]
 +
weights <- MA$weights[1,]
 +
sum(x * weights)/ sum(weights)
  
−
fit$sigma
+
# Full residual df given
−
apply(MA$M , 1, sd, na.rm=TRUE) # Not quite the same
+
fit$df.residual
  
−
plot(fit$sigma, apply(MA$M , 1, sd, na.rm=TRUE))
+
# Error calculation
 +
lmfit <- lm.fit(as.matrix(design), t(x))
 +
sqrt(sum(lmfit$res^2)/lmfit$df.residual) # unweighted linear model
  
−
# Stderr component - se(x_bar)
+
lmwfit <- lm.wfit(as.matrix(design), t(x), weights)
−
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)
 
−
 
 
−
a <- drop(fit$stdev.unscaled[,2]) * fit$sigma # Check differences here!
 
−
b <- apply(MA$M, 1, function(x,...){sd(x,...)/sqrt(length(x[!is.na(x)]))}, na.rm=TRUE)
 
−
 
 
−
plot(a, b*2)
 
−
abline(0,1)
 
−
# ========================================================================
 
−
# 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
 
fit$sigma
 +
sqrt(sum(weights * lmwfit$res^2)/lmwfit$df.residual) # weighted model
  
−
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)
 
fit <- eBayes(fit)
−
topTable(fit, adjust="none")
+
topTable(fit)
−
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))})
 

Latest revision as of 00:40, 6 June 2007

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)

Retrieved from "https://organicdesign.nz/wiki/index.php?title=LmFit-Weights.R&oldid=71167"