Ateliers R
Atelier 1
Opérations et objets de base
R peut être utilisé comme une simple calculatrice!
3 * 4## [1] 1213/3## [1] 4.3333333^3## [1] 27On peut également assigner des objets à des symboles avec l’opérateur <- et concaténer des nombres ensemble pour créer un vecteur, un objet à 1 dimension. Par la suite, on peut effectuer des opérations vectorisées.
v1 <- c(1, 2, 3, 4, 5, 6, 7, 8, 9)
v1## [1] 1 2 3 4 5 6 7 8 9v2 <- c(10, 10, 10, 10, 10, 10, 10, 10, 10)
v1 + v2## [1] 11 12 13 14 15 16 17 18 19v1 * v2 + 1## [1] 11 21 31 41 51 61 71 81 91On peut aussi coller ces deux vecteurs ensembles pour créer une matrice (matrix), un objet à 2 dimension.
m <- cbind(v1, v2) #! Colle deux vecteurs ensemble
m## v1 v2
## [1,] 1 10
## [2,] 2 10
## [3,] 3 10
## [4,] 4 10
## [5,] 5 10
## [6,] 6 10
## [7,] 7 10
## [8,] 8 10
## [9,] 9 10Une matrice est composée d’un même type de valeur, soit numérique, caractère ou logique. Plus général, un data.frame peut contenir plusieurs types de valeurs. C’est l’équivalent d’un tableau de données dans Excel.
d <- data.frame(v1, v2, z = c("a", "b", "c", "d", "e", "f", "g", "h", "i"))
d## v1 v2 z
## 1 1 10 a
## 2 2 10 b
## 3 3 10 c
## 4 4 10 d
## 5 5 10 e
## 6 6 10 f
## 7 7 10 g
## 8 8 10 h
## 9 9 10 iLes listes peuvent contenir tout type d’objets et sont des objets à une dimension.
l <- list(v1, m, d)
l## [[1]]
## [1] 1 2 3 4 5 6 7 8 9
##
## [[2]]
## v1 v2
## [1,] 1 10
## [2,] 2 10
## [3,] 3 10
## [4,] 4 10
## [5,] 5 10
## [6,] 6 10
## [7,] 7 10
## [8,] 8 10
## [9,] 9 10
##
## [[3]]
## v1 v2 z
## 1 1 10 a
## 2 2 10 b
## 3 3 10 c
## 4 4 10 d
## 5 5 10 e
## 6 6 10 f
## 7 7 10 g
## 8 8 10 h
## 9 9 10 iPrincipaux objets avec des données
- vecteur (1 dimension, 1 type de données)
- matrice (
matrix) (2 dimensions, 1 type de données) - data.frame (2 dimensions, tout type de données)
- liste (
list)(1 dimension, tout type d’objets) array(3 dimension, un type de données)
Type de données
Les principaux types de données ou de valeurs dans R sont:
class(1)
class("allo")
class(TRUE)## [1] "numeric"
## [1] "character"
## [1] "logical"Objets et assignation
Tout ce qui est contenu dans une session R est un objet qui porte un nom. Lorsque l’on modifie un objet, il faut se rappeler d’assigner le résultat afin de le conserver dans le même objet ou dans un nouvel objet. Par exemple, affichons le contenu de v1:
v1## [1] 1 2 3 4 5 6 7 8 9Maintenant, ajoutons une valeur à v1:
v1 + 100## [1] 101 102 103 104 105 106 107 108 109Le résultat s’affiche, mais l’objet n’a toujours pas changé:
v1## [1] 1 2 3 4 5 6 7 8 9L’objet v1 n’a toujours pas changé, car il faut assigner le résultat à un nouvel objet ou à l’objet v1 en “écrasant” l’ancien objet v1.
v1 <- v1 + 100
v1## [1] 101 102 103 104 105 106 107 108 109Tout ce qui est contenu dans une session R est un objet qui porte un nom.
Indexation
L’indexation est le processus par lequel on accède aux données à l’intérieur d’un objet. Pour ce faire, on utilise généralement les opérateurs [ et $ avec différentes conditions pour établir notre sélection.
Avec un vecteur
v1## [1] 101 102 103 104 105 106 107 108 109On accède au 3e élément du vecteur.
v1[3]## [1] 103On accède à tous les éléments supérieurs à 105.
v1[v1 > 105]## [1] 106 107 108 109On accède à tous les éléments supérieurs à 105 et inférieurs ou égal à 108.
v1[v1 > 105 & v1 <= 108]## [1] 106 107 108Avec un data.frame
Un data.frame est un objet en 2 dimensions qui contient des lignes et des colonnes d[i, j].
d## v1 v2 z
## 1 1 10 a
## 2 2 10 b
## 3 3 10 c
## 4 4 10 d
## 5 5 10 e
## 6 6 10 f
## 7 7 10 g
## 8 8 10 h
## 9 9 10 iOn sélectionne les lignes 4 à 6 du data.frame.
d[4:6, ]## v1 v2 z
## 4 4 10 d
## 5 5 10 e
## 6 6 10 fOn sélectionne les colonnes 1 et 3 du data.frame.
d[, c(1, 3)]## v1 z
## 1 1 a
## 2 2 b
## 3 3 c
## 4 4 d
## 5 5 e
## 6 6 f
## 7 7 g
## 8 8 h
## 9 9 iOn sélectionne la colonne nommée z.
d$z## [1] a b c d e f g h i
## Levels: a b c d e f g h iOn sélectionne toutes les lignes où v1 est supérieur à 7.
d[d$v1 > 7, ]## v1 v2 z
## 8 8 10 h
## 9 9 10 iAvec une liste
Avec les “[” simples, on accède aux éléments de la liste qui sont aussi retournés sous forme de liste.
l[1]## [[1]]
## [1] 1 2 3 4 5 6 7 8 9Avec les doubles “[[”, on accède directement à un seul élément contenu dans la liste.
l[[1]]## [1] 1 2 3 4 5 6 7 8 9On peut également accéder aux éléments de la liste par leur nom si cette liste est nommée.
l$v1## NULLPuisque la liste n’est pas nommée, la valeur NULL nous est retournée.
names(l) <- c("allo", "bonjour", "bye")
l$bye## v1 v2 z
## 1 1 10 a
## 2 2 10 b
## 3 3 10 c
## 4 4 10 d
## 5 5 10 e
## 6 6 10 f
## 7 7 10 g
## 8 8 10 h
## 9 9 10 iFonctions
On retrouve également dans R des objets qui sont des fonctions. Les fonctions peuvent être appliquées à des objets pour obtenir des résultats. On applique les fonctions aux objets en utilisant les parenthèses (), contrairement aux crochets [] qui servent à l’indexation. Par exemple, c est une fonction qui permet de coller des valeurs ensembles dans un vecteur.
temp <- c(37, 54, 79, 23, 3, 101, 4)La fonction mean permet de calculer des moyennes et la fonction sd des écarts-types.
mean(temp)
sd(temp)## [1] 43
## [1] 37.27823La fonction rep permet de répéter des valeurs. Ici, on répète le vecteur 1:2 5 fois.
rep(1:3, times = 5)## [1] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3Ceci permet d’introduire les arguments qui composent les fonctions. On pourrait aussi utiliser l’argument length.out ou each pour modfifier la façon dont les valeurs sont répétées.
rep(1:3, each = 5)
rep(1:3, length.out = 5)## [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3
## [1] 1 2 3 1 2En général, la majorité des arguments d’une fonction ont des valeurs par défaut. Par exemple, les valeurs par défaut de rep sont times = 1.
rep(1:3)## [1] 1 2 3Par défaut, le premier élément listé à l’intérieur d’une fonction est l’objet auquel on veut l’appliquer et correspond au premier argument de la fonction. Cet argument n’a généralement par besoin d’être nommé. Par contre, il est souvent préférable ou nécessaire de nommer les éléments suivants. Si ceux-ci ne sont pas nommés, l’ordre dans lequel ils sont donnés sera utilisé pour les associer aux arguments correspondant.
Opérateurs
Certaines fonctions sont représentées par des opérateurs. En voici quelques exemples.
3 == 4 # égal
3 + 4 # addition
TRUE & TRUE # 'et' logique
TRUE & FALSE # 'et' logique
3 == 4 | 3 == 5 # 'ou' logique
c(3, 4) %in% c(4, 5, 6, 7) # élément du vecteur 1 contenule vecteur 2
3:10 # crée un vecteur avec tous les entiers de 3 à 10## [1] FALSE
## [1] 7
## [1] TRUE
## [1] FALSE
## [1] FALSE
## [1] FALSE TRUE
## [1] 3 4 5 6 7 8 9 10Aide
Dès que l’on se met à jouer avec R, on est vite obligé d’aller inspecter son fichier d’aide. Ceci peut être tout aussi éclairant que frustrant! Mais c’est généralement éclairant et absolument essentiel. Pour accéder à l’aide sur la fonction rep, on n’a qu’à taper ?rep.
Section Description
Cette section contient une courte description de ce que fait la fonction.
Section Usage
Cette section montre comment la fonction est utilisée et liste également la plupart du temps les valeurs par défaut des arguments.
Section Arguments
Cette section explique l’utilisation de chacun des arguments de la fonction et donne également les valeurs par défaut.
Section Details
Cette section donne parfois des détails plus précis par rapport à l’utilisation de la fonction.
Section Value
Cette section décrit la valeur retournée par la fonction.
Section See Also
Cette section donne des suggestions concernant d’autres fonctions jouant un rôle similaire.
Section Examples
Cette section très utile donne des exemples d’utilisation de la fonction. Ces exemples peuvent être copiés et collés directement dans votre session R afin d’étudier l’utilisation de la fonction. Pour bien comprendre l’utilisation d’une fonction, cette section n’est pas à négliger!
Importation des données
Le répertoire courant (Working directory)
Lorsqu’on débute avec R, c’est souvent à cette première étape cruciale qu’on a des problèmes! Premièrement, R travaille dans un seul dossier à la fois, c’est-à-dire que la base de données doit être dans le bon dossier, appelé répertoire courant. Pour dire à R quel est ce dossier, on va dans RStudio sous l’onglet Session/Set Working Directory. On doit généralement faire cette opération à chaque début de session de travail dans R, à moins que RStudio s’ouvre par défaut au bon endroit. On peut également inspecter le répertoire courant à l’aide de la fonction getwd.
getwd()## [1] "C:/Users/rouf1703/Documents/UdeS/GitHub/ECL749/lab"Pour s’assurer que notre base de données se trouve bien à cet endroit, on peut utiliser la fonction list.files.
list.files()## [1] "lab.html" "lab.Rmd" "lab_cache" "lab_files" "rsconnect"Les données manquantes (NA)
Avant d’importer des données, il est nécessaire d’aborder la gestion des données manquantes par R. Les données manquantes sont toujours représentées par des NA, que ce soit dans les vecteurs ou dans les tableaux de données. Lorsqu’on veut importer des données avec des cellules vides, il est souvent nécessaire de remplacer les cellules vides par des NA. Parfois, il est possible de dire à une éventuelle fonction que les cellules vides correspondent à des NA.
Les fonctions d’importation
Il existe plusieurs méthodes pour importer des données dans R. La façon classique est par l’intermédiaire de la fonction read.table. Cette fonction prend des fichiers textes (.txt) qui doivent répondre à plusieurs critéres:
Il ne doit y avoir aucun espace dans les noms des différentes variables et dans les éléments de la base de données. Par exemple, « bon à modéré» va créer des problèmes. On devrait plutôt écrire “bon.à.modéré” ou “bon_à_modéré”. Encore mieux serait d’éviter les accents ou d’utiliser un code.
AUCUNE cellule ne doit être vide. Si une cellule est vide ou si on a écrit un commentaire au lieu d’un nombre, on écrit
NAdans la cellule. Autrement, R produira un message d’erreur puisqu’il considérera qu’une ligne du tableau est incomplète.Il faut généralement utiliser le point plutôt que la virgule pour les nombres décimaux (à moins de le spécifier dans les arguments)
Les messages d’erreurs sont pratiquement obligatoires lors de l’importation d’une nouvelle base de données. Il s’agit juste d’être en mesure de bien les interpréter et d’observer ces quelques règles de base.
Plusieurs autres fonctions sont disponibles selon le format des données (read.delim, read.csv, etc.). Le format .csv est souvent le plus pratique, car il réduit les risques de problèmes lors de l’importation en raison de la délimitation des cellules. Il est également possible d’utiliser le package readxl pour lire directement des fichiers .xlsx. Bien qu’étant la façon classique d’importer des données, la fonction read.table est assez capricieuse… On peut également exporter des données avec la fonction write.table.
Packages
Les packages sont des banques de fonctions destinées à des tâches spécifiques. Plusieurs packages de bases sont installés avec R, mais il existe plusieurs milliers de packages de toutes sortes. Pour installer un package donné sur son ordinateur (par exemple le package spatstat, pour les statistiques spatiales), on utilise la fonction install.packages de cette façon :
install.packages("spatstat")Dépendamment de si on travaille avec R ou RStudio, on demande parfois de sélectionner un miroir CRAN à partir duquel le téléchargement sera fait. Une fois le package installé, il faut activer le package dans notre session de travail R. Pour les packages spéciaux qui ne font pas partie des packages de base de R, on doit les activer à chaque session de travail avec la fonction library
library("spatstat")Il est fréquent d’utiliser une fonction dans R sans avoir préalablement charger le package la contenant. Dans ce cas, un message d’erreur nous dit que la fonction utilisée n’existe pas. Ceci peut également résulter d’une mauvaise orthographe de la fonction.
Exercice
Ce premier exercice consiste à importer une base de données en format .txt. Cette base de données contient des mesures morphométriques d’une espèce de bruants (un oiseau). Vous pouvez la télécharger à l’adresse suivante: https://github.com/frousseu/ECL749/tree/master/data.
Par la suite, vous devez calculer la longueur moyenne de la corde de l’aile (wingcrd) et de la masse (wt) et effectuez un graphique mettant en relation les deux mesures. La fonction plot devrait vous être utile.
Vous devez également utiliser la fonction table pour compter le nb d’observations de mâles et de femelles et apprendre la fonction boxplot pour comparer la variabilité de la masse chez les deux sexes.
Finalement, effectuez et interprétez un test de t pour comparer la masse chez les deux sexes. Ici, la fonction subset pourrait être utile, mais il est possible de s’en sortir en utilisant seulement l’indexage.
Il y a peut-être quelques petits problèmes avec la base de données…
Solution
On télécharge le fichier sparrows.txt avec le bouton de droit en appuyant sur le bouton raw. On met ce fichier à un endroit correspondant au répertoire courant ou on change le répertoire courant pour accomoder la localisation du fichier. Vous pouvez aussi utiliser le menu de RStudio Session / Set Working Directory / Choose Directory. Remarquez que la ligne setwd("chemin/vers/le/fichier.txt") va s’écrire dans la console. Le mieux serait d’inclure cette ligne dans votre script. Remarquez qu’il serait aussi possible d’écrire le chemin vers le fichier au long au lieu de seulement inscrire le nom du fichier (entre guillemets) pour la fonction read.table.
setwd("C:/Users/rouf1703/Desktop")
d <- read.table("sparrows.txt", header = TRUE)Le message d’erreur retourné nous dit qu’il manque un élément à la ligne 83 du tableau de données. Il est important de bien lire les messages d’erreur, car la plupart du temps, ils nous donnent une indication du problème. Le plus simple est probablement d’aller dans le tableau de données et de remplacer manuellement la cellule manquante par un NA.
setwd("C:/Users/rouf1703/Desktop")
d <- read.table("sparrows.txt", header = TRUE)Maintenant qu’on a les données, on peut calculer les moyennes des deux variables wingcrd et wt.
mean(d$wingcrd)## [1] NAmean(d$wt)## [1] 19.67126On on obtient un NA pour la première moyenne et un message d’erreur pour la seconde. Inspectons ce que contient le vecteur d$wingcrd pour comprendre pourquoi un NA est retourné. On voit que le vecteur contient une valeur de NA et par défaut, la fonction mean retourne un NA (?mean). Pour corriger cette situation, on doit utiliser l’argument na.rm = TRUE de la fonction mean.
mean(d$wingcrd, na.rm = TRUE)## [1] 57.63953mean(d$wt, na.rm = TRUE)## [1] 19.67126On obtient une valeur pour la longueur de l’aile, mais on obtient toujours un message d’erreur pour la masse. Inspectons le contenu de d$wt pour comprendre un peu mieux.
d$wt## [1] 20.4 21.5 20.2 21.3 16.4 19.5 16.1 22.0 19.0 19.9 20.5 19.2 19.8 20.5
## [15] 20.8 21.2 19.5 18.8 20.3 18.8 25.0 19.5 21.3 19.8 20.2 19.5 20.5 21.0
## [29] 20.4 24.2 17.2 21.2 20.8 17.0 16.5 22.5 16.4 19.1 18.1 18.6 21.9 17.2
## [43] 17.9 18.1 20.6 21.0 19.8 17.1 19.2 23.2 18.4 16.5 25.5 20.2 20.0 18.0
## [57] 18.8 21.1 25.4 20.6 17.6 16.0 16.2 21.8 20.0 20.0 20.9 19.4 18.5 19.8
## [71] 20.0 18.8 18.0 20.6 20.0 18.5 21.6 16.9 16.7 17.0 21.0 20.0 18.1 21.2
## [85] 17.6 21.2 19.0On semble avoir des valeurs, mais on voit également que le résultat indique qu’il y a 48 Levels dans notre vecteur. Ceci indique que R a transformé la colonne en facteur. Par défaut, la fonction read.table converti en facteur toute colonne qui ne semble pas être composée de données numériques ou logiques. Ceci vaut certainement pour les données de type caractère. Puisqu’une virgule se trouve quelque part dans notre série de valeurs, R a considéré que ces données étaient des facteurs et a converti le tout en conséquence. Pour comprendre un peu mieux l’utilisation des facteurs avec R, consultez la section d’aide sur la fonction ?factor. En bref, les factor sont des valeurs comportant une étiquette (la plupart du temps un caractère) et étant représentées par un entier allant de 1 au nombre d’étiquettes différentes. Pour éviter ce problème, le plus simple serait de retourner dans le tableau de données initiales et de changer la virgule pour un point pour que la colonne soit interprétée comme une colonne de valeurs numériques.
Maintenant que nous avons toutes nos valeurs, nous pouvons calculer les deux moyennes et effectuer un graphique représentant le lien entre ces deux mesures avec la fonction plot.
mean(d$wingcrd, na.rm = TRUE)## [1] 57.63953mean(d$wt, na.rm = TRUE)## [1] 19.67126plot(d$wingcrd, d$wt)
Sans surprise, les deux variables semblent positivement corrélées. Par la suite, on utilise la fonction table pour compter le nombre d’observations des deux sexes.
table(d$Sex)##
## f m
## 42 45La fonction boxplot nous permet de comparer la masse chez les deux sexes.
boxplot(wt ~ Sex, data = d)
Remarquez la forme des arguments de la fonction boxplot. Ceux-ci sont donnés sous forme de formule (?formula) et non sous forme de caractères entre guillemets. Avec l’argument data, la fonction boxplot va chercher dans le data.frame les noms de colonnes correspondants. Cette façon de spécifier les variables sous forme de formule sera beaucoup utilisée lorsque nous ferons des modèles.
Par la suite, nous devons trouver le moyen de séparer les mâles des femelles afin d’effectuer notre test de t. Pour ce faire, il est possible d’utiliser la fonction subset ou simplement d’utiliser l’indexage. Les 4 lignes suivantes produisent les mêmes objets.
m <- subset(d, subset = Sex == "m")
f <- subset(d, subset = Sex == "f")
m <- d[d$Sex == "m", ]
f <- d[d$Sex == "f", ]Un peu comme dans la fonction boxplot, il n’est pas nécessaire d’utiliser les guillemets ou le $ puisque la fonction subset cherchera dans le data.frame soumis pour trouver la colonne demandée (ici, Sex). Avec l’indexage régulier, il faut toutefois identifé où se trouve la colonne avec d$Sex. Par la suite, on effectue notre test de t.
t.test(m$wt, f$wt)##
## Welch Two Sample t-test
##
## data: m$wt and f$wt
## t = 3.808, df = 69.974, p-value = 0.0002976
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.7452093 2.3843145
## sample estimates:
## mean of x mean of y
## 20.42667 18.86190La valeur de p est inférieure à 0.05 et l’intervalle de confiance a des bornes inférieure et supérieures positives, on conclut donc que les mâles ont une masse supérieure aux femelles et que cette différence est de l’ordre de 1 ou 2 grammes.
Atelier 2
Retour
Script et console
Un script R est un fichier .R où on écrit le code qu’on veut conserver et à partir duquel on peut envoyer les différentes lignes de codes vers la console R. La console est l’endroit où sont évalué les différentes commandes et où les objets de la session R sont contenus. On peut écrire directement dans la console pour envoyer des commandes, mais il est préférable d’écrire ces commandes dans un script lorsqu’on veut conserver les commandes importantes.
Insertion de commentaires
Pour insérer des commentaires dans un script, on utilise le dièse # (ou, pour utiliser un langage plus moderne, le hashtag…). Ces lignes ne seront pas interprétées lorsque’elle seront envoyées à la console.
# Voici un commentaire qui ne sera pas interprété comme une commandeEnvironnement de travail
Il est possible de sauvegarder à la fois un script (un fichier .R) qui contient du code ou un environnement de travail (un fichier .RData, le Workspace) qui contient les objets créés dans l’environnement. En général, il est beaucoup plus utile de conserver un script avec le code important puisque celui-ci peut être relancé à nouveau et c’est celui qui contient le travail important. Sauvegarder l’environnement de travail n’est généralement utile que si on crée des objets longs à produire ou si on veut partager des objets avec d’autres personnes par exemple. Pour sauvergarder l’environnement de travail, on peut aller sous le menu de RStudio et faire Session / Save Workspace As. Il est également possible de conserver l’historique des commandes, mais ceci n’est pas très utile puisque cela conserve l’ensemble des commandes qui ont été envoyées à la console, y compris les commandes contenant des erreurs.
Importation de données
Comme vous avez pu le constater, l’importation de données avec read.table peut être ardue lorsqu’on débute. Pour éviter beaucoup de problème, une inspection préalable de la base de données est souvent nécessaire. Le message d’erreur le plus commun est probablement celui-ci.
d <- read.table("sparrows2.txt", header = TRUE)## Warning in file(file, "rt"): impossible d'ouvrir le fichier
## 'sparrows2.txt' : No such file or directory## Error in file(file, "rt"): impossible d'ouvrir la connexionCe message indique que soit on a fait une erreur dans le nom du fichier, soit le fichier ne se trouve pas dans le répertoire courant. Dans cet exemple, il s’agit d’une erreur dans le nom du fichier. Remarquez également l’utilisation de l’argument header pour spécifier que la première ligne du tableau de données contient les noms de colonnes. Un autre argument très utilisé est stringsAsFactors = FALSE qui permet de spécfier que les données de type caractère ne doivent pas être automatiquement converties en factor.
Nom du fichier / chemin
Avec la plupart des fonctions d’importation, on peut soit donner le nom du fichier (avec l’extension) si celui-ci se trouve dans le répertoire courant, soit donner le chemin complet vers où le fichier se trouve.
d <- read.table("C:/Users/rouf1703/Desktop/sparrows.txt", header = TRUE)De cette façon, on s’affranchit du répertoire courant, car on indique directement à R où aller chercher le fichier. Ceci peut être très utile lorsqu’on ne veut pas changer le répertoire courant à chaque session ou si on a des fichiers classés à différents endroits sur l’ordinateur.
Autres fonctions d’importation
L’importation de fichiers .txt est la base, mais il est souvent plus facile d’utiliser le format .csv et la fonction read.csv. Ce format cause généralement moins de problèmes avec les cellules vides et la séparation entre les différentes colonnes est souvent moins ambigue qu’avec les espaces ou les tabulations.
Il est également possible de lire directement dans un fichier Excel avec les packages readxl et xlsx, mais je vais vous laisser le soin d’étudier ces possibilités.
Inspection des données
Une fois qu’une base de données est importée, plusieurs fonctions peuvent être utilisées pour inspecter un data.frame. Les plus communes sont les fonctions head, str, summary. Voici des exemples.
head(d)## Sex wingcrd wt
## 1 m 56.5 20.4
## 2 f 60.5 21.5
## 3 m 60.0 20.2
## 4 m 61.0 21.3
## 5 f 55.0 16.4
## 6 f 58.0 19.5str(d)## 'data.frame': 87 obs. of 3 variables:
## $ Sex : Factor w/ 2 levels "f","m": 2 1 2 2 1 1 1 2 1 2 ...
## $ wingcrd: num 56.5 60.5 60 61 55 58 51.5 59 54 59 ...
## $ wt : num 20.4 21.5 20.2 21.3 16.4 19.5 16.1 22 19 19.9 ...summary(d)## Sex wingcrd wt
## f:42 Min. :51.50 Min. :16.00
## m:45 1st Qu.:56.50 1st Qu.:18.25
## Median :57.75 Median :19.80
## Mean :57.64 Mean :19.67
## 3rd Qu.:59.00 3rd Qu.:20.85
## Max. :63.50 Max. :25.50
## NA's :1La fonction head affiche les 6 premières lignes d’un data.frame. La fonction str liste les colonnes et informe sur le type de données dans chacune de celles-ci. La fonction summary calcule un résumé de l’information dans chaque colonne qui dépend du type de données.
Il est également possible d’inspecter les objets en utilisant l’onglet Environment dans RStudio et en cliquant sur les différents objets.
Les guillemets
Les guillemets sont régulièrement utilisés dans R pour indiquer des valeurs correspondant à des caractères. Par exemple, le nom d’un fichier doit être donné entre guillemets à la fonction read.table. Pour certaines fonctions (e.g. library, subset) et pour les fonctions avec des formules et un argument data, les guillemets ne sont pas nécessaires. Il peut être frustrant au début de distinguer entre ces deux possibilités, mais on finit par s’y habituer. Encore une fois, se référer au fichier d’aide!
= vs. ==
Le signe d’égalité simple = est la plupart du temps utilisé pour spécifier un argument à une fonction. Il peut parfois être utilisé pour faire des assignations au lieu de la flèche d’assignation <-, mais on recommande d’éviter cette pratique. Le double == est utilisé comme une fonction pour tester l’égalité entre des valeurs ou des objets.
Les facteurs
Étudions un peu plus les facteurs. Comme dit plus haut, un factor est un vecteur avec différents niveaux ayant une étiquette et représenter internalement par des entiers. Par défaut, les étiquettes sont classées par ordre alphabétique pour établir la représentation en entier. En général, avant de faire une analyse statistique comportant des facteurs, il faut s’assurer que nos variables catégoriques sont converties en tant que tel avec les fonctions factor ou as.factor. Ceci est particulièrement important pour les facteurs représentés par des valeurs numériques. Autrement, les modèles risquent de considérer ces variables comme numériques.
v <- c("c", "b", "b", "a", "a", "a")
v <- factor(v)
v## [1] c b b a a a
## Levels: a b cas.numeric(v)## [1] 3 2 2 1 1 1as.character(v)## [1] "c" "b" "b" "a" "a" "a"levels(v)## [1] "a" "b" "c"mean(v)## Warning in mean.default(v): argument is not numeric or logical: returning
## NA## [1] NAv + 1## Warning in Ops.factor(v, 1): '+' not meaningful for factors## [1] NA NA NA NA NA NAConsulter l’aide!
Avant d’utiliser une fonction, il faut se rappeler que c’est toujours un bon réflexe de consulter l’aide et particulièrement les exemples qui montrent comment utliser la fonction et ses arguments.
Brève introduction aux graphiques
La principale fonction pour créer des graphiques est la fonction plot. Certains arguments de cette fonction sont décrits dans le fichier d’aide de la fonction (?plot), mais la plupart des arguments sont décrits dans la page de la fonction ?par qui sert à établir les paramètres graphiques utilisés. Voici quelques paramètres (ou arguments) graphiques fréquemment utilisés. Notez la fonction par et l’argument mfrow qui sert à créer une fenêtre graphique 2 x 2. L’argument mar de la fonction par permet de modifier la largeur des marges du graphique (bas, gauche, haut et droite).
par(mfrow = c(2, 2), mar = c(5, 5, 2, 2))
plot(1:20, 1:20, pch = 1:20, col = "black", cex = 2) # graphique 1
plot(1:20, 1:20, pch = 16, col = 1:20, cex = 3) # graphique 2
plot(1:20, 1:20, pch = 1, col = "red", cex = 1:20/4) # graphique 3
plot(1:20, 1:20, xlim = c(5, 15), ylim = c(5, 15)) # graphique 4
points(1:20, 1:20 + 1, pch = 2, col = "magenta")
lines(1:20, 1:20 + 2, lwd = 2, lty = 1, col = "blue")
text(12, 6, labels = "Voici du texte", font = 2)
La fonction plot permet d’établir la fenêtre graphique. Par la suite, on peut ajouter des points, des lignes ou du texte à un graphique à l’aide des fonctions points, lines et text. Notez que ces fonctions ne permettent pas de créer ou d’ouvrir un graphique. Elles ne permettent que d’ajouter des éléments sur un graphique déjà exsitant. La fonction legend permet d’identifier des éléments sur un graphique.
par(mfrow = c(1, 1))
plot(1:10, 1:10, pch = 16, col = "red", xlim = c(0, 10), ylim = c(0, 10), xlab = "Valeurs X",
ylab = "Valeurs Y", main = "Voici un titre")
points(1:10, 1:10 - 1, pch = 16, col = "blue")
legend("topleft", pch = 16, col = c("red", "blue"), legend = c("Série 1", "Série 2"),
title = "Titre de légende")
ggplot2
Il existe trois principaux systèmes dans R pour effectuer des graphiques. Les graphiques de bases sont représentés par les fonctions plot, par et les fonctions associées. Un autre système, passablement plus compliqué est le package lattice qui est souvent utilisé par d’autres fonctions. Un autre système, très populaire, est ggplot2 qui est de plus en plus utilisé. Un des principaux avantages de ggplot2 est que ce système permet de faire rapidement des graphiques plus attrayants comparativement aux graphiques de bases. La documentation concernant ggplot2 est également très détaillée et bien expliquée pour quiconque voudrait employer ce système. Toutefois, dans ce cours, nous n’utiliserons que les graphiques de bases puisque ceux-ci sont plus souvent utilisés par d’autres packages utiles en écologie.
Exercice
Pour cet exercice, nous utiliserons une base de données ( pines.csv ) portant sur la croissance de jeunes Pins blancs (Pinus strobus) en lien avec avec l’utilisation d’un fertilisant, la déprédation par le Cerf de Virginie, la compétition avec la strate arbustive et l’espacement entre les individus. Cette base de données comporte plusieurs colonnes:
- height: hauteur de l’arbre
- width: largeur de l’arbre (maximale avec les branches)
- height_ini: hauteur initiale de l’arbre
- row: position de l’arbre selon les rangées de la plantation
- col: position de l’arbre selon les colonnes de la plantation
- deer: présence de broutage par le cerf (0 = non, 1 = oui)
- fert: utilisation d’un fertilisant (0 = non, 1 = oui)
- spacing: espacement avec les voisins (10 ou 15m)
- cover: recouvrement arbustif (0 = aucun, 1 = un peu, 2 = modéré, 3 = beaucoup)
Vous devrez réaliser plusieurs tâches et/ou répondre aux questions suivantes:
Importez la base de données.
Explorez et illustrez les liens entre la hauteur et la largeur et les différentes variables explicatives (utilisez
boxplot,mfrow).Effectuez une ANOVA afin d’étudier l’effet du broutement, de la couverture arbustive, du fertilisant et de l’espacement sur la hauteur et la largeur des individus. Pour ce faire, vous pouvez utiliser les fonctions
lmouaoven combinaison avec la fonctionsummaryInterprétez également les résultats.Réalisez une carte schématique permettant d’illustrer la plantation et les traitements (espacement et fertilisant). Pour ce faire, vous devrez utiliser les fonctions
plotetlegend. Voux aurez également besoin des différents paramètres graphiques (arguments) décrit plus haut (col,cex,pch). La tâche pourrait aussi être facilitée en créant des groupes selon les différentes possibilités.
Solution
1.
En premier, on importe la base de données. On peut visualiser son contenu avant l’importation ou après avec différentes fonctions comme head ou summary.
d <- read.csv("U:/2018/2018/Ateliers/pines.csv", header = TRUE)
head(d)d <- read.csv("https://raw.githubusercontent.com/frousseu/ECL749/master/data/pines.csv",
header = TRUE)
head(d)## height width height_ini row col deer fert spacing cover
## 1 NA NA NA 1 1 NA 0 15 0
## 2 362 162 14 1 2 1 0 15 2
## 3 442 250 17 1 3 0 0 15 1
## 4 NA NA NA 1 4 NA 0 15 0
## 5 369 176 24 1 5 0 0 15 2
## 6 365 215 22 1 6 0 0 15 12.
Par la suite, on peut explorer les liens entre les variables réponses (hauteur et largeur) et les variables explicatives. Pour ce faire, on peut utiliser le boxplot et créer deux graphiques à quatre panneaux contenant toutes les variables.
par(mfrow = c(2, 2), mar = c(4, 3, 1, 1))
boxplot(height ~ deer, data = d, xlab = "deer")
boxplot(height ~ fert, data = d, xlab = "fert")
boxplot(height ~ spacing, data = d, xlab = "spacing")
boxplot(height ~ cover, data = d, xlab = "cover")
On peut faire de même pour la largeur des arbres.
par(mfrow = c(2, 2), mar = c(4, 3, 1, 1))
boxplot(width ~ deer, data = d, xlab = "deer")
boxplot(width ~ fert, data = d, xlab = "fert")
boxplot(width ~ spacing, data = d, xlab = "spacing")
boxplot(width ~ cover, data = d, xlab = "cover")
3.
Pour effectuer une ANOVA, on peut utiliser la commande très générale lm. Avant de faire ceci, il faudra toutefois convertir nos variables numériques en facteurs, autrement, lm interprétera les différentes variables comme étant des variables numériques. Par exemple, dans la sortie suivante, on peut voir qu’il n’y a qu’une seule ligne par variable et on ne voit pas les niveaux des facteurs (1, 2 ou 3) à côté des noms de variables.
m <- lm(height ~ deer + fert + spacing + cover, data = d)
summary(m)##
## Call:
## lm(formula = height ~ deer + fert + spacing + cover, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -289.527 -40.886 8.255 50.563 188.294
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 453.269 14.356 31.574 < 2e-16 ***
## deer 18.963 5.774 3.284 0.00106 **
## fert -38.191 5.307 -7.196 1.35e-12 ***
## spacing -6.992 1.085 -6.444 1.94e-10 ***
## cover -6.821 2.436 -2.800 0.00523 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 77.84 on 858 degrees of freedom
## (137 observations deleted due to missingness)
## Multiple R-squared: 0.1137, Adjusted R-squared: 0.1096
## F-statistic: 27.52 on 4 and 858 DF, p-value: < 2.2e-16Avant de faire notre modèle, on doit donc convertir nos variables en facteur en utilisant la fonction as.factor. Ce problème ne se poserait pas si nous avions des niveaux représentés par des caractères, puisque la fonction read.csv aurait automatiquement converti ces valeurs en factor. Cela n’indique pas qu’il faut nécessairement utiliser des caractères pour désigner des facteurs, mais seulement qu’il faut porter attention aux niveaux identifiés par des entiers, ce qui est relativement commun dans les bases de données.
d$deer <- as.factor(d$deer)
d$fert <- as.factor(d$fert)
d$spacing <- as.factor(d$spacing)
d$cover <- as.factor(d$cover)
m <- lm(height ~ deer + fert + spacing + cover, data = d)
summary(m)##
## Call:
## lm(formula = height ~ deer + fert + spacing + cover, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -301.289 -40.326 8.071 53.277 187.410
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 366.819 6.583 55.718 < 2e-16 ***
## deer1 16.747 5.668 2.954 0.00322 **
## fert1 -34.710 5.219 -6.651 5.18e-11 ***
## spacing15 -35.660 5.307 -6.719 3.34e-11 ***
## cover1 21.470 7.276 2.951 0.00326 **
## cover2 21.452 7.397 2.900 0.00383 **
## cover3 -23.520 7.562 -3.110 0.00193 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 76.13 on 856 degrees of freedom
## (137 observations deleted due to missingness)
## Multiple R-squared: 0.1543, Adjusted R-squared: 0.1484
## F-statistic: 26.03 on 6 and 856 DF, p-value: < 2.2e-16On voit maintenant que les noms des coefficients ont changé puisqu’ils se rapportent maintenant aux niveaux des différents facteurs. Maintenant, si on interprète l’effet de nos variables, on voit que le broutement par le cerf à un effet positif sur la taille, mais l’utilisation du fertilisant et l’espacement entre les arbres ont tous deux un effet négatif sur la hauteur des arbres. Quant au recouvrement arbustif, des valeurs faibles ou modérées ont un effet positif comparativement à l’absence de recouvrement arbustif, mais un recouvrement élevé à un effet négatif comparativement à toutes les autres classe de recouvrement. Globablement, on peut aussi dire que nos variables expliquent une portion significative de la variance en hauteur, comme en témoigne le test du F. La proportion de variance expliquée (R2) est d’environ 15%. Comment interpréteriez vous ces résultats?
Je vous laisse explorer l’effet des variables sur la largeur des arbres.
4.
Pour construire une carte schématique de la plantation en fonction des traitements (fertilisation et espacement), nous pouvons utiliser les numéros de lignes et de colonnes de chacun des arbres. Il y a plusieurs options pour associer des symboles aux traitements. À ce stade, le plus simple est probablement de créer différents sous-ensembles et d’ajouter les points sur un graphique déjà existant. En illustrant la localisation de tous les arbres en premier, on s’assure que la fenêtre graphique couvrira l’ensemble de la plantation. Pour éviter que des points se tracent avec la fonction plot initiale, on peut soit utiliser la couleur blanche avec l’argument col = "white" ou encore utiliser l’argument type = "n". Ceci produit le graphique, mais sans tracer les points associés. On peut également utiliser la fonction legend pour identifier les différentes possibilités. Pour éviter de taper des commandes trop longue, j’utilise un vecteur leg où j’écris le texte qui composera la légende.
d1 <- d[d$fert == 0 & d$spacing == 10, ]
d2 <- d[d$fert == 0 & d$spacing == 15, ]
d3 <- d[d$fert == 1 & d$spacing == 10, ]
d4 <- d[d$fert == 1 & d$spacing == 15, ]
plot(d$col, d$row, type = "n")
points(d1$col, d1$row, col = "red", pch = 16)
points(d2$col, d2$row, col = "blue", pch = 16)
points(d3$col, d3$row, col = "red", pch = 1)
points(d4$col, d4$row, col = "blue", pch = 1)
leg <- c("fert = 0, spacing = 10", "fert = 0, spacing = 15", "fert = 1, spacing = 10",
"fert = 1, spacing = 15")
legend("bottomright", pch = c(16, 16, 1, 1), col = c("red", "blue", "red", "blue"),
legend = leg, bty = "n", title = "Traitements")
Atelier 3
Retour
L’argument sep
L’argument sep permet de spécifier aux fonctionx read.table ou read.csv quel est le séparateur utilisé dans un fichier .csv. Par défaut, sep = "," dans la fonction read.csv. Si vous spécifiez sep = ";", alors qu’en réalité vos données sont séparées par une virgule, vos données seront lues et considérées comme faisant parti d’une seule colonne. Consultez la page d’aide sur la fonction ?read.table pour plus de détails.
Vecteurs vs data.frame
Plusieurs ont été tenté de renommer les colonnes et de se créer des vecteurs indépendants pour ensuite les soumettre aux fonctions boxplot ou lm. Ceci est plutôt à éviter. Premièrement, cela implique plus de travail et plus de lignes de codes. Ensuite, en créant des vecteurs indépendants, on court le risque de perdre l’ordre selon lequel les valeurs sont emmagasinées ou d’avoir des vecteurs de longueurs différentes si on décide de supprimer les valeurs manquantes dans certains vecteurs, mais qu’on oublie de le faire dans d’autres. La meilleure façon de procéder est de faire usage de l’argument data lorsque celui-ci est disponible et donc d’utiliser notre data.frame. De plus, lorsque l’argument data est utilié, il devient inutile d’utiliser le $ pour indiquer nos variables. Par exemple, les deux lignes suivantes sont équivalentes.
boxplot(d$height ~ d$deer, data = d)
boxplot(height ~ deer, data = d)Les formules
Les formules sont beaucoup utilisées dans R, particulièrement lorsqu’on fait des modèles. Une formule s’écrit généralement de cette façon: réponse ~ var1 + var2. Le symble ~ veut souvent dire “expliqué par” ou “groupé par”. Par exemple, le premier argument à la fonction boxplot est une formule. Les arguments donnés sous cette forme sont généralement accompagnés d’un argument data.
boxplot avec num ~ num
Les boxplots sont faits pour étudier les valeurs de variables quantitatives en fonction de variables catégoriques. Si vous faites un boxplot pour deux valeurs numériques, vous obtiendrez une boîte pour chacune des valeurs uniques de la variable “explicative”, ce qui n’est pas très utile.
boxplot(height ~ width, data = d)
Transformation en facteurs
La plupart des fonctions servant à faire des modèles vont exiger que les variables catégoriques soient préalablement transformées en facteur. Dans certains cas, les variables caractères seront converties en facteurs par la fonction utilisée. Cependant, si vos facteurs sont représentés par des valeurs numériques, il est impératif de les transformer en facteurs, autrement, ils seront considérés comme des variables numériques et vos analyses seront erronées.
Mieux comprendre son modèle avec predict
La fonction predict est utilisée pour générer des prédictions à partir d’un modèle déjà existant. Par exemple, voici comment utiliser cette fonction avec un modèle fictif. Premièrement, créeons un modèle fictif avec une variable explicative numérique (x) et une variable explicative catégorique (sexe).
n <- 100
x <- runif(n, min = 0, max = 100)
val <- rep(c(0, 100), each = n/2)
sexe <- rep(c("m", "f"), each = n/2)
y <- 100 + 2 * x + val + rnorm(n, mean = 0, sd = 20)
d <- data.frame(y, x, sexe)
plot(d$x, d$y, col = ifelse(d$sexe == "m", "red", "blue"))
Voici l’estimation des paramètres du modèle avec la fonction lm.
m <- lm(y ~ x + sexe, data = d)
summary(m)##
## Call:
## lm(formula = y ~ x + sexe, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -39.947 -12.817 -0.478 10.461 41.145
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 199.13742 3.91798 50.83 <2e-16 ***
## x 2.00094 0.06633 30.17 <2e-16 ***
## sexem -101.37206 3.69129 -27.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.14 on 97 degrees of freedom
## Multiple R-squared: 0.9355, Adjusted R-squared: 0.9342
## F-statistic: 703.4 on 2 and 97 DF, p-value: < 2.2e-16Avec predict, on soumet des données à notre modèle pour générer des prédictions. Ici, on demande de faire varier x de 0 à 100 et on fixe le facteur sexe à “m”. On met ces valeurs dans un data.frame qui sera utilisé par predict pour générer les prédictions. Ensuite, on traces ces valeurs prédites à l’aide fonction lines.
dat <- data.frame(x = 0:100, sexe = "m")
p <- predict(m, newdata = dat)
plot(d$x, d$y, col = ifelse(d$sexe == "m", "red", "blue"))
lines(0:100, p, col = "red")
Exercices
Les coûts de la reproduction
Le premier exercice consiste à étudier l’effet de la reproduction sur la longévité chez les mâles de drosophiles (une espèce de mouche). Voici la base de données flies.csv. Cinq traitements ont été utilisés. Les mâles ont été mis en présence de 0, 1 ou 8 femelles et celles-ci étaient soient vierges, soient gestantes (donc non-réceptives à l’accouplement). D’autres variables ont été mesurées, notamment la longueur du thorax des mâles et la proportion du temps passé à dormir. Dans cet exercice vous devez:
Explorez les données
Déterminer si la longévité est influencée par le traitement
Vérifiez les conditions d’applications avec la fonction
plotappliquée sur votre modèle.
Solution
On importe d’abord le fichier et on effectue quelques graphiques pour explorer les données.
d <- read.csv("https://raw.githubusercontent.com/frousseu/ECL749/master/data/flies.csv")
par(mfrow = c(1, 3), mar = c(4, 4, 1, 1))
plot(d$Thorax, d$Longevity)
plot(d$Sleep, d$Longevity)
boxplot(Longevity ~ Treatment, data = d)
La longévité semble beaucoup plus faible pour les mâles exposés à plusieurs femelles réceptives. Aussi, la taille du thorax semble avoir un effet positif très fort sur la longévité. Il serait donc important d’inclure cette variable comme contrôle (même si elle n’est pas spécifiquement d’intérêt) afin de réduire le bruit dans nos observations. La proportion du temps passé à dormir ne semble pas avoir d’effet sur la longévité, mais vérifions tout de même ceci en incluant ces variables dans un modèle.
m <- lm(Longevity ~ Treatment + Thorax + Sleep, data = d)
summary(m)##
## Call:
## lm(formula = Longevity ~ Treatment + Thorax + Sleep, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.331 -7.520 -0.632 6.462 29.760
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -46.27094 10.48525 -4.413 2.27e-05 ***
## Treatment1 virgin -9.55148 2.96946 -3.217 0.00167 **
## Treatment8 pregnant 1.37939 2.97574 0.464 0.64383
## Treatment8 virgin -22.83507 2.98667 -7.646 6.13e-12 ***
## Treatmentnone -2.85763 2.97135 -0.962 0.33815
## Thorax 136.79061 12.42682 11.008 < 2e-16 ***
## Sleep -0.07738 0.05985 -1.293 0.19854
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.48 on 118 degrees of freedom
## Multiple R-squared: 0.6612, Adjusted R-squared: 0.6439
## F-statistic: 38.38 on 6 and 118 DF, p-value: < 2.2e-16On a un effet positif du thorax sur la longévité et un effet du traitement sur la longévité. Les mâles exposés à des femelles réceptives vivent moins longtemps que les mâles moins ou pas exposés et ceci est particulièrement vrai chez les mâles fortement exposés. On pourrait également comparer les traitements entre eux avec un test de comparaison multiple. Puisque nous avons utilisé la fonction lm nous ne pouvons pas prendre TukeyHSD qui ne fonctionne qu’avec aov, mais nous pouvons prendre la fonction plus générale glht dans le package multcomp.
library(multcomp)## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geysersummary(glht(m, linfct = mcp(Treatment = "Tukey")))##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = Longevity ~ Treatment + Thorax + Sleep, data = d)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 virgin - 1 pregnant == 0 -9.551 2.969 -3.217 0.014158 *
## 8 pregnant - 1 pregnant == 0 1.379 2.976 0.464 0.990400
## 8 virgin - 1 pregnant == 0 -22.835 2.987 -7.646 < 1e-04 ***
## none - 1 pregnant == 0 -2.858 2.971 -0.962 0.871640
## 8 pregnant - 1 virgin == 0 10.931 2.991 3.655 0.003451 **
## 8 virgin - 1 virgin == 0 -13.284 3.013 -4.409 0.000223 ***
## none - 1 virgin == 0 6.694 2.975 2.250 0.168841
## 8 virgin - 8 pregnant == 0 -24.214 2.976 -8.135 < 1e-04 ***
## none - 8 pregnant == 0 -4.237 2.998 -1.413 0.620259
## none - 8 virgin == 0 19.977 2.998 6.664 < 1e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)ce sont surtout les traitements où les mâles ont été mis en présence de femelles réceptives qui ont eu un effet négatif sur la longévité. Du côté des conditions d’applications de la régression, on note une certaine hétéroscédasticité, i.e. que la variance change selon les valeurs prédites. Ceci est possiblement un problème, mais nous assumerons que les conditions ne dévient pas trop de l’homoscédasticité. Par contre, on peut voir que les résidus ont une distribution normale comme en témoigne le QQ-plot.
par(mfrow = c(2, 2))
plot(m)
Pour mieux illustrer les prédictions de notre modèle, utilisons un graphique à bande avec des intervalles de confiance sur les prédictions. Je vais vous laisser étudier le code. Notez l’utilisation d’une boucle à la fin du code. Ici, i prendra toute les valeurs correspondant aux lignes du tableau de prédictions.
# on produit un tableau de valeurs à prédire
dat <- with(d, data.frame(Treatment = unique(Treatment), Thorax = mean(Thorax),
Sleep = mean(Sleep)))
# on organise nos valeurs dans un ordre logique
dat <- dat[c(2, 3, 4, 1, 5), ]
# on prédit nos valeurs avec un intervalle de confiance
p <- predict(m, newdata = dat, interval = "confidence")
# on produit un barplot
b <- barplot(p[, 1], names.arg = dat$Treatment, col = "green4", border = NA,
ylim = c(0, max(p)), ylab = "Longevity", xlab = "Treatment")
# on ajoute les bornes inférieures et supérieurs des intervalles de
# confiance
for (i in 1:nrow(dat)) {
lines(rep(b[i, 1], 2), c(p[i, 2], p[i, 3]))
}
Tendances dans la phénologie de nidification
Pour ce deuxième exercice, on veut savoir s’il y a des changements dans les dates de pontes chez une espèce d’oiseau (la Paruline jaune). Voici la base de données warblers.csv. Est-ce que ces changements sont les mêmes d’une province à l’autre? La colonne jul donne le jour julien de la date de ponte du premier oeuf (nb de jour depuis le 1er janvier. Dans cet exercice, vous devez:
Explorez les données
Déterminer s’il y a des changements dans les dates de ponte au cours des dernières années et si ces changements sont différents d’une province à l’autre.
Illustrez les prédictions de votre modèle en utilisant la fonction
predict. Pour ce faire, vous devrez utiliser l’argumentnewdataavec lequel vous soumettrez des observations à votre modèle.Est-ce que les suppositions de base sont respectées? Est-ce que vous voyez des choses à améliorer avec cette analyse?
Solution
Chargeons d’abord les données et explorons les données à l’aide d’un graphique.
d <- read.csv("https://raw.githubusercontent.com/frousseu/ECL749/master/data/warblers.csv")
plot(d$year, d$jul, type = "n")
points(d$year[d$prov == "ON"], d$jul[d$prov == "ON"], col = "darkred")
points(d$year[d$prov == "QC"], d$jul[d$prov == "QC"], col = "darkblue")
points(d$year[d$prov == "NB"], d$jul[d$prov == "NB"], col = "green4")
legend("topright", pch = 1, col = c("darkred", "darkblue", "green4"), legend = c("ON",
"QC", "NB"), bty = "n")
Pour déterminer si les changements dans la date de ponte dépendent de la province, il faut utiliser une interaction entre l’année et la province. Pour ce faire, on peut écrire notre modèle de cette façon.
m <- lm(jul ~ year * prov, data = d)
summary(m)##
## Call:
## lm(formula = jul ~ year * prov, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.172 -7.785 -1.991 5.844 53.092
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 811.2265 200.6610 4.043 5.47e-05 ***
## year -0.3270 0.1012 -3.232 0.001251 **
## provON -576.8370 203.5185 -2.834 0.004637 **
## provQC -954.5381 267.2647 -3.572 0.000363 ***
## year:provON 0.2858 0.1026 2.784 0.005410 **
## year:provQC 0.4778 0.1350 3.540 0.000409 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.22 on 2079 degrees of freedom
## Multiple R-squared: 0.07336, Adjusted R-squared: 0.07114
## F-statistic: 32.92 on 5 and 2079 DF, p-value: < 2.2e-16On a effectivement une interaction significative entre la province et l’année ce qui nous indique que les changements dans la date de ponte ne sont pas les mêmes d’une province à l’autre. Cependant, pour bien interpréter cette interaction, nous utiliserons un graphique de prédictions afin d’éviter la gymnastique difficile du calcul avec les coefficients du modèle.
dat <- expand.grid(year = min(d$year):max(d$year), prov = c("ON", "QC", "NB"))
p <- predict(m, dat)
dat$p <- p
plot(d$year, d$jul, type = "n")
points(d$year[d$prov == "ON"], d$jul[d$prov == "ON"], col = "darkred")
points(d$year[d$prov == "QC"], d$jul[d$prov == "QC"], col = "darkblue")
points(d$year[d$prov == "NB"], d$jul[d$prov == "NB"], col = "green4")
legend("topright", pch = 1, col = c("darkred", "darkblue", "green4"), legend = c("ON",
"QC", "NB"), bty = "n")
lines(p ~ year, data = dat[dat$prov == "ON", ], col = "darkred", lwd = 2)
lines(p ~ year, data = dat[dat$prov == "QC", ], col = "darkblue", lwd = 2)
lines(p ~ year, data = dat[dat$prov == "NB", ], col = "green4", lwd = 2)
Du côté des conditions d’applications, les résidus ne suivent clairement pas une distribution normale. Il n’est toutefois pas évident de déterminer ce qui devrait être fait pour corriger la situation. Une option serait peut être d’utiliser la régression quantile, qui comporte moins de conditions et qui permet d’étudier la réponse à des quantiles spécifiques au lieu de la régression qui permet d’étudier la réponse moyenne. Un autre problème évident est que nous ne prenons pas en compte la possible corrélation temporelle entre les observations. Il faudrait également ajouter des intervalles de confiance à nos courbes de prédictions.
Atelier 4
Pour cet atelier, nous nous intéresserons aux variables affectant la présence d’une espèce de grenouilles (Southern Corroboree frog) retrouvée en Australie. Pour ce faire, nous avons une base de données frogs.csv avec différentes variables explicatives. Pour cet exercice, vous devrez déterminer quelles sont les variables pertinentes à utiliser pour comprendre ou prédire l’occurence de l’espèce et vous devrez illustrer les prédictions de votre modèle à l’aide de la fonction predict et de graphiques ou encore en utilisant le package visreg. Pour mieux comprendre l’utilisation de ce package, le plus simple est de se référer à la vignette expliquant comment le package est utilisé. N’oubliez pas de faire un peu d’explorations de données avant de sélectionner vos variables explicatives.
- pres.abs: 0 = espèce absente, 1 = espèce présente
- northing: coordoonées nord
- easting: coordonnées est
- altitude: altitude en mètres
- distance: distance à la population ou à la mention de présence connue la plus près
- NoOfPools: nombre de mares de reproduction potentielles
- avrain: quantité de pluie moyenne pendant le printemps
- meanmin: température moyenne minimale pendant le printemps
- meanmax: température moyenne maximale pendant le printemps
Solution
La première étape consiste à explorer les différentes valeurs de nos variables. On peut faire cela avec la fonction summary.
d <- read.csv("https://raw.githubusercontent.com/frousseu/ECL749/master/data/frogs.csv",
sep = ";")
summary(d)## pres.abs northing easting altitude
## Min. :0.0000 Min. : 84.0 Min. : 673.0 Min. :1280
## 1st Qu.:0.0000 1st Qu.:192.0 1st Qu.: 977.8 1st Qu.:1480
## Median :0.0000 Median :222.5 Median :1023.0 Median :1580
## Mean :0.3726 Mean :228.2 Mean :1004.6 Mean :1547
## 3rd Qu.:1.0000 3rd Qu.:290.0 3rd Qu.:1086.2 3rd Qu.:1625
## Max. :1.0000 Max. :335.0 Max. :1222.0 Max. :1800
## distance NoOfPools avrain meanmin
## Min. : 250 Min. : 1.00 Min. :124.7 Min. :2.033
## 1st Qu.: 500 1st Qu.: 8.00 1st Qu.:141.7 1st Qu.:2.567
## Median : 1000 Median : 18.00 Median :148.8 Median :3.000
## Mean : 1933 Mean : 25.11 Mean :148.1 Mean :3.120
## 3rd Qu.: 2000 3rd Qu.: 32.00 3rd Qu.:155.0 3rd Qu.:3.567
## Max. :18000 Max. :232.00 Max. :198.3 Max. :4.333
## meanmax
## Min. :11.60
## 1st Qu.:12.97
## Median :13.38
## Mean :13.67
## 3rd Qu.:14.21
## Max. :15.97Par la suite, il serait préférable d’inspecter les liens entre les différentes variables explicatives potentielles. On peut s’attendre à ce que certaines soient fortement corrélées (i.e. meanmin et meanmax). Pour illustrer toutes les paires de variables, on peut utiliser les fonctions plot ou pairs appliquées à notre jeu de données. On peut également éliminer la variable réponse et les coordonnées puisque celles-ci ne feront pas partie de nos variables explicatives. Ces trois variables se trouvent dans les 3 premières colonnes de notre tableau et on peut les éliminer en utilisant l’indexage négatif.
plot(d[, -(1:3)])
À première vue, on a une corrélation très forte entre certaines variables, notamment entre l’altitude et les températures minimales et maximales. Il serait donc préférable de choisir l’une ou l’autre de ces variables pour éviter la collinéarité et la redondance dans notre modèle. On pourrait également utiliser la fonction cor pour calculer la corrélation entre chacune de nos paires de variables.
cor(d[, -(1:3)])## altitude distance NoOfPools avrain meanmin
## altitude 1.0000000 0.1820645 0.2667913 0.7780271 -0.9536610
## distance 0.1820645 1.0000000 -0.0959017 0.1484831 -0.3263706
## NoOfPools 0.2667913 -0.0959017 1.0000000 0.2437792 -0.2237348
## avrain 0.7780271 0.1484831 0.2437792 1.0000000 -0.6492240
## meanmin -0.9536610 -0.3263706 -0.2237348 -0.6492240 1.0000000
## meanmax -0.9965570 -0.2045681 -0.2633172 -0.8186997 0.9462741
## meanmax
## altitude -0.9965570
## distance -0.2045681
## NoOfPools -0.2633172
## avrain -0.8186997
## meanmin 0.9462741
## meanmax 1.0000000On a également une corrélation importante entre la pluie et l’altitude et la température maximale. Encore une fois, il pourrait être important de choisir l’une ou l’autre de ces variables. Bien que la température est certainement très importante pour les amphibiens, l’altitude est selon moi plus utile comme variable si l’intention est de faire de la prédiction sur la présence de l’espèce. En effet, cette variable est probablement plus facilement accessible à l’échelle du paysage comparativement à la température.
Maintenant, on peut également se poser la question quelle est la relation entre les variables explicatives et la probabilité de présence? Pour ce faire, on peut encore une fois utiliser la fonction plot avec une formule. L’utilisation du point . indique que nous voulons considérer toutes les variables du tableau de données.
par(mfrow = c(2, 3))
plot(pres.abs ~ ., data = d[, -(2:3)])
Ce n’est pas si évident sur les graphiques en raison des valeurs binaires, mais il semble parfois y avoir des probabilité de présence plus élevées pour des valeurs intermédiaires de nos différentes variables environnementales. Si on y réfléchit un peu, ceci n’est pas nécessairement surprenant. Imaginons une espèce fictive ayant une aire de distribution donnée et une variable climatique comme la température moyenne annuelle. Il y aura probablement une température minimale en dessous de laquelle l’espèce ne sera jamais présente ainsi qu’une température maximale au delà de laquelle l’espèce ne sera jamais rencontrée. Si on fait des inventaires à grande échelle pour trouver cette espèce, on obtiendra probablement que l’espèce est davantage rencontrée pour des valeurs intermédiaires de température moyenne annuelle. Il est possible qu’on trouve une relation linéaire positive ou négative entre la température et la présence de l’espèce, mais si on couvre un gradient de température assez étendu, on obtiendra nécessairement que l’espèce est moins abondante ou absente pour les valeurs extrêmes de température de la superficie inventoriée. Dans notre cas, il se pourrait que notre espèce de grenouille se retrouve davantage à des altitudes intermédiaires ce qui pourrait donner lieu à une relation quadratique. Il pourrait donc être judicieux d’inclure des relations quadratiques dans notre modèle.
A priori, il ne me semble pas y avoir d’interactions évidentes à inclure dans le modèle. Puisque les variables climatiques sont très fortement corrélées, je ne sélectionnerais que l’altitude. Les deux autres variables me semblent toutes deux pertinentes à inclure. On s’attend à ce que le nombre de mares NoOfPools influence positivement la probabilité de présence. Pour ce qui est de la distance, on s’attend à la même relation. En effet, il est raisonnable de penser que plus on est près d’une présence connue, plus on risque de rencontrer l’espèce. L’effet de cette variable est plus à considérer comme étant un indicateur des conditions locales que comme une variable influençant directement la présence de l’espèce. Par exemple, on a de l’information partielle sur la présence de l’espèce dans une région donnée où on veut raffiner nos connaissances de distribution de l’espèce. Elle pourrait également être liée au potentiel de dispersion de l’espèce dans un context de méta-population. Un modèle potentiel serait celui-ci:
m <- glm(pres.abs ~ NoOfPools + distance + altitude + I(altitude^2), family = "binomial",
data = d)
summary(m)##
## Call:
## glm(formula = pres.abs ~ NoOfPools + distance + altitude + I(altitude^2),
## family = "binomial", data = d)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9904 -0.8161 -0.3054 0.8766 2.8563
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.446e+01 2.793e+01 -3.382 0.000719 ***
## NoOfPools 3.180e-02 9.482e-03 3.353 0.000799 ***
## distance -6.221e-04 1.642e-04 -3.789 0.000151 ***
## altitude 1.332e-01 3.806e-02 3.500 0.000466 ***
## I(altitude^2) -4.657e-05 1.291e-05 -3.607 0.000309 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 279.99 on 211 degrees of freedom
## Residual deviance: 209.25 on 207 degrees of freedom
## AIC: 219.25
##
## Number of Fisher Scoring iterations: 6On peut voir que l’abondance de mares influence positivement la présence de l’espèce, alors que la distance à une présence connue a un effet négatif sur la probabilité de présence de l’espèce. Du côté de l’altitude, il est très difficile de décrire quel est l’effet sur la probabilité de présence sans utiliser de graphiques. On peut toutefois dire que l’effet quadratique semble justifié comme en témoigne la valeur du coefficient du terme quadratique qui est significativement différent de 0. Remarquez l’utilisation de la fonction I qui facilite la formulation du modèle et qui permet d’indiquer de la variable altitude est à la base des deux termes (linéaires et quadratiques). Cela peut nous éviter d’avoir à calculer plus de valeurs lorsque nous utiliserons la fonction predict.
Visualisons cet effet avec la fonction predict. En premier, on crée un data.frame avec des valeurs que nous allons soumettre à notre modèle. Pour ce faire, on couvre l’étendue de nos valeurs d’altitude et on fixe les autres variables à leur valeur moyenne.
# on fait varier l'altitude de la valeur min à la valeur max par bond de 10
alt <- seq(min(d$altitude), max(d$altitude), by = 10)
# on crée un newdata avec les autres variables fixées à leur valeur moyenne
dat <- data.frame(altitude = alt, NoOfPools = mean(d$NoOfPools), distance = mean(d$distance))
# on sort les prédictions du modèle en précisant que nous les voulons sous
# l'échelle de notre réponse
p <- predict(m, newdata = dat, type = "response")
# on illustre les résultats
plot(d$altitude, d$pres.abs)
lines(alt, p)
La probabilité de présence semble atteindre un pic aux alentours de 1400m. Il faut toutefois rester prudent dans l’interprétation précise de ces valeurs, car nous n’avons pas d’intervalles de confiance autour de ces courbes. On pourrait les ajouter en utilisant l’argument se.fit de predict, mais cela demanderait un peu plus de travail. Nous n’avons d’ailleurs pas encore visualisé l’effet de nos autres variables. Il faudrait produire ces graphiques pour chacune des variables. Il serait possible d’optimiser cette tâche en généralisant le code de prédictions, mais pour l’instant ceci donnerait un code un un peu plus difficile à suivre.
Pour faciliter la visualisation des résultats, nous pouvons faire appel au package visreg qui s’occupe de la majeure partie du travail. Pour comprendre son utilisation, lisez attentivement la vignette du package. Les deux premières pages vous indiquent de façon assez simple comment créer des graphiques de visualisation d’un modèle linéaire et même comment visualiser une interaction! Le package a même un site web pour décrire les fonctionnalités. Il y a notamment une section sur les GLM. Pour ce package, l’argument qui permet de spécifier que l’on veut notre prédiction sous l’échelle de la réponse s’appelle scale contrairement à type pour la fonction predict.
library(visreg)
par(mfrow = c(2, 2), mar = c(4, 3, 1, 1))
visreg(m, scale = "response", ylim = 0:1)
Certains ont utilisé les coordonnées comme variables explicatives dans leur modèle. Ceci n’est pas nécessairement mauvais, mais cela mérite une bonne réflexion. Est-ce que la relation entre la présence et les coordonnées est linéaire ou est-elle plus complexe? Est-ce que la relation entre les coordonnées et la probabilité de présence pourrait servir à prédire la présence de l’espèce dans d’autres régions? Probablement pas. Si le but est de prédire la présence de l’espèce seulement à l’intérieur de la zone étudiée, les coordonnées pourraient être utilisées. Cependant, il faudrait s’assurer que la relation est linéaire, ce qui serait assez suprenant. Si c’est le cas, il faudrait plutôt se demander s’il n’y a pas une ou des variables non prises en compte qui pourraient créer un tel gradient de probabilité de présence dans l’aire d’étude. Si on s’intéresse à comprendre l’influence des variables environnementales sur l’occurence des grenouilles à l’extérieur de la zone étudiée, la relation observée avec les coordonnées devient probablement inutile.
Malheureusement, ou plutôt heureusement, il n’y a pas nécessairment une meilleure ou une unique réponse à ce défi d’analyse de données. La mienne n’est qu’une suggestion et vous avez peut être de meilleures suggestions. L’important est de justifier nos choix et de prendre en considération les éléments les plus importants. Un élément important que j’ai ignoré est la potentielle corrélation spatiale entre les observations. Par exemple, est-ce que les observations près les unes des autres peuvent être considérées comme étant indépendantes? Aussi, est-ce qu’il y aurait d’autres variables à prendre en compte dans cette analyse? Comment les inventaires ont été réalisés? Est-ce qu’il faudrait tenir compte de la probabilité de détection de l’espèce?
Atelier 5
Les parasites chez le Colibri à gorge rubis
Pour cet exercice, vous étudierez ce qui influence la charge parasitaire chez le Colibri à gorge rubis (Archilochus colubris), une espèce retrouvée au Québec. Dans la base de données colibris.csv, vous retrouverez différentes variables comme le sexe (4 = mâle, 5 = femelle) et l’âge (1 = adulte, 2 = jeune). L’exercice consiste à déterminer quelles sont les variables qui affectent la quantité de parasites. Pour interpréter vos résultats, illustrez graphiquement votre modèle.
Solution
On importe d’abord les données et on affiche notre data.frame et un petit résumé de ce qu’il contient.
d <- read.csv("https://raw.githubusercontent.com/frousseu/ECL749/master/data/colibris.csv")
head(d)## bague sexe age masse parasites
## 1 E09801 4 1 3.2 0
## 2 E09802 4 2 3.0 0
## 3 E09808 4 1 3.5 3
## 4 E09809 4 1 3.6 0
## 5 E09810 4 1 3.6 0
## 6 E09812 5 1 3.9 0summary(d)## bague sexe age masse
## E09819 : 4 Min. :4.000 Min. :1.000 Min. :2.500
## E09822 : 3 1st Qu.:4.000 1st Qu.:1.000 1st Qu.:2.900
## E09874 : 3 Median :4.000 Median :1.000 Median :3.200
## E10196 : 3 Mean :4.466 Mean :1.204 Mean :3.273
## E10203 : 3 3rd Qu.:5.000 3rd Qu.:1.000 3rd Qu.:3.600
## E10207 : 3 Max. :5.000 Max. :2.000 Max. :5.200
## (Other):354 NA's :1
## parasites
## Min. : 0.000
## 1st Qu.: 0.000
## Median : 1.000
## Mean : 4.491
## 3rd Qu.: 7.000
## Max. :40.000
## NA's :4On voit quelques zéros dans la colonne parasites et on peut donc se demander à quel type de données nous avons à faire. On peut illustrer les valeurs de parasites avec un histogramme.
hist(d$parasites)
On a ici un nombre d’individus, avec des valeurs allant de 0 à environ 40. Il faudra donc probablement modéliser l’abondance de parasites avec une distribution de poisson. On a également plusieurs colibris dans la base de données. Une deuxième question à se poser est est-ce qu’on a plusieurs comptes de parasites pour un même individus? On peut vérifier ceci en utilisant la fonction table qui va comptabiliser le nombre de fois que chaque numéro de bague apparaît dans la base de données. En utilisant une seconde fois la fonction table, on compte le nombre de fois où chacune de ces valeurs est observée.
table(table(d$bague))##
## 1 2 3 4
## 243 48 10 1On voit que pour la majorité des colibris, on a une seule valeur, mais certains ont plusieurs valeurs (2, 3 ou 4 observations). On devrait donc par défaut se diriger vers les modèles mixtes pour traiter l’individu comme un effet aléatoire puisque les comptes de parasites sur un même individu ne sont pas indépendants. Par contre, puisque la majorité des individus ont une seule observation, cela peut créer des problèmes dans le cadre des modèles mixtes et on pourrait plutôt opter pour effectuer une moyenne ou pour une sélection aléatoire d’une seule observation pour les individus ayant plusieurs observations. On pourrait également comparer les deux approches et déterminer si nos conclusions par rapport à l’effet de nos variables sont fortement affectées par notre méthode d’analyse. Si ce n’est pas le cas, on peut choisir la méthode qui nous convient en justifiant notre choix et en exposant le fait que nos conclusions ne sont pas affectées par la méthode choisie.
On devrait également explorer les liens entre nos variables explicatives. On peut faire ceci avec des boxplot ou des graphiques standards.
par(mfrow = c(1, 2))
boxplot(masse ~ age, data = d)
boxplot(masse ~ sexe, data = d)
table(d$sexe, d$age)##
## 1 2
## 4 136 63
## 5 161 13Les femelles sont plus massives que les mâles et les adultes semblent également plus massifs que les jeunes. Cependant, il semble y avoir une disproportion de mâles par rapport aux femelles chez les jeunes.
Avant de lancer nos modèles, il faut transformer le sexe et l’âge en facteur étant donné que ces deux variables sont représentées par des valeurs numériques.
d$sexe <- as.factor(d$sexe)
d$age <- as.factor(d$age)
library(lme4)## Loading required package: Matrixm1 <- glm(parasites ~ sexe + age + masse, data = d, family = poisson)
m2 <- glmer(parasites ~ sexe + age + masse + (1 | bague), data = d, family = poisson)
summary(m1)##
## Call:
## glm(formula = parasites ~ sexe + age + masse, family = poisson,
## data = d)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.840 -2.603 -1.361 1.052 8.369
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.98469 0.25726 11.602 < 2e-16 ***
## sexe5 -0.46180 0.08081 -5.715 1.10e-08 ***
## age2 -0.19144 0.06306 -3.036 0.0024 **
## masse -0.39478 0.08789 -4.492 7.07e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2920.5 on 367 degrees of freedom
## Residual deviance: 2709.8 on 364 degrees of freedom
## (5 observations deleted due to missingness)
## AIC: 3472.8
##
## Number of Fisher Scoring iterations: 6summary(m2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: parasites ~ sexe + age + masse + (1 | bague)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 2065.3 2084.9 -1027.7 2055.3 363
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8667 -0.6979 -0.1469 0.1926 4.9784
##
## Random effects:
## Groups Name Variance Std.Dev.
## bague (Intercept) 2.612 1.616
## Number of obs: 368, groups: bague, 299
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.7788 0.6460 5.850 4.93e-09 ***
## sexe5 -0.1041 0.2626 -0.396 0.692
## age2 0.1008 0.1566 0.644 0.520
## masse -1.0082 0.2177 -4.632 3.62e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) sexe5 age2
## sexe5 0.409
## age2 0.115 0.243
## masse -0.970 -0.548 -0.202On obtient des résultats bien différents selon nos deux modèles. On a un effet négatif de la masse, mais pour l’âge et le sexe nos conclusions sont très différentes. Puisque nous avons une seule observation pour la majorité des individus, il est peut être préférable de ne pas utiliser de modèles mixtes qui requièrent plusieurs observations pour chaque niveau.
On pourrait opter pour un GLM et ne prendre qu’une seule observation aléatoirement par individu. On peut faire ceci en utilisant une combinaison de split, lapply et sample. Je vous laisse étudier le code en vous donnant quelques indications sur chaque étape.
s <- split(d, d$bague) # on sépare les données par individu dans une liste
l <- lapply(s, function(i) {
# on applique une fonction sur chaque data.frame de la liste on prend la
# ligne s'il n'y en qu'une
if (nrow(i) == 1) {
i
} else {
i[sample(1:nrow(i), 1), ] # ou on choisit aléatoirement une ligne s'il y en a plusieurs
}
})
d2 <- do.call("rbind", l) # on fusionne les différents subsets de la liste ensemble dans un nouveau data.frameDans l’objet d2, on a maintenant une seule ligne par individu.
m1sub <- glm(parasites ~ sexe + age + masse, data = d2, family = poisson)
summary(m1sub)##
## Call:
## glm(formula = parasites ~ sexe + age + masse, family = poisson,
## data = d2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.905 -2.509 -1.275 1.203 8.237
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.00779 0.28998 10.372 < 2e-16 ***
## sexe5 -0.57253 0.08878 -6.449 1.13e-10 ***
## age2 -0.19273 0.06595 -2.922 0.00348 **
## masse -0.39051 0.09831 -3.972 7.12e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2356.5 on 297 degrees of freedom
## Residual deviance: 2153.0 on 294 degrees of freedom
## (4 observations deleted due to missingness)
## AIC: 2784.2
##
## Number of Fisher Scoring iterations: 6Nos conclusions pour ce glm modifié sont très similaires. On conclut donc que les mâles et les jeunes sont davantage infectés par les parasites. Plusieurs interprétations sont possibles, mais aucune ne me semble évidente et il faudrait chercher dans la littérature pour valider ces interprétations. Quant à la masse, elle a un effet négatif sur le nombre de parasites. Ceci pourrait s’expliquer par le fait que les individus plus massifs sont aussi en meilleure santé et donc plus aptes à se défendre contre les parasites.
par(mfrow = c(2, 2))
visreg(m1sub, scale = "response", ylim = c(0, 8))
La richesse spécifique le long du littoral
Ici, on s’intéresse aux variables influençant le nb d’espèces rencontrées le long du littoral. On a la hauteur de la station d’échantillonnage par rapport à la hauteur moyenne de la marée ainsi que l’exposition qui est quantifiée à l’aide d’un index construit à partir de la “sévérité” des conditions de la station (force des vagues, longueur de la zone d’impact, pente, substrat allant de fin à grossier). La base de données se trouve à cet endroit plages.csv. Illustrez vos résultats.
Solution
Importons d’abord les données.
d <- read.csv("https://raw.githubusercontent.com/frousseu/ECL749/master/data/plages.csv")
head(d)## sample richesse exposition hauteur plage
## 1 1 11 10 0.045 1
## 2 2 10 10 -1.036 1
## 3 3 13 10 -1.336 1
## 4 4 11 10 0.616 1
## 5 5 10 10 -0.684 1
## 6 6 8 8 1.190 2Comme pour l’exercice précédent, on a une variable réponse qui correspond à un compte et qui est très asymmétrique à droite avec des valeurs de zéros.
hist(d$richesse)
range(d$richesse)## [1] 0 22On semble également avoir des comptes effectués sur plusieurs plages. Vérifions ceci.
table(d$plage)##
## 1 2 3 4 5 6 7 8 9
## 5 5 5 5 5 5 5 5 5On a cinq observations par plages et 9 plages. Il faudrait donc utiliser des modèles mixtes avec la plage en effet aléatoire. On devrait également explorer les liens entre les différentes variables. On peut inclure à la fois la richesse et les variables explicatives ce qui nous permettra d’un seul coup de vérifier la linéarité du lien avec les variables explicatives et la corrélation des variables explicatives.
plot(d[, c("richesse", "exposition", "hauteur")])
On remarque premièrement qu’il n’y a que trois valeurs d’exposition. On pourrait décider de traiter cette variable comme un facteur, mais puisque l’effet semble linéaire et qu’on peut s’imaginer que les valeurs représentent un gradient on pourrait décider de les traiter comme une variable numérique.
m <- glmer(richesse ~ hauteur + exposition + (1 | plage), data = d, family = poisson)
summary(m)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: richesse ~ hauteur + exposition + (1 | plage)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 211.0 218.2 -101.5 203.0 41
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8186 -0.5430 -0.2574 0.1204 3.9434
##
## Random effects:
## Groups Name Variance Std.Dev.
## plage (Intercept) 0.03485 0.1867
## Number of obs: 45, groups: plage, 9
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.25225 0.93145 6.712 1.92e-11 ***
## hauteur -0.50570 0.07252 -6.973 3.09e-12 ***
## exposition -0.44831 0.09279 -4.832 1.35e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) hauter
## hauteur 0.071
## exposition -0.995 -0.066On a des effets négatifs significatifs des deux variables ce qui signifie que la hauteur et l’exposition sont toutes deux associées à un moins grand nombre d’espèces. On peut visualiser l’effet de ces deux variables avec visreg.
par(mfrow = c(2, 3))
visreg(m, scale = "linear")
visreg(m, scale = "response")
La série du haut représente les prédictions du modèle sous l’échelle de la fonction de lien log et la série du bas représente les prédictions du modèle à l’échelle de la réponse. Contrairement aux GLMS, il n’y a pas d’intervalles de confiance puisque la meilleure façon de les calculer est incertaine. Ceci est expliqué dans la GLMM FAQ.
Atelier 6
Pour cet exercice, nous allons développer la débrouillardise sur les internets. Nous allons utiliser le jeu de données iris pour effectuer différentes tâches. Cette base de données contient des mesures morphométriques pour trois espèces d’iris.
Ordonnez la base de données selon l’espèce et puis selon la longueur des pétales.
Trouvez une fonction permettant de faire du clustering et tentez d’effectuer un groupement permettant de classifier les individus sur la base de leurs caractéristiques morphologiques. Est-ce que vous êtes capable d’identifier les trois espèces dans les groupes établis?
Effectuez une analyse en composantes principales afin de représenter la variation dans ce jeu de données et illustrez les résultats dans un graphique.
Solution
1.
On charge en premier les données.
d <- read.csv("https://raw.githubusercontent.com/frousseu/ECL749/master/data/iris.csv")
head(d)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.0 3.2 1.2 0.2 setosa
## 2 6.6 2.9 4.6 1.3 versicolor
## 3 4.8 3.0 1.4 0.1 setosa
## 4 7.3 2.9 6.3 1.8 virginica
## 5 4.7 3.2 1.3 0.2 setosa
## 6 5.5 3.5 1.3 0.2 setosaPour ordonner les données, on peut utiliser la fonction order. Pour illustrer son fonctionnement, on peut d’abord présenter un cas simple. Pour ce faire, créons un vecteur de trois éléments et appliquons la fonction.
x <- c("bonjour", "salut", "allo")
order(x)## [1] 3 1 2Pour chaque valeur, la fonction nous retourne à la position appropriée l’index de l’élément dans x qui devrait se retrouver à cette position. On peut donc utiliser le vecteur résultant de order pour réorganiser x avec ce nouvel indexage.
x <- x[order(x)]
x## [1] "allo" "bonjour" "salut"Pour un data.frame, order permet de trier en fonction de plusieurs colonnes successives. Il s’agit par la suite d’utiliser le résultat en indexant les lignes de notre data.frame.
d <- d[order(d$Species, d$Petal.Length), ]
head(d)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 104 4.6 3.6 1.0 0.2 setosa
## 23 4.3 3.0 1.1 0.1 setosa
## 1 5.0 3.2 1.2 0.2 setosa
## 144 5.8 4.0 1.2 0.2 setosa
## 5 4.7 3.2 1.3 0.2 setosa
## 6 5.5 3.5 1.3 0.2 setosaRemarquez que si nous avions utilisé la fonction attach, l’ordre des éléments dans notre nouveau tableau de données ne serait plus le même que dans les colonnes qui auraient été attachées précédemment. C’est une des raisons pour laquelle l’utilisation de attach est déconseillée.
2.
Ce qu’on veut ici, c’est une méthode qui nous permet de classifier les individus sur la base de leurs caractéristiques morphologiques. Sommes-nous capable de retrouver les espèces en les classifiant avec une quelconque méthode?
Pour trouver une méthode de classification, la façon la plus simple est de taper clustering in R dans Google. Le premier lien est le site Quick-R de DataCamp. On y trouve entre autres les fonctions kmeans et hclust. Testons la première. Le code est un peu compliqué puisque l’exemple cherche à identifier plusieurs groupes de façon successive. L’important dans ces exemples est la fonction kmeans. Pour son utilisation, on peut aller voir dans le fichier d’aide et dans ses exemples. Les deux arguments importants sont le tableau de variables morphologiques et le nombre de clusters désirés.
k <- kmeans(d[, -5], center = 3)
k## K-means clustering with 3 clusters of sizes 50, 38, 62
##
## Cluster means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 5.006000 3.428000 1.462000 0.246000
## 2 6.850000 3.073684 5.742105 2.071053
## 3 5.901613 2.748387 4.393548 1.433871
##
## Clustering vector:
## 104 23 1 144 5 6 43 69 79 109 115 3 9 19 51 52 66 73
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 76 102 124 131 148 149 13 17 18 27 46 50 71 85 106 130 140 142
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 145 30 61 77 91 99 111 122 87 114 136 141 64 108 29 28 38 37
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 3
## 65 7 134 25 47 63 150 20 95 101 107 126 34 55 70 58 59 89
## 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 119 117 129 22 118 139 146 31 32 53 83 94 103 135 2 96 98 35
## 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 72 112 116 127 80 128 40 137 8 39 82 92 138 26 75 110 33 42
## 3 3 3 3 3 3 2 3 2 3 3 3 3 3 3 3 3 3
## 86 10 11 21 44 68 97 120 81 123 49 54 84 113 24 45 74 12
## 3 3 2 3 2 3 3 3 2 2 2 2 2 2 2 2 2 2
## 16 56 57 88 147 48 78 100 121 125 133 105 143 36 67 41 62 90
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 4 132 15 14 60 93
## 2 2 2 2 2 2
##
## Within cluster sum of squares by cluster:
## [1] 15.15100 23.87947 39.82097
## (between_SS / total_SS = 88.4 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"Si vous étudiez le fichier d’aide (et les exemples) pour voir quelles valeurs sont retournées pas la fonction, vous verrez que les clusters sont retournés dans l’objet et qu’ils sont accessibles avec l’opérateur $. On peut donc les extraire et les comparer les clusters identifiés aux espèces à l’aide de la fonction table. On pourrait également les combiner dans un même tableau avec la fonction data.frame.
table(d$Species, k$cluster)##
## 1 2 3
## setosa 50 0 0
## versicolor 0 2 48
## virginica 0 36 14On voit que le groupe 1 semble très bien correspondre à l’espèce setosa et les deux autres groupes correspondent assez bien aux deux autres espèces, quoiqu’il y a un certain niveau de confusion entre les deux espèces. Évidemment, ici, on connaît les groupements à l’avance. Ceci est uniquement pour illustrer une méthode de classification qui permettrait éventuellement de classifier des sites ou des objets sur la base de variables quantitatives.
On pourrait également essayer la classification hiérarchique avec la fonction hclust et appliquer les différentes commandes utilisées dans l’exemple pour établir de grands groupes et illustrer les résultats.
dis <- dist(d[, -5], method = "euclidean") # on calcule la distance entre les objets
h <- hclust(dis) # on les classifie
plot(h, labels = d$Species, cex = 0.5) # on illustre les résultats en réduisant la taille du texte
groups <- cutree(h, k = 3) # on crée 3 grands groupes
rect.hclust(h, k = 3, border = "red") # on ajoute ces groupes sur l'arbre de classification
On obtient des résultats assez similaires à la méthode k-means, mais ici la classification est hiérarchique. Ceci n’est qu’un exercice de recherche de fonctions pour illustrer les méthodes de classification. Il faudrait évidemment faire des recherches plus approfondies pour identifier les meilleures méthodes adaptées à notre problème. Un point important ignoré lors de nos classification est qu’il aurait fallu s’assurer de standardiser (scaler) nos données avant la classification afin de donner le même poids à toutes les variables.
3.
Il existe plusieurs fonctions pour effectuer des PCA, mais nous allons utiliser la fonction rda du package vegan, un package classique pour les analyses multivariées. Remarquez que nous pouvons utiliser l’argument scale pour standardiser les variables afin de ne pas donner plus de poids à certaines variables dans l’analyse.
library(vegan)
m <- rda(d[, -5], scale = TRUE)
summary(m)##
## Call:
## rda(X = d[, -5], scale = TRUE)
##
## Partitioning of correlations:
## Inertia Proportion
## Total 4 1
## Unconstrained 4 1
##
## Eigenvalues, and their contribution to the correlations
##
## Importance of components:
## PC1 PC2 PC3 PC4
## Eigenvalue 2.9185 0.9140 0.14676 0.02071
## Proportion Explained 0.7296 0.2285 0.03669 0.00518
## Cumulative Proportion 0.7296 0.9581 0.99482 1.00000
##
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores: 4.940963
##
##
## Species scores
##
## PC1 PC2 PC3 PC4
## Sepal.Length 2.199 -0.89142 0.6810 -0.09290
## Sepal.Width -1.137 -2.18073 -0.2313 0.04392
## Petal.Length 2.450 -0.05785 -0.1345 0.28497
## Petal.Width 2.384 -0.15811 -0.6003 -0.18617
##
##
## Site scores (weighted sums of species scores)
##
## PC1 PC2 PC3 PC4
## 104 -0.655159 -0.193409 -0.349815 -0.055075
## 23 -0.621804 0.405731 -0.190395 0.053861
## 1 -0.521372 0.086972 0.237411 -0.473451
## 144 -0.519231 -0.784896 0.498009 -0.545836
## 5 -0.558311 0.144276 -0.046548 -0.079541
## 6 -0.482959 -0.279160 0.510263 -0.550396
## 43 -0.438794 0.986328 0.215078 -0.812492
## 69 -0.521320 -0.626044 0.005628 -0.529258
## 79 -0.603039 0.202168 -0.322000 0.186685
## 109 -0.573755 0.381540 -0.203515 0.027297
## 115 -0.539949 -0.186392 -0.036748 -0.299875
## 3 -0.523885 0.307482 0.243171 -0.006797
## 9 -0.517249 -0.206277 0.046563 -0.260934
## 19 -0.551323 0.470639 -0.152779 0.075287
## 51 -0.505256 -0.132170 0.285546 -0.236180
## 52 -0.564359 -0.272948 -0.016574 0.100692
## 66 -0.534807 -0.202559 0.134486 -0.067744
## 73 -0.577155 -0.020105 -0.353282 0.103082
## 76 -0.565430 0.162986 -0.146873 0.136884
## 102 -0.491417 0.284467 0.247065 -0.288729
## 124 -0.577663 -0.907551 0.087057 0.135146
## 131 -0.596827 -0.249932 -0.020468 0.382624
## 148 -0.488768 0.300046 0.067324 -0.393177
## 149 -0.520433 -0.003889 0.161154 -0.138390
## 13 -0.527285 -0.094163 0.093405 0.068988
## 17 -0.512107 -0.222443 0.217797 -0.028802
## 18 -0.526482 -0.421113 0.191129 0.041844
## 27 -0.521184 -0.389955 -0.168353 -0.167052
## 46 -0.512375 -0.113460 0.185222 -0.019754
## 50 -0.553384 -0.475923 -0.139672 0.105832
## 71 -0.542997 0.252085 -0.096137 0.184874
## 85 -0.432529 -0.178788 0.284827 -0.672360
## 106 -0.617518 -0.756842 -0.049736 0.642551
## 130 -0.515827 0.197911 0.267239 0.111838
## 140 -0.498269 0.194193 0.179315 -0.081352
## 142 -0.534220 -1.133542 -0.032147 -0.141647
## 145 -0.511572 -0.440410 0.282946 -0.046898
## 30 -0.563153 -0.472792 -0.060256 0.426705
## 61 -0.534939 0.142514 -0.072069 0.303510
## 77 -0.549314 -0.056156 -0.098737 0.374156
## 91 -0.460927 0.263995 0.321868 -0.122104
## 99 -0.463903 -0.199310 -0.326039 -0.496394
## 111 -0.505388 0.212903 0.078991 0.135074
## 122 -0.484377 -0.102187 -0.090949 -0.189708
## 87 -0.490158 -0.628394 -0.028400 -0.018523
## 114 -0.448346 -0.592881 0.394218 -0.171253
## 136 -0.429477 -0.036103 -0.036320 -0.423650
## 141 -0.452066 -0.172526 0.443660 -0.030613
## 64 -0.525943 -0.057918 -0.124258 0.757207
## 108 -0.504665 -0.481991 -0.261624 0.423377
## 29 -0.105716 0.651440 -0.200510 -0.560453
## 28 -0.115081 0.781677 -0.262499 0.113579
## 38 -0.085529 0.852066 -0.111439 -0.054857
## 37 -0.026023 1.119949 0.049275 -0.038572
## 65 -0.009507 0.446752 0.336540 -0.181602
## 7 -0.007877 0.185248 -0.205282 -0.305664
## 134 0.005538 0.663544 0.254377 0.091862
## 25 0.030887 0.659238 0.157947 0.026355
## 47 0.057039 0.327982 0.159240 -0.066295
## 63 0.038945 0.549667 0.182014 0.144991
## 150 0.002698 0.436329 -0.567511 0.079777
## 20 0.132741 0.744668 0.806428 -0.128185
## 95 0.112475 0.176098 0.279019 -0.318333
## 101 0.066943 0.560940 -0.094157 -0.024963
## 107 0.096227 0.740313 0.024328 -0.184351
## 126 0.079471 0.417081 0.209976 -0.018305
## 34 0.060627 0.252700 -0.096757 0.164318
## 55 0.037503 0.334244 0.318073 0.575452
## 70 0.016433 0.092625 -0.307060 0.412448
## 58 0.068149 0.361096 -0.137839 0.301050
## 59 0.021575 0.076458 -0.135826 0.644580
## 89 0.053775 0.162426 -0.164507 0.371696
## 119 0.104069 0.026709 -0.215962 -0.112474
## 117 0.165933 0.026758 0.469708 -0.121814
## 129 0.136114 0.065353 0.286073 0.055670
## 22 0.246665 0.583585 0.724840 -0.383552
## 118 0.065904 0.472623 -0.099504 0.758656
## 139 0.206727 -0.214812 0.530244 -0.294159
## 146 0.206459 -0.105828 0.497670 -0.285111
## 31 0.082712 0.082838 -0.516935 0.536803
## 32 0.156992 0.095336 -0.090423 0.102141
## 53 0.172705 -0.250912 0.099119 -0.013745
## 83 0.101341 -0.356814 -0.474560 0.307421
## 94 0.052892 0.121433 -0.700570 0.714287
## 103 0.091788 0.250350 -0.130785 0.675053
## 135 0.289305 0.684545 0.507911 -0.633200
## 2 0.219123 -0.013599 0.627822 0.083753
## 96 0.147492 -0.010517 0.021569 0.413967
## 98 0.253972 0.087948 0.419400 -0.293579
## 35 0.176292 -0.326194 -0.156879 0.216869
## 72 0.260185 -0.364152 0.720933 -0.097640
## 112 0.247656 -0.220292 0.416800 -0.104297
## 116 0.149449 0.175705 0.307394 0.768643
## 127 0.169924 0.078582 0.072304 0.461957
## 80 0.174203 -0.167343 -0.649260 0.233447
## 128 0.296723 0.032600 0.765762 -0.111247
## 40 0.293056 -0.260062 0.583421 -0.026414
## 137 0.291450 0.393839 0.387972 0.027875
## 8 0.320786 -0.139805 0.274675 -0.187320
## 39 0.250577 0.267465 -0.110903 0.515666
## 82 0.086671 0.658913 -1.039290 0.373148
## 92 0.218397 -0.007268 -0.438957 -0.014683
## 138 0.277500 0.133510 -0.136837 -0.351555
## 26 0.241097 -0.027152 -0.355646 0.024258
## 75 0.314842 0.203312 0.005717 -0.392307
## 110 0.230949 0.241270 -0.872095 -0.077799
## 33 0.297575 0.490444 -0.611056 -0.277939
## 42 0.298437 0.720202 0.281741 0.182703
## 86 0.369474 0.378378 0.027771 -0.617201
## 10 0.273433 0.294905 -0.558308 0.113580
## 11 0.322149 -0.292325 -0.299723 -0.302334
## 21 0.273433 0.294905 -0.558308 0.113580
## 44 0.449105 -0.290983 -0.136981 -1.316566
## 68 0.263105 0.123606 0.193230 0.522324
## 97 0.226858 0.010267 -0.556295 0.457110
## 120 0.346583 0.186627 -1.057168 -0.772676
## 81 0.441718 -0.163290 -0.269880 -1.091092
## 123 0.359223 -0.113540 -0.189745 -0.334038
## 49 0.375495 -0.285358 -0.672325 -0.537795
## 54 0.378471 0.177947 -0.024417 -0.163504
## 84 0.324183 -0.426723 -0.982952 -0.073239
## 113 0.437359 -0.285308 0.013345 -0.547135
## 24 0.444881 -0.176912 -0.027737 -0.410403
## 45 0.317926 -0.178255 -0.190479 0.603829
## 74 0.347477 -0.107866 -0.039419 0.435393
## 12 0.340090 0.019828 -0.172318 0.660867
## 16 0.422317 0.079062 -0.285028 -0.087139
## 56 0.283963 0.342354 0.173469 1.372029
## 57 0.475796 -0.259044 -0.451074 -0.693853
## 88 0.372231 -0.450913 -0.996072 -0.099803
## 147 0.439876 0.075343 -0.372952 -0.280329
## 48 0.401627 -0.427849 -0.314296 0.172788
## 78 0.481205 -0.384193 -0.247266 -0.470770
## 100 0.471861 -0.442722 -0.665990 -0.599972
## 121 0.441082 -0.124501 -0.416633 0.045684
## 125 0.473877 0.300209 0.414909 0.242490
## 133 0.440099 -0.237272 0.753629 0.583630
## 105 0.481876 -0.366071 -0.356097 -0.126660
## 143 0.520772 -0.237154 0.213687 -0.165894
## 36 0.461512 -0.425256 0.442283 0.612005
## 67 0.435593 -0.367295 -1.058781 0.138048
## 41 0.575248 -0.109411 0.766458 0.050239
## 62 0.660927 -0.361549 0.571731 -0.829359
## 90 0.533644 -0.810616 -0.418659 -0.293865
## 4 0.543719 -0.177257 0.686307 0.667233
## 132 0.544307 -1.108240 0.519674 0.593330
## 15 0.649852 -0.337752 0.613226 0.284182
## 14 0.572975 -1.078844 0.134671 0.767485
## 60 0.684278 -0.174546 0.902945 0.356926
## 93 0.781818 -0.007503 0.740646 -0.126664La majeure partie de la variation dans le jeu de données est capturée par le premier axe (i.e. 73%). Afin d’interpréter les résultats, il est essentiel d’illustrer les résultats dans un graphique. La construction de graphiques élégants pour illustrer les ordinations n’est pas toujours évidente et les fonctions de bases ne donnent pas toujours les plus belles illustrations. Pour ce faire, vous pouvez vous inspirer de ce que vous trouvez sur internet. Voici un article de blog qui donne différentes méthodes pour améliorer l’illustration d’une ordination à partir de rda. Bien que cela représente plus de travail, sachez qu’il est généralement possible d’extraire les coordonées de toutes les valeurs utilisées afin de produire les graphiques, ce qui donne toute la flexibilité possible pour l’illustration des résultats. Ici, nous nous contenterons de graphiques relativement simples sur lesquels nous illustrerons les différentes espèces. Remarquez que nous utilisons un vecteur de couleurs que nous indexons par la suite en se basant sur le fait que la colonne d$Species est considérée comme un factor.
col <- c("tomato4", "orange", "darkblue")
biplot(m, scaling = 1, display = c("sites", "species"), type = c("text", "points"),
col = c("white", "black"), main = "Scaling 1")
points(m, scaling = 1, col = col[as.integer(d$Species)], pch = 1)
legend("bottomright", col = col, pch = 1, legend = levels(d$Species))
biplot(m, scaling = 2, display = c("sites", "species"), type = c("text", "points"),
col = c("white", "black"), main = "Scaling 2")
points(m, scaling = 2, col = col[as.integer(d$Species)], pch = 1)
legend("bottomright", col = col, pch = 1, legend = levels(d$Species))
L’axe 1 est surtout déterminé par les dimensions des pétales et dans une moindre mesure par la longueur des sépales. Quant au second axe, il est surtout déterminé par la largeur des sépales. En consultant le premier graphique, le groupe setosa semble se distinguer de façon assez marquée par rapport aux deux autres espèces. Avec le second graphique, on constate que les dimensions des pétales sont très corrélées, alors que ce n’est pas du tout le cas pour les dimensions des sépales.
Atelier 7
Quelles questions pouvez-vous vous poser avec le jeu de données aravo du package ade4? Est-ce que des variables environnementales permettent d’expliquer les variations en composition de plantes alpines (CCA)? Réalisez et illustrer quelques analyses multivariées afin d’explorer ce jeu de données. Vous pouvez vous inspirer du site GUSTA ME (donné en référence dans les liens utiles plus bas) et de son wizard sur les analyses contraintes ou encore de l’information contenue dans ce site Analysis of community ecology data in R.
Méthodes / résultats
Articles à lire pour le cours du 20 mars.
Travail de session
Le travail de session se déroulera en deux temps. La première partie (20%) consistera en l’établissement d’une méthodologie afin de répondre à une mise en situation et la seconde partie (80%) consistera en la rédaction d’un rapport sous forme d’article scientifique à partir d’une base de données en lien avec la première partie. Dans la plupart des cas, il y aura des différences entre la méthodologie proposée pour la mise en situation et la méthodologie utilisée pour récolter les données fournies. Dans ce cas, il faudra adapter votre méthodologie décrite dans la partie 2 du travail afin de refléter la méthodologie utilisée.
Partie I
La première partie (20%) devra comprendre la méthodologie proposée pour récolter et analyser les données. Plus spécifiquement, vous devrez présenter un plan d’échantillonnage en plus d’indiquer quelles mesures seront prises sur le terrain. Une ou des méthodes d’analyses statistiques devront également être proposées. Cette partie du travail devra être remise par courriel au plus tard le vendredi 23 mars à 16h. Pour ceux désirant utiliser leur propre jeu de données, vous devrez énoncer vous même votre propre mise en situation et m’en faire part d’ici vendredi le 16 mars à 16h en me montrant également la base de données qui sera utilisée (ou un extrait).
Comment-allez vous répondre à cette mise en situation? Quel sera votre plan d’échantillonnage? Quelles seront les variables que vous utiliserez? Quelles méthodes envisagez-vous utiliser pour analyser les données récoltées? La méthodologie proposée devrait pouvoir être exposée en un maximum de 2 pages. Puisque certains éléments sont inconnus, il sera bien sûr impossible de détailler tous les aspects de votre méthodologie. Le travail proposé ici a surtout pour objectif de vous faire réfléchir à un design d’étude vous permettant de répondre à une question de nature écologique.
Vous devrez probablement vous inspirer de la littérature scientifique afin d’identifier les variables à utiliser et établir votre plan d’échantillonnage. Afin de limiter l’éventail des possibilités, il serait préférable de vous limiter à identifier un maximum de 4 ou 5 variables à utiliser. En effet, de multiples variables pourraient être mesurées, mais il s’agit surtout ici d’identifier quelques variables pertinentes et surtout d’être capable de développer et d’exposer un plan d’échantillonnage pour répondre à une question donnée. De la même façon, il n’est pas nécessaire de mentionner que vous allez inventorier des milliers de sites pour répondre à la question. Il s’agit d’utiliser un nombre de stations ou un échantillonnage suffisant en fonction de la question ou du nombre de variables étudiées. Ceci n’est pas toujours évident à évaluer et vous devrez vous fier à votre jugement pour déterminer l’étendue de votre échantillonnage.
Partie II
Pour la seconde partie du travail (80%), chaque étudiant ou étudiante recevra une base de données en lien avec ce qui a été proposé dans la première partie du travail. Elle ou il devra réaliser les analyses, les interpréter et présenter les résultats dans un bref rapport écrit qui prendra la forme d’un article scientifique comprenant les éléments suivants :
- Introduction : Courte introduction présentant et mettant en contexte la question posée et les hypothèses biologiques.
- Méthodologie : Description des méthodes et des données récoltées et description et justification des analyses statistiques effectuées.
- Résultats : Présentation des résultats accompagnée de graphiques et/ou de tableaux.
- Discussion : Courte explication et interprétation des résultats et de leur signification biologique.
- Références : Un minimum de cinq références (articles scientifiques) devront être utilisées pour mettre en contexte la problématique et/ou interpréter les résultats.
- Annexe : Code R (propre!) utilisé pour effectuer les manipulations de données, les analyses statistiques et les graphiques. Ce code devrait contenir toutes les lignes de programmation permettant de reproduire votre travail à partir de votre base de données jusqu’aux graphiques présentés. Les lignes de codes inutiles, exploratoires ou générant des erreurs devraient être éliminées. Les parties plus compliquées de votre code devraient inclure des commentaires expliquant ce qui est effectué.
Votre travail au complet ne devrait pas prendre plus de 10 pages. Ce travail est à remettre au plus tard le vendredi 20 avril à 16h.
Mises en situation
Voici les différentes mises en situation. Si vous utilisez vos propres données, inspirez-vous de ces mises en situation pour formuler la vôtre.
1. Composition du paysage, Grand-duc d’Amérique et occurrence de la Chouette rayée
Quelles sont les caractéristiques du paysage affectant la présence de la Chouette rayée (i.e. une espèce de hibou, Strix varia) dans le sud du Québec? Est-ce que le Grand-duc d’Amérique (Bubo virginianus) a une influence sur sa présence?
Les données ont été récoltées en Montérégie et en Estrie. À chaque site, trois stations d’écoute d’une durée de 10 minutes ont été effectuées au courant de la saison propice à la détection. L’espèce était considérée présente lorsqu’il y avait au moins une détection. Voici la base de données chouette et une description des variables.
- site: station d’écoute
- year: année de visite
- rayon: rayon auquel les % ont été mesurés
- for_feuil: % de forêt feuillue
- for_con: % de forêt coniférienne
- for_mix: % de forêt mixte
- for_pert: % de forêt perturbée ou secondaire
- cult_int: % de culture intensive (maïs ou soya)
- cult_ext: % de culture extensive (pâturage, céréales, cultures maraîchères)
- GDApres: Grand-duc d’Amérique présent ou absent (1/0)
- CRpres: Chouette rayée présente ou absente (1/0)
2. Richesse ichtyologique dans les récif coraliens
Quelles sont les variables influençant la richesse en poissons dans les récifs coraliens?
Ces données ont été récoltées dans différentes sections de sites dans les Caraïbes. Le nombre d’espèces de poissons a été évalué à partir de 5 heures d’observations à chaque endroit dans des moments propices. Voici la base de données poissons_coraux et une description des variables:
- nbpoissons: nb d’espèces de poissons détectées
- nbcoraux: nb d’espèces de coraux détectées
- nbgrandscoraux: nb d’espèces de coraux de grande taille ( > 50cm )
- complexité: indice de complexité du relief variant de 1 (peu complexe) à 5 (très complexe) (estimation visuelle)
- pente: pente générale de la section de récif échantillonnée (degrés)
- site: site d’échantillonnage
3. Diversité aviaire, forêt et fragmentation en milieu tropical
Quel est l’effet de la fragmentation forestière sur la richesse aviaire en milieu tropical?
Ces données ont été récoltées dans différents sites fragmentés d’une région du Costa Rica. Le nombre d’espèces d’oiseaux à l’intérieur de chaque fragment a été évalué en visitant chaque parcelle pour un minimum de 4h chacune et jusqu’à ce que le nb de nouvelles espèces rencontrées n’augmente pas pendant une période de deux heures. Les visites ont été effectuées pendant la période de reproduction et à des moments propices (i.e. tôt le matin, conditions climatiques clémentes) pour l’observation. Voici la base de données fragmentation_aviaire et une description des variables:
- alpha: nb d’espèces d’oiseaux répertoriées
- alti: altitude du site
- super: superficie du fragment (hectares)
- frag: type de fragment (small ou large)
- foret: % de superficie en forêt dans un rayon de 1000m
- vieille: % de superficie en vieille forêt dans un rayon de 1000m
- site: numéro du site
4. L’abondance des moineaux dans un contexte agricole
Quels sont les éléments du paysage affectant l’abondance du Moineau domestique (Passer domesticus) en milieu agricole dans le sud du Québec?
Ces données ont été récoltées à diverses fermes de la Montérégie et de l’Estrie. Des stations d’écoute d’une durée de 10 minutes ont effectuées à chaque ferme tôt le matin pendant les mois de mai et juin. Voici la base de données moineaux et une description des variables:
- ferme: identifiant de la ferme
- vis: numéro de la visite
- abHOSP: nombre de moineaux détectés
- int500: % de culture intensive dans un rayon de 500m (maïs ou soya)
- ext500: % de culture extensive dans un rayon de 500m (pâturage, céréales, cultures maraîchères)
- for500: % de forêt dans un rayon de 500m
- eau500: % d’eau (surface) dans un rayon de 500m
- supbat500: % de bâtiments (surface) dans un rayon de 500m
- peribat500m: périmètre total de bâtiments dans un rayon de 500m (en mètres)
- perim500km: périmètre total de bâtiments dans un rayon de 500m (en km)
5. La diversité des papillons dans les prairies naturelles
Quelles sont les caractéristiques de la végétation affectant la diversité en espèces de papillons dans les prairies naturelles?
Ces données ont été récoltées à divers sites dans les prairies américaines. À chaque site, 5 transects de 100 m de récolte active (avec un filet) et 2 nuitée de récolte passive (avec un drap et de la lumière) ont été effectuées pendant les mois de juin et juillet. Voici la base de données papillons et une description des variables:
- site: site de récolte
- habitat: type d’habitat (champs de foin ou pâturage, graminées courtes, graminés grandes, mixtes)
- batiment: % de superficie en batiment
- urbain: % de superficie en zone urbaine
- year: année
- julien: jour julien
- nbsp: nb d’espèces de papillons
Liens utiles
Ici, je tente de regrouper différentes ressources qui vous seront presque assurément utiles lors de la réalisation de vos projets. Si vous découvrez des éléments qui seraient pertinents à ajouter, vous pouvez m’en faire part. Afin de ne pas “noyer le poisson”, j’essaie de ne mentionner que les ressources les plus importantes.
Ressources en ligne
R
Quick-R: site d’introduction à R
Base R Cheat Sheet: Guide mémoire pour les fonctions de base de R
RStudio Cheat Sheets: Divers guides mémoires pour R sur différents sujets
R Spatial: Introduction à R pour l’analyse et la modélisation de données spatiales
Statistiques
GLMM FAQ: Guide pour l’utilisation des modèles mixtes avec R
GUSTA ME: GUide to STatistical Analysis in Microbial Ecology. Un guide assez complet et très user friendly sur les analyses statistiques, surtout multivariées, comprenant l’explication, l’utilisation, l’interprétation, les suppositions et des exemples de plusieurs types d’analyses multivariées. Bien que le site soit orienté vers l’écologie microbienne, il est pertinent pour les biologistes et les écologistes en général.
Analysis of community ecology data in R: Un autre site d’introduction aux analyses multivariées très utile.
Ordination Methods for Ecologists: Site abordant plusieurs sujets liés à l’ordination.
Ateliers R du CSBQ: Ateliers en ligne du CSBQ concernant plusieurs sujets liés à R (manipulation de données, programmation, ggplot) ou aux statistiques (glm, modèles mixtes, analyses multivariées, etc.)
Ressources publiées
Articles
Generalized linear mixed models: a practical guide for ecology and evolution
A protocol for conducting and presenting results of regression-type analyses
A protocol for data exploration to avoid common statistical problems
Application of multivariate statistical techniques in microbial ecology
Livres
Plusieurs des livres suivants sont disponibles à la bibliothèque et certains sont disponibles en pdf sur internet.
The R Book: Un classique très complet pour l’introduction à R et aux analyses statistiques avec R.
Introductory statistics with R
Statistics for terrified biologists
La série des Zuur: La série de Highland Statistics est très orientée sur la pratique et couvre plusieurs analyases en écologie, surtout en ce qui a trait à la modélisation (GLM, GLMM, GAM, etc.)
Numerical Ecology with R: Les versions anciennes sont disponibles en ligne en verison pdf.
Numerical Ecology: Un classique sur plusieurs types d’analyses en écologie, particulièrement les analyses multivariées. Les versions anciennes sont disponibles en ligne en verison pdf.
Forums
Quelle que soit votre question sur R ou sur les statistiques, il est fort probable que celle-ci ait déjà été posée et répondue sur l’un de ces forums. Il s’agit d’utiliser les bons mots clés et de bien explorer les questions reliées.
Stack Overflow: Programmation avec R (questions portant le tag “R”).
Cross Validated: Questions concernant les statistiques.