exercise 4

load("~/scripts/day4/tumorCellLines_anova.RData")

coln <- colnames(tumorCellLines)

coln_split <- strsplit(coln, split = "\\.")

cancerType <- sapply(coln_split, function(x) x[1])
fac <- factor(cancerType)

Define some useful functions for applying

dirtyHack <- function(list_fac) {
    return(function(x) {
        return(data.frame(expr = as.numeric(x), CancerType = list_fac))
    })
}

anovaDo <- function(x) {
    aov(expr ~ CancerType, data = x)
}

Function that has saves the factors and gives a df for every row back

givesStuff <- dirtyHack(fac)

bb <- apply(tumorCellLines, 1, givesStuff)
cc <- lapply(bb, anovaDo)

givesP <- function(x) {
    summary(x)[[1]][["Pr(>F)"]][1]
}

allp <- lapply(cc, givesP)
length(allp[as.numeric(allp) < 0.05])

## [1] 3700

p adjustemtn

bon <- p.adjust(allp, method = "bonferroni")
bh <- p.adjust(allp, method = "BH")

bon.df <- as.data.frame(bon)
bh.df <- as.data.frame(bh)

bh.df$sig <- as.factor(bh.df[, 1] < 0.05)
bon.df$sig <- as.factor(bon.df[, 1] < 0.05)

dim(bh.df[(bh.df[, 1] < 0.05), ])

## [1] 3328    2

dim(bon.df[(bon.df[, 1] < 0.05), ])

## [1] 576   2

Tow way ANOVA

library(MASS)
data(crabs)
aovFL <- aov(formula = FL ~ sp + sex + sex * sp, data = crabs)
aovRW <- aov(formula = RW ~ sp + sex + sex * sp, data = crabs)
aovCL <- aov(formula = CL ~ sp + sex + sex * sp, data = crabs)
aovCW <- aov(formula = CW ~ sp + sex + sex * sp, data = crabs)
aovBD <- aov(formula = BD ~ sp + sex + sex * sp, data = crabs)

summary(aovFL)

##              Df Sum Sq Mean Sq F value  Pr(>F)    
## sp            1    466     466   48.63 4.6e-11 ***
## sex           1      5       5    0.48  0.4913    
## sp:sex        1     81      81    8.41  0.0042 ** 
## Residuals   196   1880      10                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(aovRW)

##              Df Sum Sq Mean Sq F value  Pr(>F)    
## sp            1    131   131.4    25.3 1.1e-06 ***
## sex           1    112   112.1    21.6 6.1e-06 ***
## sp:sex        1     58    58.0    11.2 0.00099 ***
## Residuals   196   1016     5.2                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(aovCL)

##              Df Sum Sq Mean Sq F value  Pr(>F)    
## sp            1    838     838   18.58 2.6e-05 ***
## sex           1    111     111    2.46   0.118    
## sp:sex        1    293     293    6.50   0.012 *  
## Residuals   196   8843      45                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(aovCW)

##              Df Sum Sq Mean Sq F value Pr(>F)   
## sp            1    576     576   10.06 0.0018 **
## sex           1     68      68    1.19 0.2762   
## sp:sex        1    455     455    7.94 0.0053 **
## Residuals   196  11232      57                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(aovBD)

##              Df Sum Sq Mean Sq F value  Pr(>F)    
## sp            1    419     419   44.31 2.8e-10 ***
## sex           1     19      19    1.99   0.160    
## sp:sex        1     42      42    4.48   0.035 *  
## Residuals   196   1854       9                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

crabs.df <- crabs[, -3]
# crabs.df$GRP <-
# as.factor(paste(as.character(crabs.df$sp),as.character(crabs.df$sex),sep=''))

Das hier ist echt interessant, Modell an den pairwise Tukey üebrgeben

### man kann sachen übergeben
TukeyHSD(aovBD)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = BD ~ sp + sex + sex * sp, data = crabs)
## 
## $sp
##      diff   lwr   upr p adj
## O-B 2.895 2.037 3.753     0
## 
## $sex
##      diff     lwr   upr  p adj
## M-F 0.613 -0.2447 1.471 0.1603
## 
## $`sp:sex`
##           diff      lwr     upr  p adj
## O:F-B:F  3.816  2.22218  5.4098 0.0000
## B:M-B:F  1.534 -0.05982  3.1278 0.0640
## O:M-B:F  3.508  1.91418  5.1018 0.0000
## B:M-O:F -2.282 -3.87582 -0.6882 0.0015
## O:M-O:F -0.308 -1.90182  1.2858 0.9588
## O:M-B:M  1.974  0.38018  3.5678 0.0084

mi <- TukeyHSD(aovBD)
# pairwise.t.test(aovBD, c('O','B'))
# pairwise.t.test(crabs.df$CL,crabs.df$GRP)

# t<-melt(crabs.df, id=c('sp','sex'))