####因子分析####
####主成分法
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)
}
####例9.7
x<-c(1.000,
0.923, 1.000,
0.841, 0.851, 1.000,
0.756, 0.807, 0.870, 1.000,
0.700, 0.775, 0.835, 0.918, 1.000,
0.619, 0.695, 0.779, 0.864, 0.928, 1.000,
0.633, 0.697, 0.787, 0.869, 0.935, 0.975, 1.000,
0.520, 0.596, 0.705, 0.806, 0.866, 0.932, 0.943, 1.000)
names<-c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8")
R<-matrix(0,ncol=8,nrow=8,dimnames=list(names,names))
for(i in 1:8){
for(j in 1:i){
R[i,j]<-x[(i-1)*i/2+j];R[j,i]<-R[i,j]
}
}
R
fa<-factor.analy1(R,m=2);fa
E<- R-fa$loadings %*% t(fa$loadings)-diag(fa$var[,2])
sum(E^2)
####主因子法
factor.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
}
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)
}
####極大似然法
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)
}
####三法合一
factor.analy<-function(S, m=0,
d=1/diag(solve(S)), method="likelihood"){
if (m==0){
p<-nrow(S); 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)
)
}
####因子分析的計算函數####
####factnanal()函數
factanal(x, factors, data = NULL, covmat = NULL, n.obs = NA,
subset, na.action, start = NULL,
scores = c("none", "regression", "Bartlett"),
rotation = "varimax", control = NULL, ...)
#factors:因子個(gè)數;comvat:樣本的協(xié)方差矩陣或相關(guān)矩陣；scores:因子得分方法；rotation：因子旋轉
####例9.11
x<-c(1.000,
0.923, 1.000,
0.841, 0.851, 1.000,
0.756, 0.807, 0.870, 1.000,
0.700, 0.775, 0.835, 0.918, 1.000,
0.619, 0.695, 0.779, 0.864, 0.928, 1.000,
0.633, 0.697, 0.787, 0.869, 0.935, 0.975, 1.000,
0.520, 0.596, 0.705, 0.806, 0.866, 0.932, 0.943, 1.000)
names<-c("X1", "X2", "X3", "X4", "X5", "X6", "X7", "X8")
R<-matrix(0,ncol=8,nrow=8,dimnames=list(names,names))
for(i in 1:8){
for(j in 1:i){
R[i,j]<-x[(i-1)*i/2+j];R[j,i]<-R[i,j]
}
}
fa<-factanal(factors=2,covmat=R);fa
####例9.12
rt<-read.table("E:/Study/研究生/R語(yǔ)言/統計建模與Ｒ軟件/R原程序/Ch03/applicant.data")
fa<-factanal(~., factors=5, data=rt, scores="Bartlett")
fa<-factanal(~., factors=5, data=rt, scores="regression")
plot(fa$scores[, 1:2], type="n"); text(fa$scores[,1], fa$scores[,2])
plot(fa$scores[, c(1,3)], type="n"); text(fa$scores[,1], fa$scores[,3])
