TP2 ACP

Fonction acp_from_scratch calculs matriciellement

acp_from_scratch <- function(X, nombre_composants=NULL, reduire = FALSE) {
  
  n <- nrow(X)
  X <- as.matrix(X)
  
  # Data centree
  data <- centrage(X)
  
  # Data normee
  if (reduire) {
    data <- centrer_reduire(X)
  }

  #calculer la matrice de covariance matriciellement
  matrice_covariance <- (1/n) * t(data) %*% data
  
  # Calcul des valeurs propres et vecteurs propres
  eig_result <- eigen(matrice_covariance)
  valeurs_propres <- eig_result$values
  vecteurs_propres <- eig_result$vectors
  
  # Sélectionner le nombre spécifié de composants principaux
  if (!is.null(nombre_composants)) {
    vecteurs_propres <- vecteurs_propres[, 1:nombre_composants, drop = FALSE]
  }
  
  # Projection des données sur les composants principaux
  nouvelles_coordonnees <- data %*% vecteurs_propres
  
  return(list(matrice_covariance=matrice_covariance,
              nouvelles_coordonnees=nouvelles_coordonnees, 
              valeurs_propres=valeurs_propres, 
              vecteurs_propres=vecteurs_propres))
  
}

Fonction qui centre les données

centrage <- function(data) {
  n <- nrow(data)
  
  # Centrage des données
  # Définir la matrice des poids
  matrice_de_poids <- (1/n) * diag(n)
  
  # Calculer la moyenne pondérée matriciellement
  mat1_n <- matrix(rep(1, n), nrow=1)
  moyenne_ponderee <- mat1_n %*% data / n
  
  #calculer data_centree
  data_centree <- sweep(data, 2, moyenne_ponderee)
  
  return (data_centree)
}

Fonction qui normalise les données

centrer_reduire <- function(donnees) {
  n <- nrow(donnees)  # Nombre d'observations
  p <- ncol(donnees)  # Nombre de variables
  
  # Centrage des données
  donnees_centrees = centrage(donnees)
  
  # Réduction des données (calcul de l'écart-type de façon matricielle avec biais)
  # Calculer la matrice de covariance avec biais (diviser par n!)
  matrice_covariance <- (1 / n) * (t(donnees_centrees) %*% donnees_centrees)
  # Extraire l'écart-type à partir de la diagonale de la matrice de covariance
  ecarts_types <- sqrt(diag(matrice_covariance))
  
  # Réduire les données : diviser par l'écart-type
  donnees_reduites <- sweep(donnees_centrees, 2, ecarts_types, FUN = "/")
  
  return(donnees_reduites)
}

1. et 2. On a le pseudo code dans le rendu du TP1.1, testé et approuvé.

Importer les données

# Données
data_PDE19 <- read.csv("data_PDE20.csv", sep=";", dec=",")
#View(data_PDE19)
data_PDE19$Num <- NULL
data_PDE19$X <- NULL
# Convertir le data.frame en matrice
data_PDE19 <- as.matrix(data_PDE19)

1. visualisation de la matrice des donnees

1.1 afficher la dimension de la matrice

cat("Matrice X des données de dimension (n,p) =",
    "(", paste(dim(data_PDE19), collapse = ", "), ")\n\n")

## Matrice X des données de dimension (n,p) = ( 26, 8 )

1.2 Afficher la matrice avec les noms des colonnes et leurs lignes

data_PDE19_df <- as.data.frame(data_PDE19)
# Ajouter une colonne avec les numéros de ligne
data_PDE19_df$Row <- seq_len(nrow(data_PDE19_df))
# Faire de cette colonne la première colonne
data_PDE19_df <- data_PDE19_df[, c(ncol(data_PDE19_df), 1:(ncol(data_PDE19_df)-1))]

#Pander
pander(data_PDE19_df)

Row	X1	X2	X3	X4	X5	X6	X7	X8
1	303.1	24.19	0	3.29	180	8.09	360.9	120.3
2	281.9	38.59	4.29	1.06	192	10.5	353.5	117
3	277.1	34.79	0	6.85	183.8	38.89	344	114.7
4	276.4	32.43	4.14	2.04	190.8	38.53	341.2	113.9
5	253.8	39.5	3.04	1	173.8	19.33	382.1	127.4
6	243.6	34.39	2.79	3.43	166.7	27.59	391.1	130.4
7	277	34.7	0	6.85	183.8	38.8	343.9	114.7
8	294.8	28.29	1.85	1.83	182.3	10.29	360.2	120
9	303	24.2	0	3.3	180	8.1	361	120.3
10	269.4	36.89	2.99	1.03	197.7	12.59	359.5	19.82
11	283.6	28	9.3	0	186.6	13.2	359.4	119.8
12	290.3	23.2	0.8	2.34	172.4	39.4	353.6	117.9
13	285.1	25.9	0.93	7.78	180.1	39	345.7	115.2
14	265.5	40.39	0.95	5.14	184.4	38.83	348.6	116.2
15	261.9	41.49	2.33	2.89	187.3	38.69	349.1	116.3
16	274.4	29.79	6.69	0	183.6	38.89	350	116.7
17	257.9	37.2	2.96	1.1	170.9	18.89	383.3	127.8
18	238.2	29.8	2.6	0.8	166.7	14.14	410.1	136.9
19	236	33.39	5.6	0.39	166.7	15.23	407	135.7
20	247.8	36.69	5.03	1.79	166.7	27.09	386.2	128.8
21	266.6	36.4	0	2.9	166.7	30.68	365.6	121.9
22	264.8	24.19	1.19	5.6	166.7	39.69	391.9	130.7
23	235.7	24.99	0.99	4.29	166.7	39.69	392	130.7
24	239.6	47.89	0.4	4.19	166.7	39.69	376.1	125.3
25	233.1	46.59	2.29	4.83	166.7	39.69	380.1	126.7
26	241.4	34.5	5.15	0.39	166.7	39.69	383.8	127.9

1.3 Calcul des indicateurs stats classiques (mean, sd, median, quartiles etc…) pour l’ensemble des données

#SUPP cette colonne
data_PDE19_df$Row <- NULL
library(pander)

data_PDE19_stats <- data_PDE19_df

# Calculer les statistiques
stats <- sapply(data_PDE19_stats[sapply(data_PDE19_stats, is.numeric)], 
                function(x) 
                  c(Mean = mean(x, na.rm = TRUE), 
                    SD = sd(x, na.rm = TRUE), 
                    Median = median(x, na.rm = TRUE),
                    Min = min(x, na.rm = TRUE),
                    Max = max(x, na.rm = TRUE),
                    Range = diff(range(x, na.rm = TRUE)),
                    Q1 = quantile(x, 0.25, na.rm = TRUE),
                    Q3 = quantile(x, 0.75, na.rm = TRUE),
                    IQR = IQR(x, na.rm = TRUE)))


# Transformer en dataframe pour pander
stats_df <- as.data.frame(t(stats))

# Utiliser pander pour afficher le tableau
pander(stats_df)

	Mean	SD	Median	Min	Max	Range	Q1.25%	Q3.75%	IQR
X1	265.4	21.51	266	233.1	303.1	70	244.6	280.7	36.06
X2	33.4	6.877	34.45	23.2	47.89	24.69	28.07	37.12	9.05
X3	2.55	2.378	2.31	0	9.3	9.3	0.8325	3.865	3.032
X4	2.889	2.247	2.615	0	7.78	7.78	1.038	4.265	3.228
X5	176.8	9.835	176.9	166.7	197.7	31.03	166.7	183.8	17.1
X6	27.89	12.68	34.61	8.09	39.69	31.6	14.41	38.97	24.56
X7	368.5	20.36	361	341.2	410.1	68.97	350.9	383.6	32.79
X8	119	21.33	120.3	19.82	136.9	117.1	116.4	127.9	11.46

Il existe une grande disparité entre certaines variables, comme X1, X7 et X8, ce qui contribue significativement à l’inertie du modèle. Par ailleurs, les médianes sont proches des moyennes, ce qui indique qu’il n’y a pas de valeurs extrêmes notables pour une variable donnée. Cette observation est confirmée par un examen direct des données.

2. Centrage et normalisation des données

données centrées

donnees_centrees <- as.data.frame(centrage(data_PDE19))


# Ajouter une colonne avec les numéros de ligne
donnees_centrees$Row <- seq_len(nrow(donnees_centrees))
# Faire de cette colonne la première colonne
donnees_centrees <- donnees_centrees[, c(ncol(donnees_centrees), 1:(ncol(donnees_centrees)-1))]


# Mettre l'option pour avoir 3 chiffres après la virgule
panderOptions('digits', 3)
pander(donnees_centrees)

Row	X1	X2	X3	X4	X5	X6	X7	X8
1	37.6	-9.21	-2.55	0.401	3.21	-19.8	-7.55	1.38
2	16.4	5.19	1.74	-1.83	15.2	-17.4	-15	-1.95
3	11.6	1.39	-2.55	3.96	6.99	11	-24.5	-4.3
4	10.9	-0.969	1.59	-0.849	14	10.6	-27.3	-5.04
5	-11.6	6.1	0.49	-1.89	-2.98	-8.56	13.7	8.42
6	-21.9	0.991	0.24	0.541	-10.1	-0.302	22.7	11.4
7	11.6	1.3	-2.55	3.96	7	10.9	-24.5	-4.3
8	29.4	-5.11	-0.7	-1.06	5.51	-17.6	-8.25	1.05
9	37.6	-9.2	-2.55	0.411	3.22	-19.8	-7.45	1.39
10	3.93	3.49	0.44	-1.86	20.9	-15.3	-8.98	-99.1
11	18.2	-5.4	6.75	-2.89	9.82	-14.7	-9.03	0.853
12	24.9	-10.2	-1.75	-0.549	-4.38	11.5	-14.8	-1.08
13	19.6	-7.5	-1.62	4.89	3.32	11.1	-22.8	-3.72
14	0.0327	6.99	-1.6	2.25	7.61	10.9	-19.9	-2.76
15	-3.58	8.09	-0.22	0.00115	10.5	10.8	-19.4	-2.61
16	8.93	-3.61	4.14	-2.89	6.8	11	-18.5	-2.3
17	-7.55	3.8	0.41	-1.79	-5.88	-9	14.8	8.82
18	-27.2	-3.6	0.0496	-2.09	-10.1	-13.8	41.7	18
19	-29.5	-0.00923	3.05	-2.5	-10.1	-12.7	38.5	16.7
20	-17.7	3.29	2.48	-1.1	-10.1	-0.802	17.8	9.8
21	1.12	3	-2.55	0.0112	-10.1	2.79	-2.85	2.93
22	-0.657	-9.21	-1.36	2.71	-10.1	11.8	23.4	11.7
23	-29.8	-8.41	-1.56	1.4	-10.1	11.8	23.5	11.7
24	-25.9	14.5	-2.15	1.3	-10.1	11.8	7.62	6.4
25	-32.4	13.2	-0.26	1.94	-10.1	11.8	11.6	7.74
26	-24.1	1.1	2.6	-2.5	-10.1	11.8	15.3	8.97

données normés

donnees_normee <- as.data.frame(centrer_reduire(data_PDE19))

# Ajouter une colonne avec les numéros de ligne
donnees_normee$Row <- seq_len(nrow(donnees_normee))
# Faire de cette colonne la première colonne
donnees_normee <- donnees_normee[, c(ncol(donnees_normee), 1:(ncol(donnees_normee)-1))]

pander(donnees_normee)

Row	X1	X2	X3	X4	X5	X6	X7	X8
1	1.78	-1.37	-1.09	0.182	0.333	-1.59	-0.378	0.0658
2	0.779	0.77	0.746	-0.83	1.58	-1.4	-0.749	-0.0934
3	0.551	0.206	-1.09	1.8	0.725	0.884	-1.23	-0.206
4	0.518	-0.144	0.682	-0.385	1.45	0.855	-1.37	-0.241
5	-0.552	0.905	0.21	-0.857	-0.309	-0.688	0.684	0.403
6	-1.04	0.147	0.103	0.246	-1.05	-0.0243	1.14	0.546
7	0.548	0.193	-1.09	1.8	0.726	0.877	-1.23	-0.206
8	1.39	-0.758	-0.3	-0.481	0.572	-1.42	-0.413	0.05
9	1.78	-1.36	-1.09	0.187	0.334	-1.59	-0.373	0.0663
10	0.186	0.518	0.189	-0.844	2.17	-1.23	-0.45	-4.74
11	0.861	-0.801	2.89	-1.31	1.02	-1.18	-0.453	0.0408
12	1.18	-1.51	-0.751	-0.249	-0.454	0.925	-0.743	-0.0518
13	0.931	-1.11	-0.695	2.22	0.345	0.893	-1.14	-0.178
14	0.00155	1.04	-0.686	1.02	0.789	0.88	-0.996	-0.132
15	-0.17	1.2	-0.0945	0.000524	1.09	0.868	-0.971	-0.125
16	0.423	-0.535	1.78	-1.31	0.705	0.884	-0.926	-0.11
17	-0.358	0.564	0.176	-0.812	-0.609	-0.724	0.744	0.421
18	-1.29	-0.534	0.0213	-0.948	-1.04	-1.11	2.09	0.859
19	-1.4	-0.00137	1.31	-1.13	-1.05	-1.02	1.93	0.798
20	-0.838	0.488	1.06	-0.499	-1.05	-0.0645	0.892	0.468
21	0.0532	0.445	-1.09	0.00506	-1.05	0.224	-0.143	0.14
22	-0.0312	-1.37	-0.583	1.23	-1.05	0.949	1.17	0.559
23	-1.41	-1.25	-0.669	0.636	-1.05	0.949	1.18	0.559
24	-1.23	2.15	-0.922	0.591	-1.05	0.949	0.382	0.306
25	-1.53	1.96	-0.112	0.881	-1.05	0.949	0.582	0.37
26	-1.14	0.163	1.11	-1.13	-1.05	0.949	0.767	0.429

ACP from scratch CENTREE + NORMEE

result_centree <- acp_from_scratch(data_PDE19, nombre_composants=8, reduire = FALSE)
result_normee <- acp_from_scratch(data_PDE19, nombre_composants=8, reduire = TRUE)

ACP avec dudi.pca CENTREE + NORMEE

result_dudi_centree <- dudi.pca(data_PDE19,center = TRUE, scale = FALSE, nf = 8, scannf = FALSE)
result_dudi_normee <- dudi.pca(data_PDE19,center = TRUE, scale = TRUE, nf = 8, scannf = FALSE)

library(ggplot2)

# Exécution de l'ACP - centrée et normée
result_centree <- acp_from_scratch(data_PDE19, nombre_composants=8, reduire = FALSE)
result_normee <- acp_from_scratch(data_PDE19, nombre_composants=8, reduire = TRUE)

# Créer une fonction pour générer un graphique ggplot pour une paire spécifique de composantes comme demandé
generate_plot <- function(result, x, y, title) {
  data_plot <- as.data.frame(result$nouvelles_coordonnees)
  
  ggplot(data_plot, aes_string(x = paste0("V", x), y = paste0("V", y))) +
    geom_point() +
    labs(title = title, x = paste0("PC", x), y = paste0("PC", y)) +
    theme_minimal()
}

# Générer et afficher les graphiques pour chaque paire de composantes
# Pour données centrées
generate_plot(result_centree, 1, 2, "ACP Centrée: PC1 vs PC2")

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

generate_plot(result_centree, 1, 3, "ACP Centrée: PC1 vs PC3")

generate_plot(result_centree, 1, 4, "ACP Centrée: PC1 vs PC4")

generate_plot(result_centree, 2, 3, "ACP Centrée: PC2 vs PC3")

# Pour données centrées et normées
generate_plot(result_normee, 1, 2, "ACP Normée: PC1 vs PC2")

generate_plot(result_normee, 1, 3, "ACP Normée: PC1 vs PC3")

generate_plot(result_normee, 1, 4, "ACP Normée: PC1 vs PC4")

generate_plot(result_normee, 2, 3, "ACP Normée: PC2 vs PC3")

On observe qu’il semble y avoir deux groupes distincts, mais cette distinction devient moins nette au fur et à mesure que l’on avance dans les Composantes Principales (CP), car l’inertie captée diminue avec ces dernières.

Il est également notable que, dans le cas des données CENTRÉES, des valeurs extrêmes captent presque toute la variance, au détriment du reste des données. Cela souligne l’importance de NORMALISER les données, particulièrement si les variances des données originales ne sont pas similaires.

D’autres commentaires sur la différence entre les données CENTRÉES et NORMALISÉES seront abordés ultérieurement.

library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

# Exécution de l'ACP - centrée et normée
result_centree <- acp_from_scratch(data_PDE19, nombre_composants=3, reduire = FALSE)
result_normee <- acp_from_scratch(data_PDE19, nombre_composants=3, reduire = TRUE)

# Fonction pour générer un graphique Plotly 3D pour l'ACP
generate_3d_plot <- function(result, title) {
  data_plot <- as.data.frame(result$nouvelles_coordonnees)
  
  plot_ly(data = data_plot, x = ~V1, y = ~V2, z = ~V3, 
          type = 'scatter3d', mode = 'markers',
          marker = list(size = 2, opacity = 0.7)) %>%
    layout(title = title,
           scene = list(xaxis = list(title = 'PC1'),
                        yaxis = list(title = 'PC2'),
                        zaxis = list(title = 'PC3')))
}

# Générer les graphiques 3D pour les données centrées et normées
p3d_centree <- generate_3d_plot(result_centree, "ACP Centrée 3D")
p3d_normee <- generate_3d_plot(result_normee, "ACP Normée 3D")

# Afficher les graphiques
p3d_centree

p3d_normee

Comparaison des résultats DUDI et FROM SCRATCH :

Comparaison des valeurs propres CENTREE

# Mettre l'option pour avoir 5 chiffres après la virgule
panderOptions('digits', 5)
pander(t(as.matrix(result_dudi_centree$eig)))

909.47

351.32

234.58

60.747

18.652

5.2277

3.3892

0.91294

pander(t(as.matrix(result_centree$valeurs_propres)))

909.47

351.32

234.58

60.747

18.652

5.2277

3.3892

0.91294

On obtient exactement les mêmes valeurs propres entre pca_from_scratch et dudi.pca. Les valeurs propres dans une ACP reflètent la quantité de variance capturée par chaque composante principale. Dans une ACP centrée, les valeurs propres peuvent varier considérablement, en particulier si les variances des différentes variables sont très différentes. En revanche, dans une ACP normalisée, puisque toutes les variables sont mises à l’échelle pour avoir la même variance unitaire, les valeurs propres tendent à être plus petites et moins dispersées.

Comparaison des valeurs propres NORMEE

pander(t(as.matrix(result_dudi_normee$eig)))

2.8725

2.1191

1.3782

0.79864

0.51113

0.21485

0.090138

0.015465

pander(t(as.matrix(result_normee$valeurs_propres)))

2.8725

2.1191

1.3782

0.79864

0.51113

0.21485

0.090138

0.015465

On obtient exactement les mêmes valeurs propres entre pca_from_scratch et dudi.pca

Dans cette ACP NORMEE nos valeurs propres sont bien moins dispersées (ordre de grandeur 1 contre ordre 3)

Comparaison des vecteurs propres CENTREE

# Mettre l'option pour avoir 3 chiffres après la virgule
panderOptions('digits', 3)
pander(as.matrix(result_dudi_centree$c1))

	CS1	CS2	CS3	CS4	CS5	CS6	CS7	CS8
X1	0.605	0.461	-0.359	0.297	0.226	0.00608	-0.391	0.0171
X2	-0.0501	-0.146	0.193	-0.565	0.605	-0.0109	-0.504	0.0326
X3	-0.0104	-0.017	-0.023	-0.101	-0.309	0.641	-0.231	0.656
X4	0.00983	0.0242	0.0689	0.0874	0.141	-0.623	0.151	0.745
X5	0.281	-0.0862	0.0114	-0.386	-0.67	-0.412	-0.369	-0.1
X6	-0.0608	0.049	0.756	0.511	-0.0795	-0.0182	-0.39	-0.0516
X7	-0.591	-0.199	-0.503	0.314	-0.0442	-0.167	-0.478	-0.0103
X8	-0.446	0.846	0.0651	-0.258	-0.103	-0.0546	-0.0101	-0.0215

pander(as.matrix(result_centree$vecteurs_propres))

0.605	0.461	-0.359
-0.0501	-0.146	0.193
-0.0104	-0.017	-0.023
0.00983	0.0242	0.0689
0.281	-0.0862	0.0114
-0.0608	0.049	0.756
-0.591	-0.199	-0.503
-0.446	0.846	0.0651

On obtient exactement les mêmes vecteurs propres entre pca_from_scratch et dudi.pca

Comparaison des vecteurs propres NORMEE

pander(as.matrix(result_dudi_normee$c1))

	CS1	CS2	CS3	CS4	CS5	CS6	CS7	CS8
X1	0.512	0.0249	0.36	-0.0863	0.241	-0.203	0.434	0.559
X2	-0.154	-0.0609	-0.728	0.196	0.556	-0.0581	0.212	0.217
X3	-0.098	0.535	-0.189	-0.566	-0.122	0.43	0.384	0.0099
X4	0.106	-0.63	0.0505	0.094	-0.0228	0.717	0.258	-0.0185
X5	0.527	0.158	-0.237	-0.0711	0.0639	0.336	-0.669	0.266
X6	-0.095	-0.495	-0.285	-0.554	-0.387	-0.3	-0.0868	0.333
X7	-0.52	0.158	0.164	0.331	-0.227	0.212	-0.112	0.677
X8	-0.368	-0.135	0.371	-0.451	0.642	0.0941	-0.29	0.0172

pander(as.matrix(result_normee$vecteurs_propres))

0.512	0.0249	0.36
-0.154	-0.0609	-0.728
-0.098	0.535	-0.189
0.106	-0.63	0.0505
0.527	0.158	-0.237
-0.095	-0.495	-0.285
-0.52	0.158	0.164
-0.368	-0.135	0.371

CCL : On obtient exactement les mêmes vecteurs propres entre pca_from_scratch et dudi.pca OUF!

Partie 2 - qualité de l’ACP

3. Représentation de l’inertie expliquée et cumulée CENTREE et NORMEE

# Définissez d'abord vos données pour la centrée
inertie_expliquee_centree <- result_centree$valeurs_propres / sum(result_centree$valeurs_propres)
cumul_inertie_centree <- cumsum(inertie_expliquee_centree)
ylim_max_centree <- max(cumul_inertie_centree) + .1

# Puis pour la normalisée
inertie_expliquee_normee <- result_normee$valeurs_propres / sum(result_normee$valeurs_propres)
cumul_inertie_normee <- cumsum(inertie_expliquee_normee)
ylim_max_normee <- max(cumul_inertie_normee) + .1

# Définir les paramètres de la fenêtre graphique pour les subplots
par(mfrow=c(1, 2))  # 1 ligne, 2 colonnes

# Premier subplot pour la centrée
midpoints_centree <- barplot(inertie_expliquee_centree, main="CENTREE", xlab="Composantes", ylab="Inertie expliquée et cumulée", ylim=c(0, ylim_max_centree), col = "blue")
lines(midpoints_centree, cumul_inertie_centree, type="b", col="red")
abline(h = 1, col="grey", lty=2)  # Ligne pointillée pour l'asymptote

# Second subplot pour la normalisée
midpoints_normee <- barplot(inertie_expliquee_normee, main="NORMEE", xlab="Composantes", ylab="Inertie expliquée et cumulée", ylim=c(0, ylim_max_normee), col = "blue")
lines(midpoints_normee, cumul_inertie_normee, type="b", col="red")
abline(h = 1, col="grey", lty=2)  # Ligne pointillée pour l'asymptote

# Réinitialiser les paramètres de la fenêtre graphique
par(mfrow=c(1, 1))

On voit que la composante 1 est plus elevée en CENTREE qu’en NORMEE: Interpretation : l’inertie expliqué par cette composante en CENTREE reflete une variance capturé dans les données originales, pondérées selon l’importance relative des différentes variables selon leurs variances et que en NORMEE, la ou toutes les variables sont à la même echelle, cette variance capturé dans les données originales n’est plus presente.

Règle de Karlis - Saporta - Spinaki NORMEE

sueil = mean(VP) + 2sd(VP)

# Seuil calculé à partir de la moyenne et de l'écart-type des valeurs propres
karlis_threshold_value <- mean(result_normee$valeurs_propres) + 2 * sd(result_normee$valeurs_propres)
cat("Valeur de karlis_sueil plus grande que la première VP :", karlis_threshold_value, ">", result_normee$valeurs_propres[1])

## Valeur de karlis_sueil plus grande que la première VP : 3.083225 > 2.872476

La regle de Karlis - Saporta - Spinaki pour NORMEE (selon le cours) nous donne un sueil plus grand que la plus grande des VP: on ne peut pas l’utiliser ici.

Règle de Kaiser - Guttman NORMEE

sueil = 1

# Calcul des inertie expliquées et cumulées
inertie_expliquee_normee <- result_normee$valeurs_propres / sum(result_normee$valeurs_propres)
cumul_inertie_normee <- cumsum(inertie_expliquee_normee)
ylim_max_normee <- max(cumul_inertie_normee) + 0.1

# Barplot de l'inertie expliquée pour la normalisée
midpoints_normee <- barplot(inertie_expliquee_normee, main="ACP NORMEE + règles classiques", xlab="Composantes", ylab="Inertie expliquée", ylim=c(0, ylim_max_normee), col = "blue")
lines(midpoints_normee, cumul_inertie_normee, type="b", pch=19, col="red")

# Ajouter les valeurs des valeurs propres en abscisse, écriture verticale
text(x=midpoints_normee, y=-0.05, labels=round(result_normee$valeurs_propres, 2), srt=90, adj=1, xpd=TRUE, cex=0.8)

# Règle de Kaiser - Guttman
# Identifier la dernière composante avec une valeur propre >= 1 et la position de la barre correspondante
kaiser_comps <- which(result_normee$valeurs_propres >= 1)
if(length(kaiser_comps) > 0) {
  kaiser_limit <- max(kaiser_comps)
  # Utilisez les midpoints pour aligner la ligne verticale avec la barre correspondante
  abline(v=midpoints_normee[kaiser_limit], col="green", lty=2)
}
abline(v=0, col="orange", lty=2)
abline(h = .8, col="grey", lty=2)
# Ajouter une légende pour les seuils
legend("topright",inset=c(0, 0.2), legend=c("Kaiser-Guttman", "Karlis-Saporta-Spinaki"), col=c("green", "orange"), lty=2, cex=0.8)

Selon cette regle de Kaiser - Guttman, on devrait prendre les 3 premieres composantes.

Selon la règle des 80% d’inertei cumulée, 3 est celui a choisir egalement.

On ne voit pas de “Coude” apparaitre sur la courbe des inerties cumuléee.

Il semble que 3 soit un bon choix, mais cela depend de notre but évidemment.

4. Nouvelles coordonnées calculées matriciellement

Comparaison des nouvelles coordonnees CENTREE

# Mettre l'option pour avoir deux chiffres après la virgule
panderOptions('digits', 4)
pander(as.matrix(result_dudi_centree$li))

