# PCA
dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/AMSA_track_mens.csv",row.names=1)
plot(dat,pch=20)
dat.pca<-princomp(dat)
summary(dat.pca)
screeplot(dat.pca,type='l',lwd=2,pch=20,cex=4,,col='red')
abline(h=sum(summary(dat.pca)[[1]])/8,lty=2,lwd=2)

# factor analysis
dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/FA01.csv")
dat.fa<-factanal(dat,1)
dat.fa$uniquenesses
1-dat.fa$loadings[1:3]^2

dat<-read.table("http://www.matsuka.info/data_folder/tdkPCA01.txt")
dat.pca<-princomp(dat)

dat.fa<-factanal(dat,1,score="regression")
plot(dat.fa$score~dat.pca$score[,1],pch=20,cex=2,xlab="Component Score", ylab="Factor Score")
cor(dat.fa$score,dat.pca$score)
cor(rowSums(dat),dat.fa$score)

# plotting results
dat<-read.csv("http://www.matsuka.info/data_folder/tdkCFA.csv")
dat.fa<-factanal(dat,2,score="regression")
plot(dat.fa$loadings[,1],dat.fa$loadings[,2],pch=20,cex=2,col="red",
     xlab="Factor 1", ylab="Factor 2")
abline(h=0,v=0)
n.var = length(dat.fa$loadings[,1])
for (i.var  in 1:n.var){
  lines(c(0,dat.fa$loadings[i.var,1]),c(0,dat.fa$loadings[i.var,2]),
        lty=2,lwd=2,col='skyblue')
}
text(dat.fa$loadings[,1],dat.fa$loadings[,2],colnames(dat),cex=1.4)

# model comparison
dat.model1<-factanal(dat,1)
dat.model2<-factanal(dat,2)
dat.model3<-factanal(dat,3)
dat.model4<-factanal(dat,4)
print(m1res<-c(dat.model1$dof,dat.model1$PVAL))
print(m2res<-c(dat.model2$dof,dat.model2$PVAL))
print(m3res<-c(dat.model3$dof,dat.model3$PVAL))
print(m4res<-c(dat.model4$dof,dat.model4$PVAL))
1-pchisq(dat.model3$STATISTIC-dat.model4$STATISTIC,dat.model3$dof-dat.model4$dof)
dat.model2$STATISTIC
dat.model2$dof

# CFA with SEM
install.packages('sem')
library(sem)

model01=cfa(reference.indicator=FALSE)
F1:extrovert,cheerful, leadership, antisocial, talkative, motivated, hesitance, popularity

mod1<-sem(model01, cov(dat), 100)

model02=cfa(reference.indicator=FALSE)
F1: extrovert, leadership, motivated, hesitance 
F2: cheerful, antisocial, talkative, popularity

mod2<-sem(model02, cov(dat), 100)

opt <- options(fit.indices = c("RMSEA"))