#####################################
# version 0.3.3 du 19 novembre 2006 #
#####################################

Bartlett.test.sphericity <- function(X){

##########################################################
# Dziuban, C. D., & Shirkey, E. C. (1974). When is       #
# a correlation matrix appropriate for factor analysis ? #
# Psychological Bulletin, 81(6), 358-361.                #
##########################################################

# X est une matrice dont on aimerait faire l'analyse en 
# composantes principales

n <- dim(X)[[1]]
p <- dim(X)[[2]]
R <- cor(X)
khi2 <- -(n - 1 - (2*p + 5)/6)*log(det(R))
ddl <- p*(p - 1)/2
valeur.p <- pchisq(khi2, ddl, lower.tail=FALSE)
return(list(khi2=khi2, ddl=ddl, valeur.p=valeur.p))
}

KMO.test <- function(X){

##########################################################
#    Kaiser-Meyer-Olkin measure of sampling adequacy     #
#                                                        #
# Dziuban, C. D., & Shirkey, E. C. (1974). When is       #
# a correlation matrix appropriate for factor analysis ? #
# Psychological Bulletin, 81(6), 358-361.                #
##########################################################

R <- cor(X)
S <- diag(sqrt(1/diag(solve(R))))
A <- S %*% solve(R) %*% S
r2 <- sum(R^2) - sum(diag(R^2))
a2 <- sum(A^2) - sum(diag(A^2))
KMO <- r2/(r2 + a2)
conclusion <- NULL
if(0.9 <= KMO){conclusion <- "marvelous"}
if(0.8 <= KMO & KMO < 0.9){conclusion <- "meritorious"}
if(0.7 <= KMO & KMO < 0.8){conclusion <- "middling"}
if(0.6 <= KMO & KMO < 0.7){conclusion <- "mediocre"}
if(0.5 <= KMO & KMO < 0.6){conclusion <- "miserable"}
if(KMO < 0.5){conclusion <- "unacceptable"}
print(conclusion)
return(KMO)
}

# acp est le resultat d'une analyse en composantes principales

#############################
# Coordonnees des individus #
#############################

coor.ind <- function(acp, format=3){
return(round(acp$scores, format))
}
#############################
# Coordonnees des variables #
#############################

coor.var <- function(acp, format=3){
sat <- acp$loadings[] %*% diag(acp$sdev)
return(round(sat, format))
}

###########################
# Cercle des correlations #
###########################

proj.var <- function(acp, dim1=1, dim2=2, pos=3){
r <- coor.var(acp)
p <- dim(r)[1]
z <- seq(0, 2*pi, length=1000)
x <- cos(z)
y <- sin(z)

a <- attributes(r)
par(pty="s")
plot(r[, dim1], r[,dim2], type="n", xlim=c(-1, 1), ylim=c(-1, 1), xlab=paste("axe", dim1), ylab=paste("axe", dim2),
main="Cercle des correlations")

abline(v=0)
abline(h=0)
polygon(x,y)
arrows(rep(0, p), rep(0, p), r[, dim1], r[, dim2], length=0.15, angle=20)
text(r[, dim1],r[, dim2],a$dimnames[[1]],pos=pos)
}

################################
# Representation des individus #
################################

proj.ind <- function(acp, dim1=1, dim2=2, pos=3, xlim=NULL, ylim=NULL){
X <- coor.ind(acp)
x <- X[, dim1]
y <- X[, dim2]
etiquettes <- names(x)
plot(x, y, type="n", xlim=xlim, ylim=ylim, xlab=paste("Composante", dim1), ylab=paste("Composante", dim2), main="Projection des individus")
grid()
abline(v=0)
abline(h=0)
points(x, y)
text(x, y, etiquettes, pos=pos, cex=0.8)
}

###########################################
# Contribution des individus aux facteurs #
###########################################