Axis1	Axis2	Axis3	Axis4	Axis5	Axis6	Axis7	Axis8
29.24	20.15	-26.23	2.579	3.392	-1.333	0.7113	-0.2791
24.7	5.899	-10.63	-17.34	-1.903	-1.054	-1.372	-0.2002
24.73	6.456	16.88	-0.3885	0.7747	-2.807	0.8349	0.602
28.31	5.604	17.4	-3.671	-7.238	0.4854	-0.5048	-0.9967
-19.52	-2.061	-7.613	-8.231	1.843	-0.002059	-1.094	-0.669
-34.63	-4.205	-2.945	0.8816	0.2715	-0.5689	0.9644	0.7654
24.71	6.438	16.82	-0.4085	0.7076	-2.807	0.9399	0.6018
25.05	15.44	-20.59	-2.408	1.57	-0.1851	-0.1289	-0.4921
29.12	20.1	-26.24	2.578	3.366	-1.361	0.6875	-0.2756
58.47	-83.42	-14.13	5.871	0.4236	0.002492	-0.1191	0.001939
19.78	9.921	-14.3	-6.851	-6.756	3.977	0.04525	2.258
22.9	15.93	5.146	16.47	2.703	3.608	-0.04811	-1.444
27.71	11.92	11.53	8.87	-0.6275	-1.458	2.57	2.083
14.16	0.579	19.7	-6.449	0.2354	-2.374	-0.3699	-0.2051
12.35	-1.545	20.72	-9.548	-2.614	-1.394	-1.41	-1.292
18.72	6.181	13.35	1.809	-6.227	4.753	-1.009	-0.4128
-18.59	0.4923	-10.46	-4.535	3.303	0.9167	-0.8359	-0.3595
-50.99	-4.965	-21.38	-0.9477	-4.497	-2.329	1.108	-1.206
-50.18	-6.968	-17.65	-4.07	-3.609	0.3853	0.5246	0.5405
-28.58	-3.109	-2.197	-0.9118	2.11	2.8	-0.05777	1.269
-2.063	4.175	3.852	2.567	9.228	2.757	2.609	-0.7143
-22.53	7.806	-3.566	19.64	-1.33	-3.071	-6.593	0.7255
-40.24	-5.759	6.918	10.48	-7.549	-2.58	4.203	-0.8556
-27.26	-8.609	17.6	-4.864	8.597	-0.1808	-1.08	-0.2255
-34.1	-11.09	17.75	-5.285	5.527	-0.1417	-0.1631	1.266
-31.28	-5.364	10.3	4.163	-1.702	3.964	-0.4117	-0.4847

pander(as.matrix(result_centree$nouvelles_coordonnees))

29.24	20.15	-26.23
24.7	5.899	-10.63
24.73	6.456	16.88
28.31	5.604	17.4
-19.52	-2.061	-7.613
-34.63	-4.205	-2.945
24.71	6.438	16.82
25.05	15.44	-20.59
29.12	20.1	-26.24
58.47	-83.42	-14.13
19.78	9.921	-14.3
22.9	15.93	5.146
27.71	11.92	11.53
14.16	0.579	19.7
12.35	-1.545	20.72
18.72	6.181	13.35
-18.59	0.4923	-10.46
-50.99	-4.965	-21.38
-50.18	-6.968	-17.65
-28.58	-3.109	-2.197
-2.063	4.175	3.852
-22.53	7.806	-3.566
-40.24	-5.759	6.918
-27.26	-8.609	17.6
-34.1	-11.09	17.75
-31.28	-5.364	10.3

On obtient les mêmes nouvelles coordonnées par ACP_from_scratch et par DUDI.ACP On pouvais savoir cela en avance, car on a vue que les vecteurs propres étaient les mêmes, donc forcement, les nouvelles coordonnées aussi.

Qualité de la projection de l’individu Qik=1 CENTREE

result_centree <- acp_from_scratch(data_PDE19, nombre_composants=NULL, reduire = FALSE)
result_normee <- acp_from_scratch(data_PDE19, nombre_composants=NULL, reduire = TRUE)
#On prend toute les composant pour calculer la qualité de projection.

nouvelles_coordonnees <- as.matrix(result_centree$nouvelles_coordonnees)

k <- 1

# Calculer la qualité de la projection Qik pour chaque individu
qualite_projections <- apply(nouvelles_coordonnees, 1, function(Ci) sum(Ci[1:k]^2) / sum(Ci^2))

# Afficher les résultats avec deux chiffres après la virgule
panderOptions('digits', 2)
pander(t(as.matrix(qualite_projections)))

Table continues below

0.43

0.57

0.65

0.67

0.74

0.98

0.65

0.48

0.43

0.32

0.48

Table continues below

0.48

0.68

0.32

0.22

0.55

0.71

0.84

0.87

0.97

0.03

0.5

0.86

0.61

0.7

0.85

On voit des valeurs avec grande qualité de projection (comme la derniére) et d’autres médiocres (comme la 21 avec 0.03 car plus proche de l’origine notamment)

Qualité de la projection de l’individu Qik=2 CENTREE

result_centree <- acp_from_scratch(data_PDE19, nombre_composants=NULL, reduire = FALSE)
result_normee <- acp_from_scratch(data_PDE19, nombre_composants=NULL, reduire = TRUE)
#On prend toute les composant pour calculer la qualité de projection.

nouvelles_coordonnees <- as.matrix(result_centree$nouvelles_coordonnees)

k <- 2 

# Calculer la qualité de la projection Qik pour chaque individu
qualite_projections <- apply(nouvelles_coordonnees, 1, function(Ci) sum(Ci[1:k]^2) / sum(Ci^2))

# Afficher les résultats avec deux chiffres après la virgule
panderOptions('digits', 2)
pander(t(as.matrix(qualite_projections)))

Table continues below

0.64

0.61

0.69

0.69

0.75

0.99

0.69

0.67

0.64

0.98

0.61

Table continues below

0.71

0.8

0.32

0.23

0.61

0.71

0.84

0.88

0.98

0.15

0.56

0.87

0.67

0.77

0.88

En passant de Qik = 1 à Qik = 2, la qualité de la projection s’améliore. Cette amélioration est logique, car l’accès à une plus grande part de la variance, comme expliqué dans la suite du texte.

Qualité de la projection de l’individu Qik=3 CENTREE

result_centree <- acp_from_scratch(data_PDE19, nombre_composants=NULL, reduire = FALSE)
result_normee <- acp_from_scratch(data_PDE19, nombre_composants=NULL, reduire = TRUE)
#On prend toute les composant pour calculer la qualité de projection.

nouvelles_coordonnees <- as.matrix(result_centree$nouvelles_coordonnees)

k <- 3

# Calculer la qualité de la projection Qik pour chaque individu
qualite_projections <- apply(nouvelles_coordonnees, 1, function(Ci) sum(Ci[1:k]^2) / sum(Ci^2))

# Afficher les résultats avec deux chiffres après la virgule
panderOptions('digits', 2)
pander(t(as.matrix(qualite_projections)))

Table continues below

0.99

0.71

0.99

0.94

0.86

1

0.99

0.99

0.99

1

0.86

0.73

Table continues below

0.92

0.93

0.85

0.9

0.93

0.99

0.99

0.98

0.26

0.57

0.9

0.92

0.96

0.97

Pour k=3, la qualité de projection s’améliore car une plus grande proportion d’individus présente des valeurs proches de 1: Lorsque on augmente le nombre de composantes principales (de k=2 à k=3, ici), on inclut plus d’axes pour représenter nos données. Chaque axe supplémentaire est capable de capturer une partie de la variance qui n’était pas capturée par les axes précédents. En conséquence, la somme des variances capturées par les composantes principales augmente, ce qui conduit à une amélioration de la qualité de projection pour les individus. En effet d’un autre point de vue, plus il y a de dimensions dans l’espace de projection, plus il est probable que les individus soient projetés de façon à conserver leur distance relative, et donc leur variance, par rapport à l’espace original des données.

En calculant le COS2 (ce qui est la même chose mathématiquement) k=2 CENTREE

result_centree <- acp_from_scratch(data_PDE19, nombre_composants=NULL, reduire = FALSE)
result_normee <- acp_from_scratch(data_PDE19, nombre_composants=NULL, reduire = TRUE)
#On prend toute les composant pour calculer la qualité de projection.

nouvelles_coordonnees <- as.matrix(result_centree$nouvelles_coordonnees)

k <- 2 

# Calculer la qualité de la projection pour chaque individu (cos2)
qualite_cos2 <- apply(nouvelles_coordonnees, 1, function(Ci) {
  cos2_individual_axes <- Ci^2 / sum(Ci^2)
  sum(cos2_individual_axes[1:k])
})

# Afficher les résultats avec deux chiffres après la virgule
panderOptions('digits', 2)
pander(t(as.matrix(qualite_cos2)))

Table continues below

0.64

0.61

0.69

0.69

0.75

0.99

0.69

0.67

0.64

0.98

0.61

Table continues below

0.71

0.8

0.32

0.23

0.61

0.71

0.84

0.88

0.98

0.15

0.56

0.87

0.67

0.77

0.88

On obtient les mêmes valeurs de qualitées de projecton qu’avec la formule de définiton

En calculant le COS2 (ce qui est la même chose mathématiquement) k=3

nouvelles_coordonnees <- as.matrix(result_centree$nouvelles_coordonnees)
k <- 3

# Calculer la qualité de la projection pour chaque individu (cos2)
qualite_cos2 <- apply(nouvelles_coordonnees, 1, function(Ci) {
  cos2_individual_axes <- Ci^2 / sum(Ci^2)
  sum(cos2_individual_axes[1:k])
})

# Afficher les résultats avec deux chiffres après la virgule
panderOptions('digits', 2)
pander(t(as.matrix(qualite_cos2)))

Table continues below

0.99

0.71

0.99

0.94

0.86

1

0.99

0.99

0.99

1

0.86

0.73

Table continues below

0.92

0.93

0.85

0.9

0.93

0.99

0.99

0.98

0.26

0.57

0.9

0.92

0.96

0.97

On obtient les mêmes valeurs de qualitées de projecton qu’avec la formule de définiton

Comparaison des nouvelles coordonnees NORMEE

pander(as.matrix(result_dudi_normee$li))

Axis1	Axis2	Axis3	Axis4	Axis5	Axis6	Axis7	Axis8
1.8	0.2	2.2	0.92	0.57	-0.11	0.05	-0.0089
1.5	1.7	-0.6	0.04	1.3	0.31	-0.25	0.071
1.6	-2.2	-0.36	-0.074	0.19	0.39	-0.039	-0.039
1.7	0.25	-0.76	-1.4	-0.11	-0.16	-0.47	0.014
-1.1	0.93	-0.41	0.47	0.71	-0.18	-0.11	0.065
-1.9	-0.18	0.16	0.3	-0.15	0.37	0.1	-0.059
1.6	-2.2	-0.35	-0.073	0.19	0.4	-0.043	-0.046
1.4	0.94	1.3	0.44	0.67	-0.18	-0.024	0.022
1.7	0.2	2.2	0.92	0.56	-0.1	0.048	-0.0071
3.1	2.1	-2.4	2.4	-2	-0.034	0.12	0.016
1	3.1	0.32	-1.6	0.091	0.78	0.44	-0.067
0.96	-0.76	1.4	-0.7	-0.71	-1.2	0.16	0.012
1.7	-2.2	0.8	-0.51	-0.54	0.76	0.34	-0.16
0.92	-1.5	-1.2	-0.13	0.49	0.15	-0.23	0.026
0.78	-0.52	-1.6	-0.52	0.51	-0.19	-0.51	0.083
0.8	1.4	-0.47	-2.1	-0.54	-0.47	0.00083	0.016
-1.2	0.88	0.016	0.47	0.57	-0.26	0.088	0.047
-2.5	1.2	1.1	0.9	-0.15	0.13	-0.6	-0.037
-2.7	1.9	0.35	0.19	-0.07	0.47	-0.055	-0.044
-1.8	0.78	-0.3	-0.29	0.0075	0.14	0.49	-0.031
-0.48	-0.93	0.11	0.54	0.36	-0.94	0.36	-0.18
-1.1	-1.5	1.5	-0.13	-1.1	0.38	0.11	0.5
-1.9	-1.2	0.94	0.0032	-1.3	0.19	-0.65	-0.23
-1.8	-1.6	-1.6	0.64	0.67	-0.38	0.21	0.061
-2.1	-1.4	-1.7	0.24	0.38	0.3	0.38	-0.013
-2	0.7	-0.53	-1	-0.63	-0.54	0.083	-0.0069

pander(as.matrix(result_normee$nouvelles_coordonnees))

1.8	0.2	2.2	0.92	0.57	-0.11	0.05	-0.0089
1.5	1.7	-0.6	0.04	1.3	0.31	-0.25	0.071
1.6	-2.2	-0.36	-0.074	0.19	0.39	-0.039	-0.039
1.7	0.25	-0.76	-1.4	-0.11	-0.16	-0.47	0.014
-1.1	0.93	-0.41	0.47	0.71	-0.18	-0.11	0.065
-1.9	-0.18	0.16	0.3	-0.15	0.37	0.1	-0.059
1.6	-2.2	-0.35	-0.073	0.19	0.4	-0.043	-0.046
1.4	0.94	1.3	0.44	0.67	-0.18	-0.024	0.022
1.7	0.2	2.2	0.92	0.56	-0.1	0.048	-0.0071
3.1	2.1	-2.4	2.4	-2	-0.034	0.12	0.016
1	3.1	0.32	-1.6	0.091	0.78	0.44	-0.067
0.96	-0.76	1.4	-0.7	-0.71	-1.2	0.16	0.012
1.7	-2.2	0.8	-0.51	-0.54	0.76	0.34	-0.16
0.92	-1.5	-1.2	-0.13	0.49	0.15	-0.23	0.026
0.78	-0.52	-1.6	-0.52	0.51	-0.19	-0.51	0.083
0.8	1.4	-0.47	-2.1	-0.54	-0.47	0.00083	0.016
-1.2	0.88	0.016	0.47	0.57	-0.26	0.088	0.047
-2.5	1.2	1.1	0.9	-0.15	0.13	-0.6	-0.037
-2.7	1.9	0.35	0.19	-0.07	0.47	-0.055	-0.044
-1.8	0.78	-0.3	-0.29	0.0075	0.14	0.49	-0.031
-0.48	-0.93	0.11	0.54	0.36	-0.94	0.36	-0.18
-1.1	-1.5	1.5	-0.13	-1.1	0.38	0.11	0.5
-1.9	-1.2	0.94	0.0032	-1.3	0.19	-0.65	-0.23
-1.8	-1.6	-1.6	0.64	0.67	-0.38	0.21	0.061
-2.1	-1.4	-1.7	0.24	0.38	0.3	0.38	-0.013
-2	0.7	-0.53	-1	-0.63	-0.54	0.083	-0.0069

Qualité de la projection de l’individu Qik=2 NORMEE

result_centree <- acp_from_scratch(data_PDE19, nombre_composants=NULL, reduire = FALSE)
result_normee <- acp_from_scratch(data_PDE19, nombre_composants=NULL, reduire = TRUE)
#On prend toute les composant pour calculer la qualité de projection.

nouvelles_coordonnees <- as.matrix(result_normee$nouvelles_coordonnees)


k <- 2 

# Calculer la qualité de la projection Qik pour chaque individu
qualite_projections <- apply(nouvelles_coordonnees, 1, function(Ci) sum(Ci[1:k]^2) / sum(Ci^2))

# Afficher les résultats avec deux chiffres après la virgule
panderOptions('digits', 2)
pander(t(as.matrix(qualite_projections)))

Table continues below

0.34

0.71

0.96

0.5

0.69

0.93

0.96

0.55

0.34

0.48

0.76

Table continues below

0.25

0.8

0.63

0.2

0.33

0.78

0.77

0.97

0.9

0.43

0.47

0.62

0.61

0.65

0.7

Qualité de la projection de l’individu Qik=3 NORMEE

nouvelles_coordonnees <- as.matrix(result_normee$nouvelles_coordonnees)

k <- 3

# Calculer la qualité de la projection Qik pour chaque individu
qualite_projections <- apply(nouvelles_coordonnees, 1, function(Ci) sum(Ci[1:k]^2) / sum(Ci^2))

# Afficher les résultats avec deux chiffres après la virgule
panderOptions('digits', 2)
pander(t(as.matrix(qualite_projections)))

Table continues below

0.87

0.75

0.97

0.6

0.75

0.93

0.97

0.87

0.87

0.67

0.77

Table continues below

0.57

0.87

0.93

0.81

0.36

0.78

0.88

0.98

0.92

0.43

0.79

0.73

0.89

0.95

0.74

Même commentaire que pour CENTREE

En calculant le COS2 (ce qui est la même chose mathématiquement) k=2 NORMEE

nouvelles_coordonnees <- as.matrix(result_normee$nouvelles_coordonnees)
k <- 2 

# Calculer la qualité de la projection pour chaque individu (cos2)
qualite_cos2 <- apply(nouvelles_coordonnees, 1, function(Ci) {
  cos2_individual_axes <- Ci^2 / sum(Ci^2)
  sum(cos2_individual_axes[1:k])
})

# Afficher les résultats avec deux chiffres après la virgule
panderOptions('digits', 2)
pander(t(as.matrix(qualite_cos2)))

Table continues below

0.34

0.71

0.96

0.5

0.69

0.93

0.96

0.55

0.34

0.48

0.76

Table continues below

0.25

0.8

0.63

0.2

0.33

0.78

0.77

0.97

0.9

0.43

0.47

0.62

0.61

0.65

0.7

En calculant le COS2 (ce qui est la même chose mathématiquement) k=3 NORMEE

nouvelles_coordonnees <- as.matrix(result_normee$nouvelles_coordonnees)
k <- 3

# Calculer la qualité de la projection pour chaque individu (cos2)
qualite_cos2 <- apply(nouvelles_coordonnees, 1, function(Ci) {
  cos2_individual_axes <- Ci^2 / sum(Ci^2)
  sum(cos2_individual_axes[1:k])
})

# Afficher les résultats avec deux chiffres après la virgule
panderOptions('digits', 2)
pander(t(as.matrix(qualite_cos2)))

Table continues below

0.87

0.75

0.97

0.6

0.75

0.93

0.97

0.87

0.87

0.67

0.77

Table continues below

0.57

0.87

0.93

0.81

0.36

0.78

0.88

0.98

0.92

0.43

0.79

0.73

0.89

0.95

0.74

Meme resultat entre fomule du COS2 et celle du cours pour les diffenrents k Même commentaire que précédemment pour CENTREE

5. Contribution de l’individu i à l’inertie de l’axe factoriel j CENTREE et NORMEE

Contribution cas CENTREE

nouvelles_coordonnees <- as.matrix(result_centree$nouvelles_coordonnees)

valeurs_propres <- result_centree$valeurs_propres 

# Nombre total d'individus
n <- nrow(nouvelles_coordonnees)

# Calculer la contribution de chaque individu à l'inertie de chaque axe
contributions <- sweep(nouvelles_coordonnees^2, 2, valeurs_propres, FUN="/") / n

# Affichage des résultats avec deux chiffres après la virgule
panderOptions('digits', 2)
pander(contributions)

0.036	0.044	0.11	0.0042	0.024	0.013	0.0057	0.0033
On reviend	ra sur ses	contribut	ions aprés	.
# Contribu	tion cas N	ORMEE
```r nouvelles_	coordonnee	s <- as.ma	trix(resul	t_normee$n	ouvelles_c	oordonnees	)
valeurs_pr	opres <- r	esult_norm	ee$valeurs	_propres
# Nombre t n <- nrow(	otal d’ind nouvelles_	ividus coordonnee	s)
# Calculer contributi	la contri ons <- swe	bution de ep(nouvell	chaque ind es_coordon	ividu à l’ nees^2, 2,	inertie de valeurs_p	chaque ax ropres, FU	e N=“/”) / n
# Affichag panderOpti pander(con ```	e des résu ons(’digit tributions	ltats avec s’, 2) )	deux chif	fres après	la virgul	e

0.041 0.00072 0.13 0.041 0.024 0.002 0.0011 2e-04

0.03 0.054 0.0099 7.6e-05 0.13 0.017 0.027 0.013

0.033 0.088 0.0036 0.00026 0.0027 0.028 0.00066 0.0039

0.037 0.0011 0.016 0.097 0.00094 0.0046 0.093 0.00049

0.017 0.016 0.0047 0.011 0.038 0.0057 0.0056 0.01

0.047 6e-04 0.00069 0.0044 0.0017 0.024 0.0044 0.0088

0.033 0.088 0.0034 0.00025 0.0026 0.028 8e-04 0.0053

0.028 0.016 0.047 0.0094 0.034 0.0056 0.00024 0.0012

0.041 7e-04 0.13 0.041 0.024 0.0018 0.001 0.00013

0.13 0.082 0.16 0.28 0.3 2e-04 0.0066 0.00064

0.014 0.18 0.0028 0.12 0.00062 0.11 0.082 0.011

0.012 0.011 0.051 0.024 0.038 0.28 0.011 0.00033

0.039 0.09 0.018 0.013 0.022 0.1 0.051 0.06

0.011 0.042 0.042 0.00077 0.018 0.0042 0.024 0.0017

0.0082 0.005 0.074 0.013 0.02 0.0064 0.11 0.017

0.0085 0.034 0.0062 0.21 0.022 0.039 2.9e-07 6e-04

0.018 0.014 7.5e-06 0.01 0.024 0.012 0.0033 0.0055

0.086 0.026 0.034 0.039 0.0017 0.0029 0.15 0.0033

0.099 0.067 0.0033 0.0018 0.00037 0.04 0.0013 0.0048

0.046 0.011 0.0024 0.004 4.2e-06 0.0036 0.1 0.0024

0.0031 0.016 0.00037 0.014 0.0099 0.16 0.054 0.084

0.016 0.042 0.066 0.00085 0.085 0.026 0.0056 0.62

0.046 0.028 0.025 5e-07 0.13 0.0066 0.18 0.13

0.041 0.049 0.076 0.02 0.034 0.026 0.019 0.0092

0.057 0.034 0.081 0.0027 0.011 0.016 0.061 0.00041

0.056 0.0089 0.008 0.048 0.03 0.053 0.0029 0.00012 ——– ——— ——— ——— ——— ——– ——— ———

Heatmap qui affiche les individus rangés par leur contribution globale à l’inertie des axes cas CENTREE puis cas NORMEE

# Utilisez les variables données
nouvelles_coordonnees <- as.matrix(result_centree$nouvelles_coordonnees)
valeurs_propres <- result_centree$valeurs_propres
n <- nrow(nouvelles_coordonnees)

# Calculer la contribution de chaque individu à l'inertie de chaque axe
# Nous nous assurons que le résultat est une matrice
contributions <- sweep(nouvelles_coordonnees^2, 2, valeurs_propres, FUN="/") / n

# Transformer en data frame pour le traitement
contributions_df <- as.data.frame(contributions)

# Calculer la contribution totale pour chaque individu
contributions_df$Total <- rowSums(contributions_df)

# Trier le data frame par la contribution totale de manière décroissante
contributions_df <- contributions_df[order(-contributions_df$Total), ]

# Préparer les données pour la heatmap (sans la colonne Total)
contributions_for_heatmap <- contributions_df[, -ncol(contributions_df)]


# Utilisez les variables données
nouvelles_coordonnees <- as.matrix(result_normee$nouvelles_coordonnees)
valeurs_propres <- result_normee$valeurs_propres
n <- nrow(nouvelles_coordonnees)

# Calculer la contribution de chaque individu à l'inertie de chaque axe
# Nous nous assurons que le résultat est une matrice
contributions <- sweep(nouvelles_coordonnees^2, 2, valeurs_propres, FUN="/") / n

# Transformer en data frame pour le traitement
contributions_df <- as.data.frame(contributions)

# Calculer la contribution totale pour chaque individu
contributions_df$Total <- rowSums(contributions_df)

# Trier le data frame par la contribution totale de manière décroissante
contributions_df <- contributions_df[order(-contributions_df$Total), ]

# Préparer les données pour la heatmap (sans la colonne Total)
contributions_normee_for_heatmap <- contributions_df[, -ncol(contributions_df)]





# Charger la librairie 'pheatmap'
if (!requireNamespace("pheatmap", quietly = TRUE)) {
  install.packages("pheatmap")
}
library(pheatmap)


# Créer la heatmap pour les données CENTRÉES
pheatmap(contributions_for_heatmap, cluster_rows = FALSE, show_rownames = TRUE,
         cluster_cols = FALSE, display_numbers = TRUE,
         color = colorRampPalette(c("blue", "white", "red"))(100),
         main = "Contributions (Centrées)")

# La heatmap est affichée sans regroupement des lignes ou des colonnes,
# avec les noms des lignes et les nombres affichés.

La CP 2 capture la variance de l’individu 10 principalement. Ici, la somme des colonnes ne vaut pas 1 à cause des arrondis (2 chiffres) mais sinon si ( on le calcul plus tard)

# Créer la heatmap pour les données NORMALISÉES
pheatmap(contributions_normee_for_heatmap, cluster_rows = FALSE, show_rownames = TRUE,
         cluster_cols = FALSE, display_numbers = TRUE,
         color = colorRampPalette(c("blue", "white", "red"))(100),
         main = "Contributions (Normalisées)")

En NORMEE, la CP 2 ne capture plus la variance du 10 autant, car la variance du 10 est mise à la meme echelle que les autres individus. En NORMEE, la variance du 22, mise au meme niveau que les autres, est capturé grandement par le CP8, ce qui n’était pas le cas en CENTREE :

Calculer la somme des contributions pour chaque axe (Est ce egale à 1 ?) CENTREE

# Calculer la somme des contributions pour chaque axe
sum_contributions <- colSums(contributions)

# Afficher les sommes des contributions, cela devrait être égal ou très proche de 1 pour chaque axe
pander(sum_contributions)

1, 1, 1, 1, 1, 1, 1 and 1

Calculer la somme des contributions pour chaque axe (Est ce egale à 1 ?) NORMEE

# Calculer la somme des contributions pour chaque axe
sum_contributions <- colSums(contributions)

# Afficher les sommes des contributions, cela devrait être égal ou très proche de 1 pour chaque axe
pander(sum_contributions)

1, 1, 1, 1, 1, 1, 1 and 1

6. Vérification avec la fonction dudi.pca CENTREE

result_dudi_centree <- dudi.pca(data_PDE19,center = TRUE, scale = FALSE, nf = 8, scannf = FALSE)
pander(as.matrix(result_dudi_centree$c1))

	CS1	CS2	CS3	CS4	CS5	CS6	CS7	CS8
X1	0.61	0.46	-0.36	0.3	0.23	0.0061	-0.39	0.017
X2	-0.05	-0.15	0.19	-0.56	0.6	-0.011	-0.5	0.033
X3	-0.01	-0.017	-0.023	-0.1	-0.31	0.64	-0.23	0.66
X4	0.0098	0.024	0.069	0.087	0.14	-0.62	0.15	0.75
X5	0.28	-0.086	0.011	-0.39	-0.67	-0.41	-0.37	-0.1
X6	-0.061	0.049	0.76	0.51	-0.079	-0.018	-0.39	-0.052
X7	-0.59	-0.2	-0.5	0.31	-0.044	-0.17	-0.48	-0.01
X8	-0.45	0.85	0.065	-0.26	-0.1	-0.055	-0.01	-0.021

pander(as.matrix(result_centree$vecteurs_propres))

0.61	0.46	-0.36	0.3	0.23	0.0061	-0.39	0.017
-0.05	-0.15	0.19	-0.56	0.6	-0.011	-0.5	0.033
-0.01	-0.017	-0.023	-0.1	-0.31	0.64	-0.23	0.66
0.0098	0.024	0.069	0.087	0.14	-0.62	0.15	0.75
0.28	-0.086	0.011	-0.39	-0.67	-0.41	-0.37	-0.1
-0.061	0.049	0.76	0.51	-0.079	-0.018	-0.39	-0.052
-0.59	-0.2	-0.5	0.31	-0.044	-0.17	-0.48	-0.01
-0.45	0.85	0.065	-0.26	-0.1	-0.055	-0.01	-0.021

Les Vecteurs propres pour les 2 fonctions sont les mêmes, donc les plans factoriels sont aussi les mêmes: notre premier plan factoriel est donc correct en comparant avec la fonction R qui résout l’ACP dudi.pca()

6. Vérification avec la fonction dudi.pca NORMEE

result_dudi_normee <- dudi.pca(data_PDE19,center = TRUE, scale = TRUE, nf = 8, scannf = FALSE)
pander(as.matrix(result_dudi_normee$c1))

	CS1	CS2	CS3	CS4	CS5	CS6	CS7	CS8
