---
title: "Day 5 solution overview"
---

# Distance correlations

First we get the data

```{r}
load(url("http://www.biostatistics.dk/teaching/bioinformatics/data/gwasgt.rda"))
load(url("http://www.biostatistics.dk/teaching/bioinformatics/data/gwaspt.rda"))
```

1.  We pick the 21 SNPs

    ```{r}
    pick <- 791:810
    ```

2.  Computing the Pearson correlation.

    ```{r}
    library("energy")
    p.cor <- cor(genotypes[,pick])
    p.cor
    ```
    
3.  Computing the distance correlation between SNP 802 and the neighbouring SNPs.

    ```{r}
res <-     sapply(791:810, function(i) { dcor(genotypes[,802], genotypes[,i])})
    
    
    ```
    
    We can also compute the full `dcor` matrix using `apply()` and the function `pairwise_combination_indices()` from the `MESS` package:
    
    ```{r}
    res <- apply(MESS::pairwise_combination_indices(20, self=TRUE)+790, 1, function(x) {dcor(genotypes[,x[1]], genotypes[,x[2]])})
    d.cor <- diag(20)
    d.cor[lower.tri(d.cor, diag=TRUE)] <- res
    d.cor[upper.tri(d.cor, diag=TRUE)] <- res
    ```
    
    Find the largest discrepancies (absolute scale)
    
    ```{r}
    abs(p.cor)-d.cor
    # Largest discrepancy is 
    max(abs(abs(p.cor)-d.cor))
    ```
    
    
4.  Run the analysis

    ```{r}    
    correlations <- sapply(1:NCOL(genotypes), function(i) { dcor(genotypes[,i], phenotypes$BMI) })
    cor.inf <- sapply(1:NCOL(genotypes), function(i) { dcor.ttest(genotypes[,i], phenotypes$BMI) })
    ```
    
5.  Using the `dcor.ttest()` function.

    ```{r}
    res <- apply(MESS::pairwise_combination_indices(20, self=TRUE)+790, 1, function(x) {dcor.ttest(genotypes[,x[1]], genotypes[,x[2]])$p.value})
    d.cor <- diag(20)
    d.cor[lower.tri(d.cor, diag=TRUE)] <- res
    d.cor[upper.tri(d.cor, diag=TRUE)] <- res
    d.cor
    ```
    
    
6.  To add slack to the linearity assumption


# Zero-inflated and hurdle models

1.  Loading the data

    ```{r}
    library("pscl")
    load(url("http://www.biostatistics.dk/microbiome.rda"))
    head(microbiome)
    ```
    
2.  Using a zero-inflated model, i.e., where there is an excess number of zeroes. With the zero-inflated model the zeroes are due to both the binary and poisson distribution.

We need to round the outcome, `value`, to the nearest integer, if we want to model the non-excess-zeroes as poisson counts.

    ```{r}
    z1 <- zeroinfl(round(value) ~ condition + gender + age, subset=(otu==1), data=microbiome)
    summary(z1)
    ```
    
3.  Using the negative binomial

    ```{r}
    z2 <- zeroinfl(round(value) ~ condition + gender + age, 
                   dist="negbin",
                   subset=(otu==1), data=microbiome)
    summary(z2)
    ```


4.  We can specify different models by separating them with a vertical bar, `|`.

    ```{r}
    z3 <- zeroinfl(round(value) ~ condition + gender + age | condition, subset=(otu==1), data=microbiome)
    summary(z3)
    ```



5.  Using the `hurdle()` function

    ```{r}
    h1 <- hurdle(round(value) ~ condition + gender + age, subset=(otu==1), data=microbiome)
    summary(h1)
    ```
    
    
# Compositional data analysis

1.  Preparing the data and loading the libraries

    ```{r}
    library("compositions")
    exprmatrix <- matrix(microbiome$value, nrow=70)
    ```
    
    We then compute the `ilr()` function.
    
    ```{r}
    trans <- ilr(exprmatrix)
    ```
    
    
    
2.  The dimension is a 70x14 matrix. This makes sense since the 15-dimensional original data was mapped to a 14-dimensional response, and these values can be analysed n by all classical multivariate analysis tools.

    ```{r}
    dim(trans)
    ```

3.  Let's try the first example. Recall that the first 35 rows are condition A and the remaining 35 are condition B. We will code that first

    ```{r}
    condition <- rep(c("A", "B"), times=c(35, 35))
    summary(lm(trans[,1] ~ condition))
    ```



# PARAFAC


```{r warning=FALSE, message=FALSE}
library("multiway")
load(url("http://www.biostatistics.dk/parafac.rda")) # Loads
```

Let's plot the first result

```{r}
persp(x[,,1], col="lightblue")
```


1.  Start by running

    ```{r}
    res <- parafac(x, nfac=2)
    summary(res)
    ```

2.  Plotting the basis funcitons for `A`:

    ```{r}
    matplot(res$A, type="l")
    ```
    
    We can see that the two tops are captured by two of the basis functions.
    
    
3.  Approximating the matrix/tensor

    ```{r}
    appx <- fitted(res)
    ```
    
    Let's look a the difference in the first matrix slice.
    
    ```{r}
    persp(x[,,1] - appx[,,1])
    max(abs(x[,,1] - appx[,,1]))
    ```
    
    The largest difference is `r round(max(abs(x[,,1] - appx[,,1])), 2)` which is not that big.