R语言 主成分分析PCA（绘图+原理）

# install.packages("devtools")
# devtools::install_github("vqv/ggbiplot")
library(ggbiplot)
data(wine)
head(wine)
# > head(wine)
# Alcohol MalicAcid  Ash AlcAsh  Mg Phenols Flav NonFlavPhenols Proa Color  Hue   OD Proline
# 1   14.23      1.71 2.43   15.6 127    2.80 3.06           0.28 2.29  5.64 1.04 3.92    1065
# 2   13.20      1.78 2.14   11.2 100    2.65 2.76           0.26 1.28  4.38 1.05 3.40    1050
# 3   13.16      2.36 2.67   18.6 101    2.80 3.24           0.30 2.81  5.68 1.03 3.17    1185
# 4   14.37      1.95 2.50   16.8 113    3.85 3.49           0.24 2.18  7.80 0.86 3.45    1480
# 5   13.24      2.59 2.87   21.0 118    2.80 2.69           0.39 1.82  4.32 1.04 2.93     735
# 6   14.20      1.76 2.45   15.2 112    3.27 3.39           0.34 1.97  6.75 1.05 2.85    1450

wine.pca <- prcomp(wine, scale. = TRUE)
## bioplot
ggbiplot(wine.pca, obs.scale = 1, var.scale = 1,
         groups = wine.class, ellipse = TRUE, circle = TRUE) +
  scale_color_discrete(name = '') +
  theme_light()

ggbiplot(pcobj, choices = 1:2, scale = 1, pc.biplot =
TRUE, obs.scale = 1 - scale, var.scale = scale, groups =
NULL, ellipse = FALSE, ellipse.prob = 0.68, labels =
NULL, labels.size = 3, alpha = 1, var.axes = TRUE, circle
= FALSE, circle.prob = 0.69, varname.size = 3,
varname.adjust = 1.5, varname.abbrev = FALSE, ...)
pcobj # prcomp()或princomp()返回结果
choices # 选择轴，默认1：2
scale # covariance biplot (scale = 1), form biplot (scale = 0). When scale = 1, the inner product between the variables approximates the covariance and the distance between the points approximates the Mahalanobis distance.
obs.scale # 标准化观测值
var.scale # 标准化变异
pc.biplot # 兼容 biplot.princomp()
groups # 组信息，并按组上色
ellipse # 添加组椭圆
ellipse.prob # 置信区间
labels #向量名称
labels.size #名称大小
alpha  #点透明度 (0 = TRUEransparent, 1 = opaque)
circle #绘制相关环(only applies when prcomp was called with scale = TRUE and when var.scale = 1)
var.axes  #绘制变量线-菌相关
varname.size  #变量名大小
varname.adjust #标签与箭头距离 >= 1 means farther from the arrow
varname.abbrev # 标签是否缩写