X1	0.51	0.025	0.36	-0.086	0.24	-0.2	0.43	0.56
X2	-0.15	-0.061	-0.73	0.2	0.56	-0.058	0.21	0.22
X3	-0.098	0.54	-0.19	-0.57	-0.12	0.43	0.38	0.0099
X4	0.11	-0.63	0.051	0.094	-0.023	0.72	0.26	-0.018
X5	0.53	0.16	-0.24	-0.071	0.064	0.34	-0.67	0.27
X6	-0.095	-0.49	-0.29	-0.55	-0.39	-0.3	-0.087	0.33
X7	-0.52	0.16	0.16	0.33	-0.23	0.21	-0.11	0.68
X8	-0.37	-0.14	0.37	-0.45	0.64	0.094	-0.29	0.017

pander(as.matrix(result_normee$vecteurs_propres))

0.51	0.025	0.36	-0.086	0.24	-0.2	0.43	0.56
-0.15	-0.061	-0.73	0.2	0.56	-0.058	0.21	0.22
-0.098	0.54	-0.19	-0.57	-0.12	0.43	0.38	0.0099
0.11	-0.63	0.051	0.094	-0.023	0.72	0.26	-0.018
0.53	0.16	-0.24	-0.071	0.064	0.34	-0.67	0.27
-0.095	-0.49	-0.29	-0.55	-0.39	-0.3	-0.087	0.33
-0.52	0.16	0.16	0.33	-0.23	0.21	-0.11	0.68
-0.37	-0.14	0.37	-0.45	0.64	0.094	-0.29	0.017

Même commentaire : Les Vecteurs propres pour les 2 fonctions sont les mêmes, donc les plans factoriels sont aussi les mêmes: notre premier plan factoriel est donc correct en comparant avec la fonction R qui résout l’ACP dudi.pca()

7. Représentation graphique CP1 vs CP2 CENTREE

# Graphique pour CP1 vs CP2
plot(result_centree$nouvelles_coordonnees[,1], result_centree$nouvelles_coordonnees[,2], 
     xlab="CP1", ylab="CP2", main="Premier vs Deuxième axe principal CENTREE",
     type="n", xlim=c(min(c(result_centree$nouvelles_coordonnees[,1], result_dudi_centree$li[,1])),
                      max(c(result_centree$nouvelles_coordonnees[,1], result_dudi_centree$li[,1]))),
     ylim=c(min(c(result_centree$nouvelles_coordonnees[,2], result_dudi_centree$li[,2])),
            max(c(result_centree$nouvelles_coordonnees[,2], result_dudi_centree$li[,2]))))
points(result_centree$nouvelles_coordonnees[,1], result_centree$nouvelles_coordonnees[,2], col="blue", pch=16)
points(result_dudi_centree$li[,1], result_dudi_centree$li[,2], col="red", pch=17)
legend(x = 0.6, y = 0.5, legend=c("pca_from_scratch", "dudi.pca"),
       col=c("blue", "red"), pch=c(16, 17), cex=0.8, bty="n")

Les individus se superposent exactement entre ACP FROM SCRATCH et DUDI.APC, on obtient les memes résultats pour l’ACP

Représentation graphique CP1 vs CP3 CENTREE

# Graphique pour CP1 vs CP3
plot(result_centree$nouvelles_coordonnees[,1], result_centree$nouvelles_coordonnees[,3], 
     xlab="CP1", ylab="CP3", main="Premier vs Troisième axe principal CENTREE",
     type="n", xlim=c(min(c(result_centree$nouvelles_coordonnees[,1], result_dudi_centree$li[,1])),
                      max(c(result_centree$nouvelles_coordonnees[,1], result_dudi_centree$li[,1]))),
     ylim=c(min(c(result_centree$nouvelles_coordonnees[,3], result_dudi_centree$li[,3])),
            max(c(result_centree$nouvelles_coordonnees[,3], result_dudi_centree$li[,3]))))
points(result_centree$nouvelles_coordonnees[,1], result_centree$nouvelles_coordonnees[,3], col="blue", pch=16)
points(result_dudi_centree$li[,1], result_dudi_centree$li[,3], col="red", pch=17)
legend(x = 0.6, y = 0.5, legend=c("pca_from_scratch", "dudi.pca"),
       col=c("blue", "red"), pch=c(16, 17), cex=0.8, bty="n")

Représentation graphique CP1 vs CP2 NORMEE

# Graphique pour CP1 vs CP2
plot(result_normee$nouvelles_coordonnees[,1], result_normee$nouvelles_coordonnees[,2], 
     xlab="CP1", ylab="CP2", main="Premier vs Deuxième axe principal NORMEE",
     type="n", xlim=c(min(c(result_normee$nouvelles_coordonnees[,1], result_dudi_normee$li[,1])),
                      max(c(result_normee$nouvelles_coordonnees[,1], result_dudi_normee$li[,1]))),
     ylim=c(min(c(result_normee$nouvelles_coordonnees[,2], result_dudi_normee$li[,2])),
            max(c(result_normee$nouvelles_coordonnees[,2], result_dudi_normee$li[,2]))))
points(result_normee$nouvelles_coordonnees[,1], result_normee$nouvelles_coordonnees[,2], col="blue", pch=16)
points(result_dudi_normee$li[,1], result_dudi_normee$li[,2], col="red", pch=17)
legend(x = 0.6, y = 0.5, legend=c("pca_from_scratch", "dudi.pca"),
       col=c("blue", "red"), pch=c(16, 17), cex=0.8, bty="n")

7. Représentation graphique CP1 vs CP3 NORMEE

# Graphique pour CP1 vs CP3
plot(result_normee$nouvelles_coordonnees[,1], result_normee$nouvelles_coordonnees[,3], 
     xlab="CP1", ylab="CP3", main="Premier vs Troisième axe principal NORMEE",
     type="n", xlim=c(min(c(result_normee$nouvelles_coordonnees[,1], result_dudi_normee$li[,1])),
                      max(c(result_normee$nouvelles_coordonnees[,1], result_dudi_normee$li[,1]))),
     ylim=c(min(c(result_normee$nouvelles_coordonnees[,3], result_dudi_normee$li[,3])),
            max(c(result_normee$nouvelles_coordonnees[,3], result_dudi_normee$li[,3]))))
points(result_normee$nouvelles_coordonnees[,1], result_normee$nouvelles_coordonnees[,3], col="blue", pch=16)
points(result_dudi_normee$li[,1], result_dudi_normee$li[,3], col="red", pch=17)
legend(x = 0.6, y = 0.5, legend=c("pca_from_scratch", "dudi.pca"),
       col=c("blue", "red"), pch=c(16, 17), cex=0.8, bty="n")

CCL PARTIE 2 - QUALITE DE ACP

Il semble que pour nos données k=3 semble être la meilleur solution (Regle de Kaiser - Guttman + les qualites des projections des individus proche de 1 pour k=3, bine plus que k=2). Mais cela dependra de notre but : representation graphique, reduction de dimension…

Partie 3 : Etude de la forme du nuage initiale et réduction de dimension

#8. Nuage isotrope

observation graphique de l’evolution de la matrice de covariance selon n

# Valeurs de n à tester
n_values <- c( 5, 10, 20, 50, 80, 100, 200, 500, 1000, 10000)

# Listes pour stocker les résultats pour chaque valeur de n pour les données centrées et centrées réduites
results_centrees <- list()
results_normees <- list()

for (n in n_values) {
  # Nuage isotrope
  X_iso <- rnorm(n)
  Y_iso <- rnorm(n)
  Z_iso <- rnorm(n)
  
  # Normalisation des vecteurs
  norm <- sqrt(X_iso^2 + Y_iso^2 + Z_iso^2)
  X_iso <- X_iso / norm
  Y_iso <- Y_iso / norm
  Z_iso <- Z_iso / norm
  
  # ACP sur le nuage isotrope - CENTRÉ
  acp_iso_centree <- acp_from_scratch(cbind(X_iso, Y_iso, Z_iso), reduire = FALSE)
  # ACP sur le nuage isotrope - CENTRÉ RÉDUIT
  acp_iso_normee <- acp_from_scratch(cbind(X_iso, Y_iso, Z_iso), reduire = TRUE)
  
  # Stockage des matrices de covariance
  results_centrees[[as.character(n)]] <- acp_iso_centree$matrice_covariance
  results_normees[[as.character(n)]] <- acp_iso_normee$matrice_covariance
}

# Préparation des données pour le tracé
cov_data_centrees <- data.frame()
cov_data_normees <- data.frame()

for (n in n_values) {
  cov_matrix_iso_centree <- results_centrees[[as.character(n)]]
  cov_matrix_iso_normee <- results_normees[[as.character(n)]]
  
  cov_data_centrees <- rbind(cov_data_centrees, cbind(
    n = n,
    variance_X = cov_matrix_iso_centree[1, 1],
    variance_Y = cov_matrix_iso_centree[2, 2],
    variance_Z = cov_matrix_iso_centree[3, 3],
    covariance_XY = cov_matrix_iso_centree[1, 2],
    covariance_XZ = cov_matrix_iso_centree[1, 3],
    covariance_YZ = cov_matrix_iso_centree[2, 3]
  ))
  
  cov_data_normees <- rbind(cov_data_normees, cbind(
    n = n,
    variance_X = cov_matrix_iso_normee[1, 1],
    variance_Y = cov_matrix_iso_normee[2, 2],
    variance_Z = cov_matrix_iso_normee[3, 3],
    covariance_XY = cov_matrix_iso_normee[1, 2],
    covariance_XZ = cov_matrix_iso_normee[1, 3],
    covariance_YZ = cov_matrix_iso_normee[2, 3]
  ))
}

library(gridExtra)

# Création des graphiques ggplot pour les données CENTRÉES et NORMÉES
plot_centrees <- ggplot(cov_data_centrees, aes(x = n)) +
  geom_line(aes(y = variance_X, color = "Variance X")) +
  geom_line(aes(y = variance_Y, color = "Variance Y")) +
  geom_line(aes(y = variance_Z, color = "Variance Z")) +
  geom_line(aes(y = covariance_XY, color = "Covariance XY")) +
  geom_line(aes(y = covariance_XZ, color = "Covariance XZ")) +
  geom_line(aes(y = covariance_YZ, color = "Covariance YZ")) +
  geom_hline(yintercept = 1/3, linetype = "dashed", color = "blue") +
  labs(title = "Évolution COV ISO + CENTREE", y = "Valeur", color = "Légende") +
  scale_x_log10() +
  theme_minimal()

plot_normees <- ggplot(cov_data_normees, aes(x = n)) +
  geom_line(aes(y = variance_X, color = "Variance X")) +
  geom_line(aes(y = variance_Y, color = "Variance Y")) +
  geom_line(aes(y = variance_Z, color = "Variance Z")) +
  geom_line(aes(y = covariance_XY, color = "Covariance XY")) +
  geom_line(aes(y = covariance_XZ, color = "Covariance XZ")) +
  geom_line(aes(y = covariance_YZ, color = "Covariance YZ")) +
  labs(title = "Évolution COV ISO + NORMÉE", y = "Valeur", color = "Légende") +
  scale_x_log10() +
  theme_minimal()

# Utilisation de grid.arrange pour tracer les graphiques côte à côte
grid.arrange(plot_centrees, plot_normees, ncol = 2)

CAS NORMEE\ • On retrouve les résultats connus : \[ \underset{i \neq j}{\text{COV}}(X_i, X_j) = 0 \text{ et } Var(U_i)=0\] \

CAS CENTREE\ • \[\underset{i \neq j}{\text{COV}}(X_i, X_j)\] tend rapidement vers 0 avec n. Evident car ce sont 3 vecteurs gaussiens indépendants. \

• La variance des 3 VG indépependant tend vers 1/3; Démonstration: \

Considérons un vecteur gaussien tridimensionnel \(\vec{X} = (X_1, X_2, X_3)\) où \(X_1, X_2,\) et \(X_3\) sont des variables aléatoires indépendantes issues d’une distribution normale \(\mathcal{N}(0, 1)\). (CENTRAGE : soustrait que 0) La norme de \(\vec{X}\) est donnée par :

\[ \|\vec{X}\| = \sqrt{X_1^2 + X_2^2 + X_3^2} \]

La norme au carré, \(\|\vec{X}\|^2\), suit une distribution chi-carrée avec 3 degrés de liberté puisqu’elle est la somme des carrés de trois variables normales standard indépendantes.

Lorsque nous normalisons le vecteur \(\vec{X}\), nous obtenons un vecteur unitaire \(\vec{U}\) :

\[ \vec{U} = \frac{\vec{X}}{\|\vec{X}\|} \]

Les vecteurs \(\vec{U}\) sont distribués sur la surface de la sphère unité. La variance dans chaque direction pour une grande collection de tels vecteurs est impactée par la normalisation. En raison de la contrainte que la somme des carré des composantes doit être égale à 1, la variance dans chaque direction est réduite. Mathématiquement :

\[ E[\|\vec{U}\|^2] = E[U_1^2 + U_2^2 + U_3^2] = 1 \]

Puisque \(\vec{U}\) est normalisé et ses composantes sont symétriques et interchangeables, la contribution attendue de chaque composante carrée à cette somme doit être égale. Ainsi, nous avons :

\[ E[U_1^2] = E[U_2^2] = E[U_3^2] = \frac{1}{3} \]

Et la variance de chaque composante est donc :

\[ Var(U_i) = E[U_i^2] - (E[U_i])^2 = \frac{1}{3} - 0^2 = \frac{1}{3} \quad \text{pour } i = 1, 2, 3 \]

Ce qui explique pourquoi, après la normalisation des vecteurs aléatoires gaussiens indépendants sur la sphère unité, la variance dans chaque direction est de \(\frac{1}{3}\).

Obseravtion de la cascade des Valeurs propres selon n

n_values <- c( 5, 10, 20, 50, 80, 100, 200, 500, 1000, 10000)
resultats_centres <- list()
resultats_normes <- list()

# Boucle pour chaque n pour réaliser l'ACP sur les données centrées et normalisées
for (n in n_values) {
  # Nuage isotrope
  X_iso <- rnorm(n)
  Y_iso <- rnorm(n)
  Z_iso <- rnorm(n)
  
  # Normalisation des vecteurs
  norm <- sqrt(X_iso^2 + Y_iso^2 + Z_iso^2)
  X_iso <- X_iso / norm
  Y_iso <- Y_iso / norm
  Z_iso <- Z_iso / norm
  
  
  # ACP sur le nuage isotrope - CENTRÉ
  acp_iso_centre <- acp_from_scratch(cbind(X_iso, Y_iso, Z_iso), reduire = FALSE)
  # ACP sur le nuage isotrope - CENTRÉ RÉDUIT
  acp_iso_norme <- acp_from_scratch(cbind(X_iso, Y_iso, Z_iso), reduire = TRUE)
  
  # Stocker les résultats
  resultats_centres[[as.character(n)]] <- acp_iso_centre$valeurs_propres
  resultats_normes[[as.character(n)]] <- acp_iso_norme$valeurs_propres
}

# Préparation des data frames pour le tracé
donnees_cov_centrees <- data.frame(n = numeric(), valeur_propre_1 = numeric(), valeur_propre_2 = numeric(), valeur_propre_3 = numeric())
donnees_cov_normees <- data.frame(n = numeric(), valeur_propre_1 = numeric(), valeur_propre_2 = numeric(), valeur_propre_3 = numeric())

# Agréger les résultats pour le tracé
for (n in n_values) {
  donnees_cov_centrees <- rbind(donnees_cov_centrees, data.frame(n = n, valeur_propre_1 = resultats_centres[[as.character(n)]][1], valeur_propre_2 = resultats_centres[[as.character(n)]][2], valeur_propre_3 = resultats_centres[[as.character(n)]][3]))
  donnees_cov_normees <- rbind(donnees_cov_normees, data.frame(n = n, valeur_propre_1 = resultats_normes[[as.character(n)]][1], valeur_propre_2 = resultats_normes[[as.character(n)]][2], valeur_propre_3 = resultats_normes[[as.character(n)]][3]))
}

library(gridExtra)

# Création des graphiques ggplot pour les données CENTRÉES et NORMALISÉES (CENTRÉES RÉDUITES)
plot_centrees <- ggplot(donnees_cov_centrees, aes(x = n)) +
  geom_line(aes(y = valeur_propre_1, color = "Valeur propre 1")) +
  geom_line(aes(y = valeur_propre_2, color = "Valeur propre 2")) +
  geom_line(aes(y = valeur_propre_3, color = "Valeur propre 3")) +
  geom_hline(yintercept = 1/3, linetype = "dashed", color = "blue") +
  labs(title = "Évolution valeurs propres CENTRÉES", y = "Valeur", color = "Valeurs propres \n ordre décroissant") +
  scale_x_log10() +
  theme_minimal()

plot_normees <- ggplot(donnees_cov_normees, aes(x = n)) +
  geom_line(aes(y = valeur_propre_1, color = "Valeur propre 1")) +
  geom_line(aes(y = valeur_propre_2, color = "Valeur propre 2")) +
  geom_line(aes(y = valeur_propre_3, color = "Valeur propre 3")) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "blue") +
  labs(title = "Évolution valeurs propres NORMEES", y = "Valeur", color = "Valeurs propres \n ordre décroissant") +
  scale_x_log10() +
  theme_minimal()

# Utilisation de grid.arrange pour tracer les graphiques côte à côte
grid.arrange(plot_centrees, plot_normees, ncol = 2)

CAS CENTREE: • La matrice de covariance tend vers \[\frac{1}3 I_3\] et donc les 3 VP sont bien \[\frac{1}3\].

CAS NORMEE: • La matrice de covariance tend vers \[I_3\] et donc les 3 VP sont bien \[1\].

Qualité des projections des individus Qik=2 selon la formule du cours selon n CENTREE puis NORMEE

library(pander)

# Valeurs de n à itérer
n_values <- c(5, 10, 20, 50, 80, 100, 200, 500, 1000, 10000)
k <- 2  # Nombre de composantes à considérer

# Initialisation d'un data frame pour stocker les résultats
results_df <- data.frame(
  n = integer(),
  moyenne_centree = numeric(),
  ecart_type_centree = numeric(),
  moyenne_normee = numeric(),
  ecart_type_normee = numeric()
)

for (n in n_values) {
  # Génération d'un nuage isotrope
  X_iso <- rnorm(n)
  Y_iso <- rnorm(n)
  Z_iso <- rnorm(n)
  
  # Normalisation des vecteurs pour avoir des longueurs unitaires
  norm <- sqrt(X_iso^2 + Y_iso^2 + Z_iso^2)
  X_iso <- X_iso / norm
  Y_iso <- Y_iso / norm
  Z_iso <- Z_iso / norm
  
  # ACP sur le nuage isotrope - Données CENTRÉES
  acp_iso_centre <- acp_from_scratch(cbind(X_iso, Y_iso, Z_iso), nombre_composants = NULL, reduire = FALSE)
  # ACP sur le nuage isotrope - Données NORMALISÉES
  acp_iso_norme <- acp_from_scratch(cbind(X_iso, Y_iso, Z_iso), nombre_composants = NULL, reduire = TRUE)
  
  # Calcul de la qualité des projections
  qualite_projections_centrees <- apply(acp_iso_centre$nouvelles_coordonnees, 1, function(Ci) sum(Ci[1:k]^2) / sum(Ci^2))
  qualite_projections_normees <- apply(acp_iso_norme$nouvelles_coordonnees, 1, function(Ci) sum(Ci[1:k]^2) / sum(Ci^2))
  
  #view(qualite_projections_centrees)
  
  # Calcul des moyennes et des écarts types
  moyenne_centree <- mean(qualite_projections_centrees)
  ecart_type_centree <- sd(qualite_projections_centrees)
  moyenne_normee <- mean(qualite_projections_normees)
  ecart_type_normee <- sd(qualite_projections_normees)
  
  # Ajout des résultats dans le data frame
  results_df <- rbind(results_df, data.frame(n, moyenne_centree, ecart_type_centree, moyenne_normee, ecart_type_normee))
}

# Affichage du data frame avec pander pour un joli formatage
pander(results_df)

n	moyenne_centree	ecart_type_centree	moyenne_normee	ecart_type_normee
5	0.91	0.11	0.81	0.28
10	0.89	0.15	0.91	0.11
20	0.82	0.27	0.81	0.25
50	0.72	0.22	0.7	0.32
80	0.74	0.26	0.73	0.27
100	0.7	0.28	0.69	0.31
200	0.7	0.3	0.68	0.3
500	0.68	0.29	0.69	0.29
1000	0.68	0.3	0.67	0.3
10000	0.67	0.3	0.67	0.3

En CENTREE ou NORMEE, on a les memes qualités de projection de moyenne 2/3, ce qui est normale car on est en cas ISOTROPE et d’ou le 2/3, chaque direction est equivalente. Toutes les variables ont des variances egales.

Heat map des contributions de l’individu i à l’inertie de l’axe factoriel sur les 2 premiers axes selon n

library(ggplot2)
library(gridExtra)
library(reshape2)

# Fonction pour calculer les contributions
calculate_contributions <- function(data, valeurs_propres) {
  sweep(data^2, 2, valeurs_propres, FUN="/") / nrow(data)
}

# Initialisation des listes pour stocker les heatmaps
heatmap_list <- list()

# Itération sur n = 10 et n = 50
for (n in c(10, 50)) {
  # Génération d'un nuage isotrope
  X_iso <- rnorm(n)
  Y_iso <- rnorm(n)
  Z_iso <- rnorm(n)
  
  # Normalisation des vecteurs pour avoir des longueurs unitaires
  norm <- sqrt(X_iso^2 + Y_iso^2 + Z_iso^2)
  X_iso <- X_iso / norm
  Y_iso <- Y_iso / norm
  Z_iso <- Z_iso / norm
  
  # ACP sur le nuage isotrope - Données CENTRÉES
  acp_iso_centre <- acp_from_scratch(cbind(X_iso, Y_iso, Z_iso), nombre_composants = 2, reduire = FALSE)
  
  # Calcul des contributions
  contributions <- calculate_contributions(acp_iso_centre$nouvelles_coordonnees, acp_iso_centre$valeurs_propres[c(1,2)])
  
  # Préparation des données pour les heatmaps
  data_for_heatmap <- melt(contributions)
  colnames(data_for_heatmap) <- c("Individu", "Axe", "Contribution")
  
  # Création de la heatmap
  p <- ggplot(data_for_heatmap, aes(x = Individu, y = Axe, fill = Contribution)) + 
    geom_tile() + 
    scale_fill_gradient(low = "blue", high = "red") +
    labs(x = "Individu", y = "Contrib. sur 2 composantes", title = paste("CENTREE n =", n)) +
    theme_minimal()
  
  # Stockage de la heatmap dans la liste
  heatmap_list[[as.character(n)]] <- p
}

# Affichage des heatmaps côte à côte
grid.arrange(heatmap_list[["10"]], heatmap_list[["50"]], ncol = 2)

Heat map des contributions de l’individu i à l’inertie de l’axe factoriel jsur les 2 premiers axes selon n

library(ggplot2)
library(pander)


# Fonction pour calculer les contributions
calculate_contributions <- function(data, valeurs_propres) {
  sweep(data^2, 2, valeurs_propres, FUN="/") / nrow(data)
}

# Génération d'un nuage isotrope et ACP pour n = 10
n <- 100
X_iso <- rnorm(n)
Y_iso <- rnorm(n)
Z_iso <- rnorm(n)
norm <- sqrt(X_iso^2 + Y_iso^2 + Z_iso^2)
X_iso <- X_iso / norm
Y_iso <- Y_iso / norm
Z_iso <- Z_iso / norm

acp_iso_centre_10 <- acp_from_scratch(cbind(X_iso, Y_iso, Z_iso),nombre_composants = 2, reduire = FALSE)
contributions_10 <- calculate_contributions(acp_iso_centre_10$nouvelles_coordonnees, acp_iso_centre_10$valeurs_propres[c(1,2)])

# Création et affichage de la heatmap pour n = 100
data_10 <- melt(contributions_10)
ggplot(data_10, aes(Var1, Var2, fill = value)) + 
  geom_tile() + 
  scale_fill_gradient(low = "blue", high = "red") +
  labs(x = "Individu", y = "Contribution sur les 3 composantes", title = "CENTREE n=100") +
  theme_minimal()

On ne remarque pas de diffenrence sur les contribution (en terme d’amplitude) des individus entre CP1 et CP2 car en cas isotrope.

Quelques statistisques sur les contributions selon n

library(pander)

# Fonction pour calculer les contributions
calculate_contributions <- function(data, valeurs_propres) {
  sweep(data^2, 2, valeurs_propres, FUN="/") / nrow(data)
}

# Valeurs de n à itérer
n_values <- c(5, 10, 20, 50, 80, 100, 200, 500, 1000, 10000)

# Initialisation d'un data frame pour stocker les résultats
results_df <- data.frame(
  n = integer(),
  moyenne_centree_axe1 = numeric(),
  ecart_type_centree_axe1 = numeric(),
  moyenne_centree_axe2 = numeric(),
  ecart_type_centree_axe2 = numeric(),
  moyenne_normee_axe1 = numeric(),
  ecart_type_normee_axe1 = numeric(),
  moyenne_normee_axe2 = numeric(),
  ecart_type_normee_axe2 = numeric()
)

for (n in n_values) {
  # Génération d'un nuage isotrope
  X_iso <- rnorm(n)
  Y_iso <- rnorm(n)
  Z_iso <- rnorm(n)
  
  # Normalisation des vecteurs pour avoir des longueurs unitaires
  norm <- sqrt(X_iso^2 + Y_iso^2 + Z_iso^2)
  X_iso <- X_iso / norm
  Y_iso <- Y_iso / norm
  Z_iso <- Z_iso / norm
  
  # ACP sur le nuage isotrope - Données CENTRÉES
  acp_iso_centre <- acp_from_scratch(cbind(X_iso, Y_iso, Z_iso), nombre_composants = 2, reduire = FALSE)
  contributions_centree <- calculate_contributions(acp_iso_centre$nouvelles_coordonnees, acp_iso_centre$valeurs_propres[c(1,2)])
  
  # ACP sur le nuage isotrope - Données NORMALISÉES
  acp_iso_normee <- acp_from_scratch(cbind(X_iso, Y_iso, Z_iso), nombre_composants = 2, reduire = TRUE)
  contributions_normee <- calculate_contributions(acp_iso_normee$nouvelles_coordonnees, acp_iso_normee$valeurs_propres[c(1,2)])

  # Calcul des moyennes et des écarts types pour chaque axe
  moyenne_centree_axe1 <- mean(contributions_centree[,1])
  ecart_type_centree_axe1 <- sd(contributions_centree[,1])
  moyenne_centree_axe2 <- mean(contributions_centree[,2])
  ecart_type_centree_axe2 <- sd(contributions_centree[,2])
  
  moyenne_normee_axe1 <- mean(contributions_normee[,1])
  ecart_type_normee_axe1 <- sd(contributions_normee[,1])
  moyenne_normee_axe2 <- mean(contributions_normee[,2])
  ecart_type_normee_axe2 <- sd(contributions_normee[,2])
  
  # Ajout des résultats dans le data frame
  results_df <- rbind(results_df, data.frame(n, moyenne_centree_axe1, ecart_type_centree_axe1, moyenne_centree_axe2, ecart_type_centree_axe2, moyenne_normee_axe1, ecart_type_normee_axe1, moyenne_normee_axe2, ecart_type_normee_axe2))
}

# Affichage du data frame avec pander pour un joli formatage
panderOptions('digits', 2)
pander(results_df)

Table continues below

n	moyenne_centree_axe1	ecart_type_centree_axe1	moyenne_centree_axe2
5	0.2	0.29	0.2
10	0.1	0.052	0.1
20	0.05	0.042	0.05
50	0.02	0.018	0.02
80	0.013	0.0093	0.013
100	0.01	0.0078	0.01
200	0.005	0.0041	0.005
500	0.002	0.0016	0.002
1000	0.001	0.00085	0.001
10000	1e-04	8.9e-05	1e-04

Table continues below

