# R code for chapter 64 Factor analysis
# Book details:
# Health Science Statistics using R and R Commander
# by Robin Beaumont
# Paperback - ISBN 9781907904318
# http://www.amazon.co.uk/Health-Science-Statistics-using-Commander/dp/190790431X
# or direct from publisher at: 
# http://www.scionpublishing.com/
# book website: http://www.robin-beaumont.co.uk/rbook
###########################################################
#
hozdata<- read.table("http://www.robin-beaumont.co.uk/virtualclassroom/book2data/grnt_fem.dat",sep="\t", header=TRUE)
names(hozdata)

hozdatamatrix <- cor(hozdata)
hozdatamatrix
round(hozdatamatrix,  4) 

install.packages("psych")
library(psych)
describe(hozdata)

cortest.bartlett(hozdata)
kmo_result <- KMO(hozdata)
kmo_result
kmo_result$ImCov   # same as $Ais matrix in book
kmo_result$Image   # non-diagonals same as $Air matrix in book

det(hozdatamatrix) 
result1 <- eigen(hozdatamatrix)
result1

model1<- principal(hozdata, nfactors = 6, rotate = "none")
model1
plot(model1$values, type = "b")

# back to the result1 object to demonstrate
# how the model1 results are from matrix multiplication:
diag(result1$values) # adds the eigenvalues as the diagonals in a square matrix
# multiple the eigenvector matrix with the  square root of the  matrix containing 
# the eigenvalues along the diagonals. 
loadings_me <- result1$vectors %*% sqrt(diag(result1$values))
loadings_me
# to see source code to psych::principal type:
# psych::principal
# now we know we only want two factors create a new model:
model2 <- principal(hozdata, 
nfactors = 2, rotate = "none")
model2 # now returns RMSR and chi to test hypothesis 2 factors sufficient
factor.model(model2$loadings)# model produced correlations and communalities(diagonals) 
# now difference between model implied matrix and that from the raw data:
residuals <- factor.residuals(hozdatamatrix, model2$loadings)
residuals
#################
# 2017 addition:
cluster.plot(model1, cut= 0.5)
plot(model1,show.points=TRUE)
biplot(model1)
biplot(model2,labels=rownames(hozdata))
biplot(model2)
# can also produce plot with GGplot2 available on github: ggbiplot
# The points correspond to the prinipal component scores of each observation. The white circle indicates the theoretical maximum extent of the arrows. Also possible to add ellipses if you add a categorisation variable.
library(devtools)
install_github("ggbiplot", "vqv")
library(ggbiplot)
model1.pca <- prcomp(hozdata, scale. = TRUE)
model1.pca
g <- ggbiplot(model1.pca, obs.scale = 1, var.scale = 1, ellipse = TRUE, circle = TRUE)
print(g)

#################
# now factor analysis proper
model4 <- fa(hozdata, nfactors = 2,fm = "pa", rotate = "none")
model4

# rotation
model3 <- principal(hozdata, nfactors = 2,
          rotate = "varimax")
model3
model5 <-  fa(hozdata, nfactors = 2,fm = "pa", rotate = "varimax")
model5
model6 <-  fa(hozdata, nfactors = 2,fm = "pa",rotate="promax")
model6
print.psych(model6, cut = 0.3 , sort = TRUE)

###### diagram
fa.diagram(model6)

##  Obtaining factor scores for individuals 
# for the PCA analysis
model2a <- principal(hozdata, nfactors = 2, rotate = "none", scores = TRUE)
model2a

model2a$scores
head(model2a$scores, 10)
zscore_hozdata <- scale(hozdata)
head(zscore_hozdata, 10)

model2a$weights

zscore_hozdata %*% model2a$weights

plot(model2a$scores)
describe(model2a$scores)

