##PCA with Expedia destinations
##load destinations.cvs data
dest <- read.csv("../input/destinations.csv", header = T, sep=",", stringsAsFactors = F)
head(dest, 10)
apply(dest, 2, mean)
apply(dest, 2, var)
dim(dest)

require(leaps)

##I took little random sample for easy emplementation
set.seed(12345678)
mysam1  <- dest[sample(1:nrow(dest), 50, replace= FALSE),]
summary(mysam1)
## 1 forward
forward = regsubsets(d1~., mysam1, method ="forward")
summary(forward)
## its get me d35, d52,d57, d62, d89, d119, d145 forsubset to use for PCA
mysam2 <- mysam1[, c(36,53,58, 63, 90, 120, 146)]
mysam2
summary(mysam2)


## PCA with forward
pcaex1.out=prcomp(mysam2, scale=TRUE)
pcaex1.out
hist(pcaex1.out$rotation)

library(devtools)


pca_dest1 <- princomp(mysam2, cor=TRUE)
summary(pca_dest1)
loadings(pca_dest1) # pc loadings 
plot(pca_dest1,type="lines") # scree plot 
pca_dest1$scores # the principal components
biplot(pca_dest1, cex = .5)

predict(pca_dest1, 
        newdata=tail(dest, 10))
        
library(devtools)
install_github("ggbiplot", "vqv")

install.packages("vqv-ggbiplot-2623d7c.tar.gz", repos=NULL, type="source")
##Installing ggbiplot.zip from https://github.com/vqv/ggbiplot/zipball

g <- ggbiplot(pca_dest1, obs.scale = 1, var.scale = 1, 
               ellipse = TRUE, 
               circle = TRUE)
g <- g + scale_color_discrete(name = '')
g <- g + theme(legend.direction = 'horizontal', 
               legend.position = 'top')
##print(g)

## 2 backward
backward = regsubsets( d1~., mysam1, nbest = 1, nvmax = 8, method ="backward")
summary(backward)

## backward selection gave me d12, d13,d18, d22, d24, d33, d36  forsubset to use for PCA
mysam3 <- mysam1[, c(13, 14, 19, 23, 25, 34, 37)]
mysam3
summary(mysam3)
## PCA

pcaex2.out=prcomp(mysam3, scale=TRUE)
pcaex2.out
hist(pcaex2.out$rotation)

pca_dest2 <- princomp(mysam3, cor=TRUE)
summary(pca_dest2)
loadings(pca_dest2) # pc loadings 
plot(pca_dest2,type="lines") # scree plot 
pca_dest2$scores # the principal components
biplot(pca_dest2, cex = .5)

## All big data set PCA
pca_dest <- princomp(dest, cor=TRUE)
pca_dest00 <- princomp(dest, cor=TRUE, scores = T)
summary(pca_dest)
loadings(pca_dest)             # pc loadings 
pca_dest$sd^2                  # component variances
screeplot(pca_dest)               # if you would like a scree plot
plot(pca_dest,type="lines")    # scree plot in lines
pca_dest$scores             # the principal components
##biplot(pca_dest)  is to big to emplement but possible
biplot(pca_dest$sdev, pca_dest$scale)
screeplot(pca_dest)
 
pca.m4 <- princomp(~d98,
                   cor=TRUE, data=dest, scores = TRUE)
unclass(loadings(pca.m3))  # component loadings
unclass(loadings(pca.m4))
pca.m3$sd^2  # component variances
summary(pca.m3) # proportions of variance
screeplot(pca.m3) # if you would like a scree plot

##extract the component scores
first.component.scores <- pca_dest$scores[,1]
summary(first.component.scores)
length(first.component.scores)

str(pca_dest)  ## show the structure

second.component.scores <- pca_dest$scores[,2]
summary(second.component.scores)
length(second.component.scores)

# Create a correlation matrix (also called a matrix of association) object 
# of the items for PCA.
library(psych)
cor.matrix.1 <- cor(dest)
cor.matrix.1

# Script for PCA with VARIMAX rotation (an orthogonal rotation strategy).

pca.2 <- principal(r = cor.matrix.1, nfactors = 150, residuals = FALSE, rotate = "varimax")
pca.2
pca.varimax <- principal(r = cor.matrix.1, nfactors = 10, residuals = FALSE, rotate = "varimax")
pca.varimax  ## for 10 n factors

library(GPArotation)
# Script for PCA with Direct Oblimin rotation (an oblique rotation strategy).

pca.3 <- principal(r = cor.matrix.1, nfactors = 10, residuals = FALSE, rotate = "oblimin")
pca.3

##Cor matix with forward dataset
cor.matrix.2 <- cor(mysam1[,c("d35", "d52","d57", "d62", "d89", "d119", "d145")], method = "spearman")
cor.matrix.2
# Predict PCs and ggbiplot
predict(pca_dest, 
        newdata=tail(dest, 10))
ggbiplot
gd <- ggbiplot(pca_dest, obs.scale = 1, var.scale = 1, 
              ellipse = TRUE, 
              circle = TRUE)
gd <- gd + scale_color_discrete(name = '')
gd <- gd + theme(legend.direction = 'horizontal', 
               legend.position = 'top')
print(gd)