PCA主成分分析_R语言实战

作为零计算机基础，纯生物背景的实验喵，有很多技能还不会，只希望将自己学习到的知识记录下来，一方面让自己能够时常温故知新，一方面与大家分享学习内容和心得，一起进步呀。

主成分分析（principle component analysis，PCA）是常用的线性数据降维方法，通过某种线性投影，将高维的数据映射到低维的空间中，并期望在映射到的低维空间中包含的信息量最大，即使用较少的数据维度来代表（保留）较多原始数据的特性。话句话说，就是把多个变量化成少数几个主成分的方法，这些主成分能够反映原始变量的大部分信息，为原始变量的线性组合。

R中作为主成分分析最主要的函数是 princomp() 函数

princomp() 主成分分析   可以从相关阵或者从协方差阵做主成分分析

summary() 提取主成分信息

loadings() 显示主成分分析或因子分析中载荷的内容

predict() 预测主成分的值

screeplot() 画出主成分的碎石图

plot() 画出数据关于主成分的散点图和原坐标在主成分下的方向

##设置工作路径

setwd("E:/R practice result")

##查看工作路径

getwd()

##读取数据

data <- readcsv("E:/R practice file/students_datacsv",header = T)
##进行主成分分析

datapr <- princomp(data,cor = T) ##cor=T的意思是用相关系数进行主成分分析
##观察主成分的详细情况

summary(datapr,loadings = T)
##计算得到各个样本主成分的数据

pca_data <- predict(datapr)

##将文件写出，进行本地保存

writetable(pca_data,file = "E:/R practice result/pca_datatxt",rownames = T,colnames = T,quote = F,sep = "\t")
#画出主成分碎石图，主成分方差占总体方差的情况,辅助判断主成分提取个数。当折线由陡峭突然变得平稳时，陡峭到平稳对应的主成分个数即为参考提取主成分个数。X轴表示可能的因子数。在整个曲线下降的过程中，理论上会存在一点，该点处曲线出现明显的弯折，该点对应的X轴的数字即应保留的因子数。

screeplot(datapr,type="lines")

#添加划分成分的参考线

abline (v = 2 ,   col = "#900021", lwd = 2, lty = 2)
#计算前两个主成分Comp1和Comp2得分,保存到datapr2

data_pr2 <- datapr$score[,1:2]    或者data_pr2 <- pca_data[,1:2]

writetable(data_pr2,file = "E:/R practice result/data_pr2txt",rownames = F,colnames = T,quote = F,sep = "\t")
plot(data_pr2, main = "after PCA") 或者plot(pca_data[,1:2])

如果说数据挖掘是一门手艺，那么R语言就是工匠手里一种工具，要做出一件价值连城的艺术品需要先“利其器”，但更关键的是工匠能够“集百家之长，成一家之言“自成一派的创造力。
正所谓” *** 千曲而后晓声，观千剑而后识器"。建议初学者先了解一些机器学习的基础理论以及典型的应用领域实例，定下自己想要研究的方向后与行业相结合，然后再学习工具的使用。
R语言经典图书推荐：
《R in Action-Data Analysis and Graphics with R》链接：R语言实战 (豆瓣)
这本书从实用的统计研究角度，每一章节结合实际的例子讲解了R在创建数据集、绘制图形、数据管理、以及模型构建的使用方法，堪称经典。前两部分属于R基本功能介绍，第三部分以后才是精髓开始（包括了回归分析、方差分析、功效分析、广义线性模型、主成分和因子分析等统计方法详细的实例分析）。
《数据之魅-基于开源工具的数据分析》链接：数据之魅 (豆瓣)
作者是华盛顿大学理论物理学博士。这本书是数据分析系列著作的经典之一，包含大量的R语言模拟过程及结果展示，例举了很多数据分析实例和代码。吃透以后就能够对整个数据挖掘的流程有一个全方位的了解。
转载

主成分分析 因子分析 典型相关分析，三种方法的共同点主要是用来对数据降维处理的 从数据中提取某些公共部分，然后对这 些公共部分进行分析和处理。
#主成分分析 是将多指标化为少数几个综合指标的一种统计分析方法
主成分分析是一种通过降维技术把多个变量化成少数几个主成分的方法，这些主成分能够反映原始变量的大部分信息，他们通常表示为原始变量的线性组合。

