# reliability
install.packages("psych")
library("psych")
dat<-read.table("http://peach.l.chiba-u.ac.jp/course_folder/survey_data01.txt")
ca<-alpha(dat)
dat<-read.csv("http://peach.l.chiba-u.ac.jp/course_folder/fa2015.csv")
caALL<-alpha(dat[,1:10])
ca1_5<-alpha(dat[,1:10])
ca6_10<-alpha(dat[,1:10])

# w/ factor analysis
F1<-factanal(dat[,1:5],1)
F2<-factanal(dat[,6:10],1)
library(sem)
fa.model=cfa(reference.indicator=FALSE)
F1: q1,q2,q3,q4,q5
F2: q6,q7,q8,q9,q10

fa.model<-update(fa.model)
delete, F1<->F2

fa.result<-sem(fa.model, cov(dat), 300)

# IRT models
install.packages("ltm")
library(ltm)
dat<-read.table("http://peach.l.chiba-u.ac.jp/course_folder/survey_data01.txt")
descript(dat)
irt1P<-rasch(dat)
plot(irt1P)
lrt2P<-ltm(dat~z1)
plot(irt2P)

item.fit(lrt1P)
item.fit(lrt2P)
anova(irt1P,lrt2P)

# cluster analysis
cldata<-data.frame(var1=c(4,1,5,1,5), var2=c(1,5,4,3,1))
cldata.cluster=hclust(dist(cldata),method="average")
plot(cldata.cluster,cex=2)
cutree(cldata.cluster,2)

dat<-read.csv("http://matsuka.info/data_folder/tdkClust.csv", header=TRUE, row.names=1)
dat.cluster=hclust(dist(dat),method="average")
plot(dat.cluster,cex=1.5)
dat.pca=princomp(dat)
biplot(dat.pca)

dat.kmeans=kmeans(dat, centers=3)
plot(dat,col=dat.kmeans$cluster,pch=20,cex=1.75)

データ解析基礎論B　WAB07　因子分析・主成分分析
提出期限　2017.12.05 講義開始前まで

WAB07.1
このデータは、ある質問紙の回答（計10問）です。最初の５問で要因F１を、残りの５問は他の要因F２を測定する前提で作成された質問紙です。この前提が妥当に質問紙に反映されているか探索的因子分析(factanalを用いて）を行い検証してください。分析に先立ち、データを可視化してください。

WAB07.2
上記のデータを用いて主成分分析を行ってください。主成分分析の結果を用いて上記の前提に関して言及してください。

WAB07.3
上記の前提の妥当性を確認的因子分析(semを用いて）を行い検証してください。

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

データ分析基礎論B weekly assignment WAB05はここにあります。

# Split Plot Factorial (SFP)
source("http://peach.l.chiba-u.ac.jp/course_folder/tsme.txt")
dat<-read.csv("http://www.matsuka.info/data_folder/dktb3211.txt")
dat.aov<-aov(result~method*duration+Error(s+s:duration),dat)

TSME<-SPF.tsme(dat.aov,dat,"result")
dat.aov.sum<-summary(dat.aov)

# within-subjcts factor
DF_error.W = dat.aov.sum[[2]][[1]][3,1]
MS_error.W = dat.aov.sum[[2]][[1]][3,3]

q <- qtukey(0.95,4,DF_error.W)
hsd<-q*sqrt(MS_error.W/(5*2))
dat.means.duration<-tapply(dat$result,dat$duration,mean)
abs(outer(dat.means.duration,dat.means.duration,'-'))>hsd

# between-subjects factor
DF_error.B = dat.aov.sum[[1]][[1]][2,1]
MS_error.B = dat.aov.sum[[1]][[1]][2,3]

q <- qtukey(0.95,2,DF_error.B)
hsd<-q*sqrt(MS_error.B/(5*4))

dat.means.method<-tapply(dat$result,dat$method,mean)
abs(outer(dat.means.method,dat.means.method,'-'))>hsd

# randomized block factorial
dat<-read.csv("http://www.matsuka.info/data_folder/dktb3218.txt")
dat.aov<-aov(result~method+duration+method:duration +
  Error(s+method:s+duration:s+method:duration:s), data=dat)
TSME<-RBF.tsme(dat.aov, dat, "result")

dat.aov.sum<-summary(dat.aov)
# testing for method
mse  = dat.aov.sum[[2]][[1]][2,3]
qv=qtukey(0.95,2,4)
hsd=qv*sqrt(mse/(5*4))
dat.means.method=tapply(dat$result,dat$method,mean)
abs(outer(dat.means.method,dat.means.method,'-')) > hsd

# testing for duration
mse  = dat.aov.sum[[3]][[1]][2,3]
qv=qtukey(0.95,4,12)
hsd=qv*sqrt(mse/(5*2))
dat.means.duration=tapply(dat$result,dat$duration,mean)
abs(outer(dat.means.duration,dat.means.duration,'-')) > hsd