############################################################
#Fabrigar, L.R., Wegener, D.T., MacCallum, R.C., & Strahan, E.J. (1999).
#Evaluating the use of exploratory factor analysis in psychological research.
#Psychological Methods, 4, 272-299.
############################################################
library(rpsychi)
library(GPArotation)
library(psych)
library(nFactors)

#-----------------------------------#
#Breckler (1984, Study 1)
#-----------------------------------#

##(1) データの作成
x <- as.matrix(read.table("http://blue.zero.jp/yokumura/R/factor/breckler1.csv", header=F, sep=","))
x <- x[lower.tri(x)]
vnam <- c("Thurstone affect", "Mood Checklist (+)", "Mood Checklist (-)",
                                      "Action sequence", "Distance", "Thurstone behavior",
                                      "Thurstone cognition", "Semantic differential", "Listed thougts")
dat <- lower2R(x, varname=vnam)


##(2) 共通性
tab1 <- smc(dat)
tab1
range(tab1)
mean(tab1)


##(3) 固有値
dat.red <- dat
diag(dat.red) <- tab1
tab2 <- rbind(eigen(dat)$value, eigen(dat.red)$value)
round(tab2, 2)


##(4) 適合度指標
n <- 138
fit1 <- fa(r=dat, n.obs=n, nfactors=1, fm="ml")
fit2 <- fa(r=dat, n.obs=n, nfactors=2, fm="ml")
fit3 <- fa(r=dat, n.obs=n, nfactors=3, fm="ml")
fit4 <- fa(r=dat, n.obs=n, nfactors=4, fm="ml")

tab3 <- rbind(fit1$RMSEA, fit2$RMSEA, fit3$RMSEA, fit4$RMSEA)
round(tab3, 2)


##(5) 主成分負荷量
pa.unrotated <- fa(r=dat, n.obs=n, nfactors=3, fm="pa", rotate="none", max.iter=1, SMC=FALSE)
pa.or <- GPForth(loadings(pa.unrotated), method="varimax", normalize=TRUE)
pa.ob <- GPFoblq(loadings(pa.unrotated), method="quartimin", normalize=TRUE)
ord <- c(3,2,1)
round(loadings(pa.or)[,ord],2)    #Varimax
ord <- c(2,3,1)
round(loadings(pa.ob)[,ord],2)    #Direct quartimin
round(pa.ob$Phi[ord,ord],2)       #component correlation


##(6) 因子パタン
fa.unrotated <- fa(r=dat, n.obs=n, nfactors=3, fm="ml", rotate="none")
fa.or <- GPForth(loadings(fa.unrotated), method="varimax", normalize=TRUE)
fa.ob <- GPFoblq(loadings(fa.unrotated), method="quartimin", normalize=TRUE)
ord <- c(1,3,2)
round(loadings(fa.or)[,ord],2)    #Varimax
round(loadings(fa.ob)[,ord],2)    #Direct quartimin
round(fa.ob$Phi[ord,ord],2)       #factor correlation



#-----------------------------------#
#Breckler (1984, Study 2)
#-----------------------------------#
x <- as.matrix(read.table("http://blue.zero.jp/yokumura/R/factor/breckler2.csv", header=F, sep=","))
x <- x[lower.tri(x)]
vnam <- c("Thurstone affect", "Mood Checklist (+)", "Mood Checklist (-)",
                                      "Action sequence", "Distance", "Thurstone behavior",
                                      "Thurstone cognition", "Semantic differential", "Listed thougts")
dat <- lower2R(x, varname=vnam)


##(2) 共通性
tab1 <- smc(dat)
tab1
range(tab1)
mean(tab1)


##(3) 固有値
dat.red <- dat
diag(dat.red) <- tab1
tab2 <- rbind(eigen(dat)$value, eigen(dat.red)$value)
round(tab2, 2)


##(4) 適合度指標
n <- 105
fit1 <- fa(r=dat, n.obs=n, nfactors=1, fm="ml")
fit2 <- fa(r=dat, n.obs=n, nfactors=2, fm="ml")
fit3 <- fa(r=dat, n.obs=n, nfactors=3, fm="ml")
fit4 <- fa(r=dat, n.obs=n, nfactors=4, fm="ml")

tab3 <- rbind(fit1$RMSEA, fit2$RMSEA, fit3$RMSEA, fit4$RMSEA)
round(tab3, 2)


##(5) 主成分負荷量
pa.unrotated <- fa(r=dat, n.obs=n, nfactors=3, fm="pa", rotate="none", max.iter=1, SMC=FALSE)
pa.or <- GPForth(loadings(pa.unrotated), method="varimax", normalize=TRUE)
pa.ob <- GPFoblq(loadings(pa.unrotated), method="quartimin", normalize=TRUE)
ord <- c(2,1,3)
round(loadings(pa.or)[,ord],2)    #Varimax
round(loadings(pa.ob)[,ord],2)    #Direct quartimin
round(pa.ob$Phi[ord,ord],2)       #component correlation


##(6) 因子パタン
fa.unrotated <- fa(r=dat, n.obs=n, nfactors=3, fm="ml", rotate="none")
fa.or <- GPForth(loadings(fa.unrotated), method="varimax", normalize=TRUE)
fa.ob <- GPFoblq(loadings(fa.unrotated), method="quartimin", normalize=TRUE)
ord <- c(3,1,2)
round(loadings(fa.or)[,ord],2)    #Varimax
round(loadings(fa.ob)[,ord],2)    #Direct quartimin
round(fa.ob$Phi[ord,ord],2)       #factor correlation