R的功能很强大，各种包很多。但就是因为包太多，造成了很大的麻烦。不可避免的，可以做结构方程模型的包也不少，例如：sem、psych、OpenMx，lavaan等。我选择了lavaan包。原因：语法简介易懂，上手快，支持非正态、连续数据，可以处理缺失值。
lavaan包是由比利时根特大学的Yves Rosseel开发的。lavaan的命名来自于 latent variable analysis，由每个单词的前两个字母组成，la-va-an——lavaan。
为什么说它简单呢 主要是因为它的lavaan model syntax，如果你会R的回归分析，那它对你来说再简单不过了。
一、语法简介
语法一：f3~f1+f2（路径模型）
结构方程模型的路径部分可以看作是一个回归方程。而在R中，回归方程可以表示为y~ax1+bx2+c，“~”的左边的因变量，右边是自变量，“+”把多个自变量组合在一起。那么把y看作是内生潜变量，把x看作是外生潜变量，略去截距，就构成了lavaan model syntax的语法一。
语法二：f1 =~ item1 + item2 + item3（测量模型）
"=~"的左边是潜变量，右边是观测变量，整句理解为潜变量f1由观测变量item1、item2和item3表现。
语法三：item1 ~~ item1 , item1 ~~ item2
"~~"的两边相同，表示该变量的方差，不同的话表示两者的协方差
语法四：f1 ~ 1
表示截距
此外还有其它高阶的语法，详见lavaan的help文档，一般的结构方程建模分析用不到，就不再列出。
二、模型的三种表示方法
以验证性因子分析举例说明，对于如下图所示的模型：
方法一：最简化描述
只需指定最基本的要素即可，其他的由函数自动实现，对模型的控制力度最弱。只使用于函数cfa()和sem()
model<-'visual=~x1+x2+x3 textual=~x4+x5+x6 speed=~x7+x8+x9' fit <- cfa(model, data = HolzingerSwineford1939)
需要注意的是，这种指定模型的方式在进行拟合时，会默认指定潜变量的第一个测量变量的因子载荷为1，如果要指定潜变量的方差为1，可以：
modelbis <- 'visual =~ NAx1 + x2 + x3 textual =~ NAx4 + x5 + x6 speed =~ NAx7 + x8 + x9 visual ~~ 1visual textual ~~ 1textual speed ~~ 1speed'
方法二：完全描述
需要指定所有的要素，对模型控制力最强，适用于lavaan()函数，适合高阶使用者
modelfull<- ' visual =~ 1x1 + x2 +x3 textual =~ 1x4 + x5 + x6 speed =~ 1x7 + x8 +x9 x1 ~~ x1 x2 ~~ x2 x3 ~~ x3 x4 ~~ x4 x5 ~~ x5 x6 ~~ x6 x7 ~~ x7 x8 ~~ x8 x9 ~~ x9 visual ~~ visual textual ~~ textual speed ~~ speed visual ~~ textual +speed textual ~~ speed' fit <- lavaan(modelfull, data = HolzingerSwineford1939)
方法三：不完全描述
最简化和完全描述的混合版，在拟合时增加 auto 参数，适用于lavaan()函数
modelmixed<- '# latent variables visual =~ 1x1 + x2 +x3 textual =~ 1x4 + x5 + x6 speed =~ 1x7 + x8 +x9 # factor covariances visual ~~ textual + speed textual ~~ speed' fit <- lavaan(modelmixed, data = HolzingerSwineford1939, autovar = TRUE)
可以设定的参数详见help帮助文档
PS:可以在lavaan()函数里设置参数mimic="Mplus"获得与Mplus在数值和外观上相似的结果，设置mimic="EQS",输出与EQS在数值上相似的结果
三、拟合结果的查看
查看拟合结果的最简单方法是用summary()函数，例如
summary(fit, fitmeasures=TRUE)
但summary()只适合展示结果，parameterEstimates()会返回一个数据框，方便进一步的处理
parameterEstimates(fit,ci=FALSE,standardized = TRUE)
获得大于10的修正指数
MI<- modificationindices(fit) subset(MI,mi>10)
此外，还有其他的展示拟合结果的函数，功能还是蛮强大的
四、结构方程模型
（1）设定模型
model<- ' # measurement model ind60 =~ x1 + x2 +x3 dem60 =~ y1 + y2 + y3 + y4 dem65 =~ y5 + y6 + y7 + y8 # regressions dem60 ~ ind60 dem65 ~ ind60 + dem60 # redisual covariances y1 ~~ y5 y2 ~~ y4 +y6 y3 ~~ y7 y4 ~~ y8 y6 ~~ y8'
（2）模型拟合
fit <- sem(model, data = PoliticalDemocracy) summary(fit, standardized = TRUE)
（3）给回归系数设置标签
给回归系数设定标签在做有约束条件的结构方程模型时会很有用。当两个参数具有相同的标签时，会被视为同一个，只计算一次。
modelequal <- '# measurement model ind60 =~ x1 + x2 + x3 + dem60 =~ y1 + d1y2 + d2y3 + d3y4 dem65 =~ y5 + d1y6 + d2y7 + d3y8 # regressions dem60 ~ ind60 dem65 ~ ind60 + dem60 # residual covariances y1 ~~ y5 y2 ~~ y4 + y6 y3 ~~ y7 y4 ~~ y8 y6 ~~ y8'
（4）多组比较
anova(fit, fitequal)
anova()会计算出卡方差异检验
（5）拟合系数
lavaan包可以高度定制化的计算出你想要的拟合指标值，例如，我想计算出卡方、自由度、p值、CFI、NFI、IFI、RMSEA、EVCI的值
fitMeasures(fit,c("chisq","df","pvalue","cfi","nfi","ifi","rmsea","EVCI"))
（6）多组结构方程
在拟合函数里面设置 group参数即可实现，同样的可以设置groupequal参数引入等式限制
五、作图
Amos以作图化 *** 作见长，目前版本的Mplus也可以实现作图，那R语言呢，自然也是可以的，只不过是另一个包——semPlot，其中的semPaths()函数。
简单介绍一下semPaths()中的主要函数
semPaths(object, what = "paths", whatLabels, layout = "tree", ……）
（1）object：是拟合的对象，就是上文中的“fit”
（2）what：设定图中线的属性， 默认为paths,图中所有的线都为灰色，不显示参数估计值；
semPaths(fit)
若what设定为est、par，则展示估计值，并将线的颜色、粗细、透明度根据参数估计值的大小和显著性做出改变
semPaths(fit,what = "est")
若设置为stand、std，则展示标准参数估计
semPaths(fit,what = "stand")
若设置为eq、cons，则与默认path相同，如果有限制等式，被限制的相同参数会打上相同的颜色；
（3）whatLabels：设定图中线的标签
name、label、path、diagram:将边名作为展示的标签
est、par:参数估计值作为边的标签
stand、std:标准参数估计值作为边的标签
eq、cons：参数号作为标签，0表示固定参数，被限制相同的参数编号相同
no、omit、hide、invisible：隐藏标签
（4）layout:布局
主要有树状和环状两种布局，每种布局又分别有两种风格。
默认为“tree”,树状的第二种风格如下图，比第一种看起来舒服都了
semPaths(fit,layout = "tree2")
第一种环状
semPaths(fit,layout = "circle")
额，都揉成一团了！
试试第二种风格
semPaths(fit,layout = "circle2")
还好一点。如果把Rstudio默认的尺寸设计好，作图效果会更棒。
还有一种叫spring的布局，春OR泉？
semPaths(fit,layout = "spring")
看起来跟环状的很像。
详细内容可以阅读以下文献，以及相应的help文档：
[1]Rosseel Y lavaan: An R package for structural equation modeling[J] Journal of Statistical Software, 2012, 48(2): 1-36

今天被粉丝发的文章给难住了，又偷偷去学习了一下竞争风险模型，想起之前写的关于竞争风险模型的做法，真的都是皮毛哟，大家见笑了。想着就顺便把所有的生存分析的知识和R语言的做法和论文报告方法都给大家梳理一遍。
什么时候用生存分析
当你关心结局和结局发生时间的时候，就要考虑生存分析了，这种既有结局又有时间的数据叫做生存数据，英文叫做Time-to-event data 只不过因为这个方法医学上用来分析存活情况用的多，所以得名生存分析，反正你就记住一个例子，我要研究汽车发生故障，我也应该用生存分析，因为我既关心是不是有故障，我还关心用了多久（跑了多远）才出故障，就是既有time，又有event，Time-to-event data就用生存分析。
基本概念
首先是删失，对象失访了，脱落了，出现结局之前随访结束了，都叫做删失：
删失又分为左删失，区间删失和右删失，图示如下：
比如我想研究得了A病的人的生存情况，存在的所有可能情形为：
第一种，研究的开始的时候有人已经有A病，这个时候人家已经活了一段时间了，具体多久我不知道，叫做左删失；
第二种，入组随访的时候没病，中途得了A病死了，什么时候得的，没记录下来，叫区间删失；
第三种，得了A病，一直活到了研究结束还没死，叫做右删失。
你看，所有的删失情况造成的后果都是我们没法准确估计发生结局的时间，这也是其名字删失的由来，对于这类数据就需记录为删失数据。
生存分析的种类有哪些
具体的种类是为了回答具体的问题，我们做生存分析常常要回答的问题如下：
一个是描述生存情况，一个是比较，再一个就是探究影响因素。
比如我随访了很多病人，我就想知道随着时间变化这群人的生存概率是如何变化的（描述）？我就可以用简单粗暴的Kaplan-Meier method，又叫乘积极限法，基本思想就是此刻的生存概率等于上一时刻的生存概率乘以此刻的存活率。
比如我手上有如下数据：
time是随访时间，status是结局，我就可以写出如下代码拟合出总体人群的生存曲线：
fit1 <- survfit(Surv(time, status) ~ 1, data = mydata)summary（fit1）
并且得到每个时间点，病人生存概率的估计值和标准误：
如果我还恰好收集了病人的性别，我又想看看男病人和女病人生存情况是不是有不同（比较），我可以对生存曲线分组比较，代码如下：
fit2<- survfit(Surv(time, status) ~ sex, data = mydata)summary（fit2）
输出两组生存曲线比较的log-rank test的统计量：
survdiff(Surv(time, status) ~ sex, data = mydata)
并且我们还可以进行复杂分组的组间比较，大家可以翻看我之前的文章。
但是大家明白，KM法始终是在做单因素分析，而且都是做的分类变量的单因素分析，我们经常还会有的需求是考虑各种混杂的情况下去探讨影响生存时间的因素，这个时候我们就要用到The Cox regression model了。
模型形式如下：
上面的式子把h0移项到左边，等号两边同时取对数就成了一个线性模型，和广义线性模型的理解一样一样的，b就是我们的影响因素的系数，解释为lnHR的改变量，其为正就是危险因素，为负就是保护因素：
The quantities exp(bi) are called hazard ratios (HR) A value of bi greater than zero, or equivalently a hazard ratio greater than one, indicates that as the value of the ith covariate increases, the event hazard increases and thus the length of survival decreases
再观测一下上面的式子，我们拟合了预测因素对风险比例（t时刻风险比上基础风险）的模型，这个时候暗含的假设就是就是对于每个人在任何时刻风险只有一个常数，就是在所有预测因素的作用下，各个时间的风险是不变的，这个就叫做Proportional Hazards Assumption, 比例风险假设。
做COX比例风险模型的时候都有必要验证这个假设是不是满足，具体方法如下：
我们需要先做一个cox模型，拟合cox模型需要用到coxph函数，代码如下：
rescox <- coxph(Surv(time, status) ~ age + sex + wtloss, data = lung)rescox
模型输出结果里面会有预测变量的β值（coef）和其标准误，有HR（exp(coef)）值，还有预测变量的显著性检验结果：
我们将这个模型对象直接喂给coxzph，就可以得到风险比例假设的检验结果了：
testph <- coxzph(rescox)testph
可以看到，我们的p值都大于005，即都满足ph假设。
我们还可以从图上看，看scaled Schoenfeld residuals是不是和时间独立，如果独立则满足ph假设：
ggcoxzph(testph)
有竞争事件时
上面写了只有一个截距事件时生存分析的不同目的，描述，比较，和影响因素探讨，接着再来看存在竞争事件的时候该如何描述，比较，和探讨影响因素。
生存分析的另外情况就是竞争风险模型，就是time to event的event有多种，一种发生了另外一种就不可能发生了，这个就是竞争，比如好多文献中都会列举的经典例子：就是在造血干细胞移植人群中，我们关心其疾病复发风险，但是有些人因为移植死了（TRM），死了之后肯定也谈不上复发，如果你把因为移植死了的人都作为删失数据，肯定也是不对的，这个里面就会有两个竞争结局相关性造成的问题，此时我们应该用竞争风险模型来分析。
For example, disease relapse is an event of interest in studies of allogeneic hematopoietic stem cell transplantation (HSCT), as is mortality related to complications of transplantation (transplant-related mortality or TRM) Relapse and TRM are not independent in this setting because these two events are likely related to immunologic effector mechanisms following HSCT, whereby efforts to reduce TRM may adversely affect the risk of relapse; moreover, patients who die from TRM cannot be at further risk of relapse
在比例风险中我们的结局事件只有一个，我们关心的是哪些因素对事件发生的风险有影响。在竞争风险模型中我们关注的地方又变了，因为我们有好几个结局事件，这个时候我们会关心在竞争事件存在的情况下各个结局事件的累计发病函数是如何随着时间变化的，以及如何来比较不同的累计发病函数，以及如何探讨影响因素（competing risks regression analysis）。
之前写的在非竞争风险模型中累计发病率的组间比较可以用KM法，将纵坐标换一换，用log-rank检验，在竞争风险模型中我们需要用Fine and Gray提到的方法来做（Gray’s test），如果比如说我发现两组（治疗组和对照组）的累计发病风险不一样，我肯定还想探讨哪些因素影响累计发病风险，之前是用COX比例风险模型做的，在竞争结局存在的情况下我们依然是得用Fine and Gray提出的模型（Fine and Gray Model），叫做竞争风险回归模型：
Fine and Gray (6) proposed a method for direct regression modeling of the effect of covariates on the cumulative incidence function for competing risks data As in any other regression analysis, modeling cumulative incidence functions for competing risks can be used to identify potential prognostic factors for a particular failure in the presence of competing risks, or to assess a prognostic factor of interest after adjusting for other potential risk factors in the model
首先看竞争风险时候累计发病曲线的比较方法。我手上有数据如下：
其中dis是疾病类型，2分类，ftime是时间，status是结局事件，结局事件有3个水平，多的1个水平是竞争事件。现在我关心不同疾病人群各个事件累积发病曲线有无不同，我可以用cuminc函数结合plot画出各个组的累计发病曲线：
fit=cuminc (ftime, status, dis, cencode = 0) plot(fit)
fit对象的结果中还有每条曲线的比较，从比较结果可以看出两组在事件2的累积发病曲线上是有显著差异的：
上面介绍的方法相当于没有竞争风险的时候的KM法，通过上面的方法我们知道有了不同风险结局累计发病曲线的差异，继续我们会继续看影响因素，要做的就是竞争风险回归模型了，需要用到的函数就是crr。
比如我手上有如下数据，除了时间，结局还包括每个病人的像sex，age等等协变量，我想探讨说这些因素是如何影响病人某个结局的，我就可以写出一个竞争风险回归模型：
mod1 <- crr(ftime,Status,x)summary(mod1)
上面的代码中x是自变量矩阵，运行代码输出竞争风险回归模型的结果如下：
到这儿基本就写完了所有的生存分析的情况，接着再结合一篇论文看看报告方法，论文就是下面这篇，是我随意查的：
S Chen, H Sun, M Heng, X Tong, P Geldsetzer, Z Wang, P Wu, J Yang, Y Hu, C
Wang, T Bärnighausen, Factors Predicting Progression to Severe COVID-19: A Competing Risk Survival Analysis of 1753 Patients in Community Isolation in Wuhan, China, Engineering (2021), doi: >