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.
# 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
}
## 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.
# 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.
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
| 54.31 |
32.31 |
52.14 |
192 |
112.7 |
167.6 |
Table continues below
| 222.3 |
465 |
150.2 |
826.2 |
1309 |
575.1 |
Table continues below
| 100.3 |
101.5 |
74.46 |
254.3 |
Table continues below
| 517.7 |
154 |
834.9 |
1325 |
Table continues below
| 215.3 |
1652 |
1984 |
1104 |
Table continues below
| 266.7 |
538.2 |
125.5 |
546.7 |
Table continues below
| 2075 |
302.7 |
491.9 |
589.5 |
Table continues below
| 305.1 |
352.7 |
717.8 |
Table continues below
| 158.2 |
793.2 |
808.1 |
Table continues below
| 473.5 |
629.6 |
825.3 |
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
| BF2 |
15.55 |
3.441 |
14.86 |
97.98 |
33.35 |
| BF3 |
19.46 |
3.39 |
19.45 |
172.9 |
91.1 |
| CA1 |
42.2 |
9.459 |
38.91 |
218.9 |
74.67 |
| CA2 |
105.3 |
15.68 |
107.4 |
313.9 |
144.8 |
| CA3 |
53.43 |
5.744 |
52.12 |
261.2 |
63.32 |
| CA4 |
68.58 |
14.55 |
69.79 |
119.9 |
48.51 |
Table continues below
| 81.49 |
62.26 |
22.42 |
56.02 |
206.6 |
67.8 |
| 154.1 |
160.9 |
58.71 |
135.9 |
439.8 |
110.5 |
| 197.1 |
175.2 |
77.68 |
161.4 |
686.6 |
223.1 |
| 262.2 |
587.2 |
1046 |
306.2 |
2169 |
2745 |
| 255.2 |
182.2 |
83.23 |
174.5 |
668.8 |
251.9 |
| 112.8 |
144.6 |
90.45 |
130 |
640.2 |
359.3 |
Table continues below
| 190.4 |
21.4 |
5.534 |
22.18 |
34.99 |
| 409.9 |
53.41 |
34.94 |
51.1 |
106.9 |
| 634.5 |
79.8 |
23.66 |
72.72 |
133.5 |
| 1336 |
233 |
133 |
182.4 |
728.6 |
| 614.6 |
71.77 |
17.29 |
73.32 |
262.6 |
| 580.6 |
107.4 |
101.3 |
86.53 |
199.7 |
Table continues below
| 12.74 |
35.75 |
21.92 |
13.37 |
| 97.71 |
79.46 |
59.93 |
28.84 |
| 37.59 |
122.4 |
147.2 |
71.97 |
| 1127 |
480.3 |
3606 |
559.7 |
| 63.51 |
245 |
219.8 |
90.32 |
| 84.8 |
198 |
558.3 |
324 |
Table continues below
| 17.15 |
48.19 |
36.47 |
30.59 |
| 53.31 |
31.78 |
18.89 |
28.34 |
| 139.6 |
433.5 |
206.9 |
427.5 |
| 3560 |
5572 |
874.7 |
5331 |
| 215.3 |
1245 |
445 |
1355 |
| 512.2 |
1705 |
542.2 |
1687 |
Table continues below
| 21.2 |
18.17 |
11.38 |
113.8 |
| 75.79 |
69.11 |
49.74 |
317.3 |
| 133.1 |
49.15 |
131.2 |
328 |
| 855.1 |
1091 |
584.1 |
1749 |
| 180.3 |
143.9 |
123.9 |
335.9 |
| 232.6 |
186.1 |
164.7 |
332.1 |
Table continues below
| 57.35 |
87.08 |
30.82 |
4.816 |
| 176.8 |
269.9 |
55.65 |
9.739 |
| 95.82 |
293.5 |
396.8 |
317.2 |
| 4885 |
533.8 |
986.9 |
158 |
| 166.4 |
322.3 |
180.5 |
98.02 |
| 149.4 |
320.6 |
885.5 |
828.4 |
Table continues below
| 32.75 |
11.59 |
4.965 |
| 53.67 |
60.19 |
83.55 |
| 326.9 |
840.3 |
1563 |
| 963.8 |
497.9 |
172.9 |
| 191.6 |
237.1 |
96.5 |
| 580.3 |
415.2 |
699 |
Table continues below
| 12.26 |
61.35 |
43.94 |
| 39.01 |
60.92 |
33.38 |
| 184.6 |
463 |
283.4 |
| 397.3 |
1278 |
482.3 |
| 196.1 |
244.2 |
124.4 |
| 156.6 |
1772 |
619.6 |
Table continues below
| 48.22 |
77.2 |
53.02 |
| 55.35 |
147.7 |
75.86 |
| 473.5 |
380 |
264.5 |
| 1290 |
2125 |
954 |
| 211.2 |
275.9 |
120 |
| 1952 |
526.4 |
253.7 |
| 76.09 |
| 145 |
| 337.4 |
| 2169 |
| 313.3 |
| 461.1 |
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
| hiver |
32.23 |
16.81 |
25.35 |
187.3 |
90.78 |
| été |
83.43 |
23.51 |
78.35 |
198.3 |
137.3 |
Table continues below
| 170.2 |
145.6 |
79.45 |
135 |
498.5 |
260.7 |
| 150.7 |
323.5 |
693.5 |
203.8 |
1258 |
1897 |
Table continues below
| 468.7 |
56.49 |
31.95 |
55.11 |
133.6 |
| 830.6 |
158.2 |
129.7 |
119.8 |
413.5 |
Table continues below
| 102.3 |
108.3 |
110.5 |
96.08 |
| 754.9 |
278.2 |
1790 |
1572 |
Table continues below
| 70.93 |
426.5 |
547 |
105.6 |
| 650.9 |
3268 |
2038 |
1935 |
Table continues below
| 101.7 |
100.9 |
72.86 |
275.1 |
| 484.3 |
762.6 |
340.3 |
904.8 |
Table continues below
| 161 |
254.7 |
162.4 |
215.3 |
| 3138 |
386.6 |
926.5 |
644.3 |
Table continues below
| 58.49 |
280 |
820.6 |
| 840.2 |
448.6 |
547.9 |
Table continues below
| 80.98 |
202.7 |
224.3 |
| 333.4 |
1572 |
613.5 |
Table continues below
| 110.4 |
217.9 |
187.2 |
| 1492 |
1173 |
1011 |
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
| P01 |
54.84 |
182.3 |
161.9 |
679.4 |
103.1 |
| P02 |
63.49 |
281.6 |
193.9 |
743 |
142.6 |
| P03 |
54.69 |
167.1 |
157.2 |
626.1 |
97.86 |
| P04 |
88.95 |
188.2 |
218.9 |
991.3 |
170.7 |
| P05 |
54.7 |
168.9 |
177.4 |
739.2 |
89.62 |
| P06 |
45.4 |
129.9 |
109.9 |
435.3 |
59.11 |
| P07 |
52.41 |
149.8 |
126.8 |
497.5 |
80.58 |
| P08 |
37.52 |
145.3 |
155.4 |
545.5 |
56.17 |
| P09 |
42.8 |
177.8 |
146.3 |
574.9 |
58.64 |
| P10 |
59.28 |
413.4 |
2527 |
6762 |
380.5 |
| P11 |
40.54 |
214.7 |
234.4 |
866.4 |
72.04 |
| P12 |
56.88 |
176.4 |
156.2 |
600.7 |
84.32 |
| P13 |
36.32 |
148.1 |
121.1 |
438.7 |
50.61 |
| P14 |
48.79 |
147 |
155.2 |
576.4 |
68.4 |
| P15 |
51.36 |
353.9 |
386 |
1372 |
101.3 |
| P16 |
45.66 |
172.2 |
163.4 |
702.3 |
78.95 |
| P17 |
39.92 |
243.1 |
174.2 |
700.8 |
69.25 |
| P18 |
49.48 |
154 |
164.2 |
600.2 |
78.99 |
| P19 |
66.59 |
182.6 |
141.4 |
621.1 |
95.97 |
| P20 |
17.12 |
103.3 |
142.8 |
287.7 |
225.6 |
| P21 |
49.68 |
160.1 |
121.9 |
496 |
91.14 |
| P25 |
60.33 |
174.7 |
143.4 |
582.3 |
75.8 |
| P26 |
71.04 |
104.5 |
107.3 |
512 |
56.74 |
| P27 |
76.99 |
170 |
140.6 |
588.9 |
81.52 |
| P28 |
73.11 |
238.2 |
213.7 |
929 |
94.11 |
| P29 |
87.33 |
210.4 |
288.3 |
1256 |
156.7 |
| P30 |
71.9 |
141.5 |
151.6 |
717.6 |
102.9 |
| P32 |
88.18 |
249.9 |
261.3 |
1290 |
250.2 |
| P33 |
98.1 |
327.8 |
475.2 |
2369 |
197.8 |
Table continues below
| 227.2 |
1081 |
1916 |
282.6 |
| 338.6 |
1250 |
2104 |
227.7 |
| 181.3 |
975.8 |
1538 |
174.8 |
| 268.7 |
2267 |
2584 |
369.2 |
| 217.5 |
839.5 |
1868 |
205.7 |
| 186.6 |
1001 |
1609 |
282.8 |
| 172.4 |
678.3 |
1419 |
126.5 |
| 88.65 |
120.6 |
397.2 |
89.87 |
| 121.5 |
231.2 |
797.9 |
217.5 |
| 2718 |
1845 |
2520 |
2754 |
| 210.7 |
246.7 |
919.6 |
155.1 |
| 204.6 |
875.1 |
1640 |
157 |
| 103.8 |
86.3 |
664.7 |
85.43 |
| 156.7 |
676.3 |
1455 |
163.3 |
| 250.3 |
815.6 |
1691 |
342.2 |
| 212.3 |
208 |
764.3 |
118.4 |
| 166.8 |
308 |
914.7 |
149 |
| 240 |
737.1 |
1506 |
261.1 |
| 199.4 |
900.6 |
1669 |
185.3 |
| 78 |
166.2 |
529.4 |
12.3 |
| 213.4 |
730.3 |
1476 |
146.5 |
| 234.3 |
345.3 |
1617 |
325.1 |
| 141.3 |
471.5 |
1449 |
101.2 |
| 294.7 |
1731 |
3026 |
203.5 |
| 180.5 |
1221 |
2691 |
364.7 |
| 397.8 |
1795 |
3217 |
865.3 |
| 132.5 |
1941 |
3759 |
265.7 |
| 460.5 |
1562 |
3351 |
506.6 |
| 536.2 |
2263 |
4176 |
523.6 |
Table continues below
| 331.3 |
428.1 |
341.9 |
818.8 |
| 442.5 |
875.6 |
946.3 |
705.8 |
| 262.3 |
467 |
216.3 |
1046 |
| 356.1 |
456.3 |
880.8 |
1101 |
| 386.3 |
949.3 |
244.1 |
703.1 |
| 326.4 |
631 |
718.1 |
772.5 |
| 236.4 |
296.7 |
162.2 |
706.1 |
| 220 |
219.9 |
74.28 |
550.1 |
| 350.5 |
274.9 |
119 |
414.4 |
| 11327 |
510.2 |
204.8 |
410.7 |
| 371.1 |
178.1 |
94.73 |
605.2 |
| 351.2 |
428.4 |
176.5 |
757.1 |
| 202.4 |
280.8 |
137.5 |
553.4 |
| 238.5 |
387 |
181.2 |
702.6 |
| 487.6 |
333.5 |
161.3 |
712.8 |
| 360.3 |
337.3 |
133.2 |
549.6 |
| 323.3 |
379.2 |
1480 |
687 |
| 489.7 |
346.7 |
891.7 |
924.2 |
| 344.9 |
435.5 |
164.1 |
882.3 |
| 170.8 |
430.3 |
146.1 |
227 |
| 277.5 |
376 |
135.7 |
820 |
| 321 |
370.1 |
254.7 |
1159 |
| 212.3 |
643 |
201.3 |
1722 |
| 266.3 |
1422 |
1251 |
1066 |
| 266.3 |
497.6 |
359.4 |
593.2 |
| 705.7 |
764.2 |
310.2 |
1291 |
| 407 |
892.2 |
297.8 |
1814 |
| 786.8 |
1235 |
383.6 |
2183 |
| 764.4 |
640.7 |
271.6 |
921.2 |
| 359.1 |
| 529.9 |
| 1240 |
| 1784 |
| 368.9 |
| 593.8 |
| 665.6 |
| 202.1 |
| 226.5 |
| 810.2 |
| 241.6 |
| 803.4 |
| 296.2 |
| 476.7 |
| 675.5 |
| 302.8 |
| 280.6 |
| 131.7 |
| 786.8 |
| 435.8 |
| 666.6 |
| 473.6 |
| 278.3 |
| 1166 |
| 638 |
| 1492 |
| 1082 |
| 2185 |
| 981.2 |
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
| compostage |
58.68 |
226.4 |
176.1 |
707.7 |
| rural |
56.54 |
160.5 |
149.9 |
603 |
| sourceindustrie |
48.3 |
256.3 |
940 |
2637 |
| urbain |
51.81 |
210.9 |
214.7 |
851.4 |
Table continues below
| 120.7 |
276.7 |
1156 |
1999 |
| 84.38 |
186.5 |
917.9 |
1686 |
| 165.9 |
987 |
769.2 |
1372 |
| 104.9 |
225.9 |
689.1 |
1582 |
Table continues below
| 258.2 |
380.7 |
627 |
610.5 |
| 206.1 |
300.6 |
502.7 |
329.1 |
| 1063 |
4010 |
353.3 |
147.6 |
| 228.2 |
386 |
471.6 |
354.1 |
| 768.6 |
435 |
| 863.2 |
783.4 |
| 413.2 |
421.1 |
| 772.6 |
528.8 |
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
| apres |
69.55 |
25.98 |
66.07 |
214.6 |
117.3 |
| avant |
17.63 |
3.902 |
17.48 |
137.8 |
78.93 |
Table continues below
| 201.6 |
267.1 |
546.4 |
174.5 |
1032 |
1510 |
| 126.5 |
114.7 |
67.12 |
113.8 |
330.5 |
149.5 |
Table continues below
| 715.3 |
126.1 |
109.5 |
88.42 |
329.5 |
| 329.5 |
38.4 |
30.12 |
29.77 |
73.18 |
Table continues below
| 599 |
229.3 |
1164 |
1456 |
| 79.56 |
51.59 |
42.11 |
29.69 |
Table continues below
| 453.5 |
2322 |
2011 |
1533 |
| 40.27 |
39.47 |
29.22 |
30.15 |
Table continues below
| 356.6 |
618.3 |
203.4 |
681.6 |
| 50.2 |
58.14 |
36.33 |
221.9 |
Table continues below
| 2459 |
347.2 |
678 |
611.7 |
| 168.1 |
188.6 |
44.01 |
14.76 |
Table continues below
| 555.7 |
483.7 |
819.3 |
| 45.33 |
37.41 |
64.97 |
Table continues below
| 234.6 |
1097 |
779.7 |
| 23.71 |
61.12 |
38.04 |
Table continues below
| 1046 |
843.6 |
898.7 |
| 55.1 |
114.7 |
74.29 |
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 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
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 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)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:Hmisc':
##
## subplot
## 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
# 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
## 1
# 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 a 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, tres 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.
