因子分析

factor.analy1<-function(S, m){
p<-nrow(S); diag_S<-diag(S); sum_rank<-sum(diag_S)
rowname<-paste("X", 1:p, sep="")
colname<-paste("Factor", 1:m, sep="")
A<-matrix(0, nrow=p, ncol=m,
dimnames=list(rowname, colname))
eig<-eigen(S)
for (i in 1:m)
A[,i]<-sqrt(eig$values[i])*eig$vectors[,i]
h<-diag(A%*%t(A))
rowname<-c("SS loadings","Proportion Var","Cumulative Var")
B<-matrix(0, nrow=3, ncol=m,
dimnames=list(rowname, colname))
for (i in 1:m){
B[1,i]<-sum(A[,i]^2)
B[2,i]<-B[1,i]/sum_rank
B[3,i]<-sum(B[1,1:i])/sum_rank
}
method<-c("Principal Component Method")
list(method=method, loadings=A,
var=cbind(common=h, spcific=diag_S-h), B=B)

}

source("factor.analy1.R")

fa<-factor.analy1(R, m=2); fa

E<- R-fa$loadings %*% t(fa$loadings)-diag(fa$var[,2])

> sum(E^2)

or.analy2<-function(R, m, d){
p<-nrow(R); diag_R<-diag(R); sum_rank<-sum(diag_R)
rowname<-paste("X", 1:p, sep="")
colname<-paste("Factor", 1:m, sep="")
A<-matrix(0, nrow=p, ncol=m,
dimnames=list(rowname, colname))
kmax=20; k<-1; h <- diag_R-d
repeat{
diag(R)<- h; h1<-h; eig<-eigen(R)
for (i in 1:m)
A[,i]<-sqrt(eig$values[i])*eig$vectors[,i]
h<-diag(A %*% t(A))
if ((sqrt(sum((h-h1)^2))<1e-4)|k==kmax) break
k<-k+1
}
9.2 E>\g 529
rowname<-c("SS loadings","Proportion Var","Cumulative Var")
B<-matrix(0, nrow=3, ncol=m,
dimnames=list(rowname, colname))
for (i in 1:m){
B[1,i]<-sum(A[,i]^2)
B[2,i]<-B[1,i]/sum_rank
B[3,i]<-sum(B[1,1:i])/sum_rank
}
method<-c("Principal Factor Method")
list(method=method, loadings=A,
var=cbind(common=h,spcific=diag_R-h),B=B,iterative=k)

}

d<-c(0.123, 0.112, 0.155, 0.116, 0.073, 0.045, 0.033, 0.095)
> source("factor.analy2.R")

> fa<-factor.analy2(R, m=2, d); fa

E<- R-fa$loadings %*% t(fa$loadings)-diag(fa$var[,2])

> sum(E^2)

factor.analy3<-function(S, m, d){
p<-nrow(S); diag_S<-diag(S); sum_rank<-sum(diag_S)

rowname<-paste("X", 1:p, sep="")

colname<-paste("Factor", 1:m, sep="")
A<-matrix(0, nrow=p, ncol=m,
dimnames=list(rowname, colname))
kmax=20; k<-1
repeat{
d1<-d; d2<-1/sqrt(d); eig<-eigen(S * (d2 %o% d2))
for (i in 1:m)
A[,i]<-sqrt(eig$values[i]-1)*eig$vectors[,i]
A<-diag(sqrt(d)) %*% A
d<-diag(S-A%*%t(A))
if ((sqrt(sum((d-d1)^2))<1e-4)|k==kmax) break
k<-k+1

}

rowname<-c("SS loadings","Proportion Var","Cumulative Var")
B<-matrix(0, nrow=3, ncol=m,
dimnames=list(rowname, colname))
for (i in 1:m){
B[1,i]<-sum(A[,i]^2)
B[2,i]<-B[1,i]/sum_rank
B[3,i]<-sum(B[1,1:i])/sum_rank
}
method<-c("Maximum Likelihood Method")
list(method=method, loadings=A,
var=cbind(common=diag_S-d, spcific=d),B=B,iterative=k)

}

d<-c(0.123, 0.112, 0.155, 0.116, 0.073, 0.045, 0.033, 0.095)

> source("factor.analy3.R")

fa<-factor.analy3(R, m=2, d); fa

E<- R-fa$loadings %*% t(fa$loadings)-diag(fa$var[,2])

> sum(E^2)

factor.analy<-function(S, m=0,
d=1/diag(solve(S)), method="likelihood"){
if (m==0){
p<-nrow(X); eig<-eigen(S)
sum_eig<-sum(diag(S))
for (i in 1:p){
if (sum(eig$values[1:i])/sum_eig>0.70){
m<-i; break
}
}
}
source("factor.analy1.R"); source("factor.analy2.R")
source("factor.analy3.R")
switch(method,
princomp=factor.analy1(S, m),
factor=factor.analy2(S, m, d),
likelihood=factor.analy3(S, m, d)

)

}

source("factor.analy.R")
> fa<-factor.analy(R, m=2, method="princomp")

> vm1<-varimax(fa$loadings, normalize = F); vm1

> fa<-factanal(factors=2, covmat=R); fa

> fa<-factanal(factors=2, covmat=R); fa

t<-read.table("applicant.data")

> fa<-factanal(~., factors=5, data=rt, scores="regression")

plot(fa$scores[, 1:2], type="n")
> text(fa$scores[,1], fa$scores[,2])

      w=read.csv("LA.Neighborhoods.csv")

 w$density=w$Population/w$Area
    

      u=w[,-c(1,12:15)]

(a=factanal(factors=2,scale(u[,-1]),scores="regression"))
     (a$loadings) 
      (a$scores)
      
  plot(a$loadings[,1:2],type="n",ylim=c(-1.1,1.1),xlim=c(-1.1,1.2),xlab="factor1",ylab="factor2",main="loadings")    
      
      abline(h=0)
      abline(v=0)
      text(a$loadings[,1],a$loadings[,2],labels=row.names(a$loadings),cex=1)
      
      plot(a$scores[,1:2],type="n",ylim=c(-2,1.5),xlim=c(-2.5,2.5),xlab="factor1",ylab="factor2",main="factor scores")
            abline(h=0)
      abline(v=0)

      text(a$scores[,1],a$scores[,2],labels=u[,1],cex=0.7)

      

factanal(pgd,8)

model <- factanal(pgd,8)
windows(7,7)
par(mfrow=c(2,2))
plot(loadings(model)[,1],loadings(model)[,2],pch=16,xlab="Factor 1",
ylab="Factor 2",col="blue")
plot(loadings(model)[,1],loadings(model)[,3],pch=16,xlab="Factor 1",
ylab="Factor 3",col="red")
plot(loadings(model)[,1],loadings(model)[,4],pch=16,xlab="Factor 1",
ylab="Factor 4",col="green")
plot(loadings(model)[,1],loadings(model)[,5],pch=16,xlab="Factor 1",
ylab="Factor 5",col="brown")