ecart_type_centree_axe2	moyenne_normee_axe1	ecart_type_normee_axe1
0.19	0.2	0.25
0.12	0.1	0.079
0.046	0.05	0.049
0.018	0.02	0.017
0.011	0.012	0.0096
0.01	0.01	0.0078
0.0047	0.005	0.0042
0.0018	0.002	0.0017
0.00091	0.001	0.00089
9e-05	1e-04	8.9e-05

moyenne_normee_axe2	ecart_type_normee_axe2
0.2	0.17
0.1	0.15
0.05	0.044
0.02	0.014
0.012	0.011
0.01	0.01
0.005	0.0047
0.002	0.0018
0.001	0.00088
1e-04	9e-05

On retrouve bien que : Soit \(S_j\) la somme des contributions pour un \(j\) fixé :

\[ S_j = \sum_{i=1}^{n} CONTR_{ij} \]

où \(CONTR_{ij}\) est la contribution du \(i\)-ème individu au \(j\)-ème élément. La moyenne de \(S_j\), notée \(\bar{S_j}\), est alors :

\[ \bar{S_j} = \frac{S_j}{n} = \frac{1}{n} \sum_{i=1}^{n} CONTR_{ij} \]

ce qui implique que :

\[ \bar{S_j} = \frac{1}{n} \]

9. Nuage non isotrope

observation graphique de l’evolution de la matrice de covariance selon n

library(ggplot2)
library(gridExtra)

# Valeurs de n à tester
n_values <- c(5, 10, 20, 50, 80, 100, 200, 500, 1000, 10000)

# Listes pour stocker les résultats pour chaque valeur de n pour les données centrées et centrées réduites
results_centrees <- list()
results_normees <- list()

for (n in n_values) {
  # Nuage non isotrope
  X_noniso <- sort(rnorm(n))
  Y_noniso <- rnorm(n)
  Z_noniso <- rnorm(n) + atan2(X_noniso, Y_noniso)
  
  # ACP sur le nuage non isotrope - CENTRÉ
  acp_noniso_centree <- acp_from_scratch(cbind(X_noniso, Y_noniso, Z_noniso), reduire = FALSE)
  # ACP sur le nuage non isotrope - CENTRÉ RÉDUIT
  acp_noniso_normee <- acp_from_scratch(cbind(X_noniso, Y_noniso, Z_noniso), reduire = TRUE)
  
  # Stockage des matrices de covariance
  results_centrees[[as.character(n)]] <- acp_noniso_centree$matrice_covariance
  results_normees[[as.character(n)]] <- acp_noniso_normee$matrice_covariance
}

# Préparation des données pour le tracé
cov_data_centrees <- data.frame()
cov_data_normees <- data.frame()

for (n in n_values) {
  cov_matrix_noniso_centree <- results_centrees[[as.character(n)]]
  cov_matrix_noniso_normee <- results_normees[[as.character(n)]]
  
  cov_data_centrees <- rbind(cov_data_centrees, cbind(
    n = n,
    variance_X = cov_matrix_noniso_centree[1, 1],
    variance_Y = cov_matrix_noniso_centree[2, 2],
    variance_Z = cov_matrix_noniso_centree[3, 3],
    covariance_XY = cov_matrix_noniso_centree[1, 2],
    covariance_XZ = cov_matrix_noniso_centree[1, 3],
    covariance_YZ = cov_matrix_noniso_centree[2, 3]
  ))
  
  cov_data_normees <- rbind(cov_data_normees, cbind(
    n = n,
    variance_X = cov_matrix_noniso_normee[1, 1],
    variance_Y = cov_matrix_noniso_normee[2, 2],
    variance_Z = cov_matrix_noniso_normee[3, 3],
    covariance_XY = cov_matrix_noniso_normee[1, 2],
    covariance_XZ = cov_matrix_noniso_normee[1, 3],
    covariance_YZ = cov_matrix_noniso_normee[2, 3]
  ))
}

# Création des graphiques ggplot pour les données CENTRÉES et NORMÉES pour le nuage non isotrope
plot_centrees_noniso <- ggplot(cov_data_centrees, aes(x = n)) +
  geom_line(aes(y = variance_X, color = "Variance X")) +
  geom_line(aes(y = variance_Y, color = "Variance Y")) +
  geom_line(aes(y = variance_Z, color = "Variance Z")) +
  geom_line(aes(y = covariance_XY, color = "Covariance XY")) +
  geom_line(aes(y = covariance_XZ, color = "Covariance XZ")) +
  geom_line(aes(y = covariance_YZ, color = "Covariance YZ")) +
  labs(title = "Évolution covariance selon n CENTRÉE (Non Isotrope)", y = "Valeur", color = "Légende") +
  scale_x_log10() +
  theme_minimal()

plot_normees_noniso <- ggplot(cov_data_normees, aes(x = n)) +
  geom_line(aes(y = variance_X, color = "Variance X")) +
  geom_line(aes(y = variance_Y, color = "Variance Y")) +
  geom_line(aes(y = variance_Z, color = "Variance Z")) +
  geom_line(aes(y = covariance_XY, color = "Covariance XY")) +
  geom_line(aes(y = covariance_XZ, color = "Covariance XZ")) +
  geom_line(aes(y = covariance_YZ, color = "Covariance YZ")) +
  labs(title = "Évolution covariance selon n NORMÉE (Non Isotrope)", y = "Valeur", color = "Légende") +
  scale_x_log10() +
  theme_minimal()

# Utilisation de grid.arrange pour tracer les graphiques côte à côte
grid.arrange(plot_centrees_noniso, plot_normees_noniso, ncol = 2)

CAS CENTREE : • Le nuage n’est pas isotrope, on n’a pas de résulats évident.

CAS NORMEE : • Les 3 variances sont bien égales à 1.

Obseravtion de la cascade des Valeurs propres selon n

n_values <- c(5, 10, 20, 50, 80, 100, 200, 500, 1000, 10000)
resultats_centres_noniso <- list()
resultats_normes_noniso <- list()

for (n in n_values) {
  # Nuage non isotrope
  X_noniso <- sort(rnorm(n))
  Y_noniso <- rnorm(n)
  Z_noniso <- rnorm(n) + atan2(X_noniso, Y_noniso)

  # ACP sur le nuage non isotrope - CENTRÉ
  acp_noniso_centre <- acp_from_scratch(cbind(X_noniso, Y_noniso, Z_noniso), reduire = FALSE)
  # ACP sur le nuage non isotrope - CENTRÉ RÉDUIT
  acp_noniso_norme <- acp_from_scratch(cbind(X_noniso, Y_noniso, Z_noniso), reduire = TRUE)

  # Stocker les résultats pour non isotrope
  resultats_centres_noniso[[as.character(n)]] <- acp_noniso_centre$valeurs_propres
  resultats_normes_noniso[[as.character(n)]] <- acp_noniso_norme$valeurs_propres
}

donnees_cov_centrees_noniso <- data.frame(n = numeric(), valeur_propre_1 = numeric(), valeur_propre_2 = numeric(), valeur_propre_3 = numeric())
donnees_cov_normees_noniso <- data.frame(n = numeric(), valeur_propre_1 = numeric(), valeur_propre_2 = numeric(), valeur_propre_3 = numeric())

# Agréger les résultats pour non isotrope
for (n in n_values) {
  donnees_cov_centrees_noniso <- rbind(donnees_cov_centrees_noniso, data.frame(n = n, valeur_propre_1 = resultats_centres_noniso[[as.character(n)]][1], valeur_propre_2 = resultats_centres_noniso[[as.character(n)]][2], valeur_propre_3 = resultats_centres_noniso[[as.character(n)]][3]))
  donnees_cov_normees_noniso <- rbind(donnees_cov_normees_noniso, data.frame(n = n, valeur_propre_1 = resultats_normes_noniso[[as.character(n)]][1], valeur_propre_2 = resultats_normes_noniso[[as.character(n)]][2], valeur_propre_3 = resultats_normes_noniso[[as.character(n)]][3]))
}

# Création des graphiques ggplot pour les données CENTRÉES et NORMALISÉES (CENTRÉES RÉDUITES) pour non isotrope
plot_centrees_noniso <- ggplot(donnees_cov_centrees_noniso, aes(x = n)) +
  geom_line(aes(y = valeur_propre_1, color = "Valeur propre 1")) +
  geom_line(aes(y = valeur_propre_2, color = "Valeur propre 2")) +
  geom_line(aes(y = valeur_propre_3, color = "Valeur propre 3")) +
  labs(title = "Évolution vp CENTRÉES (non isotrope)", y = "Valeur", color = "Valeurs propres \n ordre décroissant") +
  scale_x_log10() +
  theme_minimal()

plot_normees_noniso <- ggplot(donnees_cov_normees_noniso, aes(x = n)) +
  geom_line(aes(y = valeur_propre_1, color = "Valeur propre 1")) +
  geom_line(aes(y = valeur_propre_2, color = "Valeur propre 2")) +
  geom_line(aes(y = valeur_propre_3, color = "Valeur propre 3")) +
  labs(title = "Évolution vp NORMÉES (non isotrope)", y = "Valeur", color = "Valeurs propres \n ordre décroissant") +
  scale_x_log10() +
  theme_minimal()

# Utilisation de grid.arrange pour tracer les graphiques côte à côte pour non isotrope
grid.arrange(plot_centrees_noniso, plot_normees_noniso, ncol = 2)

CAS CENTREE: • On a la VP rouge plus grande que les autres : L’inertie captée par la 1ere VP provient de la variance de la variable Z car Z a une variance plus imporante que les 2 autres variables d’apres la matrice de covariance. • La VP verte tend vers 1 : c’est la variable Y qui est un VG indepndant de loi N(0,1).

CAS NORMEE: • Dans le cas normée, la. variance de Z est mise sur le pied d’égalité que celle de X et Y.

Si on ne veut pas regarder les heterogeneitées entre les variables (cela dépend de notre but), alors il vaudrait mieux NORMEE les variables.

Qualité des projections des individus Qik=2 selon la formule du cours selon n CENTREE puis NORMEE

library(pander)

# Valeurs de n à itérer
n_values <- c(5, 10, 20, 50, 80, 100, 200, 500, 1000, 10000)
k <- 1  # Nombre de composantes à considérer

# Initialisation d'un data frame pour stocker les résultats
results_df <- data.frame(
  n = integer(),
  moyenne_centree_noniso = numeric(),
  ecart_type_centree_noniso = numeric(),
  moyenne_normee_noniso = numeric(),
  ecart_type_normee_noniso = numeric()
)

for (n in n_values) {
  # Nuage non isotrope
  X_noniso <- sort(rnorm(n))
  Y_noniso <- rnorm(n)
  Z_noniso <- rnorm(n) + atan2(X_noniso, Y_noniso)

  # ACP sur le nuage non isotrope - Données CENTRÉES
  acp_noniso_centre <- acp_from_scratch(cbind(X_noniso, Y_noniso, Z_noniso), nombre_composants = 2, reduire = FALSE)
  # ACP sur le nuage non isotrope - Données NORMALISÉES
  acp_noniso_norme <- acp_from_scratch(cbind(X_noniso, Y_noniso, Z_noniso), nombre_composants = 2, reduire = TRUE)
  
  # Calcul de la qualité des projections pour non isotrope
  qualite_projections_centrees_noniso <- apply(acp_noniso_centre$nouvelles_coordonnees, 1, function(Ci) sum(Ci[1:k]^2) / sum(Ci^2))
  qualite_projections_normees_noniso <- apply(acp_noniso_norme$nouvelles_coordonnees, 1, function(Ci) sum(Ci[1:k]^2) / sum(Ci^2))
  
  # Calcul des moyennes et des écarts types pour non isotrope
  moyenne_centree_noniso <- mean(qualite_projections_centrees_noniso)
  ecart_type_centree_noniso <- sd(qualite_projections_centrees_noniso)
  moyenne_normee_noniso <- mean(qualite_projections_normees_noniso)
  ecart_type_normee_noniso <- sd(qualite_projections_normees_noniso)
  
  # Ajout des résultats dans le data frame pour non isotrope
  results_df <- rbind(results_df, data.frame(n, moyenne_centree_noniso, ecart_type_centree_noniso, moyenne_normee_noniso, ecart_type_normee_noniso))
}

# Renommage des colonnes pour la clarté
names(results_df) <- c("n", "Moyenne Centrée Non-Iso", "Écart-type Centrée Non-Iso", "Moyenne Normée Non-Iso", "Écart-type Normée Non-Iso")

# Affichage du data frame avec pander pour un joli formatage
pander(results_df)

Table continues below

n	Moyenne Centrée Non-Iso	Écart-type Centrée Non-Iso
5	0.81	0.29
10	0.79	0.27
20	0.75	0.29
50	0.73	0.32
80	0.72	0.3
100	0.75	0.27
200	0.74	0.29
500	0.74	0.31
1000	0.73	0.3
10000	0.74	0.3

Moyenne Normée Non-Iso	Écart-type Normée Non-Iso
0.71	0.31
0.72	0.2
0.67	0.33
0.63	0.31
0.61	0.32
0.58	0.31
0.62	0.31
0.64	0.33
0.63	0.32
0.63	0.32

# Initialisation d'un data frame pour stocker les résultats
results_df <- data.frame(
  n = integer(),
  moyenne_centree_axe1 = numeric(),
  ecart_type_centree_axe1 = numeric(),
  moyenne_centree_axe2 = numeric(),
  ecart_type_centree_axe2 = numeric(),
  moyenne_normee_axe1 = numeric(),
  ecart_type_normee_axe1 = numeric(),
  moyenne_normee_axe2 = numeric(),
  ecart_type_normee_axe2 = numeric()
)

for (n in n_values) {
  # Génération d'un nuage non isotrope
  # Nuage non isotrope
  X_noniso <- sort(rnorm(n))
  Y_noniso <- rnorm(n)
  Z_noniso <- rnorm(n) + atan2(X_noniso, Y_noniso)

  # ACP sur le nuage non isotrope - Données CENTRÉES
  acp_noniso_centre <- acp_from_scratch(cbind(X_noniso, Y_noniso, Z_noniso), nombre_composants = 2, reduire = FALSE)
  contributions_centree <- calculate_contributions(acp_noniso_centre$nouvelles_coordonnees, acp_noniso_centre$valeurs_propres[c(1,2)])
  
  # ACP sur le nuage non isotrope - Données NORMALISÉES
  acp_noniso_normee <- acp_from_scratch(cbind(X_noniso, Y_noniso, Z_noniso), nombre_composants = 2, reduire = TRUE)
  contributions_normee <- calculate_contributions(acp_noniso_normee$nouvelles_coordonnees, acp_noniso_normee$valeurs_propres[c(1,2)])

  # Calcul des moyennes et des écarts types pour chaque axe
  moyenne_centree_axe1 <- mean(contributions_centree[,1])
  ecart_type_centree_axe1 <- sd(contributions_centree[,1])
  moyenne_centree_axe2 <- mean(contributions_centree[,2])
  ecart_type_centree_axe2 <- sd(contributions_centree[,2])
  
  moyenne_normee_axe1 <- mean(contributions_normee[,1])
  ecart_type_normee_axe1 <- sd(contributions_normee[,1])
  moyenne_normee_axe2 <- mean(contributions_normee[,2])
  ecart_type_normee_axe2 <- sd(contributions_normee[,2])
  
  # Ajout des résultats dans le data frame
  results_df <- rbind(results_df, data.frame(n, moyenne_centree_axe1, ecart_type_centree_axe1, moyenne_centree_axe2, ecart_type_centree_axe2, moyenne_normee_axe1, ecart_type_normee_axe1, moyenne_normee_axe2, ecart_type_normee_axe2))
}

La projection des individus est meilleur dans le cas NONISO que dans le cas ISO. En effet dans le cas ISO les vecteurs sont comme sur une sphere et aucune direction n’est privilégié et donc la qualité de projection est moins bonne ou ici les vecteurs ont tendencent a être orienter selon une direction, donc une qualité plus grande en NONISO que en ISO.

Heat map des contributions de l’individu i à l’inertie de l’axe factoriel sur les 2 premiers axes selon n

library(ggplot2)
library(gridExtra)
library(reshape2)

# Fonction pour calculer les contributions
calculate_contributions <- function(data, valeurs_propres) {
  sweep(data^2, 2, valeurs_propres, FUN="/") / nrow(data)
}

# Initialisation des listes pour stocker les heatmaps
heatmap_list <- list()

# Itération sur n = 10 et n = 50
for (n in c(10, 50)) {
  # Génération d'un nuage non isotrope
  # Nuage non isotrope
  X_noniso <- sort(rnorm(n))
  Y_noniso <- rnorm(n)
  Z_noniso <- rnorm(n) + atan2(X_noniso, Y_noniso)

  
  # ACP sur le nuage non isotrope - Données CENTRÉES
  acp_noniso_centre <- acp_from_scratch(cbind(X_noniso, Y_noniso, Z_noniso), nombre_composants = 2, reduire = FALSE)
  
  # Calcul des contributions
  contributions_noniso <- calculate_contributions(acp_noniso_centre$nouvelles_coordonnees, acp_noniso_centre$valeurs_propres[c(1,2)])
  
  # Préparation des données pour les heatmaps
  data_for_heatmap_noniso <- melt(contributions_noniso)
  colnames(data_for_heatmap_noniso) <- c("Individu", "Axe", "Contribution")
  
  # Création de la heatmap
  p_noniso <- ggplot(data_for_heatmap_noniso, aes(x = Individu, y = Axe, fill = Contribution)) + 
    geom_tile() + 
    scale_fill_gradient(low = "blue", high = "red") +
    labs(x = "Individu", y = "Contrib. sur 2ere composantes", title = paste("Heatmap non isotrope n =", n)) +
    theme_minimal()
  
  # Stockage de la heatmap dans la liste
  heatmap_list[[as.character(n)]] <- p_noniso
}

# Affichage des heatmaps côte à côte
grid.arrange(heatmap_list[["10"]], heatmap_list[["50"]], ncol = 2)

Heat map des contributions de l’individu i à l’inertie de l’axe factoriel sur les 2 premiers axes selon n

library(ggplot2)
library(pander)


# Fonction pour calculer les contributions
calculate_contributions <- function(data, valeurs_propres) {
  sweep(data^2, 2, valeurs_propres, FUN="/") / nrow(data)
}
n <- 100
# Nuage non isotrope
X_noniso <- sort(rnorm(n))
Y_noniso <- rnorm(n)
Z_noniso <- rnorm(n) + atan2(X_noniso, Y_noniso)


acp_noniso_centre_10 <- acp_from_scratch(cbind(X_noniso, Y_noniso, Z_noniso),nombre_composants = 2, reduire = FALSE)
contributions_10 <- calculate_contributions(acp_noniso_centre_10$nouvelles_coordonnees, acp_noniso_centre_10$valeurs_propres[c(1,2)])

# Création et affichage de la heatmap pour n = 100
data_10 <- melt(contributions_10)
ggplot(data_10, aes(Var1, Var2, fill = value)) + 
  geom_tile() + 
  scale_fill_gradient(low = "blue", high = "red") +
  labs(x = "Individu", y = "Contribution sur les 2 composantes", title = "Heatmap des contributions pour n=100") +
  theme_minimal()

Quelques statistisques sur les contributions selon n

library(pander)

# Fonction pour calculer les contributions
calculate_contributions <- function(data, valeurs_propres) {
  sweep(data^2, 2, valeurs_propres, FUN="/") / nrow(data)
}

# Valeurs de n à itérer
n_values <- c(5, 10, 20, 50, 80, 100, 200, 500, 1000, 10000)

# Initialisation d'un data frame pour stocker les résultats
results_df <- data.frame(
  n = integer(),
  moyenne_centree_axe1 = numeric(),
  ecart_type_centree_axe1 = numeric(),
  moyenne_centree_axe2 = numeric(),
  ecart_type_centree_axe2 = numeric(),
  moyenne_normee_axe1 = numeric(),
  ecart_type_normee_axe1 = numeric(),
  moyenne_normee_axe2 = numeric(),
  ecart_type_normee_axe2 = numeric()
)

for (n in n_values) {
  # Génération d'un nuage non isotrope
  # Nuage non isotrope
  X_noniso <- sort(rnorm(n))
  Y_noniso <- rnorm(n)
  Z_noniso <- rnorm(n) + atan2(X_noniso, Y_noniso)

  # ACP sur le nuage non isotrope - Données CENTRÉES
  acp_noniso_centre <- acp_from_scratch(cbind(X_noniso, Y_noniso, Z_noniso), nombre_composants = 2, reduire = FALSE)
  contributions_centree <- calculate_contributions(acp_noniso_centre$nouvelles_coordonnees, acp_noniso_centre$valeurs_propres[c(1,2)])
  
  # ACP sur le nuage non isotrope - Données NORMALISÉES
  acp_noniso_normee <- acp_from_scratch(cbind(X_noniso, Y_noniso, Z_noniso), nombre_composants = 2, reduire = TRUE)
  contributions_normee <- calculate_contributions(acp_noniso_normee$nouvelles_coordonnees, acp_noniso_normee$valeurs_propres[c(1,2)])

  # Calcul des moyennes et des écarts types pour chaque axe
  moyenne_centree_axe1 <- mean(contributions_centree[,1])
  ecart_type_centree_axe1 <- sd(contributions_centree[,1])
  moyenne_centree_axe2 <- mean(contributions_centree[,2])
  ecart_type_centree_axe2 <- sd(contributions_centree[,2])
  
  moyenne_normee_axe1 <- mean(contributions_normee[,1])
  ecart_type_normee_axe1 <- sd(contributions_normee[,1])
  moyenne_normee_axe2 <- mean(contributions_normee[,2])
  ecart_type_normee_axe2 <- sd(contributions_normee[,2])
  
  # Ajout des résultats dans le data frame
  results_df <- rbind(results_df, data.frame(n, moyenne_centree_axe1, ecart_type_centree_axe1, moyenne_centree_axe2, ecart_type_centree_axe2, moyenne_normee_axe1, ecart_type_normee_axe1, moyenne_normee_axe2, ecart_type_normee_axe2))
}

# Affichage du data frame avec pander pour un joli formatage
panderOptions('digits', 2)
pander(results_df)

Table continues below

n	moyenne_centree_axe1	ecart_type_centree_axe1	moyenne_centree_axe2
5	0.2	0.14	0.2
10	0.1	0.1	0.1
20	0.05	0.057	0.05
50	0.02	0.018	0.02
80	0.012	0.013	0.013
100	0.01	0.0089	0.01
200	0.005	0.0051	0.005
500	0.002	0.0019	0.002
1000	0.001	0.0011	0.001
10000	1e-04	1e-04	1e-04

Table continues below

ecart_type_centree_axe2	moyenne_normee_axe1	ecart_type_normee_axe1
0.23	0.2	0.1
0.14	0.1	0.094
0.055	0.05	0.052
0.024	0.02	0.02
0.017	0.013	0.013
0.014	0.01	0.0098
0.0073	0.005	0.0046
0.0031	0.002	0.0018
0.0013	0.001	0.00095
0.00015	1e-04	9.6e-05

moyenne_normee_axe2	ecart_type_normee_axe2
0.2	0.15
0.1	0.094
0.05	0.061
0.02	0.023
0.012	0.017
0.01	0.014
0.005	0.0073
0.002	0.003
0.001	0.0013
1e-04	0.00015

Les resulats sont similaires aux cas ISO, donc même commentaire.

Conclusion : Il est essentiel de NORMALISER les données pour éviter que certaines variables avec des variances élevées dominent l’ACP au détriment d’autres variables présentant des variances plus faibles. En effet, si les données originales possèdent des variances non similaires, la normalisation permet de s’assurer que chaque variable contribue équitablement à l’analyse. Ainsi, pour maximiser l’efficacité de l’ACP dans le cas de variances non similaires entre les données originales, la normalisation est une étape cruciale.

10. Points extrêmes

Code pour visualiser le nuage de points isotrope 3D avec points extremaux et nuage compacté en 3 couleurs

library(rgl)

threshold_z <- 3.5

# Non isotrope
set.seed(123) # Pour la reproductibilité
n <- 100
x <- rnorm(n)
y <- 2 * x + rnorm(n)
z <- 5 * x * rnorm(n) 
data <- data.frame(x, y, z)

# Ajouter quelques points extrêmes
extreme_points <- data.frame(x = c(10, 12, -18, -7, 7, 15), 
                             y = c(20, 24, -20, -24, 9, -30),
                             z = rnorm(6))  

data_with_extremes <- rbind(data, extreme_points)

# Supprimer certains points extrêmes sur l'axe z pour les points compacts
data_without_extremes <- data[abs(data$z) <= threshold_z, ]

# Utiliser plot3d pour afficher les points en 3D
plot3d(data_without_extremes$x, data_without_extremes$y, data_without_extremes$z, col="green", size=3)  # Points compacts
points3d(data$x, data$y, data$z, col="blue", size=3)           # Points de données initiaux
points3d(extreme_points$x, extreme_points$y, extreme_points$z, col="red", size=3)  # Points extrêmes


# Ajouter une légende si nécessaire
legend3d("topright", legend = c("Data", "Extremes", "Compacts"), col = c("blue", "red", "green"), pch = 16)

# Fenêtre interactive
rgl.open()

## Warning: 'rgl.open' is deprecated.
## Use 'open3d' instead.
## See help("Deprecated")

library(plotly)

threshold_z <- 3.5

# Non isotrope
set.seed(123) # Pour la reproductibilité
n <- 100
x <- rnorm(n)
y <- 2 * x + rnorm(n)
z <- 5 * x * rnorm(n) 
data <- data.frame(x, y, z)

# Ajouter quelques points extrêmes
extreme_points <- data.frame(x = c(10, 12, -18, -7, 7, 15), 
                             y = c(20, 24, -20, -24, 9, -30),
                             z = rnorm(6))  

data_with_extremes <- rbind(data, extreme_points)

# Supprimer certains points extrêmes sur l'axe z pour les points compacts
data_without_extremes <- data[abs(data$z) <= threshold_z, ]

# Utiliser Plotly pour afficher les points en 3D
p <- plot_ly() %>%
  add_trace(data = data, x = ~x, y = ~y, z = ~z, 
            type = 'scatter3d', mode = 'markers',
            marker = list(color = 'blue', size = 3)) %>%
  add_trace(data = extreme_points, x = ~x, y = ~y, z = ~z, 
            type = 'scatter3d', mode = 'markers',
            marker = list(color = 'red', size = 3)) %>%
  add_trace(data = data_without_extremes, x = ~x, y = ~y, z = ~z, 
            type = 'scatter3d', mode = 'markers',
            marker = list(color = 'green', size = 3)) %>%
  layout(title = "3D Plot with Extremes and Compacts")

# Afficher le graphique
p

Visualisation du nuage 3D projeté dans Z=0

library(ggplot2)
library(gridExtra)

# Tracer les points avant l'ACP pour les données avec et sans points extrêmes
plot_data <- function(data, title, color) {
  ggplot(data, aes(x = x, y = y)) +
    geom_point(color = color) +
    labs(title = title, x = "Axe X", y = "Axe Y") +
    theme_minimal()
}

# Assurez-vous d'avoir les bons dataframes ici:
# data_with_extremes et data_without_extremes

plot_data(data_with_extremes, "Données avec points extrêmes projetées sur Z=0", "red")

plot_data(data_without_extremes, "Données compactées sur Z=0", "blue")

Etude des ACP CENTREE puis NORMEE sur les différents nuages

library(ggplot2)
library(gridExtra)

set.seed(123) # Pour la reproductibilité
n <- 100
x <- rnorm(n)
y <- 2 * x + rnorm(n)
z <- 5 * x * rnorm(n) 
data <- data.frame(x, y, z)

# Ajouter quelques points extrêmes
extreme_points <- data.frame(x = c(10, 12, -18, -7, 7, 15), 
                             y = c(20, 24, -20, -24, 9, -30),
                             z = rnorm(6))  
data_with_extremes <- rbind(data, extreme_points)

# Supprimer certains points extrêmes sur l'axe z pour les points compacts
threshold_z <- 3.5
data_without_extremes <- data[abs(data$z) <= threshold_z, ]

