# R/beadCountWeights.R In limma: Linear Models for Microarray Data

```beadCountWeights <- function(y, x, design=NULL, bead.stdev=NULL, bead.stderr=NULL, nbeads=NULL, array.cv=TRUE, scale=FALSE)
#	Compute weights for BeadChips based on bead-level counts and standard deviations
#	Charity Law and Gordon Smyth
{
#	Matrix of normalized expression values
E <- as.matrix(y)
P <- nrow(E)
A <- ncol(E)

#	Matrix of non-normalized expression values on unlog scale
E.raw <- as.matrix(x)
if(nrow(E.raw) != P) stop("dimensions don't match")
if(ncol(E.raw) != A) stop("dimensions don't match")

#	Check design matrix
if(is.null(design)) design <- y\$design
if(is.null(design))
design <- matrix(1,A,1)
else {
design <- as.matrix(design)
if(mode(design) != "numeric") stop("design must be a numeric matrix")
if(nrow(design) != A) stop("row dimension of design doesn't match column dimension of data object")
}

#	Number of beads per probe
}
if(nrow(nbeads) != P) stop("dimensions don't match")
if(ncol(nbeads) != A) stop("dimensions don't match")

#	Standard deviation of beads per probe
if(is.null(y\$other\$BEAD_STDERR)) stop("Neither bead.stdev nor bead.stderr were provided (and neither could be found in the data object y). At least one is required.")
} else {
}
} else {
}
}
if(nrow(bead.stdev) != P) stop("dimensions don't match")
if(ncol(bead.stdev) != A) stop("dimensions don't match")

#	Coefficient of variation of bead-level observations
#	Array-specific or constant.
cv.array <- apply(sqrt(cv), 2, function(k) mean(k, trim=0.125, na.rm=TRUE))^2
cv.constant <- mean(sqrt(cv), trim=0.125, na.rm=TRUE)^2
if(array.cv) {cv <- rep(cv.array, each=nrow(y))}
else cv <- cv.constant

#	Predicted variance of normalized probe-level values

# 	Squared-residuals to calculate biological variance
qrX <- qr(design)
res <- qr.resid(qrX, t(E))
r2 <- (t(res))^2

# 	Leverages
Q <- qr.Q(qrX)
h <- rowSums(Q^2)
h <- matrix(rep(h, each=P), nrow=P, ncol=A)

out <- list()
out\$var.biological <- bv
out\$var.technical <- tv
out\$cv.array <- cv.array
out\$cv.constant <- cv.constant
out\$weights <- 1/(tv + bv)
out\$weights <- out\$weights/rowMeans(out\$weights)
if(scale) out\$weights <- out\$weights*rowMeans(1/tv)
out
}

#	Find the component of the between-array variance
#	not explained by bead-level variability
#	Gordon Smyth and Charity Law
#	28 July 2010. Modified 16 July 2012.
{
if(any(var.from.beads < 0)) stop("negative variances not allowed")
if(any(sq.residuals < 0)) stop("negative variances not allowed")
r2 <- as.matrix(sq.residuals)
h <- as.matrix(leverage)
P <- nrow(tv)
A <- ncol(tv)
if(nrow(r2) != P) stop("dimensions don't match")
if(ncol(r2) != A) stop("dimensions don't match")
if(nrow(h) != P) stop("dimensions don't match")
if(ncol(h) != A) stop("dimensions don't match")

# Newton's iteration
F <- rowMeans(r2/(2*tv^2)-(1-h)/(2*tv))
bv <- rep(0,length=P)
i <- (F > 0)
iter <- 0
while(any(i)) {
iter <- iter+1
if(iter > 200) {
warning("More than 200 iterations required")
return(bv)
}
Fdash <- - rowMeans(r2[i,,drop=FALSE]/(bv[i]+tv[i,,drop=FALSE])^3 - (1-h[i,,drop=FALSE])/(2*(bv[i]+tv[i,,drop=FALSE])^2))
step <- - F[i]/Fdash
bv[i] <- bv[i] + step
F[i] <- rowMeans(r2[i,,drop=FALSE]/(2*(bv[i]+tv[i,,drop=FALSE])^2) - (1-h[i,,drop=FALSE])/(2*(bv[i]+tv[i,,drop=FALSE])))
i[i] <- (step > 1e-5)
print(sum(i))
}
bv
}
```

## Try the limma package in your browser

Any scripts or data that you put into this service are public.

limma documentation built on Nov. 8, 2020, 8:28 p.m.