CSBQ Atelier R 8

0. Le modèle linéaire ... et où il échoue

1. Introduction aux GAMs

2. Plusieurs termes non linéaires

3. Interactions

4. Changer la fonction de base

5. Intro rapide aux GAMMs

6. Autres distributions

7. Les GAMs en coulisse

Aperçu

Sujets

Objectifs de l'atelier

Prérequis

• Expérience en R; assez pour être en mesure d'exécuter un script et d'examiner les données et objets R
• 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

• Que veut on dire par un modèle non-linéaire?
• Qu'est-ce qu'un modèle additif généralisé (GAM), et comment nous permet-il de traiter des relations non-linéaires?
• Comment interpréter et tracer ces relations non-linéaires
• Comment ajouter des termes d'interaction

• Utilisez la librairie mgcv pour modéliser les relations non-linéaires
• Évaluer la sortie d'un GAM afin de mieux comprendre nos données
• Utilisez 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 pourquoi ceci rend les GAMs si puissant!
• Comment modéliser la dépendance dans les données (autocorrélation, structure hiérarchique) en utilisant GAMMs

Qu'est-ce qu'un modèle linéaire?

et pourtant tant de façons pour qu'il ne le soit pas:

Alors, comment résoudre ce problème?

Il n'y a qu'une façon qu'un modèle linéaire soit correctement appliqué:

• La régression linéaire est ce que la plupart des gens apprennent avant tout en statistiques et parmi les approches les plus performantes
• Elle nous permet de modéliser une variable réponse en fonction de facteurs prédictifs, plus de l'erreur
• Elle fait cependant quatre suppositions importantes:

Pour résoudre le problème de non-linéarité, nous devons 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, en passant tout simplement une ligne entre chaque point.

Les modèles linéaires le font en trouvant la meilleure ligne de droite ajuster qui passe à travers les données.

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.

Tester la linéarité vs non-linéarité

• l'erreur est normalement distribuée
• la variance des erreurs ne change pas
• chaque erreur est indépendante des autres
• la réponse est linéaire: y = a*x + b

Voici comment nous le ferons en utilisant un terme non-linéaire avec le modèle additif généralisé (section 1.3):

Solution

Défi 1

mgcv comprend également une fonction plot qui, par défaut, nous permet de visualiser la non-linéarité du modèle (section 1.4):

Examinons un exemple. Premièrement, nous allons générer des données, et les tracer. Exécutez la section 1.1 dans le code:

Voici comment nous ajusterions le modèle en utilisant une régression linéaire (section 1.2):

> plot(gam_model)

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 (section 1.5):

Nous allons maintenait essayer cela avec d'autres données composées. Nous allons d'abord générer les données (section 1.6):

> 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)

> 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"))

> 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)

Family: gaussian

Link function: identity

Formula:

y_obs ~ x

Parametric coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 0.765606 0.021181 36.15 <2e-16 ***

x 0.157892 0.007515 21.01 <2e-16 ***

---

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) = 0.639 Deviance explained = 64%

GCV score = 0.026765 Scale est. = 0.026551 n = 250

> 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)

> linear_model = gam(y_obs~x)

> nested_gam_model = gam(y_obs~s(x)+x)

> print(anova(linear_model, nested_gam_model, test="Chisq"))

Family: gaussian

Link function: identity

Formula:

y_obs ~ s(x)

Parametric coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) 1.154422 0.006444 179.1 <2e-16 ***

---

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:

edf Ref.df F p-value

s(x) 8.317 8.872 170.9 <2e-16 ***

---

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) = 0.859 Deviance explained = 86.3%

GCV score = 0.010784 Scale est. = 0.010382 n = 250

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 ***

> 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)

Pensez à cela comme un graphique de l'effet (ou contribution) du terme de lissage.

Ou, dit autrement, comment le coefficient varie sur la gamme de xi.

Notez que lorsque nous utilisons un modèle additif, la relation entre les prédicteurs individuels et la variable dépendante est estimée par une fonction non-linéaire (un "smoother"):

g(y) = s(x1) + s(x2, x3) + ...

> summary(nested_gam_model_test)$s.table

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

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 la linéarité est justifiée pour ces données.

edf Ref.df F p-value

s(x_test) 7.602145 8.029057 294.0944 0

Comme nous allons le voir dans la diapositive suivante, l'avantage d'un GAM est que la forme optimale de la non-linéarité est déterminée automatiquement

