Analyse en composantes principales (ACP) en R
L'analyse en composantes principales est un classique des statistiques. Très populaire, elle est souvent utilisée en R via le packages ADE4.
En R :
Cette fonction a été implémentée au sein d’un cours d’analyse des données en M1 Data Science et n’a pas été vérifiée depuis. Une erreur n’est donc pas impossible. J’ai constaté des outputs identiques à ceux de l’ACP d’ADE4 pour différents datasets. Si jamais vous remarquez un problème, n’hésitez pas à m’envoyer un mail à [email protected].
ACP <- function(m, mode = 3, nbDim = 2, inertieMin=80, centerReduc=TRUE, df=FALSE, aff = TRUE){
# Dimension de la matrice
nR <- nrow(m)
nC <- ncol(m)
# Si m est un dataframe, nous voulons le transformer en matrice
if (df){
m = as.matrix(m)
}
# Si nous souhaitons centrer et réduire la matrice m
if (centerReduc){
mMeans <- colMeans(m)
mMeans <- matrix(mMeans, nrow=nR, ncol=length(mMeans), byrow=TRUE)
m <- m - mMeans
# L'écart-type utilise une division par (n-1) et non par (n) sur R.
# C'est la raison pour laquel nous la calculons manuellement sans utiliser "sd(x)"
SD <- sqrt(colMeans((m - colMeans(m))^2))
m <- sweep(m,2, SD,'/')
sq <- 1/nR * t(m) %*% m
}
else {
sq <- cov(m)
}
# Calcul des valeurs et vecteurs propres (norme euclydienne = 1)
eig = eigen(sq, symmetric = TRUE)
eigVal = eig$val # Valeurs propres
eigVect = eig$vectors # Vecteurs propre
# On reverse la somme cumulé des valeurs propres
cumEigValLoss = rev(cumsum(rev(eigVal)))
cumEigVal <- cumsum(eigVal) # Valeurs propres cumulés
inertieTotal <- sum(eigVal) # Inertie total du nuage de point
nbVP <- length(eigVal) # Nombre de valeurs propres
inertieMoyenne <- inertieTotal/nbVP # Inertie moyenne ( = 1 si centré)
cumEigValPercent <- cumEigVal/inertieTotal*100 # Inertie cumulé en %.
if (mode == 1){ # Règle de Kaiser
# On séléctionne les vecteurs/valeurs propres uniquement > Inertie moyenne
selEigVal <- eigVal[eigVal >= inertieMoyenne]
selEigVect <- eigVect[,1:length(selEigVal)]
}
else if (mode == 2){ # Nombre de dimension
# On prend littéralement nbDim dimensions
selEigVect <- eigVect[,1:nbDim]
selEigVal <- eigVal[1:nbDim]
}
else { # Règle de l'inertie minimum
nbSelVal <- length(cumEigValPercent[cumEigValPercent < inertieMin])+1
selEigVal <- eigVal[1:nbSelVal]
selEigVect <- eigVect[,1:nbSelVal]
}
nbSelVal <- length(selEigVal) # Nombre de valeurs propres séléctionnées
sumSelVal <- sum(selEigVal) # Inertie gardé
loss = 100 - sumSelVal/inertieTotal*100 # Perte d'inertie
###########################
# INDIVIDUS
# Composantes principales : Coordonnées des individus
cooInd = m %*% selEigVect
# Contribution des individus pour chaque composante
contribInd <- 100/nR * sweep(cooInd^2,2,selEigVal, '/')
# Qualité de représentation des individus pour chaque composante
qualInd <- 100 * sweep(cooInd^2,1,t(rowSums(m^2)), '/') * sign(cooInd)
###########################
# VARIABLES
# Coordonnées des variables
cooVar = (selEigVect %*% diag(sapply(selEigVal, sqrt)))
# Contribution des variables pour chaque composante
contribVar <- 100 * sweep(cooVar^2,2,selEigVal, '/')
# Qualité de représentation des varariables pour chaque composante
qualVar <- 10^3 * sweep(cooVar^2,1,t(colSums(m^2)), '/') * sign(cooVar)
###########################
)
}
newcoo <-ACP(A,df = TRUE, mode=2, reduc=TRUE, aff=TRUE, center=TRUE)
Alternative
Alternativement, vous pouvez retrouver le code utilisé par ADE4 sur GitHub dans les fichiers :