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 :

Sources