7A-Multivariate.qmd

---
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r, include=FALSE}
set.seed(42)
```

# Multivariate Statistics

Multivariate statistics help us analyze datasets that have many variables at once. Instead of looking at each variable separately, we look at how they relate to each other and to the samples in our data. In ecology, for example, this helps us understand patterns in species composition or environmental variables. In this class, we’ll explore two common multivariate methods:

-   **Ordination** (to reduce complexity and visualize patterns) and

-   **Clustering** (to group similar observations).

## Ordination:

### PCA

We'll start with **Principal Component Analysis (PCA)**. This is a method to reduce the number of variables while keeping most of the variation in the data. It's helpful for visualizing and summarizing complex datasets.The focus for the PCA is in the relationship among variables!

We'll use the classic `iris` dataset, which has flower measurements for three iris species.

```{r}
head(iris)
pairs(iris, col = iris$Species)
```

The `pairs()` plot shows the relationships between all pairs of variables, and how they differ among the species.

We run PCA using the four numeric columns and **scale the variables to give them equal weight.** The summary shows how much variance is explained by each principal component. The "cumulative proportion" shows how much of the total variation is captured as we include more components.

```{r}
pca = prcomp(iris[, 1:4], scale = T) # always set scale = T
# when data is very skewed --> better transform e.g. log

summary(pca)
# standard deviation^2 is variance!!!
# cum prop of PC2 is the variance that is visualized in a biplot
```

We can visualize the variance explained by each component. The scree plot shows the absolute variance, while the barplot shows the proportion explained by each component. These plots help us decide how many principal components to keep. Often, we look for a small number of components that explain most of the variation.

```{r}
# absolute variance of each component
plot(pca) # see row1 of the summary(pca): (sd)^2 = variance

# rel variance of each component
barplot(summary(pca)$importance[2, ], 
        ylab="proportion of variance explained") # displays % of variance explained by PCs
```

The **biplot** shows how the samples are positioned (numbers) in the reduced space (based on the first two components), and how the variables contribute to those components (red arrows).

```{r}
biplot(pca) # displays PC1 and PC2 AND rotation (vectors) of the different variables AND observations
abline(v=0, lty="dashed") # adding lines in 0,0 to see better
abline(h=0, lty="dashed")
```

### NMDS

Now let’s try **Non-metric Multidimensional Scaling (NMDS)**. This is another type of ordination that works with distance (or dissimilarity) between samples (observations) instead of raw data. It’s often used in ecology to explore patterns in species composition across sites.

The focus on the NMDS is in the observations - the distance between observations given their properties (the variables).

We'll use the `dune` dataset from the `vegan` package, which is a community dataset for plants in dunes. Each **row is a site**, and each **column is a species**. The values are species counts.

```{r}
library(vegan)
#?vegan
data("dune")
str(dune) 
#?dune
summary(dune)
```

First, let's chose the distance/dissimilarity metric. We use here a **Bray-Curtis dissimilarity matrix** (default in the `metaMDS`function), which is commonly used for community data. It measures how different the species compositions are between sites.

The `metaMDS` function runs NMDS and tries to represent the distances between sites in 2D space as accurately as possible.

```{r}
NMDS = metaMDS(dune)

NMDS # gives information on NMDS: distance measure, stress (should be low)
```

The STRESS (Standard Residuals Sum of Squares), is a measure of how the sites positions in the bidimensional configuration deviates from the original distances (from the distance matrix). I can be used as a measure of how adequate is the analysis for the dataset. As a rule of thumb, we have:

-   stress \< 0.05 great representation;

-   stress \< 0.01 good ordination

-   stress \< 0.2 reazonable ordination

-   stress \> 0.2 - be suspicious

-   stress \>= 0.3 ordination is abitrary (do not interpret it)

The NMDS plot shows how similar or different the sites are in terms of species composition. Sites that are close together have more similar species compositions.

```{r}
ordiplot(NMDS, type = "t") #"t" = text for the species
```

::: callout-caution
Why we should be careful when interpreting patterns in ordination plots? Because it may find patterns even if the relationship between variables doesn't exist:

```{r}
set.seed(123)
random = data.frame(pollution = rnorm(30),
                    temperature = rnorm(30),
                    moisture = rnorm(30),
                    tourists = rnorm(30),
                    wind = rnorm(30),
                    dogs = rnorm(30))
head(random)
pca = prcomp(random, scale = T)
biplot(pca)
summary(pca) # similar variance on all axes
```
:::



## Constrained Ordination

A constrained ordination is similar to a PCA or NMDS, but considers additional predictors for the ordination.

Here an example, using the dune dataset (see via ?dune, str(dune)) together with the dune.env environmental data

```{r}
data(dune.env)
RDA = rda(dune ~ as.numeric(Manure) + as.numeric(Moisture), 
         data = dune.env)
plot(RDA)

summary(RDA)


