Loading presentation...

Present Remotely

Send the link below via email or IM

Copy

Present to your audience

Start remote presentation

  • Invited audience members will follow you as you navigate and present
  • People invited to a presentation do not need a Prezi account
  • This link expires 10 minutes after you close the presentation
  • A maximum of 30 users can follow your presentation
  • Learn more about this feature in our knowledge base article

Do you really want to delete this prezi?

Neither you, nor the coeditors you shared it with will be able to recover it again.

DeleteCancel

Make your likes visible on Facebook?

Connect your Facebook account to Prezi and let your likes appear on your timeline.
You can change this under Settings & Account at any time.

No, thanks

CSBQ Atelier R 7

No description
by

CSBQ QCBS

on 6 September 2016

Comments (0)

Please log in to add your comment.

Report abuse

Transcript of CSBQ Atelier R 7

Call:
glm(formula = Faramea.occidentalis ~ Elevation, family = poisson,
data = faramea)

Deviance Residuals:
Min 1Q Median 3Q Max
-3.3319 -2.7509 -1.5451 0.1139 11.3995

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.7687001 0.1099136 16.092 < 2e-16 ***
Elevation -0.0027375 0.0006436 -4.253 2.11e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 414.81 on 42 degrees of freedom
Residual deviance: 388.12 on 41 degrees of freedom
AIC: 462.01
Spécifiez un modèle de régression logistique avec le jeu de données mites.csv
Centre de la Science de la Biodiversité du Québec

Série d'ateliers R

Atelier 7: Modèles linéaires généralisés (mixtes)
Site web:
http://qcbs.ca/wiki/r/atelier7

Comment marche un GLM Poisson?
spécifie des probabilités pour valeurs entières seulement
(ne comptez pas 7.56 arbres)
probabilité d'observer une valeur négative est nulle
(ne comptez pas -7 arbres)
la relation entre la moyenne et la variance permet de manipuler des données hétérogènes
(variance plus grande pour les grandes valeurs)
La distribution de Poisson
Distribution Poisson est appropriée pour modéliser ces données car :
Un GLM avec une distribution de Poisson considère que les valeurs Y ont été générées par une distribution de Poisson dont la moyenne et la variance sont égales à µ:
Donc seul le paramètre µ est nécessaire pour calculer la probabilité de différentes valeurs de Y; il suffit d'utiliser l'équation de la distribution Poisson:
Rappel : un GLM est composé de deux parties, le prédicteur linéaire et la fonction de lien.
Étape 1.

On suppose que les Yi suivent une distribution de Poisson de paramètre (moyenne et variance) µi
GLM Poisson sous R
Trois étapes
µi
correspond au nombre attendu d'individus
Comment marche un GLM Poisson?
Étape 2.

On spécifie le prédicteur linéaire comme dans un modèle linéaire
Étape 3.

La fonction de lien entre la moyenne de Yi et la partie systématique est un
log
ou
glm.poisson = glm(Faramea.occidentalis~Elevation, data=faramea, family=poisson)
Permet de spécifier le type de distribution et la fonction de lien (log)
summary(glm.poisson)
Un GLM Poisson se spécifie à l'aide de la fonction
glm()
Tout comme avec
lm()
vous pouvez accéder au résumé du modèle à l'aide de la fonction
summary()
Estimation des paramètres
La déviance
Dans notre modèle, les paramètres à estimer sont l'ordonnée à l'origine (alpha) et le coefficient de régression de l'élévation (Beta)
Le maximum de vraisemblance est utilisé pour estimer les paramètres du modèle
Résumé du modèle
La déviance résiduelle et approx. la différence entre la vraisemblance d'un modèle saturé (n paramètres pour chaque observation) et le modèle complet (p paramètres):
Res dev
= 2[log(L(saturé)) - log(
L
(complet)]
Donc,
Res dev
~ df = n - p =
Res df
...
ce qui n'est pas le cas ici:
388.12 >> 41
La surdispersion
Quand la déviance résiduelle est supérieure au nombre de degrés de liberté résiduels le modèle est surdispersé
Se produit lorsque la variance dans les données est plus grande que la moyenne. Dans ce cas la distribution de Poisson n'est plus appropriée (beaucoup de zéros, covariables manquantes, etc.)
Solutions
2.
Choisir une autre distribution: la
binomiale négative
1.
Corriger la surdispersion en utilisant un
GLM

quasi-Poisson
GLM Quasi-Poisson
La variance du modèle tient compte de la surdispersion en ajoutant le paramètre de surdispersion a l'équation :
sera estimé avant les paramètres. Corriger pour la surdispersion ne va pas affecter l'estimation des paramètres, mais leur significativité. En effet, les écarts-types des paramètres seront multipliés par

le
prédicteur linéaire
et la
fonction de lien
restent les mêmes
Certaines p-values marginalement significatives peuvent devenir non significatives!
La surdispersion
1
5
GLM Poisson
10
>15-20
GLM Quasi-Poisson
Binomiale négative
Ajuster un GLM quasi-Poisson sous R
glm.quasipoisson <- glm(Faramea.occidentalis~Elevation, data=faramea, family=quasipoisson)
Spécifie la nouvelle famille
Call:
glm(formula = Faramea.occidentalis ~ Elevation, family = quasipoisson,
data = faramea)

Deviance Residuals:
Min 1Q Median 3Q Max
-3.3319 -2.7509 -1.5451 0.1139 11.3995

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.768700 0.439233 4.027 0.000238 ***
Elevation -0.002738 0.002572 -1.064 0.293391
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 15.96936)

