在高维度数据不断出现的现在,输入特征的数量呈指数形式迅速增长,形成了所谓的维数灾难。今天,我们将重点讲讲主成分分析,从原始数据集中找出一个更小的,但能最大程度保留原来大部分信息的变量集合。
所谓主成分分析,即Principle Component Analysis (PCA),其实就是寻找主成分的过程,且通常被认为是一种特殊的非监督学习算法,可以对复杂或多变的变量进行处理分析。一般而言,成分可以认为是特征的规范化线性组合。其中,第一主成分是能够最大程度解释数据中方差的特征线性组合,而第二主成分是在方向上与第一主成分垂直的条件下,最大程度解释数据中方差的另一种线性组合。其后的每一个主成分都遵循相同的原则。此外,主成分得分是对于每个主成分和每个观测变量计算而得到的。
下面,我们一起来看下如何进行主成分分析并建立模型。
1.R包的安装与读取
1.R包的安装与读取
rm(list = ls()) #清空环境变量
options(stringsAsFactors = F)
###1. R包的安装与读取
if(!require(psych))install.packages("psych")
if(!require(ggplot2))install.packages("ggplot2")
if(!require(ggpubr))install.packages("ggpubr")
if(!require(ROCR))install.packages("ROCR")
library(psych) #PCA packages
library(ggplot2)
library(ggpubr)
library(ROCR)
psych包常被用于进行主成分分析,包括:判断主成分的个数;提取主成分;获取主成分得分;以及列出主成分方程,解释主成分意义。
2.数据的读取与预处理
2.1 读取表达数据
rt <- read.table("exp.txt", header=T, sep=" ", check.names=F, row.names = 1)
exp <- t(log2(rt+1))
str(exp)
与之前的模型构建一样,首先使用read.table()函数读取表达数据,并通过结构函数str()检查数据。
2.2 数据标准化
exp.scale <- scale(exp)
exp.cor <- cor(exp.scale)
cor.plot(exp.cor)
在分析之前,我们需要对数据进行标准化,使数据的均值为0,标准差为1。同时,完成标准化后,使用psych包提供的cor.plot()函数,将输入基因之间的相关性进行可视化展示。
![e8d61c282155e3ab3df361c1e00773e2.png](https://i-blog.csdnimg.cn/blog_migrate/3ccb75381c9856d7df58f14b5107c209.jpeg)
2.3 读取临床数据
cli <- read.table("cliData.txt",header=T,sep=" ",check.names=F, row.names = 1)
head(cli)
str(cli)
随后,将提前清洗准备好的临床数据进行读取。
![f8b61f130096144ddeeadd97729a53ed.png](https://i-blog.csdnimg.cn/blog_migrate/f41b49f891a890715e8a62b7606b44bb.png)
3.主成分分析
3.1 主成分抽取
pca <- principal(exp.scale, rotate = "none")
通过psych包的principal()函数抽取主成分。首先,我们将rotate的参数设置为none,得到其中的各个成分。
plot(pca$values, type = "b", ylab = "Eigenvalues", xlab = "Component")
接着,为了确定要保留的成分的数量,使用碎石图来评估能解释大部分数据方差的主成分,其中x轴表示主成分的数量,y轴表示相应的特征。在碎石图中,需要找出使变化率降低的那个点,也就是常说的“肘点”或者弯曲点。从图中可以看出,4个主成分是比较令人信服的。
![5b4e6e98695abb134e7686e6a6c87026.png](https://i-blog.csdnimg.cn/blog_migrate/b6459e5b7f3bdbf1e19b7b53d2c953b9.png)
3.2 正交旋转与解释
正交旋转又称为“方差最大法”,而旋转的意义是使变量在某个主成分上的载荷最大化,减少主成分之间的相关性,从而有助于对主成分的解释。
pca.rotate <- principal(exp.scale, nfactors = 4, rotate = "varimax")
pca.rotate
在输出的结果中,4个主成分的变量载荷分别标注为RC1到RC4,而SS loading中的值是每个主成分的特征值。同时,Cumulative Var结果表示这4个旋转后的主成分可以解释86%的全部方差。一般而言,我们选择的主成分至少应该解释大于70%的全部方差。可见,4个主成分已经可以用于后续的建模分析。
![2a97f0d84def2f410567562085af3192.png](https://i-blog.csdnimg.cn/blog_migrate/2c06ff2f86bea343d9934a3e353f5e2a.jpeg)
3.3 根据主成分建立因子得分
pca.scores <- data.frame(pca.rotate$scores)
head(pca.scores)
检查旋转和的主成分载荷,并将其作为每个患者的因子得分。而且,这些得分说明了每个患者与旋转后的主成分之间的相关程度。
![7bdcd38e6cd54cacbf84d0276cd14562.png](https://i-blog.csdnimg.cn/blog_migrate/bbb8a2408382d14d618a3d9aa3138fd2.jpeg)
exp.pca <-cbind(cli, pca.scores)
最后,将得分与临床信息进行合并。
4.模型构建与评价
4.1 回归分析
#通过lm()函数建立线性模型
lm_all <- lm(fustat ~ RC1+RC2+RC3+RC4, data = exp.pca)
summary(lm_all)
首先,将所有的因子作为输入,然后查看结果摘要。结果显示,整体模型在统计学上是具有显著性的,P值为0.0348。而且,其中RC1和RC4具有显著性。
![027ed1cde703a96e253a417f423d7def.png](https://i-blog.csdnimg.cn/blog_migrate/d1f9e86b6ababfdf75878848ea8107ad.jpeg)
#使用RC1和RC4变量建立线性模型
lm_2 <- lm(fustat ~ RC1+RC4, data = exp.pca)
summary(lm_2)
随后,将不显著的主成分因子排出模型,只保留RC1和RC4,再一次构建模型。可以发现,整体模型的统计学显著性得到了提升,P值为0.0057。
![68d71715ffce8a5c25072d9909fdbd4c.png](https://i-blog.csdnimg.cn/blog_migrate/e951e23e5c44ef349abeb08391ec9a91.jpeg)
4.2 模型评价
exp.pca$predict <- lm_2$fitted.values
ggboxplot(exp.pca,
x = "fustat", y = "predict",
color = "fustat", palette = "jco",
add = "jitter") +
stat_compare_means()
提取模型得到的预测值,且箱型图结果显示两组之间具有显著性差异,P值为0.0035。
![faab530569efbf761233b1e527464970.png](https://i-blog.csdnimg.cn/blog_migrate/165325871a2d74e6622f44173b6b8f97.png)
4.3 ROC分析
pred <- prediction(exp.pca$predict, exp.pca$fustat)
ROC <- performance(pred,"tpr","fpr")
auc <- performance(pred,"auc")@y.values[[1]]
auc
#[1] 0.7101449
plot(ROC,
col="red", #曲线的颜色
xlab="False positive rate", ylab="True positive rate", #x轴和y轴的名称
lty=1,lwd=3,
main=paste("AUC=", auc))
abline(0, 1, lty=2, lwd=3) #绘制对角线
同时,ROC曲线分析结果显示,AUC值为0.71,表明其模型的预测能力较为一般。
![f9e93c29c214a18bdd7b458179baed5c.png](https://i-blog.csdnimg.cn/blog_migrate/b1ea19d0481e0be31a423ef5530a957a.png)
好了,PCA的相关内容就介绍到这里了,简单的讲呢,就是使用一些相关的变量构造一个新的变量。在今后的文章分析中,遇到其他模型不理想时,或者其他高维度数据时候,不妨试试主成分分析来进行分析~~~
![f3850f4c52f8485fc0bef6dd7695325e.png](https://i-blog.csdnimg.cn/blog_migrate/be14c4d7e78aa56b1ad4f2c952beeddf.jpeg)
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
![d04ddeb88d03f6f501e51b7390bc077b.png](https://i-blog.csdnimg.cn/blog_migrate/ef3eeb419eb34e17e44e7727f1e54cbd.png)
本文首发于“挑圈联考”微信公众号
转载请注明:解螺旋·临床医生科研成长平台