```

* important part at the top
* variance explained by the two variables = prop constrained = 37.09%
* how much is explained by each RDA = see importance of components prop explained
* PCs are the unconstrained axes
* species scores = coordinates of species in the plot
* site scores = coordinates of sites in the plot
* biplot scores = coordinates of environmental variable vectors

```{r}

barplot(summary(RDA)$cont$importance[2, ],  las = 2,
        col = c(rep ('red', 2), 
                rep ('black', length(summary(RDA)$cont$importance[2, ])-2)),
        ylab="proportion of variance explained") 
```

displays % of variance explained by PCs

## Cluster analysis

Cluster analysis groups similar observations. It can be useful when we want to classify sites or samples based on their characteristics (variables).

There are two main types of clustering:

-   **Hierarchical clustering**: builds a tree (dendrogram) of nested clusters.

-   **Non-hierarchical clustering** (like K-means): partitions the data into a fixed number of clusters.

| Feature | Hierarchical | K-means (Non-hierarchical) |
|-----------------------|-------------------------|-------------------------|
| Number of groups | Decided after tree is built | Must be chosen before running |
| Input/ data type | Distance matrix | Raw (standardized) data |
| Result | Dendrogram (tree) | Cluster assignments |
| Flexibility | Can inspect different levels | Fixed number of groups |
| Sensitive to scaling? | Yes | Yes |
| Good for | Exploring nested group structure | Partitioning into defined clusters |

: Comparing hiearchical and non-hiearchical clustering

Hierarchical clustering gives a full tree of nested groups, which is great for exploring data. K-means is simpler and faster, but you need to choose the number of groups first.

### Hierarchical clustering

Let's pretend we don't know which species the iris individuals belong. We will use a hierarchical clustering technique to see how the individuals will be clustered by their morphological variables. We plot the results and then color individuals by their known species to see if the clustering could provide reliable classifications

For clustering you also need a dissimilarity matrix, here we are using the default from `dist`, the euclidian distance.

```{r}
library(cluster) # clustering
# example for distance matrix
dist(iris[1:3, 1:4]) # creates a distance matrix (comparison of all possible sample pairs)
```

Clustering and plotting.

```{r}
hc = hclust(dist(iris[, 1:4]))
plot(hc)
```

THe default `plot.hclust` is ugly and don't allow to color the tips (individuals). For that, we will use the package `ape`for phylogenetic analyses to get a pretty dendogram.

```{r}
library(ape)
new.hc <- as.phylo(hc, method = "ward.D2") # converts object type "hcclust" into type "phylo"
plot(new.hc, tip.color = as.numeric(iris$Species)) 
# change plotting type:
plot(as.phylo(hc), tip.color = as.numeric(iris$Species), type = "fan")
```

The clustering algorithm did a good job in classifying the 3 species. Only 4 individuals (71,73,8,107) were wrongly attributed to a cluster formed mostly by another species. However, we see that the only species that had an isolated branch for itself was the setosa (black color), the other species had individuals spread in different braches.

````{=html}
<!---
MELINA: I hided it because it would be too much to explain
```{r, eval=F}
# lets try another clustering algorithm
data(animals)
str(animals)

# Agglomerative Nesting 
aa <- agnes(animals)
plot(aa, which.plots = 2) #which.plots: plots only plot 2
# first is banner plot...
```
--->
````

### Non-hiearchical clustering

The most common method for non-hierarchical clustering is **K-means clustering**, which partitions the data into a fixed number of groups. K-means clustering requires you to choose the number of clusters in advance. It works by assigning observations to the nearest cluster center, then updating the centers.

Note: K-means uses raw data (not a distance matrix), so we need to prepare the data differently.

```{r}
set.seed(123) # choice of first k centers is random and depends on (random) seed
cl = kmeans(iris[, 1:4], centers = 3) # centers = number of clusters

cl$cluster
table(cl$cluster) # see the cluster assignment
as.numeric(iris$Species)
```

Let's compare the clustering assignmet with the real species data:

```{r}
table(clusters = cl$cluster, real = iris$Species)
```

We see again that the `setosa` species had all its individuals correctly clustered. The species `versicolor` and `virginica` presented some overlap.

Let's see the individuals that were wrongly attributed to another species in the pca plot (now using the function `ordiplot` from the `vegan` package):

```{r}
wrong.class <- data.frame(cluster = cl$cluster, species = iris$Species)
wrong.class$wrong <- "black" # creating column to color
# color red for the wrong classification
wrong.class$wrong[wrong.class$species=="virginica" & 
                           wrong.class$cluster != 3 |
                           wrong.class$species=="versicolor" & 
                           wrong.class$cluster != 2] <- "red"

ordiplot(pca, display="sites",) |>
    points("site", col=wrong.class$wrong, pch=16)
abline(v=0, lty="dashed") # adding lines in 0,0 to see better
abline(h=0, lty="dashed")
```

This was just a brief overview of clustering methods. [Here](https://bookdown.org/brittany_davidson1993/bookdown-demo/cluster-analysis.html#cluster-analysis) is a good online source for you to learn more.