Ce document a pour but d'introduire les praticiens et scientifiques en science du sol à l'analyse de données avec R. Il s'agit d'une création communautaire évolutive.
Laissons le New York Times, chez qui on utilise R, en faire l’éloge : « R est […] le nom d’un langage de programmation utilisé par un nombre grandissant d’analystes de données dans les milieux corporatifs et académiques. Il devient la lingua franca notamment puisque l’exploration de données est à son âge d’or, que ce soit pour déterminer le prix d’une publicité, découvrir de nouveaux médicaments ou ajuster des modèles financiers.[…] Mais R a rapidement gagné en popularité parce que les statisticiens, les ingénieurs et les scientifiques n’ayant pas de compétences en programmation le trouvent facile à utiliser. »
L’environnement d’analyse de données R est un langage de programmation à source libre (donc gratuit) particulièrement adapté pour effectuer des calculs statistiques et générer des graphiques.
Un tableur est un tableau d’informations sur lesquelles il est possible d’effectuer des calculs. Les logiciels tableurs les plus utilisé sont Excel, distribué avec la suite Microsoft Office, et Calc, distribué gratuitement avec la suite LibreOffice. Les tableurs sont largement utilisés pour effectuer des calculs grâce à leur facilité à prendre en main. Toutefois, ils deviennent rapidement lourds à utiliser lorsque les calculs deviennent complexes où sont effectués sur un grand nombre de données. À la suite d’audits effectués en Europe, le European Spreadsheet Risks Interest Group (EuSpRI) a découvert que 90% des tableurs contenaient de sérieuses erreurs de calcul (site internet de EuSpRI).
“Asok, selon mon tableur, tu as fait un travail épouvantable. | Peut-être que votre tableur est mal conçu et qu’il ne capte pas la complexité du monde réel. | Et n’oublions pas la quasi certitude que votre formule pointe sur la mauvaise cellule. – Les nombres ne mentent pas.” crédit: dilbert.com
“Le MBA versus la vielle sorcière folle. – Je ne sais qui croire. | Les tableurs ne mentent pas. Mais les excréments de chauve-souris non plus. | Rappelez moi encore qui a ruiné l’économie. Étaient-ce les sorcières?” crédit: dilbert.com
Les tableurs sont néanmoins d’une grande utilité, notamment pour présenter esthétiquement et de manière claire de l’information tabulée : petits budgets et soumissions, brefs bilans hydrauliques, brève planification quinquennale, etc. Ils sont aussi utiles pour effectuer des opérations simples et rapides sur un petit nombre de données : sommations, moyennes, écarts-type, filtre et classement de données et graphiques.
Mais ils sont à utiliser avec précaution pour traiter un grand nombre de données ou effectuer des opérations statistiques complexes. Dans ce cas, R est préférable.
R fonctionne avec des lignes de commande. Sa courbe d’apprentissage est plutôt abrupte, mais une fois que ses principes de base sont maîtrisés (disons, après un ou deux mois d’utilisation), l’analyste s’ouvre à un large éventail de possibilités : ANOVA, tests de normalité, analyse en composante principale, partitionnement de données, analyse compositionnelle, méta-analyses, création de graphiques de grande qualité, etc.
Pour interagir avec R, la méthode classique et peu pratique est l’utilisation d’un interprétateur de commande (shell).
Note: “La liberté c'est la liberté de dire 2+2=4” George Orwell
Sur Windows, RGUI s’installe automatiquement avec R. RGUI permet d’utiliser une feuille de calcul au lieu d’entrer une commande à la fois. Toutefois, l’interface est plutôt rigoureuse.
Il est suggéré ici d’utiliser RStudio, une interface libre pour R développée récemment, téléchargeable pour Linux, Windows et MacOS à rstudio.org.
RStudio offre une aire pour la feuille de calcul (haut-gauche, avec coloration syntaxique), une aire pour la console (bas-gauche), une aire présentant les variables que vous avez créé (haut-droit) et une aire pour les graphiques (bas-droit). Vous observerez que d’autres onglets sont disponibles: navigation dans vos fichiers, gestionnaire de modules et interface d’aide.
Aux librairies inclues dans le module de base, des milliers d’autres peuvent être installées gratuitement depuis internet. Par exemple, pour les calculs écologiques, vous apprécierez les modules vegan et ade4. Pour l’analyse compositionnelle, vous utiliserez le module compositions. Le module nortest vous offrira des tests de normalité et le module outliers permettra de détecter les valeurs aberrantes. Les pédologues pourraient utiliser le module aqp pour, entre autres, la classification des taxons. Pour les méta-analyses, vous choisirez entre meta et metafor. Pour effectuer des graphiques avancés de grande qualité, vous vous intéresserez à ggplot2. Pour la conception expérimentale en agriculture, un étudiant péruvien a mis à la disposition de quiconque le module agricolae. Les modules peuvent être installés et mis à jour par ligne de commande ou via RStudio (onglet Packages), à condition de posséder une connexion internet. Nous verrons cela plus tard.
La gratuité et l’ouverture font de R un outil de plus en plus utilisé. Des centaines de documents pour débutants ou avancés ont été mis disponibles sur le web par de généreux contributeurs. On répond à la plupart de vos questions en quelques secondes en interrogeant Google, ou pour obtenir un résultat probablement plus précis, une adaptation de Google spécialement conçue pour R: RSeek.org. Par exemple, si vous vous demandez comment effectuer une analyse en composantes principales, dans Google, cherchez R pca ou R acp et vous trouverez des documents PDF et des forums de discussion qui vous guideront dans vos analyses. R introduction pourrait aussi être intéressant. De telles recherches sont aussi utiles pour vérifier si un module de R a été créé pour effectuer une fonction qui vous intéresse, par exemple : R discriminant analysis.
Pour découvrir davantage les possibilités de R, visitez RPubs, un site de publication de documents R, ou R-Bloggers, un site aggrégateur de plusieurs blogues à propos de R, de son développement et de ses capacités.
Robin Hankin: I'd say that without a tool like R you cannot learn statistics.
David Whiting: I believe Fisher and a few others managed to get by without it.
Peter Dalgaard: But think how far they could have got with R!
– Robin Hankin, David Whiting and Peter Dalgaard (on teaching/learning statistics with R)
R-help (December 2004)
2 + 2
## [1] 4
a <- c(1, 2, 3, 4, 5)
b <- c(2, 4, 5, 6, 1)
a + b
## [1] 3 6 8 10 6
R est un langage de programmation orienté objet. Il vous permet d’assigner des fonctions, des tableaux, des matrices, des vecteurs, des listes, des variables, etc. à des variables représentées par des chaînes de caractères. Les fonctions que vous utilisez sont des chaînes de caractères auxquelles une série de commandes ont déjà été assignées, comme la fonction length qui retourne la longueur d’un vecteur:
length(a)
## [1] 5
Un vecteur peut s’écrire de différentes façons:
c(2, 4, 5, "janvier", TRUE)
## [1] "2" "4" "5" "janvier" "TRUE"
1:10
## [1] 1 2 3 4 5 6 7 8 9 10
seq(from = -3, to = 1, by = 0.5)
## [1] -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0
rep("carotte", times = 3)
## [1] "carotte" "carotte" "carotte"
Une matrice est un tableau à plus d’une ligne, dont tous les éléments sont du même type. En effet, il y a différents types d’éléments: chiffres (int ou numeric), caractères (character), date, facteur, logique, etc.
Si votre tableau contient des colonnes de différentes natures, par exemple, une colonne de date, une colonne de cultivar (facteur), et des colonnes de concentrations d’éléments, on parlera d’un dataframe. Il s’agit du format que R attribuera sans doute aux tableaux que vous importerez.
Une liste est une série d’objets. On peut y insérer des matrices, des dataframes, des variables, etc. Les fonctions utilisent parfois des types particuliers. Pour obtenir le type, on utilise la fonction class.
a <- c(1, 2, 3, 4)
b <- "ananas"
c <- data.frame(col1 = c(1, 5, 7, 9), col2 = c("a", "b", "c", "d"), col3 = rep("poisson",
4))
maListe <- list()
maListe[[1]] <- a
maListe[[2]] <- b
maListe[[3]] <- c
maListe
## [[1]]
## [1] 1 2 3 4
##
## [[2]]
## [1] "ananas"
##
## [[3]]
## col1 col2 col3
## 1 1 a poisson
## 2 5 b poisson
## 3 7 c poisson
## 4 9 d poisson
class(maListe)
## [1] "list"
class(maListe[[1]])
## [1] "numeric"
class(maListe[[2]])
## [1] "character"
class(maListe[[3]])
## [1] "data.frame"
Notez que le = peut être utilisé en tout temps à la place du <-, mais il est préférable d’un point de vue de lecture du code de réserver la flèche <- pour l’égalité et le = pour attribuer des objets aux arguments d’une fonction.
Avant d’importer des données, vous devrez vous assurer que votre tableau est bien arrangé. Dans Excel, vous pouvez fusionner des cellules pour attribuer un titre (ou une unité de mesure) à une série de colonne. Vous pouvez aussi avoir plusieurs tableaux dans une même feuille. Dans R, une colonne ne peut pas avoir plus d’une seule entête. Par exemple, dans Excel, vous pouvez travailler avec ces tableaux.
| Culture | Annee | N.per | P.per | K.per | Ca.per | Mg.per |
|---|---|---|---|---|---|---|
| PanetRave | 2007 | 3,13 | 0,36 | 5,33 | 1,96 | 0,36 |
| PanetRave | 2008 | 2,94 | 0,36 | 4,91 | 1,73 | 0,32 |
| PanetRave | 2009 | 3,05 | 0,36 | 5,19 | 1,81 | 0,42 |
| ChouMarseille | 2007 | 5,47 | 0,68 | 4,38 | 0,90 | 0,22 |
| ChouMarseille | 2008 | 5,46 | 0,68 | 4,29 | 0,85 | 0,19 |
| ChouMarseille | 2009 | 5,47 | 0,68 | 4,29 | 0,87 | 0,20 |
| PommeMer | 2007 | 3,00 | 0,47 | 4,12 | 0,65 | 0,26 |
| PommeMer | 2008 | 3,04 | 0,49 | 4,11 | 0,69 | 0,27 |
| PommeMer | 2009 | 3,03 | 0,50 | 4,06 | 0,70 | 0,26 |
| Courgine | 2007 | 4,25 | 0,27 | 5,77 | 1,80 | 0,30 |
| Courgine | 2008 | 4,30 | 0,33 | 5,79 | 1,74 | 0,30 |
| Courgine | 2009 | 4,50 | 0,50 | 6,01 | 1,60 | 0,30 |
Il est préférable d’enregistrer le tableau en format csv (pour comma separated values, même si le séparateur de colonne par défaut dans Excel est un point-virgule). Ce format de garde aucun formatage, ni aucun calcul. Il s’agit d’un fichier texte à sa plus simple expression. Si vous ouvrez le fichier csv dans un bloc note, vous aurez :
Culture;Annee;N.per;P.per;K.per;Ca.per;Mg.per
PanetRave;2007;3,13;0,36;5,33;1,96;0,36
PanetRave;2008;2,94;0,36;4,91;1,73;0,32
PanetRave;2009;3,05;0,36;5,19;1,81;0,42
ChouMarseille;2007;5,47;0,68;4,38;0,90;0,22
ChouMarseille;2008;5,46;0,68;4,29;0,85;0,19
ChouMarseille;2009;5,47;0,68;4,29;0,87;0,20
PommeMer;2007;3,00;0,47;4,12;0,65;0,26
PommeMer;2008;3,04;0,49;4,11;0,69;0,27
PommeMer;2009;3,03;0,50;4,06;0,70;0,26
Courgine;2007;4,25;0,27;5,77;1,80;0,30
Courgine;2008;4,30;0,33;5,79;1,74;0,30
Courgine;2009;4,50;0,50;6,01;1,60;0,30
Pour importer le tableau dans R, vous utiliserez la fonction setwd pour indiquer votre dossier de travail, puis read.csv pour charger le fichier.
# setwd('C://Users//MonNom//Desktop//donnees')# sous Windows
# setwd('/home/monnom/Bureua/donnees') # sous Linux
data <- read.csv("/home/serge-etienne/Documents/etudes/ulaval/RPubs/rSciencesSol/legumes_fictifs.csv",
header = TRUE, sep = ";", dec = ",")
Dans la fonction, j’ai ajouté les arguments header=TRUE pour signifier que la première ligne est le titre de la colonne, non pas une donnée, sep=';' pour signifier que le séparateur de colonne est un point-virgule et dec=',' pour signifier que la décimale est une virgule. Pour exécuter des lignes dans RStudio, sélectionnez les lignes désirées et appuyez sur Run (ou Ctrl + Enter). Visionnez maintenant vos données dans RStudio en cliquant sur le titre approprié dans la fenêtre Workspace.
Pour appeler une cellule du tableau :
data[6, 4]
## [1] 0.38
Pour appeler les lignes 1 à 5 sur une ou plusieurs colonnes:
data[1:5, 2]
## [1] 2007 2007 2007 2007 2007
data$Annee[1:5]
## [1] 2007 2007 2007 2007 2007
data[1:5, c(1, 3, 5)]
## Culture N K
## 1 PanetRave 3.13 5.33
## 2 PanetRave 3.03 4.94
## 3 PanetRave 3.01 5.08
## 4 PanetRave 3.53 5.44
## 5 PanetRave 3.96 6.10
data[1:5, -c(1, 2)]
## N P K Ca Mg
## 1 3.13 0.36 5.33 1.96 0.36
## 2 3.03 0.37 4.94 2.33 0.40
## 3 3.01 0.36 5.08 1.76 0.36
## 4 3.53 0.36 5.44 2.13 0.37
## 5 3.96 0.38 6.10 2.76 0.50
Pour sélectionner seulement les données des pommes:
data[data$Espece == "PommeMer", ] # le double egale signifie une operation booleenne
## [1] Culture Annee N P K Ca Mg
## <0 rows> (or 0-length row.names)
Et effectuer des statistiques de base:
apply(data[, -c(1, 2)], 2, mean) # données data[,-c(1,2)], 2: effectuer par colonne (1: par lignes), la fonction mean
## N P K Ca Mg
## 4.6131 0.5614 5.1443 1.4418 0.2790
apply(data[, -c(1, 2)], 2, sd) # idem pour sd, fonction pour calculer l'écart-type
## N P K Ca Mg
## 1.8583 0.2287 1.0608 0.6547 0.1097
colMeans(data[, -c(1, 2)]) # moyennes par colonne, commande raccourcie
## N P K Ca Mg
## 4.6131 0.5614 5.1443 1.4418 0.2790
Note 1: effectuer des moyennes et des écarts-type sur des compositions n’a aucun sens.
Note 2: le # signifie que le reste de la ligne est un commentaire que R n’interprétera pas.
Pour installer puis charger un paquet spécifique :
# install.packages('compositions', dep=TRUE) # dep=TRUE pour installer les
# dependances
library(compositions)
## Loading required package: rgl
## Loading required package: tensorA
## Attaching package: 'tensorA'
## The following object(s) are masked from 'package:base':
##
## norm
## Loading required package: robustbase
## Loading required package: energy
## Loading required package: boot
## Attaching package: 'boot'
## The following object(s) are masked from 'package:robustbase':
##
## salinity
## Loading required package: MASS
## Welcome to compositions, a package for compositional data analysis. Find
## an intro with "? compositions"
## Attaching package: 'compositions'
## The following object(s) are masked from 'package:stats':
##
## cor, cov, dist, var
## The following object(s) are masked from 'package:base':
##
## %*%, scale, scale.default
Pour effectuer un graphique, plot est une fonction générique qui appellera le graphique selon le type de données.
class(data[, 3:5])
## [1] "data.frame"
plot(data[, 3:5])
composition <- acomp(data[, 3:5])
class(composition)
## [1] "acomp"
plot(composition) # lance en fait la fonction plot.acomp
plot(composition, pch = 2) # pour changer le type de point
Les boucles et les conditions sont parfois très pratiques pour éviter d’écrire les mêmes lignes de code en répétition.
Ou pour générer des graphiques
par(mfrow = c(2, 2)) # les graphiques apparaîtront dans une matrice 2 par 2
for (i in 1:length(levels(data$Culture))) {
plot(acomp(data[data$Culture == levels(data$Culture)[i], -c(1, 2, 6, 7)]),
pch = 21, bg = i)
}
L’interface de RStudio vous permet d’exporter vos graphiques en formats bitmap ou vectoriel. Pour l’édition sur le pouce, cliquez sur Export > Copy to clipboard, copy as Meta File, ou enregistrer en format SVG. Collez le graphique dans Inkscape (logiciel libre d’édition vectorielle), ou ouvrez le fichier SVG.
Pour terminer, notez qu’il est tout à fait normal d’obtenir des erreurs lorsque l’on effectue des opérations. Lisez bien ce que R vous indique comme erreur. À force d’utiliser R, vous éviterez des erreurs banales et saurez reconnaître rapidement la source des erreurs pour les corriger rapidement. Vous saurez aussi comment cherche dans le moteur d’aide ou sur internet pour trouver réponse à vos questions.
On suppose que vous connaissez les bases de l’analyse compositionnelle et leurs transformations, ainsi que leur importance en sciences du sol. John Aitchison (2003) présente l’analyse compositionnelle dans le document PDF A Concise Guide to Compositional Data Analysis.
Le document ‘‘compositions’’: A unified R package to analyze compositional data, de K. Gerald van den Boogaart et Raimon Tolosana-Delgado (2008), est une ressource qui facilite la prise en main des fonctions offertes par le module compositions.
Voici une liste des fonctions les plus courantes.
On vous présente ici un exemple d’utilisation du module compositions, inspiré du document de van den Boogaart et Tolosana-Delgado (2008). Nous traiterons de données hydrogéochimiques collectées dans différents cours d’eau.
On suppose que vous avez déjà installé le module compositions.
library("compositions") # charger les fonctions du module compositions
data(Hydrochem) # charger des donnees de R
Donnees <- Hydrochem[, c("Cl", "SO4", "HCO3")] # selectionner seulement les colonnes a traiter
Composition <- acomp(Donnees) # fermer le simplex et declarer que les donnees sont compositionnelles
plot(Composition, col = Hydrochem$River) # diagramme ternaire, couleur du point selon la riviere
isoPortionLines(by = 0.2, lwd = 0.2) # grille
isoProportionLines(by = 0.5, lwd = 0.5, labs = F) # lignes d’isoproportion
Vous obtiendrez un diagramme ternaire :
image008.png
Nous avons déclaré le tableau Composition comme étant un tableau de données compositionnelles :
class(Composition)
## [1] "acomp"
Pour calculer des statistiques sur les données, nous pouvons utiliser des appels de fonctions génériques. R reconnaîtra que ces données sont compositionnelles et effectuera les calculs statistiques en conséquence.
mean(Composition) # moyenne
## Cl SO4 HCO3
## 0.2262 0.3267 0.4471
## attr(,"class")
## [1] acomp
cov(Composition) # matrice de variance-covariance
## Cl SO4 HCO3
## Cl 0.5394 -0.19776 -0.34165
## SO4 -0.1978 0.16205 0.03571
## HCO3 -0.3416 0.03571 0.30594
cor(Composition) # matrice de correlations
## Cl SO4 HCO3
## Cl 1.0000 -0.6689 -0.8410
## SO4 -0.6689 1.0000 0.1604
## HCO3 -0.8410 0.1604 1.0000
Pour effectuer les statistiques sur une rivière en particulier, il faut déclarer de nouveau la composition avec la fonction acomp (le filtre retire la classe) :
mean(acomp(Composition[Hydrochem$River == "Cardener", ]))
## Cl SO4 HCO3
## 0.3027 0.2650 0.4324
## attr(,"class")
## [1] acomp
Si l’on travaille avec des ratios de logarithmes additifs, on peut calculer la transformation, effectuer les calculs statistiques désirés, puis inverser la transformation.
transf <- alr(Composition) # transformation alr
moyenne <- mean(transf) # moyenne sur les alr
alrInv(moyenne) # transformation inverse de la moyenne des alr en moyenne de la composition
## Cl SO4
## 0.2262 0.3267 0.4471
## attr(,"class")
## [1] acomp
Les transformations ilr demande une partition binaire séquentielle (SBP). Elle n’est pas obligatoire: si vous n’en spécifiez pas, R en générera une automatiquement. Toutefois, bien que la SBP n’ait pas d’influence sur les résultats de vos statistiques une fois que vous inversez vos résultats pour retourner dans le simplex, une SBP basée sur une théorie peut vous aider à interpréter vos résultats, comme nous allons le voir dans le chapitre suivant.
Pour charger une SBP, vous faites de même que pour un tableau de données :
sbp1 <- read.csv("/home/serge-etienne/Documents/etudes/ulaval/RPubs/rSciencesSol/sbpHydrochem.csv",
header = TRUE, sep = ";", dec = ".")
Vous devez transformer votre SBP en matrice psi, qui tient compte des balances et de leur poids pour assurer l’orthogonalité. Puis, vous calculez les balances ilr. La dernière ligne, qui lance la fonction CoDaDendrogram, fait apparaître un dendrogramme compositionnel.
psi1 <- gsi.buildilrBase(t(sbp1))
balances1 <- ilr(Composition, V = psi1)
CoDaDendrogram(Composition, V = psi1)
Le dendrogramme compositionnel s’interprète comme suit :
Voici un dendrogramme plus complexe. Essayez de l’interpréter (les boxplots sont petits)!
sbp2 <- read.csv("/home/serge-etienne/Documents/etudes/ulaval/RPubs/rSciencesSol/sbp2Hydrochem.csv",
header = TRUE, sep = ";", dec = ".")
psi2 <- gsi.buildilrBase(t(sbp2))
CoDaDendrogram(acomp(Hydrochem[, 6:19]), V = psi2)
La transformation ilr:
La régression est utilisée pour analyser la relation entre une variable (ou un vecteur) et une ou plusieurs autres. La régression linéaire dans le simplex peut être utilisée pour modéliser un vecteur de données compositionnelles en fonction d’une variable extérieure. Pour plus d’information, consultez Tolosana-Delgado et van den Boogart (2011).
Nous allons mettre en relation la distribution granulométrique (sable, limon, argile) de sédiments arctiques en fonction de la profondeur. Ces données sont fournies avec le module compositions :
library(compositions)
data(ArcticLake) # charger les donnees
ArcticLake[1:5, ] # montrer les lignes 1 a 5
## sand silt clay depth
## 1 77.5 19.5 3.0 10.4
## 2 71.9 24.9 3.2 11.7
## 3 50.7 36.1 13.2 12.8
## 4 52.2 40.9 6.6 13.0
## 5 70.0 26.5 3.5 15.7
On peut explorer les données en créant un graphique du pourcentage en fonction de la profondeur :
plot(x = ArcticLake[, 1], y = ArcticLake[, 4], xlab = "%", ylab = "Profondeur (m)",
pch = 21, bg = 2, ylim = c(110, 0))
points(x = ArcticLake[, 2], y = ArcticLake[, 4], pch = 22, bg = 3)
points(x = ArcticLake[, 3], y = ArcticLake[, 4], pch = 23, bg = 4)
legend(x = 65, y = 80, legend = c("Sable", "Limon", "Argile"), pch = c(21, 22,
23), pt.bg = c(2, 3, 4))
La relation n’est évidemment pas linéaire. À première vue, ou pourrait appliquer sur chacune des classes granulométrique, une relation de croissance ou décroissance à saturation (tendance exponentielle avec convergence vers un maximum ou un minimum asymptotique). Ces relations s’appliqueraient toutefois sur des données non normales et ne tiendraient pas compte de l’interrelation entre les classes granulométriques. On passera donc par des données transformées :
comp <- acomp(ArcticLake[, 1:3]) # fermer le simplex colonnes 1 a 3
balances <- ilr(comp) # creer les balances ilr avec sbp par defaut
depth <- ArcticLake[, 4] # vecteur de profondeur: colonne 4
Nos balances étant créées, on peut les explorer. Dans notre cas, le log en base 10 de la profondeur est utilisé :
plot(y = log10(depth), x = balances[, 1], xlim = c(-4, 3), pch = 24, bg = 5)
points(y = log10(depth), x = balances[, 2], pch = 25, bg = 6)
Nous pouvons modéliser linéairement les balances en fonction du log10 de la profondeur, grâce à la fonction lm :
ilr.lm <- lm(balances ~ log10(depth))
Nous obtenons le diagnostique suivant :
summary(ilr.lm)
## Response Y1 :
##
## Call:
## lm(formula = Y1 ~ log10(depth))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.8874 -0.2865 0.0355 0.2604 0.9074
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.459 0.394 -8.77 1.4e-10 ***
## log10(depth) 2.681 0.243 11.02 3.0e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.444 on 37 degrees of freedom
## Multiple R-squared: 0.766, Adjusted R-squared: 0.76
## F-statistic: 121 on 1 and 37 DF, p-value: 3.04e-13
##
##
## Response Y2 :
##
## Call:
## lm(formula = Y2 ~ log10(depth))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6149 -0.3257 0.0537 0.2995 1.2279
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.921 0.643 -9.21 4.1e-11 ***
## log10(depth) 3.609 0.396 9.10 5.6e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.723 on 37 degrees of freedom
## Multiple R-squared: 0.691, Adjusted R-squared: 0.683
## F-statistic: 82.9 on 1 and 37 DF, p-value: 5.59e-11
Vous obtenez le R², ainsi que le résultat du test d’hypothèse sur les paramètres (Pr) pour
Au seuil 0,05, nous rejetons les hypothèses nulles pour l’origine et la pente pour les régressions sur les deux balances. Pour plus de détails sur les diagnostiques référez-vous au Chapitre 6 du livre Introductory statistics with R, à Tolosana-Delgado et van den Boogart (2011) ou à toute autre source expliquant la fonction lm. Nous obtenons les origines et les pentes pour chacune des deux balances:
coefficients(ilr.lm)
## [,1] [,2]
## (Intercept) -3.459 -5.921
## log10(depth) 2.681 3.609
Et la prédiction:
plot(x = balances[, 1], y = log10(depth), xlim = c(-4, 3), pch = 24, bg = 5)
points(x = balances[, 2], y = log10(depth), pch = 25, bg = 6)
lines(x = predict(ilr.lm)[, 1], y = log10(depth), lwd = 2, col = 5)
lines(x = predict(ilr.lm)[, 2], y = log10(depth), lwd = 2, col = 6)
Les coefficients sont sous forme de balance. Ils peuvent être transposés dans le simplex en inversant la transformation ilr.
coefficients <- ilrInv(coef(ilr.lm))
intercept <- acomp(coefficients[1, ]) #intercept (perturbation)
pente <- acomp(coefficients[2, ]) #log(depth) (puissance)
La commande acomp défini la ligne comme étant de classe compositionnelle. Lorsqu’un vecteur est de class acomp, les opérations dans le simplex (perturbation ⊕ et powering ⊗ sont effectuées en appelant respectivement + et *.
logDepthPred <- log10(1:500) # generer un vecteur de profondeur
compPred <- intercept + pente * logDepthPred # calculer avec pert. et pow. la rel. lineraire dans le simplex
Notez qu’il revient au même d’effectuer la prédiction en ilr et de convertir la prédiction de l’ilr en compositionnel par après.
ilr1Pred <- coef(ilr.lm)[1, 1] + coef(ilr.lm)[2, 1] * logDepthPred
ilr2Pred <- coef(ilr.lm)[1, 2] + coef(ilr.lm)[2, 2] * logDepthPred
ilrPred <- data.frame(ilr1Pred, ilr2Pred)
compPred <- ilrInv(ilrPred)
La régression peut être présentée sur graphique dans le diagramme ternaire ou en % versus profondeur.
# % versus profondeur
plot(y = depth, x = comp[, 1] * 100, ylab = "Profondeur (m)", xlab = "%", pch = 21,
bg = 2, ylim = c(110, 0), xlim = c(0, 100))
points(y = depth, x = comp[, 2] * 100, pch = 22, bg = 3)
points(y = depth, x = comp[, 3] * 100, pch = 23, bg = 4)
lines(y = 10^logDepthPred, x = compPred[, 1] * 100, col = 2)
lines(y = 10^logDepthPred, x = compPred[, 2] * 100, col = 3)
lines(y = 10^logDepthPred, x = compPred[, 3] * 100, col = 4)
grid()
legend(y = 50, x = 60, legend = c("Sable", "Limon", "Argile"), pch = 21:23,
pt.bg = c(2, 3, 4), cex = 0.8)
# diagramme ternaire
plot(comp, pch = 21, bg = 7) # points
lines(compPred, col = 3, lwd = 2) # regression lineaire
isoPortionLines()
Tolosana-Delgado et van den Boogart (2011). Linear models with compositions in R, dans Compositional data analysis : Theory and applications, wiley, 378 pages.