---
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/modernstatistics/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]
```



# Principal component analysis


```{r}
library("MASS")
data(biopsy)
names(biopsy)
predictors <- biopsy[complete.cases(biopsy),2:10]
fit <- prcomp(predictors, scale=TRUE)
summary(fit)
fit
```

1.  The `rotation` element is the interesting part. 

2.  (and 3 + 4)

```{r}
plot(fit)
biplot(fit)
```

Only one real PC needed (from the elbow plot).  The biplot shows that one variable drives PC2 while all of them have essentially the same effect on PC1. Not surprising as they are all more and more critical indicators of tumor features.


5.  Fitting using a formula model

```{r}
fit2 <- prcomp(~ V1 + V2 + V3 + V3 + V5 + V6 + V7 + V8 + V9, 
              data=biopsy, na.action=na.omit)
print(fit2)
```


6.  

```{r}
cor(predictors)
```




# Principal component regression


We already computed the PCs above so we can just use those. The computed values can either be extracted as the `x` element or by using the `predict()` function.

1.  The outcome is binary so we will use a logistic regression model. Trying with the first four PCs.

    ```{r}
    y <- biopsy$class[complete.cases(biopsy)]
    pcr <- glm(y ~ fit$x[,1] + fit$x[,2] + fit$x[,3] + fit$x[,4],
               family=binomial)
    summary(pcr)
    ```
    

2.  PC 1 provides a lot of info about the outcome, and the coefficient is negative which means that at increase for the PC value will lower the odds of a malignant outcome. However, recall that the rotation matrix for the first principal component only had negative weights so when that is combined we get that an increase in any predictor will result in an increase in log odds of a malignant tumor.  
    We also see that even though the PCs are ranked from highest variation to lowest variation their p-values are not necessarily ranked in the same order.
    
3.  Primarily prediction as the components are difficult to interpret.

4.  
    ```{r}
    pcr2 <- glm(y ~ fit$x[,1] + fit$x[,2],
               family=binomial)
    summary(pcr2)
    ```
    
    We see that there is hardly any difference for the predictors. The reason why the predictors are not completely identical to the fit above where we had 4 predictors is because of the non-identity link function.
    
5.  We will need to see how the predictability improves with number of components.     