install.packages("scatterplot3d")
library(scatterplot3d)
scatterplot3d(model2a$scores)
############
##  Obtaining factor scores for individuals 
# for the FA analysis simple in the scores attribute, two columns,as we requested two factors:
model6$scores
str(model6$scores)
# num [1:73, 1:2] -0.7 0.178 2.083 -0.454 -0.266 ...
# - attr(*, "dimnames")=List of 2
#  ..$ : NULL
# ..$ : chr [1:2] "PA1" "PA2"
# 
 attr(model6$scores,  "dimnames")
# [[1]]
# NULL
#
# [[2]]
# [1] "PA1" "PA2"
### an see each column seperately:
 model6$scores[,1] # = PA1
 model6$scores[,2] # = PA2
##########################################
##########################################
#####################
# using a matrix in place of raw data:
library(psych)
hozdatamatrix <-  matrix(c(
1.000, .696, .743, .369, .184, .230,
.696, 1.000, .724, .335, .179, .367, 
.743, .724, 1.000, .326, .211, .343, 
.369, .335, .326, 1.000, .492, .492, 
.184, .179,.211, .492, 1.000,  .483, 
.230, .367, .343, .492, .483, 1.000), ncol = 6, byrow = TRUE)
# now give the columns and rows names
colnames(hozdatamatrix)<- c("wordmean",	"sentence",	"paragrap",	"lozenges",	"cubes", "visperc")
rownames(hozdatamatrix)<- c("wordmean",	"sentence",	"paragrap",	"lozenges",	"cubes", "visperc")
######### you can NOT use the standard cor function i.e. cor(hozdatamatrix)
#  With a correlation matrix as it produces correlations of the correlations.
# you can produce a correlation plot by using the cor.plot function
# the darker the shading for the cell the higher the correlation
cor.plot(hozdatamatrix)
# According to the psych package manual you should be able to use the 
# coor.p() function below to obtain p values for the associated correlation matrix
# but all you get is the original correlations back not the p values,
# a nice email from  the author professor Revelle (package developer)
# said this bug will be fixed in the next version (1.4.3) and sent a bug fix below 
##  given below I have called it corr.p_temp
corr.p_temp <-function(r, n, adjust="holm") {
cl <- match.call()
if(missing(n)) stop("The number of subjects must be specified")
t <- (r*sqrt(n-2))/sqrt(1-r^2)
p <- 2*(1 - pt(abs(t),(n-2)))
p[p>1] <- 1
if (adjust !="none") {
if(isSymmetric(p)) {sym <- TRUE
lp <- upper.tri(p)   #the case of a symmetric matrix
   pa <- p[lp]
   pa <- p.adjust(pa,adjust)
   p[upper.tri(p,diag=FALSE)] <- pa
 } else {
p[] <- p.adjust(p ,adjust)   #the case of an asymmetric matrix
 sym <- FALSE}
} 
result <- list(r = r,n=n,t=t,p=p,sym=sym, Call=cl)
class(result) <- c("psych", "corr.test")
return(result)
}

aresult1<- corr.p_temp(r = hozdatamatrix, n= 73, adjust="holm");  aresult1
aresult2<-corr.p_temp (r = hozdatamatrix, n= 73);                 aresult2
aresult3<- corr.p_temp (r = hozdatamatrix, n= 73, adjust="none"); aresult3
aresult4<- corr.p_temp (r = hozdatamatrix, n= 73, adjust="holm"); aresult4
det(hozdatamatrix) #  gives same answer as using the raw data 0.07374609
eigen(hozdatamatrix) provides the eigenvalues 
# To carry out a PCA analysis using a correlation matrix need to 
#  tell the principal function how many observations formed the
#  basis of the correlations specifying a value for the  the n.obs parameter
model1<- principal(hozdatamatrix, nfactors = 6, n.obs = 73,  rotate = "none")
model1
# to carry out a factor analysis using a correlation matrix
# adapt the fa function in a similar way:
modelb <- fa(hozdatamatrix, nfactors = 2, fm = "pa", n.obs = 73, rotate = "none")
modelb
##############################
##### example of how the eigenvalues and iegenvector are calculated:
result1 <- eigen(hozdatamatrix) # provides the eigenvalues/vectors 
result1
# show the structure:
 str(result1)