# Fonction pour tracer les résultats de l'ACP
plot_acp <- function(acp_results, title) {
  df <- as.data.frame(acp_results$nouvelles_coordonnees)
  ggplot(df, aes(x = V1, y = V2)) +
    geom_point() +
    labs(title = title, x = "PC1", y = "PC2") +
    theme_minimal()
}

# Perform PCA and plot for each dataset
# ACP centree
acp <- acp_from_scratch(as.matrix(data), nombre_composants = 2, reduire = FALSE)
acp_with_extremes <- acp_from_scratch(as.matrix(data_with_extremes), nombre_composants = 2, reduire = FALSE)
acp_without_extremes <- acp_from_scratch(as.matrix(data_without_extremes), nombre_composants = 2, reduire = FALSE)

plot1 <- plot_acp(acp, "iso CENT")
plot2 <- plot_acp(acp_with_extremes, "extrêmes CENT")
plot3 <- plot_acp(acp_without_extremes, "compacté CENT")

# ACP normee
acp_norm <- acp_from_scratch(as.matrix(data), nombre_composants = 2, reduire = TRUE)
acp_with_extremes_norm <- acp_from_scratch(as.matrix(data_with_extremes), nombre_composants = 2, reduire = TRUE)
acp_without_extremes_norm <- acp_from_scratch(as.matrix(data_without_extremes), nombre_composants = 2, reduire = TRUE)

plot4 <- plot_acp(acp_norm, "iso NORM")
plot5 <- plot_acp(acp_with_extremes_norm, "extrêmes NORM")
plot6 <- plot_acp(acp_without_extremes_norm, " compacté NORMEE")


plots <- list(plot1, plot4, plot2, plot5, plot3, plot6)
g <- gridExtra::marrangeGrob(plots, nrow = 2, ncol = 3)
grid::grid.newpage()
grid::grid.draw(g)

## [[1]]
## NULL

Observation:

Lorsque on execute ce code et observons les graphiques, on peut comparer les effets des points extrêmes sur la projection des individus en ACP centrée et normalisée. Avec des points extrêmes, les axes principaux peuvent largemement être influencés par ces points, alors que sans ces points, la structure intrinsèque des données pourrait être mieux représentée. Ils captent de la variance “pour rien”. En passant de CENTREE a NORMEE, on obtient une ACP avec une echelle différente pour les nouvelles coordonnées. Mais la structure reste trés similaire. Plus on prend des points exremes, plus les changements d’echelles seront grand et inversement avec le compactage, jusqu’a avoir quasi la meme echelle entre CENTREE et NORMEE si le compactage est important, car alors peut de variance restante a capter dans les donnees originales par rapport a la ponderation aux autres variances.

Conclusion:

L’ajout ou la suppression de points extrêmes peut avoir un impact significatif sur les résultats de l’ACP, en particulier sur la direction des premiers axes. Cela souligne l’importance de bien comprendre et prétraiter nos données avant de réaliser une ACP. Centrer ou normaliser les données peut aider à atténuer l’effet des points extrêmes, mais cela dépend également de la nature et de la distribution des données, comme on le verra en partie données environement.

Partie 4:Etude de la forme du nuage initiale sur la réduction de dimension dans les deux espaces

Mettre en oeuvre l’ACP Rn DUAL dans le cas ISO et cas NONISO

Fonction qui calcul matriciellement l’ACP DUAL FROM SCRATCH

acp_dual_from_scratch <- function(X, nombre_composants=NULL, reduire= FALSE) {
  
  n <- nrow(X)
  X <- as.matrix(X)
  
  # Data centree
  data <- centrage(X)
  
  # Data normee
  if (reduire) {
    data <- centrer_reduire(X)
  }
  
  
  
  # Calculer la matrice de Gram
  G <- data %*% t(data)
  
  # Calcul des valeurs propres et vecteurs propres pour G
  eig_result <- eigen(G)
  valeurs_propres <- eig_result$values
  vecteurs_propres_G <- eig_result$vectors
  
  # Sélectionner le nombre spécifié de composants principaux
  if (!is.null(nombre_composants)) {
    vecteurs_propres_G <- vecteurs_propres_G[, 1:nombre_composants, drop = FALSE]
  }
  
  # Projection des données sur les composants principaux
  nouvelles_coordonnees <- t(data) %*% vecteurs_propres_G
  
  return(list(matrice_gram=G,
              nouvelles_coordonnees=nouvelles_coordonnees, 
              valeurs_propres=valeurs_propres, 
              vecteurs_propres=vecteurs_propres_G))
  
}

1 et 2 : Mettre en oeuvre et comparer les resulats de nos 2 ACP classique et dual sur nos 2 cas ISO et NONISO

n <- 900
# Génération d'un nuage isotrope
X_iso <- rnorm(n)
Y_iso <- rnorm(n)
Z_iso <- rnorm(n)

# Normalisation des vecteurs pour avoir des longueurs unitaires
norm <- sqrt(X_iso^2 + Y_iso^2 + Z_iso^2)
X_iso <- X_iso / norm
Y_iso <- Y_iso / norm
Z_iso <- Z_iso / norm

data_isotrope <- data.frame(X_iso, Y_iso, Z_iso)


# Générer des données pour le cas non isotrope
X <- sort(rnorm(n))
Y <- rnorm(n)
Z <- rnorm(n) + atan2(X, Y)

data_non_isotrope <- data.frame(X, Y, Z)

# Appliquer l'ACP classique et duale sur les deux jeux de données.
res_acp_isotrope_classique <- acp_from_scratch(as.matrix(data_isotrope), reduire = TRUE)
res_acp_non_isotrope_classique <- acp_from_scratch(as.matrix(data_non_isotrope), reduire = TRUE)

res_acp_isotrope_dual <- acp_dual_from_scratch(as.matrix(data_isotrope), reduire = TRUE)
res_acp_non_isotrope_dual <- acp_dual_from_scratch(as.matrix(data_non_isotrope), reduire = TRUE)

plot(res_acp_isotrope_classique$nouvelles_coordonnees[,1], res_acp_isotrope_classique$nouvelles_coordonnees[,2], 
     main="ACP Classique - Isotrope", xlab="CP1", ylab="CP2", asp=1)

plot(res_acp_non_isotrope_classique$nouvelles_coordonnees[,1], res_acp_non_isotrope_classique$nouvelles_coordonnees[,2], 
     main="ACP Classique - Non Isotrope", xlab="CP1", ylab="CP2", asp=1)

plot(res_acp_isotrope_dual$nouvelles_coordonnees[,1], res_acp_isotrope_dual$nouvelles_coordonnees[,2], 
     main="ACP Duale - Isotrope", xlab="CP1", ylab="CP2", asp=1)

plot(res_acp_non_isotrope_dual$nouvelles_coordonnees[,1], res_acp_non_isotrope_dual$nouvelles_coordonnees[,2], 
     main="ACP Duale - Non Isotrope", xlab="CP1", ylab="CP2", asp=1)

library(plotly)

# Fonction pour générer un graphique Plotly 3D pour l'ACP
generate_3d_plot <- function(result, title) {
  ind_coords <- result$nouvelles_coordonnees
  
  p <- plot_ly() %>%
    add_trace(x = ~ind_coords[, 1], y = ~ind_coords[, 2], z = ~ind_coords[, 3],
              type = 'scatter3d', mode = 'markers',
              marker = list(size = 2, opacity = 0.8, color = 'darkblue')) %>%
    layout(title = title,
           scene = list(xaxis = list(title = 'PC1'),
                        yaxis = list(title = 'PC2'),
                        zaxis = list(title = 'PC3')))
  return(p)
}

# Générer les graphiques 3D pour chaque cas
p3d_acp_isotrope_classique <- generate_3d_plot(res_acp_isotrope_classique, "ACP Classique - Isotrope")
p3d_acp_non_isotrope_classique <- generate_3d_plot(res_acp_non_isotrope_classique, "ACP Classique - Non Isotrope")
p3d_acp_isotrope_dual <- generate_3d_plot(res_acp_isotrope_dual, "ACP Duale - Isotrope")
p3d_acp_non_isotrope_dual <- generate_3d_plot(res_acp_non_isotrope_dual, "ACP Duale - Non Isotrope")

# Afficher les graphiques
p3d_acp_isotrope_classique

p3d_acp_non_isotrope_classique

p3d_acp_isotrope_dual

p3d_acp_non_isotrope_dual

On obtient une quasi sphére du fait de l’isotropie des données choisies.

vérification avec PCA

# Appliquer l'ACP classique et duale sur les deux jeux de données.
res_acp_isotrope_classique_verif <- PCA(as.matrix(data_isotrope))

res_acp_non_isotrope_classique_verif <- PCA(as.matrix(data_non_isotrope))

Les resultats des 2 cas NON ISO des dex fonctions sont les mêmes; Mais dans le cas ISO la fonction PCA “évite” de placer sur les axes les points alors que notre fonction FROM SCRATCH non, il n’y a pas de raison mathématique de faire cela. C’est surement pour une raison “informatique” ou pratique.

Avec 900 individus et 3 variables, l’ACP dual est plus efficace en termes de calcul car elle évite de calculer une grande matrice de covariance et sa décomposition en VP, qui seraient de taille 900×900. À la place, on travaille avec une matrice de gram G de taille 3×3, ce qui est beaucoup plus gérable en termes de calcul.

3.

\[\begin{equation} \mathbf{v}_{\alpha} = \frac{1}{\sqrt{\lambda_{\alpha}}} \mathbf{Xu}_{\alpha} \end{equation}\]

\[\begin{equation} \mathbf{u}_{\alpha} = \frac{1}{\sqrt{\lambda_{\alpha}}} \mathbf{X}'\mathbf{v}_{\alpha} \end{equation}\]

4. Valider les résultats sur vos exemples

n<- 20
# Générer des données pour le cas non isotrope
X <- sort(rnorm(n))
Y <- rnorm(n)
Z <- rnorm(n) + atan2(X, Y)

data_non_isotrope <- data.frame(X, Y, Z)

# Data normee
data_non_isotrope <- centrer_reduire(as.matrix(data_non_isotrope))
  
# Calculer la matrice de Gram
G <- data_non_isotrope %*% t(data_non_isotrope)
M = t(data_non_isotrope) %*% data_non_isotrope
  
# Calcul des valeurs propres et vecteurs propres pour G
eigenG <- eigen(G)
lambdaG <- eigenG$values
vecteurs_propresG <- eigenG$vectors

eigenM <- eigen(M)
lambdaM <- eigenM$values
vecteurs_propresM <- eigenM$vectors

# On a bien les mêmes premiers valeurs propres
print(lambdaM)

## [1] 37.656686 14.114260  8.229053

print(lambdaG)

##  [1]  3.765669e+01  1.411426e+01  8.229053e+00  5.138273e-15  4.229563e-16
##  [6]  3.386689e-16  1.664530e-16  1.609501e-16  1.236492e-16  6.462368e-18
## [11] -2.240617e-17 -5.697641e-17 -2.078910e-16 -2.599617e-16 -4.704905e-16
## [16] -5.342331e-16 -9.797849e-16 -1.211067e-15 -1.626362e-15 -1.895188e-15

On a bien les mêmes premieres valeurs propres

On vérifie les formule de passsage mentionnées :

# Extrayons les valeurs et vecteurs propres pour la vérification

lambda <- lambdaM[1]
v_alphaM <- eigen(G)$vectors[,1]
u_alphaG <- eigen(M)$vectors[,1]

# data_non_iso (centré) matrice nxp
# G matrice nxn

v_alpha_calcule <- 1/sqrt(lambda) * data_non_isotrope %*% u_alphaG
u_alpha_calcule <- 1/sqrt(lambda) * t(data_non_isotrope) %*% v_alphaM

# Créer un DataFrame
df_v <- data.frame(v_alpha_calcule, v_alphaM)
df_u <- data.frame(u_alpha_calcule, u_alphaG)


# Imprimer le DataFrame
pander(df_v)

v_alpha_calcule	v_alphaM
-0.24	-0.24
-0.34	-0.34
-0.16	-0.16
-0.31	-0.31
-0.11	-0.11
-0.12	-0.12
-0.2	-0.2
-0.15	-0.15
-0.11	-0.11
-0.23	-0.23
0.24	0.24
0.36	0.36
0.024	0.024
0.01	0.01
0.18	0.18
0.29	0.29
0.41	0.41
0.075	0.075
0.1	0.1
0.27	0.27

pander(df_u)

	u_alpha_calcule	u_alphaG
X	0.58	0.58
Y	-0.51	-0.51
Z	0.63	0.63

On voit donc que u et v sont bien égaux a nos u et v calculés par relation de passage.

5.

CONCLUSION intermediaire : ACP Classique :

Utilisation courante : Lorsque le nombre de variables est relativement petit par rapport au nombre d’observations. Avantage : Plus efficace computationnellement dans des situations où le nombre de variables est beaucoup moins important que le nombre d’observations. Application typique : analyses où les variables sont d’un intérêt principal, par exemple, réduire la dimensionnalité d’un ensemble de caractéristiques dans un modèle de machine learning comme vu ce matin en cours. ## ACP Duale : Utilisation courante : Lorsque le nombre d’observations est beaucoup moins important que le nombre de variables. Avantage : Plus efficace computationnellement dans des situations où les calculs sur un nombre moins important d’observations sont préférables (par exemple, des données à haute dimensionnalité comme les données génomiques). Application typique : Dans les cas où les observations elles-mêmes sont d’un intérêt principal, comme dans certaines pour les gènes.

Quelques résultats ACP Rn et Rp :

Fonction pour imprimer jusqu’à 5 premières valeurs propres et les 5 premières coordonnées des vecteurs propres

# Fonction pour imprimer jusqu'à 3 premières valeurs propres et les 3 premières coordonnées des vecteurs propres
print_first_10_eigen_info <- function(result, title) {
  len_val_propres <- min(length(result$valeurs_propres), 3)
  len_vect_propres <- min(ncol(result$vecteurs_propres), 3)
  
  cat("=== ", title, " ===\n")
  cat("Premières Valeurs propres:\n")
  print(result$valeurs_propres[1:len_val_propres])
  cat("3 premières coordonnées des premiers Vecteurs propres :\n")
  if(len_vect_propres > 0){
    print(head(result$vecteurs_propres, 3)[, 1:len_vect_propres])
  } else {
    cat("Aucun vecteur propre disponible.\n")
  }
  cat("\n")
}

# Imprimer les 5 premières valeurs propres et vecteurs propres pour chaque ACP
print_first_10_eigen_info(res_acp_isotrope_classique, "ACP Classique - Isotrope")

## ===  ACP Classique - Isotrope  ===
## Premières Valeurs propres:
## [1] 1.0452412 1.0005029 0.9542559
## 3 premières coordonnées des premiers Vecteurs propres :
##            [,1]       [,2]       [,3]
## [1,]  0.7089907 0.01372909  0.7050842
## [2,]  0.2252462 0.94303205 -0.2448565
## [3,] -0.6682787 0.33241851  0.6655085

print_first_10_eigen_info(res_acp_non_isotrope_classique, "ACP Classique - Non Isotrope")

## ===  ACP Classique - Non Isotrope  ===
## Premières Valeurs propres:
## [1] 1.6195385 0.9990138 0.3814478
## 3 premières coordonnées des premiers Vecteurs propres :
##             [,1]       [,2]       [,3]
## [1,]  0.70674767 0.01481237  0.7073106
## [2,] -0.04832455 0.99845644  0.0273766
## [3,]  0.70581334 0.05352882 -0.7063726

print_first_10_eigen_info(res_acp_isotrope_dual, "ACP Duale - Isotrope")

## ===  ACP Duale - Isotrope  ===
## Premières Valeurs propres:
## [1] 940.7171 900.4526 858.8303
## 3 premières coordonnées des premiers Vecteurs propres :
##             [,1]        [,2]        [,3]
## [1,]  0.03819541 -0.02755634  0.03563111
## [2,] -0.03098348  0.04259255  0.01994804
## [3,] -0.02591993 -0.01720188 -0.04915552

print_first_10_eigen_info(res_acp_non_isotrope_dual, "ACP Duale - Non Isotrope")

## ===  ACP Duale - Non Isotrope  ===
## Premières Valeurs propres:
## [1] 1457.5846  899.1124  343.3030
## 3 premières coordonnées des premiers Vecteurs propres :
##             [,1]        [,2]        [,3]
## [1,] -0.08473914  0.01707746 -0.06731584
## [2,] -0.05803718 -0.01117835 -0.10924922
## [3,] -0.08024318 -0.02717637 -0.06158008

6. ProposONS la projection de points supplémentaires (espace sur R𝑝 ) et de points colonnes (variables) dans ce nouvel espace R𝑛

Bleu : Représente les données originales avant toute modification ou extension.

Rouge : Représente les données après la première extension, projeté par la matrice des données BLEU

Vert : Représente les données qui sont obtenues en centrant et réduisant l’ensemble de données étendu.

# Générer deux vecteurs gaussiens et utiliser ces vecteurs pour créer de nouvelles variables ou pour étendre l'échantillon.
n <- 900
# Génération d'un nuage isotrope
X_iso <- rnorm(n)
Y_iso <- rnorm(n)
Z_iso <- rnorm(n)

# Normalisation des vecteurs pour avoir des longueurs unitaires
norm <- sqrt(X_iso^2 + Y_iso^2 + Z_iso^2)
X_iso <- X_iso / norm
Y_iso <- Y_iso / norm
Z_iso <- Z_iso / norm

data_isotrope <- data.frame(X_iso, Y_iso, Z_iso)
data_isotrope <- centrer_reduire(as.matrix(data_isotrope))

V1 <- rnorm(n)
V2 <- rnorm(n)
V3 <- rnorm(n)
norm2 <- sqrt(V1^2 + V2^2 + V3^2)

# Normalisation des vecteurs pour avoir des longueurs unitaires
V1 <- V1 / norm2
V2 <- V2 / norm2
V3 <- V3 / norm2

# Ajouter ces vecteurs comme nouvelles colonnes à nos données
data_isotrope_etendu_variable <- as.data.frame(data_isotrope)
data_isotrope_etendu_variable$V1 <- V1
data_isotrope_etendu_variable$V2 <- V2
data_isotrope_etendu_variable$V3 <- V3
data_isotrope_etendu_variable2 <- data_isotrope_etendu_variable
data_isotrope_etendu_variable2 <- centrer_reduire(as.matrix(data_isotrope_etendu_variable2))
data_isotrope_etendu_variable$X_iso <- NULL
data_isotrope_etendu_variable$Y_iso <- NULL
data_isotrope_etendu_variable$Z_iso <- NULL
#print(data_isotrope_etendu_variable)
#print(dim(data_isotrope_etendu_variable))
data_isotrope_etendu_variable <- centrer_reduire(as.matrix(data_isotrope_etendu_variable))


# data non isotrope

X <- sort(rnorm(n))
Y <- rnorm(n)
Z <- rnorm(n) + atan2(X, Y)

data_non_isotrope <- data.frame(X, Y, Z)
data_non_isotrope <- centrer_reduire(as.matrix(data_non_isotrope))

data_non_isotrope_etendu_variable <- as.data.frame(data_non_isotrope)
data_non_isotrope_etendu_variable$V1 <- X + Z
data_non_isotrope_etendu_variable$V2 <- Y + X
data_non_isotrope_etendu_variable$V3 <- Z - X - Y
data_non_isotrope_etendu_variable2 <- data_non_isotrope_etendu_variable
data_non_isotrope_etendu_variable2 <- centrer_reduire(as.matrix(data_non_isotrope_etendu_variable2))
data_non_isotrope_etendu_variable$X <- NULL
data_non_isotrope_etendu_variable$Y <- NULL
data_non_isotrope_etendu_variable$Z <- NULL
data_non_isotrope_etendu_variable <- centrer_reduire(as.matrix(data_non_isotrope_etendu_variable))
#print(dim(data_non_isotrope_etendu_variable))





# On trouve les vecteurs propres SANS LES NV DONNEES
M_iso = t(as.matrix(data_isotrope)) %*% as.matrix(data_isotrope)
#print(dim(G_iso))
eig_result_iso <- eigen(M_iso)
vecteurs_propres_M_pre_extension_isotrope <- eig_result_iso$vectors
#print(dim(vecteurs_propres_G_pre_extension_isotrope))
#print(dim(data_isotrope_etendu_variable))

M_noniso = t(as.matrix(data_non_isotrope)) %*% as.matrix(data_non_isotrope)
eig_result_noniso <- eigen(M_noniso)
vecteurs_propres_M_pre_extension_non_isotrope <- eig_result_noniso$vectors


M_noniso2 = t(as.matrix(data_non_isotrope_etendu_variable2)) %*% as.matrix(data_non_isotrope_etendu_variable2)
eig_result_noniso2 <- eigen(M_noniso2)
vecteurs_propres_M_pre_extension_non_isotrope2 <- eig_result_noniso2$vectors

M_iso2 = t(as.matrix(data_isotrope_etendu_variable2)) %*% as.matrix(data_isotrope_etendu_variable2)
eig_result_iso2 <- eigen(M_iso2)
vecteurs_propres_M_pre_extension_isotrope2 <- eig_result_iso2$vectors


#print(dim(data_isotrope_etendu_variable))
#print(dim(vecteurs_propres_M_pre_extension_isotrope))

#print(dim(data_non_isotrope_etendu_variable))
#print(dim(data_non_isotrope_etendu_variable))
#print(dim(vecteurs_propres_M_pre_extension_non_isotrope))

# Projection des données sur les composants principaux
nouvelles_coordonnees_iso2 <- as.matrix(data_isotrope_etendu_variable) %*% vecteurs_propres_M_pre_extension_isotrope
nouvelles_coordonnees_non_iso2 <- as.matrix(data_non_isotrope_etendu_variable) %*% vecteurs_propres_M_pre_extension_non_isotrope

nouvelles_coordonnees_iso3 <- as.matrix(data_isotrope_etendu_variable2) %*% vecteurs_propres_M_pre_extension_isotrope2
nouvelles_coordonnees_non_iso3 <- as.matrix(data_non_isotrope_etendu_variable2) %*% vecteurs_propres_M_pre_extension_non_isotrope2

nouvelles_coordonnees_iso1 <- as.matrix(data_isotrope) %*% vecteurs_propres_M_pre_extension_isotrope
nouvelles_coordonnees_non_iso1 <- as.matrix(data_non_isotrope) %*% vecteurs_propres_M_pre_extension_non_isotrope 

# Créer un data frame avec vos données pour l'ensemble isotrope et non isotrope
data_iso <- data.frame(
  X = c(nouvelles_coordonnees_iso1[, 1], nouvelles_coordonnees_iso2[, 1], nouvelles_coordonnees_iso3[, 1]),
  Y = c(nouvelles_coordonnees_iso1[, 2], nouvelles_coordonnees_iso2[, 2], nouvelles_coordonnees_iso3[, 2]),
  Cas = rep(c("Cas 1", "Cas 2", "Cas 3"), times = c(nrow(nouvelles_coordonnees_iso1), nrow(nouvelles_coordonnees_iso2), nrow(nouvelles_coordonnees_iso3)))
)

data_non_iso <- data.frame(
  X = c(nouvelles_coordonnees_non_iso1[, 1], nouvelles_coordonnees_non_iso2[, 1], nouvelles_coordonnees_non_iso3[, 1]),
  Y = c(nouvelles_coordonnees_non_iso1[, 2], nouvelles_coordonnees_non_iso2[, 2], nouvelles_coordonnees_non_iso3[, 2]),
  Cas = rep(c("Cas 1", "Cas 2", "Cas 3"), times = c(nrow(nouvelles_coordonnees_non_iso1), nrow(nouvelles_coordonnees_non_iso2), nrow(nouvelles_coordonnees_non_iso3)))
)

# Créer le graphique pour l'ensemble isotrope avec toutes les coordonnées
ggplot(data_iso, aes(x = X, y = Y, color = Cas, shape = Cas)) +
  geom_point(size = 1) +
  labs(
    x = "Composante Principale 1",
    y = "Composante Principale 2",
    title = "ACP Isotrope"
  ) +
  scale_color_manual(values = c("blue", "red", "green")) +
  scale_shape_manual(values = c(16, 17, 18)) +
  theme_minimal() +
  theme(legend.position = "topright") +
  guides(color = guide_legend(title = "Cas"))

# Créer le graphique pour l'ensemble non isotrope avec toutes les coordonnées
ggplot(data_non_iso, aes(x = X, y = Y, color = Cas, shape = Cas)) +
  geom_point(size = 1) +
  labs(
    x = "Composante Principale 1",
    y = "Composante Principale 2",
    title = "ACP Non Isotrope"
  ) +
  scale_color_manual(values = c("blue", "red", "green")) +
  scale_shape_manual(values = c(16, 17, 18)) +
  theme_minimal() +
  theme(legend.position = "topright") +
  guides(color = guide_legend(title = "Cas"))

Partie ENVIRONEMENT

# Importer les données de l'énoncé
data <- read_excel("TP4_covC1234_DS19_20.xlsx")

# Supprimer les lignes avec des valeurs manquantes et des 0
data_noNA <- na.omit(data)
numeric_columns <- sapply(data_noNA, is.numeric)
data_clean <- data_noNA[apply(data_noNA[, numeric_columns], 1, function(row) all(row != 0)), ]
data_clean$indice <- NULL
#data_clean <- cbind(data_clean, ID=c(1:nrow(data_clean)))

#Sauvegarder les données nettoyées dans un nouveau fichier Excel
#write_xlsx(data_clean, "data_clean_ACP.xlsx")

# Sélectionner uniquement les colonnes numériques
data_numeric <- data_clean[sapply(data_clean, is.numeric)]
data_numeric <- data.frame(data_numeric)

On represente les box plot des 14 variables avec en rouge les limites des valeurs aberrantes selon la regle des 1,5 Q1 Q3

# Fonction pour calculer les bornes des valeurs aberrantes  
calculate_bounds <- function(column) {
  Q1 <- quantile(column, 0.25)
  Q3 <- quantile(column, 0.75)
  IQR <- Q3 - Q1
  lower_bound <- Q1 - 1.5 * IQR
  upper_bound <- Q3 + 1.5 * IQR
  return(c(lower_bound, upper_bound))
}

# Initialiser une liste pour stocker les bornes pour chaque colonne
outlier_bounds_list <- list()

# Créer un boxplot pour chaque variable numérique avec des lignes pour les limites des valeurs aberrantes
plot_list <- list()

for (col_name in names(data_numeric)) {
# Ajouter des backticks pour les noms de colonnes qui commencent par un chiffre ou contiennent des caractères spéciaux (au cas par cas)
  fixed_col_name <- ifelse(grepl("^[0-9]|[^[:alnum:]_]", col_name), paste0("`", col_name, "`"), col_name)

  bounds <- calculate_bounds(data_numeric[[col_name]])
  
  # Stocker les bornes dans la liste
  outlier_bounds_list[[col_name]] <- bounds
  
  p <- ggplot(data_numeric, aes_string(y = fixed_col_name)) +
    geom_boxplot() +
    geom_hline(yintercept = bounds[1], colour = "red", linetype = "dashed") +
    geom_hline(yintercept = bounds[2], colour = "red", linetype = "dashed") +
    labs(title = col_name, y = col_name)
  plot_list[[col_name]] <- p
}

# Afficher les boxplots
do.call("grid.arrange", c(plot_list, ncol = 4))

Pré traitement des données : On va garder la méthode 1 !

Methode 1, on garde nos data avec juste les 0 et les NA supprimées.

On importe les données en enlevant les valeurs NA et les 0.

En effet, pour les 0, les capteurs renverraient donc des mesures de cocnentration de 0 en moyenne, alors que les concentrations pour d’autres mêmes localisations et/ou d’autres même périodes et/ou d’autres même type et/ou campagne donne une cocnentration de plusieurs centaines ou plusieurs miliers en unité approprié. Il semble donc que ce soit des problèmes qui viennent des capteurs, des conditions environementales et étalonnages etc.. On choisi de supprimer les lignes avec 0, et on aura tout de meme assez de points pour realiser l’ACP.

Methode 2, on suprime ce qu’on considère comme outliers.

• Le problème est que selon une colonne, on a des valeurs mesurées en hiver, été, loin de l’usine de composatge polluante, plus proche etc… Donc un outlier peut ne pas l’être et inversemment. • Dans cette methode 2, on va retirer les éléments considérés comme outliers selon la règle des 1,5 Q1Q3, mais cela va aboutir à certains problèmes que l’on soulignera dans la suite.

Mise en oeuvre méthode 2:

# Function to identify if a value is an outlier based on the bounds
is_outlier <- function(x, bounds) {
  return(x < bounds[1] | x > bounds[2])
}

# Assuming 'data_clean' is your original dataset with numeric and non-numeric columns

# Apply functions only to numeric columns and save the names of the numeric columns
numeric_columns <- names(select_if(data_clean, is.numeric))

# Flag outliers within the dataset
data_with_outliers_flagged <- data_clean %>%
  mutate(across(all_of(numeric_columns), ~if_else(is_outlier(.x, calculate_bounds(.x)), TRUE, FALSE), .names = "outlier_{.col}"))

# Count the total number of outliers for each row and create a column 'is_outlier_row' to indicate rows with any outlier
data_with_outliers_flagged <- data_with_outliers_flagged %>%
  rowwise() %>%
  mutate(is_outlier_row = any(c_across(starts_with("outlier_")), na.rm = TRUE)) %>%
  ungroup()

# Get indices of rows to remove (which have outliers)
rows_to_remove <- which(data_with_outliers_flagged$is_outlier_row)

# Print the indices of rows with outliers
print(rows_to_remove)

##  [1]  19  35  36  38  44  48  49  50  51  52  53  54  55  56  57  58  59  60  61
## [20]  62  63  64  65  66  72  84  86  87  88  95 101 104 107

# Remove rows with any outlier
data_noVA <- data_with_outliers_flagged %>%
  filter(!is_outlier_row) %>%
  select(-starts_with("outlier_"), -is_outlier_row)

data_numeric_noVA <- data_noVA[sapply(data_noVA, is.numeric)]
data_numeric_noVA <- data.frame(data_numeric_noVA)

# Print out the number of rows removed
num_rows_removed <- length(rows_to_remove)
print(paste("Number of rows removed:", num_rows_removed))

## [1] "Number of rows removed: 33"

On voit que les lignes supprimées sont assez regroupés (par SAISON, TYPE etc…). Ceci suggère que les outliers supprimées n’en sont surement pas et que en fait on retrouve que : tout les hiver ne se ressemblent pas, tout les été aussi. Il y a des mesures localisé tout proche de l’usine, d’autres éloigné etc.. Ceci impacte grandement les données au sein même d’une colonne.

# Créer un boxplot pour chaque variable numérique avec des lignes pour les limites des valeurs aberrantes
plot_list <- list()

for (col_name in names(data_numeric_noVA)) {
# Ajouter des backticks pour les noms de colonnes qui commencent par un chiffre ou contiennent des caractères spéciaux (au cas par cas)
  fixed_col_name <- ifelse(grepl("^[0-9]|[^[:alnum:]_]", col_name), paste0("`", col_name, "`"), col_name)

  bounds <- calculate_bounds(data_numeric[[col_name]])
  p <- ggplot(data_numeric_noVA, aes_string(y = fixed_col_name)) +
    geom_boxplot() +
    geom_hline(yintercept = bounds[1], colour = "red", linetype = "dashed") +
    geom_hline(yintercept = bounds[2], colour = "red", linetype = "dashed") +
    labs(title = col_name, y = col_name)
  plot_list[[col_name]] <- p
}

# Afficher les boxplots
do.call("grid.arrange", c(plot_list, ncol = 4))

library(corrplot)

## corrplot 0.92 loaded

library(dplyr)

# Définir la mise en page pour afficher deux graphiques côte à côte
par(mfrow = c(1, 2))

# Calculer la matrice de corrélation et créer un graphique de corrélation
cor(data_numeric) %>% 
  corrplot::corrplot()

# Calculer la matrice de corrélation et créer un second graphique de corrélation
# avec une méthode de coloration différente et une commande de regroupement hiérarchique
cor(data_numeric) %>% 
  corrplot::corrplot(method = "color", order = "hclust")

on voit que les variables correles sont 13 et 14 ane, et les 10 ane, BTm, X1M2PA, E et X. On peut penser que ce sont peut etre des mêmes machines (ou plus generalement mêmes causes) qui emettent donc ces composant. Les correlations sont toutes psoitives ! Ce qui est logique, un uisne en fonctionnement normale va émetre de tout les composant et non pas certain et pas d’autres en fonction du temps et en moyenne.

1) Le traitement statistique des données