• le degré de courbure est déterminé en utilisant une validation croisée pour éviter un sur-ajustage

Oui la non-linéarité est justifiée. Les degrés de libertés effectifs (edf) sont » 1. (nous reviendrons sur ce point plus tard...)

Centre de la science de la biodiversité du Québec

Série d'ateliers R

Ou, encore une fois, en utilisant vis.gam():

> plot(smooth_interact,page=1)

Nous pouvons également utiliser la fonction de vis.gam pour créer un graphique 2D

L'argument »by« vous permet de varier un terme non-linéaire entre les différents niveaux d'un facteur (section 3.1):

Nous allons maintenant examiner l'interaction non linéaire entre deux termes quantitatifs, x1 et x2 (section 3.2):

> plot(categorical_interact,page=1)

> anova(two_smooth_model,smooth_interact,test="Chisq")

Comment modéliser des interactions avec les GAMs?

Nous allons utiliser ANOVA pour déterminer si les courbes diffèrent significativement entre les niveaux du facteur.

> vis.gam(smooth_interact,view=c("x1","x2"),theta=40, n.grid=500, border=NA)

> vis.gam(categorical_interact,view=c("x2","x0"),theta=40,n.grid=500, border=NA)

> smooth_interact = gam(y~x0+s(x1,x2),data= gam_data)

> smooth_interact_summary = summary(smooth_interact)

> print(smooth_interact_summary$s.table)

> anova(two_smooth_model, categorical_interact,test="Chisq")

> 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)

Il y a deux façons de modéliser une interaction entre deux variables:

• si une variable est quantitative, et l'autre est qualitative: en utilisant l'argument »by« s(x, by=facteur)
• si les deux variables sont quantitatives: en mettant les deux termes sous une fonction non linéaire s(x1,x2).

Commençons par un modèle de base, avec un terme non-linéaire (x1) et facteur qualitatif (X0, avec 4 niveaux) (section 2.2).

Analysis of Deviance Table

Model 1: y ~ x0 + s(x1) + s(x2)

Model 2: y ~ x0 + s(x1, x2)

Resid. Df Resid.Dev Df Deviance Pr(>Chi)

1 385.73 1839.5

2 370.08 1710.3 15.644 129.2 0.02804 *

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> plot(two_term_model)

Plusieurs variables

Quelques mots sur l'edf ...

Modélisez x2 avec une fonction non linéaire (section 2.4):

Maintenant, ajoutez une seconde variable non linéaire (section 2.3):

Défi 2

edf Ref.df F p-value

s(x1,x2) 25.91803 28.42892 36.39357 9.315854e-165

Le graphique du modèle illustre bien la forme de ce terme non linéaire:

> plot(two_smooth_model,page=1)

Comme avant, nous pouvons utiliser une ANOVA pour tester si le terme non-linéaire est nécessaire (section 2.5):

Solution

Analysis of Deviance Table

Model 1: y ~ x0 + s(x1) + s(x2)

Model 2: y ~ x0 + s(x1) + s(x2, by = x0)

Resid. Df Resid. Dev Df Deviance Pr(>Chi)

1 385.73 1839.5

2 368.67 1740.9 17.06 98.608 0.2347

Solution

> basic_model = gam(y~x0+s(x1), data= gam_data)

Avec les GAMs, il est facile d'ajouter des termes non-linéaires et linéaires dans un seule modèle, ou plusieurs termes non-linéaires, ou même des interactions non-linéaires. Pour cette section, nous allons utiliser un ensemble de données généré automatiquement par mgcv (section 2.1):

> 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)

> plot(basic_model)

> 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)

edf Ref.df F p-value

s(x1) 2.401350 2.996664 67.18567 6.325091e-38

s(x2):x01 6.606775 7.708292 21.40363 2.572371e-27

s(x2):x02 6.436261 7.571201 19.96801 5.469274e-25

s(x2):x03 5.467172 6.618716 28.73636 2.938870e-32

s(x2):x04 6.422564 7.574388 26.43148 2.092391e-33

A partir de là, nous pouvons regarder les tableaux de coefficients de ce modèle:

> anova(basic_model,two_term_model,two_smooth_model, test="Chisq")

> gam_data = gamSim(eg=5)

> basic_summary = summary(basic_model)

> print(basic_summary$p.table)

> anova(two_smooth_model,three_term_model,test="Chisq")

