Here are few slides explaining the problems of multiple hypothesis testing and the concept of ANOVA: analysis of variance.

## 2.1. Multiple hypothesis testing

Let’s run a small test described in the lecture: generate a random dataset with 1000 “genes” for 6 “samples”. 3 samples we denote as “A”, 3 - as “B”

library(pheatmap)

ng = 1000  ## n genes
nr = 3     ## n replicates

## random matrix
X = matrix(rnorm(ng*nr*2),nrow=ng, ncol=nr*2)
rownames(X) = paste0("gene",1:ng)
colnames(X) = c(paste0("A",1:nr),paste0("B",1:nr))
pheatmap(X,fontsize_row = 1,cellwidth=20,main="Random expression matrix",scale="row")

## calculate p-values
pv = double(ng)
names(pv) = rownames(X)
for (i in 1:ng){
pv[i] = t.test(X[i,1:nr],X[i,(nr+1):(2*nr)])$p.value } pheatmap(X[pv<0.05,],cellwidth=20,main="Significant genes by p-value",scale="row") ## Significant genes by p-value: print(names(pv)[pv<0.05]) ## [1] "gene6" "gene37" "gene115" "gene144" "gene198" "gene255" "gene258" ## [8] "gene279" "gene305" "gene361" "gene427" "gene470" "gene478" "gene490" ## [15] "gene508" "gene530" "gene540" "gene552" "gene571" "gene612" "gene623" ## [22] "gene627" "gene636" "gene651" "gene685" "gene700" "gene715" "gene726" ## [29] "gene727" "gene750" "gene823" "gene849" "gene850" "gene882" "gene886" ## [36] "gene934" "gene935" "gene946" "gene975" "gene985" ## Correction for multiple testing using Benjamini-Hochberg's FDR and Holm's FWER fdr = p.adjust(pv,"fdr") fwer = p.adjust(pv,"holm") ## Significant genes by FDR: print(names(fdr)[fdr<0.05]) ## character(0) data.frame(PV=pv[pv<0.05],FDR=fdr[pv<0.05],FWER=fwer[pv<0.05]) ## PV FDR FWER ## gene6 0.023475903 0.9767980 1.000000 ## gene37 0.039105373 0.9767980 1.000000 ## gene115 0.011645949 0.8318535 1.000000 ## gene144 0.010088350 0.8318535 1.000000 ## gene198 0.037773811 0.9767980 1.000000 ## gene255 0.005340796 0.8318535 1.000000 ## gene258 0.016426336 0.9767980 1.000000 ## gene279 0.000234788 0.2347880 0.234788 ## gene305 0.029846842 0.9767980 1.000000 ## gene361 0.006277142 0.8318535 1.000000 ## gene427 0.032706265 0.9767980 1.000000 ## gene470 0.022784497 0.9767980 1.000000 ## gene478 0.041778884 0.9767980 1.000000 ## gene490 0.006138305 0.8318535 1.000000 ## gene508 0.029061558 0.9767980 1.000000 ## gene530 0.007618467 0.8318535 1.000000 ## gene540 0.045535434 0.9767980 1.000000 ## gene552 0.014186959 0.9457973 1.000000 ## gene571 0.002258944 0.8318535 1.000000 ## gene612 0.011405607 0.8318535 1.000000 ## gene623 0.040115231 0.9767980 1.000000 ## gene627 0.045401580 0.9767980 1.000000 ## gene636 0.044793932 0.9767980 1.000000 ## gene651 0.024583983 0.9767980 1.000000 ## gene685 0.008105792 0.8318535 1.000000 ## gene700 0.035901391 0.9767980 1.000000 ## gene715 0.003663151 0.8318535 1.000000 ## gene726 0.037740048 0.9767980 1.000000 ## gene727 0.038503746 0.9767980 1.000000 ## gene750 0.035975886 0.9767980 1.000000 ## gene823 0.011430813 0.8318535 1.000000 ## gene849 0.009887362 0.8318535 1.000000 ## gene850 0.048993315 0.9767980 1.000000 ## gene882 0.047366697 0.9767980 1.000000 ## gene886 0.039158580 0.9767980 1.000000 ## gene934 0.005424578 0.8318535 1.000000 ## gene935 0.030403728 0.9767980 1.000000 ## gene946 0.028513979 0.9767980 1.000000 ## gene975 0.024840923 0.9767980 1.000000 ## gene985 0.029278682 0.9767980 1.000000 ## 2.2. Analysis of Variance (ANOVA) ## load data Dep = read.table("http://edu.modas.lu/data/txt/depression2.txt",header=T,sep="\t",as.is=FALSE) str(Dep) ## 'data.frame': 120 obs. of 3 variables: ##$ Location  : Factor w/ 3 levels "Florida","New York",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $Health : Factor w/ 2 levels "bad","good": 2 2 2 2 2 2 2 2 2 2 ... ##$ Depression: int  3 7 7 3 8 8 8 5 5 2 ...

### 2.2.1. One-way (one-factor) ANOVA

We start with one-way ANOVA. Let’s consider only healthy respondents in Dep dataset. ANOVA is performed in R: aov(formula, data)

• formula - “Values ~ Factor1”, where “Values”, “Factor1” should be replaced by the actual column names in the data frame data