CTR.ind <- function(acp, nb.comp=2, fraction=10000){
X2 <- acp$scores^2
contr <- fraction*(X2 %*% diag(1/apply(X2, 2, sum)))
dimnames(contr)[[2]] <- paste("Comp.", 1:(dim(X2)[2]), sep="")
return(round(contr[, 1:nb.comp]))
}

###########################################
# Contribution des variables aux facteurs #
###########################################

CTR.var <- function(acp, nb.comp=2, fraction=10000){
X2 <- coor.var(acp)^2
contr <- fraction*(X2 %*% diag(1/apply(X2, 2, sum)))
dimnames(contr)[[2]] <- paste("Comp.", 1:(dim(X2)[2]), sep="")
return(round(contr[, 1:nb.comp]))
}

###########################################
# Qualite de representation des individus #
###########################################

QLT.ind <- function(acp, nb.comp=2, fraction=10000){
X2 <- acp$scores^2
D2 <- apply(X2, 1, sum)
qual <- diag(1/D2) %*% X2
dimnames(qual) <- list(dimnames(X2)[[1]], paste("Comp.", 1:(dim(X2)[2]), sep=""))
return(round(10000*cbind(qual[, 1:nb.comp], somme=apply(qual[, 1:nb.comp], 1, sum))))
}

###########################################
# Qualite de representation des variables #
###########################################

QLT.var <- function(acp, nb.comp=2, fraction=10000){
X2 <- (acp$loadings[] %*% diag(acp$sdev))^2
D2 <- apply(X2, 1, sum)
qual <- diag(1/D2) %*% X2
dimnames(qual) <- list(dimnames(X2)[[1]], paste("Comp.", 1:(dim(X2)[2]), sep=""))
return(round(10000*cbind(qual[, 1:nb.comp], somme=apply(qual[, 1:nb.comp], 1, sum))))
}

#########################
# Inertie des individus #
#########################

Inertie.ind <- function(acp, format.distance=3, fraction.inertie=10000){
X2 <- acp$scores^2
D2 <- apply(X2, 1, sum)
return(cbind(distance=round(sqrt(D2), format.distance), inertie=round(fraction.inertie*D2/sum(D2))))
}

#############################
# Individus supplementaires #
#############################

suppl.ind <- function(acp, adj, format=3){
if(is.vector(adj)){print("Les individus supplementaires sont les lignes d'une MATRICE."); stop}
n <- dim(adj)[1]
coor <- (as.matrix(adj) - (rep(1, n) %*% t(acp$center))) %*% diag(1/acp$scale) %*% acp$loadings
return(round(coor, format))
}

#############################
# Variables supplementaires #
#############################

suppl.var <- function(acp, adj, format=3){
if(is.vector(adj)){n <- length(adj)}
if(is.matrix(adj)){n <- dim(adj)[1]}
adj.cr <- centrer.reduire(adj)
score.cr <- centrer.reduire(acp$score)
coor <- (1/n)*t(adj.cr) %*% score.cr
return(round(coor, format))
}

############################################
# Coordonnees des individus apres rotation #
############################################

coor.ind.rot <- function(acp, nb.comp, format=3){
cor.v <- acp$loadings[] %*% diag(acp$sdev)
rot <- varimax(cor.v[, 1:nb.comp])
coor <- acp$scores[, 1:nb.comp] %*% rot$rotmat
return(round(coor, format))
}

############################################
# Coordonnees des variables apres rotation #
############################################

coor.var.rot <- function(acp, nb.comp, format=3){
cor.v <- acp$loadings[] %*% diag(acp$sdev)
rot <- varimax(cor.v[, 1:nb.comp])
return(round(rot$loadings[], format))
}

##########################################
# Cercle des correlations apres rotation #
##########################################

