Hierarchical Models Exercises
Exercises
Load the following data (you can install it from Bioconductor) and extract the data matrix using exprs
:
library(Biobase)
library(SpikeInSubset)
data(rma95)
y < exprs(rma95)
This dataset comes from an experiment in which RNA was obtained from the same background pool to create six replicate samples. Then RNA from 16 genes were artificially added in different quantities to each sample. These quantities (in picoMolars) and gene IDs are stored here:
pData(rma95)
These quantities where the same in the first three arrays and in the last three arrays. So we define two groups like this:
g < factor(rep(0:1,each=3))
and create an index of which rows are associated with the artificially added genes:
spike < rownames(y) %in% colnames(pData(rma95))

Only these 16 genes are diferentially expressed since the six samples differ only due to sampling (they all come from the same background pool of RNA).
Perform a ttest on each gene using the
rowttest
function.What proportion of genes with a pvalue < 0.01 (no multiple comparison correction) are not part of the artificially added (false positive)?
 Now compute the within group sample standard deviation for each gene (you can use group 1). Based on the pvalue cutoff, split the genes into true positives, false positives, true negatives and false negatives. Create a boxplot comparing the sample SDs for each group. Which of the following best describes the boxplot?
 A) The standard deviation is similar across all groups.
 B) On average, the true negatives have much larger variability.
 C) The false negatives have larger variability.
 D) The false positives have smaller standard deviation.

In the previous two questions, we observed results consistent with the fact that the random variability associated with the sample standard deviation leads to tstatistics that are large by chance.
The sample standard deviation we use in the ttest is an estimate and with just a pair of triplicate samples, the variability associated with the denominator in the ttest can be large.
The following steps perform the basic
limma
analysis. We specifycoef=2
because we are interested in the difference between groups, not the intercept. TheeBayes
step uses a hierarchical model that provides a new estimate of the gene specific standard error.library(limma) fit < lmFit(y, design=model.matrix(~ g)) colnames(coef(fit)) fit < eBayes(fit)
Here is a plot of the original, new, hierarchical models based estimate versus the sample based estimate:
sampleSD = fit$sigma posteriorSD = sqrt(fit$s2.post)
Which best describes what the hierarchical model does?
 A) Moves all the estimates of standard deviation closer to 0.12.
 B) Increases the estimates of standard deviation to increase t.
 C) Decreases the estimate of standard deviation.
 D) Decreases the effect size estimates.

Use these new estimates of standard deviation in the denominator of the ttest and compute pvalues. You can do it like this:
library(limma) fit = lmFit(y, design=model.matrix(~ g)) fit = eBayes(fit) ##second coefficient relates to diffences between group pvals = fit$p.value[,2]
What proportion of genes with a pvalue < 0.01 (no multiple comparison correction) are not part of the artificially added (false positive)?
Compare to the previous volcano plot and notice that we no longer have small pvalues for genes with small effect sizes.