##数据预处理
df<-read.csv(file="重庆板块数据(20210210).csv")         
dim(df)
head(df)
X<-df[,-1]
X1<-scale(X)
p<-ncol(X)
n<-nrow(X)

#数据可视化
#install.packages("aplpack")
library(aplpack)
rownames(X1)<-df[,1]
faces(X1,face.type = 2)#脸谱图

##主成分分析
#install.packages("psych")
library(psych)
#主成分个数确定
fa.parallel(X,fa="pc",n.iter=100,show.legend=FALSE,xlab="主成分数")
abline(h=1)
mtext("主成分数", side=1,line=2)
mtext("主成分的特征值", side=2,line=2)

#主成分
pca.fit<-principal(X,nfactors=4,rotate="none",score = TRUE)
pca.fit
par(mfrow=c(1,2))
biplot(pca.fit,choose = 1:2,main="")
biplot(pca.fit,choose = 3:4,main="")
#因子载荷
pca.fit$loadings 
-pca.fit$loadings%*%diag(1/sqrt(pca.fit$values[1:4])) #princomp函数得到的因子载荷矩阵
#主成分得分
pca.fit$scores
X1%*%pca.fit$loadings%*%diag(1/pca.fit$values[1:4])   #计算方法
#综合得分
f.score<-pca.fit$scores%*%pca.fit$Vaccounted[2,]
#排名
result<-data.frame(pca.fit$scores,f.score)
rownames(result)<-df[,1]
result<-result[order(result$f.score,decreasing = TRUE),]
result$order<-1:n
round(result,4)


##主成分旋转
rpca.fit<-principal(X,nfactors=4,rotate="varimax",score = TRUE)
rpca.fit
par(mfrow=c(1,2))
biplot(rpca.fit,choose = 1:2,main = "")
biplot(rpca.fit,choose = 3:4,main="")
#因子载荷
rpca.fit$loadings
#主成分得分
rpca.fit$scores
#综合得分
rf.score<-rpca.fit$scores%*%rpca.fit$Vaccounted[2,]
#排名
rresult<-data.frame(rpca.fit$scores,rf.score)
rownames(rresult)<-df[,1]
rresult<-rresult[order(rresult$rf.score,decreasing = TRUE),]
rresult$order<-1:n
round(rresult,4)


##稀疏主成分分析
#install.packages("elasticnet")
library(elasticnet)
Omiga<-function(X,K0,lambda1,lambda2){
  N1<-length(lambda1)
  N2<-length(lambda2)
  p<-ncol(X)
  omiga<-matrix(0,N1,N2)
  for (i in 1:N1) {
    for (j in 1:N2){
      spca.fit<-spca(cor(X),K=K0,type="Gram",sparse="penalty",lambda = lambda1[i],para=rep(lambda2[j],K0))
      if(length(which(rowSums(spca.fit$loadings)!=0))==p){
        omiga[i,j]<-sum(spca.fit$pev)+length(which(spca.fit$loadings==0))/(p*K0)
      }
    }
  }
  return(omiga)
}
lambda1<-c(0,0.001,0.01,0.1,0.5,1,2,5,10)
lambda2<-c(0,0.001,0.01,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,2,3,4,5)
omiga<-Omiga(X,4,lambda1,lambda2)
ij<-which(omiga==max(omiga),arr.ind = TRUE)
lamda1<-lambda1[ij[1]]
lamda2<-lambda2[ij[2]]
spca.fit<-spca(cor(X),K=4,type="Gram",sparse="penalty",trace=TRUE,lambda=lamda1,para=rep(lamda2,4))
spca.fit$pev
sum(spca.fit$pev)
spca.fit$loadings

#转化为principal函数一致的主成分得分
score<-X1%*%spca.fit$loadings%*%diag(1/sqrt(pca.fit$values[1:4]))  
#综合得分
sf.score<-score%*%spca.fit$pev
#排名
sresult<-data.frame(score,sf.score)
rownames(sresult)<-df[,1]
sresult<-sresult[order(sresult$sf.score,decreasing = TRUE),]
sresult$order<-1:n
round(sresult,4)