Calcul de la variabilité pour l’ensemble des données numeriques

data_numeric <- data.frame(data_numeric)
library(readxl)
library(dplyr)
library(ggplot2)
library(pander)

# Calcul de la variabilité pour l'ensemble des données numeriques
variability_overall <- data_numeric %>%
  summarise(across(where(is.numeric), list(moyenne = mean, ecart_type = sd, mediane = median)))

pander(variability_overall)

Table continues below

B_moyenne	B_ecart_type	B_mediane	T_moyenne	T_ecart_type	T_mediane
54	32	52	192	113	168

Table continues below

E_moyenne	E_ecart_type	E_mediane	X_moyenne	X_ecart_type	X_mediane
222	465	150	826	1309	575

Table continues below

X9_ane_moyenne	X9_ane_ecart_type	X9_ane_mediane	X10_ane_moyenne
100	102	74	254

Table continues below

X10_ane_ecart_type	X10_ane_mediane	X13_ane_moyenne	X13_ane_ecart_type
518	154	835	1325

Table continues below

X13_ane_mediane	X14_ane_moyenne	X14_ane_ecart_type	X14_ane_mediane
215	1652	1984	1104

Table continues below

X1_M_2_PA_moyenne	X1_M_2_PA_ecart_type	X1_M_2_PA_mediane	BTM_moyenne
267	538	126	547

Table continues below

BTM_ecart_type	BTM_mediane	FormicAcid_moyenne	FormicAcid_ecart_type
2075	303	492	590

Table continues below

FormicAcid_mediane	aceticacid_moyenne	aceticacid_ecart_type
305	353	718

Table continues below

aceticacid_mediane	NonaDecanoicAc_moyenne	NonaDecanoicAc_ecart_type
158	793	808

Table continues below

NonaDecanoicAc_mediane	Tot_OcNoDecana_moyenne	Tot_OcNoDecana_ecart_type
474	630	825

Tot_OcNoDecana_mediane
324

On remarque que les ecart-types sont trés important et donc on comprend que selon les variables qualitatives (SAISON, Campagnes etc…) on va avoir des mesures trés différnentes. Parfois, on a des moyennes et des médiannes trés différentes, c qui traduit la présence de valeures extrêmes dans les données, comme on l’a deja vu.

Calcul de la variabilité pour chaque campagne

# Calcul de la variabilité pour chaque campagne
campagnes_variability <- data_clean %>%
  group_by(Campagne) %>%
  summarise(across(where(is.numeric), list(moyenne = mean, ecart_type = sd, mediane = median)))
pander(campagnes_variability)

Table continues below

Campagne	B_moyenne	B_ecart_type	B_mediane	T_moyenne	T_ecart_type
BF2	16	3.4	15	98	33
BF3	19	3.4	19	173	91
CA1	42	9.5	39	219	75
CA2	105	16	107	314	145
CA3	53	5.7	52	261	63
CA4	69	15	70	120	49

Table continues below

T_mediane	E_moyenne	E_ecart_type	E_mediane	X_moyenne	X_ecart_type
81	62	22	56	207	68
154	161	59	136	440	110
197	175	78	161	687	223
262	587	1046	306	2169	2745
255	182	83	175	669	252
113	145	90	130	640	359

Table continues below

X_mediane	9_ane_moyenne	9_ane_ecart_type	9_ane_mediane	10_ane_moyenne
190	21	5.5	22	35
410	53	35	51	107
634	80	24	73	133
1336	233	133	182	729
615	72	17	73	263
581	107	101	87	200

Table continues below

10_ane_ecart_type	10_ane_mediane	13_ane_moyenne	13_ane_ecart_type
13	36	22	13
98	79	60	29
38	122	147	72
1127	480	3606	560
64	245	220	90
85	198	558	324

Table continues below

13_ane_mediane	14_ane_moyenne	14_ane_ecart_type	14_ane_mediane
17	48	36	31
53	32	19	28
140	433	207	427
3560	5572	875	5331
215	1245	445	1355
512	1705	542	1687

Table continues below

1_M_2_PA_moyenne	1_M_2_PA_ecart_type	1_M_2_PA_mediane	BTM_moyenne
21	18	11	114
76	69	50	317
133	49	131	328
855	1091	584	1749
180	144	124	336
233	186	165	332

Table continues below

BTM_ecart_type	BTM_mediane	FormicAcid_moyenne	FormicAcid_ecart_type
57	87	31	4.8
177	270	56	9.7
96	293	397	317
4885	534	987	158
166	322	180	98
149	321	885	828

Table continues below

FormicAcid_mediane	aceticacid_moyenne	aceticacid_ecart_type
33	12	5
54	60	84
327	840	1563
964	498	173
192	237	96
580	415	699

Table continues below

aceticacid_mediane	NonaDecanoicAc_moyenne	NonaDecanoicAc_ecart_type
12	61	44
39	61	33
185	463	283
397	1278	482
196	244	124
157	1772	620

Table continues below

NonaDecanoicAc_mediane	Tot_OcNoDecana_moyenne	Tot_OcNoDecana_ecart_type
48	77	53
55	148	76
474	380	264
1290	2125	954
211	276	120
1952	526	254

Tot_OcNoDecana_mediane
76
145
337
2169
313
461

On voit que les ecart-types sont plus petits que précédemment, meme si restent elevés donc on peut penser que, au sein d’une même campagne, les variables qualitatives comme la SAISON jouent un role. Les moyennes sont aussi plus proches des medianes, il y a moins de valeurs extrêmes au sein des campagnes.

Calcul de la variabilité pour les saisons été VS hiver

# Calcul de la variabilité pour les saisons
saisons_variability <- data_clean %>%
  group_by(SAISON) %>%
  summarise(across(where(is.numeric), list(moyenne = mean, ecart_type = sd, mediane = median)))
pander(saisons_variability)

Table continues below

SAISON	B_moyenne	B_ecart_type	B_mediane	T_moyenne	T_ecart_type
hiver	32	17	25	187	91
été	83	24	78	198	137

Table continues below

T_mediane	E_moyenne	E_ecart_type	E_mediane	X_moyenne	X_ecart_type
170	146	79	135	499	261
151	324	693	204	1258	1897

Table continues below

X_mediane	9_ane_moyenne	9_ane_ecart_type	9_ane_mediane	10_ane_moyenne
469	56	32	55	134
831	158	130	120	413

Table continues below

10_ane_ecart_type	10_ane_mediane	13_ane_moyenne	13_ane_ecart_type
102	108	111	96
755	278	1790	1572

Table continues below

13_ane_mediane	14_ane_moyenne	14_ane_ecart_type	14_ane_mediane
71	427	547	106
651	3268	2038	1935

Table continues below

1_M_2_PA_moyenne	1_M_2_PA_ecart_type	1_M_2_PA_mediane	BTM_moyenne
102	101	73	275
484	763	340	905

Table continues below

BTM_ecart_type	BTM_mediane	FormicAcid_moyenne	FormicAcid_ecart_type
161	255	162	215
3138	387	926	644

Table continues below

FormicAcid_mediane	aceticacid_moyenne	aceticacid_ecart_type
58	280	821
840	449	548

Table continues below

aceticacid_mediane	NonaDecanoicAc_moyenne	NonaDecanoicAc_ecart_type
81	203	224
333	1572	613

Table continues below

NonaDecanoicAc_mediane	Tot_OcNoDecana_moyenne	Tot_OcNoDecana_ecart_type
110	218	187
1492	1173	1011

Tot_OcNoDecana_mediane
149
711

On voit que clairement, il y une difference entre l’été et l’hiver, en été les moyennes sont bien plus élévées, ce qui veut dire que les conditions en été sont favorables a une augmentation des concentrations : temperature, humidité, activité des usines plus importantes ….. On va observer dans la suite par l’ACP ceci.

Calcul de la variabilité pour les localisations

# Calcul de la variabilité pour les localisations
localisation_variability <- data_clean %>%
  group_by(Localisation) %>%
  summarise(across(where(is.numeric), list(moyenne = mean)))

pander(localisation_variability)

Table continues below

Localisation	B_moyenne	T_moyenne	E_moyenne	X_moyenne	9_ane_moyenne
P01	55	182	162	679	103
P02	63	282	194	743	143
P03	55	167	157	626	98
P04	89	188	219	991	171
P05	55	169	177	739	90
P06	45	130	110	435	59
P07	52	150	127	498	81
P08	38	145	155	545	56
P09	43	178	146	575	59
P10	59	413	2527	6762	381
P11	41	215	234	866	72
P12	57	176	156	601	84
P13	36	148	121	439	51
P14	49	147	155	576	68
P15	51	354	386	1372	101
P16	46	172	163	702	79
P17	40	243	174	701	69
P18	49	154	164	600	79
P19	67	183	141	621	96
P20	17	103	143	288	226
P21	50	160	122	496	91
P25	60	175	143	582	76
P26	71	104	107	512	57
P27	77	170	141	589	82
P28	73	238	214	929	94
P29	87	210	288	1256	157
P30	72	142	152	718	103
P32	88	250	261	1290	250
P33	98	328	475	2369	198

Table continues below

10_ane_moyenne	13_ane_moyenne	14_ane_moyenne	1_M_2_PA_moyenne
227	1081	1916	283
339	1250	2104	228
181	976	1538	175
269	2267	2584	369
218	840	1868	206
187	1001	1609	283
172	678	1419	126
89	121	397	90
121	231	798	217
2718	1845	2520	2754
211	247	920	155
205	875	1640	157
104	86	665	85
157	676	1455	163
250	816	1691	342
212	208	764	118
167	308	915	149
240	737	1506	261
199	901	1669	185
78	166	529	12
213	730	1476	146
234	345	1617	325
141	472	1449	101
295	1731	3026	204
180	1221	2691	365
398	1795	3217	865
133	1941	3759	266
461	1562	3351	507
536	2263	4176	524

Table continues below

BTM_moyenne	FormicAcid_moyenne	aceticacid_moyenne	NonaDecanoicAc_moyenne
331	428	342	819
443	876	946	706
262	467	216	1046
356	456	881	1101
386	949	244	703
326	631	718	773
236	297	162	706
220	220	74	550
351	275	119	414
11327	510	205	411
371	178	95	605
351	428	177	757
202	281	138	553
239	387	181	703
488	334	161	713
360	337	133	550
323	379	1480	687
490	347	892	924
345	435	164	882
171	430	146	227
277	376	136	820
321	370	255	1159
212	643	201	1722
266	1422	1251	1066
266	498	359	593
706	764	310	1291
407	892	298	1814
787	1235	384	2183
764	641	272	921

Tot_OcNoDecana_moyenne
359
530
1240
1784
369
594
666
202
226
810
242
803
296
477
676
303
281
132
787
436
667
474
278
1166
638
1492
1082
2185
981

on voit que les moyennes varient grandement pour les diffenrentes localisations. Donc on voit bien que selon que l’on soit proche ou loin des usines, en milieu rurale ou urbaion etc.. cela a une gande importance sur les mesures, et donc il semble que les produits toxiques que l’on mesure soient bien non uniformes geographiquement.

Calcul de la variabilité pour les TYPE

# Calcul de la variabilité pour les localisations
localisation_variability <- data_clean %>%
  group_by(TYPE) %>%
  summarise(across(where(is.numeric), list(moyenne = mean)))

pander(localisation_variability)

Table continues below

TYPE	B_moyenne	T_moyenne	E_moyenne	X_moyenne
compostage	59	226	176	708
rural	57	161	150	603
sourceindustrie	48	256	940	2637
urbain	52	211	215	851

Table continues below

9_ane_moyenne	10_ane_moyenne	13_ane_moyenne	14_ane_moyenne
121	277	1156	1999
84	187	918	1686
166	987	769	1372
105	226	689	1582

Table continues below

1_M_2_PA_moyenne	BTM_moyenne	FormicAcid_moyenne	aceticacid_moyenne
258	381	627	610
206	301	503	329
1063	4010	353	148
228	386	472	354

NonaDecanoicAc_moyenne	Tot_OcNoDecana_moyenne
769	435
863	783
413	421
773	529

On voit que les moyennes sont soient similaires (comme B pour les 4 TYPES), soient trés differentes (comme E 5x plus grande pour sourceindustrie que pour les 3 autres TYPES). On voit donc ici que certains composants capté viennent surement des usines.

Calcul de la variabilité pour avant et après ouverture du site.

# Marquer les périodes avant et après activité
data <- data_clean %>%
  mutate(Periode = case_when(
    str_detect(Campagne, "BF") ~ "avant",
    TRUE ~ "apres"
  ))

# Calcul de la variabilité pour avant et après ouverture du site
periode_variability <- data %>%
  group_by(Periode) %>%
  summarise(across(where(is.numeric), list(moyenne = mean, ecart_type = sd, mediane = median)))
pander(periode_variability)

Table continues below

Periode	B_moyenne	B_ecart_type	B_mediane	T_moyenne	T_ecart_type
apres	70	26	66	215	117
avant	18	3.9	17	138	79

Table continues below

T_mediane	E_moyenne	E_ecart_type	E_mediane	X_moyenne	X_ecart_type
202	267	546	175	1032	1510
127	115	67	114	330	149

Table continues below

X_mediane	9_ane_moyenne	9_ane_ecart_type	9_ane_mediane	10_ane_moyenne
715	126	109	88	330
329	38	30	30	73

Table continues below

10_ane_ecart_type	10_ane_mediane	13_ane_moyenne	13_ane_ecart_type
599	229	1164	1456
80	52	42	30

Table continues below

13_ane_mediane	14_ane_moyenne	14_ane_ecart_type	14_ane_mediane
454	2322	2011	1533
40	39	29	30

Table continues below

1_M_2_PA_moyenne	1_M_2_PA_ecart_type	1_M_2_PA_mediane	BTM_moyenne
357	618	203	682
50	58	36	222

Table continues below

BTM_ecart_type	BTM_mediane	FormicAcid_moyenne	FormicAcid_ecart_type
2459	347	678	612
168	189	44	15

Table continues below

FormicAcid_mediane	aceticacid_moyenne	aceticacid_ecart_type
556	484	819
45	37	65

Table continues below

aceticacid_mediane	NonaDecanoicAc_moyenne	NonaDecanoicAc_ecart_type
235	1097	780
24	61	38

Table continues below

NonaDecanoicAc_mediane	Tot_OcNoDecana_moyenne	Tot_OcNoDecana_ecart_type
1046	844	899
55	115	74

Tot_OcNoDecana_mediane
473
112

On voit clairement la differnence AVANT et APRES la mise en place du site sur les moyennes des mesures des composants; on retrouvera ce resultat dans notre analyse ACP.

library(dplyr)
library(Hmisc)

## 
## Attaching package: 'Hmisc'

## The following object is masked from 'package:plotly':
## 
##     subplot

## The following objects are masked from 'package:dplyr':
## 
##     src, summarize

## The following objects are masked from 'package:base':
## 
##     format.pval, units

# Split the data by season
data_by_season <- split(data_clean, data_clean$SAISON)

# Calculate the correlation matrix only for the winter season
corr_hiver <- lapply(data_by_season["hiver"], function(df) {
  # Ensure we only use the relevant columns
  relevant_df <- df[, c("B", "T", "E", "X", "9_ane", "10_ane", "13_ane", "14_ane", "1_M_2_PA", "BTM", "FormicAcid", "aceticacid", "NonaDecanoicAc", "Tot_OcNoDecana")]
  rcorr(as.matrix(relevant_df), type = "pearson")
})

# To print the correlation matrix for winter
print(corr_hiver)

## $hiver
##                   B    T    E    X 9_ane 10_ane 13_ane 14_ane 1_M_2_PA  BTM
## B              1.00 0.59 0.42 0.65  0.59   0.72   0.78   0.80     0.60 0.39
## T              0.59 1.00 0.82 0.86  0.79   0.80   0.62   0.57     0.60 0.70
## E              0.42 0.82 1.00 0.92  0.56   0.52   0.47   0.39     0.62 0.65
## X              0.65 0.86 0.92 1.00  0.66   0.63   0.66   0.58     0.74 0.68
## 9_ane          0.59 0.79 0.56 0.66  1.00   0.77   0.58   0.43     0.51 0.72
## 10_ane         0.72 0.80 0.52 0.63  0.77   1.00   0.77   0.75     0.62 0.66
## 13_ane         0.78 0.62 0.47 0.66  0.58   0.77   1.00   0.89     0.72 0.61
## 14_ane         0.80 0.57 0.39 0.58  0.43   0.75   0.89   1.00     0.68 0.38
## 1_M_2_PA       0.60 0.60 0.62 0.74  0.51   0.62   0.72   0.68     1.00 0.67
## BTM            0.39 0.70 0.65 0.68  0.72   0.66   0.61   0.38     0.67 1.00
## FormicAcid     0.41 0.24 0.14 0.34  0.30   0.23   0.43   0.23     0.26 0.24
## aceticacid     0.19 0.24 0.14 0.27  0.18   0.18   0.26   0.08     0.24 0.19
## NonaDecanoicAc 0.51 0.37 0.26 0.46  0.55   0.31   0.51   0.37     0.33 0.32
## Tot_OcNoDecana 0.47 0.41 0.33 0.45  0.42   0.32   0.37   0.31     0.23 0.28
##                FormicAcid aceticacid NonaDecanoicAc Tot_OcNoDecana
## B                    0.41       0.19           0.51           0.47
## T                    0.24       0.24           0.37           0.41
## E                    0.14       0.14           0.26           0.33
## X                    0.34       0.27           0.46           0.45
## 9_ane                0.30       0.18           0.55           0.42
## 10_ane               0.23       0.18           0.31           0.32
## 13_ane               0.43       0.26           0.51           0.37
## 14_ane               0.23       0.08           0.37           0.31
## 1_M_2_PA             0.26       0.24           0.33           0.23
## BTM                  0.24       0.19           0.32           0.28
## FormicAcid           1.00       0.75           0.40           0.22
## aceticacid           0.75       1.00           0.23          -0.01
## NonaDecanoicAc       0.40       0.23           1.00           0.69
## Tot_OcNoDecana       0.22      -0.01           0.69           1.00
## 
## n= 62 
## 
## 
## P
##                B      T      E      X      9_ane  10_ane 13_ane 14_ane 1_M_2_PA
## B                     0.0000 0.0006 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000  
## T              0.0000        0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000  
## E              0.0006 0.0000        0.0000 0.0000 0.0000 0.0001 0.0016 0.0000  
## X              0.0000 0.0000 0.0000        0.0000 0.0000 0.0000 0.0000 0.0000  
## 9_ane          0.0000 0.0000 0.0000 0.0000        0.0000 0.0000 0.0004 0.0000  
## 10_ane         0.0000 0.0000 0.0000 0.0000 0.0000        0.0000 0.0000 0.0000  
## 13_ane         0.0000 0.0000 0.0001 0.0000 0.0000 0.0000        0.0000 0.0000  
## 14_ane         0.0000 0.0000 0.0016 0.0000 0.0004 0.0000 0.0000        0.0000  
## 1_M_2_PA       0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000         
## BTM            0.0015 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0022 0.0000  
## FormicAcid     0.0008 0.0598 0.2692 0.0073 0.0175 0.0701 0.0004 0.0741 0.0414  
## aceticacid     0.1434 0.0642 0.2666 0.0333 0.1681 0.1550 0.0391 0.5366 0.0640  
## NonaDecanoicAc 0.0000 0.0027 0.0387 0.0001 0.0000 0.0137 0.0000 0.0030 0.0085  
## Tot_OcNoDecana 0.0001 0.0008 0.0083 0.0002 0.0006 0.0114 0.0030 0.0138 0.0760  
##                BTM    FormicAcid aceticacid NonaDecanoicAc Tot_OcNoDecana
## B              0.0015 0.0008     0.1434     0.0000         0.0001        
## T              0.0000 0.0598     0.0642     0.0027         0.0008        
## E              0.0000 0.2692     0.2666     0.0387         0.0083        
## X              0.0000 0.0073     0.0333     0.0001         0.0002        
## 9_ane          0.0000 0.0175     0.1681     0.0000         0.0006        
## 10_ane         0.0000 0.0701     0.1550     0.0137         0.0114        
## 13_ane         0.0000 0.0004     0.0391     0.0000         0.0030        
## 14_ane         0.0022 0.0741     0.5366     0.0030         0.0138        
## 1_M_2_PA       0.0000 0.0414     0.0640     0.0085         0.0760        
## BTM                   0.0593     0.1302     0.0112         0.0255        
## FormicAcid     0.0593            0.0000     0.0013         0.0792        
## aceticacid     0.1302 0.0000                0.0776         0.9412        
## NonaDecanoicAc 0.0112 0.0013     0.0776                    0.0000        
## Tot_OcNoDecana 0.0255 0.0792     0.9412     0.0000

Pour analyser la matrice des correlations, on cherche les valeurs proche de 1 ou -1: On peut noter E et X, T et X, 13_ane et 14_ane notamment. De plus, les correlations sont positives, ce qui est logique car si on a un composant emis, alors les autres composant sont aussi emis des usines en questions et donc la correlation est positive.

# Split the data by season
data_by_season <- split(data_clean, data_clean$SAISON)

# Calculate the correlation matrix only for the summer season
correlations_ete <- lapply(data_by_season["été"], function(df) {
  # Ensure we only use the relevant columns
  relevant_df <- df[, c("B", "T", "E", "X", "9_ane", "10_ane", "13_ane", "14_ane", "1_M_2_PA", "BTM", "FormicAcid", "aceticacid", "NonaDecanoicAc", "Tot_OcNoDecana")]
  rcorr(as.matrix(relevant_df), type = "pearson")
})

# To print the correlation matrix for summer
print(correlations_ete)

