__AVIS IMPORTANT__ Depuis l'automne 2021, ce wiki a été discontinué et n'est plus activement développé. Tout le matériel mis à jour et les annonces pour la série d'ateliers R du CSBQ se trouvent maintenant sur le [[https://r.qcbs.ca/fr/workshops/r-workshop-08/|site web de la série d'ateliers R du CSBQ]]. Veuillez mettre à jour vos signets en conséquence afin d'éviter les documents périmés et/ou les liens brisés. Merci de votre compréhension, Vos coordonnateurs de la série d’ateliers R du CSBQ. ======= Ateliers R du CSBQ ======= [[http://qcbs.ca/fr/|{{:logo_text.png?nolink&500|}}]] Cette série de [[r|10 ateliers]] guide les participants à travers les étapes requises afin de maîtriser le logiciel R pour une grande variété d’analyses statistiques pertinentes en recherche en biologie et en écologie. Ces ateliers en libre accès ont été créés par des membres du CSBQ à la fois pour les membres du CSBQ et pour la grande communauté d’utilisateurs de R. //Le contenu de cet atelier a été révisé par plusieurs membres du CSBQ. Si vous souhaitez y apporter des modifications, veuillez SVP contacter les coordonnateurs actuels de la série, listés sur la page d'accueil// AVIS IMPORTANT: MISES À JOUR MAJEURES **Mise à jour de mars 2021:** Ce wiki n'est plus activement développé ou mis à jour. Le matériel mis à jour pour la série d'ateliers R du CSBQ est maintenant hébergé sur [[https://github.com/QCBSRworkshops/workshop08|la page GitHub]] des ateliers R du CSBQ. Le matériel disponible inclut; - La [[https://qcbsrworkshops.github.io/workshop08/pres-fr/workshop08-pres-fr.html|présentation Rmarkdown]] pour cet atelier; - Un [[https://qcbsrworkshops.github.io/workshop08/book-fr/workshop08-script-fr.R|script R]] qui suit la présentation - //*en construction*//. - [[https://qcbsrworkshops.github.io/workshop08/book-fr/index.html|Le matériel écrit]] qui accompagne la présentation en format bookdown - //*en construction*//. ====== Atelier 8 : Modèles additifs généralisés ====== Développé par : Eric Pedersen et Zofia Taranu Révision par : Cédric Frenette Dussault **Résumé:** L'objectif de l'atelier d'aujourd'hui sera d'examiner ce que nous entendons par un modèle non-linéaire et comment les GAMs (modèles additifs généralisés) nous permettent de modéliser les relations non-linéaires. Nous examinerons également comment tracer et interpréter ces relations non-linéaires, comment ajouter des interactions, comment prendre en compte la non-indépendance des données (//e.g.// erreurs autocorrélées) et comment inclure des effets aléatoires en se basant sur les ateliers précédents. Enfin, nous allons brièvement aborder la mécanique derrière le fonctionnement des GAMs. **Lien vers la nouvelle [[https://qcbsrworkshops.github.io/workshop08/workshop08-fr/workshop08-fr.html|présentation Rmarkdown]]** Lien vers la présentation Prezi associée : [[https://prezi.com/ddkjfsxlr-za/| Prezi]] Télécharger le script R et les données pour cette leçon: - [[http://qcbs.ca/wiki/_media/gam_e.r|Script]] - [[http://qcbs.ca/wiki/_media/other_dist.csv|Données]] ===== Objectifs de l'atelier ===== - Utiliser la librairie //mgcv// pour modéliser les relations non linéaires - Évaluer la sortie d'un GAM afin de mieux comprendre nos données - Utiliser des tests pour déterminer si nos relations correspondent à des modèles non linéaires ou linéaires - Ajouter des interactions non linéaires entre les variables explicatives - Comprendre l'idée d'une fonction de base (//basis function//) et la raison pour laquelle ça rend les GAMs si puissants ! - Comment modéliser la dépendance dans les données (autocorrélation, structure hiérarchique) en utilisant les GAMMs **Prérequis:** Expérience du logiciel R (assez pour être en mesure d'exécuter un script et d'examiner les données et les objets dans R) et une connaissance de base de la régression simple (vous devez savoir ce qu'on entend par une régression linéaire et une ANOVA). ===== 0. Le modèle linéaire ... et où il échoue ===== Qu'est-ce qu'un modèle linéaire ? La régression linéaire est ce que la plupart des gens apprennent avant tout en statistiques et est parmi les méthodes les plus performantes. Elle nous permet de modéliser une variable réponse en fonction de facteurs prédictifs et d'une erreur résiduelle. Tel que vu dans l'atelier sur les [[http://qcbs.ca/wiki/r_atelier4|modèles linéaires]], la régression fait cependant quatre suppositions importantes : - l'erreur est distribuée normalement - la variance des erreurs est constante - chaque erreur est indépendante des autres (homoscédasticité) - la réponse est linéaire: {y} = {β_0} + {β_1}{x} Il n'y a qu'une façon pour qu'un modèle linéaire soit correctement appliqué : {{::graphic0.1.png?350|}} et pourtant tant de façons pour qu'il ne le soit pas : {{::graphic0.2.png|}} Alors, comment résoudre ce problème ? Nous devons premièrement savoir ce que le modèle de régression cherche à faire : * ajuster une ligne qui passe au milieu des données, * faire cela sans sur-ajuster les données, c'est-à-dire en passant une ligne entre chaque point. Les modèles linéaires le font en trouvant la meilleure ligne droite qui passe à travers les données. En revanche, les modèles additifs font cela en ajustant une courbe à travers les données, mais tout en contrôlant le degré de courbure de la ligne //(plus d'information sur cela plus bas)//. ===== 1. Introduction aux GAMs ===== Examinons un exemple. Premièrement, nous allons générer des données et les représenter graphiquement. library(ggplot2) set.seed(10) n = 250 x = runif(n,0,5) y_model = 3*x/(1+2*x) y_obs = rnorm(n,y_model,0.1) data_plot = qplot(x, y_obs) + geom_line(aes(y=y_model)) + theme_bw() print(data_plot) {{::graphic1.1.jpg?350|}} Si nous modélisions cette relation par une régression linéaire, les résultats ne respecteraient pas les suppositions énumérées ci-dessus. Commençons par modéliser une régression en utilisant la méthode des moindres carrés en utilisant la fonction ''gam()'' de la librairie //mgcv// - donc en tant que modèle linéaire //(nous verrons plus bas comment utiliser cette fonction pour spécifier un terme non linéaire).// library(mgcv) linear_model = gam(y_obs~x) model_summary=summary(linear_model) print(model_summary) data_plot = data_plot+ geom_line(colour="red", aes(y=fitted(linear_model))) print(data_plot) Nous pouvons constater à partir du sommaire que notre modèle linéaire explique une grande partie de la variance (R²_(adj) = 0.639). Toutefois, les graphiques de diagnostic des résidus du modèle montrent que l'écart type ne suit pas une distribution normale et que la variance n'est pas homoscédastique. De plus, il reste un patron non-linéaire important. Essayons maintenant de résoudre ce problème en ajustant les données avec un terme non linéaire. Nous reviendrons sur ceci un peu plus tard, mais brièvement, les GAMs sont une forme non paramétrique de la régression où le beta*x_i d'une régression linéaire est remplacé par une fonction de lissage des variables explicatives, f (x_i), et le modèle devient : {y_i} = f(x_{i}) + {ε_i} où y_i est la variable réponse, x_i est la covariable, et //f// est la fonction lissage. Étant donné que la fonction de lissage f(x_i) est non linéaire et locale, l'ampleur de l'effet de la variable explicative peut varier en fonction de la relation entre la variable et la réponse. Autrement dit, contrairement à un coefficient fixe beta, la fonction //f// peut changer tout au long du gradient x_i. Le degré de lissage de //f// est contrôlée en utilisant une régression pénalisée qui est déterminée automatiquement à l'aide d'une méthode de validation croisée généralisée (GCV) de la librairie //mgcv// (Wood 2006). Avec ''gam()'' les termes non linéaires sont spécifiés par des expressions de la forme: ''s(x)''. gam_model = gam(y_obs~s(x)) summary(gam_model) data_plot = data_plot + geom_line(colour="blue",aes(y=fitted(gam_model))) print(data_plot) La variance expliquée par notre modèle a augmenté de 20% (R²_(adj) = 0.859) et quand on compare l'ajustement du modèle linéaire (rouge) au modèle non-linéaire (bleu), il est clair que l'ajustement de ce dernier est relativement meilleur. {{::graphic1.3.jpg?350|}} La librairie //mgcv// comprend également une fonction ''plot'' qui, par défaut, nous permet de visualiser la non-linéarité du modèle. plot(gam_model) Comment utilisons-nous les GAMs pour savoir si un modèle linéaire est suffisant pour modéliser nos données ? Nous pouvons utiliser les fonctions ''gam()'' et ''anova()'' pour tester formellement si une hypothèse de linéarité est justifiée. Nous devons simplement le configurer de sorte que notre modèle non-linéaire est emboîté dans notre modèle linéaire; c'est à dire, nous devons créer un objet qui inclut à la fois ''x'' (linéaire) et ''s(x)'' (non-linéaire). En utilisant la fonction ''anova()'', on vérifie si l'ajout de ''s(x)'' au modèle avec seulement ''x'' comme covariable est justifié par les données. linear_model = gam(y_obs~x) nested_gam_model = gam(y_obs~s(x)+x) print(anova(linear_model, nested_gam_model, test="Chisq")) Le terme non linéaire est significatif: Analysis of Deviance Table Model 1: y_obs ~ x Model 2: y_obs ~ s(x) + x Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 248.00 6.5846 2 240.68 2.4988 7.3168 4.0858 < 2.2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ---- **DÉFI 1** Nous allons maintenant essayer cela avec d'autres données générées aléatoirement. Nous allons d'abord générer les données. Ensuite, nous allons ajuster un modèle linéaire et un GAM à la relation entre ''x_test'' et ''y_test_obs''. Quels sont les degrés de libertés effectifs du terme non-linéaire ? Déterminez si l'hypothèse de linéarité est justifiée pour ces données. n <- 250 x_test <- runif(n,-5,5) y_test_fit <- 4*dnorm(x_test) y_test_obs <- rnorm(n,y_test_fit, 0.2) ++++ Réponse au défi 1 | data_plot <- qplot(x_test, y_test_obs) + geom_line(aes(y=y_test_fit))+ theme_bw() print(data_plot) linear_model_test <- gam(y_test_obs~x_test) nested_gam_model_test <- gam(y_test_obs~s(x_test)+x_test) print(anova(linear_model_test, nested_gam_model_test, test="Chisq")) summary(nested_gam_model_test)$s.table Analysis of Deviance Table Model 1: y_test_obs ~ x_test Model 2: y_test_obs ~ s(x_test) + x_test Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 248.0 81.09 2 240.5 7.46 7.5012 73.629 < 2.2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 edf Ref.df F p-value s(x_test) 7.602145 8.029057 294.0944 0 Réponse: Oui la non-linéarité est justifiée. Les degrés de libertés effectifs (edf) sont >> 1. {{::challenge_1_soln.jpg?350|}} ++++ ===== 2. Plusieurs termes non linéaires ===== Avec les GAMs, il est facile d'ajouter des termes non linéaires et linéaires dans un seul modèle, plusieurs termes non linéaires ou même des interactions non linéaires. Dans cette section, nous allons utiliser un ensemble de données générées automatiquement par //mgcv//. gam_data = gamSim(eg=5) head(gam_data) Nous allons voir comment nous pouvons prédire la variable réponse y en fonction des autres variables. Commençons par un modèle de base comprenant un terme non linéaire (x1) et un facteur qualitatif (X0 avec 4 niveaux). basic_model = gam(y~x0+s(x1), data= gam_data) basic_summary = summary(basic_model) print(basic_summary$p.table) print(basic_summary$s.table) plot(basic_model) Ici, la sortie de ''p.table'' fournit le tableau de résultats pour chaque terme paramétrique et le tableau ''s.table'' nous donne les résultats du terme non linéaire. Notez que pour le second tableau, la courbure du terme non linéaire ''s(X1)'' est indiquée par le paramètre edf (degrés de libertés effectifs); plus la valeur de l'edf est élevée, plus la non-linéarité est forte. Une valeur élevée (8-10 ou plus) signifie que la courbe est fortement non linéaire, alors qu'une courbe avec un edf égal à 1 est une ligne droite. En revanche, dans la régression linéaire, les degrés de libertés du //modèle// sont équivalents au nombre de paramètres libres non redondants p dans le modèle (et les degrés de libertés //résiduels// sont égaux à n-p). Nous reviendrons plus tard sur le concept d'edf. print(basic_summary$p.table) Estimate Std. Error t value Pr(>|t|) (Intercept) 8.550030 0.3655849 23.387258 1.717989e-76 x02 2.418682 0.5165515 4.682364 3.908046e-06 x03 4.486193 0.5156501 8.700072 9.124666e-17 x04 6.528518 0.5204234 12.544629 1.322632e-30 > print(basic_summary$s.table) edf Ref.df F p-value s(x1) 1.923913 2.406719 42.84268 1.076396e-19 Dans notre modèle de base, l'edf du terme non linéaire s(x1) est ~ 2, ce qui indique une courbe non linéaire. Le graphique du modèle illustre bien la forme de ce terme non linéaire : {{::graphic2.2.jpg?425|}} Nous pouvons ajouter un second terme x2, mais spécifier une relation linéaire avec Y (//i.e.// les GAMs peuvent inclure à la fois des termes linéaires et non linéaires dans le même modèle). Ce nouveau terme linéaire x2 sera présenté dans le tableau ''p.table'', pour lequel une estimation du coefficient de régression sera indiquée. Dans le tableau ''s.table'', nous retrouvons encore une fois le terme non linéaire s(x1) et son paramètre de courbure. two_term_model <- gam(y~x0+s(x1)+x2, data=gam_data) two_term_summary <- summary(two_term_model) print(two_term_summary$p.table) print(two_term_summary$s.table) Pour évaluer si la relation entre y et x2 est non linéaire, on peut modéliser x2 avec une fonction non linéaire. Tel que vu auparavant, nous pouvons utiliser une ANOVA pour tester si le terme non linéaire est nécessaire. two_smooth_model <- gam(y~x0+s(x1)+s(x2), data=gam_data) two_smooth_summary <- summary(two_smooth_model) print(two_smooth_summary$p.table) print(two_smooth_summary$s.table) plot(two_smooth_model,page=1) {{::graphic2.4.jpg?600|}} Lorsqu'il y a plus d'une variable d'incluse dans le modèle, comme ci-dessus, la réponse ajustée peut-être partitionnée entre les contributions de chaque variable. Ici, nous pouvons évaluer l'effet de chaque variable où l'axe des ordonnées représente la contribution (effet) de chaque covariable à la réponse ajustée centrée sur 0. Si l'intervalle de confiance chevauche zéro pour certaines valeurs de x, cela indique que l'effet est non significatif. Lorsque la contribution varie selon l'axe x, un changement de cette variable cause un changement de la variable réponse. anova(basic_model,two_term_model,two_smooth_model, test="Chisq") Analysis of Deviance Table Model 1: y ~ x0 + s(x1) Model 2: y ~ x0 + s(x1) + x2 Model 3: y ~ x0 + s(x1) + s(x2) Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 394.08 5231.6 2 393.10 4051.3 0.97695 1180.2 < 2.2e-16 *** 3 385.73 1839.5 7.37288 2211.8 < 2.2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Le meilleur modèle est le modèle avec deux fonctions non linéaires pour x1 et pour x2. ---- **DÉFI 2** Créez deux nouveaux modèles avec la variable x3 : un modèle avec x3 comme paramètre linéaire et un autre modèle avec x3 avec un paramètre non linéaire. Utilisez des graphiques, les tables des coefficients et la fonction ''anova()'' afin de déterminer s'il est nécessaire d'inclure x3 dans le modèle. ++++ Réponse au défi 2| three_term_model <- gam(y~x0+s(x1)+s(x2)+x3, data=gam_data) three_smooth_model <- gam(y~x0+s(x1)+s(x2)+s(x3), data=gam_data) three_smooth_summary <- summary(three_smooth_model) print(three_smooth_summary$p.table) print(three_smooth_summary$s.table) plot(three_smooth_model,page=1) # edf = 1 -> le terme est donc linéaire. anova(two_smooth_model,three_term_model,test="Chisq") # le terme x3 n'est pas significatif ++++ ===== 3. Interactions ===== Il y a deux façons de modéliser une interaction entre deux variables: * si une variable est quantitative et l'autre est qualitative, on utilise l'argument ''by'' -> s(x, by=facteur), * si les deux variables sont quantitatives, on inclut les deux termes sous une même fonction non linéaire -> s(x1, x2). L'argument ''by'' permet de faire varier un terme non linéaire selon les différents niveaux d'un facteur. Nous allons examiner ceci en utilisant notre variable qualitative ''x0'' et examiner si la non-linérité de ''s(x2)'' varie selon les différents niveaux de ''x0''. Pour déterminer si les courbes diffèrent significativement entre les niveaux du facteur, nous allons utiliser une ANOVA sur l'interaction. categorical_interact <- gam(y~x0+s(x1)+s(x2,by=x0),data=gam_data) categorical_interact_summary <- summary(categorical_interact) print(categorical_interact_summary$s.table) plot(categorical_interact,page=1) # ou nous pouvons utiliser la fonction vis.gam où theta représente la rotation du plan x-y vis.gam(categorical_interact,view=c("x2","x0"),theta=40,n.grid=500,border=NA) anova(two_smooth_model, categorical_interact,test="Chisq") {{::graphic3.1b.png?350|}} Nous pouvons constater à partir du graphique que les formes des termes non linéaires sont comparables entre les quatre niveaux de ''x0''. L'ANOVA le confirme également (déviance = 98,6, p = 0,2347). Ensuite, nous allons examiner l'interaction non linéaire entre deux termes quantitatifs, ''x1'' et ''x2''. Cette fois-ci, l'argument ''by'' est supprimé. smooth_interact <- gam(y~x0+s(x1,x2),data=gam_data) smooth_interact_summary <- summary(smooth_interact) print(smooth_interact_summary$s.table) plot(smooth_interact,page=1,scheme=3) # plot(smooth_interact,page=1,scheme=1) donne un graphique comparable à vis.gam() vis.gam(smooth_interact,view=c("x1","x2"),theta=40,n.grid=500,border=NA) anova(two_smooth_model,smooth_interact,test="Chisq") {{::graphic3.2b.png?350|}} L'interaction entre ''s(x1)'' et ''s(x2)'' est significative et le graphique en deux dimensions illustre très bien cette interaction non linéaire. La relation entre y et x1 change en fonction de la valeur de x2. Vous pouvez changez la valeur de l'argument ''theta'' pour tourner l'axe du graphique. Si vous prévoyez exécuter un grand nombre de graphiques, supprimez l'argument ''n.grid = 500'', car ceci fait appel à des calculs intensifs et ralentit R. ===== 4. Changer la fonction de base ===== Sans entrer dans le détail, sachez qu'il est possible de modifier le modèle de base que nous avons vu avec : * des fonctions plus complexes en modifiant la fonction de base (par exemple, cyclique), * d'autres distributions : tout ce que vous pouvez faire avec un GLM (tel que spécifier l'argument ''family'') est possible avec les GAMs, * des modèles à effets mixtes en utilisant la fonction //gamm// ou la fonction //gamm4// de la librairie //gamm4//. Nous allons d'abord jeter un coup d’œil au changement de la fonction de base puis une introduction rapide aux autres distributions et les GAMMs (modèles additifs généralisés à effets mixtes) suivra. Commençons par regarder un cas où modifier la fonction de base peut être utile, soit avec des données cycliques. Imaginez que vous avez une série temporelle de données climatiques, divisées en mesures mensuelles, et que vous voulez déterminer s'il y a une tendance de température annuelle. Nous allons utiliser la série temporelle de température de Nottingham pour cette section : data(nottem) n_years <- length(nottem)/12 nottem_month <- rep(1:12, times=n_years) nottem_year <- rep(1920:(1920+n_years-1),each=12) nottem_plot <- qplot(nottem_month,nottem, colour=factor(nottem_year), geom="line") + theme_bw() print(nottem_plot) En utilisant les données //nottem//, nous avons créé trois nouveaux vecteurs ; ''n_years'' correspond au nombre d'années de données (20 ans), ''nottem_month'' est un codage qualitatif pour les 12 mois de l'année, pour chaque année échantillonnée (série de 1 à 12, répétée 20 fois) et ''nottem_year'' est une variable où l'année correspondant à chaque mois est fournie. Données mensuelles des années 1920 à 1940: {{::graphic4.1.jpg?500|}} Pour modéliser cela, nous devons utiliser ce qu'on appelle un "spline cubique cyclique", ou ''cc'', pour modéliser les effets du mois et de l'année. year_gam <- gam(nottem~s(nottem_year)+s(nottem_month, bs="cc")) summary(year_gam)$s.table plot(year_gam,page=1, scale=0) {{::graphic4.2.jpg?700|}} Il y a une hausse d'environ 1 - 1,5ºC au cours de la série, mais au cours d'une année, il y a une variation d'environ 20ºC. Les données réelles varient autour de ces valeurs prédites et ceci représente donc la variance inexpliquée. Ici, nous pouvons voir l'un des avantages très intéressants de l'utilisation des GAMs. Nous pouvons soit tracer la surface réponse (valeurs prédites) ou les termes (contribution de chaque covariable) tel qu'indiqué ci-haut. Vous pouvez imaginer ce dernier en tant qu'une illustration de la variation des coefficients de régression et comment leur contribution (ou taille de leur effet) varie au fil du temps. Dans le premier graphique, nous voyons que les contributions positives de la température sont survenues après 1930. Sur des échelles de temps plus longues, en utilisant par exemple des données paléolimnologiques, d'autres ([[http://www.aslo.info/lo/toc/vol_54/issue_6_part_2/2529.pdf|Simpson & Anderson 2009;. Fig 3c]]) ont utilisé des GAMs pour tracer la contribution (effet) de la température sur la composition d'algues dans les lacs afin d'illustrer comment les contributions significatives ont seulement eu lieu au cours de deux périodes extrêmement froides (c'est-à-dire, la contribution est importante lorsque les intervalles de confiance ne recoupent pas la valeur de zéro à environ 300 et 100 ans AVJC). Cela a permis aux auteurs de non seulement déterminer combien de variance est expliquée par la température au cours des derniers siècles, mais aussi de repérer dans le temps cet effet significatif. Si cela vous intéresse, le code pour tracer soit la surface de réponse (''type = "response"'') ou les termes (''type = "terms"'') est disponible ci-dessous. Lorsque les termes sont sélectionnés, vous obtiendrez la même figure que celle ci-dessus. ++++ Graphique de contribution vs réponse ajustée| pred<-predict(year_gam, type = "terms", se = TRUE) I<-order(nottem_year) plusCI<-I(pred$fit[,1] + 1.96*pred$se[,1]) minusCI<-I(pred$fit[,1] - 1.96*pred$se[,1]) xx <- c(nottem_year[I],rev(nottem_year[I])) yy <- c(plusCI[I],rev(minusCI[I])) plot(xx,yy,type="n",cex.axis=1.2) polygon(xx,yy,col="light grey",border="light grey") lines(nottem_year[I], pred$fit[,1][I],lty=1) abline(h=0) ++++ ===== 5. Intro rapide aux GAMMs ===== ==== La non-indépendance des données ==== Lorsque les observations ne sont pas indépendantes, les GAMs peuvent être utilisés soit pour incorporer: * une structure de corrélation pour modéliser les résidus autocorrélés (autorégressif (AR), moyenne mobile (MA), ou une combinaison des deux (ARMA)) * des effets aléatoires qui modélisent l'indépendance entre les observations d'un même site. En plus de changer la fonction de base, nous pouvons aussi complexifier le modèle en intégrant une structure d'auto-corrélation (ou même des effets mixtes) en utilisant les fonctions //gamm// ou //gamm4//. Pour commencer, nous allons jeter un coup d’œil au premier cas ; un modèle avec autocorrélation temporelle dans les résidus. Ré-examinons le modèle de la température de Nottingham; nous allons vérifier si les résidus sont corrélés en faisant appel à la fonction (partielle) d'autocorrélation. par(mfrow=c(1,2)) acf(resid(year_gam), lag.max = 36, main = "ACF") pacf(resid(year_gam), lag.max = 36, main = "pACF") {{::graphic5.1.jpg?700|}} Les graphiques des fonctions d'autocorrélation suggèrent qu'un modèle AR de faible ordre est nécessaire (avec un ou deux intervalles de temps décalés), donc nous pouvons évaluer deux modèles; ajouter un AR(1) ou un AR(2) au modèle de la température de Nottingham et évaluer le meilleur avec une ANOVA. year_gam <- gamm(nottem~s(nottem_year)+s(nottem_month, bs="cc")) year_gam_AR1 <- gamm(nottem~s(nottem_year)+s(nottem_month, bs="cc"), correlation = corARMA(form = ~ 1|nottem_year, p = 1)) year_gam_AR2 <- gamm(nottem~s(nottem_year)+s(nottem_month, bs="cc"), correlation = corARMA(form = ~ 1|nottem_year, p = 2)) anova(year_gam$lme,year_gam_AR1$lme,year_gam_AR2$lme) Model df AIC BIC logLik Test L.Ratio p-value year_gam$lme 1 5 1109.908 1127.311 -549.9538 year_gam_AR1$lme 2 6 1101.218 1122.102 -544.6092 1 vs 2 10.689206 0.0011 year_gam_AR2$lme 3 7 1101.598 1125.962 -543.7988 2 vs 3 1.620821 0.2030 Le modèle avec la structure AR(1) prévoit une augmentation significative comparativement au premier modèle (LRT = 10,69, p = 0,0011), mais il y a très peu d'intérêt à considérer le modèle AR(2) (LRT = 1,62, b = 0,203). ==== Modélisation avec effets mixtes ==== Comme nous l'avons vu dans la section précédente, ''bs'' spécifie la fonction de base sous-jacente. Pour les facteurs aléatoires (origine et pente linéaire), nous utilisons ''bs = "re"'' et pour les pentes aléatoires non linéaires, nous utilisons ''bs = "fs"''. Trois types d'effets aléatoires différents sont possibles lors de l'utilisation des GAMMs (où //fac// représente une variable qualitative utilisée pou l'effet aléatoire et //x0// est un effet quantitatif fixe) : * //interceptes aléatoires// ajustent la hauteur des termes du modèle avec une valeur constante de pente : s(fac, bs="re") * //pentes aléatoires// ajustent la pente d'une variable explicative numérique: s(fac, x0, bs="re") * //surfaces lisses aléatoires// ajustent la tendance d'une prédiction numérique de façon non linéaire: s(x0, fac, bs="fs", m=1) où l'argument m = 1 met une plus grande pénalité au lissage qui s'éloigne de 0, ce qui entraîne un retrait vers la moyenne. Nous examinerons d'abord un GAMM avec un interception aléatoire. Tel que vu précédemment, nous allons utiliser ''gamSim()'' pour générer un ensemble de données, cette fois-ci avec une composante d'effet aléatoire. Ensuite, nous construirons un modèle avec un intercepte aléatoire en utilisant //fac// comme facteur aléatoire. # Générez des données gam_data2 <- gamSim(eg=6) head(gam_data2) # Faites rouler un modèle avec intercepte aléatoire gamm_intercept <- gam(y ~ s(x0) + s(fac, bs="re"), data=gam_data2) summary(gamm_intercept) Notez le terme aléatoire dans le tableau. Vous pouvez le visualiser: plot(gamm_intercept, select=2) # select=2 parce que le terme aléatoire se trouve sur la 2e ligne du tableau sommaire. Une fonction de traçage vraiment intéressante que nous allons maintenant utiliser est le ''plot_smooth'' de la librairie //itsadug//. Contrairement au graphique par défaut ''plot.gam'', cette fonction présente l'effet additionné du GAMM avec l'option de ne pas inclure les courbes aléatoires dans le graphique. Ici, nous allons premièrement tracer l'effet combiné de x0 (sans les niveaux de l'effet aléatoire) et ensuite une courbe pour les quatre niveaux de //fac//: par(mfrow=c(1,2), cex=1.1) plot_smooth(gamm_intercept, view="x0", rm.ranef=TRUE, main="intercept + s(x1)", rug=FALSE) plot_smooth(gamm_intercept, view="x0", cond=list(fac="1"), main="... + s(fac)", col='orange', ylim=c(8,21), rug=FALSE) plot_smooth(gamm_intercept, view="x0", cond=list(fac="2"), add=TRUE, col='red') plot_smooth(gamm_intercept, view="x0", cond=list(fac="3"), add=TRUE, col='purple') plot_smooth(gamm_intercept, view="x0", cond=list(fac="4"), add=TRUE, col='turquoise') {{::graphic5.3.jpg?650|}} Ensuite, nous allons générer et tracer un modèle avec une pente aléatoire : gamm_slope <- gam(y ~ s(x0) + s(x0, fac, bs="re"), data=gam_data2) summary(gamm_slope) plot_smooth(gamm_slope, view="x0", rm.ranef=TRUE, main="intercept + s(x0)", rug=FALSE) plot_smooth(gamm_slope, view="x0", cond=list(fac="1"), main="... + s(fac)", col='orange',ylim=c(7,22), rug=FALSE) plot_smooth(gamm_slope, view="x0", cond=list(fac="2"), add=TRUE, col='red') plot_smooth(gamm_slope, view="x0", cond=list(fac="3"), add=TRUE, col='purple') plot_smooth(gamm_slope, view="x0", cond=list(fac="4"), add=TRUE, col='turquoise') {{::graphic5.4.jpg?650|}} Nous allons maintenant inclure à la fois un intercepte et une pente aléatoires. gamm_int_slope <- gam(y ~ s(x0) + s(fac, bs="re") + s(fac, x0, bs="re"), data=gam_data2) summary(gamm_int_slope) plot_smooth(gamm_int_slope, view="x0", rm.ranef=TRUE, main="intercept + s(x0)", rug=FALSE) plot_smooth(gamm_int_slope, view="x0", cond=list(fac="1"), main="... + s(fac) + s(fac, x0)", col='orange', ylim=c(7,22), rug=FALSE) plot_smooth(gamm_int_slope, view="x0", cond=list(fac="2"), add=TRUE, col='red', xpd=TRUE) plot_smooth(gamm_int_slope, view="x0", cond=list(fac="3"), add=TRUE, col='purple', xpd=TRUE) plot_smooth(gamm_int_slope, view="x0", cond=list(fac="4"), add=TRUE, col='turquoise', xpd=TRUE) {{::graphic5.5.jpg?650|}} Notez que les pentes aléatoires sont statique : plot(gamm_int_slope, select=3) # select=3 parce que la pente aléatoire se trouve sur la 3e ligne du tableau sommaire. Enfin, nous allons examiner un modèle avec une surface lisse aléatoire. gamm_smooth <- gam(y ~ s(x0) + s(x0, fac, bs="fs", m=1), data=gam_data2) summary(gamm_smooth) Ici, si les pentes aléatoires variaient selon x0, nous auront vue des courbe variable pour chaque niveau : plot(gamm_smooth, select=1) # select=1 parce que le terme se trouve sur la 1e ligne du tableau sommaire. Finalement, tous ces modèles mixes peuvent être compare en utilisant la fonction anova() pour trouver le meilleur modèle. ===== 6. Autres distributions ===== Pour vous donner un bref aperçu de l'utilisation des GAMs lorsque la variable réponse ne suit pas une distribution normale ou que les données sont des abondances ou proportions (par exemple, distribution Gamma, binomiale, Poisson, binomiale négative), l'exemple qui suit utilise un ensemble de données où une répartition binomiale sera nécessaire, y compris une modélisation d'une relation non linéaire. La variable réponse représente le nombre de succès (l'événement a eu lieu) en fonction des défaillances au cours d'une expérience. gam_data3 <- read.csv("other_dist.csv") summary(gam_data3) str(gam_data3) 'data.frame': 514 obs. of 4 variables: $ prop : num 1 1 1 1 0 1 1 1 1 1 ... $ total: int 4 20 20 18 18 18 20 20 20 20 ... $ x1 : int 550 650 750 850 950 650 750 850 950 550 ... $ fac : Factor w/ 4 levels "f1","f2","f3",..: 1 1 1 1 1 1 1 1 1 1 ... ''prop'' est la variable réponse, égal à la proportion de //succès / (succès + échecs)//. Notez qu'il existe de nombreux cas où la proportion est égal à 1 ou 0 qui indique que les résultats ont toujours été des succès ou des échecs, respectivement, à ce moment mesuré durant l'expérience.\\ ''x1'' est le temps écoulé depuis le début de l'expérience (variable explicative).\\ ''total'' représente le nombre de //succès + échecs// observé au moment x1i de l'expérience.\\ ''fac'' est un facteur qui code pour l'essai 1 à 4 de l'expérience (nous n'utiliserons pas cette variable dans cette section). Commençons par la visualisation des données. Nous sommes intéressés par le nombre de succès par rapport aux échecs à mesure que x1 augmente. Étant donné qu'il y a des mesures répétées pour la valeur de x1 (essais 1 à 4, avec nombreuses observations par essai), nous pouvons d'abord présenter la proportion de succès en moyenne par boîte de temps (x1): emptyPlot(range(gam_data3$x1), c(0,1), h=.5, main="Probability of successes", ylab="Probability",xlab="x1") avg <- aggregate(prop ~ x1, data=gam_data3, mean, na.rm=TRUE) lines(avg$x1, avg$prop, col="orange",lwd=2) {{::graphic1_otherdist.jpg?450|}} Notez comment la probabilité de succès augmente avec x1. D'après vous, est-ce que cette tendance est linéaire ou non linéaire? Nous allons tester cela en utilisant un GAM logistique (nous utilisons une distribution ''binomiale'' puisque la variable réponse représente des proportions). prop_model <- gam(prop~ s(x1), data=gam_data3, weights=total, family="binomial") prop_summary <- summary(prop_model) print(prop_summary$p.table) print(prop_summary$s.table) plot(prop_model) Estimate Std. Error z value Pr(>|z|) (Intercept) 1.173978 0.02709613 43.32641 0 edf Ref.df Chi.sq p-value s(x1) 4.591542 5.615235 798.9407 2.027751e-164 Qu'est-ce que l'ordonnée représente dans ce modèle? * Rappel : le modèle utilise le nombre de succès vs échecs pour calculer le //logit//, qui est le logarithme du rapport entre les succès et échecs: logit = log({N_{success}/N_{failures}}) - Si succès = échecs, le rapport est de 1 et le logit est 0 (log (1) = 0).\\ - Si les succès ont un nombre plus grand que les échecs, le ratio est supérieur à 1 et le logit a une valeur positive (par exemple, log(2) = 0,69).\\ - Si les succès ont un nombre plus petit que les échecs, le ratio est inférieur à 1 et le logit a une valeur négative (par exemple, log(0,5) = -0.69).\\ Donc, l'ordonnée est le //logit//, et indique s'il y a en moyenne plus de succès que d'échecs. Ici, l'estimé est positif ce qui signifie, qu'en moyenne, il n'y a plus de succès que d'échecs. Qu'est-ce que le terme de lissage indique? * Ceci représente la façon dont les chances de succès vs échecs changent sur l'échelle de x1 (l'échelle du temps dans cet exemple). Donc, puisque l'edf > 1, la proportion de succès augmente plus rapidement au fil du temps (si par exemple, la réponse représente le nombre d'individus de l'espèce A vs l'espèce B et que nous augmentons la concentration des nutriments au fil du temps, ces résultats indiqueront que l'espèce A est de plus en plus observée alors que les concentrations de nutriments approchent de l'optimum de cette espèce au cours de l'expérience). === Visualiser la tendance au fil du temps === Enfin, nous allons voir les différentes façons de représenter ces relations graphiquement. par(mfrow=c(1,2)) plot(prop_model, select=1, scale=0, shade=TRUE) abline(h=0) plot_smooth(prop_model, view="x1",main="") (diff <- find_difference(out$fv$fit, out$fv$CI, xVals=out$fv$x1)) addInterval(0, lowVals=diff$start, highVals = diff$end, col='red', lwd=2) abline(v=c(diff$start, diff$end), lty=3, col='red') text(mean(c(diff$start, diff$end)), 2.1, "sign. more \n success", col='red', font=3) {{::graphic2_otherdist.jpg?800|}} Quels renseignements ces graphiques nous apportent-ils vis à vis les succès et échecs ? * **graphique de gauche**: contribution (ou effet partiel si nous avions plus qu'une variable explicative) au fil du temps. La valeur logit augmente, donc les succès augmentent et les échecs diminuent. * **graphique de droite**: valeurs ajustées, ordonnée incluse (somme des effets si nous avions plus d'une variable explicative dans le modèle). Nous voyons ici que la valeur logit est estimée près de zéro au début de l'expérience ; cela signifie qu'il y a des quantités égales de succès et d'échecs. Peu à peu, les succès augmentent et à environ x1=400 il y a beaucoup plus de succès que d'échecs (l'effet est significativement différent de zéro). Nous avons également montré comment nous pouvons utiliser le graphique pour déterminer à quelle valeur de x1 cela se produit. Enfin, pour nous aider à interpréter les résultats, nous pouvons re-transformer l'effet sur une échelle de proportions avec la fonction ''plot_smooth'' de la librairie //itsadug//: par(mfrow=c(1,1)) plot_smooth(prop_model, view="x1", main="", transform=plogis, ylim=c(0,1)) abline(h=.5, v=diff$start, col='red', lty=2) {{::graphic3_otherdist.jpg?450|}} Comme nous l'avons déjà vu avec le graphique précédent des valeurs logits, nous voyons qu'à approximativement x1=400 la proportion de succès augmente de façon significative au-dessus de 0,5. ===== 7. Les GAMs en coulisse ====== Nous allons maintenant prendre quelques minutes pour regarder comment fonctionnent les GAMs. Commençons en considérant d'abord un modèle qui contient une fonction lisse d'une covariable, x_i : {y_i} = f(x_{i}) + {ε_i} Pour estimer la fonction //f//, nous avons besoin de représenter l'équation ci-dessus de manière à ce qu'elle devienne un modèle linéaire. Cela peut être fait en choisissant une base, b_i(x), définissant l'espace des fonctions dont //f// est un élément: f(x) = sum{i=1}{q}{b_{i}(x)β_{i}} **Un exemple simple : une base polynomiale** Supposons que //f// est considéré comme un polynôme d'ordre 4, de sorte que l'espace des polynômes d'ordre 4 et moins contiennent //f//. Une base de cet espace serait alors : b_{1}(x)=1, b_{2}(x)=x, b_{3}(x)=x^{2}, b_{4}(x)=x^{3}, b_{5}(x)=x^{4} et f(x) devient: f(x) = β_{1} + x_{i}β_{2} + x^{2}_{i}β_{3} + x^{3}_{i}β_{4}(x) + x^{4}_{i}β_{5} et le modèle complet devient: y_{i} = β_{1} + x_{i}β_{2} + x^{2}_{i}β_{3} + x^{3}_{i}β_{4}(x) + x^{4}_{i}β_{5} + ε_{i} Chaque fonction de base est multipliée par un paramètre à valeur réelle, β_i, et est ensuite additionnée pour donner la courbe finale **f(x)**. {{::polynomial_basis_example.png?600|}} En faisant varier le coefficient β_i, on peut faire varier la forme de f(x) pour produire une fonction polynomiale d'ordre 4 ou moins. **Un autre exemple : une base de spline cubique** Un spline cubique est une courbe construite à partir de sections d'un polynôme cubique reliées entre elles de sorte qu'elles sont continues en valeur. Chaque section du spline a des coefficients différents. {{::cubic_spline_fr.png?550|}} Voici une représentation d'une fonction lisse utilisant une base de régression spline cubique de rang 5 avec des nœuds situés à incréments de 0,2: {{::graphic6.1.jpg?300|}} Dans cet exemple, les nœuds sont espacés uniformément à travers la gamme des valeurs observées de x. Le choix du degré de finesse du modèle est pré-déterminé par le nombre de noeuds, qui était arbitraire. Y a-t-il une meilleure façon de sélectionner les emplacements des nœuds? **Contrôler le degré de lissage avec des splines de régression pénalisés** Au lieu de contrôler le lissage (non linéarité) en modifiant le nombre de nœuds, nous gardons celle-ci fixée à une taille un peu plus grande que raisonnablement nécessaire et on contrôle le lissage du modèle en ajoutant une pénalité sur le niveau de courbure. Donc, plutôt que d'ajuster le modèle en minimisant (comme avec la méthode des moindres carrés) : ||y - XB||^{2} il peut être modélisé en minimisant : ||y - XB||^{2} + {lambda}int{0}{1}{[f^{''}(x)]^{2}dx} Quant lambda tend vers infty , le modèle devient linéaire. Pour la sélection du meilleur paramètre de lissage, lambda, on utilise une approche de validation croisée. Si lambda est trop élevé, les données seront trop lissées et si elle est trop faible, les données ne seront pas assez lissées. Idéalement, il serait bon de choisir une valeur lambda de sorte que le //f prédit// est aussi proche que possible du //f observé//. Un critère approprié pourrait être de choisir lambda pour minimiser : M = 1/n sum{i=1}{n}{(hat{f_{i}} - f_{i})^{2}} Étant donné que //f// est inconnu, M est estimé en utilisant une technique de validation croisée généralisée qui laisse de côté, à chaque tour, une donnée et estime la capacité moyenne des modèles, construits sur les données restantes, de prédire la donnée qui a été mise de côté (pour plus de détails, consultez Wood 2006). **Illustration du principe de validation croisée** {{::illustration_of_smooth_sel.png?600|}} Dans le premier panneau, la courbe correspond à un ajustement faible par rapport aux données et ne fait pas mieux avec le point manquant. Dans le troisième panneau, la courbe ajuste le bruit aussi bien que le signal et la variabilité supplémentaire induite l'amène à prédire la donnée manquante plutôt mal. Dans le deuxième panneau, cependant, nous voyons que l'ajustement de la courbe du signal sous-jacent est très bien, le lissage passe à travers le bruit et la donnée manquante est raisonnablement bien prédite. {{::gcv.png?600|}} **Note supplémentaire sur les degrés de liberté effectifs (edf)** Combien de degrés de liberté a un GAM ? Au lieu de fournir la sortie de la validation croisée en termes de lambda (un paramètre qui détermine la complexité de l'ajustement), la fonction ''gam()'' utilise un terme appelé les degrés de liberté effectifs (edf), de manière cohérente à quantifier la complexité du modèle (pour justifier notre intuition que les degrés de liberté offrent une manière cohérente de quantifier la complexité du modèle). Parce que le nombre de paramètres libres des splines de lissage (tel que les GAMs) est souvent difficile à définir, les edf sont liés à lambda, où l'effet de la pénalité est de réduire les degrés de libertés. Donc, si le nombre de noeuds est arbitrairement réglé à k = 10, k-1 définit la limite supérieure des degrés de libertés associés à un terme de lissage. Ce nombre diminue alors que la pénalité lambda augmente jusqu'à ce que le meilleur modèle soit trouvé par validation croisée. ===== Références ===== Il existe beaucoup d'information sur les GAMs... Simon Wood, l'auteur de la librairie //mgcv//, a un site très utile avec des conférences et des notes introductives sur la façon d'utiliser les GAMs : * http://people.bath.ac.uk/sw283/mgcv/ Il a aussi écrit un livre, //Generalized Additive Models: An Introduction with R//, que nous avons utilisé comme référence pour cet atelier. Le matériel de cet atelier a également été obtenu à partir des blogs et des tutoriels suivants : * http://www.fromthebottomoftheheap.net/blog/ * http://www.sfs.uni-tuebingen.de/~jvanrij/Tutorial/GAMM.html * http://www.sfs.uni-tuebingen.de/~jvanrij/LSA2015/AnswersLab2.html Enfin, les pages d'aide, disponibles via ''?gam'' dans R sont une excellente ressource.