Nous reviendrons sur l'interprétation de l'edf plus tard dans l'atelier, mais pour le moment garder en tête que le plus que l'edf est élevé, plus que 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 égale à 1 est une ligne droite.

En revanche, dans la régression linéaire, les degrés de libertés du modèle est équivalent au nombre de paramètres libres non redondants, p, dans le modèle (et les degrés de libertés résiduels sont égales à n-p)

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.

Estimate Std. Error t value Pr(>|t|)

(Intercept) 8.937862 0.2217506 40.305927 2.373755e-140

x02 2.008045 0.3137690 6.399756 4.518133e-10

x03 3.832496 0.3143049 12.193562 3.758930e-29

x04 6.041521 0.3145299 19.208098 3.520507e-58

Additive model + factor

Estimate Std. Error t value Pr(>|t|)

(Intercept) 11.400658 0.4177614 27.289879 9.396207e-93

x02 2.314405 0.4552138 5.084216 5.723467e-07

x03 4.487653 0.4543299 9.877520 1.063008e-20

x04 6.596149 0.4585778 14.383925 5.468771e-38

x2 -5.825948 0.5436671 -10.716021 1.114046e-23

Visuellement, nous avons vu que la forme de termes lisses étaient largement comparables entre les quatre niveaux de x0. Le test anova confirme ceci (Déviance = 98,6, p = 0,2347).

> 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$s.table)

> plot(three_smooth_model,page=1)

> head(gam_data)

Créer deux nouveaux modèles, avec x3 comme un terme linéaire et non linéaire.

Utilisez les graphiques, les tableaux de coefficients et la fonction anova() pour déterminer si nous devons ajouter x3 a notre modèle.

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

*plot(smooth_interact,page=1,scheme=1) donne un graphique comparable

Model 1: y ~ x0 + s(x1) + s(x2)

Model 2: y ~ x0 + s(x1) + s(x2) + x3

Resid. Df Resid. Dev Df Deviance Pr(>Chi)

1 385.73 1839.5

2 384.73 1839.3 0.99997 0.18818 0.8427

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

> print(two_smooth_summary$s.table)

> print(two_term_summary$s.table)

y x0 x1 x2 x3

1 4.723147 1 0.02573032 0.70706571 0.69248543

2 8.886671 2 0.83272144 0.84997218 0.88974095

3 11.196905 3 0.66302652 0.88025265 0.08469529

4 10.886068 4 0.11126873 0.80087554 0.15109792

5 12.270534 1 0.87969756 0.37692184 0.51467778

6 9.020910 2 0.12441532 0.05154493 0.86526950

> print(basic_summary$s.table)

edf Ref.df F p-value

s(x1) 2.542296 3.170441 67.75922314 1.573017e-39

s(x2) 7.731424 8.584355 80.90188472 1.745178e-119

s(x3) 1.000000 1.000000 0.02039697 8.865087e-01

edf Ref.df F p-value

s(x1) 1.900864 2.377544 49.85908 1.938884e-22

edf Ref.df F p-value

s(x1) 1.923913 2.406719 42.84268 1.076396e-19

edf Ref.df F p-value

s(x1) 2.546757 3.175726 68.08822 9.272491e-40

s(x2) 7.726989 8.582003 81.42428 2.845063e-120

Le terme x3 n'est pas significatif. Nous devons le supprimer du modele.

Nous allons voir comment nous pouvons prédire la variable réponse, y, en fonction des autres variables.

La courbure du terme non-linéaire, s(X1), est indiquée par le paramètre edf (degrés de libertés effectifs).

Le edf est plus grand que 1, donc le terme est donc linéaire. Evaluons maintenant la relation linéaire entre y et x3; est-elle significative?

Atelier 8: Modèles additifs généralisés

Modèle avec erreurs autocorréles

Modélisation mixte

Les graphiques ACF et pACF suggèrent qu'un modèle AR de faible ordre est nécessaire (avec une ou deux intervalles de temps décalés), donc nous pouvons évaluer deux modèles; AR(1) et AR(2) (section 5.2)

Comment faire face à la non-indépendance

Ici, si les pentes aléatoires variaient selon x0, nous auront vue des courbe variable pour chaque niveau:

> gamm_intercept <- gam(y ~ s(x0) + s(fac, bs="re"), data=gam_data2)

> summary(gamm_intercept)$s.table