Null deviance: 414.81 on 42 degrees of freedom
Residual deviance: 388.12 on 41 degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 10
summary(glm.quasipoisson)
Les écart-types des paramètres sont multipliés par :
Mêmes estimé, mais :
Pas d'AIC
GLM Binomiale négative
Ajuster une binomiale négative sous R
Une distribution binomiale négative est favorisée quand la surdispersion est forte
La distribution binomiale négative n'est pas dans la fonction
glm()
,
il faut installer et charger la librairie MASS
install.packages("MASS")
glm.negbin <- glm.nb(Faramea.occidentalis~Elevation, data=faramea)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.369226 0.473841 5.00 5.73e-07 ***
Elevation -0.007038 0.002496 -2.82 0.00481 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.2593) family taken to be 1)

Null deviance: 41.974 on 42 degrees of freedom
Residual deviance: 36.343 on 41 degrees of freedom
AIC: 182.51

Number of Fisher Scoring iterations: 1


Theta: 0.2593
Std. Err.: 0.0755

2 x log-likelihood: -176.5090
Représenter le modèle final
Étape 1
Représenter les données et utiliser les estimations des paramètres pour représenter le modèle.
Utilisez
summary()
pour obtenir les paramètres
Étape 2
Utilisez les écarts-types pour construire l'intervalle de confiance.
C'est une combinaison de deux distributions (
Poisson
et
Gamma
)
Les
Yi
suivent une distribution de Poisson dont la moyenne, μ suit une distribution Gamma!
La distribution a
deux paramètres

µ
et
k
.
k
contrôle pour la dispersion (plus la dispersion est forte, plus
k
est petit)
GLM avec des données d'abondance
GLM avec variable binaire
Variables binaires
Présence/Absence d'une espèce
Présence/Absence d'une maladie
Succès/Échec d'observer un comportement
Survie/Mort d'un organisme
*** Présence/Absence ~ Environnement
*** Régression logistique ou modèle logit
Variables binaires
Dans R, on code une variable binaire avec 1 et 0:
Site Présence
A 1
B 0
C 1
D 1
E 0
F 1
L'espèce est présente !
L'espèce est absente !
GLMMs
Pourquoi être normale?
Régression logistique
La fonction
glm()
logit.reg <- glm(formula, data, family)
Variables binaires
Pas distribué normalement !
Variables binaires
Les valeurs prédites peuvent se trouver hors de l'intervalle [0,1] avec lm() :
Permet de spécifier une distribution de probabilité
ET
une fonction de lien
Distribution de probabilité
La distribution de Bernoulli :
E(
Y
) =
p
Var(
Y
) =
p
* (1-
p
)
La fonction de lien
Pour un modèle linéaire de base :
μ = Xβ
où μ représente les valeurs prédites par le modèle
X est la matrice du modèle (
i.e.
les variables explicatives)
β est le vecteur des paramètres estimés
(
i.e.
intercept et pente)
Xβ se nomme le
prédicteur linéaire
La fonction de lien
μ = Xβ
est vrai seulement pour la distribution normale
Si ce n'est pas le cas, on utilise
une transformation sur μ
g
(μ) = Xβ
La fonction de lien
g
(μ) = log
Pour les données binaires, on utilise la transformation
logit
.
1-μ
μ
1) La cote (μ / 1-μ)
2) Transformation log
La fonction de lien
g
(μ) = log
1-μ
μ
1) La cote (μ / 1-μ)
2) Transformation log
La cote met les valeurs prédites
sur une échelle de 0 à +Inf
La transformation log met les valeurs prédites sur une échelle de -Inf à +Inf
Les valeurs prédites transformées sont reliées
linéairement
au prédicteur linéaire
μ = variables prédites
Exercice 1
setwd(~/workshop6)
mites <- read.csv("mites.csv")
str(mites)
str(mites)
'data.frame': 70 obs. of 9 variables:
$ Galumna : int 8 3 1 1 2 1 1 1 2 5 ...
$ pa : int 1 1 1 1 1 1 1 1 1 1 ...
$ totalabund: int 140 268 186 286 199 209 162 126 123 166 ...
$ prop : num 0.05714 0.01119 0.00538 0.0035 0.01005 ...
$ SubsDens : num 39.2 55 46.1 48.2 23.6 ...
$ WatrCont : num 350 435 372 360 204 ...
$ Substrate : Factor w/ 7 levels "Barepeat","Interface",..: 4 3 2 4 4 4 4 2 3 4 ...
$ Shrub : Factor w/ 3 levels "Few","Many","None": 1 1 1 1 1 1 1 2 2 2 ...
$ Topo : Factor w/ 2 levels "Blanket","Hummock": 2 2 2 2 2 2 2 1 1 2 ...
Exercice 1
Spécifiez un modèle de la présence de
Galumna
sp. en fonction du contenu en eau du sol et de la topographie :
logit.reg <- glm(pa ~ WatrCont + Topo, data=mites,
family = binomial(link = "logit"))

summary(logit.reg)
Call:
glm(formula = pa ~ WatrCont + Topo, family = binomial, data = mites)

Deviance Residuals:
Min 1Q Median 3Q Max
-2.0387 -0.5589 -0.1594 0.4112 2.0252

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.464402 1.670622 2.672 0.007533 **
WatrCont -0.015813 0.004535 -3.487 0.000489 ***
TopoHummock 2.090757 0.735348 2.843 0.004466 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 91.246 on 69 degrees of freedom
Residual deviance: 48.762 on 67 degrees of freedom
AIC: 54.762

