install.packages(c("psych", "Rcmdr","FactoMineR","GPArotation"),contriburl = "http://cran.stat.ucla.edu/bin/windows/contrib/2.6/")  #Not needed if you already ran this for the notes example
library(pcaMethods)  #Change the repository to include Bioconductor (BioC Software) where this package can be found. Only if you want to estimate missing
library(psych)

data(bfi)  #load the data
bfimod <- pca(bfi, nPcs = 5, method = "bpca")  #estimate missing values
newbfi= data.frame(bfimod@completeObs)  #create completed dataset
bartlett.test(newbfi)
FAmod2a=factor.pa(newbfi,5, rotate="varimax") #Perform standard fair MaxLik factor analysis
FAmod2a
FAmod2b=factor.pa(newbfi,5, rotate="oblimin") #Allow correlated factors
FAmod2b

cluster.plot(FAmod2a)
cluster.plot(FAmod2b)

#compare to previous dendrogram from cluster analysis


#test for hierarchical structure
FAmod2c=omega(newbfi,5, pc = "pa", digits=2) 
FAmod2c


#Omega (hierarchical) is an estimate of the general factor saturation of a test. A higher value would indicate a general factor
#Alpha is Cronbach's
#loadings are then provided for the hierarchical factor g and the rest, with default suppressing those less than .1;
#The next to last column is the sum of the squared loadings (even if not visible), i.e. the communality, and the last is the shared variance with a unique factor, i.e. 1-communality
#After that you get an orthogonal rotation solution w/o g
#Next is an oblique w/o the g factor
#Inter-factor correlations
#Oblique factor loadings on the g factor

#Omega Total (not available in R 2.6, but is in newer versions) an estimate of the total reliability of a test
#G.6 (not available in R 2.6, but is in newer versions) is one of Guttman's measures, the amount of variance in each item that can be accounted for the linear regression of all of the other items (the squared multiple correlation or smc)
#omega.graph(FAmod1b)  #Can't run as it requires another graphing program.  Just letting you know it's available.




#Final example if time to better demonstrate hierarchical structure. Library data imported from SPSS dataset via R-commander with the name Library

library=read.table("http://www.unt.edu/rss/class/mike/data/library.txt",header=T)
libdat=library[,3:14]
FAmod3=factor.pa(libdat, nfactors=3, rotate = "promax")  #how might you interpret the factors?
FAmod3
cluster.plot(FAmod3)

fa=omega(libdat, 3, pc = "pa")
fa  #looks like we do have a general factor of 'quality

#questions
#1 Willingness to help others
#2 Giving users individual attention
#3 Employees who deal with users in a caring fashion
#4 Employees who are consistently courteous

#5 A haven for quiet and solitude
#6 A meditative place
#7 A contemplative environment
#8 Space that facilitates quiet study

#9 Comprehensive print collections
#10 Complete runs of journal titles
#11 Interdisciplinary library needs being addressed
#12 Timely document deliver/interlibrary loan