MSA课程作业代码记录

比较主成分分析和因子分析的步骤和异同，选用的是R语言自带的ISLR包中的College数据集，菜鸟萌新，代码部分来自老师指导和网络，这里只是个课程作业记录，侵删

library(ISLR)
library(psych)
data4<-data.frame(College)
head(data4)
summary(data4)
主成分分析
#KMO检验
data4<-data4[,-1]
KMO(cor(data4))
#Bartlett检验
cortest.bartlett(cor(data4), n=777)
#cattell碎石检验#
fa.parallel(data4,fa=“pc”,n.iter = 100,show.legend = T,main=“Cattell碎石检验”,ylabel=“特征值”)
abline(2,0)#发现有2个主成分#
#提取主成分#
pc2<-principal(data4,nfactors=2)
pc2
#作出主成分散点图#
data4.pr <- princomp(data4, cor = TRUE)
biplot(data4.pr)
#获取每个变量在主成份上的得分#
pc<-principal(data4,nfactors=2,score=TRUE)
head(pcKaTeX parse error: Expected 'EOF', got '#' at position 9: scores) #̲获取主成分得分系数 round…weights),2)
因子分析
#测试数据整理#
options(digits=2)
data4cov<-as.matrix(data4)
n<-nrow(data4cov)
mx<-diag(1,n)-matrix(1,n,n)/n
covA<-t(data4cov)%%mx%%data4cov/(n-1);covA
covariances<-covA
covariances
correlations<-cov2cor(covariances)
correlations
#极大似然因子分析
fa.parallel(correlations,n.obs=777,fa=“both”,n.iter=100,main=“Scree plots with parallel analysis”)
#提取3个因子，使用最大方差法旋转
fit <- factanal(data4, 2, rotation=“varimax”) # 第2个参数是提取的因子个数
print(fit, digits=2, cutoff=0.3, sort=TRUE) # 输出结果
fa<-fa(correlations,nfactors=2,rotate=“none”,fm=“pa”)
fa
#用正交旋转提取因子
fa.varimax<-fa(correlations,nfactors=2,rotate=“varimax”,fm=“pa”)
fa.varimax
#用斜交旋转提取因子
install.packages(“GPArotation”)
library(GPArotation)
fa.promax<-fa(correlations,nfactors=2,rotate=“promax”,fm=“pa”)
fa.promax
#绘制正交结果的图形
factor.plot(fa.varimax,labels=rownames(fa.varimaxKaTeX parse error: Expected 'EOF', got '#' at position 47: …x,simple=TRUE) #̲绘制斜交结果的图形 facto…loadings))
fa.diagram(fa.promax,simple=TRUE)
#求因子的相关系数矩阵
fa <- factanal(~., factors =2, data = data4, score = “regression”, rotation = “varimax”)
faKaTeX parse error: Expected 'EOF', got '#' at position 13: correlation #̲求因子得分 fa.varima…weights

箭头图
library(ISLR)
library(psych)
data4<-data.frame(College)
data4<-data4[,-1]

rm(list=ls())
library(ggplot2)
pca <- prcomp(data4,scale. = T,rank=2,retx=T)
#取PC1和PC2的比例
xlab <- paste(“PC1”,"(",round((summary(pca))$importance[2,1]\ast 100,1),"ylab<-paste("PC2","(",round((summary(pca))$importance[2,2]*100,1),"%)",sep="")
x<-“PC1”
y<-“PC2”
data_x <- data.frame(varnames=rownames(pca$x),pca$x) #为方便取用数据，将pca结果放在一个数据库里面
plot_1 <- ggplot(data_x, aes(PC1,PC2))+geom_point(aes(color=varnames),size=3)+
coord_equal(ratio=1)+xlab(xlab)+ylab(ylab) # 这里先画出点，coord_equal(ratio=1) 将X轴和y轴比例设置为一样的
data_rotation <- data.frame(obsnames=row.names(pca$rotation),pca$rotation)

#获取箭头缩放比例
mult <- min(
(max(data_x[,y]) - min(data_x[,y])/(max(data_rotation[,y])-min(data_rotation[,y]))),
(max(data_x[,x]) - min(data_x[,x])/(max(data_rotation[,x])-min(data_rotation[,x])))
)

#设置箭头坐标
data_2 <- transform(data_rotation,
v1 = mult * (get(x)),
v2 = mult * (get(y))
)
#添加箭头
plot_1<-plot_1+geom_segment(data=data_2,aes(x=0,y=0,xend=v1,yend=v2),arrow=arrow(length=unit(0.2,“cm”)),alpha=0.75)
#添加箭头名称
plot_1<-plot_1+geom_text(data=data_2,aes(v1,v2,label=obsnames),size=3,nudge_x=-0.05,nudge_y=-0.01)
#对图形结果进行修饰
plot_1 <- plot_1+scale_color_discrete(guide=guide_legend(title=“stage type”))+theme_bw()+
theme(plot.background=element_blank(),panel.background=element_blank(),
panel.grid.minor=element_blank(),panel.grid.major=element_blank(),
axis.title=element_text(color=“black”,size=15),
axis.text=element_text(size=15))+guides(color=F)
plot_1