## $été
##                    B     T     E     X 9_ane 10_ane 13_ane 14_ane 1_M_2_PA
## B               1.00  0.63  0.19  0.30  0.17   0.24   0.76   0.73     0.30
## T               0.63  1.00  0.68  0.78  0.64   0.66   0.69   0.69     0.71
## E               0.19  0.68  1.00  0.98  0.71   0.98   0.31   0.27     0.95
## X               0.30  0.78  0.98  1.00  0.71   0.96   0.39   0.36     0.94
## 9_ane           0.17  0.64  0.71  0.71  1.00   0.70   0.45   0.43     0.68
## 10_ane          0.24  0.66  0.98  0.96  0.70   1.00   0.34   0.31     0.96
## 13_ane          0.76  0.69  0.31  0.39  0.45   0.34   1.00   0.96     0.39
## 14_ane          0.73  0.69  0.27  0.36  0.43   0.31   0.96   1.00     0.35
## 1_M_2_PA        0.30  0.71  0.95  0.94  0.68   0.96   0.39   0.35     1.00
## BTM             0.12  0.56  0.98  0.94  0.65   0.98   0.21   0.17     0.94
## FormicAcid      0.16  0.06  0.02  0.02  0.05   0.02   0.06   0.05     0.02
## aceticacid      0.16  0.04 -0.03 -0.03 -0.02  -0.02   0.11   0.05    -0.02
## NonaDecanoicAc -0.15 -0.34 -0.26 -0.28 -0.31  -0.21  -0.38  -0.31    -0.25
## Tot_OcNoDecana  0.65  0.59  0.17  0.24  0.47   0.20   0.72   0.67     0.27
##                  BTM FormicAcid aceticacid NonaDecanoicAc Tot_OcNoDecana
## B               0.12       0.16       0.16          -0.15           0.65
## T               0.56       0.06       0.04          -0.34           0.59
## E               0.98       0.02      -0.03          -0.26           0.17
## X               0.94       0.02      -0.03          -0.28           0.24
## 9_ane           0.65       0.05      -0.02          -0.31           0.47
## 10_ane          0.98       0.02      -0.02          -0.21           0.20
## 13_ane          0.21       0.06       0.11          -0.38           0.72
## 14_ane          0.17       0.05       0.05          -0.31           0.67
## 1_M_2_PA        0.94       0.02      -0.02          -0.25           0.27
## BTM             1.00       0.01      -0.03          -0.20           0.09
## FormicAcid      0.01       1.00       0.64           0.03           0.14
## aceticacid     -0.03       0.64       1.00          -0.14           0.10
## NonaDecanoicAc -0.20       0.03      -0.14           1.00          -0.23
## Tot_OcNoDecana  0.09       0.14       0.10          -0.23           1.00
## 
## n= 47 
## 
## 
## P
##                B      T      E      X      9_ane  10_ane 13_ane 14_ane 1_M_2_PA
## B                     0.0000 0.1971 0.0436 0.2540 0.0995 0.0000 0.0000 0.0417  
## T              0.0000        0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000  
## E              0.1971 0.0000        0.0000 0.0000 0.0000 0.0350 0.0657 0.0000  
## X              0.0436 0.0000 0.0000        0.0000 0.0000 0.0064 0.0117 0.0000  
## 9_ane          0.2540 0.0000 0.0000 0.0000        0.0000 0.0016 0.0026 0.0000  
## 10_ane         0.0995 0.0000 0.0000 0.0000 0.0000        0.0187 0.0320 0.0000  
## 13_ane         0.0000 0.0000 0.0350 0.0064 0.0016 0.0187        0.0000 0.0073  
## 14_ane         0.0000 0.0000 0.0657 0.0117 0.0026 0.0320 0.0000        0.0167  
## 1_M_2_PA       0.0417 0.0000 0.0000 0.0000 0.0000 0.0000 0.0073 0.0167         
## BTM            0.4363 0.0000 0.0000 0.0000 0.0000 0.0000 0.1489 0.2458 0.0000  
## FormicAcid     0.2779 0.6909 0.8973 0.9157 0.7436 0.8862 0.6859 0.7556 0.9099  
## aceticacid     0.2840 0.7804 0.8590 0.8618 0.8973 0.8947 0.4576 0.7292 0.9055  
## NonaDecanoicAc 0.3058 0.0197 0.0753 0.0530 0.0324 0.1596 0.0079 0.0312 0.0945  
## Tot_OcNoDecana 0.0000 0.0000 0.2563 0.1086 0.0010 0.1817 0.0000 0.0000 0.0680  
##                BTM    FormicAcid aceticacid NonaDecanoicAc Tot_OcNoDecana
## B              0.4363 0.2779     0.2840     0.3058         0.0000        
## T              0.0000 0.6909     0.7804     0.0197         0.0000        
## E              0.0000 0.8973     0.8590     0.0753         0.2563        
## X              0.0000 0.9157     0.8618     0.0530         0.1086        
## 9_ane          0.0000 0.7436     0.8973     0.0324         0.0010        
## 10_ane         0.0000 0.8862     0.8947     0.1596         0.1817        
## 13_ane         0.1489 0.6859     0.4576     0.0079         0.0000        
## 14_ane         0.2458 0.7556     0.7292     0.0312         0.0000        
## 1_M_2_PA       0.0000 0.9099     0.9055     0.0945         0.0680        
## BTM                   0.9500     0.8584     0.1883         0.5528        
## FormicAcid     0.9500            0.0000     0.8247         0.3349        
## aceticacid     0.8584 0.0000                0.3484         0.5029        
## NonaDecanoicAc 0.1883 0.8247     0.3484                    0.1157        
## Tot_OcNoDecana 0.5528 0.3349     0.5029     0.1157

Même commentaire que precedemment, les correlations sont principalement positives sauf pour NonaDecanoicAc et aceticacid ou les correlations sont negatives (mais faibles). En hiver il y a aussi beaucoup plus de correlation tres proche de 1 qu’en été, donc VAR1 = cst*VAR2, par exemple X et E.// Donc l’effet de cause commune des usines ou de l’environnement hivernale sur les mesures est bien plus fort en hiver qu’en été.

# Split the data by 'Periode'
data_by_periode <- split(data, data$Periode)

# Calculate the correlation matrix for the 'avant' period using rcorr
correlations_avant <- lapply(data_by_periode["avant"], function(df) {
  # Select only the relevant numeric columns
  relevant_df <- df[, c("B", "T", "E", "X", "9_ane", "10_ane", "13_ane", "14_ane", "1_M_2_PA", "BTM", "FormicAcid", "aceticacid", "NonaDecanoicAc", "Tot_OcNoDecana")]
  # Use rcorr for a more detailed correlation analysis
  rcorr(as.matrix(relevant_df), type = "pearson")
})

# Print the correlation matrix for the 'avant' period
print(correlations_avant)

## $avant
##                    B    T     E     X 9_ane 10_ane 13_ane 14_ane 1_M_2_PA  BTM
## B               1.00 0.39  0.45  0.53  0.43   0.38   0.49   0.09     0.51 0.61
## T               0.39 1.00  0.76  0.76  0.93   0.91   0.74   0.04     0.36 0.72
## E               0.45 0.76  1.00  0.98  0.67   0.58   0.64  -0.15     0.64 0.73
## X               0.53 0.76  0.98  1.00  0.70   0.61   0.67  -0.15     0.66 0.78
## 9_ane           0.43 0.93  0.67  0.70  1.00   0.98   0.89   0.13     0.40 0.80
## 10_ane          0.38 0.91  0.58  0.61  0.98   1.00   0.86   0.13     0.35 0.78
## 13_ane          0.49 0.74  0.64  0.67  0.89   0.86   1.00   0.23     0.53 0.82
## 14_ane          0.09 0.04 -0.15 -0.15  0.13   0.13   0.23   1.00     0.14 0.09
## 1_M_2_PA        0.51 0.36  0.64  0.66  0.40   0.35   0.53   0.14     1.00 0.82
## BTM             0.61 0.72  0.73  0.78  0.80   0.78   0.82   0.09     0.82 1.00
## FormicAcid      0.43 0.69  0.74  0.76  0.71   0.66   0.73  -0.14     0.41 0.64
## aceticacid      0.21 0.87  0.47  0.48  0.93   0.95   0.78   0.11     0.18 0.64
## NonaDecanoicAc -0.17 0.31  0.05  0.04  0.45   0.46   0.48   0.37     0.03 0.22
## Tot_OcNoDecana  0.27 0.68  0.53  0.52  0.78   0.76   0.84   0.16     0.27 0.55
##                FormicAcid aceticacid NonaDecanoicAc Tot_OcNoDecana
## B                    0.43       0.21          -0.17           0.27
## T                    0.69       0.87           0.31           0.68
## E                    0.74       0.47           0.05           0.53
## X                    0.76       0.48           0.04           0.52
## 9_ane                0.71       0.93           0.45           0.78
## 10_ane               0.66       0.95           0.46           0.76
## 13_ane               0.73       0.78           0.48           0.84
## 14_ane              -0.14       0.11           0.37           0.16
## 1_M_2_PA             0.41       0.18           0.03           0.27
## BTM                  0.64       0.64           0.22           0.55
## FormicAcid           1.00       0.65           0.09           0.67
## aceticacid           0.65       1.00           0.48           0.75
## NonaDecanoicAc       0.09       0.48           1.00           0.58
## Tot_OcNoDecana       0.67       0.75           0.58           1.00
## 
## n= 32 
## 
## 
## P
##                B      T      E      X      9_ane  10_ane 13_ane 14_ane 1_M_2_PA
## B                     0.0295 0.0098 0.0020 0.0141 0.0325 0.0044 0.6152 0.0031  
## T              0.0295        0.0000 0.0000 0.0000 0.0000 0.0000 0.8164 0.0423  
## E              0.0098 0.0000        0.0000 0.0000 0.0005 0.0000 0.4055 0.0000  
## X              0.0020 0.0000 0.0000        0.0000 0.0002 0.0000 0.4071 0.0000  
## 9_ane          0.0141 0.0000 0.0000 0.0000        0.0000 0.0000 0.4865 0.0217  
## 10_ane         0.0325 0.0000 0.0005 0.0002 0.0000        0.0000 0.4647 0.0521  
## 13_ane         0.0044 0.0000 0.0000 0.0000 0.0000 0.0000        0.1985 0.0019  
## 14_ane         0.6152 0.8164 0.4055 0.4071 0.4865 0.4647 0.1985        0.4406  
## 1_M_2_PA       0.0031 0.0423 0.0000 0.0000 0.0217 0.0521 0.0019 0.4406         
## BTM            0.0002 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.6281 0.0000  
## FormicAcid     0.0137 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.4294 0.0185  
## aceticacid     0.2452 0.0000 0.0060 0.0051 0.0000 0.0000 0.0000 0.5347 0.3153  
## NonaDecanoicAc 0.3649 0.0872 0.7802 0.8091 0.0094 0.0076 0.0060 0.0384 0.8597  
## Tot_OcNoDecana 0.1374 0.0000 0.0017 0.0024 0.0000 0.0000 0.0000 0.3954 0.1298  
##                BTM    FormicAcid aceticacid NonaDecanoicAc Tot_OcNoDecana
## B              0.0002 0.0137     0.2452     0.3649         0.1374        
## T              0.0000 0.0000     0.0000     0.0872         0.0000        
## E              0.0000 0.0000     0.0060     0.7802         0.0017        
## X              0.0000 0.0000     0.0051     0.8091         0.0024        
## 9_ane          0.0000 0.0000     0.0000     0.0094         0.0000        
## 10_ane         0.0000 0.0000     0.0000     0.0076         0.0000        
## 13_ane         0.0000 0.0000     0.0000     0.0060         0.0000        
## 14_ane         0.6281 0.4294     0.5347     0.0384         0.3954        
## 1_M_2_PA       0.0000 0.0185     0.3153     0.8597         0.1298        
## BTM                   0.0000     0.0000     0.2279         0.0011        
## FormicAcid     0.0000            0.0000     0.6344         0.0000        
## aceticacid     0.0000 0.0000                0.0058         0.0000        
## NonaDecanoicAc 0.2279 0.6344     0.0058                    0.0005        
## Tot_OcNoDecana 0.0011 0.0000     0.0000     0.0005

Correlations positives comme on pouvait s’y attendre

# Calculate the correlation matrix for the 'apres' period using rcorr
correlations_apres <- lapply(data_by_periode["apres"], function(df) {
  # Select only the relevant numeric columns
  relevant_df <- df[, c("B", "T", "E", "X", "9_ane", "10_ane", "13_ane", "14_ane", "1_M_2_PA", "BTM", "FormicAcid", "aceticacid", "NonaDecanoicAc", "Tot_OcNoDecana")]
  # Use rcorr for a more detailed correlation analysis
  rcorr(as.matrix(relevant_df), type = "pearson")
})

# Print the correlation matrix for the 'apres' period
print(correlations_apres)

## $apres
##                    B     T     E     X 9_ane 10_ane 13_ane 14_ane 1_M_2_PA
## B               1.00  0.30  0.22  0.33  0.36   0.30   0.82   0.82     0.38
## T               0.30  1.00  0.61  0.69  0.48   0.57   0.44   0.42     0.59
## E               0.22  0.61  1.00  0.98  0.69   0.98   0.33   0.29     0.94
## X               0.33  0.69  0.98  1.00  0.71   0.96   0.43   0.40     0.94
## 9_ane           0.36  0.48  0.69  0.71  1.00   0.70   0.55   0.53     0.70
## 10_ane          0.30  0.57  0.98  0.96  0.70   1.00   0.38   0.36     0.95
## 13_ane          0.82  0.44  0.33  0.43  0.55   0.38   1.00   0.96     0.45
## 14_ane          0.82  0.42  0.29  0.40  0.53   0.36   0.96   1.00     0.43
## 1_M_2_PA        0.38  0.59  0.94  0.94  0.70   0.95   0.45   0.43     1.00
## BTM             0.16  0.49  0.98  0.94  0.64   0.98   0.24   0.20     0.93
## FormicAcid      0.42 -0.07  0.08  0.11  0.22   0.10   0.32   0.32     0.14
## aceticacid     -0.03  0.05 -0.02 -0.02 -0.04  -0.02   0.02  -0.04    -0.01
## NonaDecanoicAc  0.44 -0.33 -0.06 -0.02  0.12   0.01   0.22   0.30     0.05
## Tot_OcNoDecana  0.70  0.38  0.21  0.29  0.54   0.25   0.78   0.73     0.34
##                  BTM FormicAcid aceticacid NonaDecanoicAc Tot_OcNoDecana
## B               0.16       0.42      -0.03           0.44           0.70
## T               0.49      -0.07       0.05          -0.33           0.38
## E               0.98       0.08      -0.02          -0.06           0.21
## X               0.94       0.11      -0.02          -0.02           0.29
## 9_ane           0.64       0.22      -0.04           0.12           0.54
## 10_ane          0.98       0.10      -0.02           0.01           0.25
## 13_ane          0.24       0.32       0.02           0.22           0.78
## 14_ane          0.20       0.32      -0.04           0.30           0.73
## 1_M_2_PA        0.93       0.14      -0.01           0.05           0.34
## BTM             1.00       0.07      -0.02          -0.03           0.13
## FormicAcid      0.07       1.00       0.41           0.41           0.33
## aceticacid     -0.02       0.41       1.00          -0.08          -0.01
## NonaDecanoicAc -0.03       0.41      -0.08           1.00           0.24
## Tot_OcNoDecana  0.13       0.33      -0.01           0.24           1.00
## 
## n= 77 
## 
## 
## P
##                B      T      E      X      9_ane  10_ane 13_ane 14_ane 1_M_2_PA
## B                     0.0074 0.0528 0.0033 0.0013 0.0085 0.0000 0.0000 0.0006  
## T              0.0074        0.0000 0.0000 0.0000 0.0000 0.0000 0.0002 0.0000  
## E              0.0528 0.0000        0.0000 0.0000 0.0000 0.0037 0.0103 0.0000  
## X              0.0033 0.0000 0.0000        0.0000 0.0000 0.0001 0.0004 0.0000  
## 9_ane          0.0013 0.0000 0.0000 0.0000        0.0000 0.0000 0.0000 0.0000  
## 10_ane         0.0085 0.0000 0.0000 0.0000 0.0000        0.0007 0.0014 0.0000  
## 13_ane         0.0000 0.0000 0.0037 0.0001 0.0000 0.0007        0.0000 0.0000  
## 14_ane         0.0000 0.0002 0.0103 0.0004 0.0000 0.0014 0.0000        0.0001  
## 1_M_2_PA       0.0006 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001         
## BTM            0.1686 0.0000 0.0000 0.0000 0.0000 0.0000 0.0347 0.0741 0.0000  
## FormicAcid     0.0002 0.5561 0.5020 0.3574 0.0532 0.3825 0.0047 0.0047 0.2115  
## aceticacid     0.8230 0.6734 0.8599 0.8895 0.7399 0.8317 0.8589 0.7112 0.9045  
## NonaDecanoicAc 0.0000 0.0037 0.6118 0.8338 0.3161 0.9627 0.0559 0.0090 0.6419  
## Tot_OcNoDecana 0.0000 0.0006 0.0718 0.0105 0.0000 0.0302 0.0000 0.0000 0.0023  
##                BTM    FormicAcid aceticacid NonaDecanoicAc Tot_OcNoDecana
## B              0.1686 0.0002     0.8230     0.0000         0.0000        
## T              0.0000 0.5561     0.6734     0.0037         0.0006        
## E              0.0000 0.5020     0.8599     0.6118         0.0718        
## X              0.0000 0.3574     0.8895     0.8338         0.0105        
## 9_ane          0.0000 0.0532     0.7399     0.3161         0.0000        
## 10_ane         0.0000 0.3825     0.8317     0.9627         0.0302        
## 13_ane         0.0347 0.0047     0.8589     0.0559         0.0000        
## 14_ane         0.0741 0.0047     0.7112     0.0090         0.0000        
## 1_M_2_PA       0.0000 0.2115     0.9045     0.6419         0.0023        
## BTM                   0.5652     0.8851     0.7967         0.2611        
## FormicAcid     0.5652            0.0002     0.0002         0.0031        
## aceticacid     0.8851 0.0002                0.4825         0.9144        
## NonaDecanoicAc 0.7967 0.0002     0.4825                    0.0331        
## Tot_OcNoDecana 0.2611 0.0031     0.9144     0.0331

Correlations positives comme on pouvait s’y attendre

#Mettre en œuvre une ACP avec analyse et interprétation des résultats (avec représentativité de l’ACP )

Avant d’appliquer l’ACP sur les donnees :

• Taille de l’echantillon : méthode ou on supprime les outliers : regle des 5 respectée car il y a plus de 5x14 = 70 individus > 76 individus (33 outliers + le reste de NA et 0) • Taille de l’échantillon : pas de problème //

• Test de KMO : corrélations partielles entre variables petites ou pas ?

DatabaseACP <- data_numeric
psych::KMO(DatabaseACP)

## Kaiser-Meyer-Olkin factor adequacy
## Call: psych::KMO(r = DatabaseACP)
## Overall MSA =  0.81
## MSA for each item = 
##              B              T              E              X         X9_ane 
##           0.89           0.78           0.76           0.83           0.90 
##        X10_ane        X13_ane        X14_ane      X1_M_2_PA            BTM 
##           0.83           0.77           0.79           0.95           0.71 
##     FormicAcid     aceticacid NonaDecanoicAc Tot_OcNoDecana 
##           0.79           0.41           0.73           0.85

Indice KMO de 81% soit “bon”.

• Test de Bartlett : vérifier si les variables sont corrélées dans l’ensemble de données et s’il est donc approprié de procéder à une réduction de dimension par ACP.

BARTLETT(DatabaseACP, use="complete.obs")

## ℹ 'x' was not a correlation matrix. Correlations are found from entered raw data.

## 
## ✔ The Bartlett's test of sphericity was significant at an alpha level of .05.
##   These data are probably suitable for factor analysis.
## 
##   𝜒²(91) = 2500.57, p < .001

pvalue < 5% : rejeter l’hypothèse nulle que la matrice de corrélations est égale à la matrice d’identité. Cela signifie qu’au moins une des variables dans l’ensemble de données est corrélée avec une autre, validant ainsi l’utilisation de l’ACP pour ces données.

Mise en oeuvre de l’ACP

ACP_test <- prcomp(DatabaseACP, scale=TRUE)
ACP_test

## Standard deviations (1, .., p=14):
##  [1] 2.70681870 1.73818904 1.09412628 0.97579595 0.63723734 0.60484105
##  [7] 0.53269577 0.45874442 0.32750804 0.23348479 0.19616364 0.16100819
## [13] 0.08790143 0.05018823
## 
## Rotation (n x k) = (14 x 14):
##                       PC1          PC2         PC3         PC4          PC5
## B              0.26059276  0.346941039 -0.05046069 -0.06361327  0.384982120
## T              0.25351256 -0.087691536 -0.15082432  0.53518465  0.517089466
## E              0.30694518 -0.308049822  0.07467670 -0.06458986 -0.006744262
## X              0.33101409 -0.234569833  0.03770182 -0.01890689  0.098402330
## X9_ane         0.31344764 -0.002847728 -0.05067527 -0.01555145 -0.473761943
## X10_ane        0.32309128 -0.253789394  0.07226998 -0.10807600  0.018561309
## X13_ane        0.28016754  0.271718969 -0.24718835  0.15838355 -0.142891709
## X14_ane        0.27982428  0.305537225 -0.22260329  0.06652200  0.045449542
## X1_M_2_PA      0.33337282 -0.196964975  0.04902381 -0.09944134  0.015225240
## BTM            0.28249484 -0.340073991  0.13035292 -0.16471834 -0.063824617
## FormicAcid     0.17065334  0.323538215  0.48312446 -0.11304139 -0.210320684
## aceticacid     0.05860446  0.147236587  0.71594690  0.49513392  0.003283061
## NonaDecanoicAc 0.14184506  0.345639570  0.15806705 -0.58817498  0.362406669
## Tot_OcNoDecana 0.24695427  0.300258602 -0.24267846  0.15053953 -0.384071078
##                        PC6          PC7           PC8          PC9        PC10
## B              -0.09722808  0.061157289  0.0717957832  0.776498341  0.14022292
## T               0.47147984  0.178628854 -0.0434742217 -0.213191263 -0.11061726
## E              -0.04508246  0.052506841  0.0226953582 -0.090104408  0.32888030
## X              -0.03521193  0.064118389 -0.0004483977 -0.093513560  0.41357362
## X9_ane          0.53568955 -0.325544413 -0.4530021342  0.257196220 -0.03964463
## X10_ane        -0.10138289 -0.023424679  0.0270756129  0.066361745  0.01714446
## X13_ane        -0.43593549 -0.043476676 -0.1810387132 -0.270834073  0.09858360
## X14_ane        -0.31043668 -0.055780580 -0.3590776935 -0.172301116 -0.12515554
## X1_M_2_PA      -0.15075265 -0.008059655  0.1879206854  0.053477579 -0.80372234
## BTM            -0.14048316 -0.009601323  0.1121931071  0.004888223  0.06593168
## FormicAcid      0.12767421  0.726763454 -0.1348922153 -0.094593559 -0.04978094
## aceticacid     -0.14624948 -0.432971976  0.0746420316  0.007977544  0.02457947
## NonaDecanoicAc  0.26120376 -0.360514415  0.0825817951 -0.379921458  0.01563885
## Tot_OcNoDecana  0.18152078 -0.013272911  0.7391181748 -0.070057356  0.08946652
##                        PC11        PC12          PC13          PC14
## B              -0.046708303  0.09071283 -7.332675e-02  0.0130989048
## T               0.121097291  0.12721100 -7.394214e-02 -0.0373716923
## E              -0.196513300 -0.03508278 -3.904337e-01  0.6980141443
## X              -0.501943033 -0.26098780  3.995641e-01 -0.3988641763
## X9_ane         -0.060073364  0.06386319 -4.192434e-03 -0.0300395668
## X10_ane         0.639880597 -0.01514251  5.800653e-01  0.2270345634
## X13_ane        -0.090883416  0.64373485  1.031508e-01  0.0004915407
## X14_ane         0.173305557 -0.65430599 -1.998353e-01 -0.0211241147
## X1_M_2_PA      -0.352579573  0.02212975  5.347407e-02  0.0542875542
## BTM             0.334386470  0.15817580 -5.355291e-01 -0.5440014116
## FormicAcid      0.028821272  0.02294182  1.482229e-02 -0.0096725770
## aceticacid     -0.011407325 -0.04220537 -5.872887e-03  0.0099890654
## NonaDecanoicAc -0.000593936  0.09165122  5.406512e-05  0.0116909820
## Tot_OcNoDecana  0.064428910 -0.15234024 -9.646162e-03  0.0026770442

On peut ici voir quelle variable a le plus d’influence sur quelle CP: B a le plus d’influence sur PC9 T a une grande influence sur PC4 …

Combien de CP retient-on ?

get_eig(ACP_test)

##         eigenvalue variance.percent cumulative.variance.percent
## Dim.1  7.326867469      52.33476763                    52.33477
## Dim.2  3.021301131      21.58072236                    73.91549
## Dim.3  1.197112315       8.55080225                    82.46629
## Dim.4  0.952177739       6.80126956                    89.26756
## Dim.5  0.406071424       2.90051017                    92.16807
## Dim.6  0.365832695       2.61309068                    94.78116
## Dim.7  0.283764788       2.02689134                    96.80805
## Dim.8  0.210446443       1.50318888                    98.31124
## Dim.9  0.107261514       0.76615367                    99.07740
## Dim.10 0.054515149       0.38939392                    99.46679
## Dim.11 0.038480175       0.27485839                    99.74165
## Dim.12 0.025923638       0.18516884                    99.92682
## Dim.13 0.007726662       0.05519044                    99.98201
## Dim.14 0.002518858       0.01799184                   100.00000

• 1ere règle : on retient les inerties > 1. Ce qui donne ici qu’on retient les 3 premiers PC. • 2eme règle : on retient au seuil 80% d’inertie expliquée. Ce qui donne aussi les 3 premieres PC. • 3eme règle : on cherche les coudes dans le scree plot. Ici, un coude est visible pour la 2eme ou 3eme PC.

fviz_eig(ACP_test, addlabels=TRUE)

CCL : On retient les 3 premiers composantes principales CP. q=3.

On a extrait 80% des informations avec 3 composantes au lieu de 100% avec 14.

res.pca <- PCA(data_numeric, graph = FALSE)
fviz_pca_var(res.pca, col.var = "cos2", axes = 1:2)

#res.pca_noNA <- PCA(data_noNA, quali.sup = which(names(data_noNA) %in% c("TYPE", "SAISON", "Campagne", "Localisation")), graph = FALSE)
#fviz_pca_var(res.pca, col.var = "cos2", axes = 1:2)

DatabaseACP <- data_numeric

ACP_3factNONE <- principal(DatabaseACP, nfactors=3, 
                           rotate="none")

ACP_3factVAR <- principal(DatabaseACP, nfactors=3, 
                           rotate="varimax") #rotation Varimax visant à simplifier l'interprétation des facteurs.

loadings(ACP_3factVAR)

## 
## Loadings:
##                RC1    RC2    RC3   
## B               0.186  0.873  0.259
## T               0.620  0.357       
## E               0.988              
## X               0.959  0.224       
## X9_ane          0.668  0.521       
## X10_ane         0.965  0.174       
## X13_ane         0.288  0.888       
## X14_ane         0.254  0.920       
## X1_M_2_PA       0.925  0.271       
## BTM             0.975              
## FormicAcid             0.490  0.751
## aceticacid                    0.838
## NonaDecanoicAc         0.599  0.421
## Tot_OcNoDecana  0.187  0.868       
## 
##                  RC1   RC2   RC3
## SS loadings    5.690 4.309 1.546
## Proportion Var 0.406 0.308 0.110
## Cumulative Var 0.406 0.714 0.825

#loadings(ACP_3factNONE)

CP1 corrélé à : T, E, X, 10_ane, 1_M_2_PA et BTM # 6 variable PC1 CP2 corrélé à : B, 9_ane, 13_ane, 14_ane, NonaDecanoicAc et Tot_OcNoDecana. # 6 variables PC2 CP3 corrélé à : FormicAcid et aceticacid. # 2 variables PC3

Ceux proches de 1 : mieux expliquées par cette PC (Ex: aceticacid pour PC3)

On peut deja en dire que les acides sont corrélées grandement avec PC3. Donc on peut interpreter que la PC3 représente les acides.

PC1 represente les solvants pétroliers (à cycle benzeique notamment)

PC2 represente les hydrocarbures et leurs dérivés.

On va tracé en 2D (ie on garde que 2 CP pour meilleur visibilité avec R)

DatabaseACP <- data_numeric
DatabaseACP <- cbind(DatabaseACP, ACP_3factVAR$scores)
#view(DatabaseACP)

plot(ACP_3factVAR$loadings, type="n", xlim=c(-1.1,1.1), ylim=c(-1,1))
abline(v=0, col="indianred") # Ajout Ligne verticale
abline(h=0, col="indianred") # Ajout Ligne horizontale
text(ACP_3factVAR$loadings, labels=names(DatabaseACP), cex=1)
title(main = "Representation des charges en 2D pour les 2ère PC ") 

On voit bien les 3 composantes : les solvant, les hydrocarbures et les acides

# Assurer que 'ACP_3factVAR$loadings' contient bien les données des composantes principales de votre ACP.
# Le code ci-dessous suppose que 'ACP_3factVAR' est un objet contenant les résultats de l'ACP avec les loadings.