proj.var.rot <- function(acp, nb.comp, dim1=1, dim2=2, pos=3){
# nb.comp est le nombre de composantes retenues
# dim1 et dim2 definissent le plan de projection
r <- coor.var(acp, format=20)
rot <- varimax(r[, 1:nb.comp])$loadings
p <- dim(rot)[1]

a <- attributes(r)
par(pty="s")
plot(rot[, dim1], rot[,dim2], type="n", xlim=c(-1, 1), ylim=c(-1, 1), xlab=paste("axe", dim1), ylab=paste("axe", dim2),
main="Projection des variables apres rotation")

abline(v=0)
abline(h=0)
symbols(0, 0, circles=1, inches=FALSE, add=TRUE)
arrows(rep(0, p), rep(0, p), rot[, dim1], rot[, dim2], length=0.15, angle=20)
text(rot[, dim1],rot[, dim2],a$dimnames[[1]],pos=pos)
}

###############################################
# Representation des individus apres rotation #
###############################################

proj.ind.rot <- function(acp, nb.comp, dim1=1, dim2=2, pos=3, xlim=NULL, ylim=NULL){
X <- coor.ind.rot(acp, nb.comp)
x <- X[, dim1]
y <- X[, dim2]
etiquettes <- names(x)
plot(x, y, type="n", xlim=xlim, ylim=ylim, xlab=paste("axe", dim1), ylab=paste("axe", dim2), main="Projection des individus apres rotation")
grid()
abline(v=0)
abline(h=0)
points(x, y)
text(x, y, etiquettes, pos=pos, cex=0.8)
}

##########################################################
# Contribution des individus aux facteurs apres rotation #
##########################################################

CTR.ind.rot <- function(acp, nb.comp, fraction=10000){
X2 <- coor.ind.rot(acp, nb.comp, format=20)^2
contr <- fraction*(X2 %*% diag(1/apply(X2, 2, sum)))
dimnames(contr)[[2]] <- paste("Axe", 1:(dim(X2)[2]), sep="")
return(round(contr[, 1:nb.comp]))
}

##########################################################
# Contribution des variables aux facteurs apres rotation #
##########################################################

CTR.var.rot <- function(acp, nb.comp, fraction=10000){
X2 <- coor.var.rot(acp, nb.comp, format=20)^2
contr <- fraction*(X2 %*% diag(1/apply(X2, 2, sum)))
dimnames(contr)[[2]] <- paste("Axe", 1:(dim(X2)[2]), sep="")
return(round(contr[, 1:nb.comp]))
}

##########################################################
# Qualite de representation des individus apres rotation #
##########################################################

QLT.ind.rot <- function(acp, nb.comp, fraction=10000){
X2 <- acp$scores^2
D2 <- apply(X2, 1, sum)
Y2 <- coor.ind.rot(acp, nb.comp, format=20)^2
qual <- diag(1/D2) %*% Y2
dimnames(qual) <- list(dimnames(Y2)[[1]], paste("Axe", 1:(dim(Y2)[2]), sep=""))
return(round(10000*cbind(qual[, 1:nb.comp], somme=apply(qual[, 1:nb.comp], 1, sum))))
}

##########################################################
# Qualite de representation des variables apres rotation #
##########################################################

QLT.var.rot <- function(acp, nb.comp=2, fraction=10000){
X2 <- (acp$loadings[] %*% diag(acp$sdev))^2
D2 <- apply(X2, 1, sum)
Y2 <- coor.var.rot(acp, nb.comp, format=20)^2
qual <- diag(1/D2) %*% Y2
dimnames(qual) <- list(dimnames(Y2)[[1]], paste("Axe", 1:(dim(Y2)[2]), sep=""))
return(round(10000*cbind(qual[, 1:nb.comp], somme=apply(qual[, 1:nb.comp], 1, sum))))
}

############################################
# Individus supplementaires apres rotation #
############################################

