---
title: "Principle Components Analysis"
author: "Ryan G. McClarren"
date: "November 3, 2014"
output:
  html_document:
    toc: yes
---

Assume you have a lot of correlated data: what to do?

```{r, include=FALSE}
library(MASS)
library(ggplot2)
library(GGally)
library(reshape)
library(xtable)
```

#A simple 4-dimensional example

First, let's generate multivariate normal random samples from the 4-D space.  We have to create a covariance matrix and make sure it is positive definite:
```{r}
Ndims = 4
Nsamples = 25
Sigma <- matrix(runif(Ndims^2,min = 0.1,max=10),Ndims,Ndims)
Sigma <- Sigma %*% t(Sigma)
samples <- data.frame(mvrnorm(n=Nsamples, runif(4,min=-2,max=2), Sigma))
```

Looking at the data with ggpairs shows that the dimensions are correlated with each other:
```{r}
ggpairs(samples)
```

The next thing we do is compute the singular value decomposition of the data matrix, in this case, <tt> samples.</tt> It is useful to make the columns of the matrix have a zero mean to help interpret the singular values.

```{r}
#create mean 0 columns
norm_samples <- samples
columnMeans <- 0
for (i in 1:Ndims){
  columnMeans[i] <- mean(norm_samples[,i])
  norm_samples[,i] <- norm_samples[,i] - columnMeans[i]  
}

#compute SVD
sv_decomp <- svd(norm_samples)
svDF = data.frame(singular_values = sv_decomp$d, sv = 1:Ndims)
svDF$cumulativeFrac <- sv_decomp$d/sum(sv_decomp$d)
ggplot(svDF,aes(x=as.factor(sv),y=singular_values, label=paste(round(cumulativeFrac*100,2),"%"))) + geom_bar(stat="identity") + scale_x_discrete("Singular Value #") + scale_y_continuous("Singular Value") + geom_text(color="red",  aes(y=(sv_decomp$d)*.5))
```

From the figure, we see that the first singular value represents `r round(100*max(svDF$cumulativeFrac),2)`% of the cumulative sum of the singular values. One way to interpret this is the that `r round(100*max(svDF$cumulativeFrac),2)`% of the variability of the data is represented by a singular principle component.  The sum of the first two singular values is `r round(100*sum(svDF$cumulativeFrac[1:2]),2)`% of the cumulative sum of the singular values.  That is, about `r round(100*sum(svDF$cumulativeFrac[1:2]),2)`% of the total variability in the data is represented by two components.  

The upshot of this is that if we replace the 4-d space with an appropriate 2-d space, we can capture most of the affect of the 4-d space.

##One PC Reconstuction
Reconstructing the data using only one principle component requires using the singular value decomposition formula with all of the singular values, besides the first, being set to zero.

```{r}
#one_var reconstruction
normData1 <- sv_decomp$u %*% diag(c(sv_decomp$d[1],rep(0,Ndims-1))) %*% t(sv_decomp$v)
origScale1 <- data.frame(normData1)
for (i in 1:Ndims){
  origScale1[,i] <- origScale1[,i] + columnMeans[i]  
}

plotDF <- cbind(samples,origScale1)
colnames(plotDF) <- c("O1","O2","O3","O4","X1","X2","X3","X4")

ggplot(plotDF,aes(x=O1,y=X1)) + geom_point() + geom_smooth(method="lm",se=F) + scale_x_continuous("Original Data") + scale_y_continuous("1 Principle Component Reconstruction") + ggtitle(paste("Dimension 1, RMS Error:", round(sqrt(mean((plotDF$O1 - plotDF$X1)^2)),2)))
ggplot(plotDF,aes(x=O2,y=X2)) + geom_point() + geom_smooth(method="lm",se=F) + scale_x_continuous("Original Data") + scale_y_continuous("1 Principle Component Reconstruction") + ggtitle(paste("Dimension 2, RMS Error:", round(sqrt(mean((plotDF$O2 - plotDF$X2)^2)),2)))
ggplot(plotDF,aes(x=O3,y=X3)) + geom_point() + geom_smooth(method="lm",se=F) + scale_x_continuous("Original Data") + scale_y_continuous("1 Principle Component Reconstruction") + ggtitle(paste("Dimension 3, RMS Error:", round(sqrt(mean((plotDF$O3 - plotDF$X3)^2)),2)))
ggplot(plotDF,aes(x=O4,y=X4)) + geom_point() + geom_smooth(method="lm",se=F) + scale_x_continuous("Original Data") + scale_y_continuous("1 Principle Component Reconstruction") + ggtitle(paste("Dimension 4, RMS Error:", round(sqrt(mean((plotDF$O4 - plotDF$X4)^2)),2)))
```

##Two PC Reconstuction
We can do much better using 2 principle components:

```{r}
#one_var reconstruction
normData2 <- sv_decomp$u %*% diag(c(sv_decomp$d[1:2],rep(0,Ndims-2))) %*% t(sv_decomp$v)
origScale2 <- data.frame(normData2)
for (i in 1:Ndims){
  origScale2[,i] <- origScale2[,i] + columnMeans[i]  
}

plotDF <- cbind(samples,origScale2)
colnames(plotDF) <- c("O1","O2","O3","O4","X1","X2","X3","X4")

ggplot(plotDF,aes(x=O1,y=X1)) + geom_point() + geom_smooth(method="lm",se=F) + scale_x_continuous("Original Data") + scale_y_continuous("2 Principle Component Reconstruction") + ggtitle(paste("Dimension 1, RMS Error:", round(sqrt(mean((plotDF$O1 - plotDF$X1)^2)),2)))
ggplot(plotDF,aes(x=O2,y=X2)) + geom_point() + geom_smooth(method="lm",se=F) + scale_x_continuous("Original Data") + scale_y_continuous("2 Principle Component Reconstruction") + ggtitle(paste("Dimension 2, RMS Error:", round(sqrt(mean((plotDF$O2 - plotDF$X2)^2)),2)))
ggplot(plotDF,aes(x=O3,y=X3)) + geom_point() + geom_smooth(method="lm",se=F) + scale_x_continuous("Original Data") + scale_y_continuous("2 Principle Component Reconstruction") + ggtitle(paste("Dimension 3, RMS Error:", round(sqrt(mean((plotDF$O3 - plotDF$X3)^2)),2)))
ggplot(plotDF,aes(x=O4,y=X4)) + geom_point() + geom_smooth(method="lm",se=F) + scale_x_continuous("Original Data") + scale_y_continuous("2 Principle Component Reconstruction") + ggtitle(paste("Dimension 4, RMS Error:", round(sqrt(mean((plotDF$O4 - plotDF$X4)^2)),2)))
```

