Formation Mixed model
nicolas Saby
du 29/09
Prise de notes pendant la formation R modèles mixte avec Sophie ? Anastats.
1 Modèle linéaire classique et MLM
Un travail essentiel consiste à représenter ses données graphiquement. Voici des exemples
library(car)## Warning: package 'car' was built under R version 3.2.2## Package SparseM (0.97) loaded.
## To cite, see citation("SparseM")setwd(dir = "E:/nsaby/cours_fac/Modeles_mixtes")
longi <- read.table(file="Exemples/Engrais.txt",
header = TRUE,
dec = ".",
sep = "\t")
summary(longi)## Plant Traitement Temps Mesure Groupe
## P1 : 3 Trt1:60 Sem1:40 Min. :28.0 Trt1-Sem1:20
## P10 : 3 Trt2:60 Sem2:40 1st Qu.:46.8 Trt1-Sem2:20
## P11 : 3 Sem3:40 Median :53.0 Trt1-Sem3:20
## P12 : 3 Mean :52.6 Trt2-Sem1:20
## P13 : 3 3rd Qu.:60.0 Trt2-Sem2:20
## P14 : 3 Max. :77.0 Trt2-Sem3:20
## (Other):102with(longi,table(Traitement,Plant,Temps))## , , Temps = Sem1
##
## Plant
## Traitement P1 P10 P11 P12 P13 P14 P15 P16 P17 P18 P19 P2 P20 P21 P22 P23
## Trt1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0
## Trt2 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1
## Plant
## Traitement P24 P25 P26 P27 P28 P29 P3 P30 P31 P32 P33 P34 P35 P36 P37 P38
## Trt1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## Trt2 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1
## Plant
## Traitement P39 P4 P40 P5 P6 P7 P8 P9
## Trt1 0 1 0 1 1 1 1 1
## Trt2 1 0 1 0 0 0 0 0
##
## , , Temps = Sem2
##
## Plant
## Traitement P1 P10 P11 P12 P13 P14 P15 P16 P17 P18 P19 P2 P20 P21 P22 P23
## Trt1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0
## Trt2 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1
## Plant
## Traitement P24 P25 P26 P27 P28 P29 P3 P30 P31 P32 P33 P34 P35 P36 P37 P38
## Trt1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## Trt2 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1
## Plant
## Traitement P39 P4 P40 P5 P6 P7 P8 P9
## Trt1 0 1 0 1 1 1 1 1
## Trt2 1 0 1 0 0 0 0 0
##
## , , Temps = Sem3
##
## Plant
## Traitement P1 P10 P11 P12 P13 P14 P15 P16 P17 P18 P19 P2 P20 P21 P22 P23
## Trt1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0
## Trt2 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1
## Plant
## Traitement P24 P25 P26 P27 P28 P29 P3 P30 P31 P32 P33 P34 P35 P36 P37 P38
## Trt1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## Trt2 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1
## Plant
## Traitement P39 P4 P40 P5 P6 P7 P8 P9
## Trt1 0 1 0 1 1 1 1 1
## Trt2 1 0 1 0 0 0 0 0Un graphique boxplot par traitement
## Warning: package 'ggplot2' was built under R version 3.2.2
On observe une forte hétérogénéité entre les individus. Elle est plus forte que l’effet du facteur.
Commencons par construire le modèle a effet fixe uniquement Vérifions la variabilité par groupe
on prend groupe car on s’interesse au résidus des 2 facteurs TRaitement et Temps
utilise l’option qui calcule à partir de la médiane
Intérêt de ce test est que on a pas besoin d’avoir des distributions gaussiennes
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 0.39 0.86
## 114## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 5 0.39 0.86
## 114On peut examiner la normalité des distributions par groupe avec la fonction qqmath
qqmath(~Mesure|Groupe, type = c("p","r") , data = longi)
Construisons le modèle a effets fixes. calculons l’analyse de variance Vérifions également la normalité des résidus
LE dernier graphique présente la variabilité issu de l’individu , les résidus par individus ne sont pas du tout centré. on confirme l’effet aléatoire de l’infividu
modfix <- lm( Mesure ~ Temps*Traitement, data = longi)
Anova(modfix)## Anova Table (Type II tests)
##
## Response: Mesure
## Sum Sq Df F value Pr(>F)
## Temps 109 2 0.59 0.56
## Traitement 21 1 0.22 0.64
## Temps:Traitement 42 2 0.23 0.80
## Residuals 10594 114shapiro.test(modfix$residuals)##
## Shapiro-Wilk normality test
##
## data: modfix$residuals
## W = 0.99, p-value = 0.4hist(modfix$residuals)
hist(rstandard(modfix))
par(mfrow=c(2,2))
plot(modfix)
par(mfrow=c(1,1))
plot(rstandard(modfix)~longi$Plant)
abline(0,0,lty=2)
Attention, discussion sur la déclaration des effets aléatoire, il y a toujours le 1, les effets aléatoires sont toujours après le | Dans la synthaxe ~1+Temps|Plant, Temps est fixe mais fait varier la pente et l’intercept par plante.
Dans ce contexte on prend la partie aléatoire en compte mais on ne la discute pas. On regarde que la partie fixe. l’interêt est de retirer la variabilité inhérente aux individus
library(nlme)
modmix1 <- lme(Mesure~Traitement*Temps,
random=~1|Plant, data = longi, method="REML")
AIC(modmix1)## [1] 731.6modmix2 <- update(modmix1, random = ~1+Temps|Plant)
modmix3 <- update(modmix1, random = ~1+Temps|Plant)
anova(modmix1,modmix2)## Model df AIC BIC logLik Test L.Ratio p-value
## modmix1 1 8 731.6 753.5 -357.8
## modmix2 2 13 725.5 761.1 -349.8 1 vs 2 16.07 0.0066Anova(modmix2)## Analysis of Deviance Table (Type II tests)
##
## Response: Mesure
## Chisq Df Pr(>Chisq)
## Traitement 1.24 1 0.27
## Temps 23.49 2 7.9e-06 ***
## Traitement:Temps 3.79 2 0.15
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(modmix2)## Linear mixed-effects model fit by REML
## Data: longi
## AIC BIC logLik
## 725.5 761.1 -349.8
##
## Random effects:
## Formula: ~1 + Temps | Plant
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 9.535 (Intr) TmpsS2
## TempsSem2 2.120 0.274
## TempsSem3 2.842 -0.514 -0.436
## Residual 1.850
##
## Fixed effects: Mesure ~ Traitement * Temps
## Value Std.Error DF t-value p-value
## (Intercept) 51.25 2.1718 76 23.598 0.0000
## TraitementTrt2 0.20 3.0714 38 0.065 0.9484
## TempsSem2 2.50 0.7529 76 3.321 0.0014
## TempsSem3 0.35 0.8638 76 0.405 0.6865
## TraitementTrt2:TempsSem2 -0.40 1.0647 76 -0.376 0.7082
## TraitementTrt2:TempsSem3 2.30 1.2216 76 1.883 0.0636
## Correlation:
## (Intr) TrtmT2 TmpsS2 TmpsS3 TT2:TS2
## TraitementTrt2 -0.707
## TempsSem2 0.065 -0.046
## TempsSem3 -0.463 0.327 0.061
## TraitementTrt2:TempsSem2 -0.046 0.065 -0.707 -0.043
## TraitementTrt2:TempsSem3 0.327 -0.463 -0.043 -0.707 0.061
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.629385 -0.341479 -0.007788 0.363696 1.686953
##
## Number of Observations: 120
## Number of Groups: 40plot(residuals(modmix2)~longi$Plant)
longi$predicted <- predict(modmix2)
ggplot(longi, aes(Mesure,predicted)) + geom_point() + facet_wrap(~Traitement) 
2 Modèle linéaire généralisé
On reprend ici le travail avec des données 0/1. on teste avec et sans effet aléatoire Labo. Par contre, on ne peut pas atteindre l’interêt de l’insertion de l’effet aléatoire. Il faut regarder les résidus. regarder la variance estimée. mais
Aubergines <- read.table(file="Exemples/Aubergines.txt",
header = TRUE,
dec = ".",
sep = "\t")
summary(longi)## Plant Traitement Temps Mesure Groupe
## P1 : 3 Trt1:60 Sem1:40 Min. :28.0 Trt1-Sem1:20
## P10 : 3 Trt2:60 Sem2:40 1st Qu.:46.8 Trt1-Sem2:20
## P11 : 3 Sem3:40 Median :53.0 Trt1-Sem3:20
## P12 : 3 Mean :52.6 Trt2-Sem1:20
## P13 : 3 3rd Qu.:60.0 Trt2-Sem2:20
## P14 : 3 Max. :77.0 Trt2-Sem3:20
## (Other):102
## predicted
## Min. :28.8
## 1st Qu.:47.2
## Median :51.9
## Mean :52.6
## 3rd Qu.:59.7
## Max. :74.8
## Table<-xtabs(~Developpement+Taille, data=Aubergines) # Réalise la table de contingence du nombre de cas en fonction de la taille
plot(Table, main="Nombre de cas en fonction de la taille et le ")# Réalise un graphique de la table de contingence créée ci-dessus
modglm <- glm(Developpement ~ Taille*Temperature, family=binomial(logit), data=Aubergines) # Réalise le modèle avec effets fixes.
modglm2 <- glm(Developpement ~ Taille, family=binomial(logit), data=Aubergines) # Réalise le modèle avec effets fixes.
Anova(modglm) # Affiche les résultats## Analysis of Deviance Table (Type II tests)
##
## Response: Developpement
## LR Chisq Df Pr(>Chisq)
## Taille 13.25 2 0.0013 **
## Temperature 1.18 2 0.5540
## Taille:Temperature 3.41 4 0.4916
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(modglm,type="III") # avec le type III, le calcul met plus de poids sur l'interaction ## Analysis of Deviance Table (Type III tests)
##
## Response: Developpement
## LR Chisq Df Pr(>Chisq)
## Taille 3.56 2 0.17
## Temperature 1.09 2 0.58
## Taille:Temperature 3.41 4 0.49AIC(modglm) # Donne la valeur du critère d'information d'Akaike du modèle.## [1] 244Table<-xtabs(~Developpement+Labo+Taille, data=Aubergines) # Réalise la table de contingence du nombre de cas en fonction de la taille par laboratoire.
plot(Table, main="Nombre de cas en fonction de la taille par labo")# Réalise un graphique de la table de contingence créée ci-dessus
library(lmerTest) # Charge la bibliothèque.## Warning: package 'lmerTest' was built under R version 3.2.2## Loading required package: Matrix
## Loading required package: lme4
## Loading required package: Rcpp## Warning: package 'Rcpp' was built under R version 3.2.2##
## Attaching package: 'lme4'
##
## The following object is masked from 'package:nlme':
##
## lmList
##
##
## Attaching package: 'lmerTest'
##
## The following object is masked from 'package:lme4':
##
## lmer
##
## The following object is masked from 'package:stats':
##
## stepmodglmer <- glmer(Developpement ~ Taille*Temperature + (1|Labo), family=binomial(logit), data=Aubergines) # Réalise le modèle avec effets mixtes.
modglmer2 <-glmer(Developpement ~ Taille + (1|Labo), family=binomial(logit), data=Aubergines) # Réalise le modèle avec effets mixtes.
summary(modglmer2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Developpement ~ Taille + (1 | Labo)
## Data: Aubergines
##
## AIC BIC logLik deviance df.resid
## 238.6 251.3 -115.3 230.6 172
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3693 -1.0761 0.0443 0.9293 1.5584
##
## Random effects:
## Groups Name Variance Std.Dev.
## Labo (Intercept) 0 0
## Number of obs: 176, groups: Labo, 8
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.629 0.310 2.03 0.04230 *
## Taillemoyenne -0.482 0.381 -1.27 0.20539
## Taillepetite -1.516 0.443 -3.42 0.00063 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tllmyn
## Taillemoynn -0.813
## Taillepetit -0.698 0.568Anova(modglmer) # Affiche les résultats## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Developpement
## Chisq Df Pr(>Chisq)
## Taille 11.03 2 0.004 **
## Temperature 1.11 2 0.574
## Taille:Temperature 3.37 4 0.499
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(modglmer2) # Affiche les résultats## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Developpement
## Chisq Df Pr(>Chisq)
## Taille 12.4 2 0.0021 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(modglmer) # Donne la valeur du critère d'information d'Akaike du modèle.## [1] 246plot(residuals(modglmer),cex=1.5)
points(residuals(modglm),col=2)
Un exercice avec PPCB
On peut faire la courbe ROC - package dédié, on calcul l’air sous la courbe pour valider l’ajustement.
ppcb <- read.table(file="Exercices/PPCB.txt",
header = TRUE,
dec = ".",
sep = "\t")
summary(ppcb)## Troupeau Période Incidence
## trp_01 : 4 p_1:15 NON:22
## trp_03 : 4 p_2:14 OUI:34
## trp_04 : 4 p_3:14
## trp_05 : 4 p_4:13
## trp_06 : 4
## trp_07 : 4
## (Other):32library(lattice)
levelplot(as.numeric(Incidence) ~ Période + Troupeau , data = ppcb)
Table<-xtabs(~Incidence+Période+Troupeau, data=ppcb) # Réalise la table de contingence du nombre de cas en fonction de la taille
plot(Table, main="Nombre de cas en fonction Troupeau+Période ")# Réalise un graphique de la table de contingence créée ci-dessus
modglm <- glm(Incidence ~ Période, family=binomial(logit), data=ppcb) # Réalise le modèle avec effets fixes.
Anova(modglm) # Affiche les résultats## Analysis of Deviance Table (Type II tests)
##
## Response: Incidence
## LR Chisq Df Pr(>Chisq)
## Période 13 3 0.0046 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(modglm,type="III") # avec le type III, le calcul met plus de poids sur l'interaction ## Analysis of Deviance Table (Type III tests)
##
## Response: Incidence
## LR Chisq Df Pr(>Chisq)
## Période 13 3 0.0046 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modglmer <- glmer(I(Incidence == "OUI") ~ Période + (1|Troupeau), family=binomial(logit), data=ppcb) # Réalise le modèle avec effets mixtes.
VarCorr(modglmer)## Groups Name Std.Dev.
## Troupeau (Intercept) 0summary(modglmer)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: I(Incidence == "OUI") ~ Période + (1 | Troupeau)
## Data: ppcb
##
## AIC BIC logLik deviance df.resid
## 72.0 82.2 -31.0 62.0 51
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.742 -0.809 0.267 0.745 1.265
##
## Random effects:
## Groups Name Variance Std.Dev.
## Troupeau (Intercept) 0 0
## Number of obs: 56, groups: Troupeau, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.64 1.04 2.55 0.0108 *
## Périodep_2 -2.05 1.18 -1.75 0.0811 .
## Périodep_3 -2.93 1.17 -2.51 0.0122 *
## Périodep_4 -3.11 1.18 -2.63 0.0085 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:## Warning: abbreviate utilisé avec des caractères non ASCII## (Intr) Pérd_2 Pérd_3
## Périodep_2 -0.880
## Périodep_3 -0.887 0.780
## Périodep_4 -0.876 0.771 0.777Anova(modglmer) # Affiche les résultats## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: I(Incidence == "OUI")
## Chisq Df Pr(>Chisq)
## Période 8.23 3 0.042 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Conclusion, l’effet aléatoire ne fonctionne pas… à creuser sur l’ajusteùent
3 Analyse des plans en blocks complet
MAintenant on s’interesse aux données par block. Discussion sur la variabilité des blocks. Qu’est-ce qu’un block, comment l’insérer dans le modèle, en fixe, en aléatoire. Cela dépend de la question. Si on veut avoir une réponse moyenne ou si on veut prendre en compte dans le modèle la variabilité dans le modèle
library(nlme) # Charge la bibliothèque
In <- read.table(file="Exemples/Incubateur.txt",
header = TRUE,
dec = ".",
sep = "\t")
modlmeComp<-lme(PrctCroiss~Traitement, random=~1|Exp, data=In) # Réalise le modèle avec effets mixtes.
Anova(modlmeComp)# Affiche les résultats## Analysis of Deviance Table (Type II tests)
##
## Response: PrctCroiss
## Chisq Df Pr(>Chisq)
## Traitement 6.76 3 0.08 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ici un effet expérimentateur et dispositif (serre). On peut prendre en compte l’effet un à un ou avec interaction
library(nlme) # Charge la bibliothèque
choux <- read.table(file="Exemples/Choux.txt",
header = TRUE,
dec = ".",
sep = "\t")
boxplot(mesure ~ var, data=choux)
modfix <- lm(mesure ~ var, data=choux) # Réalise le modèle avec effets fixes seuls.
Anova(modfix) # Affiche les résultats.## Anova Table (Type II tests)
##
## Response: mesure
## Sum Sq Df F value Pr(>F)
## var 2282 3 2.8 0.051 .
## Residuals 11936 44
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modfix.resid<-rstandard(modfix) # Stocke les résidus du modèle dans un objet.
plot(modfix.resid~choux$Serre, main="variation selon la serre")# Réalise le graphique des résidus selon les modalités de 'Serre'.
abline(0,0,lty=2) # Ajoute une droite horizontale passant par le 0 des ordonnées.
plot(modfix.resid~choux$Exp, main="variation selon l'expérimentateur")# Réalise le graphique des résidus selon les modalités de 'Exp'.
abline(0,0,lty=2) # Ajoute une droite horizontale passant par le 0 des ordonnées.
# examinons les interactions possibes
xyplot(mesure~var, group = Serre, data = choux, type=c("p"), col=c(1,2))
library(ggplot2)
ggplot(choux,aes(var,mesure,col=Serre)) + geom_point() 
interaction.plot(choux$var,choux$Exp,choux$mesure)
interaction.plot(choux$var,choux$Serre,choux$mesure)
# Il y a un effet experimentateur fort
xyplot(mesure ~ var | Exp,
group = Serre,
type = c("p","a"),
data=choux)
library(lmerTest) # Charge la bibliothèque.
modlmer1<-lmer(mesure~var+(1|Serre), data=choux) # Réalise le modèle avec l'effet 'Serre' comme variable aléatoire uniquement et sans interaction avec l'effet 'var'.
modlmer2<-lmer(mesure~var+(1|Exp), data=choux) # Réalise le modèle avec l'effet 'Exp' comme variable aléatoire uniquement et sans interaction avec l'effet 'var'.
modlmer3<-lmer(mesure~var+(1|Serre)+(1|Exp),
data=choux) # Réalise le modèle avec les effets 'Exp' et 'Serre' comme variables aléatoires et sans interaction avec l'effet 'var'.
AIC(modlmer1, modlmer2, modlmer3) # Donne le critère d'Akaike pour les trois modèles. ## df AIC
## modlmer1 6 392.0
## modlmer2 6 392.9
## modlmer3 7 393.5anova(modlmer1, modlmer2, modlmer3) # Donne le critère d'Akaike pour les trois modèles. ## refitting model(s) with ML (instead of REML)## Warning: convergence code 3 from bobyqa: bobyqa -- a trust region step
## failed to reduce q## Data: choux
## Models:
## object: mesure ~ var + (1 | Serre)
## ..1: mesure ~ var + (1 | Exp)
## ..2: mesure ~ var + (1 | Serre) + (1 | Exp)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 6 412 424 -200 400
## ..1 6 413 424 -200 401 0.00 0 1.00
## ..2 7 414 427 -200 400 0.81 1 0.37modlmer4<-lmer(mesure~var+(var|Exp)+(1|Serre),
data=choux) # Réalise le modèle avec les effets 'Exp' et 'Serre' comme variables aléatoires et avec interaction entre l'effet 'Exp' et l'effet 'var'.
AIC(modlmer4) # Donne le critère d'Akaike du modèle.## [1] 385.5Anova(modlmer4) # Donne les résultats du modèle.## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: mesure
## Chisq Df Pr(>Chisq)
## var 16.8 3 0.00078 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmer4.resid <- residuals(modlmer4) # Stocke les résidus du modèle dans un objet.
qqnorm(lmer4.resid)
qqline(lmer4.resid) # Permet de vérifier graphiquement la normalité des résidus.
shapiro.test(lmer4.resid) # Teste la normalité des résidus.##
## Shapiro-Wilk normality test
##
## data: lmer4.resid
## W = 0.99, p-value = 0.9predites<-predict(modlmer4) # Stocke les valeurs prédites par le modèle dans un objet.
library(lattice) # Charge la bibliothèque.
xyplot(choux$mesure~predites, type=c("p","r"),
main="mesures observées et mesures prédites")# Réalise le graphique des valeurs observées en fonction des valeurs prédites par le modèle.
modfix.resid<-rstandard(modfix) # Stocke les résidus du modèle dans un objet.
plot(modfix.resid~choux$Serre, main="variation selon la serre")# Réalise le graphique des résidus selon les modalités de 'Serre'.
abline(0,0,lty=2) # Ajoute une droite horizontale passant par le 0 des ordonnées.
plot(modfix.resid~choux$Exp, main="variation selon l'expérimentateur")# Réalise le graphique des résidus selon les modalités de 'Exp'.
Pour aller plus loin, on pourrait regarder ?
4 Analyse des plans avec effets emboîtés
Vive les mycoses Test ici des différents modèles avec emboitement On ne trouve pas de transformation
library(nlme) # Charge la bibliothèque
mycose <- read.table(file="Exemples/Mycose.txt",
header = TRUE,
dec = ".",
sep = "\t")
summary(mycose) # Donne un résumé du jeu de données.## sujet semaine jour traitement mesure
## S1 : 6 Sem1:45 lundi :30 drogue :48 Min. :293
## S10 : 6 Sem2:45 mercredi:30 placebo:42 1st Qu.:313
## S11 : 6 vendredi:30 Median :326
## S12 : 6 Mean :324
## S13 : 6 3rd Qu.:336
## S14 : 6 Max. :345
## (Other):54boxplot(mesure ~ traitement+semaine,mycose)
qqmath(~ ((mesure^6.3)-1)/(6.3)|traitement,data=mycose,type=c("p","r"))
qqmath(~ (mesure)|traitement,data=mycose,type=c("p","r"))
xyplot(mesure ~ traitement | semaine+jour,
data = mycose,
type=c("p","a"))
xyplot(mesure ~ semaine:jour | traitement, group = sujet,
data = mycose,
type=c("p","a"))
var.test(mesure ~ traitement,alternative='two.sided', conf.level=.95, data=mycose) # Réalise le test d'homogénéité des variances entre les deux groupes.##
## F test to compare two variances
##
## data: mesure by traitement
## F = 1.1, num df = 47, denom df = 41, p-value = 0.7
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.6179 2.0522
## sample estimates:
## ratio of variances
## 1.133mod.fix<-lm(mesure~traitement, data=mycose) # Réalise le modèle avec effets fixes uniquement.
Anova(mod.fix) # Donne le résultat du modèle.## Anova Table (Type II tests)
##
## Response: mesure
## Sum Sq Df F value Pr(>F)
## traitement 3221 1 23.6 5.2e-06 ***
## Residuals 12038 88
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1resid.fix<-residuals(mod.fix) # Stocke les résidus du modèle dans un objet.
boxplot(resid.fix~mycose$sujet, main ="Résidus selon 'sujet'") # Réalise le graphique des résidus du modèle en fonction des modalités 'Sujet'.
boxplot(resid.fix~mycose$semaine, main ="Résidus selon 'semaine'") # Réalise le graphique des résidus du modèle en fonction des modalités 'Semaine'.
boxplot(resid.fix~mycose$jour, main ="Résidus selon 'jour'") # Réalise le graphique des résidus du modèle en fonction des modalités 'Jour'.
m1<-lmer(mesure~traitement+(1|sujet),data=mycose) # Réalise le modèle avec la variable 'sujet' comme variable aléatoire, sans interaction avec l'effet 'traitement'. Ce modèle est le meilleur trouvé pour ces données.
Anova(m1) # Donne le résultat du modèle.## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: mesure
## Chisq Df Pr(>Chisq)
## traitement 5.41 1 0.02 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1resid.m1 <- residuals(m1) # Stocke les résidus du modèle dans un objet.
qqnorm(resid.m1, main="Résidus du modèle m1")
qqline(resid.m1) # Réalise un graphique des résidus pour estimer graphiquement leur normalité.
shapiro.test(resid.m1) # Réalise le test de normalité des résidus.##
## Shapiro-Wilk normality test
##
## data: resid.m1
## W = 0.96, p-value = 0.006boxplot(resid.m1~mycose$sujet) # Réalise un graphique des résidus du modèle 'm1' en fonction des modalités 'Sujet'.
mycose$pred.m1<-predict(m1) # Stocke les valeurs prédites par le modèle dans une nouvelle colonne du jeu de données 'mycose'.
mycose$tps <- as.numeric(as.factor(paste(mycose$semaine,mycose$jour,sep="-")))
mcomplet<-lmer(mesure~traitement*tps+(1+tps|sujet),data=mycose) # ce modeèle n'est pas mieuxIl y a un gros pb de non normalité à traiter par des transformations
5 Exercices
5.1 Exercice 1
data("Wheat")
summary(Wheat)## Tray Moisture fertilizer DryMatter
## 3 : 4 Min. :10.0 Min. :2.0 Min. : 1.98
## 1 : 4 1st Qu.:17.5 1st Qu.:3.5 1st Qu.: 5.24
## 2 : 4 Median :25.0 Median :5.0 Median : 8.12
## 4 : 4 Mean :25.0 Mean :5.0 Mean : 8.45
## 5 : 4 3rd Qu.:32.5 3rd Qu.:6.5 3rd Qu.:10.81
## 6 : 4 Max. :40.0 Max. :8.0 Max. :15.71
## (Other):24str(Wheat)## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame': 48 obs. of 4 variables:
## $ Tray : Ord.factor w/ 12 levels "3"<"1"<"2"<"4"<..: 2 2 2 2 3 3 3 3 1 1 ...
## $ Moisture : num 10 10 10 10 10 10 10 10 10 10 ...
## $ fertilizer: num 2 4 6 8 2 4 6 8 2 4 ...
## $ DryMatter : num 3.35 4.32 4.56 5.88 4.04 ...
## - attr(*, "formula")=Class 'formula' length 3 DryMatter ~ fertilizer | Tray
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## - attr(*, "labels")=List of 2
## ..$ y: chr "Dry matter"
## ..$ x: chr "Fertlizer"
## - attr(*, "units")=List of 2
## ..$ y: chr "(%)"
## ..$ x: chr "(mg)"boxplot(DryMatter~Tray,Wheat)
xyplot(DryMatter ~ fertilizer | Moisture, group = Tray,
data = Wheat,
type=c("p","a"))
xyplot(DryMatter ~ fertilizer, group = Tray,
data = Wheat,
type=c("p","a"))
boxplot(Wheat$DryMatter)
# un modèle simple
mod.fix<-lm(DryMatter ~ fertilizer, data=Wheat) # Réalise le modèle avec effets fixes uniquement.
residualPlot(mod.fix)
resid.fix<-residuals(mod.fix) # Stocke les résidus du modèle dans un objet.
boxplot(resid.fix~Wheat$Tray, main ="Résidus selon 'sujet'") # Réalise le graphique des résidus du modèle en fonction des modalités 'Sujet'.
m1<-lmer(DryMatter ~ fertilizer+(1|Tray),data=Wheat) # Réalise le modèle avec la variable 'sujet' comme variable aléatoire, sans interaction avec l'effet 'traitement'. Ce modèle est le meilleur trouvé pour ces données.
AIC(mod.fix) # Donne le résultat du modèle.## [1] 238.4AIC(m1) # Donne le résultat du modèle.## [1] 202.6resid.m1 <- residuals(m1) # Stocke les résidus du modèle dans un objet.
boxplot(resid.m1~Wheat$Tray, main ="Résidus selon 'sujet'") # Réalise le graphique des résidus du modèle en fonction des modalités 'Sujet'.
qqnorm(resid.m1, main="Résidus du modèle m1")
qqline(resid.m1) # Réalise un graphique des résidus pour estimer graphiquement leur normalité.
shapiro.test(resid.m1) # Réalise le test de normalité des résidus.##
## Shapiro-Wilk normality test
##
## data: resid.m1
## W = 0.98, p-value = 0.5boxplot(resid.m1~Wheat$Tray) # Réalise un graphique des résidus du modèle 'm1' en fonction des modalités 'Sujet'.
# toujours pas bon
m2 <- lmer(DryMatter ~ fertilizer+(fertilizer|Tray),data=Wheat) # Réalise le modèle avec la variable 'sujet' comme variable aléatoire, sans interaction avec l'effet 'traitement'. Ce modèle est le meilleur trouvé pour ces données.
AIC(m2) # Donne le résultat du modèle.## [1] 185.5resid.m2 <- residuals(m2) # Stocke les résidus du modèle dans un objet.
boxplot(resid.m2~Wheat$Tray, main ="Résidus selon 'sujet'") # Réalise le graphique des résidus du modèle en fonction des modalités 'Sujet'.
qqnorm(resid.m2, main="Résidus du modèle m1")
qqline(resid.m2) # Réalise un graphique des résidus pour estimer graphiquement leur normalité.
shapiro.test(resid.m2) # Réalise le test de normalité des résidus.##
## Shapiro-Wilk normality test
##
## data: resid.m2
## W = 0.97, p-value = 0.2il reste un effet de variance qui augmente car la fonction est plutôt exponentielle il faudrait prendre en compte dans le modèle avec SSlogis see it demain
5.2 Exercice 2
## to be done6 Modèle non linéaire à effets fixes
On part sur un exemple
Les relations non linéaires disponibles sont décrites dans l’annexe 3
data("CO2")
library(car)
scatterplot(uptake ~conc , data = CO2 , reg.line = FALSE)
On remarque bien une reation non linéaire. On va ajuster une fonction asymptotique. 2 étape, on calcule d’abord les valeurs initiales puis on ajuste. Cela se fait sur les paramètres de la fonction non linéaire
vecStart <- getInitial(uptake ~ SSasymp( input = conc , Asym, R0 , lrc ), CO2) # Bien mettre les noms des paramètres comme défini , ici Asym, R0, lrc
nl1 <- nls(uptake ~ SSasymp(conc , Asym , R0 , lrc ) , start = vecStart , data = CO2)
CO2$fitted <- predict(nl1)
ggplot( CO2 , aes( conc, uptake , col= Treatment) ) + geom_point() +
geom_point(aes(y=fitted), col = 3, size= 5) +
geom_line(aes(y = fitted))
# regardons la qualité du modèle
AIC(nl1)## [1] 600.5hist(residuals(nl1))
qqnorm(residuals(nl1))
plot(nl1) # forte hétéroscedasticité
Effet du paramètre traitement, il faut séparer le jeu de données en 2 et ajuster un modèle par sous jeu
d’abord une représentation graphique
Vincent parle de la méthode gnls
ggplot( CO2 , aes( conc, uptake , col= Treatment) ) + geom_point() +
geom_smooth(method = "gam", formula = y ~ poly(x,2), size = 1) +
geom_point(aes(y=fitted), col = 3, size= 5) +
geom_line(aes(y = fitted))
CO2.CH<-CO2[CO2$Treatment=="chilled",] # Crée un sous-fichier ne contenant que les valeurs "chilled"
CO2.NCH<-CO2[CO2$Treatment=="nonchilled",] # Crée un sous-fichier ne contenant que les valeurs "non chilled"
nlmCH<- nls(uptake ~ SSasymp(conc, Asym, R0,
lrc), data = CO2.CH) # Réalise le modèle pour les données « chilled »
nlmNCH<-update(nlmCH, data=CO2.NCH) # Réalise le modèle pour les données « non-chilled »
summary(nlmCH) ; summary(nlmNCH) # affiche les résultats des paramètres du modèle.##
## Formula: uptake ~ SSasymp(conc, Asym, R0, lrc)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Asym 28.388 2.374 11.96 1.3e-14 ***
## R0 -11.786 23.079 -0.51 0.61
## lrc -4.718 0.579 -8.14 6.1e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.36 on 39 degrees of freedom
##
## Number of iterations to convergence: 0
## Achieved convergence tolerance: 1.21e-07##
## Formula: uptake ~ SSasymp(conc, Asym, R0, lrc)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Asym 36.554 1.340 27.28 <2e-16 ***
## R0 -23.811 16.304 -1.46 0.15
## lrc -4.610 0.251 -18.35 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.52 on 39 degrees of freedom
##
## Number of iterations to convergence: 0
## Achieved convergence tolerance: 7.06e-08t.test(fitted(nlmCH),fitted(nlmNCH)) # compare les valeurs prédites par les 2 groupes avec un test t.##
## Welch Two Sample t-test
##
## data: fitted(nlmCH) and fitted(nlmNCH)
## t = -4.4, df = 75, p-value = 3e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -9.938 -3.781
## sample estimates:
## mean of x mean of y
## 23.78 30.64boxplot(fitted(nlmCH), fitted(nlmNCH), names=c("CH","NCH"), main="Valeurs prédites par le modèle") # affiche les valeurs prédites sous forme de box-plot.
CO2$fittedCH <- predict(nlmCH)
CO2$fittedNCH<- predict(nlmNCH)
ggplot( CO2 , aes( conc, uptake , col= Treatment) ) + geom_point() +
geom_point(aes(y=fitted), col = 3, size= 5) +
geom_line(aes(y = fitted)) +
geom_line(aes(y = fittedCH ) , col = 4 ) +
geom_line(aes(y = fittedNCH ), col = 5 )
library(nlme) # charge la bibliothèque.
res.traitement<-nlsList(uptake ~ SSasymp(conc, Asym, R0, lrc)|Treatment, start=vecStart, data = CO2) # fournit les paramètres du modèle selon les modalités de la colonne 'Traitement'.
intervals(res.traitement, level=0.95) # Donne les intervalles de confiance à 95% de ces paramètres.## , , Asym
##
## lower est. upper
## nonchilled 32.84 36.55 40.27
## chilled 24.51 28.39 32.27
##
## , , R0
##
## lower est. upper
## nonchilled -69.02 -23.81 21.40
## chilled -49.50 -11.79 25.92
##
## , , lrc
##
## lower est. upper
## nonchilled -5.306 -4.610 -3.913
## chilled -5.665 -4.718 -3.771lots<-nlsList(uptake ~ SSasymp(conc, Asym, R0, lrc)| Plant, start=vecStart, data = CO2) # calcule les paramètres du modèle selon les lots de la colonne 'Plant'.
lots # affiche les résultats## Call:
## Model: uptake ~ SSasymp(conc, Asym, R0, lrc) | Plant
## Data: CO2
##
## Coefficients:
## Asym R0 lrc
## Qn1 38.14 -34.276 -4.381
## Qn2 42.87 -29.656 -4.666
## Qn3 44.23 -37.627 -4.486
## Qc1 36.43 -9.901 -4.862
## Qc3 40.68 -11.542 -4.945
## Qc2 39.82 -51.535 -4.464
## Mn3 28.48 -17.371 -4.592
## Mn2 32.13 -29.044 -4.466
## Mn1 34.08 -8.813 -5.065
## Mc2 13.56 -1.982 -4.561
## Mc3 18.54 -136.114 -3.465
## Mc1 21.79 2.449 -5.142
##
## Degrees of freedom: 84 total; 48 residual
## Residual standard error: 1.798plot(intervals(lots)) # représente graphiquement les intervalles de confiance autour des paramètres du modèle.
vecStart.mic<-getInitial(uptake~SSmicmen(conc, Vm, K), data=CO2) # Permet d'obtenir les valeurs initiales des paramètres du modèle et de les stocker dans un vecteur.
print(vecStart.mic) # Affiche le contenu du vecteur## Vm K
## 39.86 140.33nl1<- nls(uptake ~ SSasymp(conc, Asym, R0, lrc), start=vecStart, data = CO2) # réalise le modèle avec la fonction 'asymp'.
nl2<-nls(uptake~SSmicmen(conc, Vm, K), data=CO2) # réalise le modèle avec la fonction de Micaelis-Menten.
AIC(nl1,nl2) # affiche les valeurs du critère d'Akaike pour les 2 modèles## df AIC
## nl1 4 600.5
## nl2 3 603.6anova(nl1,nl2) # réalise la comparaison des 2 modèles## Analysis of Variance Table
##
## Model 1: uptake ~ SSasymp(conc, Asym, R0, lrc)
## Model 2: uptake ~ SSmicmen(conc, Vm, K)
## Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
## 1 81 5691
## 2 82 6044 -1 -353 5.03 0.028 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 16.1 Essai linéaire avec poly
NE fonctionne pas, pris exemple dans le pinhiero
fm1Dial.lme <-
lme(uptake ~ ((conc) + (conc)^2 + (conc)^4 + (conc)^3) * Treatment,
CO2, ~ conc + conc^2 )
fm1Dial.lme <-
lme(uptake ~ poly(conc,1) * Treatment,
CO2, ~ conc + conc^2 )
fm1Dial.lme## Linear mixed-effects model fit by REML
## Data: CO2
## Log-restricted-likelihood: -268.4
## Fixed: uptake ~ poly(conc, 1) * Treatment
## (Intercept) poly(conc, 1)
## 30.64 53.45
## Treatmentchilled poly(conc, 1):Treatmentchilled
## -6.86 -11.29
##
## Random effects:
## Formula: ~conc + conc^2 | Plant
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 4.270093 (Intr)
## conc 0.006891 0.906
## Residual 5.674525
##
## Number of Observations: 84
## Number of Groups: 12anova( fm1Dial.lme )## numDF denDF F-value p-value
## (Intercept) 1 70 130.50 <.0001
## poly(conc, 1) 1 70 37.48 <.0001
## Treatment 1 10 2.19 0.1697
## poly(conc, 1):Treatment 1 70 0.52 0.4721plot( augPred(fm1Dial.lme), grid = T )
plot( fm1Dial.lme, resid(.) ~ conc, abline = 0 )
fm2Dial.lme <- update( fm1Dial.lme,
weights = varPower(form = ~ conc) )
anova( fm1Dial.lme, fm2Dial.lme )## Model df AIC BIC logLik Test L.Ratio p-value
## fm1Dial.lme 1 8 552.9 572.0 -268.4
## fm2Dial.lme 2 9 465.1 486.5 -223.5 1 vs 2 89.84 <.0001plot( fm2Dial.lme, resid(., type = "p") ~ conc,
abline = 0 )
plot(fm2Dial.lme, resid(.) ~ conc|Treatment, abline = 0)
fm3Dial.lme <- update(fm2Dial.lme,
weights=varPower(form = ~ conc | Treatment))
anova( fm2Dial.lme, fm3Dial.lme )## Model df AIC BIC logLik Test L.Ratio p-value
## fm2Dial.lme 1 9 465.1 486.5 -223.5
## fm3Dial.lme 2 10 467.0 490.8 -223.5 1 vs 2 0.1023 0.7491fm4Dial.lme <- update(fm2Dial.lme,
weights = varConstPower(form = ~ conc))
anova( fm2Dial.lme, fm4Dial.lme )## Model df AIC BIC logLik Test L.Ratio p-value
## fm2Dial.lme 1 9 465.1 486.5 -223.5
## fm4Dial.lme 2 10 467.1 490.9 -223.5 1 vs 2 6.1e-06 0.998plot( augPred(fm2Dial.lme), grid = T )
7 Modèles non linéaires à effets mixtes
plot(CO2) # réalise un graphique par plant.
library(nlme) # charge la bibliothèque.
nlmCO2<-nlme(uptake~SSasymp(conc,Asym,R0,lrc),
data = CO2,
fixed=Asym + R0 +lrc ~ 1 ,
random=Asym ~1 | Plant,
start = c(Asym = 32.46, R0 = -17.82,
lrc = -4.65)) # réalise le modèle
summary(nlmCO2) # affiche les résultats du modèle## Nonlinear mixed-effects model fit by maximum likelihood
## Model: uptake ~ SSasymp(conc, Asym, R0, lrc)
## Data: CO2
## AIC BIC logLik
## 445.4 457.6 -217.7
##
## Random effects:
## Formula: Asym ~ 1 | Plant
## Asym Residual
## StdDev: 8.94 2.365
##
## Fixed effects: Asym + R0 + lrc ~ 1
## Value Std.Error DF t-value p-value
## Asym 32.94 2.663 70 12.37 0.0000
## R0 -10.89 3.311 70 -3.29 0.0016
## lrc -4.81 0.080 70 -60.04 0.0000
## Correlation:
## Asym R0
## R0 0.060
## lrc -0.099 -0.895
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.9543 -0.5434 -0.0698 0.5257 2.2890
##
## Number of Observations: 84
## Number of Groups: 12CO2$pred<-predict(nlmCO2) # Stocke les valeurs prédites par le modèle dans une nouvelle colonne du jeu de données.
plot(nlmCO2) # Permet d'observer graphiquement les résidus en fonction des valeurs prédites.
qqnorm(residuals(nlmCO2)); qqline(residuals(nlmCO2)) # Permet d'observer graphiquement la normalité des résidus.
shapiro.test(residuals(nlmCO2)) # Réalise le test de normalité des résidus.##
## Shapiro-Wilk normality test
##
## data: residuals(nlmCO2)
## W = 0.98, p-value = 0.2nl1 ; nlmCO2 # Affiche les résultats des deux modèles.## Nonlinear regression model
## model: uptake ~ SSasymp(conc, Asym, R0, lrc)
## data: CO2
## Asym R0 lrc
## 32.46 -17.82 -4.65
## residual sum-of-squares: 5691
##
## Number of iterations to convergence: 0
## Achieved convergence tolerance: 1.45e-07## Nonlinear mixed-effects model fit by maximum likelihood
## Model: uptake ~ SSasymp(conc, Asym, R0, lrc)
## Data: CO2
## Log-likelihood: -217.7
## Fixed: Asym + R0 + lrc ~ 1
## Asym R0 lrc
## 32.943 -10.893 -4.813
##
## Random effects:
## Formula: Asym ~ 1 | Plant
## Asym Residual
## StdDev: 8.94 2.365
##
## Number of Observations: 84
## Number of Groups: 12AIC(nl1, nlmCO2) # Affiche le critère d'Akaike des deux modèles.## df AIC
## nl1 4 600.5
## nlmCO2 5 445.4nlm2CO2<-update(nlmCO2,random=Asym+R0 ~1|Plant) # réalise le modèle avec les effets aléatoires affectant Asym et R0
nlm3CO2<-update(nlmCO2,random=Asym+R0+lrc ~1|Plant) # réalise le modèle avec les effets aléatoires affectant tous les paramètres du modèle
anova(nlmCO2, nlm2CO2, nlm3CO2) # Compare les trois modèles obtenus.## Model df AIC BIC logLik Test L.Ratio p-value
## nlmCO2 1 5 445.4 457.6 -217.7
## nlm2CO2 2 7 417.0 434.0 -201.5 1 vs 2 32.41 <.0001
## nlm3CO2 3 10 420.9 445.2 -200.4 2 vs 3 2.13 0.5464anova(nlm2CO2, nlm3CO2)# Compare les deux modèles ayant des AIC voisins.## Model df AIC BIC logLik Test L.Ratio p-value
## nlm2CO2 1 7 417.0 434.0 -201.5
## nlm3CO2 2 10 420.9 445.2 -200.4 1 vs 2 2.127 0.5464nlm4CO2<-update(nlmCO2,fixed=Asym + R0 +lrc ~ Treatment + Type, start = c(Asym = c(32.46,32.46,32.46), R0 = c(-17.82,-17,-17),
lrc = c(-4.65,-4,-4))) # réalise le modèle avec les effets
anova(nlmCO2, nlm4CO2)# Compare les deux modèles ayant des AIC voisins.## Model df AIC BIC logLik Test L.Ratio p-value
## nlmCO2 1 5 445.4 457.6 -217.7
## nlm4CO2 2 11 400.5 427.2 -189.2 1 vs 2 56.91 <.0001anova(nlmCO2, nlm2CO2, nlm3CO2,nlm4CO2) # Compare les trois modèles obtenus.## Model df AIC BIC logLik Test L.Ratio p-value
## nlmCO2 1 5 445.4 457.6 -217.7
## nlm2CO2 2 7 417.0 434.0 -201.5 1 vs 2 32.41 <.0001
## nlm3CO2 3 10 420.9 445.2 -200.4 2 vs 3 2.13 0.5464
## nlm4CO2 4 11 400.5 427.2 -189.2 3 vs 4 22.38 <.0001comp<-compareFits(coef(nlm2CO2), coef(nlm3CO2)) # Compare les paramètres des 2 modèles par plant
comp # Affiche les résultats## , , Asym
##
## coef(nlm2CO2) coef(nlm3CO2)
## Qn1 39.26 38.69
## Qn2 42.60 42.61
## Qn3 45.01 44.44
## Qc1 35.75 35.70
## Qc3 39.26 39.45
## Qc2 40.06 40.25
## Mn3 28.62 28.60
## Mn2 32.61 32.48
## Mn1 32.02 32.46
## Mc2 13.70 13.61
## Mc3 19.61 19.63
## Mc1 20.90 20.92
##
## , , R0
##
## coef(nlm2CO2) coef(nlm3CO2)
## Qn1 -22.6312 -29.838
## Qn2 -28.6294 -31.010
## Qn3 -29.0995 -35.493
## Qc1 -20.7704 -23.634
## Qc3 -24.9732 -26.788
## Qc2 -27.6318 -27.434
## Mn3 -14.2961 -15.300
## Mn2 -18.1093 -20.386
## Mn1 -19.2061 -17.119
## Mc2 0.3434 1.865
## Mc3 -3.9687 -5.218
## Mc1 -6.4121 -5.789
##
## , , lrc
##
## coef(nlm2CO2) coef(nlm3CO2)
## Qn1 -4.65 -4.502
## Qn2 -4.65 -4.634
## Qn3 -4.65 -4.534
## Qc1 -4.65 -4.618
## Qc3 -4.65 -4.661
## Qc2 -4.65 -4.672
## Mn3 -4.65 -4.633
## Mn2 -4.65 -4.604
## Mn1 -4.65 -4.736
## Mc2 -4.65 -4.648
## Mc3 -4.65 -4.634
## Mc1 -4.65 -4.670plot(comp) # Représente graphiquement les paramètres pour les 2 modèles sur le même graphique.
plot(nlm4CO2) # en non linéaire, ce"
plot( augPred(nlm4CO2), grid = T ) # graphique des prédictions du MMixte non linéaire
predict(nlm4CO2)## Qn1 Qn1 Qn1 Qn1 Qn1 Qn1 Qn1 Qn2 Qn2 Qn2
## 13.112 28.756 34.881 38.012 39.243 39.485 39.523 14.611 30.722 37.029
## Qn2 Qn2 Qn2 Qn2 Qn3 Qn3 Qn3 Qn3 Qn3 Qn3
## 40.253 41.521 41.770 41.810 16.199 32.805 39.306 42.629 43.936 44.193
## Qn3 Qc1 Qc1 Qc1 Qc1 Qc1 Qc1 Qc1 Qc2 Qc2
## 44.233 11.782 24.131 30.058 33.829 35.832 36.424 36.582 13.645 26.725
## Qc2 Qc2 Qc2 Qc2 Qc2 Qc3 Qc3 Qc3 Qc3 Qc3
## 33.004 36.999 39.121 39.748 39.915 13.447 26.449 32.690 36.662 38.771
## Qc3 Qc3 Mn1 Mn1 Mn1 Mn1 Mn1 Mn1 Mn1 Mn2
## 39.394 39.561 11.709 22.023 26.808 29.744 31.227 31.636 31.737 12.179
## Mn2 Mn2 Mn2 Mn2 Mn2 Mn2 Mn3 Mn3 Mn3 Mn3
## 22.669 27.538 30.524 32.033 32.449 32.551 10.156 19.881 24.394 27.163
## Mn3 Mn3 Mn3 Mc1 Mc1 Mc1 Mc1 Mc1 Mc1 Mc1
## 28.561 28.947 29.042 11.224 15.444 17.752 19.443 20.536 20.955 21.114
## Mc2 Mc2 Mc2 Mc2 Mc2 Mc2 Mc2 Mc3 Mc3 Mc3
## 8.005 10.756 12.260 13.362 14.074 14.348 14.451 10.726 14.720 16.903
## Mc3 Mc3 Mc3 Mc3
## 18.503 19.537 19.934 20.084
## attr(,"label")
## [1] "Fitted values (umol/m^2 s)"8 Exercices Non linéaire
8.1 Exercices 1
On procède d’abord à une selection des effets fixes sur la base des AIC et ensuite on cherche les meilleurs effets aléatoires.
data("Orange")
library(lattice)
library(nlme)
library(ggplot2)
library(car)
library(lme4)
xyplot(circumference ~ age, group = Tree, Orange,
type=c("p","a"))
hist(Orange$age)
hist(Orange$circumference)
vecStart <- getInitial(circumference ~ SSasymp( input = age , Asym, R0 , lrc ), Orange) # Bien mettre les noms des paramètres comme défini , ici Asym, R0, lrc
nl1 <- nls(circumference ~ SSasymp(age , Asym , R0 , lrc ) , start = vecStart , data = Orange)
plot(residuals(nl1))
Orange$fitted <- predict(nl1)
ggplot( Orange , aes( age, circumference , col= Tree) ) + geom_point() +
geom_point(aes(y=fitted), col = 3, size= 5) +
geom_line(aes(y = fitted))
AIC(nl1)## [1] 326.5nl2<-nlme(circumference ~ SSasymp(age , Asym , R0 , lrc ),
data = Orange,
fixed=Asym + R0 +lrc ~ 1 ,
random= lrc + R0 ~1 | Tree,
start = c(Asym = c(559.2826), R0 = c( 10.81873), lrc = c(-8.344031))) # réalise le modèle
AIC(nl2)##
## 287.2plot(nl2)
plot( augPred(nl2), grid = T ) # graphique des prédictions du MMixte non linéaire
cor(Orange$circumference,nl2$fitted[,2])^2## [1] 0.9768# essai avec SSgompertz
vecStart <- getInitial(circumference ~ SSgompertz( x = age , Asym, b2,b3 ), Orange) # Bien mettre les noms des paramètres comme défini , ici Asym, R0, lrc
nl1.g <- nls(circumference ~ SSgompertz( x = age , Asym, b2,b3 ) , start = vecStart , data = Orange)
plot(residuals(nl1.g))
nl2.g <-nlme(circumference ~ SSgompertz( x = age , Asym, b2,b3 ) ,
data = Orange,
fixed=Asym + b2 +b3 ~ 1 ,
random= Asym ~1 | Tree,
start = c(Asym = c(223.6533447), R0 = c( 2.5735171), lrc = c(0.9984465))) # réalise le modèle
AIC(nl2.g)##
## 276.2plot(nl2.g)
plot( augPred(nl2.g), grid = T ) # graphique des prédictions du MMixte non linéaire
library(nlme)
Orange$fix <- "1"
orange.modmix <- lme( circumference ~age,
random = ~ 1 + age | Tree,
data = Orange,
method = "ML")
orange.modmix1 <- lme( circumference ~age,
random = ~ 1 | fix,
data = Orange,
method = "ML")
plot( augPred(orange.modmix), grid = T ) # graphique des prédictions du MMixte non linéaire
anova(orange.modmix,orange.modmix1)## Model df AIC BIC logLik Test L.Ratio p-value
## orange.modmix 1 6 288.9 298.2 -138.4
## orange.modmix1 2 4 327.0 333.2 -159.5 1 vs 2 42.09 <.0001# modmix1 <- lmer(circumference~ age + (age|Tree), data = Orange, REML = FALSE)
AIC(orange.modmix)## [1] 288.9# cor(Orange$circumference,predict(modmix1))^2
cor(Orange$circumference,nl2$fitted[,2])^2## [1] 0.97688.2 exercice 2
data("ChickWeight")
library(lattice)
library(nlme)
library(ggplot2)
library(car)
library(lme4)
plot(ChickWeight)
xyplot(weight ~ Time | Diet, group = Chick, ChickWeight,
type=c("p","a"))
hist(Orange$age)
hist(Orange$circumference)
nlm2<-nlme(weight~SSlogis(Time,Asym, xmid, scal), data=ChickWeight,
fixed=Asym+ xmid+ scal ~1,
random=Asym +xmid ~1|Chick,
start=c(Asym=323.81,xmid=15.28,scal=7.82))
nlm2## Nonlinear mixed-effects model fit by maximum likelihood
## Model: weight ~ SSlogis(Time, Asym, xmid, scal)
## Data: ChickWeight
## Log-likelihood: -2217
## Fixed: Asym + xmid + scal ~ 1
## Asym xmid scal
## 285.06 11.53 6.69
##
## Random effects:
## Formula: list(Asym ~ 1, xmid ~ 1)
## Level: Chick
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## Asym 133.104 Asym
## xmid 4.644 0.915
## Residual 7.741
##
## Number of Observations: 578
## Number of Groups: 50