##数据预处理
df<-read.csv(file="新冠检测概念数据(20210219).csv")         
dim(df)
head(df)
x<-df[,1]
X<-df[,-1]
X1<-scale(X)
p<-ncol(X)
n<-nrow(X)
R<-cor(X)

#数据可视化
color <- c('#8DD3C7', '#FFFFB3', '#BEBADA', '#FB8072', '#80B1D3', '#FDB462', '#B3DE69', '#FCCDE5', '#BC80BD', '#CCEBC5', 'gray') 
stars(X1,draw.segments = TRUE,col.segments = color[1:11], key.loc= c(22,10),labels = x)

#因子分析检验
#install.packages("psych")
library(psych)
cortest.bartlett(R,n)

library(MASS)
kmo<-function(data){
  x <- cor(as.matrix(data))
  ix <- ginv(x)
  S2 <- diag(diag((ix^-1)))
  AIS <- S2%*%ix%*%S2                      # anti-image covariance matrix
  IS <- x+AIS-2*S2                         # image covariance matrix
  Dai <- sqrt(diag(diag(AIS)))
  IR <- ginv(Dai)%*%IS%*%ginv(Dai)         # image correlation matrix
  AIR <- ginv(Dai)%*%AIS%*%ginv(Dai)       # anti-image correlation matrix
  a <- apply((AIR - diag(diag(AIR)))^2, 2, sum)
  AA <- sum(a)
  b <- apply((x - diag(nrow(x)))^2, 2, sum)
  BB <- sum(b)
  MSA <- b/(b+a)                        # indiv. measures of sampling adequacy
  
  AIR <- AIR-diag(nrow(AIR))+diag(MSA)  # Examine the anti-image of the
  # correlation matrix. That is the
  # negative of the partial correlations,
  # partialling out all other variables.
  
  kmo <- BB/(AA+BB)                     # overall KMO statistic
  
  # Reporting the conclusion
  if (kmo >= 0.00 && kmo < 0.50){
    test <- 'The KMO test yields a degree of common variance
unacceptable for FA.'
  } else if (kmo >= 0.50 && kmo < 0.60){
    test <- 'The KMO test yields a degree of common variance miserable.'
  } else if (kmo >= 0.60 && kmo < 0.70){
    test <- 'The KMO test yields a degree of common variance mediocre.'
  } else if (kmo >= 0.70 && kmo < 0.80){
    test <- 'The KMO test yields a degree of common variance middling.'
  } else if (kmo >= 0.80 && kmo < 0.90){
    test <- 'The KMO test yields a degree of common variance meritorious.'
  } else {
    test <- 'The KMO test yields a degree of common variance marvelous.'
  }
  
  ans <- list(  overall = kmo,
                report = test,
                individual = MSA,
                AIS = AIS,
                AIR = AIR )
  return(ans)
  
}    
kmo(X)$overall

##因子分析
#因子个数的确定
fa.parallel(X,n.obs = n,fa="fa",n.iter=100,show.legend=FALSE)
abline(h=0,lty=2)
mtext("因子数", side=1,line=2)
mtext("主成分因子特征值", side=2,line=2)

#极大似然法提取公因子
fa.fit<-fa(X, nfactors=2, rotate = "none", fm = "ml", scores="regression" )  #改变fm还可以采用主轴迭代,加权最小二乘等方法
fa.fit  
#与factanal函数得到的有区别,这里是Standardized loadings (pattern matrix) based upon correlation matrix

##因子正交旋转
fa.varimax <- fa(X, nfactors = 2, rotate = "varimax",fm = "ml", scores="regression" )
fa.varimax

#因子斜交旋转
fa.promax <- fa(X, nfactors = 2, rotate = "Promax",fm = "ml", scores="regression" )
fa.promax  
#注意斜交fa.promax$loading是标准化的回归系数,并非是相关系数.
#计算斜交的因子载荷矩阵
fsm <- function(oblique) {
  if (class(oblique)[2]=="fa" & is.null(oblique$Phi)) {
    warning("Object doesn't look like oblique EFA")
  } else {
    P <- unclass(oblique$loading)
    F <- P %*% oblique$Phi
    colnames(F) <- c("PA1", "PA2")
    return(F)
  }
}
fsm(fa.promax)

par(mfrow=c(1,2))
factor.plot(fa.varimax,labels = rownames(fa.promax$loadings),
            pch=14,main="",xlim=c(-0.3,1),ylim=c(-0.4,0.9),xlab="ML2(正交)")
factor.plot(fa.promax,labels = rownames(fa.promax$loadings), 
            pch=14,main="",xlim=c(-0.4,1),ylim=c(-0.4,0.9),xlab="ML2(斜交)")

fa.diagram(fa.varimax,simple =FALSE,main = "正交",digits=2)
fa.diagram(fa.promax,simple = FALSE,main = "斜交",digits=2)

par(mfrow=c(1,2))
plot(fa.varimax$scores[,2],fa.varimax$scores[,1],pch=3,cex=0.5,
     xlab = "因子1(正交)",ylab = "因子2")
text(fa.varimax$scores[,2],fa.varimax$scores[,1],labels(df[,1]))

plot(fa.promax$scores[,2],fa.promax$scores[,1],pch=4,cex=0.5,
     xlab = "因子1(斜交)",ylab = "因子2")
text(fa.promax$scores[,2],fa.promax$scores[,1],labels(df[,1]))

##稀疏因子分析
#install.packages("fanc")
library(fanc)
X<-as.matrix(X)
sfa.fit<-fanc(X,factors=2,type="MC")
sfa.fit
plot(sfa.fit)
rho<-select(sfa.fit,criterion="BIC",gamma=1.96)$rho
rho

#因子载荷
sfa.out<-out(sfa.fit,rho,gamma=1.96)
sfa.out
sfa.out$loadings

#回归法计算稀疏因子得分
A<-as.matrix(sfa.out$loadings)
B<- t(A) %*%solve(R)
sfa.scores<-X1%*%t(B)
sfa.scores

#因子得分
scores<-data.frame(df[,1],fa.fit$scores,fa.varimax$scores,fa.promax$scores,sfa.scores)
scores

plot(sfa.scores,pch=3,cex=0.5, xlab = "因子1",ylab = "因子2")
text(sfa.scores,labels(df[,1]))