```{r echo=FALSE}
# Set up chunk for all slides
knitr::opts_chunk$set(
  fig.width = 4,
  fig.height = 4,
  fig.align = "center",
  out.width = "60%",
  code.line.numbers = FALSE,
  fig.retina = 3,
  echo = TRUE,
  message = FALSE,
  warning = FALSE,
  cache = FALSE,
  dev.args = list(pointsize = 11)
)
```

```{r}
#| echo: true
#| code-fold: true
#| code-summary: "Load the libraries and avoid conflicts"
# Load libraries used everywhere
library(tidyverse)
library(tidymodels)
library(conflicted)
library(colorspace)
library(patchwork)
library(MASS)
library(randomForest)
library(gridExtra)
library(GGally)
library(geozoo)
library(mulgar)
conflicts_prefer(dplyr::filter)
conflicts_prefer(dplyr::select)
conflicts_prefer(dplyr::slice)
conflicts_prefer(palmerpenguins::penguins)
conflicts_prefer(tourr::flea)
```

```{r}
#| echo: false
# Set plot theme
theme_set(theme_bw(base_size = 14) +
   theme(
     aspect.ratio = 1,
     plot.background = element_rect(fill = 'transparent', colour = NA),
     plot.title.position = "plot",
     plot.title = element_text(size = 24),
     panel.background = element_rect(fill = 'transparent', colour = NA),
     legend.background = element_rect(fill = 'transparent', colour = NA),
     legend.key = element_rect(fill = 'transparent', colour = NA)
   )
)
```

## `r emo::ji("target")` Objectives

The goal for this week is for you to learn and practice visualising high-dimensional data. 

## `r emo::ji("wrench")` Preparation 

- Complete the quiz
- Do the reading related to week 2


## Exercises: 

Open your project for this unit called `iml.Rproj`. 

#### 1. The sparseness of high dimensions

Randomly generate data points that are uniformly distributed in a hyper-cube of 3, 5 and 10 dimensions, with 500 points in each sample, using the `cube.solid.random` function of the `geozoo` package. What differences do we expect to see? Now visualise each set in a grand tour and describe how they differ, and whether this matched your expectations? 

The code to generate and view the cubes is:

```{r eval=FALSE}
#| code-fold: true
#| code-summary: "Code to generate the data and show in a tour"
library(tourr)
library(geozoo)
set.seed(1234)
cube3 <- cube.solid.random(3, 500)$points
cube5 <- cube.solid.random(5, 500)$points
cube10 <- cube.solid.random(10, 500)$points

animate_xy(cube3, axes="bottomleft")
animate_xy(cube5, axes="bottomleft")
animate_xy(cube10, axes="bottomleft")
```


::: unilur-solution

Each of the projections has a boxy shape, which gets less distinct as the dimension increases. 

As the dimension increases, the points tend to concentrate in the centre of the plot window, with a smattering of points in the edges. 
:::

#### 2. Detecting clusters

For the data sets, `c1`, `c3` from the `mulgar` package, use the grand tour to view and try to identify structure (outliers, clusters, non-linear relationships). 

```{r eval=FALSE}
#| code-fold: true
#| code-summary: "Code to show in a tour"
animate_xy(c1)
animate_xy(c3)
```

::: unilur-solution

The first data set `c1` has 6 clusters, 4 small ones, and two big ones. The two big ones look like planes because they have no variation in some dimensions.

The second data set `c3` has a triangular prism shape, which itself is divided into several smaller triangular prisms. It also has several dimensions with no variation, because the points collapse into a line in some projections.

:::

#### 3. Effect of covariance

Examine 5D multivariate normal samples drawn from populations with a range of variance-covariance matrices.  (You can use the `mvtnorm` package to do the sampling, for example.) Examine the data using a grand tour. What changes when you change the correlation from close to zero to close to 1?  Can you see a difference between strong positive correlation and strong negative correlation?

```{r}
#| code-fold: true
#| code-summary: "Code to generate the samples"
library(mvtnorm)
set.seed(501)

s1 <- diag(5)
s2 <- diag(5)
s2[3,4] <- 0.7
s2[4,3] <- 0.7
s3 <- s2
s3[1,2] <- -0.7
s3[2,1] <- -0.7

s1
s2
s3

set.seed(1234)
d1 <- as.data.frame(rmvnorm(500, sigma = s1))
d2 <- as.data.frame(rmvnorm(500, sigma = s2))
d3 <- as.data.frame(rmvnorm(500, sigma = s3))
```

::: unilur-solution

