Exploratory data analysis 3rd
Dr. Ning-ning ZHOU, Yunnan university
2023-10-12
统计术语
方差：衡量数据分布的离散程度。它描述的是一组数据与其平均值之间的偏离程度。
协方差： 描述的是变量间如何一起变化。协方差的值如果为正值，则说明两者是正相关的，结果为负值就说明负相关的，如果为0，也是就是统计上说的“相互独立”。它与相关性（Correlation）的主要区别是，协方差未被标准化，取值可以从“-∞ ~ ∞”间的任何数值，而相关性是标准化的，相关系数的取值范围是”-1 ~ 1”。
协方差矩阵: 将所有变量之间的协方差关系以矩阵的形式表现，即为协方差矩阵。
主成分分析（Principal Component Analysis, PCA）：是一种数据降维技巧，它能将大量相关变量转换为一组很少的不相关变量，这些无关变量称为主成分。例如，PCA可以将30个相关的环境变量转化为5个无关的变量，并且尽可能的保留原始数据集的信息。
PCA的目标是用一组较少的不相关变量代替大量相关变量，同时尽可能保留初始变量的信息，这些推导所得的变量称为主成分，它们是观测变量的线性组合。第一主成分可以等效地定义为使投影数据的方差最大化的方向。降维的直观感受如下图：
安装R包
已安装过的，无需重复安装
install.packages("ggplot2")
install.packages("ggpubr")
install.packages("cowplot")
install.packages("psych")
install.packages("devtools")
install.packages("factoextra")加载数据
data("iris")
head(iris[c(1, 2, 51, 52, 101, 102), ])## 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
## 51 7.0 3.2 4.7 1.4 versicolor
## 52 6.4 3.2 4.5 1.5 versicolor
## 101 6.3 3.3 6.0 2.5 virginica
## 102 5.8 2.7 5.1 1.9 virginicastr(iris)## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...4 探索性数据分析
4.8 PCA的基本原理
视频教程推荐： StatQuest with Josh Starmer@YouTube: https://www.youtube.com/watch?v=FgakZw6K1QQ&t=67s
交通数据小旭学长@bilibili: https://www.bilibili.com/video/BV1E5411E71z/?spm_id_from=333.337.search-card.all.click&vd_source=106e54383f94ae586c91138e171409f1
4.9 PCA分析的通用方法
在R语言中执行PCA有两种通用方法：
一种是Spectral
decomposition，它检查变量间的协方差/相关性；
另一种是Singular value
decomposition，它检查个体间的协方差/相关性。
4.9.1 主成分分析
4.9.1.1 可用函数
目前来讲，主要有5个包用于PCA的计算：
（1）stats默认加载包，提供procomp()与princomp()函数。
（2）psych包，提供principal()函数。
（3）FactoMineR包，提供PCA()函数。
（4）ade4包duli.pca()函数。
（5）vegan包的rda()函数。
不同的包计算出的PCA结果会有一点差别。
4.9.1.2 数据输入
可以输入原始数据矩阵或者相关系数矩阵到prcomp()函数中，若输入原始数据，相关系数矩阵将会自动计算，在计算前请确保数据中没有缺失值。下面的小节中将分别介绍这两类数据的分析方法。
4.9.1.3 计算主成分
1. 逐步拆解法
library(psych)数据预处理
train.scale<-scale(iris[,-5]) #`scale()`函数对数据进行标准化处理，即将每个变量减去其均值，再除以其标准差，以确保每个变量在相同的尺度上进行分析
nhl.cor <- cor(train.scale) #计算标准化后数据的相关系数矩阵。`cor()`函数用于计算相关系数。
cor.plot(nhl.cor) #绘制相关系数矩阵的图形
主成分分析
pc.iris1 <- principal(train.scale, rotate="none")
plot(pc.iris1$values, type="b", ylab="Eigenvalues", xlab="Component")
2. 输入原始数据的方法
pc.iris2 <- prcomp(iris[,1:4], scale = TRUE) #PCAattributes(pc.iris2) #查看pca分析输出对象的参数## $names
## [1] "sdev" "rotation" "center" "scale" "x"
##
## $class
## [1] "prcomp"3. 结果解读
sdev（Standard Deviation）：表示每个主成分的标准差。它反映了对应主成分方向上的数据变化程度。sdev越大，说明主成分方向上的数据变化越大，对总体数据的解释能力也就越强。
rotation：表示主成分与原始变量之间的线性关系。rotation矩阵描述了如何将原始变量投影到主成分空间中。每一列代表一个主成分，它们是原始变量的线性组合。
center：表示原始数据在每个变量上的均值。通过减去center的值，可以使得原始数据在各个维度上的均值为零。
scale：表示原始数据在每个变量上的标准差。通过除以scale的值，可以使得原始数据在各个维度上的标准差为1。
x：表示经过PCA转换后的数据。x矩阵的每一列代表一个主成分，它们是原始数据经过rotation变换后得到的新的特征变量。
————————–
根据attributes(pc.iris2)给出的各种参数信息，我们可以单独查看某一参数。
例如：
pc.iris2$sdev #返回主成分的标准差## [1] 1.7083611 0.9560494 0.3830886 0.1439265结果解读
第一个主成分的标准差为1.7083611。
第二个主成分的标准差为0.9560494。
第三个主成分的标准差为0.3830886。
第四个主成分的标准差为0.1439265。
标准差是衡量数据分散程度的一种指标，它描述了在该主成分方向上数据的变化程度。标准差越大，说明该主成分方向上的数据变化越大。在主成分分析中，标准差可以帮助我们理解每个主成分方向上对总体数据的贡献度大小，以及选择保留多少主成分的依据。通常来说，我们会选择标准差较大的主成分进行保留，因为它们能够很好地解释数据的变化。
————————–
pc.iris2$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该结果是主成分分析后，原始变量与每个主成分之间的相关系数（旋转矩阵）。具体解读如下：
第一个主成分（PC1）与原始变量之间的相关系数为：
Sepal.Length: 0.5210659
Sepal.Width: -0.2693474
Petal.Length:
0.5804131
Petal.Width: 0.5648565
第二个主成分（PC2）与原始变量之间的相关系数为：
Sepal.Length: -0.37741762
Sepal.Width: -0.92329566
Petal.Length: -0.02449161
Petal.Width: -0.06694199
第三个主成分（PC3）与原始变量之间的相关系数为：
Sepal.Length:
0.7195664
Sepal.Width: -0.2443818
Petal.Length: -0.1421264
Petal.Width: -0.6342727
第四个主成分（PC4）与原始变量之间的相关系数为：
Sepal.Length:
0.2612863
Sepal.Width: -0.1235096
Petal.Length: -0.8014492
Petal.Width: 0.5235971
这些相关系数表示了原始变量在每个主成分方向上的贡献程度。它们代表了原始变量与主成分之间的线性关系。相关系数的绝对值越大，表示原始变量在相应主成分中的贡献越大。这些相关系数可以帮助我们理解每个主成分所代表的特征或模式，并且可以用于解释主成分分析的结果。
————————–
pc.iris2$center #返回数据集各变量（列）的均值## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 5.843333 3.057333 3.758000 1.199333结果解读
Sepal.Length的均值为5.843333
Sepal.Width的均值为3.057333
Petal.Length的均值为3.758000
Petal.Width的均值为1.199333
————————–
pc.iris2$scale #返回标准化后的数据集各变量（列）的标准差## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.8280661 0.4358663 1.7652982 0.7622377pc.iris2$scale返回的结果，每个值表示对应变量在进行PCA分析时的标准差。
对于Sepal.Length来说，标准差为0.8280661。
对于Sepal.Width来说，标准差为0.4358663。
对于Petal.Length来说，标准差为1.7652982。
对于Petal.Width来说，标准差为0.7622377。
————————–
pc.iris2$x # 返回经过PCA转换后的数据，即主成分。## PC1 PC2 PC3 PC4
## [1,] -2.25714118 -0.478423832 0.127279624 0.024087508
## [2,] -2.07401302 0.671882687 0.233825517 0.102662845
## [3,] -2.35633511 0.340766425 -0.044053900 0.028282305
## [4,] -2.29170679 0.595399863 -0.090985297 -0.065735340
## [5,] -2.38186270 -0.644675659 -0.015685647 -0.035802870
## [6,] -2.06870061 -1.484205297 -0.026878250 0.006586116
## [7,] -2.43586845 -0.047485118 -0.334350297 -0.036652767
## [8,] -2.22539189 -0.222403002 0.088399352 -0.024529919
## [9,] -2.32684533 1.111603700 -0.144592465 -0.026769540
## [10,] -2.17703491 0.467447569 0.252918268 -0.039766068
## [11,] -2.15907699 -1.040205867 0.267784001 0.016675503
## [12,] -2.31836413 -0.132633999 -0.093446191 -0.133037725
## [13,] -2.21104370 0.726243183 0.230140246 0.002416941
## [14,] -2.62430902 0.958296347 -0.180192423 -0.019151375
## [15,] -2.19139921 -1.853846555 0.471322025 0.194081578
## [16,] -2.25466121 -2.677315230 -0.030424684 0.050365010
## [17,] -2.20021676 -1.478655729 0.005326251 0.188186988
## [18,] -2.18303613 -0.487206131 0.044067686 0.092779618
## [19,] -1.89223284 -1.400327567 0.373093377 0.060891973
## [20,] -2.33554476 -1.124083597 -0.132187626 -0.037630354
## [21,] -1.90793125 -0.407490576 0.419885937 0.010884821
## [22,] -2.19964383 -0.921035871 -0.159331502 0.059398340
## [23,] -2.76508142 -0.456813301 -0.331069982 0.019582826
## [24,] -1.81259716 -0.085272854 -0.034373442 0.150636353
## [25,] -2.21972701 -0.136796175 -0.117599566 -0.269238379
## [26,] -1.94532930 0.623529705 0.304620475 0.043416203
## [27,] -2.04430277 -0.241354991 -0.086075649 0.067454082
## [28,] -2.16133650 -0.525389422 0.206125707 0.010241084
## [29,] -2.13241965 -0.312172005 0.270244895 0.083977887
## [30,] -2.25769799 0.336604248 -0.068207276 -0.107918349
## [31,] -2.13297647 0.502856075 0.074757996 -0.048027970
## [32,] -1.82547925 -0.422280389 0.269564311 0.239069476
## [33,] -2.60621687 -1.787587272 -0.047070727 -0.228470534
## [34,] -2.43800983 -2.143546796 0.082392024 -0.048053409
## [35,] -2.10292986 0.458665270 0.169706329 0.028926042
## [36,] -2.20043723 0.205419224 0.224688852 0.168343905
## [37,] -2.03831765 -0.659349230 0.482919584 0.195702902
## [38,] -2.51889339 -0.590315163 -0.019370918 -0.136048774
## [39,] -2.42152026 0.901161067 -0.192609402 -0.009705907
## [40,] -2.16246625 -0.267981199 0.175296561 0.007023875
## [41,] -2.27884081 -0.440240541 -0.034778398 0.106626042
## [42,] -1.85191836 2.329610745 0.203552303 0.288896090
## [43,] -2.54511203 0.477501017 -0.304745527 -0.066379077
## [44,] -1.95788857 -0.470749613 -0.308567588 0.176501717
## [45,] -2.12992356 -1.138415464 -0.247604064 -0.150539117
## [46,] -2.06283361 0.708678586 0.063716370 0.139801160
## [47,] -2.37677076 -1.116688691 -0.057026813 -0.151722682
## [48,] -2.38638171 0.384957230 -0.139002234 -0.048671707
## [49,] -2.22200263 -0.994627669 0.180886792 -0.014878291
## [50,] -2.19647504 -0.009185585 0.152518539 0.049206884
## [51,] 1.09810244 -0.860091033 0.682300393 0.034717469
## [52,] 0.72889556 -0.592629362 0.093807452 0.004887251
## [53,] 1.23683580 -0.614239894 0.552157058 0.009391933
## [54,] 0.40612251 1.748546197 0.023024633 0.065549239
## [55,] 1.07188379 0.207725147 0.396925784 0.104387166
## [56,] 0.38738955 0.591302717 -0.123776885 -0.240027187
## [57,] 0.74403715 -0.770438272 -0.148472007 -0.077111455
## [58,] -0.48569562 1.846243998 -0.248432992 -0.040384912
## [59,] 0.92480346 -0.032118478 0.594178807 -0.029779844
## [60,] 0.01138804 1.030565784 -0.537100055 -0.028366154
## [61,] -0.10982834 2.645211115 0.046634215 0.013714785
## [62,] 0.43922201 0.063083852 -0.204389093 0.039992104
## [63,] 0.56023148 1.758832129 0.763214554 0.045578465
## [64,] 0.71715934 0.185602819 0.068429700 -0.164256922
## [65,] -0.03324333 0.437537419 -0.194282030 0.108684396
## [66,] 0.87248429 -0.507364239 0.501830204 0.104593326
## [67,] 0.34908221 0.195656268 -0.489234095 -0.190869932
## [68,] 0.15827980 0.789451008 0.301028700 -0.204612265
## [69,] 1.22100316 1.616827281 0.480693656 0.225145511
## [70,] 0.16436725 1.298259939 0.172260719 -0.051554138
## [71,] 0.73521959 -0.395247446 -0.614467782 -0.083006045
## [72,] 0.47469691 0.415926887 0.264067576 0.113189079
## [73,] 1.23005729 0.930209441 0.367182178 -0.009911322
## [74,] 0.63074514 0.414997441 0.290921638 -0.273304557
## [75,] 0.70031506 0.063200094 0.444537765 0.043313222
## [76,] 0.87135454 -0.249956017 0.471001057 0.101376117
## [77,] 1.25231375 0.076998069 0.724727099 0.039556002
## [78,] 1.35386953 -0.330205463 0.259955701 0.066604931
## [79,] 0.66258066 0.225173502 -0.085577197 -0.036318171
## [80,] -0.04012419 1.055183583 0.318506304 0.064571834
## [81,] 0.13035846 1.557055553 0.149482697 -0.009371129
## [82,] 0.02337438 1.567225244 0.240745761 -0.032663020
## [83,] 0.24073180 0.774661195 0.150707074 0.023572390
## [84,] 1.05755171 0.631726901 -0.104959762 -0.183354200
## [85,] 0.22323093 0.286812663 -0.663028512 -0.253977520
## [86,] 0.42770626 -0.842758920 -0.449129446 -0.109308985
## [87,] 1.04522645 -0.520308714 0.394464890 0.037084781
## [88,] 1.04104379 1.378371048 0.685997804 0.136378719
## [89,] 0.06935597 0.218770433 -0.290605718 -0.146653279
## [90,] 0.28253073 1.324886147 -0.089111491 0.008876070
## [91,] 0.27814596 1.116288852 -0.094172116 -0.269753497
## [92,] 0.62248441 -0.024839814 0.020412763 -0.147193289
## [93,] 0.33540673 0.985103828 0.198724011 0.006508757
## [94,] -0.36097409 2.012495825 -0.105467721 0.019505467
## [95,] 0.28762268 0.852873116 -0.130452657 -0.107043742
## [96,] 0.09105561 0.180587142 -0.128547696 -0.229191812
## [97,] 0.22695654 0.383634868 -0.155691572 -0.132163118
## [98,] 0.57446378 0.154356489 0.270743347 -0.019794366
## [99,] -0.44617230 1.538637456 -0.189765199 0.199278855
## [100,] 0.25587339 0.596852285 -0.091572385 -0.058426315
## [101,] 1.83841002 -0.867515056 -1.002044077 -0.049085303
## [102,] 1.15401555 0.696536401 -0.528389994 -0.040385459
## [103,] 2.19790361 -0.560133976 0.202236658 0.058986583
## [104,] 1.43534213 0.046830701 -0.163083761 -0.234982858
## [105,] 1.86157577 -0.294059697 -0.394307408 -0.016243853
## [106,] 2.74268509 -0.797736709 0.580364827 -0.101045973
## [107,] 0.36579225 1.556289178 -0.983598122 -0.132679346
## [108,] 2.29475181 -0.418663020 0.649530452 -0.237246445
## [109,] 1.99998633 0.709063226 0.392675073 -0.086221779
## [110,] 2.25223216 -1.914596301 -0.396224508 0.104488870
## [111,] 1.35962064 -0.690443405 -0.283661780 0.107500284
## [112,] 1.59732747 0.420292431 -0.023108991 0.058136869
## [113,] 1.87761053 -0.417849815 -0.026250468 0.145926073
## [114,] 1.25590769 1.158379741 -0.578311891 0.098826244
## [115,] 1.46274487 0.440794883 -1.000517746 0.274738504
## [116,] 1.58476820 -0.673986887 -0.636297054 0.191222383
## [117,] 1.46651849 -0.254768327 -0.037306280 -0.154811637
## [118,] 2.41822770 -2.548124795 0.127454475 -0.272892966
## [119,] 3.29964148 -0.017721580 0.700957033 0.045037725
## [120,] 1.25954707 1.701046715 0.266643612 -0.064963167
## [121,] 2.03091256 -0.907427443 -0.234015510 0.167390481
## [122,] 0.97471535 0.569855257 -0.825362161 0.027662914
## [123,] 2.88797650 -0.412259950 0.854558973 -0.126911337
## [124,] 1.32878064 0.480202496 0.005410239 0.139491837
## [125,] 1.69505530 -1.010536476 -0.297454114 -0.061437911
## [126,] 1.94780139 -1.004412720 0.418582432 -0.217609339
## [127,] 1.17118007 0.315338060 -0.129503907 0.125001677
## [128,] 1.01754169 -0.064131184 -0.336588365 -0.008625505
## [129,] 1.78237879 0.186735633 -0.269754304 0.030983849
## [130,] 1.85742501 -0.560413289 0.713244682 -0.207519953
## [131,] 2.42782030 -0.258418706 0.725386035 -0.017863520
## [132,] 2.29723178 -2.617554417 0.491826144 -0.210968943
## [133,] 1.85648383 0.177953334 -0.352966242 0.099675959
## [134,] 1.11042770 0.291944582 0.182875741 -0.185721512
## [135,] 1.19845835 0.808606364 0.164173760 -0.487849130
## [136,] 2.78942561 -0.853942542 0.541093785 0.294893130
## [137,] 1.57099294 -1.065013214 -0.942695700 0.035486875
## [138,] 1.34179696 -0.421020154 -0.180271551 -0.214702016
## [139,] 0.92173701 -0.017165594 -0.415434449 0.005220919
## [140,] 1.84586124 -0.673870645 0.012629804 0.194543500
## [141,] 2.00808316 -0.611835930 -0.426902678 0.246711805
## [142,] 1.89543421 -0.687273065 -0.129640697 0.468128374
## [143,] 1.15401555 0.696536401 -0.528389994 -0.040385459
## [144,] 2.03374499 -0.864624030 -0.337014969 0.045036251
## [145,] 1.99147547 -1.045665670 -0.630301866 0.213330527
## [146,] 1.86425786 -0.385674038 -0.255418178 0.387957152
## [147,] 1.55935649 0.893692855 0.026283300 0.219456899
## [148,] 1.51609145 -0.268170747 -0.179576781 0.118773236
## [149,] 1.36820418 -1.007877934 -0.930278721 0.026041407
## [150,] 0.95744849 0.024250427 -0.526485033 -0.162533529结果解读
x矩阵的每一列代表一个主成分，它们是原始数据经过rotation变换后得到的新的特征变量。
————————–
我们也可以使用print()函数和summary()查看pca对象的主要信息,例如：
print(pc.iris2) ## 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结果解读
首先，给出了四个标准的值，分别对应四个PC。与pc.iris2$sdev输出的结果相同。
然后是一个旋转矩阵（rotation matrix），与pc.iris2$rotation输出的结果相同。它描述了变量（Sepal.Length、Sepal.Width、Petal.Length和Petal.Width）与主成分之间的关系。每个主成分都由一组权重（PC1, PC2, PC3, PC4）来表示，这些权重表明了变量对于该主成分的贡献程度。
————————–
summary(pc.iris2)## 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结果解读
1）“Standard
deviation”表示PC1、PC2、PC3和PC4的标准差分别为1.7084、0.9560、0.38309和0.14393。
2）“Proportion of
Variance”表示每个主成分解释的方差比例。在这里，PC1解释了数据方差的大约72.96%，PC2解释了22.85%，PC3解释了3.67%，PC4解释了0.52%。这意味着PC1包含了最多的信息，其次是PC2，依此类推。
3）“Cumulative
Proportion”给出了累积方差比例。它表示前n个主成分解释的方差比例之和。在这里，前两个主成分（PC1和PC2）累计解释了约95.81%的方差，前三个主成分累计解释了约99.48%，而所有四个主成分累计解释了100%的方差。换句话说，PC1解释了变异的72.96%，PC1+PC2解释了变异的95.81%，以此类推。
————————–
4.9.2 判断主成分的个数
判断PCA中需要多少个主成分的一般准则：
1）根据先验经验和理论知识判断；
2）根据变量方差积累值的阈值判断；
3）通过检验变量k✖️k的相关系数矩阵来判断保留的主成分。
最常见的是基于特征值的方法,PC1与最大的特征值相关，PC2与第二大特征值相关，依次类推。
Kaiser-Harri准则建议保留特征值>1的主成分，一般特征值小于1的成分所解释的方差比包含在单个变量中的方差更小。注意这只是一般准则，不是绝对准则。
你还可以进行模拟，依据与初始矩阵相同大小的随机数据矩阵来判断要提取的特征值，若基于真实数据的某个特征值大于一组随机数据矩阵相应的平均特征值，那么主成分可以保留，该方法称为平行分析。
根据变异方差积累值，并没有标准准则，这个值代表了你能解释变异的程度，当然是能解释的越多越好，比如你解释了80%的变异比解释了70%要好，这分别代表了数据的20%和30%的变异你没有解释清楚。
4.9.3 可视化
1）碎石图：显示每个主成分解释的方差百分比。
library(ggplot2)
library(factoextra)fviz_eig(pc.iris2)
————————–
2）Graph of individuals(样本图)
具有相似特征的样本被分组在一起。
library(gridExtra)
#总体样本
fviz_pca_ind(pc.iris2,
col.ind = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = FALSE
)
代码解读
col.ind = “cos2”,
表示使用个体在主成分上的贡献度(原始数据中的各个样本individual)来对散点图的颜色进行着色，以反映各个个体在主成分上的解释度量。
————————–
3）Graph of variables（变量图)
正相关变量指向图的同一侧。负相关变量指向图表的相反两侧。
fviz_pca_var(pc.iris2,
col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = FALSE # Avoid text overlapping
)
代码解读
col.var = “contrib”
是使用变量在主成分上的贡献度来（原始数据中各个变量）来对散点图的颜色进行着色，以反应变量中主成分上的解释度量。
————————–
4）Biplot of individuals and variables
方法一：fviz_pca_biplot()函数
fviz_pca_biplot(pc.iris2, repel = FALSE,
col.ind = "cos2", # Individuals color
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
col.var = "#00AFBB" # Variables color
)
方法二：ggbiplot()函数
library(devtools)install_github("vqv/ggbiplot")
# 需要先library(devtools)，才能进行安装，安装过程中会提示你是否需要更新，可以选择3（不更新）library(ggbiplot)
ggbiplot(pc.iris2,
obs.scale = 1,
var.scale = 1,
groups = iris$Species,
ellipse = TRUE,
circle = TRUE,
ellipse.prob = 0.75) +
scale_color_discrete(name = 'Species') +
theme(legend.direction = 'vertical',
legend.position = 'right')
代码解读
使用了ggbiplot库来可视化主成分分析的结果：
obs.scale = 1和var.scale = 1表示对观测值和变量进行比例尺标准化。
groups = iris$Species将观测点按照物种标签进行分组。
ellipse = TRUE和circle = TRUE表示在图中绘制椭圆和圆形。
ellipse.prob = 0.75表示正太概率中椭圆的大小，可以尝试更改数值感受图的变化，取值0-1之间。
scale_color_discrete(name = ’’)：这一行代码用于设置颜色比例尺的名称为Species。
theme(legend.direction = ‘horizontal’, legend.position = ‘top’)：用于设置图例的方向为垂直方向，位置位于右侧。
结果解读
1)
PC1解释了大约73%的变异性，PC2解释了约22.9%的变异性。
箭头之间越接近表示相关性越高。例如，萼片宽度与其他变量之间的相关性很弱。
解读这个图的另一种方式是，PC1与花瓣长度、花瓣宽度和萼片长度呈正相关，而与萼片宽度呈负相关。PC2与萼片宽度呈高度负相关。
Biplot在主成分分析中是一个重要工具，可以理解数据集中发生了什么。
————————–
4.10 PCA分析进阶玩法
4.10.1 主成分的正交性
在上一节中，我们讲到对于具有多重共线性情况的数据，直接基于原始数据集进行模型预测可能是错误的。那么基于主成分数据pc.iris2，看看多重共线性问题是否得到解决？
library(psych)
pairs.panels(pc.iris2$x,
gap=0,
bg = c("red", "yellow", "blue")[iris$Species],
pch=21)
代码解读
pairs.panels()函数被用于可视化主成分分析的结果。
pc.iris2$x是主成分分析的结果，其中包含了转换后的主成分数据。
结果解读
现在相关系数为零，因此我们可以摆脱多重共线性问题。
————————–
4.10.2 主成分预测
前面我们知道使用主成分数据（pc.iris2）可以解决多重共线性问题，因为我们可以用找个数据进行下一步的预测模型分析。
首先，我们需要把将iris数据集的前四列的原始数据转换为主成分分析的结果，操作如下：
pd <- predict(pc.iris2, iris) #predict()函数将原始数据转换为主成分分析的结果。
pd <- data.frame(pd, iris[5]) #将物种列信息添加到 pd 中因为我们的因变量factor有3个水平(即3个种)，所以我们将使用多元逻辑回归模型（Multinomial Logistic Regression Model）。
#install.packages("nnet")
library(nnet)
pd$Species <- relevel(pd$Species, ref = "setosa")#把"setosa"设置为参考类别
mymodel <- multinom(Species~PC1+PC2, data = pd)## # weights: 12 (6 variable)
## initial value 164.791843
## iter 10 value 25.921971
## iter 20 value 24.637111
## iter 30 value 24.515626
## iter 40 value 24.513212
## iter 40 value 24.513212
## final value 24.513212
## converged代码解读
multinom()函数对训练数据进行多元逻辑回归模型的拟合。
Species~PC1+PC2表示使用PC1和PC2来预测三个种（Species）。
结果解读
首先，代码对数据集中的”Species”变量进行了重新水平化（relevel）操作，将”setosa”设置为参考类别。
接下来，通过使用multinom函数构建了一个多项逻辑回归模型（mymodel），其中自变量为”PC1”和”PC2”，因变量为”Species”。
在模型拟合的过程中，迭代进行了多次，最终收敛到一个最优解。
converged意味着模型已经找到了最佳参数配置，可以用于预测新的数据。这个模型可以用来根据给定的”PC1”和”PC2”值预测鸢尾花的种类，并且通过比较这两个主成分的权重，可以了解它们对种分类的贡献程度。
————————–
summary(mymodel)## Call:
## multinom(formula = Species ~ PC1 + PC2, data = pd)
##
## Coefficients:
## (Intercept) PC1 PC2
## versicolor 9.125730 12.77516 2.853780
## virginica 2.728464 18.56054 3.396886
##
## Std. Errors:
## (Intercept) PC1 PC2
## versicolor 108.4973 67.86478 82.44931
## virginica 108.5069 67.87592 82.45103
##
## Residual Deviance: 49.02642
## AIC: 61.02642结果解读
系数（Coefficients）部分给出了每个预测变量的回归系数。在这里，拟合出两个种（versicolor和virginica）的回归系数。例如，对于versicolor种，回归方程为：Species = 9.125730 + 12.77516 * PC1 + 2.853780 * PC2。类似地，对于virginica种，回归方程为：Species = 2.728464 + 18.56054 * PC1 + 3.396886 * PC2。这些回归系数用于计算特定PC1和PC2值时，归类到每个种的概率。
标准误差（Std. Errors）部分给出了每个回归系数的标准误差。标准误差用于衡量回归系数的稳定性。
残差（Residual Deviance）是模型的拟合优度指标，度量了模型与观测数据之间的偏差。较小的残差值表示模型更好地拟合了数据。
赤池信息量准则 (AIC, Akaike information criterion)，是衡量统计模型拟合优良性的一种标准，是由日本统计学家赤池弘次创立和发展的。赤池信息量准则建立在熵的概念基础上。 AIC越小，模型越好，通常选择AIC最小的模型.
综上所述，这些结果提供了有关多项逻辑回归模型的参数估计、标准误差、拟合优度和模型选择的信息。可以使用这些结果来解释主成分对鸢尾花种类预测的影响以及模型的准确性。
下面给出了模型，我们只使用了前两个主要成分，因为大多数信息都在第一个成分中可用。
————————–
4.10.3 混淆矩阵和错误分类错误
让我们看看数据集中的正确分类和错误分类。
pd2 <- predict(mymodel, pd)
tab <- table(pd2, iris$Species)
tab##
## pd2 setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 43 5
## virginica 0 7 45代码解读
pd2 <- predict(mymodel,
pd)是使用训练好的PCA模型 mymodel 对数据集 iris
进行预测。这将为每个样本在主成分空间中分配一个预测标签。
tab <-
table(pd2, iris\(Species) 则是根据预测标签 pd2
和实际类别 iris\)Species 创建了一个频数表（contingency
table）tab，用于展示预测结果和实际类别之间的对应关系。
结果解读
在预测结果中，有50个样本被分类为setosa类别，并且都正确预测；
在预测结果中，有43个样本被分类为versicolor类别，并且都正确预测；
在预测结果中，有45个样本被分类为virginica类别，并且都正确预测；
在预测结果中，有37个样本被错误地分类为virginica类别，它们实际属于versicolor类别；
在预测结果中，有5个样本被错误地分类为versicolor类别，它们实际属于virginica类别。
计算Misclassification error
1 - sum(diag(tab))/sum(tab)## [1] 0.08代码解读
通过计算预测结果与实际类别不一致的样本数占总样本数的比例，可以评估模型的错误分类率（misclassification
rate）。 结果解读
错误分类误差为
8%，即约有8%的样本被错误地分类到了其他类别。
4.10.4 PCA绘图进阶
library(ggplot2)
library(ggpubr)
library(cowplot)步骤1:自定义绘图函数
函数代码源自：Molania, Ramyar, et al. “Removing unwanted variation from large-scale RNA sequencing data with PRPS.” Nature Biotechnology 41.1 (2023): 82-95.
此处自定义设置了.pca()和.scatter.density.pc()函数，只需要复制运行即可。
.pca <- function(data, is.log) {
if (is.log)
data <- data
else
data <- log2(data + 1)
svd <- base::svd(scale(
x = t(data),
center = TRUE,
scale = FALSE
))
percent <- svd$d ^ 2 / sum(svd$d ^ 2) * 100
percent <- sapply(
seq_along(percent),
function(i) round(percent[i], 1))
return(
list(
sing.val = svd,
variation = percent)
)
}
.scatter.density.pc <- function(
pcs,
pc.var,
group.name,
group,
color,
strokeSize,
pointSize,
strokeColor,
alpha
){
pair.pcs <- utils::combn(ncol(pcs), 2)
pList <- list()
for(i in 1:ncol(pair.pcs)){
if(i == 1){
x <- pair.pcs[1,i]
y <- pair.pcs[2,i]
p <- ggplot(mapping = aes(
x = pcs[,x],
y = pcs[,y],
fill = group)) +
xlab(paste0('PC', x, ' (', pc.var[x], '%)')) +
ylab(paste0('PC', y, ' (', pc.var[y], '%)')) +
geom_point(
aes(fill = group),
pch = 21,
color = strokeColor,
stroke = strokeSize,
size = pointSize,
alpha = alpha) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 5)) +
scale_y_continuous(breaks = scales::pretty_breaks(n = 5)) +
theme(
legend.position = "right",
panel.background = element_blank(),
axis.line = element_line(colour = "black", size = 1.1),
legend.background = element_blank(),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14),
legend.key = element_blank(),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 10),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14)) +
guides(fill = guide_legend(override.aes = list(size = 4))) +
scale_fill_manual(name = group.name, values = color)
le <- ggpubr::get_legend(p)
}else{
x <- pair.pcs[1,i]
y <- pair.pcs[2,i]
p <- ggplot(mapping = aes(
x = pcs[,x],
y = pcs[,y],
fill = group)) +
xlab(paste0('PC', x, ' (',pc.var[x], '%)')) +
ylab(paste0('PC', y, ' (',pc.var[y], '%)')) +
geom_point(
aes(fill = group),
pch = 21,
color = strokeColor,
stroke = strokeSize,
size = pointSize,
alpha = alpha) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 5)) +
scale_y_continuous(breaks = scales::pretty_breaks(n = 5)) +
theme(
panel.background = element_blank(),
axis.line = element_line(colour = "black", size = 1.1),
legend.position = "none",
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 10),
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14)) +
scale_fill_manual(values = color, name = group.name)
}
p <- p + theme(legend.position = "none")
xdens <- cowplot::axis_canvas(p, axis = "x")+
geom_density(
mapping = aes(
x = pcs[,x],
fill = group),
alpha = 0.7,
size = 0.2
) +
theme(legend.position = "none") +
scale_fill_manual(values = color)
ydens <- cowplot::axis_canvas(
p,
axis = "y",
coord_flip = TRUE) +
geom_density(
mapping = aes(
x = pcs[,y],
fill = group),
alpha = 0.7,
size = 0.2) +
theme(legend.position = "none") +
scale_fill_manual(name = group.name, values = color) +
coord_flip()
p1 <- insert_xaxis_grob(
p,
xdens,
grid::unit(.2, "null"),
position = "top"
)
p2 <- insert_yaxis_grob(
p1,
ydens,
grid::unit(.2, "null"),
position = "right"
)
pList[[i]] <- ggdraw(p2)
}
pList[[i+1]] <- le
return(pList)
}步骤2:PCA高级绘图
#PCA分析
pca.ncg <- .pca(data = iris[,1:4], is.log = FALSE)
#使用函数 .scatter.density.pc 生成一个散点图,并赋予给ph。
ph <- .scatter.density.pc(pcs = pca.ncg$sing.val$v[,1:3],
pc.var =pca.ncg$variation,
strokeColor = 'gray30',
strokeSize = .2,
pointSize = 2,
alpha = .6,
group.name = "Species",
group = iris$Species,
color = c("red","blue","green"))
#可视化
do.call(gridExtra::grid.arrange, c(ph, ncol= 2))