Tracer les effets combinés de x0 (sans les niveaux de l'effet aléatoire), et ensuite une courbe pour tous les quatre niveaux de fac:

Notez que les pentes aléatoires sont statiques :

Enfin, nous allons examiner un modèle avec une surface lisse aléatoire (section 5.6).

> lme.ex <- lme(y ~ 1,data=gam_data2,~1|fac,method="REML")

> gamm.ex <- gam(y ~ s(fac,bs="re"),data=gam_data2,method="REML")

Ensuite, nous allons générer et tracer un modèle avec une pente aléatoire (section 5.4).

> plot(year_gam,page=1, scale=0)

Nous examinerons d'abord un GAMM avec un interception aléatoire.

Comme avant, nous allons utiliser gamSim() pour générer automatiquement un ensemble de données, cette fois avec une composante d'effet aléatoire, ensuite rouler un modèle avec un origine aléatoire en utilisant fac comme facteur aléatoire (section 5.3)

Pour commencer, regardons les résidus du modèle de température Nottingham en utilisant la fonction d'autocorrélation (partielle) (section 5.1)

Nous allons maintenant inclure à la fois un origine de droite et une pente aléatoire (section 5.5).

Tracer les effets combinés de x0 (sans les niveaux de l'effet aléatoire), et ensuite une courbe pour tous les quatre niveaux de fac:

Paramètres des fonctions lisses

Données mensuelles des années 1920 à 1940

Terme de lissage

Comment changer la fonction de base?

plot(gamm_smooth, select=2)

fac 4

Tracer les effets combinés de x0 (sans les niveaux de l'effet aléatoire), et ensuite une courbe pour tous les quatre niveaux de fac:

plot(gamm_int_slope, select=3)

Pour modéliser cela, nous devons utiliser ce qu'on appelle une spline cubique cyclique, ou «cc», pour modéliser les effets de mois ainsi qu'un terme pour l'année (section 4.2):

edf Ref.df F p-value

s(x0) 3.044725 3.776643 2.811258 2.835663e-02

s(fac) 2.960269 3.000000 95.146996 2.696335e-54

charac. saisonnier

(spline cyclique cubique)

tendance

> gamm_smooth <- gam(y ~ s(x0)

+ s(x0, fac, bs="fs", m=1), data=gam_data2)

> summary(gamm_smooth)

Les fonctions lisses ont beaucoup de paramètres qui pourraient changer leur comportement. Les paramètres les plus souvent utilisés sont les suivants :

intervals(lme.ex)

> gamm_slope <- gam(y ~ s(x0) + s(x0, fac, bs="re"), data=gam_data2)

> summary(gamm_slope)

fac 3

> par(mfrow=c(1,2))

> acf(resid(year_gam), lag.max = 36, main = "ACF")

> pacf(resid(year_gam), lag.max = 36, main = "pACF")

Pour modéliser une surface lisse ou non linéaire, nous pouvons utiliser trois fonctions lisses différentes :

> gam_data2 <- gamSim(eg=6)

> gamm_int_slope <- gam(y ~ s(x0) + s(fac, bs="re")

+ s(fac, x0, bs="re"), data=gam_data2)

> summary(gamm_int_slope)

Regardons un cas où changer la fonction de base peut être utile:

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 si il y a une tendance de température annuelle.

Nous allons utiliser la série temporelle de température de Nottingham pour cette section (section 4.1):

> 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)

fac 2

4 term additive + random effect

> year_gam = gam(nottem~s(nottem_year) + s(nottem_month, bs = "cc"))

> summary(year_gam)$s.table

Random Effects:

Level: fac

lower est. upper

sd((Intercept)) 1.903369 4.267269 9.567026

Et faire un beau graphique:

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))
• effets aléatoires qui modélise l'indépendance entre les observations d'un même site

C'est à dire, en plus de changer la fonction de base, nous pouvons aussi ajouter de la complexité du modèle en intégrant une structure d'auto-corrélation ou des effets mixtes en utilisant les librairies gamm ou gamm4.

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 se caractérisent lors de l'utilisation GAMMs:

• origines aléatoires ajuste l'hauteur des termes du modèle avec une valeur constante: s(fac, bs="re")
• pentes aléatoires ajuste la pente d'une variable explicative numérique: s(fac, x0, bs="re")
• smooths aléatoires ajuste 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 l'éloignement de 0, ce qui entraîne un retrait vers la moyenne.

> head(gam_data)

fac 1

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)

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')

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')

