setwd("C:/Dropbox/Lavoro/Didattica/Course Statistics for Economics and Business/Datasets_Codes")

#####################
### Example Passito
#####################

# import dataset
X=read.csv("Passito.csv",header=TRUE,sep=";")
head(X)
# remove unuseful variables (columns)
X=X[,-c(1,2,3,4,5)]
head(X)
dim(X)

# perform factor analysis with no rotation and 5 factors
F=factanal(X,factors=5,rotation="none",scores="regression")
# compute correlations between the original variables 
R=F$correlation
R
# compute eigenvalues and eigenfactors of the matrix of correlations
e=eigen(R)
e$values
# create the scree plot to represent the eigenvalues
barplot(rep(1,ncol(X))),col="yellow",ylim=c(0,max(e$values)+2))
barplot(e$values,main="Scree plot",ylab="eigenvalue",add=TRUE)
F

# perform factor analysis with varimax rotation and 5 factors
F=factanal(X,factors=5,rotation="varimax",scores="regression")
F

# perform factor analysis with promax rotation and 5 factors
F=factanal(X,factors=5,rotation="promax",scores="regression")
F
# Factor loadings
L=F$loadings
L
# variances of unique factors (proportions of variance of the original observed variables not explained by the factors)
U=F$uniquenesses
U
# communalities (proportions of variance of the original observed variables explained by the factors)
Comm=1-U
Comm
# rotation matrix
G=F$rotmat
G
# factor scores (factor values for each statistical unit (country) )
FScores=F$scores
FScores
# p-value of the test on the hypothesis that the considered factrs are sufficient
p=F$PVAL
p

# Principal Component Analysis
PC=prcomp(X,retx =TRUE,center=TRUE,scale.=TRUE)
e_values=PC$sdev^2
e_values
e_vectors=PC$rotation
e_vectors

# YScores=PC$x
YScores=as.matrix(X)%*%e_vectors

ifelse(abs(e_vectors)>0.3,e_vectors,0) 

#####################
### Example Mall
#####################


X=read.csv("Mall.csv",header=TRUE,sep=";")
head(X)
dim(X)

F=factanal(X,factors=2,rotation="none",scores="regression")
R=F$correlation
R
e=eigen(R)
e$values
barplot(rep(1,ncol(X)),col="yellow",ylim=c(0,max(e$values)+2))
barplot(e$values,main="Scree plot",ylab="eigenvalue",add=TRUE)
F

F=factanal(X,factors=2,rotation="varimax",scores="regression")
F

F=factanal(X,factors=2,rotation="promax",scores="regression")
F
L=F$loadings
L
U=F$uniquenesses
U
Comm=1-U
Comm
R=F$correlation
R
G=F$rotmat
G
FScores=F$scores
FScores
plot(FScores[,1],FScores[,2])


#####################
### Example Eating Habits
#####################

data=read.csv("EatingHabits.csv",header=TRUE,sep=";")
head(data)
X=data[,-c(1,2,15)]
dim(X)

F=factanal(X,factors=5,rotation="none",scores="regression")
R=F$correlation
R
plot(X)
e=eigen(R)
e$values
barplot(rep(1,ncol(X)),col="yellow",ylim=c(0,max(e$values)+2))
barplot(e$values,main="Scree plot",ylab="eigenvalue",add=TRUE)
F

F=factanal(X,factors=4,rotation="varimax",scores="regression")
F

F=factanal(X,factors=4,rotation="promax",scores="regression")
F
L=F$loadings
L
ifelse(L>0.5,L,0)
U=F$uniquenesses
U
Comm=1-U
Comm
R=F$correlation
R
G=F$rotmat
G
FScores=F$scores
FScores
class(FScores)
cor(FScores)
plot(as.data.frame(FScores))

data$Countries
FScores[data$Countries=="Italy",]

p=F$PVAL
p


PC=prcomp(X,retx =TRUE,center=TRUE,scale.=TRUE)
e_values=PC$sdev^2
e_values
e_vectors=PC$rotation
e_vectors
YScores=PC$x
YScores