library(psych)
library(ggplot2)
PC <- principal(wine, nfactors=2, rotate ="none")
pc <- data.frame(PC$scores)
p  <- ggplot(pc, aes(x=PC1, y=PC2,color=wine.class )) +
  geom_point(size=4,alpha=0.5)+
  theme(axis.text= element_text(size=20))+
  theme(panel.grid.major =element_blank(),
        panel.grid.minor = element_line(color="steelblue"),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"))+
  stat_ellipse(lwd=1,level = 0.8)

data(wine)
wine.pca <- prcomp(wine, scale. = TRUE)
## bioplot
TB= data.frame(wine.pca$x)
TB$class = wine.class

p0 <- ggbiplot(wine.pca, obs.scale = 1, var.scale = 1,
               groups = wine.class, ellipse = TRUE, circle = TRUE) +
  scale_color_discrete(name = '') +
  theme_light()+theme(legend.position = 'none',
                      axis.title=element_blank())

p1 <- ggplot() +
  geom_density(data=TB, aes(x=PC1,fill=class),alpha=0.5) +
  theme_map() +
  theme(legend.position = 'none')

p2 <- ggplot() +
  geom_density(data=TB, aes(x=PC2,fill=class),alpha=0.5)  +
  theme_map() +theme(legend.position = 'none') +
  coord_flip()  +  scale_y_reverse()

p3 <- ggplot(TB, aes(x=class,y=PC1)) +
  geom_tufteboxplot(aes(color=class),median.type = "line",
                    hoffset = 0, width = 3) +
  coord_flip()+ theme_map()+
  theme(legend.position = 'none')

p4 <- ggplot(TB, aes(x=class,y=PC2)) +
  geom_tufteboxplot(aes(color=class),median.type = "line",
                    hoffset = 0, width = 3) +
  theme_map()+
  theme(legend.position = 'none')


TB= data.frame(wine.pca$x)
TB$class = wine.class

GGlay <-
  "#BBBBBB#
CAAAAAAE
CAAAAAAE
CAAAAAAE
CAAAAAAE
CAAAAAAE
CAAAAAAE
#DDDDDD#"
p0+p1+p2+p3+p4+plot_layout(design = GGlay)

library(dplyr)
#特征分解
dat_eigen<-scale(iris[,-5],scale=T)%>%cor()%>%eigen()
#特征值提取
dat_eigen$values 
#特征向量提取
dat_eigen$vectors
#求loadings=eigen vector*sqrt(eigen value)，与princomp不同
#主成分载荷表示各个主成分与原始变量的相关系数。
sweep(dat_eigen$vectors,2,sqrt(dat_eigen$values),"*")
#将中心化的变量矩阵得到每个观测值的得分
scale(iris[,-5],scale=T)%*%dat_eigen$vectors%>%head()

iris.pca<-prcomp(iris[,-5],scale=T,rank=4,retx=T) #相关矩阵分解
#retx表四返回score，scale表示要标准化
summary(iris.pca) #方差解释度
iris.pca$sdev #特征值的开方
iris.pca$rotation #特征向量，回归系数
iris.pca$x #样本得分score

data(wine) #三种葡萄酿造的红酒品质分析数据集
wine.pca<-princomp(wine,cor=T,scores=T) 
#默认方差矩阵(cor=F),改为cor=T则结果与prcomp相同
summary(wine.pca) #各主成份的SVD值以及相对方差
wine.pca$loading #特征向量，回归系数
wine.pca$score
screenplot(wine.pca) #方差分布图
biplot(wine.pca,scale=F) #碎石图,直接把x与rotation绘图，而不标准化

model_pca<-psych::principal(iris[,-5],nfactors=4,rotate="none")
model_pca$values # 特征值=sdev^2
# 此处loadings与手动计算相同=prcomp的rotation*sdev 
model_pca%>%.$loadings #载荷，不是特征向量
#此处score=prcomp的score/sdev
model_pca$scores[1:5,] #此处为因子得分，不是主成分得分
model_pca$weights #此处为上面loadings/特征值，也称成份得分系数或者因子系数

prcomp函数会返回主成分的标准差、特征向量和主成分构成的新矩阵。
不同主成分对数据差异的贡献和主成分与原始变量的关系。
1. 主成分的平方为为特征值,其含义为每个主成分可以解释的数据差异，计算方式为 eigenvalues = (pca$sdev)^2
2. 每个主成分可以解释的数据差异的比例为 percent_var = eigenvalues*100/sum(eigenvalues)
3. 可以使用summary(pca)获取以上两条信息。

这两个信息可以判断主成分分析的质量：
    成功的降维需要保证在前几个为数不多的主成分对数据差异的解释可以达到80-90%。

指导选择主成分的数目：
  1.  选择的主成分足以解释的总方差大于80% (方差比例碎石图)
  2. 从前面的协方差矩阵可以看到，自动定标(scale)的变量的方差为1 (协方差矩阵对角线的值)。待选择的主成分应该是那些方差大于1的主成分，即其解释的方差大于原始变量（特征值碎石图，方差大于1，特征值也会大于1，反之亦然）。

鉴定核心变量和变量间的隐性关系:
    原始变量与主成分的相关性 Variable correlation with PCs (var.cor) = loadings * sdev
    原始数据对主成分的贡献度 var.cor^2 / (total var.cor^2)