```{r}
#| eval: false
animate_xy(d1)
animate_xy(d2)
animate_xy(d3)
```

The points in data `d1` are pretty spread in every projection. 
For the data `d2`, `d3` have some projections where the data is concentrated along a line. This should be seen to be when variables 3 and 4 are contributing to the projection in `d2`, and when variables 1, 2, 3, 4 contributing to the projection in `d3`.
:::

#### 4. Principal components analysis on the simulated data

🧐 For data sets `d2` and `d3` what would you expect would be the number of PCs suggested by PCA?

👨🏽‍💻👩‍💻Conduct the PCA. Report the variances (eigenvalues), and cumulative proportions of total variance, make a scree plot, and the PC coefficients.

🤯Often, the selected number of PCs are used in future work. For both `d3` and `d4`, think about the pros and cons of using 4 PCs and 3 PCs, respectively. 

::: unilur-solution

Thinking about it: In `d2` there is strong correlation between variables 3 and 4, which means probably only 4PC s would be needed. In `d3` there is strong correlation also between variables 1 and 2 which would mean that only 3 PCs would be needed.

```{r}
d2_pca <- prcomp(d2, scale=TRUE)
d2_pca

d2_pca$sdev^2/5

mulgar::ggscree(d2_pca, q=5)
```

Four PCs explain 93% of the variation. PC1 is the combination of variables 3 and 4, which captures this reduced dimension.

```{r}
d3_pca <- prcomp(d3, scale=TRUE)
d3_pca

d3_pca$sdev^2/5

mulgar::ggscree(d3_pca, q=5)
```

Three PCs explain 88% of the variation, and the last two PCs have much smaller variance than the others. PC 1 and 2 are  combinations of variables 1, 2, 3 and 4, which captures this reduced dimension, and PC 3 is primarily variable 5.

The PCs are awkward combinations of the original variables. For `d2`, it would make sense to use PC1 (or equivalently and equal combination of V3, V4), and then keep the original variables V1, V2, V5. 

For `d3` it's harder to make this call because the first two PCs are combinations of four variables. Its hard to see from this that the ideal solution would be to use an equal combination of V1, V2, and equal combination of V3, V4 and V5 on its own. 

Often understanding the variance that is explained by the PCs is hard to interpret.
:::

#### 5. PCA on cross-currency time series

The `rates.csv` data has 152 currencies relative to the USD for the period of Nov 1, 2019 through to Mar 31, 2020. Treating the dates as variables, conduct a PCA to examine how the cross-currencies vary, focusing on this subset: ARS, AUD, BRL, CAD, CHF, CNY, EUR, FJD, GBP, IDR, INR, ISK, JPY, KRW, KZT, MXN, MYR, NZD, QAR, RUB, SEK, SGD, UYU, ZAR. 

```{r}
#| message: false
rates <- read_csv("https://raw.githubusercontent.com/numbats/iml/master/data/rates_Nov19_Mar20.csv") |>
  select(date, ARS, AUD, BRL, CAD, CHF, CNY, EUR, FJD, GBP, IDR, INR, ISK, JPY, KRW, KZT, MXN, MYR, NZD, QAR, RUB, SEK, SGD, UYU, ZAR)
```

a. Standardise the currency columns to each have mean 0 and variance 1. Explain why this is necessary prior to doing the PCA or is it? Use this data to make a time series plot overlaying all of the cross-currencies. 
```{r}
#| label: rates-std
#| code-fold: true
#| code-summary: "Code to standardise currencies"
#| eval: false
library(plotly)
rates_std <- rates |>
  mutate_if(is.numeric, function(x) (x-mean(x))/sd(x))
rownames(rates_std) <- rates_std$date
p <- rates_std |>
  pivot_longer(cols=ARS:ZAR, 
               names_to = "currency", 
               values_to = "rate") |>
  ggplot(aes(x=date, y=rate, 
             group=currency, label=currency)) +
    geom_line() 
ggplotly(p, width=400, height=300)
```

::: unilur-solution
```{r}
#| label: rates-std
#| echo: false
```

It isn't necessary to standardise the variables before using the `prcomp` function because we can set `scale=TRUE` to have it done as part of the PCA computation. However, it is useful to standardise the variables to make the time series plot where all the currencies are drawn. This is useful for interpreting the principal components.
:::

b. Conduct a PCA. Make a scree plot, and summarise proportion of the total variance. Summarise these values and the coefficients for the first five PCs, nicely. 