# List of 2
#  $ values : num [1:6] 3.099 1.348 0.549 0.488 0.282 ...
#  $ vectors: num [1:6, 1:6] -0.451 -0.46 -0.467 -0.381 -0.299 ...
result1$values # gives eigenvalues
result1$vectors # gives eigenvectors
# The above value correspond to those produced by principal() but the 
vectors do not equal loading in the output
# this is because:
# The eigen vectors are rescaled by the sqrt of the 
# eigen values to produce the component loadings more typical in factor analysis.
#########
# you can see the code for the principal() by typing into the r console window
#  psych::principal
# has: loadings <- eigens$vectors %*% sqrt(diag(eigens$values))
# converting to our result1:
loadings_sd <- result1$vectors %*% sqrt(diag(result1$values))
loadings_sd
##################################################
#############################################
 ######################
 ## Allocating subject to a particular factor, for graphical purposes:
### test example first 
 # test idea of creating a new variable which indicates which of pA1 or pA2 is larger for a subject:
 a <- c(5, 7, 6, 9, 4,5,3,10,67,-5,-7,0,5, 2)
 b <- c(5, 5, 9, 9, 6,4,9,4,  1, 7 ,4,4,8,-5)
  c <- c(1, 1,1, 1, 1,1,1,1,  1, 1, 1,1,1,-1)
 tempdf = data.frame(a,b,c)
 tempdf$newvar <- ifelse(tempdf$a >  tempdf$b, 1,   0 )  # if column a value> col. b put 1 else put 0
head(tempdf) 
   a b c newvar
# 1 5 5 1      0
# 2 7 5 1      1
# 3 6 9 1      0
# 4 9 9 1      0
# 5 4 6 1      0
# 6 5 4 1      1
#  works
#########################
# apply above technique to the results:
newvar <- ifelse( model6$scores[,1] >   model6$scores[,2], 1,   2 )  
 # if PA1 value> PAb put 1 else put 2
newvar

############################
### now we have a classifier we can make use of it:
# several ways:
# 1. using R car package:
## need to create a new datafram with the necessary columns:
### biplot uses the x$scores columns for the x and y axises 
# but we want two ellipses based upon the two subsets of data defined by the category
tempdfbiplot <-  data.frame(model6$scores[,1], model6$scores[,2], newvar)
# using the R core function plot:
plot(tempdfbiplot[,1],tempdfbiplot[,2],xlim=c(-3,3), ylim=c(-3,3))
library(car) # needed for dataEllipse function in the package:
dataEllipse(tempdfbiplot[,1],tempdfbiplot[,2], 
groups= factor(tempdfbiplot[,3]) , 
add=TRUE,  levels = c(0.95), fill.alpha=0.1)
# 2. using psych package:
biplot(model6) # points all one colour
biplot(model6, group= newvar) # colours the points on the plot
# the psych package has a ellipses function, but you dont seem to be able to define subgroups of variables,
# workaround seems to be selecting each subgroup of rows in turn:  
# n specifies sigmal level to show 1= 1 sigma; 2 shows both 1 and 2 levels 
ellipses(tempdfbiplot[(tempdfbiplot$newvar==1),1:2], col="blue", add=TRUE, data=FALSE, n=2, smooth=FALSE)
ellipses(tempdfbiplot[(tempdfbiplot$newvar==2),1:2], col="red", add=TRUE, data=FALSE, n=2, smooth=FALSE)

########################
We can also do the same thing with the ggbiplot:
library(ggbiplot)
g <- ggbiplot(model1.pca, obs.scale = 1, var.scale = 1, groups=newvar, ellipse = TRUE,
 ellipse.prob = .95) # not how the ellipses get bigger as you increase the CI interval)
print(g)

#############################################################  end


DF <- matrix(sample(1:9,9),ncol=3,nrow=3)
DF <- as.data.frame.matrix(DF)
DF