---
title: "Day 1 solution overview"
---

# Multiple testing

1.  Out of 21500 test we'd expect around `0.05*21500` to be significant by chance.

4.  (and 3, 5 and 6)

```{r warning=FALSE, message=FALSE}
library("MESS")
data(superroot2)

pval <- by(superroot2,
           superroot2$gene,
           FUN=function(dat) {
              anova(lm(log(signal) ~ array + color + plant,
                       data=dat))[3,5]} )

library("ggplot2")
library("tidyverse")
DF <- data.frame(pval=as.vector(pval))

DF %>% ggplot(aes(x=pval)) + geom_histogram(col="lightblue")

# Remember that some are missing

p1 <- p.adjust(as.vector(pval), method="bonferroni")
p2 <- p.adjust(as.vector(pval), method="holm")
p3 <- p.adjust(as.vector(pval), method="BH") # FDR

head(sort(pval))
head(sort(p1))
head(sort(p2))
head(sort(p3))
```

9.  Note that one of the tests resulted in a `NA` so there are actually only 21499 tests undertaken. If this number is used as the denominator then the numbers match.

10. Due to the nature of the BH approach. If an earlier p-value had a lower adjusted fdr-value then that value is used for lower ranked adjusted fdr values.

# Penalized regression

```{r warning=FALSE, message=FALSE}
library("glmnet")
load(url("http://www.biostatistics.dk/teaching/advtopicsA/data/lassodata.rda"))

m1 <- glmnet(genotype, phenotype)

```

3.  Ensure that the penalty is the same across predictors such that it does not depend on scale.

4.

```{r}
m2 <- cv.glmnet(genotype, phenotype)
plot(m2)
lambda <- m2$lambda.1se
```

5.  

```{r}
ests <- coef(m1, s=lambda)
ests[ests != 0]
```


6.  

```{r}
m3 <- glmnet(genotype, phenotype, family="binomial")
m4 <- cv.glmnet(genotype, phenotype, family="binomial")
plot(m3)
plot(m4)
ests2 <- coef(m3, s=m4$lambda.1se)
ests2[ests2 != 0]
```



7.  

```{r}
m5 <- glmnet(genotype, phenotype, alpha=0, family="binomial")
m6 <- cv.glmnet(genotype, phenotype, alpha=0, family="binomial")
plot(m5)
plot(m6)
ests2 <- coef(m5, s=m6$lambda.1se)  # Gets them all since none are zero
```


9.  

```{r}
m7 <- glmnet(genotype, phenotype, alpha=.67, family="binomial")
m8 <- cv.glmnet(genotype, phenotype, alpha=.67, family="binomial")
plot(m7)
plot(m8)
ests2 <- coef(m7, s=m8$lambda.1se)
ests2[ests2 != 0]
```


