# 一、R基本操作# 1、将数据文件mydata1.txt按照以下要求整理成标准形式。#(1)读入数据文件mydata.txt命名为insurance。insurance<-read.table("mydata1.txt")head(insurance)dim(insurance)#192个数据#(2)将insurance转换为3列的矩阵。insurance<-matrix(insurance$V1,nrow = 64,ncol = 3)#nrow =192/3=64insurance#(3)将insurance转换为数据框。insurance<-as.data.frame(insurance)class(insurance)#(4)将列名命名为"District", "Holders"和"Claims"。names(insurance)<-c("District", "Holders","Claims")insurance#(5)随机无放回抽取50行数据。sub<-insurance[sample(1:nrow(insurance),50),]#无放回不用设置replacesub#(6)将抽样数据写入result1.txt。write.table(sub,"result1.txt",row.names = FALSE)####################################################################### 2、将数据文件mydata2.txt按照以下要求整理成标准形式。#(1)读入数据文件mydata2.txt命名为iris。iris<-read.table("mydata2.txt")head(iris)dim(iris)#600个数据#(2)将iris转换为4列的矩阵。iris<-matrix(iris$V1,nrow = 150,ncol = 4)#nrow =600/3=150iris#(3)将iris转换为数据框。iris<-as.data.frame(iris)class(iris)#(4)将列名命名为"Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width"。names(iris)<-c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")iris#(5)随机无放回抽取90行数据。sub<-iris[sample(1:nrow(iris),90),]#无放回不用设置replacesub#(6)将抽样数据写入result2.txt。write.table(sub,"result2.txt",row.names = FALSE)####################################################################### 3.将数据文件data.csv按照以下要求进行数据预处理。#(1)读入数据文件data.csv命名为nhanes2。nhanes2<-read.csv("data.csv")#(2) 载入缺失值处理所需要的包。install.packages("lattice")install.packages("MASS")install.packages("nnet")library(lattice)library(MASS)library(nnet)#(3) 判断nhanes2是否存在缺失值。sum(is.na(nhanes2))#(4) 利用插补法处理chl变量的缺失值。sub=which(is.na(nhanes2[,4])==TRUE)#在数据集中chl变量是第4列，返回nhanes2数据集中第4列为NA的行dataTR<-nhanes2[-sub,]#将第4列不为NA的数存入数据集dataTRdataTE<-nhanes2[sub,]#将第4列为NA的数存入数据集dataTE中dataTE[,4]<-sample(dataTR[,4],length(dataTE[,4]),replace = T)#在非缺失值中简单抽样dataTE #(5) 将插补法处理后的数据写入result3.txt。write.table(dataTE,"result3.txt",row.names = FALSE)#############################################################################################################################################二、函数调用#1、测得某班学术X（身高(cm)）与Y（体重(kg)）的数据如下,试画出散点图,建立线性回归方程，并作进一步分析。# (1) 建立数据集，并画出散点图，考察数据点的分布趋势，看是否呈直线条状分布。x1<-c(171,175,159,155,152,158,154,164,168,166,159,164)#身高y1<-c(57,64,41,38,35,44,41,51,57,49,47,46)#体重#构建数据集model <- data.frame(x1,y1) #探索性分析-做散点图查看数据的分布情况：plot(x1,y1)# (2)进行回归分析,列出回归方程，画拟合线，并对结果进行解读。# 用lm()函数构建模型lm.reg<-lm(y1~ x1)# 添加回归曲线查看拟合效果 abline(lm.reg) #模型解读summary(lm.reg)# (3)对回归系数进行假设检验。anova(lm.reg) # 回归模型的方差分析summary(lm.reg) #回归系数t检验：提取模型计算结果，其中有t检验的结果# (4)对回归模型进行诊断。#模型检验对方程进行进一步检验，以检查回归方程是否满足模型的先验条件及模型的稳健性。par(mfrow=c(2,2))#画布分面plot(lm.reg)#结果解读：#1.左上图：残差与拟合图，理论上散点应该散乱的分布在横线两侧；#2.右上图：正太Q-Q图，用于检验因变量的正太分布性，若服从正太分布，则散点应分布在一条直线线#3.左下图：齐方差检验，若满足其方差，则散点在水平线周围随机分布#4.右下图：独立性检验，即一个样本是否会影响另一个样本##################################################################2、研究某抗心律失常药对电刺激狗右心室致颤阙的影响,实验测得狗静脉注射不同剂量的抗心律失常药与右心室致颤阙的数据如下,试画出散点图,建立线性回归方程，并作进一步分析。# (1) 建立数据集，并画出散点图，考察数据点的分布趋势，看是否呈直线条状分布。x <- c(1,3,5,7,9)y <- c(8.03, 14.97, 19.23, 27.83, 36.23)#构建数据集model <- data.frame(x,y) #探索性分析-做散点图查看数据的分布情况：plot(model)#画散点图# (2)进行回归分析,列出回归方程，画拟合线，并对结果进行解读。# 用lm()函数构建模型fm <- lm(y ~ x)#建立回归模型fm# 添加回归曲线查看拟合效果abline(fm)# 添加回归曲线至散点图 #模型解读summary(fm)# (3)对回归系数进行假设检验。anova(fm) # 回归模型的方差分析summary(fm) # 提取模型计算结果，其中有t检验的结果# (4)对回归模型进行诊断。#模型检验对方程进行进一步检验，以检查回归方程是否满足模型的先验条件及模型的稳健性。par(mfrow=c(2,2))#画布分面plot(fm)#结果解读：#1.左上图：残差与拟合图，理论上散点应该散乱的分布在横线两侧；#2.右上图：正太Q-Q图，用于检验因变量的正太分布性，若服从正太分布，则散点应分布在一条直线线#3.左下图：齐方差检验，若满足其方差，则散点在水平线周围随机分布#4.右下图：独立性检验，即一个样本是否会影响另一个样本################################################################### 3、countries数据集含有69个国家和地区的出生率与死亡率。# (1) 请使用K-均值聚类将样本点聚为3个类别。countries=read.csv("countries.csv")head(countries)#查看前6行names(countries)=c("country","birth","death")#修改变量名称var=as.character(countries$country)#将变量country转为字符型并赋值给varfor(i in 1:69) row.names(countries)[i]=var[i]#将数据集的行名命名为国家名称km1=kmeans(countries[,-1],center=3)#用kmeans算法对countries数据集进行聚类# (2) 输出聚类结果及各类别的中心点坐标。km1$cluster#获取类别km1$centers#获取中心点坐标# (3) 绘制聚类结果将中心点以星号标识。#画出聚为四类的类别图，标注中心点。plot(countries[,-1],pch=c(1,2,3))#将中心点用星号标示出来points(km1$centers,pch=8,col="red")#对中心点添加标注legend(km1$centers[1,1],km1$centers[1,2],"Center_1",bty="n",xjust=0.5,cex=0.8)legend(km1$centers[2,1],km1$centers[2,2],"Center_2",bty="n",xjust=0.5,cex=0.8)legend(km1$centers[3,1],km1$centers[3,2],"Center_3",bty="n",xjust=0.5,cex=0.8)# (4) 判断与中国大陆同属于一个类别的国家和地区有哪些。cluster_CHINA=km1$cluster[which(countries$country=="CHINA")]which(km1$cluster==cluster_CHINA)###############################################################################################################################三、数据分析# 1、使用arules软件包中的Groceries数据集，该数据集是某一食品杂货店一个月的真实交易数据，使用R完成以下要求：(软件包：arules；数据集：Groceries； 函数：apriori())# (1)利用apriori()函数进行关联分析，支持度为0.01，置信度为0.5。install.packages("arules")library(arules)data("Groceries")rules0<-apriori(Groceries,parameter=list(support=0.01,confidence=0.5))inspect(rules0[1:10])# (2)利用sort()函数按照支持度排序。rules.sorted_sup<-sort(rules0,by="support")inspect(rules.sorted_sup[1:5])# (3)捆绑销售：寻找蛋黄酱（mayonnaise）的捆绑商品。（supp=0.001,conf=0.1,minlen=2, maxlen=6）rules1=apriori(Groceries,parameter=list(minlen=2,maxlen=6,supp=0.001,conf=0.1),appearance=list(rhs="mayonnaise",default="lhs"))inspect(rules1)# (4)查看销量最高的商品。itemsets_apr=apriori(Groceries,parameter=list(supp=0.001,target="frequent itemsets"),control=list(sort=-1))inspect(itemsets_apr[1:5])# (5)适合捆绑销售的商品。（supp=0.001,minlen=2, maxlen=3）itemsets_apr1=eclat(Groceries,parameter=list(supp=0.001,minlen=2,maxlen=3,target="frequent itemsets"),control=list(sort=-1))inspect(itemsets_apr1[1:5])# (6)关联规则的可视化（support=0.001,con=0.5）install.packages("arulesViz")library(arulesViz)rules5=apriori(Groceries,parameter=list(support=0.002,con=0.5))rules5plot(rules5)######################################################################## 2、根据breast-cancer-wisconsin.csv威斯康星州乳腺癌数据集，通过对数据的分析，提取出关键特征来判断乳腺癌患病情况。(软件包：rpart；函数：rpart()。)# (1)属性名依次设置为"编号","肿块厚度","肿块大小","肿块形状","边缘黏附","单个表皮细胞大小","细胞核大小","染色质","细胞核常规","有丝分裂","类别"),并将类别为2的设为"良性",为4的设为"恶性"。install.packages("rpart")library(rpart)install.packages("rpart.plot")library(rpart.plot)#############加载数据breast.cancer<-read.csv('breast-cancer-wisconsin.csv',header=F)head(breast.cancer)#数据整理names(breast.cancer)=c("编号","肿块厚度","肿块大小","肿块形状","边缘黏附","单个表皮细胞大小","细胞核大小","染色质","细胞核常规","有丝分裂","类别")breast.cancer$类别[breast.cancer$类别==2]="良性"breast.cancer$类别[breast.cancer$类别==4]="恶性"head(breast.cancer)# (2)抽取训练数据集为原数据的70%，测试数据集取30%。#数据预处理(分层抽样,划分训练集和测试集)#分别计算良性和恶性组中应抽取测试集样本数，记为a，ba=round(0.3*sum(breast.cancer$类别=="良性"))b=round(0.3*sum(breast.cancer$类别=="恶性"))ab #输出a,b值install.packages("sampling")library(sampling)#使用strata函数对数据集中的“分组油耗”变量进行分层抽样sub=strata(breast.cancer,stratanames="类别",size=c(b,a),method="srswor")sub #所抽出的所有测试集样本信息#生成训练集train1和测试集test1train1=breast.cancer[-sub$ID_unit,]test1=breast.cancer[sub$ID_unit,]nrow(train1)nrow(test1) #显示训练集和测试集的行数，检查两者比例是否为7：3# (3) minsplit=5，建立决策树。#CART建立分类树formula_cla=类别~肿块厚度+肿块大小+肿块形状+边缘黏附+单个表皮细胞大小+细胞核大小+染色质+细胞核常规+有丝分裂cla1=rpart(formula_cla,train1,method="class",minsplit=5)#cla1# (4)选择cp=0.05来剪枝。######修改cp的值cla2=rpart(formula_cla,train1,method="class",minsplit=5,cp=0.05)cla2# (5)画出type为2和4的树图。rpart.plot(cla1,type=2)#修改typerpart.plot(cla1,type=4) # (6)测试数据进行预测，并输出混淆矩阵，给出模型准确率为。#预测pre1=predict(cla1,test1,type="class")pre1table(test1$类别,pre1)#获取混淆矩阵#计算样本错误率error1<-sum(as.numeric(pre1!=test1$类别))/nrow(test1)error1#################################################################### 3、美国科罗拉多州某加油站连续 57 天的OVERSHORTS序列“OVERSHORTS.csv”# (1) 判断该序列的平稳性与纯随机性。# （时序图检验、白噪声检验）install.packages("fUnitRoots")install.packages("TSA")install.packages("forecast")install.packages("zoo")library(fUnitRoots)library(TSA)library(forecast)library(zoo)#读取数据c<-read.csv("OVERSHORTS.csv")#转换为时间序列overshort<-ts(c$overshort,start = 1)#平稳性，纯随机（白噪声检验）## 绘制序列的时间序列图plot.ts(overshort, xlab = "time", ylab = "prop")##对序列做单位根检验unitrootTest(overshort)##对序列做白噪声检验Box.test(overshort, lag = 1, type = "Ljung-Box")# (2) 如果序列平稳且非白噪声，选择适当模型拟合该序列的发展。(10分)# （模型的识别、参数估计（模型显著性、模型参数的显著性））#模型识别##观察自相关，偏自相关图，模型定阶par(mfrow=c(1,2))acf(overshort)###衰减到零是突然的，所以自相关系数1阶截尾pacf(overshort)### 衰减到零不是突然的，所以偏相关系数托尾# 推荐模型为 MA(1)##或者对序列进行模型识别，自动定阶auto.arima(overshort)# 推荐模型为 MA(1)#参数估计###模型检验x.fit<-arima(overshort,order=c(0,0,1),method="ML")x.fit##对残差x.fit$residual进行白噪声检验for(i in 1:2) print(Box.test(x.fit$residual,lag=6*i))##P>0.05,接受原假设，即残差为白噪声，所以拟合模型显著有效####参数检验###模型参数的显著性检验t1<--0.8477/0.1206pt(t1,df=56,lower.tail=T) ###p<0.05参数显著非零t0<--4.7942/1.0253pt(t0,df=56,lower.tail=T) ###p<0.05参数显著非零# (3) 利用拟合模型，预测该加油站未来5天的OVERSHORTS。(10分)# （模型预测、绘制预测图）####模型预测c<-read.csv("OVERSHORTS.csv")x<-ts(c$overshort,start=1)x.fit<-arima(x,order=c(0,0,1))x.fitx.fore<-forecast(x.fit,h=5)#预测x.foreplot(x.fore)###############################################################4、使用是survival软件包中的“pbc”数据集，该数据集记录的是肝硬化数据， 使用R完成一下要求：(软件包：survival；数据集：pbc； 函数：Surv()、survfit()、survdiff()、coxph()、cox.zph(), 将答案保存在“姓名.doc”文件中。)# (1)生成生存分析对象，拟合生存曲线模型。install.packages("survival") #安装survival包library(survival) #加载survival包#使用survival包自带的“pbc”数据集为例（418*20） data("pbc")str(pbc)head(pbc)#生成生存分析对象Sur_Obj<-Surv(pbc$time,pbc$status)Sur_Obj#拟合曲线模型model<-survfit(Sur_Obj~1) summary(model)# (2)两种方法绘制生存曲线。plot(model,ylab = "生存率",xlab="天")#用survminer进行漂亮的展示install.packages("survminer")library(survminer) ggsurvplot(model, data = pbc)# (3)进行单因素比较分析，并进行结果解释。#survdiff(formula)函数进行log-rank检验。survdiff(Sur_Obj~pbc$trt) #trt是分组条件# (4)考虑年龄，性别以及trt是否会影响肝硬化的生存时间,进行多因素分析Cox模型的建立，并进行结果解释。coxmodel<-coxph(Sur_Obj~pbc$age+pbc$sex+pbc$bili)coxmodel# (5)模型诊断——PH检验。zphmodel<-cox.zph(coxmodel)zphmodel############################################################### 5、life.csv为50位急性淋巴细胞白血病病人的数据，包括：入院治疗时取得外辕血中细胞数X1，淋巴结浸润等级X2，出院后有无巩固治疗X3(1表示有巩固治疗，0表示无巩固治疗)；随访后，变量Y=0表示生存期在1年以内，Y=1表示生存时间在1年以上，使用R完成一下要求：(函数：glm()，predict()。)# (1)建立全变量logistic回归，对模型结果进行解释。life<-read.csv("life.csv")#建立全变量logistic回归glm.sol<-glm(Y~X1+X2+X3, family=binomial, data=life)#回归模型解读summary(glm.sol)# (2)预测当X1=5,X2=2,X3=0时，y的概率是多少？pre<-predict(glm.sol, data.frame(X1=5,X2=2,X3=0))p<-exp(pre)/(1+exp(pre))p# (3)预测当X1=5,X2=2,X3=1时，y的概率是多少？(6分)pre<-predict(glm.sol, data.frame(X1=5,X2=2,X3=1))p<-exp(pre)/(1+exp(pre))p# (4)对回归模型参数进行检验，用step()函数做变量筛选。step(glm.sol)glm.new<-glm(Y~X2+X3, family=binomial, data=life)summary(glm.new)# (5)对筛选后的变量进行建模，预测。pre<-predict(glm.new, data.frame(X2=2,X3=0))p<-exp(pre)/(1+exp(pre))ppre<-predict(glm.new, data.frame(X2=2,X3=1))p<-exp(pre)/(1+exp(pre))p