# Tracer RC1 contre RC3
plot(ACP_3factVAR$loadings[, c(1, 3)], type="n", xlim=c(-1.1, 1.1), ylim=c(-1, 1))
abline(v=0, col="indianred") # Ajout d'une ligne verticale à x=0
abline(h=0, col="indianred") # Ajout d'une ligne horizontale à y=0

# Ajouter les noms des variables aux chargements correspondants pour RC1 et RC3
text(ACP_3factVAR$loadings[, 1], ACP_3factVAR$loadings[, 3], labels=names(DatabaseACP), cex=1)

# Titre du graphique
title(main = "PC1 vs PC3")

Ici on voti bien sur le 3eme axe les acides.

library(plotly)

# Définition des étiquettes des individus
var_labels <- c("B", "T", "E", "X", "9_ane", "10_ane", "13_ane", "14_ane", 
                "1_M_2_PA", "BTM", "FormicAcid", "aceticacid", "NonaDecanoicAc", 
                "Tot_OcNoDecana")


# Utiliser plotly pour une visualisation 3D interactive
plot_ly(x = ~ACP_3factVAR$loadings[, 1], 
        y = ~ACP_3factVAR$loadings[, 2], 
        z = ~ACP_3factVAR$loadings[, 3], 
        text = ~var_labels, 
        type = 'scatter3d', 
        mode = 'markers+text',
        marker = list(size = 5, color = 'blue')) %>%
  layout(title = "Charges ACP: PC1 vs PC2 vs PC3",
         scene = list(xaxis = list(title = 'PC1'),
                      yaxis = list(title = 'PC2'),
                      zaxis = list(title = 'PC3')))

on peut ici aussi voir les 3 composantes solvants, hydrocarbures et acides.

#Et voici les individus ou l’on prend en compte la rotation

library(FactoMineR)
library(plotly)
library(stats)

# Réaliser une PCA avec FactoMineR
res.pca <- PCA(data_numeric, graph = FALSE)

# Appliquer la rotation Varimax sur les charges des composantes
rotation <- varimax(res.pca$var$coord)
rotated_scores <- res.pca$ind$coord %*% rotation$rotmat  # Les scores après rotation Varimax

# Créer un graphique 3D avec Plotly en utilisant les scores tournés
plot_ly(x = rotated_scores[, 1],  # Premier axe principal après rotation
        y = rotated_scores[, 2],  # Deuxième axe principal après rotation
        z = rotated_scores[, 3],  # Troisième axe principal après rotation
        type = 'scatter3d', 
        mode = 'markers',
        marker = list(size = 5, color = 'blue', opacity = 0.5)) %>%
  layout(title = "PCA avec rotation Varimax: Vue 3D des scores des individus",
         scene = list(xaxis = list(title = 'PC1'),
                      yaxis = list(title = 'PC2'),
                      zaxis = list(title = 'PC3')))

Sans rotation, c’est plus facile de travailler sans car moins de calculs matricielle cara on n’a aps de fonction toute faite (non trouvée qui donne les coordonnées)

fviz_pca_ind (res.pca, 
    pointsize = "cos2",
    pointshape = 21,
    fill = "#E7B800",
    repel = TRUE ,
    axes = c(1, 2),
    #xlim = c(-4, 4),
    #ylim = c(-2, 4)
    )

#En zoomant un peu plus

fviz_pca_ind (res.pca, 
    pointsize = "cos2",
    pointshape = 21,
    fill = "#E7B800",
    repel = TRUE ,
    axes = c(1, 2),
    xlim = c(-4, 4),
    ylim = c(-2, 4))

On semble voir des groupes/classes, esseyons de voir cela

#install.packages("rgl")
library(rgl)
 
# Supposons que ACP_3factVAR est un objet contenant les résultats d'une ACP avec les loadings
# et que DatabaseACP contient les noms de vos variables.

# Ouvrir une nouvelle fenêtre de graphique 3D
open3d()

## glX 
##   3

# Tracer les points dans l'espace 3D
plot3d(ACP_3factVAR$loadings[, 1],  # Charges de la première composante principale
       ACP_3factVAR$loadings[, 2],  # Charges de la deuxième composante principale
       ACP_3factVAR$loadings[, 3],  # Charges de la troisième composante principale
       col = "blue",                 # Couleur des points
       size = 3)                     # Taille des points

# Ajouter des étiquettes aux points
text3d(ACP_3factVAR$loadings[, 1],
       ACP_3factVAR$loadings[, 2],
       ACP_3factVAR$loadings[, 3],
       text = names(DatabaseACP),
       adj = c(1.5, 1.5))

# Ajouter des axes au graphique (facultatif)
axes3d(edges = c("x--", "y--", "z--"), 
       col = "black", 
       tick = TRUE, 
       labels = TRUE)


# Titre du graphique
title3d(main = "3D Plot of Principal Component Loadings")

Afficher les individus et les les variables dans un meme graphique

library(FactoMineR)
library(factoextra)

res.pca <- PCA(data_numeric, graph = FALSE)
fviz_pca_biplot(res.pca)

# Extraire les scores des composantes pour les 3 premières composantes principales
scores <- res.pca$ind$coord[, 1:3]

# Créer un graphique 3D avec rgl
plot3d(scores[,1], scores[,2], scores[,3], type = "s", col = "blue", size = 1)

# Ajouter des étiquettes aux points
text3d(scores[,1], scores[,2], scores[,3], texts = 1:nrow(scores), col = "red")

# Ajouter des axes avec des étiquettes
axes3d(col = "black", box = FALSE)
title3d("CP 1", "C2 2", "CP 3")

# Vous pouvez personnaliser l'apparence avec diverses fonctions 'rgl'
rglwidget() # Si vous utilisez RStudio Viewer ou un environnement HTML pour intégrer le graphique 3D

# Charger les packages nécessaires
library(FactoMineR)
library(factoextra)
library(plotly)

# Assurez-vous que vous avez les données numériques correctes pour l'ACP
# Remplacez 'data_numeric' par votre propre jeu de données numériques
res.pca <- PCA(data_numeric, graph = FALSE)

# Maintenant, nous pouvons accéder aux coordonnées des individus et aux charges des variables
ind_coords <- res.pca$ind$coord  # Coordonnées des individus
var_coords <- res.pca$var$coord  # Charges des variables

# Définition des étiquettes pour les variables comme précédemment
var_labels <- c("B", "T", "E", "X", "9_ane", "10_ane", "13_ane", "14_ane", 
                "1_M_2_PA", "BTM", "FormicAcid", "aceticacid", "NonaDecanoicAc", 
                "Tot_OcNoDecana")

# Générer des étiquettes pour les individus de 1 à n
ind_labels <- as.character(1:nrow(ind_coords))

# Visualisation 3D avec plotly
p <- plot_ly() %>%
  add_trace(x = var_coords[, 1], y = var_coords[, 2], z = var_coords[, 3],
            text = var_labels, type = 'scatter3d', mode = 'markers+text',
            marker = list(size = 5, color = 'red'), name = "Variables") %>%
  add_trace(x = ind_coords[, 1], y = ind_coords[, 2], z = ind_coords[, 3],
            text = ind_labels, type = 'scatter3d', mode = 'markers',
            marker = list(size = 3, color = 'blue'), name = "Individus") %>%
  layout(title = "ACP: Charges des variables et positions des individus",
         scene = list(xaxis = list(title = 'PC1'),
                      yaxis = list(title = 'PC2'),
                      zaxis = list(title = 'PC3')))

# Afficher le graphique
p

On revoit les groupes de tout à l’heure, on va maitnenant labeliser :

library(plotly)
library(FactoMineR)
library(dplyr)
library(stringr)


data_cleanACP <- as.data.frame(data_clean)
# Ajouter la colonne 'Avant_Apres' à 'data_cleanACP'
data_cleanACP <- data_cleanACP %>%
  mutate(Avant_Apres = ifelse(str_detect(Campagne, "BF"), "avant", "apres"))

# Coordonnées des individus
ind_coords <- res.pca$ind$coord

# Création de vecteurs de couleurs pour les individus en fonction de 'Avant_Apres'
data_cleanACP$color <- ifelse(data_cleanACP$Avant_Apres == "avant", "blue", "green")

# Séparation des individus en deux groupes : avant et après
ind_avant <- ind_coords[data_cleanACP$Avant_Apres == "avant", ]
ind_apres <- ind_coords[data_cleanACP$Avant_Apres == "apres", ]

# Création du graphique Plotly sans les variables
p <- plot_ly() %>%
  add_trace(x = ind_avant[, 1], y = ind_avant[, 2], z = ind_avant[, 3],
            type = 'scatter3d', mode = 'markers',
            marker = list(size = 3, color = 'blue', opacity = 0.7), name = "Avant") %>%
  add_trace(x = ind_apres[, 1], y = ind_apres[, 2], z = ind_apres[, 3],
            type = 'scatter3d', mode = 'markers',
            marker = list(size = 3, color = 'green', opacity = 0.7), name = "Apres") %>%
  layout(title = "ACP: avant (bleu) VS apres (vert)",
         scene = list(xaxis = list(title = 'PC1'),
                      yaxis = list(title = 'PC2'),
                      zaxis = list(title = 'PC3')),
         legend = list(x = 0.1, y = 0.9))

# Afficher le graphique
p

On voit clairement les 3 groupes AVANT et APRES. Le groupe Avant est plus petit, recentré et compacte: Cela veut dire qu’avant usine les composant n’avait pas de variations (variance non captée) par rapport a aprés, donc il y a un impact significatif des usines. Cela soulgine nos resultats statistique précédents avant / apres.

library(plotly)
library(FactoMineR)
library(dplyr)
library(stringr)

# Supposons que res.pca est déjà défini et que data_cleanACP est votre dataframe avec les résultats PCA

# Coordonnées des individus
ind_coords <- res.pca$ind$coord

# Création de vecteurs de couleurs pour les individus en fonction de 'SAISON'
data_cleanACP$color <- ifelse(data_cleanACP$SAISON == "été", "orange", "blue")

# Séparation des individus en deux groupes : été et hiver
ind_ete <- ind_coords[data_cleanACP$SAISON == "été", ]
ind_hiver <- ind_coords[data_cleanACP$SAISON == "hiver", ]

# Création du graphique Plotly sans les variables
p <- plot_ly() %>%
  add_trace(x = ind_ete[, 1], y = ind_ete[, 2], z = ind_ete[, 3],
            type = 'scatter3d', mode = 'markers',
            marker = list(size = 3, color = 'orange', opacity = 0.7), name = "Été") %>%
  add_trace(x = ind_hiver[, 1], y = ind_hiver[, 2], z = ind_hiver[, 3],
            type = 'scatter3d', mode = 'markers',
            marker = list(size = 3, color = 'blue', opacity = 0.7), name = "Hiver") %>%
  layout(title = "ACP: Été VS Hiver",
         scene = list(xaxis = list(title = 'PC1'),
                      yaxis = list(title = 'PC2'),
                      zaxis = list(title = 'PC3')),
         legend = list(x = 0.1, y = 0.9))

# Afficher le graphique
p

De même, il y a clairement 2 groupes été et hiver. Celui été est moins compacte. Donc les conditions environnemental en été jouent en faveur d’une concentration plus elevé (moins de vent, temperature, secheresse etc..) Et on retrouve ce constat dans notre vie de tout les jours par exemple les pics de polution a Paris sont en été !!!! Il a plusde variance dans les composant en été qu’en hiver.

library(plotly)
library(FactoMineR)
library(dplyr)
library(stringr)

# Supposons que res.pca est déjà défini et que data_cleanACP est votre dataframe avec les résultats PCA

# Coordonnées des individus
ind_coords <- res.pca$ind$coord

# Définir une palette de couleurs pour les différents types
type_colors <- c("urbain" = "red", "sourceindustrie" = "blue", "rural" = "green", "compostage" = "orange")

# Création du graphique Plotly pour les individus en utilisant une boucle pour ajouter une trace par type
p <- plot_ly() %>%
  layout(title = "ACP par type",
         scene = list(xaxis = list(title = 'PC1'),
                      yaxis = list(title = 'PC2'),
                      zaxis = list(title = 'PC3')))

# Ajouter une trace pour chaque type
for (type in names(type_colors)) {
  # Filtrer les données par type
  type_data <- data_cleanACP[data_cleanACP$TYPE == type, ]
  type_coords <- ind_coords[rownames(ind_coords) %in% rownames(type_data), ]
  
  # Ajouter la trace pour ce type spécifique
  p <- p %>%
    add_trace(x = type_coords[, 1], y = type_coords[, 2], z = type_coords[, 3],
              type = 'scatter3d', mode = 'markers',
              marker = list(size = 3, color = type_colors[type], opacity = 0.7),
              name = type)  # Le nom de la trace pour la légende
}

# Afficher le graphique
p

on n’observe pas vraiment de groupe, on ne peut pas plus en dire sur le TYPE.

library(FactoMineR)
library(factoextra)
library(dplyr)
library(stringr)

data_cleanACP <- as.data.frame(data_clean)
#view(data_cleanACP)


# Effectuer l'ACP en spécifiant la variable 'Localisation' comme variable supplémentaire (quali.sup)
res.pca <- PCA(data_cleanACP, quali.sup = which(names(data_cleanACP) %in% c("TYPE", "SAISON", "Campagne", "Localisation", "Avant_Apres")), graph = FALSE)

# Utiliser fviz_pca_biplot pour créer un biplot qui affiche les individus et les variables
# Avec un habillage selon la variable 'Localisation'
fviz_pca_ind(res.pca, 
                label = "none",
                habillage="SAISON",
                palette = c("été" = "red", "hiver" = "blue"),
                addEllipses = TRUE,  # Ajouter des ellipses de confiance autour des groupes
                ellipse.level = 0.95  # Niveau de confiance pour les ellipses
)

On remarque que la valeur en bas a droite tire l’ellipsoide dans sa direction.

Est une anomalie ou un outlier, trés difficile à dire.

library(FactoMineR)
library(factoextra)
library(dplyr)
library(stringr)

data_cleanACP <- as.data.frame(data_clean)
# Ajouter la colonne 'Avant_Apres' à 'data_cleanACP'
data_cleanACP <- data_cleanACP %>%
  mutate(Avant_Apres = ifelse(str_detect(Campagne, "BF"), "avant", "apres"))


#view(data_cleanACP)


# Effectuer l'ACP en spécifiant la variable 'Localisation' comme variable supplémentaire (quali.sup)
res.pca <- PCA(data_cleanACP, quali.sup = which(names(data_cleanACP) %in% c("TYPE", "SAISON", "Campagne", "Localisation", "Avant_Apres")), graph = FALSE)

# Utiliser fviz_pca_biplot pour créer un biplot qui affiche les individus et les variables
# Avec un habillage selon la variable 'Localisation'
fviz_pca_ind(res.pca, 
                label = "none",
                habillage="Avant_Apres",
                palette = c("avant" = "blue", "apres" = "green"),
                addEllipses = TRUE,  # Ajouter des ellipses de confiance autour des groupes
                ellipse.level = 0.95  # Niveau de confiance pour les ellipses
)

De même, la valeur extreme capte vers elle l’ellipsoide verte.

library(FactoMineR)
library(factoextra)
library(dplyr)
library(stringr)

# Assurez-vous que data_cleanACP est déjà chargé et préparé avant cette partie du code
res.pca <- PCA(data_cleanACP, quali.sup = which(names(data_cleanACP) %in% c("TYPE", "SAISON", "Campagne", "Localisation", "Avant_Apres", "color")), graph = FALSE)

# Utiliser fviz_pca_ind pour créer un biplot qui affiche les individus
fviz_pca_ind(res.pca, 
             label = "none",
             habillage = "Localisation",
             #palette = c("été" = "red", "hiver" = "blue"), # Décommentez et ajustez selon vos données
             #addEllipses = TRUE,  # Ajouter des ellipses de confiance autour des groupes
             #ellipse.level = 0.95,  # Niveau de confiance pour les ellipses
             xlim = c(-4, 4),
             ylim = c(-2, 4)
)

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '28'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '29'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '26'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '27'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '28'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '29'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '26'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '27'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '26'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '27'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '28'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '29'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '26'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '26'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '27'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '27'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '28'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '28'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '29'

## Warning in grid.Call.graphics(C_points, x$x, x$y, x$pch, x$size): unimplemented
## pch value '29'

# Selon la localisation, on ne pas plus en dire contrairement au cas ete hiver ou avant apres.

library(FactoMineR)
library(factoextra)
library(dplyr)
library(stringr)

data_cleanACP <- as.data.frame(data_cleanACP)
#view(data_cleanACP)


# Effectuer l'ACP en spécifiant la variable 'Localisation' comme variable supplémentaire (quali.sup)
res.pca <- PCA(data_cleanACP, quali.sup = which(names(data_cleanACP) %in% c("TYPE", "SAISON", "Campagne", "Localisation", "Avant_Apres")), graph = FALSE)

# Utiliser fviz_pca_biplot pour créer un biplot qui affiche les individus et les variables
# Avec un habillage selon la variable 'Localisation'
fviz_pca_ind(res.pca, 
                label = "none",
                habillage="TYPE",
                #palette = c("été" = "red", "hiver" = "blue"),
                addEllipses = TRUE,  # Ajouter des ellipses de confiance autour des groupes
                ellipse.level = 0.95,  # Niveau de confiance pour les ellipses
                xlim = c(-4, 4),
                ylim = c(-2, 4)
)

On a ici un probleme dans le biplot : en effet, le graphique des variables et individus est trés difficilement interpretable car certaines données se comportent comme des valeurs extrêmes ou aberrantes par rapport au reste des données. Donc les points des invidus par rapport aux varibales sont peu interpretable car les point ont tendance à etre proche de 0 donc a cause de ces valeurs extrêmes, comme l’individu 59 en bas a droite.

#Quelques autres tests statistique :

Database <- data_clean
Database <- cbind(Database, ACP_3factVAR$scores)
#view(DatabaseACP)
#install.packages("car")
#install.packages("data.table")
library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:psych':
## 
##     logit

## The following object is masked from 'package:purrr':
## 
##     some

## The following object is masked from 'package:dplyr':
## 
##     recode

library(data.table)

## 
## Attaching package: 'data.table'

## The following objects are masked from 'package:lubridate':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year

## The following object is masked from 'package:purrr':
## 
##     transpose

## The following objects are masked from 'package:reshape2':
## 
##     dcast, melt

## The following objects are masked from 'package:dplyr':
## 
##     between, first, last

t.test(Database$RC1~Database$SAISON,
       alternative="two.sided", var.equal=TRUE)

## 
##  Two Sample t-test
## 
## data:  Database$RC1 by Database$SAISON
## t = 0.96158, df = 107, p-value = 0.3384
## alternative hypothesis: true difference in means between group été and group hiver is not equal to 0
## 95 percent confidence interval:
##  -0.1974978  0.5695770
## sample estimates:
##   mean in group été mean in group hiver 
##          0.10582070         -0.08021892

# pvalue > 5% : Il n'y a pas de difference entre hiver et été pour la composante 1

t.test(Database$RC2~Database$SAISON,
       alternative="two.sided", var.equal=TRUE)

## 
##  Two Sample t-test
## 
## data:  Database$RC2 by Database$SAISON
## t = 11.249, df = 107, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group été and group hiver is not equal to 0
## 95 percent confidence interval:
##  1.218769 1.740224
## sample estimates:
##   mean in group été mean in group hiver 
##           0.8415485          -0.6379481

# difference significative entre hiver et été composante 2

t.test(Database$RC3~Database$SAISON,
       alternative="two.sided", var.equal=TRUE)

## 
##  Two Sample t-test
## 
## data:  Database$RC3 by Database$SAISON
## t = 3.9483, df = 107, p-value = 0.0001412
## alternative hypothesis: true difference in means between group été and group hiver is not equal to 0
## 95 percent confidence interval:
##  0.3568829 1.0766163
## sample estimates:
##   mean in group été mean in group hiver 
##           0.4076924          -0.3090572

# difference significative entre hiver et été composante 3

pvalue > 5% : Il n’y a pas de difference entre hiver et été pour la CP 1 difference significative entre hiver et été CP 2 difference significative entre hiver et été CP 3

Donc entre l’été et l’hiver, il y a les hydrocarbures et les acides qui sont statistiquement de moyenne different. Pour les solvant ce n’est pas prouvé par la statistique.

library(dplyr)
library(stringr)
Database <- data_clean
Database <- cbind(Database, ACP_3factVAR$scores)
# Ajout de la nouvelle colonne 'Avant_Apres' à 'Database'
Database <- Database %>%
  mutate(Avant_Apres = case_when(
    str_detect(Campagne, "BF") ~ "avant",
    TRUE ~ "apres"
  ))

# Afficher les premières lignes de Database pour vérifier
#head(Database)


t.test(Database$RC1~Database$Avant_Apres,
       alternative="two.sided", var.equal=TRUE)

## 
##  Two Sample t-test
## 
## data:  Database$RC1 by Database$Avant_Apres
## t = 1.143, df = 107, p-value = 0.2556
## alternative hypothesis: true difference in means between group apres and group avant is not equal to 0
## 95 percent confidence interval:
##  -0.1762851  0.6564274
## sample estimates:
## mean in group apres mean in group avant 
##           0.0704796          -0.1695915

t.test(Database$RC2~Database$Avant_Apres,
       alternative="two.sided", var.equal=TRUE)

## 
##  Two Sample t-test
## 
## data:  Database$RC2 by Database$Avant_Apres
## t = 6.5406, df = 107, p-value = 2.158e-09
## alternative hypothesis: true difference in means between group apres and group avant is not equal to 0
## 95 percent confidence interval:
##  0.8140972 1.5221989
## sample estimates:
## mean in group apres mean in group avant 
##           0.3429425          -0.8252055

t.test(Database$RC3~Database$Avant_Apres,
       alternative="two.sided", var.equal=TRUE)

## 
##  Two Sample t-test
## 
## data:  Database$RC3 by Database$Avant_Apres
## t = 4.1708, df = 107, p-value = 6.189e-05
## alternative hypothesis: true difference in means between group apres and group avant is not equal to 0
## 95 percent confidence interval:
##  0.4288727 1.2058713
## sample estimates:
## mean in group apres mean in group avant 
##           0.2399624          -0.5774096

pvalue > 5% : Il n’y a pas de difference entre avant et apres pour la CP 1 difference significative entre avant et apres CP 2 difference significative entre avant et apres CP 3

Donc entre avant et aprés, il y a les hydrocarbures et les acides qui sont statistiquement de moyenne different. Pour les solvant ce n’est pas prouvé par la statistique.

On retrouve le même resultat avec été et hiver.

Test ANOVA

#install.packages("ez")
#install.packages("DescTools")
library(ez)
library(DescTools)

## 
## Attaching package: 'DescTools'

## The following object is masked from 'package:data.table':
## 
##     %like%

## The following object is masked from 'package:car':
## 
##     Recode

## The following objects are masked from 'package:psych':
## 
##     AUC, ICC, SD

## The following objects are masked from 'package:Hmisc':
## 
##     %nin%, Label, Mean, Quantile

Database <- data_clean
Database <- cbind(Database, ACP_3factVAR$scores)
Database <- cbind(Database, ID=c(1:nrow(Database)))
#Database <- as.matrix(Database)
#view(Database)
#Database$ID.1 <- NULL

ezANOVA(data = Database,
        dv = .(X),
        wid = .(ID),
        between = .(Campagne),
        detailed = TRUE, type = 3)

## Warning: Converting "ID" to factor for ANOVA.

## Warning: Converting "Campagne" to factor for ANOVA.

## Warning: Data is unbalanced (unequal N per group). Make sure you specified a
## well-considered value for the type argument to ezANOVA().

## Coefficient covariances computed by hccm()

## $ANOVA
##        Effect DFn DFd      SSn       SSd         F            p p<.05       ges
## 1 (Intercept)   1 103 66683871 140954831 48.727941 2.944440e-10     * 0.3211534
## 2    Campagne   5 103 44212390 140954831  6.461469 2.908251e-05     * 0.2387701
## 
## $`Levene's Test for Homogeneity of Variance`
##   DFn DFd      SSn       SSd        F          p p<.05
## 1   5 103 15177066 127757447 2.447196 0.03871207     *

PostHocTest(aov(X~Campagne, data=Database), 
            which="Campagne", method="hsd", conf.level=.95) 

## 
##   Posthoc multiple comparisons of means : Tukey HSD 
##     95% family-wise confidence level
## 
## $Campagne
##                diff     lwr.ci    upr.ci    pval    
## BF3-BF2   233.13656  -970.3307 1436.6039 0.99318    
## CA1-BF2   479.98524  -760.5205 1720.4909 0.87050    
## CA2-BF2  1962.83219   789.4313 3136.2331 6.1e-05 ***
## CA3-BF2   462.15050  -778.3552 1702.6562 0.88754    
## CA4-BF2   433.53266  -653.4910 1520.5563 0.85539    
## CA1-BF3   246.84867  -956.6186 1450.3160 0.99112    
## CA2-BF3  1729.69563   595.5224 2863.8688 0.00033 ***
## CA3-BF3   229.01393  -974.4534 1432.4812 0.99372    
## CA4-BF3   200.39609  -844.1610 1244.9532 0.99348    
## CA2-CA1  1482.84696   309.4460 2656.2479 0.00506 ** 
## CA3-CA1   -17.83474 -1258.3404 1222.6710 1.00000    
## CA4-CA1   -46.45258 -1133.4762 1040.5710 1.00000    
## CA3-CA2 -1500.68170 -2674.0826 -327.2808 0.00437 ** 
## CA4-CA2 -1529.29954 -2539.0696 -519.5294 0.00038 ***
## CA4-CA3   -28.61784 -1115.6415 1058.4058 1.00000    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#install.packages("ez")
#install.packages("DescTools")
library(ez)
library(DescTools)

Database <- data_clean
Database <- cbind(Database, ACP_3factVAR$scores)
Database <- cbind(Database, ID=c(1:nrow(Database)))
#Database <- as.matrix(Database)
#view(Database)
#Database$ID.1 <- NULL

ezANOVA(data = Database,
        dv = .(X),
        wid = .(ID),
        between = .(TYPE),
        detailed = TRUE, type = 3)

## Warning: Converting "ID" to factor for ANOVA.

## Warning: Converting "TYPE" to factor for ANOVA.

## Warning: Data is unbalanced (unequal N per group). Make sure you specified a
## well-considered value for the type argument to ezANOVA().

## Coefficient covariances computed by hccm()

## $ANOVA
##        Effect DFn DFd      SSn       SSd         F            p p<.05       ges
## 1 (Intercept)   1 105 71885754 162894305 46.336820 6.366546e-10     * 0.3061834
## 2        TYPE   3 105 22272917 162894305  4.785631 3.637463e-03     * 0.1202854
## 
## $`Levene's Test for Homogeneity of Variance`
##   DFn DFd      SSn       SSd       F           p p<.05
## 1   3 105 20716638 146579302 4.94669 0.002979762     *

PostHocTest(aov(X~TYPE, data=Database), 
            which="TYPE", method="hsd", conf.level=.95) 

## 
##   Posthoc multiple comparisons of means : Tukey HSD 
##     95% family-wise confidence level
## 
## $TYPE
##                                  diff     lwr.ci    upr.ci   pval    
## rural-compostage            -104.6715 -1283.9118 1074.5689 0.9956    
## sourceindustrie-compostage  1929.4623   215.6771 3643.2475 0.0208 *  
## urbain-compostage            143.6931 -1043.6522 1331.0383 0.9890    
## sourceindustrie-rural       2034.1337   627.7128 3440.5547 0.0015 ** 
## urbain-rural                 248.3645  -423.0137  919.7428 0.7691    
## urbain-sourceindustrie     -1785.7692 -3198.9927 -372.5457 0.0072 ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

On voit ici quelque chose de trés interessant car l’ACP ne montre pas grand chose pour les location ou les TYPE, mais ici on voit que les differences significatives sont TOUTES avec sourceindustrie dedans et sourceindustrie et dans TOUT les tests significatifs. Donc il semble bien que l’industrie ait un impacte non negligeable sur son entourage en terme d’émission de composant.