Code : Tout sélectionner
library(rgl)
library(MASS)
library(ade4)
co12 <- 0.86 ; co13 <- 0.87 ; co23 <- 0.99
n <- 47 # nombre de points
Sigma <- matrix(c(1.0, co12, co13, co12, 1.0, co23, co12, co23, 1.0), 3, 3)
set.seed(1) # pour reproductibilité
dta <- MASS::mvrnorm(n = n, rep(0, 3), Sigma, empirical = TRUE)
acp <- ade4::dudi.pca(dta, scannf = FALSE, nf = 3)
# Sanity check : les axes de l'ACP sont bien tous orthogonaux :
X <- as.matrix(acp$c1) # axes de l'ACP de norme unité
zapsmall(X %*% t(X))
# V1 V2 V3
# V1 1 0 0
# V2 0 1 0
# V3 0 0 1
#
# Ajout des axes de la base cannonique
#
mkaxes <- function(dta){
noms <- paste("Axe", 1:3)
for(i in 1:3){
p1 <- rep(0, 3) ; p1[i] <- max(dta[ , i])
arrow3d(p0 = c(0, 0, 0), p1 = p1, s = 1/20)
text3d(1.2*p1[1], 1.2*p1[2], 1.2*p1[3], noms[i])
}
}
#
# Ajout des axes de l'ACP
#
mkaxesell <- function(dta, scalefactor = 6){
nomaxex <- c("Longueur", "Largeur", "Épaisseur")
for(i in 1:3){
p1 <- scalefactor*dta[ , i] ; p0 <- -1*p1
arrow3d(p0 = p0, p1 = p1, s = 1/20, col = "red")
text3d(1.2*p1[1], 1.2*p1[2], 1.2*p1[3], nomaxex[i], col = "red")
}
}
scalefactor <- 2
mycol <- grey(0.7)
plot3d(dta, type = "p", box = FALSE, axes = FALSE, xlab = "", ylab = "", zlab = "")
mkaxes(dta)
mkaxesell(acp$c1, scalefactor)
plot3d(ellipse3d(cor(dta)), col = grey(0.9), add = TRUE, alpha = 0.5)
for(i in -scalefactor:scalefactor)
for(j in -scalefactor:scalefactor)
for(z in -scalefactor:scalefactor){
lines3d(c(i, i), c(j, j), c(-z, z), col = mycol)
lines3d(c(-i, i), c(j, j), c(z, z), col = mycol)
lines3d(c(i, i), c(-j, j), c(z, z), col = mycol)
}
sessionInfo()
# R version 4.0.3 (2020-10-10)
# Platform: x86_64-apple-darwin17.0 (64-bit)
# Running under: macOS Catalina 10.15.7
#
# Matrix products: default
# BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
#
# locale:
# [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
#
# attached base packages:
# [1] stats graphics grDevices utils datasets methods base
#
# other attached packages:
# [1] ade4_1.7-15 MASS_7.3-53 rgl_0.103.5
#
# loaded via a namespace (and not attached):
# [1] Rcpp_1.0.4.6 knitr_1.30 magrittr_1.5 xtable_1.8-4 R6_2.4.1
# [6] rlang_0.4.6 fastmap_1.0.1 tools_4.0.3 webshot_0.5.2 xfun_0.19
# [11] miniUI_0.1.1.1 htmltools_0.4.0 crosstalk_1.1.0.1 yaml_2.2.1 digest_0.6.25
# [16] shiny_1.4.0.2 later_1.0.0 htmlwidgets_1.5.1 promises_1.1.0 manipulateWidget_0.10.1
# [21] evaluate_0.14 mime_0.9 rmarkdown_2.5 compiler_4.0.3 jsonlite_1.6.1
# [26] httpuv_1.5.2
Faites tourner le graphique pour mettre l'axe 1 horizontal orienté vers la droite, l'axe 2 vertical orienté vers le haut et l'axe 3 perpendiculaire au plan, pas de problème : les axes 1 et 2 sont bien à angle droit. Faites tourner le graphique pour mettre l'axe 1 horizontal orienté vers la droite, l'axe 3 vertical orienté vers le haut et l'axe 2 perpendiculaire au plan, pas de problème : les axes 1 et 3 sont bien à angle droit. Faites tourner le graphique pour mettre l'axe 2 horizontal orienté vers la droite, l'axe 3 vertical orienté vers le haut et l'axe 1 perpendiculaire au plan, pas de problème : les axes 2 et 3 sont bien à angle droit. De cette petite expérience je déduis que rgl est parfaitement capable de représenter correctement des axes orthogonaux.
Essayez maintenant de faire l'ACP à la main en mettant l'axe de la longueur horizontal orienté vers la droite, l'axe de la largeur vertical orienté vers le haut et l'axe de l'épaisseur perpendiculaire au plan. Problème : l'axe de la longueur et de la largeur ne sont pas à angle droit, l'axe de la largeur penche vers la gauche.
Qu'est-ce qui peut bien clocher ici ?
Amicalement,
jean lobry