Number of Fisher Scoring iterations: 6
summary(logit.reg)
Défi 5
En utilisant le jeu de données 'bacteria', spécifiez un modèle de la présence de
H. influenzae
en fonction du traitement et de la semaine de test. Commencez avec un modèle saturé et trouvez le modèle le plus parcimonieux.
install.packages("MASS")
library(MASS)
data(bacteria)
str(bacteria)
?bacteria
Solution
model.bact1 <- glm(y ~ trt * week,
data = bacteria, family = binomial)
model.bact2 <- glm(y ~ trt + week,
data = bacteria, family = binomial)
anova(model.bact1, model.bact2, model.bact3,
test = "LRT")
model.bact3 <- glm(y ~ week, data = bacteria,
family = binomial)
Analysis of Deviance Table

Model 1: y ~ trt * week
Model 2: y ~ trt + week
Model 3: y ~ trt
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 214 203.12
2 216 203.81 -2 -0.6854 0.709836
3 217 210.72 -1 -6.9140 0.008552 **
En se basant sur ces résultats, on doit choisir le modèle #3 comme celui représentant le mieux les données.
Interpréter la sortie dans R
Regardez les coefficients du modèle 'logit.reg' :
summary(logit.reg)$coefficients
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.46440 1.67062 2.67230 0.00753
WatrCont -0.01581 0.00454 -3.48673 0.00049
TopoHummock 2.09076 0.73535 2.84322 0.00447
Interpréter la sortie dans R
On a utilisé une transformation logit !
Pour interpréter les coefficients du modèle, on doit utiliser une fonction "inverse" :
La fonction exponentielle pour obtenir les cotes : e
x
La fonction logit inverse pour obtenir les probabilités :
logit = 1 / (1 + 1/e )
-1
x
Interpréter la sortie dans R
Pour le contenu en eau sur l'échelle des cotes :
exp(logit.reg$coefficient[2])
WatrCont
0.9843118
Pour le contenu en eau sur l'échelle
des probabilités :
1 / (1 + 1/exp(logit.reg$coefficient[2]))
WatrCont
0.4960469
Pouvoir prédictif et ajustement du modèle
Le pseudo-R, un concept analogue du R pour les modèles estimés par maximisation de la vraisemblance :
Déviance nulle - déviance résid.
Déviance nulle
pseudo-R =
2
2
2
Variance expliquée par le modèle
Exercice 2
Calculez le pseudo-R du modèle 'logit.reg' :
2
pseudoR2 <- (logit.reg$null.deviance –
logit.reg$deviance) / logit.reg$null.deviance
pseudoR2
[1] 0.91
Pouvoir prédictif et ajustement du modèle
Le test de Hosmer-Lemeshow :
Similaire à un test de khi carré
Dans le paquet 'binomTools'
Faites le test de Hosmer-Lemeshow :
HLtest(Rsq(logit.reg))
Une valeur non significative indique un modèle adéquat !
Défi 6
En utilisant le modèle créé avec le jeu de données 'bacteria', évaluez le pouvoir prédictif et l'ajustement de ce modèle.
Comment faire pour améliorer le pouvoir explicatif du modèle ?
1)
2)
Solution
null.d <- model.bact2$null.deviance
resid.d <- model.bact2$deviance
bact.pseudoR2 <- (null.d -resid.d) / null.d
1)
2)
Ajouter des variables explicatives
pertinentes
pourrait certainement augmenter le pouvoir explicatif du modèle.
HLtest(Rsq(model.bact2)
Modélisation des données d'abondance
faramea <- read.csv(‘faramea.csv’, header = TRUE)
str(faramea)
Nous allons maintenant utiliser un nouveau jeux de données pour illustrer les GLMs avec données d'abondance.
Le nombre d'arbres de l'espèce
Faramea occidentalis
a été compté dans 43 quadrats sur l'île de Barro Colorada (Panama). Des données environnementales, comme l'élévation et la précipitation ont aussi été mesurées mesurées .
Examinons maintenant à quoi ressemble la distribution du nombre d'arbres par transect.
Étant donné cette spécificité propore aux données de dénombrement, la
distribution Poisson
semble le meilleur choix pour modéliser ces données...
Commençons par ajuster ce type de modèle pour modéliser l'abondance de
Faramea occidentalis
en fonction de l'élévation:
Coefficient de régression de l'élévation
Ordonnée à l'origine
Solution
# GLM Poisson
glm.p <- glm(Galumna~WatrCont+SubsDens, data=mites, family = poisson)

# GLM quasi-Poisson
glm.qp <- update(glm.p,family = quasipoisson)

# sélection de modèle
drop1(glm.qp, test = "Chi")
# ou
glm.qp2 <- glm(Galumna~WatrCont, data=mites, family = quasipoisson)
anova(glm.qp2, glm.qp, test="Chisq")
Déviance résiduelle
Degrés de liberté résiduels
glm.quasipoisson <- update(glm.poisson,family=quasipoisson)
Créez un nouveau GLM à l'aide de la famille '
quasipoisson
' ou actualisez le modèle précédent :
= 4!
0.0006436*4
null.model <- glm(Faramea.occidentalis~1, data=faramea, family=quasipoisson)

anova(null.model, glm.quasipoisson, test="Chisq")
Testons l'effet de l'élévation par une analyse de déviance :
Model 1: Faramea.occidentalis ~ 1
Model 2: Faramea.occidentalis ~ Elevation
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 42 414.81
2 41 388.12 1 26.686 0.1961
summary(glm.negbin)
k
summary(glm.negbin)$coefficients[1,1]
summary(glm.negbin)$coefficients[2,1]
summary(glm.negbin)$coefficients[1,2]
summary(glm.negbin)$coefficients[2,2]
Représenter le modèle final
L'effet de l'élévation sur le nombre de
F. occidentalis
est significatif, mais cette relation est faible.
Challenge
Utilisez le jeux de données
mites
!
Modélisez l'abondance de l'espèce
Galumna
en fonction des caractéristiques du substrat (son contenu en eau [
WatrCont
] et sa densité [
SubsDens
])
mites <- read.csv("mites.csv", header = TRUE)
N.S.
Limite sup
Limite inf
Faut-il contrôler pour la surdispersion ?
Quelles variables explicatives ont un effet significatif ?
Sélectionnez le meilleur modèle !

Sélection pas à pas en retirant à chaque fois une variable et en comparant le modèle emboité au modèle complet :
Conseils
Sélection de modèle avec des GLMs

Spécifiez un modèle emboîté manuellement, appelez le MyGLM2, et utilisez la fonction
anova()
:
drop1(MyGLM, test = "Chi")
anova(MyGLM, MyGLM2, test = "Chi")
setwd("~/Desktop")
mites <- read.csv('mites.csv')
head(mites)
str(mites)
Limites des modèles linéaires (mixtes)
Charger les données et appliquer un ML
Le jeu de données chargé contient une partie du jeu de données classique des mites oribatidées:





Objectif:
Modéliser l'abondance ($abund), la présence ($pa), et la proportion de Galumna ($prop) en fonction des 5 paramètres environnementaux.
Explorer les relations
Pouvons-nous détecter des relations entre
Galumna
et les variables environnementales?
plot(mites)
Une relation négative entre Galumna et le contenu d'eau du sol?
par(mfrow=c(1,3))

plot(Galumna ~ WatrCont, data = mites, xlab = 'Water content', ylab='Abundance')

boxplot(WatrCont ~ pa, data = mites, xlab='Presence/Absence', ylab = 'Water content')

plot(prop ~ WatrCont, data = mites, xlab = 'Water content', ylab='Proportion')

par(mfrow=c(1,1))


Utiliser des modèles linéaires pour voir si l'abondance, la présence/absence, et la proportion de
Galumna
varie en fonction du contenu d'eau.
lm.abund <- lm(Galumna ~ WatrCont, data = mites)
summary(lm.abund)
lm.pa <- lm(pa ~ WatrCont, data = mites)
summary(lm.pa)
lm.prop <- lm(prop ~ WatrCont, data = mites)
summary(lm.prop)
Un effet significatif dans tous les modèles! Attends une seconde...
plot(Galumna ~ WatrCont, data = mites)
abline(lm.abund)
plot(lm.abund)
Encore pire pour les autres modèles:
#Proportion
plot(prop ~ WatrCont, data = mites)
abline(lm.prop)
plot(lm.prop)

#Présence/Absence
plot(pa ~ WatrCont, data = mites)
abline(lm.pa)
plot(lm.pa)
Très commun en écologie que les suppositions des modèles linéaires ne soient pas respectées.
Raison principale pourquoi nous avons souvent besoin des modèles linéaires généralisé (GLMs)



Revisitons les hypothèses du modèle linéaire...
Équation du modèle linéaire :

y =
β
+
β
x +
ε


où:
yi = valeur estimée pour la variable réponse
β0 = l'ordonnée à l'origine de la droite
β1 = pente
xi = valeur de la variable indépendante
εi = résidus du modèle obtenus d'une distribution normale dont la moyenne varie mais dont la variance est constante
**
**Point clé:
les résidus (distances entre les observations et la droite de régression) peuvent être prédits en échantillonnant de façon aléatoire une distribution normales
Rappel:
la distribution normale a deux paramètres: μ (la moyenne) et σ (la variance):
i
i
i
0
1
2
Une autre équation de représenter un modèle linéaire est:

yi ~ N(μ = β0 + β1xi, σ ),

Ce qui signifie que yi est échantillonner d'une distribution normale ayant les paramètres μ (dont la valeur dépend de xi) et σ (qui garde la même valeur pour toutes les valeurs de Y)



Essayons de prédire l'abondance de
Galumna
en fonction du contenu d'eau en utilisant le modèle linéaire que nous avons vus plus tôt...
2
2
Nous avons besoin de coefficients de régression (betas) et de sigma:


coef(lm.abund)
summary(lm.abund)$sigma
Quels sont les paramètres de la distribution normale utilisée pour modéliser y lorsque le contenu d'eau = 300?
yi ~ N(μ = β0 + β1xi, σ )
2
μ = 3.44 + (-0.006 x 300) = 1.63
σ = 1.51
2
Lorsque y = 300, les résidus devraient suivre une distribution normale avec μ = 1.63 et σ = 1.51.
Lorsque y = 400, les résidus devraient suivre une distribution normale avec μ = 1.02 and σ = 1.51, etc.

Graphiquement, notre modèle ressemble à:
Distribution des données biologiques
Problèmes
:
σ n'est pas homogène mais la fonction lm nous contraint d'utiliser toujours la même valeur de σ
les valeurs estimées devraient être des nombres entiers
Les statisticiens ont développé une foule de lois de probabilité (distributions) correspondant à divers types de données

Une distribution donne la probabilité d'observer chaque issue possible d'une expérience ou échantillonage (e.g. abund = 8
Galumna
)

Les distributions peuvent être
discrètes
(que des nombres entiers) ou
continues
(incluent aussi des fractions)

Toutes les distributions ont des
paramètres
qui déterminent la forme (e.g. μ et σ pour la distribution normale)
L'abondance de
Galumna
suit une distribution discrète (que des nombres entiers)

Une loi utile pour modéliser les données d'abondance est la loi de Poisson:
une distribution discrète avec un seul paramètre, λ (lambda), qui détermine la moyenne et la variance de la distribution:
Galumna
semble suivre une loi de Poisson avec une faible valeur de λ:
hist(mites$Galumna)
mean(mites$Galumna)
Présence-absence par contre prend une autre forme:
Inclut que des 0 et des 1, de telle sorte que la loi de Poisson n'est appropriée pour cette variable
hist(mites$pa)
Distribution “Bernoulli”
n'inclut que deux issues possibles dans son ensemble: succès (1) ou échec (0)
N'a qu'un paramètre: p, la probabilité de succès
Distribution binomiale:
Lorsqu'il y a plusieurs épreuves (chacune avec un succès/échec), la loi de Bernoulli se transforme en loi binomiale
Inclue le paramètre additionel n, le nombre d'épeuves
Prédit la probabilité d'observer une certaine proportion de succès, p, sur le nombre total d'épreuves, n
Distribution binomiale:
utilisée pour modéliser des données lorsque le nombre de succès est donné par un nombre entier, et lorsque le nombre d'épreuves, n, est connu

Différence principale avec la loi de Poisson:
L'étendue de la loi binomiale a une limite supérieure, n.
Par conséquent, est asymétrique et décalée à gauche lorsque p est faible, mais décalée à droite lorsque p est élevé.
vs.
Retournons à notre problème... pouvons changer la distribution des résidus (εi) de normale à Poisson:

yi ~ Poisson(λ = β0 + β1xi)

Le problème est résolu!
1. λ varie selon x (contenu d'eau)
donc la variance des résidus variera aussi selon x, et nous venons de nous défaire de la supposition d'homogénéité de la variance!

2. Les valeurs estimées seront des nombres entiers plutôt que des nombres décimaux

3. Ce modèle ne prédira jamais de valeurs négatives (Poisson est toujours strictement positif)
Ce modèle est
presqu'un
MLG de Poisson, qui ressemble à ça:
Les probabilités (en orange) sont maintenant des nombres entiers, et la variance que la moyenne de la distribution déclinent lorsque λ diminue avec le contenu d'eau.
Structure dans le jeu de données ou corrélation entre les observations peut entraîner une
dépendance entre les observations
échantillonnés à partir des mêmes sites ou points dans le temps

Effets aléatoires:
Un échantillon de la population; les sujets que vous avez échantillonné par hasard
Explique la variation de la variable réponse

Effets fixe:
Reproductible; serait le même dans toutes études
Explique la moyenne de la variable réponse

Révision: modèles linéaires mixtes
Sont une extension des GLMs qui tiennent compte de cette structure supplémentaire des données.

Nous pouvons préciser ces modèles en suivant des étapes similaires à celles introduites dans cet atelier et l'atelier sur les LMMs.
Les modèles linéaires généralisés mixtes (GLMM)
Estimations de retrait
Lancer un GLMM sous R
dat.tf <- read.csv("Banta_TotalFruits.csv")
str(dat.tf)
Chargez les données
Arabidopsis
'Banta_TotalFruits.csv' :
L'effet de la disponibilité de nutriments et d'herbivores (
effets fixes
)

sur la production de fruits d'
Arabidopsis thaliana
(arabette des dames) a été évaluée en mesurant 625 plantes à travers neuf populations différentes, constituées chacune de 2 à 3 génotypes (
effets aléatoires
).
# popu facteur avec un niveau pour chaque population
# gen facteur avec un niveau pour chaque génotype
# nutrient facteur avec niveaux bas ( = 1) ou élevé (value = 8)
# amd facteur précisant l'absence ou la présence d'herbivorie
# total.fruits nombre entier indiquant le nombre de fruits par plante
dat.tf <- transform(dat.tf,
X=factor(X),
gen=factor(gen),
rack=factor(rack),
amd=factor(amd,levels=c("unclipped","clipped")),
nutrient=factor(nutrient,label=c("Low","High")))

Exploration des données
Avant de commencer, nous allons nettoyer les données (rendre nombres entiers en facteurs, ajuster niveaux d'herbivorie (amd) et renommer les niveaux de nutriments) :
require(lme4)
require(coefplot2)
require(reshape)
require(ggplot2)
require(plyr)
require(gridExtra)
require(emdbook)
source("glmm_funs.R")
Exploration de la variance
dat.tf <- within(dat.tf,
{

# genotype x nutrient x herbivory
gna <- interaction(gen,nutrient,amd)
gna <- reorder(gna, total.fruits, mean)

# population x nutrient x herbivory
pna <- interaction(popu,nutrient,amd)
pna <- reorder(pna, total.fruits, mean)
})
Pour illustrer l'hétérogénéité de la variance, nous allons d'abord créer des boîtes à moustaches (boxplots) de la variable réponse par rapport aux différents facteurs environnementaux
Créons de nouvelles variables qui représentent toutes les combinaisons de
nutriments
x
herbivorie
x
facteur aléatoire
# Boîte à moustaches total fruits vs nouvelle variable (génotype x nutriments x herbivorie)
ggplot(data = dat.tf, aes(factor(x = gna),y = log(total.fruits + 1))) +
geom_boxplot(colour = "skyblue2", outlier.shape = 21,
outlier.colour = "skyblue2") +
theme_bw() + theme(axis.text.x=element_text(angle=90)) +
stat_summary(fun.y=mean, geom="point", colour = "red")
* De même, la boîte à moustaches total des fruits vs population x nutriments x herbivorie montre une grande quantité d'hétérogénéité entre les populations.
Distribution de la variable réponse
La variable réponse constitue des données d'abondance, donc nous devons choisir une
distribution de Poisson
(variance est égale à la moyenne)
Cependant, comme nous le verrons, la variance de chaque groupe augmente beaucoup plus rapidement que prévu...
Choisir une distribution
Il existe donc une importante hétérogénéité parmi la variance de chaque groupe, même lorsque la variable réponse est transformé

Si nous représentons graphiquement les
écarts vs moyennes par groupes
(génotype x nutriment x herbivorie), on voit que la distribution de Poisson est la moins appropriée (écart augmentent beaucoup plus vite que la moyenne)
Poisson GLMM
Compte tenu de la relation moyenne-variance, nous avons besoin d'un modèle avec surdispersion. Commençons avec un modèle de Poisson.

Pour lancer un GLMM dans R, nous faisons appel à la fonction
glmer()
:
mp1 <- glmer(total.fruits ~ nutrient*amd + rack + status +
(1|popu)+
(1|gen),
data=dat.tf, family="poisson")
Spécifie la famille
Spécifie les effets aléatoires
Nous pouvons vérifier la surdispersion en utilisant la fonction
overdisp_fun
qui divise la déviance des résidus (résidus de Pearson) par les degrés de liberté des résidus et test si le rapport est plus grand que 1
Vérification de la surdispersion
# Surdispersion?
overdisp_fun(mp1)
Les effets aléatoires sont souvent appelés des
estimations de rétrécissement
parce qu'ils représentent une moyenne pondérée des données et de l'ajustement global (effet fixe)

Le retrait vers l'ajustement global est plus sévère si la variabilité intra-groupe est grande par rapport à la variabilité inter-groupes
Révision: LMMs
Exploration de la variance
factor(genotype x nutrient x clipping)
> overdisp_fun(mp1)
chisq ratio p logp
15755.86833 25.57771 0.00000 -6578.47027
Ratio est » 1
Comme on s'y attendait, nous devons modéliser une distribution différente où la variance augmente plus rapidement que la moyenne
GLMM binomiale négative
(Poisson-gamma)
La distribution binomiale négative satisfait la supposition que la
variance est proportionnelle au carré de la moyenne
mnb1 <- glmer.nb(total.fruits ~ nutrient*amd + rack + status +
(1|popu)+
(1|gen),
data=dat.tf, control=glmerControl(optimizer="bobyqa"))
# Control spécifie la façon dont nous optimisons les valeurs des paramètres
# Surdispersion?
overdisp_fun(mnb1)
Le rapport est maintenant beaucoup plus près de 1 mais la valeur p < 0,05
GLMM Poisson-lognormal
Une autre option
(que nous n’avons pas encore vue)
est la distribution
Poisson-lognormal

Une distribution Poisson-lognormal avec moyenne µ et variance lognormal σ^2 est donnée par:

var(y) = µ + µ^2 [exp(σ^2) - 1]

En revanche, nous avons vu que la distribution
binomiale négative
(Poisson-gamma) était donnée par:

var(y) = µ + µ^2/k

Pour modéliser un modèle Poisson-lognormal, nous devons modéliser les εi
a priori
avec une distribution lognormal, afin que σ^2 dépende du niveau de regroupement sélectionné (individues, génotype ou population)
Ici, nous allons évaluer un modèle avec la σ^2 regroupé au niveau individuel (1|X)
GLMM Poisson-lognormal
mpl1 <- glmer(total.fruits ~ nutrient*amd + rack + status +
(1|X) +
(1|popu)+
(1|gen),
data=dat.tf, family="poisson", control=glmerControl(optimizer="bobyqa"))

overdisp_fun(mpl1)
chisq ratio p logp
1.775361e+02 2.886766e-01 1.000000e+00 -3.755374e-73
Rapport maintenant conforme avec notre critère
Intercepts aléatoires
Évaluation des intercepts pour
population
et
genotype
avec
drop1 et l'approche LRT :

# popu seulement
mpl1.popu <- glmer(total.fruits ~ nutrient*amd + rack + status +
(1|X) +
(1|popu),
data=dat.tf, family="poisson", control=glmerControl(optimizer="bobyqa"))

# gen seulement
mpl1.gen <-glmer(total.fruits ~ nutrient*amd + rack + status +
(1|X) +
(1|gen),
data=dat.tf, family="poisson", control=glmerControl(optimizer="bobyqa"))

# Avec AICc
ICtab(mpl1,mpl1.popu, type = c("AICc"))
* Nous testons l'hypothèse nulle d'une variance = 0. Mais nous ne pouvons pas avoir un écart négatif. Donc, nous testons le paramètre à la limite de sa région réalisable et par conséquent, la valeur de p rapportée est environ le double de ce qu'elle devrait être.
Représentation graphique des paramètres du modèle
Une représentation graphique des paramètres du modèle peut être obtenu en utilisant la fonction
coefplot2()

:
# Paramètres de la variance
coefplot2(mpl1,ptype="vcov",intercept=TRUE)

# Effets fixes
coefplot2(mpl1,intercept=TRUE)
Graphiques des intercepts aléatoires
Nous pouvons aussi extraire les effets aléatoires en utilisant la fonction
ranef()
et les tracer en utilisant la fonction
dotplot()
Il y a une variabilité régionale parmi les populations:
Populations espagnoles (SP) ont des valeurs plus élevées que les populations suédoises (SW) et néerlandaises (NL)
Différence entre les génotypes semble largement induite par génotype 34
*La fonction
grid.arrange()

a été utilisé pour faire ce graphique (voir scripte).
Modèle final
Nous pouvons aussi utiliser les fonctions
drop1()
et
dfun()
pour évaluer nos effets fixes (dfun convertit les valeurs AIC retournées par drop1 en valeurs ΔAIC)
(dd_LRT <- drop1(mpl1,test="Chisq"))
(dd_AIC <- dfun(drop1(mpl1)))
Df dAIC
<none> 0.000
rack 1 55.083
status 2 1.612
nutrient:amd 1 1.444
Fort effet de
rack
(dAIC = 55.08 si on enlève cette variable)
Effets de
status
et de l'
interaction
sont faibles (dAIC < 2)
Commençons par
enlever l'interaction non significative
afin de tester les effets principaux de nutriments et d’herviborie
Modèle final
anova(mpl1,mpl2)

(dd_LRT2 <- drop1(mpl2,test="Chisq"))
(dd_AIC2 <- dfun(drop1(mpl2)))
Df dAIC
<none> 0.000
nutrient 1 135.688
amd 1 10.218
rack 1 54.231
status 2 1.286
Fort effets de nutriments et d’herbivorie
(grand changement d'AIC si supprimé)
Notre modèle final inclut l’effet fixe de nutriments, d’herbivorie, la variable nuisance de rack, l’effet aléatoire au niveau de l'observation (1 | X) et la variation de fruits par populations et génotypes
DÉFI 7
En utilisant l’ensemble de données
inverts
(temps de développement larvaire (PLD) de 74 espèces d'invertébrés et vertébrés marins élevés à différentes températures et temps), répondez aux questions suivantes:

Quel est l'effet du type d'alimentation et du climat (
effets fixes
) sur PLD?
Est-ce que cette relation varie selon les taxons (
effets aléatoires
)?
Quelle est la
meilleure famille de distribution
pour ces données?
GLM et données de proportions
Les données de proportions peuvent être analysées par régression logistique.
Si on mesure un nombre d'occurrences et qu'on connaît la taille d'échantillon, on obtient des données de proportions !
Supposons qu'on mesure la prévalence d'une maladie sur dix cerfs dans 10 populations différentes:
x cerfs infectés
10 cerfs
}
toujours entre
0 et 1 !
Exercice 3
Dans R, on doit spécifier le nombre de fois qu'un événement s'est produit et le nombre de fois qu'un événement ne s'est pas produit :
> prop.reg <- glm(cbind(Galumna, totalabund - Galumna) ~ Topo + WatrCont, data = mites, family = binomial)

> summary(prop.reg)
Call:
glm(formula = cbind(Galumna, totalabund - Galumna) ~ WatrCont +
Topo, family = binomial, data = mites)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.4808 -0.9699 -0.6327 -0.1798 4.1688

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.288925 0.422109 -7.792 6.61e-15 ***
WatrCont -0.005886 0.001086 -5.420 5.97e-08 ***
TopoHummock 0.578332 0.274928 2.104 0.0354 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 140.702 on 69 degrees of freedom
Residual deviance: 85.905 on 67 degrees of freedom
AIC: 158.66
Puis charger les librairies et fonctions nécessaires :
NB = binomiale negative


QP = quasi-Poisson



loess = Locally weighted
regression smoothing
Barres d'erreur ne sont visibles que pour les effets fixes parce glmer ne nous donne pas d'informations sur l'incertitude des effets aléatoires.
variance: effets aléatoires
coefficient: effets fixes
DÉFI 7

1) Pour comparer un glmer avec un glm, vous devez extraire et de comparer les probabilités de chaque modèle manuellement - vous ne pouvez pas utiliser
anova()
:
NL1 <- -logLik(glmer.mod)
# glmer.mod est le nom de votre modèle glmer
NL0 <- -logLik(glm.mod)
# glm.mod est le nom de votre modèle glm
devdiff <- 2*(NL0-NL1)
dfdiff <- attr(NL1,"df")-attr(NL0,"df")
pchisq(devdiff,dfdiff,lower.tail=FALSE)
CONSEILS
DÉFI 7

2) Pour exécuter la fonction
drop1()
sur un GLMM NB, vous devez spécifier thêta et exécuter glmer (..., la famille = negative.binomial) au lieu de glmer.nb (...)
# Obtenir thêta du modèle NB:
theta.mnb1 <- theta.md(inverts$PLD, fitted(mnb1), dfr = df.residual(mnb1))
# Utilisez glmer au lieu de glmer.nb
mnb1 <- glmer(..., family=negative.binomial(theta=theta.mnb1))
# Maintenant drop1 et dfun fonctionnent
(dd_LRT <- drop1(mnb1,test="Chisq")
(dd_AIC <- dfun(drop1(mnb1)))
CONSEILS
library(binoTools)
On peut coder le modèle directement avec les proportions :
> prop.reg2 <- glm(prop ~ Topo + WatrCont, data = mites, family = binomial, weights = totalabund)

Intercepts aléatoires
Maintenant que nous avons la distribution appropriée, nous pouvons tester l'importance des intercepts aléatoires (pour
population
et
génotype
)

Les modèles emboîtés peuvent être comparés en utilisant soit
l'approche AICc (comme dans l'atelier 5), ou
l'approche fréquentiste (drop1), où l'importance de chaque terme est évaluée en utilisant
anova()

et le test du rapport de vraisemblance (LRT)*
Intercepts aléatoires
Évaluation des intercepts pour
population
et
génotype
avec
AICc
:
p <0,05 pour lest deux tests de vraisemblance (LRT), donc le modèle avec les deux effets aléatoires (mpl1) est sélectionné.
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
mpl1.gen 9 5031.5 5071.5 -2506.8 5013.5
mpl1 10 5015.4 5059.8 -2497.7 4995.4 18.177 1 2.014e-05 ***
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
mpl1.popu 9 5017.4 5057.4 -2499.7 4999.4
mpl1 10 5015.4 5059.8 -2497.7 4995.4 4.0639 1 0.04381 *
Modèle final
Nous terminons cette section sur les GLMMs avec l'évaluation des effets fixes

Nous décrivons d'abord les modèles potentiels et les comparons en utilisant AICc
*
:
mpl2 <- update(mpl1, . ~ . - rack)
# modèle sans rack
mpl3 <- update(mpl1, . ~ . - status)
# modèle sans status
mpl4 <- update(mpl1, . ~ . - amd:nutrient)
# modèle sans interaction amd:nutrient
ICtab(mpl1,mpl2,mpl3,mpl4, type = c("AICc"))
dAICc df
mpl1 0.0 10
mpl4 1.4 9
mpl3 1.5 8
mpl2 55.0 9
Modèle final
mpl2 <- update(mpl1, . ~ . - amd:nutrient)


mpl3 <- update(mpl2, . ~ . - rack)
# modèle sans rack ni l’interaction
mpl4 <- update(mpl2, . ~ . - status)
# modèle sans status ni l’interaction
mpl5 <- update(mpl2, . ~ . - nutrient)
# modèle sans nutrient ni l’interaction
mpl6 <- update(mpl2, . ~ . - amd)
# modèle sans amd ni l’interaction
ICtab(mpl2,mpl3,mpl4,mpl5,mpl6, type = c("AICc"))
Ou avec
drop1
:
Supprimer l'interaction et évaluer les effets nutriments et herbivorie, rack et status en utilisant AICc
dAICc df
mpl2 0.0 9
mpl4 1.2 7
mpl6 10.2 8
mpl3 54.2 8
mpl5 135.6 8
anova(mpl1,mpl1.popu)
dAICc df
mpl1 0.0 10
mpl1.popu 2.0 9
mpl1.gen 16.1 9
anova(mpl1,mpl1.gen)
Meilleur modèle est mpl1!
*
Nous ne couvrons pas tous les modèles possibles ci-dessus, cependant, l'interaction amd: nutriments ne peut être évaluée que si amd et nutriments sont présents dans le modèle.
Explorer les relations
Pouvons-nous détecter des relations entre
Galumna
et les variables environnementales?
Explorer les relations
Tester la linéarité
Utiliser des modèles linéaires pour voir si l'abondance, la présence/absence, et la proportion de
Galumna
varie en fonction du contenu d'eau.
Tester la linéarité
Tester la linéarité
Tester la linéarité
Suppositions des modèles linéaires
Suppositions des modèles linéaires
Distribution normale des résidus
Distribution normale des résidus
Prédiction du modèle
Prédiction du modèle
Distribution des données biologiques
Distribution des données biologiques
Distribution des données biologiques
Distribution des données biologiques
Distribution des données biologiques
Distribution des données biologiques
Distribution des données biologiques
Distribution des données biologiques
= Xβ
Modélisation des données d'abondance
GLM avec une distribution de Poisson
Ajuster une binomiale négative sous R
data.frame': 43 obs. of 9 variables:
...
$ Precipitation : num 2530 2530 2993 3072 3007 ...
$ Elevation : int 120 120 20 100 180 180 40 30 60 50 ...
$ Age : int 3 3 2 3 1 1 2 2 1 3 ...
...
$ Faramea.occidentalis: int 14 7 0 0 2 1 0 0 2 0 ...
Importez les données 'faramea.csv' dans R.
hist(faramea$Faramea.occidentalis, breaks=seq(0,45,1), main="", col="grey")
Que des valeurs entières et positives
plot(faramea$Elevation, faramea$Faramea.occidentalis)
Y ~ Poisson(µ)

E(Y) = Var(Y) = µ
La transformation log permet aux valeurs d'abondance (0 à l'infini) d'être distribuées de moins l'infini à l'infini.
Res dev
est une statistique qui suit un distribution Chi-carré avec df = df_sat - df_complet = n - p. En moyenne, ce Chi-carré est égale a df = n - p.
Pour donner un bref aperçu des GLMM, nous allons examiner une étude des données d'abondance parvenant de différentes populations (possiblement un manque d'indépendance).
[1] 0.4655937
Calculez le ROC (recursive operating curve) du modèle 'logit.reg' :
library("ROCR")
ROC
# 70 échantillons de mousses & mites
# 5 variables environnementales, abondance de la mite ”Galumna sp.”, et abondance totale de mites.
(vos données sont correctes;mais le modèle n'est pas adéquat pour vos données)
Nous avons besoin de coefficients de régression (betas) et de sigma:


coef(lm.abund)
summary(lm.abund)$sigma
Quels sont les paramètres de la distribution normale utilisée pour modéliser y lorsque le contenu d'eau = 300?
Prédiction du modèle
yi ~ N(μ = β0 + β1xi, σ )
2
2
2
2
2
2
Nous pouvons utiliser la loi de Bernoulli pour calculer la probabilité d'obtenir l'issue “Galumna présent” (1) vs. “Galumna absent” (0)
Full transcript