#Correcting Co-linearity with Principle Components
Given data that are nearly perfectly correlated, simple linear regression models can fall on their face. Here we create data that are nearly collinear, and then produce a dependent variable that has the simple formula $$y = x_1 + x_2 + 3.$$ 
```{r}
dHW <- data.frame(x1=0,x2=0,y=0)
dHW[1:20,"x1"] <- rnorm(20)
dHW[1:20,"x2"] <- rnorm(20,mean=dHW$x1,sd=0.01)
dHW[1:20,"y"] <- rnorm(20,mean=dHW$x1 + dHW$x2 + 3)
ggplot(dHW, aes(x=x1, y = x2)) + geom_point()
```

Running R's standard linear model on this gives something pretty bad:
```{r}
summary(simpMod <- lm(y ~ x1 + x2, data = dHW))
```

In other words we get $$y = `r coef(simpMod)[2]` x_1 + `r coef(simpMod)[3]` x_2 + `r coef(simpMod)[1]`.$$
The estimates of the standard error of the coefficients gives that neither of the coefficients for $x_1$ or $x_2$ are significant. 

We can fix this with principle components.  First, we compute the SVD,
```{r}
svdLin <- svd(dHW[,c("x1","x2")])
```
Then we multiply the original data by the matrix $V$:
```{r}
print(svdLin$v)
newData <- dHW
newData[,c("x1","x2")] <- as.matrix(dHW[,c("x1","x2")]) %*% svdLin$v
```
Finally, we build the model on the new variables
```{r}
summary(pcMod <- lm(y ~ x1 + x2, data = newData))
```
Now $x_1$ is the original value of $x_1$ multiplied by `r svdLin$v[1,1]` plus the original $x_2$ multiplied by `r svdLin$v[2,1]`, and the modified $x_2$ is the original value of $x_1$ multiplied by `r svdLin$v[1,2]` plus the original $x_2$ multiplied by `r svdLin$v[2,2]`.  We can solve for the coefficient in the original basis by solving the system of equations $$ `r svdLin$v[1,1]` x_1 + (`r svdLin$v[2,1]`) x_2 = `r coef(pcMod)[2]`,$$ $$ `r svdLin$v[1,2]` x_1 + (`r svdLin$v[2,2]`) x_2 = 0.$$ Here we have set the second coefficient to zero because it was not significant. 

The solution to this equation gives that the coefficient for $x_1$ in the original basis is `r solve((svdLin$v), c(coef(pcMod)[2],0))[1]` and the coefficient for $x_2$ in the original basis is `r solve((svdLin$v), c(coef(pcMod)[2],0))[2]` or 
$$y = `r solve((svdLin$v), c(coef(pcMod)[2],0))[1]` x_1 + `r solve((svdLin$v), c(coef(pcMod)[2],0))[2]` x_2 + `r coef(pcMod)[1]`.$$
This model is much closer to the true model that we used to generate the data.

#How should we measure success
Consider the following data from a course.
```{r}
grades <- read.csv("PostExam2Grades.csv")
summary(grades)
ggplot(grades,aes(x=Exam.1,y=Exam.2, color = Quiz.1 + Quiz.2 + Quiz.3, size = HW.1 + HW.2 + HW.3 + HW.4)) + geom_point() + scale_colour_gradient2("Quiz Total",low="red",mid="orange",high="blue",midpoint=20) + scale_x_continuous("Exam 1 Score", lim=c(40,100)) + scale_y_continuous("Exam 2 Score") + scale_size(range=c(1,12), "Homework Total") 
```

There are 9 total grades, grouped into 3 categories: homework, quizzes and exams.  Let's see how the different grades vary together:
```{r}
ggpairs(grades)
```


Let's throw a principle components analysis at this. We will use the built in R function <tt> prcomp<\tt>, and scale the data because the grades have different scales, e.g., quizzes are out of 10 and exams are out of 100. 
```{r}
print(gradePC <- prcomp(grades,scale=T))
summary(gradePC)
```

Based on the principle components analysis, there are only 5-7 grades hidden in the data. We can look at the rotation matrix to interpret the data.  For instance, PC1 is a measure of the overall performance of the student (it is a weighted sum of all the grades). Therefore, the higher this score, the higher the overall score of the student. 

PC2 is negative for exams and quizzes and positive for homeworks.  Students with a high PC2 outperformed on homework relative to exams and quizzes; a negative PC2 score would indicate they rose to the occasion on quizzes and exams.  The presence of a measure such as this indicates that there was a measureable partition of the class that had a different type of performance on homework versus quizzes and exams. 

The other PC's can be interpreted also.  The benefit of this kind of approach is that it makes the story the data is telling, a bit more clear because you have transformed to a uncorrelated basis.  As we saw in the figures above, There are strong correlations in the data.