ANOVA tells whether a factor is significant, but tells nothing with treatments (factor levels) are different! Use TukeyHSD() to perform post-hoc analysis and get corrected p-values for the contrasts.

DepGH = Dep[Dep$Health == "good",] ## build 1-factor ANOVA model res1 = aov(Depression ~ Location, DepGH) summary(res1) ## Df Sum Sq Mean Sq F value Pr(>F) ## Location 2 78.5 39.27 6.773 0.0023 ** ## Residuals 57 330.4 5.80 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 TukeyHSD(res1) ## Tukey multiple comparisons of means ## 95% family-wise confidence level ## ## Fit: aov(formula = Depression ~ Location, data = DepGH) ## ##$Location
##                         diff        lwr       upr     p adj
## New York-Florida         2.8  0.9677423 4.6322577 0.0014973
## North Carolina-Florida   1.5 -0.3322577 3.3322577 0.1289579
## North Carolina-New York -1.3 -3.1322577 0.5322577 0.2112612

### 2.2.2. Two-way (two-factor) ANOVA

For two-way ANOVA the call is similar aov(formula, data). Please note that factor columns in the data frame data should be of class "factor" (newer R version accepts "character"). If factors are "numeric", aov will build regression instead.

• formula - "Values ~ Factor1 + Factor2 + Factor1*Factor2“, where”Values“,”Factor1" and “Factor2” should be replaced by the actual column names in the data frame data

• Interaction "Factor1*Factor2" can be estimated in case, if you table contains replicates (values measured for the same conditions)

If needed, ANOVA can be followed by TukeyHSD().

## build 2-factor ANOVA model (with interaction)
res2 = aov(Depression ~  Location + Health + Location*Health, Dep)
summary(res2)
##                  Df Sum Sq Mean Sq F value Pr(>F)
## Location          2   73.8    36.9   4.290  0.016 *
## Health            1 1748.0  1748.0 203.094 <2e-16 ***
## Location:Health   2   26.1    13.1   1.517  0.224
## Residuals       114  981.2     8.6
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(res2)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
##
## Fit: aov(formula = Depression ~ Location + Health + Location * Health, data = Dep)
##
## $Location ## diff lwr upr p adj ## New York-Florida 1.850 0.2921599 3.4078401 0.0155179 ## North Carolina-Florida 0.475 -1.0828401 2.0328401 0.7497611 ## North Carolina-New York -1.375 -2.9328401 0.1828401 0.0951631 ## ##$Health
##               diff       lwr       upr p adj
## good-bad -7.633333 -8.694414 -6.572252     0
##
## \$Location:Health
##                                         diff         lwr       upr     p adj
## New York:bad-Florida:bad                0.90  -1.7893115  3.589311 0.9264595
## North Carolina:bad-Florida:bad         -0.55  -3.2393115  2.139311 0.9913348
## Florida:good-Florida:bad               -8.95 -11.6393115 -6.260689 0.0000000
## New York:good-Florida:bad              -6.15  -8.8393115 -3.460689 0.0000000
## North Carolina:good-Florida:bad        -7.45 -10.1393115 -4.760689 0.0000000
## North Carolina:bad-New York:bad        -1.45  -4.1393115  1.239311 0.6244461
## Florida:good-New York:bad              -9.85 -12.5393115 -7.160689 0.0000000
## New York:good-New York:bad             -7.05  -9.7393115 -4.360689 0.0000000
## North Carolina:good-New York:bad       -8.35 -11.0393115 -5.660689 0.0000000
## Florida:good-North Carolina:bad        -8.40 -11.0893115 -5.710689 0.0000000
## New York:good-North Carolina:bad       -5.60  -8.2893115 -2.910689 0.0000003
## North Carolina:good-North Carolina:bad -6.90  -9.5893115 -4.210689 0.0000000
## New York:good-Florida:good              2.80   0.1106885  5.489311 0.0361494
## North Carolina:good-Florida:good        1.50  -1.1893115  4.189311 0.5892328
## North Carolina:good-New York:good      -1.30  -3.9893115  1.389311 0.7262066

TukeyHSD() output may be quite long. If you wish saving it to a text file for later analysis - use sink() function

## start writing logs to the file
sink("tukey.txt")
## print output
TukeyHSD(res2)
## close the file and write it to the current folder
sink() 

## Exercise 2

1. Dataset teeth2 contains the result of an experiment, conducted to measure the effect of various doses of vitamin C on the tooth growth (model animal – Guinea pigs). Vitamin C and orange juice were given to animals in different quantities. Using 2-way ANOVA study the effects of vitamin-C and orange juice. Transform Dose to a factor.

factor(), aov(), summary()

1. Each of six garter snakes (Thamnophis radix) was observed during the presentation of petri dishes containing solutions of different chemical stimuli. The number of tongue flicks during a 5-minute interval of exposure was recorded. The three petri dishes were presented to each snake in random order. See data in “snakes” table. Does solution affect tongue flicks? Should we correct the results for effect of difference b/w snakes (1- or 2-factor ANOVA)? Perform post hoc analysis for the significant factor. Use snakes dataset. It should be transformed before being used - see melt of reshape2 package.

str(), melt(Data,id=c("Snake")), aov, summary

 By Petr Nazarov