主成分分析

主成分分析（(Principal Component Analysis，PCA）是一种数据降维技巧，它能将大量相关变量转化为一组很少的不相关变量，这些无关变量称为主成分（原来变量的线性组合）。整体思想就是化繁为简，抓住问题关键，也就是降维思想。

主成分分析法是通过恰当的数学变换，使新变量——主成分成为原变量的线性组合，并选取少数几个在变差总信息量中比例较大的主成分来分析事物的一种方法。主成分在变差信息量中的比例越大，它在综合评价中的作用就越大。

因子分析

探索性因子分析法（Exploratory Factor Analysis，EFA）是一系列用来发现一组变量的潜在结构的方法。它通过寻找一组更小的、潜在的或隐藏的结构来解释已观测到的、显式的变量间的关系。

PCA与EFA模型间的区别

参见图14-1。主成分（PC1和PC2）是观测变量（X1到X5）的线性组合。形成线性组合的权重都是通过最大化各主成分所解释的方差来获得，同时还要保证个主成分间不相关。相反，因子（F1和F2）被当做是观测变量的结构基础或“原因”，而不是它们的线性组合。

R的基础安装包提供了PCA和EFA的函数，分别为princomp()和factanal()。

最常见的分析步骤

(1)数据预处理。PCA和EFA都根据观测变量间的相关性来推导结果。用户可以输入原始数据矩阵或者相关系数矩阵到principal()和fa()函数中。若输入初始数据，相关系数矩阵将会被自动计算，在计算前请确保数据中没有缺失值。

(2)选择因子模型。判断是PCA（数据降维）还是EFA（发现潜在结构）更符合你的研究目标。如果选择EFA方法，你还需要选择一种估计因子模型的方法（如最大似然估计）。

(3)判断要选择的主成分/因子数目。

(4)选择主成分/因子。

(5)旋转主成分/因子。

(6)解释结果。

(7)计算主成分或因子得分。

PCA的目标是用一组较少的不相关变量代替大量相关变量，同时尽可能保留初始变量的信息，这些推导所得的变量称为主成分，它们是观测变量的线性组合。如第一主成分为：

它是k个观测变量的加权组合，对初始变量集的方差解释性最大。第二主成分也是初始变量的线性组合，对方差的解释性排第二，同时与第一主成分正交（不相关）。后面每一个主成分都最大化它对方差的解释程度，同时与之前所有的主成分都正交。理论上来说，你可以选取与变量数相同的主成分，但从实用的角度来看，我们都希望能用较少的主成分来近似全变量集。

主成分与原始变量之间的关系

（1）主成分保留了原始变量绝大多数信息。

（2）主成分的个数大大少于原始变量的数目。

（3）各个主成分之间互不相关。

（4）每个主成分都是原始变量的线性组合。

数据集USJudgeRatings包含了律师对美国高等法院法官的评分。数据框包含43个观测，12个变量。

用来判断PCA中需要多少个主成分的准则：

根据先验经验和理论知识判断主成分数；

根据要解释变量方差的积累值的阈值来判断需要的主成分数；

通过检查变量间k × k的相关系数矩阵来判断保留的主成分数。

最常见的是基于特征值的方法。每个主成分都与相关系数矩阵的特征值相关联，第一主成分与最大的特征值相关联，第二主成分与第二大的特征值相关联，依此类推。

Kaiser-Harris准则建议保留特征值大于1的主成分，特征值小于1的成分所解释的方差比包含在单个变量中的方差更少。Cattell碎石检验则绘制了特征值与主成分数的图形。这类图形可以清晰地展示图形弯曲状况，在图形变化最大处之上的主成分都可保留。最后，你还可以进行模拟，依据与初始矩阵相同大小的随机数据矩阵来判断要提取的特征值。若基于真实数据的某个特征值大于一组随机数据矩阵相应的平均特征值，那么该主成分可以保留。该方法称作平行分析。

图形解读：线段和x符号组成的图（蓝色线）：特征值曲线；

红色虚线：根据100个随机数据矩阵推导出来的平均特征值曲线；

绿色实线：特征值准则线（即：y=1的水平线）

判别标准：特征值大于平均特征值，且大于y=1的特征值准则线，被认为是可保留的主成分。根据判别标准，保留1个主成分即可。

fa.parallel函数学习

fa.parallel(data,n.obs=,fa=”pc”/”both”,n.iter=100,show.legend=T/F)

data：原始数据数据框；

n.obs：当data是相关系数矩阵时，给出原始数据（非原始变量）个数，data是原始数据矩阵时忽略此参数；

fa：“pc”为仅计算主成分，“fa”为因子分析，“both”为计算主成分及因子；

n.iter：模拟平行分析次数；

show.legend：显示图例。

principal(r, nfactors = , rotate = , scores = )

r：相关系数矩阵或原始数据矩阵；

nfactors：设定主成分数（默认为1）；

rotate：指定旋转的方法，默认最大方差旋转（varimax）。

scores：设定是否需要计算主成分得分（默认不需要）。

PC1栏包含了成分载荷，指观测变量与主成分的相关系数。如果提取不止一个主成分，那么还将会有PC2、PC3等栏。成分载荷（component loadings）可用来解释主成分的含义，解释主成分与各变量的相关程度。

h2栏为成分公因子方差，即主成分对每个变量的方差解释度。

u2栏为成分唯一性，即方差无法被主成分解释的部分（1-h2）。

SS loadings包含了与主成分相关联的特征值，其含义是与特定主成分相关联的标准化后的方差值，即可以通过它来看90%的方差可以被多少个成分解释，从而选出主成分（即可使用nfactors=原始变量个数来把所有特征值查出，当然也可以直接通过eigen函数对它的相关矩阵进行查特征值）。

Proportion Var表示每个主成分对整个数据集的解释程度。

Cumulative Var表示各主成分解释程度之和。

Proportion Explained及Cumulative Proportion分别为按现有总解释方差百分比划分主成分及其累积百分比。

结果解读：第一主成分（PC1）与每个变量都高度相关，也就是说，它是一个可用来进行一般性评价的维度。ORAL变量99.1%的方差都可以被PC1来解释，仅仅有0.91%的方差不能被PC1解释。第一主成分解释了11个变量92%的方差。

结果解读：通过碎石图可以判定选择的主成分个数为2个。

结果解读：从结果Proportion Var： 0.58和0.22可以判定，第一主成分解释了身体测量指标58%的方差，而第二主成分解释了22%，两者总共解释了81%的方差。对于高度变量，两者则共解释了其88%的方差。

旋转是一系列将成分载荷阵变得更容易解释的数学方法，它们尽可能地对成分去噪。旋转方法有两种：使选择的成分保持不相关（正交旋转），和让它们变得相关（斜交旋转）。旋转方法也会依据去噪定义的不同而不同。最流行的正交旋转是方差极大旋转，它试图对载荷阵的列进行去噪，使得每个成分只是由一组有限的变量来解释（即载荷阵每列只有少数几个很大的载荷，其他都是很小的载荷）。 结果列表中列的名字都从PC变成了RC，以表示成分被旋转。

当scores = TRUE时，主成分得分存储在principal()函数返回对象的scores元素中。

如果你的目标是寻求可解释观测变量的潜在隐含变量，可使用因子分析。

EFA的目标是通过发掘隐藏在数据下的一组较少的、更为基本的无法观测的变量，来解释一

组可观测变量的相关性。这些虚拟的、无法观测的变量称作因子。（每个因子被认为可解释多个

观测变量间共有的方差，因此准确来说，它们应该称作公共因子。）

其中 是第i个可观测变量（i = 1…k）， 是公共因子（j = 1…p），并且p<k。 是 变量独有的部分（无法被公共因子解释）。 可认为是每个因子对复合而成的可观测变量的贡献值。

碎石检验的前两个特征值（三角形）都在拐角处之上，并且大于基于100次模拟数据矩阵的特征值均值。对于EFA，Kaiser-Harris准则的特征值数大于0，而不是1。

结果解读：PCA结果建议提取一个或者两个成分，EFA建议提取两个因子。

fa(r, nfactors=, n.obs=, rotate=, scores=, fm=)

 r是相关系数矩阵或者原始数据矩阵；

 nfactors设定提取的因子数（默认为1）；

 n.obs是观测数（输入相关系数矩阵时需要填写）；

 rotate设定旋转的方法（默认互变异数最小法）；

 scores设定是否计算因子得分（默认不计算）；

 fm设定因子化方法（默认极小残差法）。

与PCA不同，提取公共因子的方法很多，包括最大似然法（ml）、主轴迭代法（pa）、加权最小二乘法（wls）、广义加权最小二乘法（gls）和最小残差法（minres）。统计学家青睐使用最大似然法，因为它有良好的统计性质。

结果解读：两个因子的Proportion Var分别为0.46和0.14，两个因子解释了六个心理学测试60%的方差。

结果解读：阅读和词汇在第一因子上载荷较大，画图、积木图案和迷宫在第二因子上载荷较大，非语言的普通智力测量在两个因子上载荷较为平均，这表明存在一个语言智力因子和一个非语言智力因子。

正交旋转和斜交旋转的不同之处。

对于正交旋转，因子分析的重点在于因子结构矩阵（变量与因子的相关系数），而对于斜交旋转，因子分析会考虑三个矩阵：因子结构矩阵、因子模式矩阵和因子关联矩阵。

因子模式矩阵即标准化的回归系数矩阵。它列出了因子预测变量的权重。因子关联矩阵即因子相关系数矩阵。

图形解读：词汇和阅读在第一个因子（PA1）上载荷较大，而积木图案、画图和迷宫在第二个因子（PA2）上载荷较大。普通智力测验在两个因子上较为平均。

与可精确计算的主成分得分不同，因子得分只是估计得到的。它的估计方法有多种，fa()函数使用的是回归方法。

R包含了其他许多对因子分析非常有用的软件包。FactoMineR包不仅提供了PCA和EFA方法，还包含潜变量模型。它有许多此处我们并没考虑的参数选项，比如数值型变量和类别型变量的使用方法。FAiR包使用遗传算法来估计因子分析模型，它增强了模型参数估计能力，能够处理不等式的约束条件，GPArotation包则提供了许多因子旋转方法。最后，还有nFactors包，它提供了用来判断因子数目的许多复杂方法。

主成分分析

1.数据导入

数据结构：对10株玉米进行了生物学性状考察，考察指标有株高，穗位，茎粗，穗长，秃顶，穗粗，穗行数，行粒数。

结果解读：选择2个主成分即可保留样本大量信息。

3.提取主成分

结果解读：主成分1可解释44%的方差，主成分2解释了26%的方差，合计解释了70%的方差。

4.获取主成分得分

5.主成分方程

PC1 = 0.27 株高 - 0.04 穗位 + 0.29 茎粗 - 0.01 穗长 - 0.21 秃顶 - 0.13 穗粗 + 0.16 穗行数 + 0.24 行粒数

PC2 = -0.01 株高 + 0.36 穗位 - 0.10 茎粗 + 0.41 穗长 - 0.08 秃顶 + 0.43 穗粗 - 0.15 穗行数 + 0.01 行粒数

图形解读：此图反映了变量与主成分的关系，三个蓝点对应的RC2值较高，点上的标号2，4，6对应变量名穗位，穗长，穗粗，说明第2主成分主要解释了这些变量，与这些变量相关性强；黑点分别对应株高，茎粗，穗行数，行粒数，说明第一主成分与这些变量相关性强，第一主成分主要解释的也是这些变量，而5号点秃顶对于两个主成分均没有显示好的相关性。

因子分析

图解：可以看到需要提取4个因子。

2.提取因子

结果解读：因子1到4解释了80%的方差。

3.获取因子得分

图解：可以看出，因子1和因子2的相关系数为0.4，行粒数，株高，茎粗，秃顶在因子1的载荷较大，穗长，穗位在因子2上的载荷较大；因子3只有穗行数相关，因子4只有穗粗相关。

参考资料：

---