suppl.ind.rot <- function(acp, adj, nb.comp=2, format=3){
if(is.vector(adj)){print("Les individus supplementaires sont les lignes d'une MATRICE."); stop}
n <- dim(adj)[1]
coor <- (as.matrix(adj) - (rep(1, n) %*% t(acp$center))) %*% diag(1/acp$scale) %*% acp$loadings
cor.v <- acp$loadings[] %*% diag(acp$sdev)
rot <- varimax(cor.v[, 1:nb.comp])
coor <- coor[, 1:nb.comp] %*% rot$rotmat
return(round(coor, format))
}

############################################
# Variables supplementaires apres rotation #
############################################

suppl.var.rot <- function(acp, adj, nb.comp=2, format=3){
if(is.vector(adj)){n <- length(adj)} else {n <- dim(adj)[1]}
adj.cr <- centrer.reduire(as.matrix(adj))
score.cr <- centrer.reduire(acp$score)
coor <- (1/n)*t(adj.cr) %*% score.cr
cor.v <- acp$loadings[] %*% diag(acp$sdev)
rot <- varimax(cor.v[, 1:nb.comp])
coor <- coor[, 1:nb.comp] %*% rot$rotmat
return(round(coor, format))
}

####################################################
# Adjonction d'individus dans le plan (dim1, dim2) #
####################################################

adj.ind <- function(coor, etiquettes=NULL, dim1=1, dim2=2, col="red", pch=19, cex=1, pos=3){

# coor est le vecteur (la matrice) des coordonnees de l'individu (des individus)
# que l'on souhaite positionner dans l'espace factorielle

if(is.vector(coor)){points(coor[dim1], coor[dim2], col=col, pch=pch, cex=cex)
text(coor[dim1], coor[dim2], etiquettes, col=col, pos=pos)
} else{
points(coor[, dim1], coor[, dim2], col=col, pch=pch, cex=cex)
text(coor[, dim1], coor[, dim2], etiquettes, col=col, pos=pos)}
}

######################################################
# Adjonction de variables dans l'espace (dim1, dim2) #
######################################################

adj.var <- function(coor, etiquettes=NULL, dim1=1, dim2=2, col="red", pos=3){

# coor est le vecteur (la matrice) des coordonnees de la variable (des variables)
# que l'on souhaite positionner dans l'espace factorielle

if(is.vector(coor)){arrows(0, 0, coor[dim1], coor[dim2], col=col, length=0.15, angle=20)
text(coor[dim1], coor[dim2], etiquettes, col=col, pos=pos)
} else{
arrows(0, 0, coor[, dim1], coor[, dim2], col=col, length=0.15, angle=20)
text(coor[, dim1], coor[, dim2], etiquettes, col=col, pos=pos)}
}

##############################################
# Determination du nombre de valeurs propres #
##############################################

valeurs.propres.significatives <- function(dim1, dim2, p, n.iter=500){
lambda <- NULL
for(i in 1:n.iter){
acp <- princomp(matrix(rnorm(dim1 * dim2), ncol= dim2), cor=TRUE, scores=FALSE)
lambda <- rbind(lambda, acp$sdev^2)
}
q <- round(apply(lambda, 2, quantile, p),2)
names(q) <- paste("comp", 1:dim2)
return(q)
}

######################
# Centrer et reduire #
######################

centrer.reduire <- function(data){
if(is.data.frame(data)){X <- as.matrix(data)}else{X <- data}
if(is.vector(X)){n <- length(X)
X.cr <- (X - mean(X))/( sqrt((n - 1)/n) *sd(X))}
if(is.matrix(X)){
n <- dim(X)[1]
p <- dim(X)[2]
if(p == 1){X.cr <- (X - mean(X))/( sqrt((n - 1)/n) *sd(X))} else{
m <- apply(X, 2, mean)
X.c <- X - (rep(1, n) %*% t(m))
s <- sqrt(apply(X.c^2, 2, mean))
X.cr <- X.c %*% diag(1/s)
dimnames(X.cr) <- dimnames(X)}}
return(X.cr)
}