Nous n'allons pas rentrer trop en détails sur ceci dans cet atelier, mais prenez note qu'il est possible de modifier le modèle de base que nous avons vu avec:

• des fonctions plus compliquées en changeant la fonction de base (par exemple, cyclique)
• des modèles à effets mixtes, en utilisant les librairies gamm ou gamm4
• d'autres distributions: tout ce que vous pouvez faire avec un GLM (tel que spécifier l'argument "family") est faisable avec les GAMs

Nous allons d'abord jeter un coup d'oeil au changement de la fonction de base, suivit d'une intro rapide aux GAMMs, puis un exemple de comment spécifier la distribution des données...

edf Ref.df F p-value

s(nottem_year) 2.333879 2.906998 2.528043 0.06007933

s(nottem_month) 4.923943 8.000000 390.029032 0.00000000

k: nombre de «nœuds» - détermine la limite supérieure du nombre de fonctions de base utilisées pour construire la courbe

• # fonctions de base est représenté par l'edf
• par défaut, la limite supérieure de k pour s() est ~9 et de 5 pour chaque dimension de te() et ti() (k devrait être < #données uniques)

d - spécifie quels variables d'une interaction qui se trouvent sur la même échelle lorsque qu'on utilise te() ou ti()

• Par exemple, te(temps, largeur, hauteur, d=c(1,2)), indique que la largeur et la hauteur sont sur la même échelle, mais que temps ne l'est pas

bs - spécifie la fonction de base sous-jacente

• pour s() on utilise "tp" (thin plate regression spline) par défaut, et pour te() et ti() on utilise la base "cr" (cubic regression spline) par défaut.

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.68921 0.0011

year_gam_AR2$lme 3 7 1101.598 1125.962 -543.7988 2 vs 3 1.62082 0.2030

gam.vcomp(gamm.ex)

s():

• pour modéliser un terme lisse 1-dimensionnel, ou pour modéliser une interaction entre des variables mesurées sur la même échelle

te():

• pour modéliser une surface d'interaction 2- ou n-dimensionnel entre des variables qui ne sont pas sur la même échelle. Comprend les effets «principaux».

ti():

• pour modéliser une surface d'interaction 2- ou n-dimensionnel qui ne comprend pas les effets «principaux».

y x0 x1 x2 x3 [...] fac

1 4.925515 0.04879793 0.273053853 0.59921214 0.18689231 1

2 10.352629 0.73410097 0.312996866 0.99805171 0.03375583 2

3 13.527934 0.32637799 0.165864279 0.52143350 0.99192673 3

4 15.511255 0.33131391 0.009219346 0.07767554 0.61888402 4

5 11.455404 0.58911290 0.469944164 0.48423685 0.22028845 1

6 10.124126 0.04054477 0.009816532 0.42261209 0.28148881 2

> plot_smooth(gamm_smooth, view="x0", rm.ranef=TRUE,

main="intercept + s(x0)", rug=FALSE)

> plot_smooth(gamm_smooth, view="x0", cond=list(fac="1"),

main="... + s(x0, fac)", col='orange',

ylim=c(7,21), rug=FALSE)

> plot_smooth(gamm_smooth, view="x0", cond=list(fac="2"),

add=TRUE, col='red', xpd=TRUE)

> plot_smooth(gamm_smooth, view="x0", cond=list(fac="3"),

add=TRUE, col='purple', xpd=TRUE)

> plot_smooth(gamm_smooth, view="x0", cond=list(fac="4"),

add=TRUE, col='turquoise', xpd=TRUE)

> 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")

> print(nottem_plot)

Standard deviations and 0.95 confidence intervals:

std.dev lower upper

s(fac) 4.267269 1.904181 9.562950

scale 3.925177 3.661113 4.208288

Le modèle AR(1) prévoit une augmentation significative de l'ajustement du modèle, mais il y a très peu d'amélioration en allant au modèle AR(2).

Fonctions non linéaire pour la tendance et les caractéristiques saisonniers

Finalement, tous ces modèles mixes peuvent être compare en utilisant la fonction anova() pour trouver le meilleur modèle.

Il y a une hausse d'environ 1 - 1,5ºC au cours de la série, mais dans 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.

Site web: http://qcbs.ca/wiki/r_atelier8

Visualiser la tendance au fil du temps

Enfin, pour nous aider à interpréter les résultats, nous pouvons re-transformer l'effet à des proportions avec la fonction plot_smooth() (section 6.4).

Nombre de succès relatifs aux échecs en fonction de x1 (section 6.2)

Nous allons tester ceci en utilisant un GAM logistique (section 6.3).

Comment définir la distribution pour des données non normales ou d'abondance?

Enfin, nous allons voir les différentes façons de représenter cet relation graphiquement (section 6.4).

emptyPlot(range(gam_data3$x1), c(0,1), h=.5)

avg <- aggregate(prop ~ x1, data=gam_data3, mean, na.rm=TRUE)

lines(avg$x1, avg$prop, col="orange",lwd=2)

Rappel : le modèle utilise le nombre de succès vs échecs pour calculer le logit, qui est la logarithme du rapport entre les succès et échecs:

- succès = échecs: le rapport est de 1 et le logit est 0 (log (1) = 0).

- succès > échecs: le ratio est supérieur à 1 et le logit a une valeur positive (par exemple, log(2) = 0,69).

- succès < échecs, le ratio est inférieur à 1 et le logit a une valeur négative (par exemple, log(0,5) = -0.69).

Valeurs ajustées, ordonné inclus :

Succès et échecs égale jusqu'à ~400

> 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)

