---
title: "Multivariate Demos"
output: html_notebook
---

This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook. When you execute code within the notebook, the results appear beneath the code.

First, let's get the session set up. These two commands will clear out any leftover variables from a previous session and load a library we're using.

```{r}
rm(list = ls())
library(mvtnorm)
```

Next, let's set a couple of simulation parameters. We are going to be working with two variables to keep it simple, and the relationship between them can be defined with a single correlation value which is set using the `input_correlation` variable. The default value is 0.8, but after going through this the first time, you can play with this to see how it affects the answer. What is the allowable range for this value?

The second variable `sample_size` lets you set the sample size for the simulation, which affects the precision of the results. The default is 1,000 but you can change this as you like as well.

```{r}
input_correlation <- 0.8
sample_size <- 1000
```

The following code sets up the parameter inputs for the bivariate normal distribution. `x.mean` is a vector of zeroes with length 2 - one for each dimension of the bivariate normal. `x.cov` is a 2x2 (again, dimension 2) matrix where there the values on the diagonal are all 1 (because the variance of each variable is 1) and the off-diagonals have our correlation value. *If we used values on the diagonal that are not 1, then the off-diagonal values have to be covariance, not correlation.*

```{r}
x.mean <- c(0, 0)
x.cov <- matrix(c(1, input_correlation, input_correlation, 1), nrow = 2)
```

The following code generates a sample from the bivariate normal distribution we just specified using the `rmvnorm` function from the `mvtnorm` library. Note that the result is a matrix with 2 columns and a number of rows equal to the sample size.

```{r}
Z <- rmvnorm(sample_size, mean = x.mean, sigma = x.cov)  # note that the result is a matrix with 2 columns and "sample_size" rows
```

Let's look at the sample correlation values using the three measures we've talked about. The code uses matrix notation to get one value because the `cor` function returns a symmetric matrix, and we want the off-diagonal values. For large samples, Kendall's tau takes a little while to compute.

```{r}
cor(Z, method = "pearson")[1,2]
cor(Z, method = "spearman")[1,2]
cor(Z, method = "kendall")[1,2]
```

Because the data are fairly linear, the Pearson correlation and Spearman's rho have similar values. Kendall's tau is a bit lower, which is not unusual. How close was the Pearson correlation to your input value?

Now, let's look at a very simple scatter plot of the data. The line below uses R's base graphics.

```{r}
plot(Z, xlab = "X1", ylab = "X2", main = "Bivariate Normal Sample")
```

You should notice that the values cluster around the origin (0, 0) because the mode of the bivariate normal distribution is right at the mean vector we specified earlier. Depending on the value for the correlation coefficient you chose, the plot will either seem strongly related (as in the default case of 0.8) or may behave differently. If you change the `input_correlation` value above, re-run all the chunks of code above to update the sample and the plot using the "Run All Chunks Above" button (grey downward arrow with green line button) in this chunk.

Now, let's visualize the distribution that we took the samples from. We will represent this bivariate normal distribution with some contour lines and place it over our sampled points to see how they line up.

```{r}
z1.seq <- seq(min(Z[,1]), max(Z[,1]), length.out=50)
z2.seq <- seq(min(Z[,2]), max(Z[,2]), length.out=50)
density_values <- dmvnorm(expand.grid(z1.seq, z2.seq), x.mean, x.cov)
dim(density_values) <- c(length(z1.seq), length(z2.seq))
plot(Z, xlab = "X1", ylab = "X2", main = "Bivariate Normal Sample")
contour(z1.seq, z2.seq, density_values, col="#A60F2D", drawlabels=FALSE, nlevels=10, lwd=1, lty = 2, add = TRUE)
```

Due to the way we set up our bivariate normal distribution (mean vector with zeroes, ones on the diagonal of the covariance matrix) we know that our two variables should have a standard normal distribution. (Recall that a standard normal is normal with mean = 0 and variance = 1.) Let's compare the expected distribution to the sampled ones.

First we compute the standard normal density using the `dnorm` function for a range of values and plot that in black. Then, we use the `density` function to compute a kernel density estimate of the distribution of the two samples in blue and red.

```{r}
z.ref <- seq(-5, 5, length.out = 100)
plot(z.ref, dnorm(z.ref, mean = 0, sd = 1), type = "l", lwd = 3, lty = 2, main = "Marginal Densities", xlab = "Value", ylab = "f(x)")
lines(density(Z[,1]), col = "blue", lwd = 2)
lines(density(Z[,2]), col = "red", lwd = 2)
```

Now, let's look at a regression line for X1 predicting the value of X2 using simple linear regression. The `lm` function fits a linear model using some funky notation we will see more of in a later workshop. `abline` is an easy way to add a fitted regression model to an existing plot. The `summary` function produces the standard regression output information.

```{r}
Z.df <- data.frame(Z)
Z.lm <- lm(X2 ~ X1, data = Z.df)
plot(Z, xlab = "X1", ylab = "X2", main = "Bivariate Normal Sample with Regression Line")
abline(Z.lm)
summary(Z.lm)
```

The regression summary shows three interesting things:

1.  The coefficient of the regression model is close to the value of the correlation coefficient you specified.

2.  The intercept of the model is close to zero because the value of X2 when X1 is zero should also be zero.

3.  The model r-squared is close to the value of your chosen correlation coefficient, squared. The reason these are close and not exact is due to sample error. If you choose a larger sample size, the sample values get closer to the true values.

There is a close relationship between the bivariate normal distribution and simple linear regression.

Finally, let's look at you you would use a bivariate normal distribution as part of a Gaussian copula. The marginal distributions for our two variables are standard normal as we saw earlier. Copulas let us use any probability distribution for our margins, so we have to transform the variables from standard normal to our target distributions. The first step is to get into probability space by using the CDF of the standard normal distribution with the `pnorm` function. Let's look at what the data look like before and after this transformation.

```{r}
P <- pnorm(Z, mean = 0, sd = 1)
plot(Z, xlab = "X1", ylab = "X2", main = "Bivariate Normal Sample")
plot(P, xlab = "F(X1)", ylab = "F(X2)", main = "Bivariate Normal Sample in Probability Space")
```

Let's look at how the correlation measures change after transformation.

```{r}
cor(P, method = "pearson")[1,2]
cor(P, method = "spearman")[1,2]
```

Notice how the Spearman correlation didn't change, but the Pearson correlation decreased. Spearman doesn't care about monotonic transformations of the data because it is based on rank.

Finally, let's transform these margins into other distributions. We will use a gamma(1, 3) distribution for X1 and beta(2, 2)distribution for X2. We use the quantile function for these distributions (`qbeta` and `qgamma`, respectively) to take the CDF values and turn them into the target distributions.

```{r}
y1 <- qgamma(P[,1], 1, 3)
y2 <- qbeta(P[,2], 2, 2)
plot(y1, y2, main = "Bivariate Normal Sample with Transformed Margins")
```

We can also look at the marginal densities for these two variables.

```{r}
plot(density(y1), type = "l", col = "blue", main = "Transformed Marginal Densities")
lines(density(y2), col = "red")
```

Notice how they are no longer standard normal distributions. Even if you're not familiar with the shape of the beta or gamma distributions, you can see that our transformations have an effect!