# load 24 psych test
psychtest<-read.csv("24psychtests.csv",header = TRUE)
# choose "Grant-White school"
psychtestGRANT <- subset(psychtest,select = c(V1:V2, V5:V26), group == "GRANT")
names(psychtestGRANT)[23]<-paste("V3")
names(psychtestGRANT)[24]<-paste("V4") 

# Maximum Likelihood Factor Analysis
# entering raw data and extracting 3 factors
fit <- factanal(psychtestGRANT, 3, rotation="none")
print(fit, digits=2, cutoff=.3, sort=FALSE)
p = 24
k = 3
df = p*(p+1)/2 - p*k -p + k*(k-1)/2
# three factors - rejected (small p-value) 
fit <- factanal(psychtestGRANT, 5, rotation="none")
print(fit, digits=2, cutoff=.3, sort=FALSE)

# plot factor 1 by factor 2 
load <- fit$loadings[,1:2] 
par(mfrow=c(2,2))
plot(load,type="n",main="MLE, no rotation") # set up plot 
text(load,labels=names(psychtestGRANT),cex=.7)

# with varimax rotation 
fit <- factanal(psychtestGRANT, 5, rotation="varimax")
print(fit, digits=2, cutoff=.3, sort=FALSE)
# plot factor 1 by factor 2 
load <- fit$loadings[,1:2] 
plot(load,type="n",main="MLE,varimax") # set up plot 
text(load,labels=names(psychtestGRANT),cex=.7) # add variable names


#Principal Component Factor Estimation
pc<-prcomp(psychtestGRANT,scale=TRUE) #correlation PCA
sum(pc$sdev > 1)
pc$sdev
# five factors to be used (Kaiser's rule)
pcfactorloadings <- pc$rotation[,1:5]
print(pc$rotation[,1:5],digits = 2, cutoff=.3)
load <- pcfactorloadings[,1:2]
plot(load,type="n",main="PCF,no rotation") # set up plot 
text(load,labels=names(psychtestGRANT),cex=.7) # add variable names
 

rot<-varimax(pcfactorloadings, normalize = FALSE)
print(rot$loadings,digits = 2, cutoff=.3)

load <- rot$loadings[,1:2] 
plot(load,type="n",main="PCF,varimax") # set up plot 
text(load,labels=names(psychtestGRANT),cex=.7) # add variable names