```{r}
#| label: rates-pca
#| code-fold: true
#| code-summary: "Code to do PCA and screeplot"
#| eval: false
rates_pca <- prcomp(rates_std[,-1], scale=FALSE)
mulgar::ggscree(rates_pca, q=24)
options(digits=2)
summary(rates_pca)
```

```{r}
#| label: rates-smry
#| code-fold: true
#| code-summary: "Code to make a nice summary"
#| eval: false
# Summarise the coefficients nicely
rates_pca_smry <- tibble(evl=rates_pca$sdev^2) |>
  mutate(p = evl/sum(evl), 
         cum_p = cumsum(evl/sum(evl))) |> 
  t() |>
  as.data.frame()
colnames(rates_pca_smry) <- colnames(rates_pca$rotation)
rates_pca_smry <- bind_rows(as.data.frame(rates_pca$rotation),
                            rates_pca_smry)
rownames(rates_pca_smry) <- c(rownames(rates_pca$rotation),
                              "Variance", "Proportion", 
                              "Cum. prop")
rates_pca_smry[,1:5]
```

::: unilur-solution
```{r}
#| label: rates-pca
#| echo: false
```

```{r}
#| label: rates-smry
#| echo: false
```

- The first two principal components explain 85% of the total variation. 
- PC1 is a combination of all of the currencies except for CHF, EUR, JPY, QAR. 
- PC2 is a combination of CHF, EUR, JPY.
:::

c. Make a biplot of the first two PCs. Explain what you learn.

```{r}
#| label: rates-biplot
#| code-fold: true
#| code-summary: "Biplot code"
#| eval: false
library(ggfortify)
autoplot(rates_pca, loadings = TRUE, 
         loadings.label = TRUE) 
```

::: unilur-solution
```{r}
#| label: rates-biplot
#| echo: false
```

- Most of the currencies contribute substantially to PC1. Only three contribute strongly to PC2: CHF, JPY, EUR. Similar to what is learned from the summary table (made in b).
- The pattern of the points is most unusual! It has a curious S shape. Principal components are supposed to be a random scattering of values, with no obvious structure. This is a very strong pattern.
:::

d. Make a time series plot of PC1 and PC2. Explain why this is useful to do for this data. 

```{r}
#| label: rates-PC-ts
#| code-fold: true
#| code-summary: "Code to plot PCs"
#| eval: false
rates_pca$x |>
  as.data.frame() |>
  mutate(date = rates_std$date) |>
  ggplot(aes(x=date, y=PC1)) + geom_line()

rates_pca$x |>
  as.data.frame() |>
  mutate(date = rates_std$date) |>
  ggplot(aes(x=date, y=PC2)) + geom_line()

```

::: unilur-solution
```{r}
#| label: rates-PC-ts
#| echo: false
```

- Because there is a strong pattern in the first two PCs, it could be useful to understand if this is related to the temporal context of the data. 
- Here we might expect that the PCs extract the main temporal patterns. We see this is the case.
- PC1 reflects the large group of currencies that greatly increase in mid-March.
- PC2 reflects the few currencies that decrease at the start of March. 

Note that: increase here means that the value of the currency declines relative to the USD and a decrease indicates stronger relative to the USD. *Is this correct?*
:::

e. You'll want to drill down deeper to understand what the PCA tells us about the movement of the various currencies, relative to the USD, over the volatile period of the COVID pandemic. Plot the first two PCs again, but connect the dots in order of time. Make it interactive with plotly, where the dates are the labels. What does following the dates tell us about the variation captured in the first two principal components?

```{r}
#| label: rates-PC-interactive
#| code-fold: true
#| code-summary: "Code to use interaction of the PC plot"
#| eval: false
library(plotly)
p2 <- rates_pca$x |>
  as.data.frame() |>
  mutate(date = rates_std$date) |>
  ggplot(aes(x=PC1, y=PC2, label=date)) +
    geom_point() +
    geom_path()
ggplotly(p2, width=400, height=400)
```

::: unilur-solution
```{r}
#| label: rates-PC-interactive
#| echo: false
```

The pattern in PC1 vs PC2 follows time. Prior to the pandemic there is a tangle of values on the left. Towards the end of February, when the world was starting to realise that COVID was a major health threat, there is a dramatic reaction from the world currencies, at least in relation to the USD. Currencies such as EUR, JPY, CHF reacted first, gaining strength relative to USD, and then they lost that strength. Most other currencies reacted later, losing value relative to the USD.
:::

#### 6. Write a simple question about the week's material and test your neighbour, or your tutor. 

## `r emo::ji("wave")` Finishing up

Make sure you say thanks and good-bye to your tutor. This is a time to also report what you enjoyed and what you found difficult.