Nous allons maintenant voir un bref aperçu de comment d'utiliser les GAMs lorsque la variable réponse ne suit pas une distribution normale ou que les données soient des abondances ou proportions, utilisant un ensemble de données où un répartition binomiale sera nécessaire (mesurant le nombre de succès et d'échecs d'un événement donné; section 6.1)

Resources

Ici, l'estimé est positif ce qui signifie, qu'en moyenne, il n'y a plus de succès que d'échecs.

llustration du principe de validation croisée

Regardons les GAMs de plus près

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:

Un exemple simple : une base polynomiale

Contrôler le degré de lissage avec des splines de régression pénalisés

Un autre exemple : une base de spline cubique

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).

D'après vous, est-ce que cette tendance est linéaire ou non linéaire?

Selection de lambda: principe de validation croisée

Contribution/ effet partiel :

Au fil du temps, la valeur logit augmente (succès augmentent et les échecs diminuent).

Note supplémentaire sur les degrés de liberté effectifs (edf)

régression spline pénalisé modeliser sur les données restantes (º)

observation omise (•)

Logit

> gam_data3 <- read.csv("other_dist.csv")

> summary(gam_data3)

> str(gam_data3)

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, xi :

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

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.

Il existe beaucoup d'information sur les GAMs ... ce fut juste un aperçu.

Pouvez tester la colinéarité des variables et la sélection de modèle:

http://distancesampling.org/workshops/duke-spatial-2015/practicals/3-advanced-dsms-solutions.html

• utilisez "select = TRUE" lors de l'ajustement du modèle (pénalité suppl. sur chaque terme lisse et élimine les relations linéaires)

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 :

Au lieu de contrôler le lissage 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) :

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 :

'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 ...

spline

Qu'est-ce que l'ordonné représente dans ce modèle?

Qu'est-ce que le terme de lissage indique?

et f(x) devient :

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.

Qu'est-ce que ces graphiques nous renseignent vis à vis les succès et échecs?

Simon Wood, l'auteur de la librairie mgcv, a un site très utile sur les GAMs :

http://people.bath.ac.uk/sw283/mgcv/

lI a aussi écrit le livre, Generalized Additive Models: An Introduction with R, que nous avons utilisé comme référence pour cet atelier.

Combien de degrés de liberté a un GAM ?

Au lieu de fournir lambda, gam() utilise les degrés de liberté effectifs (edf) pour quantifier la complexité du modèle (justifie notre intuition que les degrés de liberté quantifient la complexité).

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.

Pour estimer 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, bi(x), définissant l'espace des fonctions dont f est un élément:

il peut être modélisé en minimisant :

nœuds

et le modèle complet devient :

É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é.

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.

Ajustement faible par rapport aux données et ne fait pas mieux avec le point manquant

Ajustement du signal sous-jacent est très bien, le lissage passe à travers le bruit et la donnée manquante est bien prédite.

Le matériel de cet atelier a également été obtenu à partir des liens 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.

Ajuste le bruit aussi bien que le signal et la variabilité induite l'amène à prédire la donnée manquante plutôt mal.

Quant lambda tend vers infinité, le modèle devient linéaire.

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.

Chaque section du spline a des coefficients différents.

Y a-t-il une meilleure façon de sélectionner les emplacements des nœuds?