灰色预测模型代码_生信审稿人最常问的验证！临床预测模型中的PCA主成分分析！这点你注意到了没！（附代码）...

在高维度数据不断出现的现在，输入特征的数量呈指数形式迅速增长，形成了所谓的维数灾难。今天，我们将重点讲讲主成分分析，从原始数据集中找出一个更小的，但能最大程度保留原来大部分信息的变量集合。

所谓主成分分析，即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()函数，将输入基因之间的相关性进行可视化展示。

2.3 读取临床数据

cli <- read.table("cliData.txt",header=T,sep="  ",check.names=F, row.names = 1) 
head(cli) 
str(cli)

随后，将提前清洗准备好的临床数据进行读取。

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个主成分是比较令人信服的。

3.2 正交旋转与解释

正交旋转又称为“方差最大法”，而旋转的意义是使变量在某个主成分上的载荷最大化，减少主成分之间的相关性，从而有助于对主成分的解释。

pca.rotate <- principal(exp.scale, nfactors = 4, rotate = "varimax")
pca.rotate

在输出的结果中，4个主成分的变量载荷分别标注为RC1到RC4，而SS loading中的值是每个主成分的特征值。同时，Cumulative Var结果表示这4个旋转后的主成分可以解释86%的全部方差。一般而言，我们选择的主成分至少应该解释大于70%的全部方差。可见，4个主成分已经可以用于后续的建模分析。

3.3 根据主成分建立因子得分

pca.scores <- data.frame(pca.rotate$scores)
head(pca.scores)

检查旋转和的主成分载荷，并将其作为每个患者的因子得分。而且，这些得分说明了每个患者与旋转后的主成分之间的相关程度。

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具有显著性。

#使用RC1和RC4变量建立线性模型 
lm_2 <- lm(fustat ~ RC1+RC4, data = exp.pca) 
summary(lm_2)

随后，将不显著的主成分因子排出模型，只保留RC1和RC4，再一次构建模型。可以发现，整体模型的统计学显著性得到了提升，P值为0.0057。

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。

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，表明其模型的预测能力较为一般。

好了，PCA的相关内容就介绍到这里了，简单的讲呢，就是使用一些相关的变量构造一个新的变量。在今后的文章分析中，遇到其他模型不理想时，或者其他高维度数据时候，不妨试试主成分分析来进行分析~~~

欢迎大家关注解螺旋生信频道-挑圈联靠公号~

—END—

本文首发于“挑圈联考”微信公众号

转载请注明：解螺旋·临床医生科研成长平台