# Estimating array weights for simulated data # Figure 7 # Simulate data set.seed(20081978) a <- 10 sdgene <- 1 arraygammas <- c(seq(-1,-0.2, by=0.2), seq(0.2,1, by=0.2)) design <- matrix(1, a, 1) n <- 10000 mu <- rnorm(n, sd=1) M <- rep(0, n*a) dim(M) <- c(n, a) for(i in 1:a) { M[,i] <- rnorm(n, mean=0, sd=sqrt(exp(arraygammas[i]))) } # Estimate array weights using varying numbers of genes # Gene-by-gene algorithm array.wts.10000 <- arrayWeights(M, design) array.wts.1000 <- arrayWeights(M[1:1000,], design) array.wts.100 <- arrayWeights(M[1:100,], design) # Full algorithm array.wts.full.10000 <- arrayWeights(M, design, method="reml") array.wts.full.1000 <- arrayWeights(M[1:1000,], design, method="reml") array.wts.full.100 <- arrayWeights(M[1:100,], design, method="reml") # Calculate estimated array parameters (log scale) array.gammas.10000 <- -log(array.wts.10000) array.gammas.1000 <- -log(array.wts.1000) array.gammas.100 <- -log(array.wts.100) array.gammas.full.10000 <- -log(array.wts.full.10000) array.gammas.full.1000 <- -log(array.wts.full.1000) array.gammas.full.100 <- -log(array.wts.full.100) # Root mean square error for array parameters estimates rms.100 <- sqrt(mean((array.gammas.100-arraygammas)^2)) rms.1000 <- sqrt(mean((array.gammas.1000-arraygammas)^2)) rms.10000 <- sqrt(mean((array.gammas.10000-arraygammas)^2)) rms.approx.full.100 <- sqrt(mean((array.gammas.100-array.gammas.full.100)^2)) rms.approx.full.1000 <- sqrt(mean((array.gammas.1000-array.gammas.full.1000)^2)) rms.approx.full.10000 <- sqrt(mean((array.gammas.10000-array.gammas.full.10000)^2)) round(c(rms.100, rms.1000, rms.10000),2) # [1] 0.16 0.07 0.01 round(c(rms.approx.full.100, rms.approx.full.1000, rms.approx.full.10000),2) # [1] 0.17 0.03 0.01 # Plot estimated versus actual array gammas plot(arraygammas, array.gammas.100, pch=1, col=gray(0.1), xlab=expression(paste("Actual ",gamma['array'], sep="")), ylab=expression(paste("Estimated ", hat(gamma)['array'], sep="")), asp=1, ylim=c(-1.05, 1.2), xlim=c(-1.05, 1.2)) points(arraygammas, array.gammas.1000, pch=2, col=gray(0.35)) points(arraygammas, array.gammas.10000, pch=3, col=gray(0.6)) abline(0,1, col="black", lwd=1, lty=2) legend(-1.05,1.2, legend=c("100 genes", "1,000 genes", "10,000 genes"), col=gray(c(0.1,0.35,0.6)), pch=seq(1:3))