PCA
Description for PCA notebook.
Notebook Contents
This notebook covers:
- Topic 1
- Topic 2
- Topic 3
Use the buttons above to download the notebook or open it in your preferred environment.
📓 Notebook Preview
M9a: PCA example using hubway data¶
library(mvtnorm) # col.norm
library(tilting) # col.norm
library(ggplot2)
library("np") #npreg: you may need to install it
hpdata <- read.csv("./data/dataset-hubway2.csv", header=TRUE, sep = ",", quote="\"",dec = ".")
hpdata$TOWN <- NULL
hpdata$year <- NULL #zero-variance
#write.table(hpdata,"./output/dataset-hubway3.csv", sep = ",",dec = ".")
colnames(hpdata)
# PCA Scree plot without standardizing data
hpca_cov <- prcomp(hpdata, scale=FALSE)
# We can extract the information summarized above (and much more)
# from the attributes of the object hpca_cov
standard_deviation_of_each_component <- hpca_cov$sdev
var_per_dim <- standard_deviation_of_each_component^2
var_tot <- sum(var_per_dim)
var_tot
var_per_dim/var_tot
var_prop <- var_per_dim / sum(var_per_dim)
var_prop
cum_var <- cumsum(var_prop)
cum_var
plot(cum_var,xlab="Principal component",
ylab="Proportion of variance explained", ylim=c(0,1), type='b')
apply(hpdata, 2, mean)
apply(hpdata, 2, var)
From now on, we use only the correlation-matrix-based PCA (standardized).
First, we describe with plots and comment on how the first two principal components relate to the original columns in both covariance-matrix-based PCA and correlation-matrix-based (standardized) PCA. In both cases, is it possible to give a name to each of the two principal components?
hpca_cor <- prcomp(hpdata, scale=TRUE)
standard_deviation_of_each_component <- hpca_var$sdev
var_per_dim <- standard_deviation_of_each_component^2
var_tot <- sum(var_per_dim)
var_prop <- var_per_dim / sum(var_per_dim)
cum_var <- cumsum(var_prop)
plot(cum_var,xlab="Principal component",
ylab="Proportion of variance explained", ylim=c(0,1), type='b')
eigenvectors <- hpca_cor$rotation
col.norm(eigenvectors)
eigenvectors
# Let us plot the contribution of the original dimension to the 1st PCA
PC_contr <- eigenvectors[,c("PC1")]
# PC_contr
# We order by the magnitude of the contribution
# We use the - sign because we want a descending order
ord <- order( -abs(PC_contr) )
PC_contr <- PC_contr[ord]
#PC_contr
# We just select the 5 highest contributing dimensions (highest loadings)
PC_contr <- PC_contr[1:5]
PC_contr
barplot(PC_contr, main="Contribution to the 1st component", xlab="Original Dimensions")
Since the household size-related variables have negative loadings, we can describe the first principal component as the "Small low-car-usage household"
# Second principal component vector
PC_contr <- eigenvectors[,c("PC2")]
# We order by the magnitude of the contribution
# We use the - sign because we want a descending order
ord <- order( -abs(PC_contr) )
PC_contr <- PC_contr[ord]
# We just select the 5 highest contributing dimensions
PC_contr <- PC_contr[1:5]
options(repr.plot.width=12, repr.plot.height=8)
barplot(PC_contr, main="Contribution to the 2nd component",xlab="Original Dimensions")
Based on the value of the important loadings, the second component can be described as "Bike Usage" or "Bike unavailability"
scores <- hpca_cor$x
par(mfrow = c(1,2))
options(repr.plot.width=16, repr.plot.height=8)
plot(scores[,1], hpdata$num_of_pickups, xlab="Small low-car-usage household propensity",
ylab="num_of_pickups",cex=1.2,cex.lab=1.2)#,xlim=c(-1000,0))
plot(scores[,2], hpdata$num_of_pickups, xlab="Bike unavailability",
ylab="num_of_pickups",cex=1.2,cex.lab=1.2)
We show with plots and comments if the variable num of bikes can be explained by the first two principal components. We visualize in the same plot the first 2 components and some indication of the number of pickups.
options(repr.plot.width=15, repr.plot.height=10)
df <- data.frame(SmallHH=scores[,1], Unavailability=scores[,2],
Num.pickups=hpdata$num_of_pickups)
#quantile(df$pickups, 0.25)
df$Pickups.yes <- cut(df$Num.pickups, c(-Inf, 0, Inf)) # create category for pickups (yes) or no pickups
#df <- df[which(hpdata$hour>=8 & hpdata$hour<=10),]
ggplot(df, aes(x=SmallHH, y=Unavailability) ) + geom_point( aes(size=Num.pickups, color=Pickups.yes),alpha=0.75 ) + theme_gray(base_size = 25)
Comment on your observations.