## chisq test
chisq.test(c(72,23,16,49),p=rep(40,4),rescale.p=T)
M=matrix(c(52,48,8,42),nrow=2)
chisq.test(M,correct=F)

## PCA
Mdat<-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)

dist = c(100,200,400,800,1500,5000,10000,42195)
log.dist=log(dist)

plot(log.dist, dat.pca$loadings[,1],pch=20,col='red',
  cex=2,ylab="PC loadings (1st)",xlab="Log(distance)")
plot(log.dist, dat.pca$loadings[,2],pch=20,col='red',
  cex=5,ylab="PC loadings (2nd)",xlab="Log(distance)")

### EFA
dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/FA01.csv")
dat.fa<-factanal(dat,1)

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

dat<-read.csv("http://www.matsuka.info/data_folder/tdkCFA.csv")
dat.fa<-factanal(dat,2,rotation="varimax",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)
text(dat.fa$loadings[,1],dat.fa$loadings[,2],colnames(dat),cex=1.4)

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)

AIC1=2*(36-dat.model1$dof)+dat.model1$STATISTIC
AIC2=2*(36-dat.model2$dof)+dat.model2$STATISTIC
AIC3=2*(36-dat.model3$dof)+dat.model3$STATISTIC
AIC4=2*(36-dat.model4$dof)+dat.model4$STATISTIC

### 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), nrow(dat))
summary(mod1)

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"))