降维与度量学习
最后发布时间 : 2023-10-10 23:08:01
浏览量 :
学习资料
- R语言进行PCA(主成分分析)
- Principal Component Analysis through Singular Value Decomposition
- Principal Component Methods in R: Practical Guide
使用函数prcomp计算主成分
# 加载 iris 数据集
data(iris)
iris数据集总共有5列，其中前四列是特征，最后一列是标签
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
提取前四列作为PCA的输入，并为每一个样本生成一个行名
# 提取特征变量
x <- iris[, 1:4]
rownames(x) <- paste("sample",1:nrow(x),sep = "")
Sepal.Length Sepal.Width Petal.Length Petal.Width
sample1 5.1 3.5 1.4 0.2
sample2 4.9 3.0 1.4 0.2
sample3 4.7 3.2 1.3 0.2
sample4 4.6 3.1 1.5 0.2
sample5 5.0 3.6 1.4 0.2
sample6 5.4 3.9 1.7 0.4
使用函数prcomp执行pca
# 执行主成分分析
pca <- prcomp(x, scale= TRUE,center = TRUE)
打印pca变量，可以看到分析结果包括主成分的标准差和Rotation矩阵，其中Rotation矩阵的行是特征，列是主成分
稍后可以看到
- 主成分的标准差是：相关系数矩阵特征值的标准差
Rotation矩阵是：相关系数矩阵特征向量
Standard deviations (1, .., p=4):
[1] 1.7083611 0.9560494 0.3830886 0.1439265
Rotation (n x k) = (4 x 4):
PC1 PC2 PC3 PC4
Sepal.Length 0.5210659 -0.37741762 0.7195664 0.2612863
Sepal.Width -0.2693474 -0.92329566 -0.2443818 -0.1235096
Petal.Length 0.5804131 -0.02449161 -0.1421264 -0.8014492
Petal.Width 0.5648565 -0.06694199 -0.6342727 0.5235971
提取降维后的主成分矩阵
res1 <- as.data.frame(pca$x)
PC1 PC2 PC3 PC4
sample1 -2.257141 -0.4784238 0.12727962 0.024087508
sample2 -2.074013 0.6718827 0.23382552 0.102662845
sample3 -2.356335 0.3407664 -0.04405390 0.028282305
sample4 -2.291707 0.5953999 -0.09098530 -0.065735340
sample5 -2.381863 -0.6446757 -0.01568565 -0.035802870
sample6 -2.068701 -1.4842053 -0.02687825 0.006586116
这个主成分矩阵是如何计算得到？
scale(x, center = TRUE, scale = TRUE) %*% as.matrix(pca$rotation)
可以看到是中心化后的特征变量(150 * 4)乘以rotation矩阵(4 * 4)，也就是说，只要计算得到rotation矩阵就可以得到主成分矩阵
> scale(x, center = TRUE, scale = TRUE)
Sepal.Length Sepal.Width Petal.Length Petal.Width
sample1 -0.89767388 1.01560199 -1.33575163 -1.3110521482
sample2 -1.13920048 -0.13153881 -1.33575163 -1.3110521482
> as.matrix(pca$rotation)
PC1 PC2 PC3 PC4
Sepal.Length 0.5210659 -0.37741762 0.7195664 0.2612863
Sepal.Width -0.2693474 -0.92329566 -0.2443818 -0.1235096
Petal.Length 0.5804131 -0.02449161 -0.1421264 -0.8014492
Petal.Width 0.5648565 -0.06694199 -0.6342727 0.5235971
提取主成分的方差贡献率
summary(pca)
- 其中第一行是相关系数矩阵特征值的标准差
- 第二行是相关系数矩阵特征值的比例，也就是主成分的方差贡献率，也可以通过第一行计算得到
pca$sdev^2 / sum(pca$sdev^2) - 第三行是累计比例：0.7296 +0.2285 =0.9581
Importance of components:
PC1 PC2 PC3 PC4
Standard deviation 1.7084 0.9560 0.38309 0.14393
Proportion of Variance 0.7296 0.2285 0.03669 0.00518
Cumulative Proportion 0.7296 0.9581 0.99482 1.00000
使用相关系数计算Rotation矩阵
首先计算相关系数矩阵
cor_mat <- cor(x)
Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 1.0000000 -0.1175698 0.8717538 0.8179411
Sepal.Width -0.1175698 1.0000000 -0.4284401 -0.3661259
Petal.Length 0.8717538 -0.4284401 1.0000000 0.9628654
Petal.Width 0.8179411 -0.3661259 0.9628654 1.0000000
求解相关系数矩阵的特征值和特征向量
rs_mat <- eigen(cor_mat)
可以看到其中特征向量矩阵与Rotation矩阵是相同的
eigen() decomposition
$values
[1] 2.91849782 0.91403047 0.14675688 0.02071484
$vectors
[,1] [,2] [,3] [,4]
[1,] 0.5210659 -0.37741762 0.7195664 0.2612863
[2,] -0.2693474 -0.92329566 -0.2443818 -0.1235096
[3,] 0.5804131 -0.02449161 -0.1421264 -0.8014492
[4,] 0.5648565 -0.06694199 -0.6342727 0.5235971
计算主成分的方差贡献率
对特征值取平方根，即为函数prcomp结果的Standard deviation
sqrt(rs_mat$values)
[1] 1.7083611 0.9560494 0.3830886 0.1439265
计算特征值的比例，即为主成分的方差贡献率
rs_mat$values/sum(rs_mat$values)
[1] 0.729624454 0.228507618 0.036689219 0.005178709
使用奇异值分解(SVD)计算Rotation矩阵
# 对特征变量进行中心化
x_centered <- scale(x, center = TRUE, scale = TRUE)
svd_result <- svd(x_centered)
# 提取奇异值、左奇异向量和右奇异向量
singular_values <- svd_result$d
left_singular_vectors <- svd_result$u
right_singular_vectors <- svd_result$v
可以看多SVD的右矩阵与Rotation矩阵相同
> singular_values (1 * 4)
[1] 20.853205 11.670070 4.676192 1.756847
> head(left_singular_vectors) (150*4)
[,1] [,2] [,3] [,4]
[1,] -0.10823953 -0.04099580 0.027218646 0.013710648
[2,] -0.09945776 0.05757315 0.050003401 0.058435855
[3,] -0.11299630 0.02920003 -0.009420891 0.016098333
[4,] -0.10989710 0.05101939 -0.019457133 -0.037416661
[5,] -0.11422046 -0.05524180 -0.003354363 -0.020379051
[6,] -0.09920300 -0.12718049 -0.005747892 0.003748828
> right_singular_vectors (4*4)
[,1] [,2] [,3] [,4]
[1,] 0.5210659 -0.37741762 0.7195664 0.2612863
[2,] -0.2693474 -0.92329566 -0.2443818 -0.1235096
[3,] 0.5804131 -0.02449161 -0.1421264 -0.8014492
[4,] 0.5648565 -0.06694199 -0.6342727 0.5235971
计算主成分的方差贡献率
singular_values^2/sum(singular_values^2)
[1] 0.729624454 0.228507618 0.036689219 0.005178709
主成分可视化
使用函数biplot
biplot(pca)
使用ggplot
library(ggplot2) # 加载ggplot2库
summ1 <- summary(pca) # 对pca进行汇总统计
xlab1 <- paste0("PC1(",round(summ1$importance[2,1]*100,2),"%)") # 设置x轴标签
ylab1 <- paste0("PC2(",round(summ1$importance[2,2]*100,2),"%)") # 设置y轴标签
ggplot(data = res1, aes(x = PC1, y = PC2, color = iris$Species)) + # 创建ggplot对象，设置数据和映射
stat_ellipse(aes(fill = iris$Species), type = "norm", geom = "polygon", alpha = 0.25, color = NA) + # 添加置信椭圆
geom_point(size = 3.5) + # 添加散点图
labs(x = xlab1, y = ylab1, color = "Condition", title = "PCA Scores Plot") + # 设置标签和标题
guides(fill = "none") + # 隐藏填充颜色的图例
theme_bw() + # 设置主题为白色背景
scale_fill_manual(values = c("purple", "orange", "pink")) + # 设置填充颜色的手动比例尺
scale_colour_manual(values = c("purple", "orange", "pink")) + # 设置颜色的手动比例尺
theme(plot.title = element_text(hjust = 0.5, size = 15),
axis.text = element_text(size = 11),
axis.title = element_text(size = 13),
legend.text = element_text(size = 11),
legend.title = element_text(size = 13),
plot.margin = unit(c(0.4, 0.4, 0.4, 0.4), 'cm')) # 设置图表的主题样式
PCA Loadings
pca$rotation
因子负荷量是一个矩阵，每一列代表一个主成分，每一行代表一个原始变量。因子负荷量的元素表示了每个原始变量对应主成分的权重或相关性。
PC1 PC2 PC3 PC4
Sepal.Length 0.5210659 -0.37741762 0.7195664 0.2612863
Sepal.Width -0.2693474 -0.92329566 -0.2443818 -0.1235096
Petal.Length 0.5804131 -0.02449161 -0.1421264 -0.8014492
Petal.Width 0.5648565 -0.06694199 -0.6342727 0.5235971