成分分析（Principal Component Analysis，PCA）， 将多个变量通过线性变换以选出较少个数重要变量的一种多元统计分析方法。又称主分量分析。当在研究中， 经常会遇到多变量的情况。那么问题来了，如果在分析数据的时候，每个变量都去加以分析那必定会使我们的数据分析变得很困难，而且加大了分析难度和工作量。但是如果贸然的删除变量的话，必定会影响结果的有效性和准确性。那遇到这种问题我们该怎么办呢？？PCA就派上用场了，一般我们取第一主成分（PC1）第二主成分（PC2）就能解决问题了，有必要的时候可以加上PC3。
还有，如果研究对象是生物群体的话，会遇到具有群体结构的情况，这样如果不考虑群体结构的话直接分析，结果会不准确或有效性不高，所以面对此类群体，首先要查看群体分层情况，看看两个群体或者多个群体之间是否存在群体分层。此时PCA会起到关键性作用，可以用来检测群体之间是否存在遗传差异和群体分层。
好了闲话少说，开始正经内容。
>setwd("C:/Users/gaom/Desktop")#打开文件所在路径，并将文件所在目录作为工作目录。
>data<-read.table(file= "test_data.txt",header= T,sep = "\t")#读取数据。
用r语言做pca分析有连个函数可以调用，princomp函数可以从相关矩阵或者从协方差矩阵做主成分分析；prcomp 函数可以从奇异值分解的算法做主成分分析。上述两个函数都可以做主成分分析，不过相对的函数，需要的数据格式也不相同。
其中R-mode简单来说就是数据数量要大于变量（因素）数量；而Q-mode正与其相反，简单来讲就是变量数量要多余数据数量。
OK，了解了上述函数，接着往下进行(首先介绍princomp函数)
>pca=princomp(data,cor=TRUE) #调用函数--princomp，#当cor=TRUE表示用样本的相关矩阵R做主成分分析当cor=FALSE表示用样本的协方差阵S做主成分分析。
>summary(pca)
>pc=loadings(pca)
>PCA1 = pca$loadings[,1]
>PCA2 = pca$loadings[,2]
>PCA3 = pca$loadings[,3]
接着再介绍一下（prcomp函数的用法）
>data1=t(as.matrix(data)) #由于两个函数对数据的模型不同所以需要转置一下。
>pca =prcomp(data1, scale= TRUE) #scale= TRUE对数据进行标准化
>summary(pca）
PCA1=pca$x[,1]
PCA2=pca$x[,2]
PCA3=pca$x[,3]
上述两个函数都能够得到pca1，pca2，和pca3。接下来画图
首先可以利用pca1和pca2画出二维的图。
>plot(PCA1,PCA2)
>group = factor(c(rep("HIGH",500 ), rep("LOW",400))) #对群体进行分组，我这里就是high组有500个体，low组有400个体。
>colour_group<-rainbow(length(unique(group)))
>colour<-colour_group[as.numeric(factor(group))]
>plot(PCA1,PCA2,col = colour ,pch = c(20,20)[group],cex=2,main="PCA-PLOT")
title("test") #这样就完成了二维pca的图，可以自己调整颜色和形状。
还可以用ggplot2这个r语言的package来画图，相对来讲ggplo2画出的图片要美观一些。
>library(ggplot2)
>group2<-data.frame(group)
>pca_reuslt<-as.data.frame(pca$x)
>pca_reuslt<-cbind(pca_reuslt,group2)
>p=ggplot(pca_reuslt)+geom_point(aes(x=pca_reuslt[,1],y=pca_reuslt[,2],color=pca_reuslt$group,shape= pca_reuslt$group),size=5)
>p<-p+theme(legend.title=element_blank())+labs(x="PCA1",y="PCA2")
>p
这样就画好了pca的分布图了，可以直观的看出群体之间的分层状况。
要是pca1，pca2，不足以满足需求，比如二维图不能分开某些群体，或者群体之间的分层不明显。就可以利用pca3，画出一个空间坐标轴，做出一个3d的图片，这3d图片可以利用plot3D这个package来画。
>library(plot3D)
>library(scatterplot3d)
>par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
>scatterplot3d(pca_reuslt[,1:3],pch=20 ,color=colour,angle=45,main= "test_3D",cex.symbols= 2,mar=c(5.1, 4.1, 4.1, 8.1))
>legend("right",legend= group,col =colour,pch = 20,bg="white",xpd=TRUE,inset= -0.5)
这样3D的图片就画好了。3D图片不仅看起来更加直观而且能分开之前分层不明显的群体。
希望此篇自学笔记，能够给需要用r语言做pca的同志们带来一点点的帮助。以上的代码均为自己手打，可能会出现一些标点的错误，使用时注意。还有以上纯属自己的个人理解和看法，以及做法，如有不规范的或者错误的地方，望大神们指正。
网友评论