From 3142f2eb847b6dc5a0040b68b592baef1aa8aace Mon Sep 17 00:00:00 2001 From: Daniel Schoenig Date: Wed, 10 Mar 2021 17:47:21 -0500 Subject: [PATCH 1/4] Other distributions with real data; model checking --- pres-en/workshop08-pres-en.Rmd | 562 ++++++++++++++++++++++++--------- pres-fr/workshop08-pres-fr.Rmd | 553 +++++++++++++++++++++++--------- 2 files changed, 817 insertions(+), 298 deletions(-) diff --git a/pres-en/workshop08-pres-en.Rmd b/pres-en/workshop08-pres-en.Rmd index fad4499..35bd929 100644 --- a/pres-en/workshop08-pres-en.Rmd +++ b/pres-en/workshop08-pres-en.Rmd @@ -77,8 +77,8 @@ install.packages(c('ggplot2', 'itsadug', 'mgcv')) 2. Introduction to GAM 3. GAM with multiple smooth terms 4. Interactions -5. Changing basis -6. Other distributions +5. Generalization of the additive model +6. Changing basis 7. Quick intro to GAMM --- @@ -291,7 +291,7 @@ data_plot <- data_plot + data_plot ``` -.comment[Note: as opposed to one fixed coefficient, $\beta$ in linear regression, the smooth function can continually change over the range of the predictor x] +Note: as opposed to one fixed coefficient, $\beta$ in linear regression, the smooth function can continually change over the range of the predictor $x$ --- exclude:true @@ -695,9 +695,9 @@ There is a pronounced difference in bioluminescence between the seasons. basic_summary$s.table ``` -The `edf` shown in the `s.table` are the **effective degrees of freedom** (EDF) of the the smooth term –- essentially, more EDF imply more complex, wiggly splines. +The `edf` shown in the `s.table` are the **effective degrees of freedom** (EDF) of the the smooth term - essentially, more EDF imply more complex, wiggly splines. -- A value close to 1 tend to be close to a linear term. +- A value close to 1 tends to be close to a linear term. - A higher value means that the spline is non-linear. @@ -719,6 +719,7 @@ An upper bound on the **EDF** is determined by the basis dimension $k$ for each In practice, the exact choice of $k$ is arbitrary, but it should be **large enough** to accomodate a sufficiently complex smooth function. +We will talk about choosing $k$ in section 5. --- # GAM with multiple variables @@ -974,6 +975,7 @@ AIC(two_smooth_model, three_term_model, test = "Chisq") --- + class: inverse, center, middle ## 4. Interactions @@ -1095,272 +1097,364 @@ AIC(two_smooth_model, smooth_interact) The model with the interaction between `s(SampleDepth)` and `s(RelativeDepth)` have a lower AIC and the 2D plot nicely illustrates this non-linear interactions, where `Sources` is lower at high values of `SampleDepth` but mid to high values of `RelativeDepth`. + + --- class: inverse, center, middle -## 5. Changing basis +## 5. Generalization of the additive model --- + # Generalization of the additive model The basic additive model can be extended in the following way: -1. Using different kinds of basis functions, -2. Using other distributions for the response variable with the `family` argument (just as in a GLM), -3. Using different kinds of random effects to fit mixed effect models. +1. Using other **distributions** for the response variable with the `family` argument (just as in a GLM), +2. Using different kinds of **basis functions**, +3. Using different kinds of **random effects** to fit mixed effect models. We will now go over these aspects. --- -# Other smooth functions -To model a non-linear smooth variable or surface, smooth functions can be built in different ways: +# Generalized additive models -`s()` ![:faic](arrow-right) for modeling a 1-dimensional smooth or for modeling interactions among variables measured using the same unit and the same scale +We have so far worked with simple (Gaussian) additive models, the non-linear equivalent to a linear model. -`te()` ![:faic](arrow-right) for modeling 2- or n-dimensional interaction surfaces of variables that are not on the same scale. Includes main effects. +-- -`ti()` ![:faic](arrow-right) for modeling 2- or n-dimensional interaction surfaces that do not include the main effects. +But what can we do if: +- the observations of the response variable do **not follow a normal distribution**? +- the **variance is not constant**? (heteroscedasticity) + +-- + +.comment[These cases are very common!] + +Just like generalized linear models (GLM), we can formulate **generalized** additive models to deal with these issues --- -# Parameters of smooth functions -The smooth functions have several parameters that can be set to change their behavior. The most often-used parameters are : +# Generalized additive models -`k` ![:faic](arrow-right) basis dimension - - determines the upper bound of the number of base functions used to build the curve. - - constrains the wigglyness of a smooth. - - the default $k$ is 10 for `s()`, and 5 for `te()` and `ti()`. - - k should be < the number of unique data points. - - the comlexity (i.e. non-linearity) of a smooth function in a fitted model is reflected by its effective degrees of freedom (**EDF**). +Lets recall the interaction model for the bioluminescence data: + +```{r} +smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth), data = isit) + +summary(smooth_interact)$p.table + +summary(smooth_interact)$s.table +``` --- -# Parameters of smooth functions -The smooth functions have several parameters that can be set to change their behavior. The most often-used parameters are : +# GAM model checking -`d` ![:faic](arrow-right) specifies that predictors in the interaction are on the same scale or dimension (only used in `te()` and `ti()`). - - For example, in `te(Time, width, height, d=c(1,2))`, indicates that `width` and `height` are one the same scale, but not `Time`. +As for a GLM, it is essential to check whether we correctly specifified the model, especially the *distribution* for the response. -`bs` ![:faic](arrow-right) specifies the type of underlying base functions. - - the default for `s()` is `tp` (thin plate regression spline) and for `te()` and `ti()` is `cr` (cubic regression spline). +**We have two verify:** - +1. The choice of basis dimension `k`. +2. The residuals plots (just as for a GLM). ---- -# Example smooth for cyclical data -Cyclical data is a good example where changing basis is useful: you want the predictor to match at the ends. +-- -Let's use a time series of climate data, with monthly measurements, and see if there's a temporal trend in yearly temperature. +
+Useful functions included in `mgcv`: -```{r, eval = FALSE} -data(nottem) # Nottingham temperature time series -n_years <- length(nottem)/12 -nottem_month <- rep(1:12, times = n_years) -nottem_year <- rep(1920:(1920 + n_years - 1), each = 12) -qplot(nottem_month, nottem, colour = factor(nottem_year), geom = "line") + - theme_bw() -``` +- `k.check()` performs a basis dimension check. +- `gam.check()` produces residual plots (and also calls `k.check()`). --- -# Example smooth for cyclical data -```{r, echo = FALSE, fig.height=8, fig.width=12} -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) -qplot(nottem_month, nottem, colour = factor(nottem_year), geom = "line") + - theme_bw() -``` +# GAM model checking ---- -# Example smooth for cyclical data +##### First step: -We can model both the cyclic change of temperature across months and the non-linear trend through years, using a cyclical cubic spline, or `cc`, for the month variable and a regular smooth for the year variable. +Have we chosen `k` large enough? -```{r, fig.height=4, fig.width=8} -year_gam <- gam(nottem ~ s(nottem_year) + s(nottem_month, bs = "cc")) -summary(year_gam)$s.table +.comment[The default for smooth interactions is `k = 30`] + +```{r} +k.check(smooth_interact) ``` +-- + +The **EDF are very close to** `k`, this could be problematic. + --- -# Example smooth for cyclical data -```{r, fig.height = 5, fig.width = 8} -plot(year_gam, page = 1, scale = 0) +# GAM model checking + +Let's refit the model with a larger `k`: + +```{r} +smooth_interact_k30 <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), + data = isit) +summary(smooth_interact_k30)$p.table +summary(smooth_interact_k30)$s.table ``` -There is about 1-1.5 degree rise in temperature over the period, but within a given year there is about 20 degrees variation in temperature, on average. The actual data vary around these values and that is the unexplained variance. +-- ---- -class: inverse, center, middle +Is `k` large enough this time? + +```{r} +k.check(smooth_interact_k30) +``` + +-- -## 6. Other distributions +Looks better, let's replace the old model: +```{r} +smooth_interact <- smooth_interact_k30 +``` --- -# GAM using other distributions -Let's now take a look on how to use GAMs when the response variable does not follow a normal distributions and is either count or proportion data (e.g., Gamma, binomial, Poisson, negative binomial). +# GAM model checking -We will use an example dataset where a binomial distribution is needed; the response variable represents the number of successes vs failures over the course of an experiment. +##### Second Step: -```{r} -gam_data3 <- read.csv("data/other_dist.csv") -str(gam_data3) +Let's look at the residual plots, using `gam.check()`: + +```{r, eval = FALSE} +par(mfrow = c(2,2)) # Show all 4 plots side by side +gam.check(smooth_interact) ``` - +.comment[In addition to the plots, `gam.check()` also provides the output of `k.check()`.] --- -# GAM using other distributions -```{r, fig.height=5, fig.width=5} -plot(range(gam_data3$x1), c(0,1), type = "n", - main = "Probability of successes over time", - ylab = "Probability", xlab = "x1 (time)") -abline(h = 0.5) +# GAM model checking -avg <- aggregate(prop ~ x1, data = gam_data3, mean) -lines(avg$x1, avg$prop, col = "orange", lwd = 2) +```{r, echo=FALSE, results = FALSE, fig.height=6, fig.width=8} +par(mfrow=c(2,2), mar = c(4,4,2,1.1), oma =c(0,0,0,0)) +gam.check(smooth_interact) ``` +-- +
+ +.alert[Pronounced heteroscedasticity and patterns in the residuals] + + +??? + +- These plots are a little different than those produced by `plot` for a linear model (e.g. no leverage plot) +- Participants should already be familiar with residual plots (they are explaned more in detail in workshops 4 and 6) +- Obvious problem: heteroscedasticity, visible in residuals vs. linear predictor. +- Another problem: a few strong outliers (visible in QQ plot, response vs. fitted) + --- -# GAM using other distributions -We will test if this trend is linear or not using a logistic GAM (we use a binomial family distribution given that our response is proportion data). +# Other distributions -```{r, warning=F} -prop_model <- gam(prop ~ s(x1), data = gam_data3, weights = total, family = "binomial") -prop_summary <- summary(prop_model) -``` - +For our interaction model, we need a probability distribution that allows the **variance to increase with the mean**. -- -.comment[What does the intercept represent in this model?] +One family of distributions that has this property and that works well in a GAM is the **Tweedie** family. -```{r} -prop_summary$p.table +A common link function for *Tweedie* distributions is the $log$. + +-- + +
+As in a GLM, we can use the `family = ` argument in `gam()` to fit models with other distributions (including distributions such as `binomial`, `poisson`, `gamma` etc.). + +To get an overview of families available in `mgcv`: +```{r, eval = FALSE} +?family.mgcv ``` + + + +--- + +# Challenge 3 ![:cube]() + +1. Fit a new model `smooth_interact_tw` with the same formula as the `smooth_interact` model but with a distribution from the *Tweedie* family (instead of the normal distribution) and `log` link function. You can do so by using `family = tw(link = "log")` inside `gam()`. +2. Check the choice of `k` and the residual plots for the new model. +3. Compare `smooth_interact_tw` with `smooth_interact`. Which one would you choose? + -- -.comment[What does the smooth term indicate?] +
+ +.comment[Hint:] ```{r} -prop_summary$s.table +# Normal distribution +smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), + data = isit) +# Tweedie family with log link +smooth_interact_tw <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), + family = tw(link = "log"), + data = isit) +``` + + +--- + +# Challenge 3 - Solution ![:cube]() + +Fit the model: + +```{r} +smooth_interact_tw <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), + family = tw(link = "log"), + data = isit) +summary(smooth_interact_tw)$p.table +summary(smooth_interact_tw)$s.table ``` +-- + +Check the basis dimension: +```{r} +k.check(smooth_interact_tw) +``` --- -# GAM using other distributions -.comment[What does the intercept represent in this model?] +# Challenge 3 - Solution ![:cube]() -```{r, echo = FALSE} -prop_summary$p.table +Residual plots: + +```{r, eval = FALSE} +par(mfrow=c(2,2)) +gam.check(smooth_interact_tw) ``` -Recall that the model uses the count data to calculate the logit, which is the log odds ratio between successes and failures: -.small[ -- If successes = failures, the ratio = 1 and the logit is 0 (log(1) = 0). -- If successes > failures, the ratio > 1 and the logit has a positive value (log(2) = 0.69). -- If successes < failures, the ratio < 1 and the logit has a negative value (log(.5) = -0.69). -] --- +```{r, echo = FALSE, results = FALSE, fig.height = 6, fig.width = 8} +par(mfrow=c(2,2), mar = c(4,4,2,1.1), oma =c(0,0,0,0)) +gam.check(smooth_interact_tw) +``` -> Here, the estimated intercept coefficient is positive, which means that there are more successes than failures overall. + +??? +- The residuals look much better, but it is clear that something is missing from the model. This could be a spatial affect (longtitude and latitude), or a random effect (e.g. based on `Station`). --- -# GAM using other distributions -.comment[What does the smooth term indicate?] +# Challenge 3 - Solution ![:cube]() -```{r, echo = FALSE} -prop_summary$s.table +Compare the models: + +```{r} +AIC(smooth_interact, smooth_interact_tw) ``` -This represents how the log odds of successes vs failures change over time (x1). +.comment[AIC allows us to compare models that are based on different distributions!] -- -> As the **EDF** > 1, the proportion of successes increases faster over time +Using a *Tweedie* instead of a *normal* distribution greatly improves our model! -```{r, fig.height=4.3, fig.width=4.3, eval = FALSE} -plot(prop_model) -``` -```{r, fig.height=4.3, fig.width=4.3, echo = FALSE} -par(mar = c(4,4,0,0)) -plot(prop_model) -``` +--- +class: inverse, center, middle + +## 6. Changing basis --- -# Visualizing the trend over time +# Other smooth functions -There are different ways this relationship can be represented graphically: +To model a non-linear smooth variable or surface, smooth functions can be built in different ways: -- **partial effects** are the isolated effects of one particular predictor or interaction. If you visualize your GAM model with `plot()`, you get the partial effects. -- **summed effects** are the predicted response measures for a given value or level of predictors. If you visualize your GAM model with `itsadug::plot_smooth()`, you get the summed effects. +`s()` ![:faic](arrow-right) for modeling a 1-dimensional smooth or for modeling interactions among variables measured using the same unit and the same scale +`te()` ![:faic](arrow-right) for modeling 2- or n-dimensional interaction surfaces of variables that are not on the same scale. Includes main effects. + +`ti()` ![:faic](arrow-right) for modeling 2- or n-dimensional interaction surfaces that do not include the main effects. --- -# Visualizing the trend over time +# Parameters of smooth functions -What do these plots tell us about successes vs failures? +The smooth functions have several parameters that can be set to change their behavior. The most often-used parameters are : -```{r, echo = FALSE, fig.height=4.5, fig.width=9} -library(itsadug) -par(mfrow=c(1,2), mar = c(4,4,0,0)) -plot(prop_model, select = 1, scale = 0, shade=TRUE) -abline(h=0) +`k` ![:faic](arrow-right) basis dimension + - determines the upper bound of the number of base functions used to build the curve. + - constrains the wigglyness of a smooth. + - the default $k$ is 10 for `s()`, and 5 for `te()` and `ti()`. + - k should be < the number of unique data points. + - the comlexity (i.e. non-linearity) of a smooth function in a fitted model is reflected by its effective degrees of freedom (**EDF**). -out <- plot_smooth(prop_model, view = "x1",main = "", print.summary = F) -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, "successes > failures", col = 'red', font = 2) -``` +--- +# Parameters of smooth functions +The smooth functions have several parameters that can be set to change their behavior. The most often-used parameters are : -.pull-left[ -**Contribution / partial effect** +`d` ![:faic](arrow-right) specifies that predictors in the interaction are on the same scale or dimension (only used in `te()` and `ti()`). + - For example, in `te(Time, width, height, d=c(1,2))`, indicates that `width` and `height` are one the same scale, but not `Time`. -Over time the log odds increases, so over time successes increase and failures decrease.] +`bs` ![:faic](arrow-right) specifies the type of underlying base functions. + - the default for `s()` is `tp` (thin plate regression spline) and for `te()` and `ti()` is `cr` (cubic regression spline). -.pull-right[ -**Fitted values, summed effect, intercept included** + -Equal amounts of successes and failures up to x1=400. -] +--- +# Example smooth for cyclical data + +Cyclical data is a good example where changing basis is useful: you want the predictor to match at the ends. + +Let's use a time series of climate data, with monthly measurements, and see if there's a temporal trend in yearly temperature. + +```{r, eval = FALSE} +data(nottem) # Nottingham temperature time series +n_years <- length(nottem)/12 +nottem_month <- rep(1:12, times = n_years) +nottem_year <- rep(1920:(1920 + n_years - 1), each = 12) +qplot(nottem_month, nottem, colour = factor(nottem_year), geom = "line") + + theme_bw() +``` - --- -# Visualizing the trend over time +# Example smooth for cyclical data -Lastly, to help interpret the results, we could transform the summed effects back to proportions with the function `itsadug::plot_smooth()`: +```{r, echo = FALSE, fig.height=8, fig.width=12} +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) +qplot(nottem_month, nottem, colour = factor(nottem_year), geom = "line") + + theme_bw() +``` -```{r, echo = -1, fig.height=4, fig.width==4} -par(mar=c(3.8,4,0,0)) -plot_smooth(prop_model, view = "x1", main = "", - transform = plogis, ylim = c(0,1), print.summary = F) -abline(h = 0.5, v = diff$start, col = 'red', lty = 2) +--- +# Example smooth for cyclical data + +We can model both the cyclic change of temperature across months and the non-linear trend through years, using a cyclical cubic spline, or `cc`, for the month variable and a regular smooth for the year variable. + +```{r, fig.height=4, fig.width=8} +year_gam <- gam(nottem ~ s(nottem_year) + s(nottem_month, bs = "cc")) +summary(year_gam)$s.table ``` -As in the logit plot, the proportion of successes increases above 0.5 at x1=400. +--- +# Example smooth for cyclical data + +```{r, fig.height = 5, fig.width = 8} +plot(year_gam, page = 1, scale = 0) +``` + +There is about 1-1.5 degree rise in temperature over the period, but within a given year there is about 20 degrees variation in temperature, on average. The actual data vary around these values and that is the unexplained variance. - --- @@ -1709,3 +1803,171 @@ class: inverse, center, bottom # Thank you for attending this workshop! ![:scale 50%](images/qcbs_logo.png) + +--- +class: inverse, center, middle + +# Additional example with other distributions + + +--- +# GAM using other distributions + +Let's now take a look on how to use GAMs when the response variable does not follow a normal distributions and is either count or proportion data (e.g., Gamma, binomial, Poisson, negative binomial). + +We will use an example dataset where a binomial distribution is needed; the response variable represents the number of successes vs failures over the course of an experiment. + +```{r} +gam_data3 <- read.csv("data/other_dist.csv") +str(gam_data3) +``` + + + +--- +# GAM using other distributions + +```{r, fig.height=5, fig.width=5} +plot(range(gam_data3$x1), c(0,1), type = "n", + main = "Probability of successes over time", + ylab = "Probability", xlab = "x1 (time)") +abline(h = 0.5) + +avg <- aggregate(prop ~ x1, data = gam_data3, mean) +lines(avg$x1, avg$prop, col = "orange", lwd = 2) +``` + +--- +# GAM using other distributions + +We will test if this trend is linear or not using a logistic GAM (we use a binomial family distribution given that our response is proportion data). + +```{r, warning=F} +prop_model <- gam(prop ~ s(x1), data = gam_data3, weights = total, family = "binomial") +prop_summary <- summary(prop_model) +``` + + + +-- + +.comment[What does the intercept represent in this model?] + +```{r} +prop_summary$p.table + +``` + +-- + +.comment[What does the smooth term indicate?] + +```{r} +prop_summary$s.table +``` + + + +--- +# GAM using other distributions + +.comment[What does the intercept represent in this model?] + +```{r, echo = FALSE} +prop_summary$p.table +``` + +Recall that the model uses the count data to calculate the logit, which is the log odds ratio between successes and failures: +.small[ +- If successes = failures, the ratio = 1 and the logit is 0 (log(1) = 0). +- If successes > failures, the ratio > 1 and the logit has a positive value (log(2) = 0.69). +- If successes < failures, the ratio < 1 and the logit has a negative value (log(.5) = -0.69). +] + +-- + +> Here, the estimated intercept coefficient is positive, which means that there are more successes than failures overall. + +--- +# GAM using other distributions + +.comment[What does the smooth term indicate?] + +```{r, echo = FALSE} +prop_summary$s.table +``` + +This represents how the log odds of successes vs failures change over time (x1). + +-- + +> As the **EDF** > 1, the proportion of successes increases faster over time + +```{r, fig.height=4.3, fig.width=4.3, eval = FALSE} +plot(prop_model) +``` + +```{r, fig.height=4.3, fig.width=4.3, echo = FALSE} +par(mar = c(4,4,0,0)) +plot(prop_model) +``` + + +--- +# Visualizing the trend over time + +There are different ways this relationship can be represented graphically: + +- **partial effects** are the isolated effects of one particular predictor or interaction. If you visualize your GAM model with `plot()`, you get the partial effects. +- **summed effects** are the predicted response measures for a given value or level of predictors. If you visualize your GAM model with `itsadug::plot_smooth()`, you get the summed effects. + + +--- +# Visualizing the trend over time + +What do these plots tell us about successes vs failures? + +```{r, echo = FALSE, fig.height=4.5, fig.width=9} +library(itsadug) +par(mfrow=c(1,2), mar = c(4,4,0,0)) +plot(prop_model, select = 1, scale = 0, shade=TRUE) +abline(h=0) + +out <- plot_smooth(prop_model, view = "x1",main = "", print.summary = F) +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, "successes > failures", col = 'red', font = 2) +``` + + + +.pull-left[ +**Contribution / partial effect** + +Over time the log odds increases, so over time successes increase and failures decrease.] + +.pull-right[ +**Fitted values, summed effect, intercept included** + +Equal amounts of successes and failures up to x1=400. +] + + +--- +# Visualizing the trend over time + +Lastly, to help interpret the results, we could transform the summed effects back to proportions with the function `itsadug::plot_smooth()`: + +```{r, echo = -1, fig.height=4, fig.width==4} +par(mar=c(3.8,4,0,0)) +plot_smooth(prop_model, view = "x1", main = "", + transform = plogis, ylim = c(0,1), print.summary = F) +abline(h = 0.5, v = diff$start, col = 'red', lty = 2) +``` + +As in the logit plot, the proportion of successes increases above 0.5 at x1=400. + + + diff --git a/pres-fr/workshop08-pres-fr.Rmd b/pres-fr/workshop08-pres-fr.Rmd index 2e5bbbd..2af16ac 100644 --- a/pres-fr/workshop08-pres-fr.Rmd +++ b/pres-fr/workshop08-pres-fr.Rmd @@ -78,9 +78,9 @@ install.packages(c('ggplot2', 'itsadug', 'mgcv')) 2. Introduction aux GAMs 3. GAM avec plusieurs termes non-linéaires 4. Interactions -5. Changer la fonction de base -6. Intro rapide aux GAMMs -7. Autres distributions +5. Généralisation du modèle additif +6. Changer la fonction de base +7. Intro rapide aux GAMMs --- # Objectifs d'apprentissage @@ -1039,7 +1039,8 @@ Le modèle avec l'intéraction entre `s(SampleDepth)` et `s(RelativeDepth)` a un --- class: inverse, center, middle -## 5. Changer la fonction de base +# 5. Généralisation du modèle additif + --- @@ -1047,257 +1048,349 @@ class: inverse, center, middle Le modèle additif de base peut être étendu de plusieurs façons : -1. Utiliser de différents types de fonctions de base, -2. Utiliser d'autres distributions pour la variable de réponse avec l'argument `famille` (comme dans un GLM), -3. Utilisation de différents types d'effets aléatoires pour ajuster des modèles à effets mixtes. +1. Utilisation de différents types d'effets aléatoires pour ajuster des modèles à effets mixtes. +2. Utiliser de différents types de fonctions de base, +3. Utiliser d'autres distributions pour la variable de réponse avec l'argument `famille` (comme dans un GLM), Nous allons maintenant examiner ces aspects. + --- -# Fonctions lisses -Pour modéliser une surface lisse ou non-linéaire, nous pouvons construire des fonctions lisses de différentes manières: +# Modèle additif généralisé -`s()` ![:faic](arrow-right) pour modéliser un terme lisse 1-dimensionnel, ou pour modéliser une intéraction entre des variables mesurées sur la *même échelle* +Jusqu'à présent, nous avons utilisé des modèles additifs simples (gaussiens), l'équivalent non linéaire d'un modèle linéaire. -`te()` ![:faic](arrow-right) 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()` ![:faic](arrow-right) pour modéliser une surface d'interaction 2- ou n-dimensionnel *qui ne comprend pas les effets principaux*. +Mais que pouvons-nous faire si : +- les observations de la variable de réponse ne **suivent pas une distribution normale** ? +- la **variance n'est pas constante** ? (hétéroscédasticité) ---- -# Paramètres des fonctions lisses +-- -Les fonctions lisses ont beaucoup de paramètres qui pourraient changer leur comportement. Les paramètres les plus souvent utilisés sont les suivants : +.comment[Ces situations se produisent fréquemment !] -`k` ![:faic](arrow-right) dimensions de base - - détermine la limite supérieure du nombre de fonctions de base utilisées pour construire la courbe. - - contraint l'ondulation d'une fonction lisse. - - par défaut, $k$ est 10 pour `s()`, et 5 pour chaque dimension de `te()` and `ti()`. - - k devrait être < au nombre de données uniques. - - la comlexité (ou la non-linéarité) d'une fonction lisse dans un modèle ajusté est reflétée par ses degrés de liberté effectifs (**EDF**) +Tout comme les modèles linéaires généralisés (GLM), nous pouvons formuler des modèles additifs **généralisés** pour répondre à ces problèmes. --- -# Paramètres des fonctions lisses -Les fonctions lisses ont beaucoup de paramètres qui pourraient changer leur comportement. Les paramètres les plus souvent utilisés sont les suivants : +# Modèle additif généralisé -`d` ![:faic](arrow-right) spécifie quelles variables d'une intéraction se trouvent sur la même échelle lorsqu'on utilise `te()` and `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. +Rappelons le modèle d'interaction pour les données de bioluminescence : -`bs` ![:faic](arrow-right) spécifie la fonction de base sous-jacente. - - pour `s()` on utilise `tp` (*thin plate regression spline*) et pour `te()` et `ti()` on utilise la base `cr` (*cubic regression spline*). +```{r} +smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth), data = isit) - +summary(smooth_interact)$p.table + +summary(smooth_interact)$s.table +``` --- -# Exemple avec des données cycliques -Les données cycliques sont un bon exemple où changer la base est utile : vous voulez que le prédicteur corresponde à la fin. +# Validation d'un GAM -Nous allons utiliser une série temporelle de données climatiques, divisées en mesures mensuelles, afin de déterminer s'il y a une tendance de température annuelle. +Comme pour un GLM, il est essentiel de vérifier si nous avons correctement spécifié le modèle, en particulier la *distribution* de la réponse. + +**Il faut vérifier :** + +1. Le choix des dimensions de base `k`. +2. Les tracés des résidus (comme pour un GLM). + + +-- + +
+Fonctions incluses dans `mgcv` : + +- `k.check()` effectue une vérification des dimensions de base. +- `gam.check()` produit des tracés de résidus (et fournit également la sortie de `k.check()`. -```{r, eval = FALSE} -data(nottem) # Nottingham temperature time series -n_years <- length(nottem)/12 -nottem_month <- rep(1:12, times = n_years) -nottem_year <- rep(1920:(1920 + n_years - 1), each = 12) -qplot(nottem_month, nottem, colour = factor(nottem_year), geom = "line") + - theme_bw() -``` --- -# Exemple avec des données cycliques -```{r, echo = F, fig.height=8, fig.width=12} -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) -qplot(nottem_month, nottem, colour = factor(nottem_year), geom = "line") + - theme_bw() +# Validation d'un GAM + +##### Première étape : + +Avons-nous choisi `k` assez grand ? + +.comment [Le défaut pour les interactions lisses est `k = 30`]. + +```{r} +k.check(smooth_interact) ``` +-- + +Les **EDF se rapprochent beaucoup de** `k`, cela pourrait être problématique. + --- -# Exemple avec des données cycliques -Nous pouvons modéliser le changement cyclique de température à travers les mois et la tendance non-linéaire à travers les années, en utilisant une spline cubique, ou `cc` pour modéliser les effets de mois ainsi qu'un terme non-linéaire pour la variable année. +# Validation d'un GAM -```{r, fig.height=4, fig.width=8} -year_gam <- gam(nottem ~ s(nottem_year) + s(nottem_month, bs = "cc")) -summary(year_gam)$s.table +Refaisons le modèle avec un `k` plus grand : + +```{r} +smooth_interact_k30 <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), + data = isit) +summary(smooth_interact_k30)$p.table +summary(smooth_interact_k30)$s.table ``` ---- -# Exemple avec des données cycliques +-- -```{r, fig.height=5, fig.width=8} -plot(year_gam, page = 1, scale = 0) +Est-ce que `k` est assez grand ? + +```{r} +k.check(smooth_interact_k30) ``` -Il y a une hausse d'environ 1-1.5 degrés au cours de la série, mais au cours d'une année, il y a une variation d'environ 20 degrés. Les données réelles varient autour de ces valeurs prédites et ceci représente donc la variance inexpliquée. ---- -class: inverse, center, middle +-- + +C'est mieux, remplaçons l'ancien modèle : -## 6. Autres distributions +```{r} +smooth_interact <- smooth_interact_k30 +``` --- -# GAM avec d'autres distributions -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). +# Validation d'un GAM -Nous allons utiliser un exemple de données où une répartition binomiale sera nécessaire; 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. +##### Deuxième étape : -```{r} -gam_data3 <- read.csv("data/other_dist.csv") -str(gam_data3) +Regardons les tracés des résidus, en utilisant `gam.check()` : + +```{r, eval = FALSE} +par(mfrow = c(2,2)) # Afficher les 4 tracés à la fois +gam.check(smooth_interact) ``` - +.comment [En plus des tracés, `gam.check()` fournit également la sortie de `k.check()`]. --- -# GAM avec d'autres distributions -```{r, fig.height=4, fig.width=4} -plot(range(gam_data3$x1), c(0,1), type = "n", - main = "Probabilités de succès dans le temps", - ylab = "Probabilité", xlab = "x1 (temps)") -abline(h = 0.5) +# Validation d'un GAM -avg <- aggregate(prop ~ x1, data=gam_data3, mean) -lines(avg$x1, avg$prop, col = "orange", lwd = 2) +```{r, echo=FALSE, results = FALSE, fig.height=6, fig.width=8} +par(mfrow=c(2,2), mar = c(4,4,2,1.1), oma =c(0,0,0,0)) +gam.check(smooth_interact) ``` +-- +
+ +.alert[Hétéroscédasticité marquée et tendances dans les résidus] + + +??? + +- Ces tracés sont un peu différents de ceux produits par `plot` pour un modèle linéaire (par exemple, pas de tracé de levier). +- Les participants devraient déjà être familiarisés avec les graphiques de résidus (ils sont expliqués plus en détail dans les ateliers 4 et 6). +- Problème : hétéroscédasticité, visible dans "residuals vs. linear predictor". +- Autre problème : quelques observations extrêmes(visibles dans le tracé QQ, et "response vs. fitted"). + + --- -# GAM avec d'autres distributions -Nous allons tester si cette tendance est linéaire ou non avec un GAM logistique (nous utilisons une famille de distributions binomiales parce que nous avons des données de proportion). +# Autres distributions -```{r, warning=F} -prop_model <- gam(prop ~ s(x1), data = gam_data3, weights = total, family = "binomial") -prop_summary <- summary(prop_model) -``` - +Pour notre modèle d'interaction, nous avons besoin d'une distribution de probabilité qui permet à la **variance d'augmenter avec la moyenne**. -- -.comment[Qu'est ce que représente l'intercepte dans ce modèle?] +Une famille de distributions qui possède cette propriété et qui fonctionne bien dans un GAM est la famille **Tweedie**. -```{r} -prop_summary$p.table -``` +Une fonction de liaison commune pour les distributions *Tweedie* est le $log$. -- -.comment[Qu'est ce que le terme de lissage indique?] +
-```{r} -prop_summary$s.table +Comme dans un GLM, nous pouvons utiliser l'argument `family = ` dans `gam()` pour ajuster des modèles avec d'autres distributions (y compris des distributions telles que `binomial`, `poisson`, `gamma` etc). + +Pour en savoir plus sur les familles disponibles dans `mgcv` : +```{r, eval = FALSE} +?family.mgcv ``` + --- -# GAM avec d'autres distributions -```{r, echo = FALSE} -prop_summary$p.table +# Défi 3 ![:cube]() + +1. Ajuster un nouveau modèle `smooth_interact_tw` avec la même formule que le modèle `smooth_interact` mais avec une distribution de la famille *Tweedie* (au lieu de la distribution normale) et `log` comme fonction de liaison. Pour ce faire, on peut utiliser `family = tw(link = "log")` dans `gam()`. +2. Vérifier le choix de `k` et les tracés de résidus pour le nouveau modèle. +3. Comparer `smooth_interact_tw` avec `smooth_interact`. Lequel est meilleur ? + +-- + +
+ +.comment[Indice :] + +```{r} +# Distribution normale +smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), + data = isit) +# Tweedie avec log comme fonction de liaison +smooth_interact_tw <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), + family = tw(link = "log"), + data = isit) ``` -.comment[Que représente l'intercepte dans ce modèle?] -**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 : +--- -.small[ -- Si succès = échecs, le rapport = 1 et le logit est de 0 (log(1) = 0). -- Si succès > échecs, le rapport > 1 et le logit a une valeur positive (log(2) = 0.69). -- Si succès < échecs, le rapport < 1 et le logit a une valeur négative (log(.5) = -0.69). -] +# Défi 3 - Solution ![:cube]() + +Ajuster le modèle : + +```{r} +smooth_interact_tw <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), + family = tw(link = "log"), + data = isit) +summary(smooth_interact_tw)$p.table +summary(smooth_interact_tw)$s.table +``` -- -> Ici, l'estimé est positif ce qui signifie, qu'en moyenne, il y a plus de succès que d'échecs. +Vérifier les dimensions de base : +```{r} +k.check(smooth_interact_tw) +``` --- -# GAM avec d'autres distributions -```{r, echo = FALSE} -prop_summary$s.table +# Challenge 3 - Solution ![:cube]() + +Tracé des résidus : + +```{r, eval = FALSE} +par(mfrow=c(2,2)) +gam.check(smooth_interact_tw) ``` -.comment[Qu'est ce que le terme de lissage indique?] -Cela représente comment le ratio de succès vs échecs change sur l'échelle de $x1$. +```{r, echo = FALSE, results = FALSE, fig.height = 6, fig.width = 8} +par(mfrow=c(2,2), mar = c(4,4,2,1.1), oma =c(0,0,0,0)) +gam.check(smooth_interact_tw) +``` --- -> Puisque les **EDF** > 1, la proportion des succès augmente plus rapidement avec $x1$ +??? +- Les résidus semblent meilleurs, mais il est évident que le modèle manque quelque chose. Il pourrait s'agir d'un effet spatial (longitude et latitude), ou d'un effet aléatoire (par exemple basé sur `Station`). -```{r, fig.height=4.3, fig.width=4.3, eval = FALSE} -plot(prop_model) -``` +--- -```{r, fig.height=4.3, fig.width=4.3, echo = FALSE} -par(mar = c(4,4,0,0)) -plot(prop_model) +# Défi 3 - Solution ![:cube]() + +Comparer les modèles: + +```{r} +AIC(smooth_interact, smooth_interact_tw) ``` +.comment[L'AIC nous permet de comparer des modèles qui sont basés sur des distributions différentes !] + +-- + +Utiliser une distribution *Tweedie* au lieu d'une distribution *normale* améliore beaucoup notre modèle ! + + --- -# Visualiser la tendance au fil du temps +class: inverse, center, middle -Il y a différente façon de représenter cette relation graphiquement : +## 6. Changer la fonction de base -- **Contribution/effet partiel** correspond aux effets isolés d'une interaction ou prédiction particulière. Si vous visualisez votre modèle GAM avec `plot()`, vous obtenez les effets partiels. -- **effets additionnés** correspond aux mesures réponse prédites pour une valeur ou niveau donné de prédicteurs. Si vous visualisez votre GAM avec `itsadug::plot_smooth()`, vous obtenez les effets additionnés. +--- +# Fonctions lisses + +Pour modéliser une surface lisse ou non-linéaire, nous pouvons construire des fonctions lisses de différentes manières: + +`s()` ![:faic](arrow-right) pour modéliser un terme lisse 1-dimensionnel, ou pour modéliser une intéraction entre des variables mesurées sur la *même échelle* + +`te()` ![:faic](arrow-right) 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()` ![:faic](arrow-right) pour modéliser une surface d'interaction 2- ou n-dimensionnel *qui ne comprend pas les effets principaux*. --- -# Visualiser la tendance au fil du temps +# Paramètres des fonctions lisses -Que nous disent ces graphes sur les succès et échecs? +Les fonctions lisses ont beaucoup de paramètres qui pourraient changer leur comportement. Les paramètres les plus souvent utilisés sont les suivants : -```{r, echo = FALSE, fig.height = 4.5, fig.width = 9} -library(itsadug) -par(mfrow=c(1,2), mar = c(4,4,0,0)) -plot(prop_model, select = 1, scale = 0, shade=TRUE) -abline(h=0) +`k` ![:faic](arrow-right) dimensions de base + - détermine la limite supérieure du nombre de fonctions de base utilisées pour construire la courbe. + - contraint l'ondulation d'une fonction lisse. + - par défaut, $k$ est 10 pour `s()`, et 5 pour chaque dimension de `te()` and `ti()`. + - k devrait être < au nombre de données uniques. + - la comlexité (ou la non-linéarité) d'une fonction lisse dans un modèle ajusté est reflétée par ses degrés de liberté effectifs (**EDF**) -out <- plot_smooth(prop_model, view = "x1",main = "", print.summary = F) -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, "succès > échecs", col = 'red', font = 2) -``` +--- +# Paramètres des fonctions lisses + +Les fonctions lisses ont beaucoup de paramètres qui pourraient changer leur comportement. Les paramètres les plus souvent utilisés sont les suivants : +`d` ![:faic](arrow-right) spécifie quelles variables d'une intéraction se trouvent sur la même échelle lorsqu'on utilise `te()` and `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. -.pull-left[ -**Contribution / effets partiels** +`bs` ![:faic](arrow-right) spécifie la fonction de base sous-jacente. + - pour `s()` on utilise `tp` (*thin plate regression spline*) et pour `te()` et `ti()` on utilise la base `cr` (*cubic regression spline*). -La valeur logit augmente, donc les succès augmentent et les échecs diminuent.] + -.pull-right[ -**Valeurs ajustées, effets additionnés, intercepte inclu** +--- +# Exemple avec des données cycliques -Quantités égales de succès et d'échecs jusqu'à $x1 = 400$. -] +Les données cycliques sont un bon exemple où changer la base est utile : vous voulez que le prédicteur corresponde à la fin. + +Nous allons utiliser une série temporelle de données climatiques, divisées en mesures mensuelles, afin de déterminer s'il y a une tendance de température annuelle. + +```{r, eval = FALSE} +data(nottem) # Nottingham temperature time series +n_years <- length(nottem)/12 +nottem_month <- rep(1:12, times = n_years) +nottem_year <- rep(1920:(1920 + n_years - 1), each = 12) +qplot(nottem_month, nottem, colour = factor(nottem_year), geom = "line") + + theme_bw() +``` - --- -# Visualiser la tendance au fil du temps +# Exemple avec des données cycliques -Enfin, pour nous aider à interpréter les résultats, nous pouvons re-transformer l'effet sur une échelle de proportions avec la fonction `itsadug::plot_smooth()` : +```{r, echo = F, fig.height=8, fig.width=12} +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) +qplot(nottem_month, nottem, colour = factor(nottem_year), geom = "line") + + theme_bw() +``` -```{r, echo = -1, fig.height=3.5, fig.width==4} -par(mar=c(3.8,4,0,0)) -plot_smooth(prop_model, view = "x1", main = "", - transform = plogis, ylim = c(0,1), print.summary = F) -abline(h = 0.5, v = diff$start, col = 'red', lty = 2) +--- +# Exemple avec des données cycliques + +Nous pouvons modéliser le changement cyclique de température à travers les mois et la tendance non-linéaire à travers les années, en utilisant une spline cubique, ou `cc` pour modéliser les effets de mois ainsi qu'un terme non-linéaire pour la variable année. + +```{r, fig.height=4, fig.width=8} +year_gam <- gam(nottem ~ s(nottem_year) + s(nottem_month, bs = "cc")) +summary(year_gam)$s.table ``` -Comme précédemment, la proportion de succès augmente au-dessus de 0.5 à $x1 = 400$. +--- +# Exemple avec des données cycliques - +```{r, fig.height=5, fig.width=8} +plot(year_gam, page = 1, scale = 0) +``` +Il y a une hausse d'environ 1-1.5 degrés au cours de la série, mais au cours d'une année, il y a une variation d'environ 20 degrés. Les données réelles varient autour de ces valeurs prédites et ceci représente donc la variance inexpliquée. --- @@ -1656,3 +1749,167 @@ class: inverse, center, bottom
![:scale 50%](images/qcbs_logo.png) + +--- +class: inverse, center, middle + +# Exemple supplémentaire avec d'autres distributions + +--- +# GAM avec d'autres distributions + +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). + +Nous allons utiliser un exemple de données où une répartition binomiale sera nécessaire; 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. + +```{r} +gam_data3 <- read.csv("data/other_dist.csv") +str(gam_data3) +``` + + + +--- +# GAM avec d'autres distributions + +```{r, fig.height=4, fig.width=4} +plot(range(gam_data3$x1), c(0,1), type = "n", + main = "Probabilités de succès dans le temps", + ylab = "Probabilité", xlab = "x1 (temps)") +abline(h = 0.5) + +avg <- aggregate(prop ~ x1, data=gam_data3, mean) +lines(avg$x1, avg$prop, col = "orange", lwd = 2) +``` + +--- +# GAM avec d'autres distributions + +Nous allons tester si cette tendance est linéaire ou non avec un GAM logistique (nous utilisons une famille de distributions binomiales parce que nous avons des données de proportion). + +```{r, warning=F} +prop_model <- gam(prop ~ s(x1), data = gam_data3, weights = total, family = "binomial") +prop_summary <- summary(prop_model) +``` + + + +-- + +.comment[Qu'est ce que représente l'intercepte dans ce modèle?] + +```{r} +prop_summary$p.table +``` + +-- + +.comment[Qu'est ce que le terme de lissage indique?] + +```{r} +prop_summary$s.table +``` + +--- +# GAM avec d'autres distributions + +```{r, echo = FALSE} +prop_summary$p.table +``` + +.comment[Que représente l'intercepte dans ce modèle?] + +**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 : + +.small[ +- Si succès = échecs, le rapport = 1 et le logit est de 0 (log(1) = 0). +- Si succès > échecs, le rapport > 1 et le logit a une valeur positive (log(2) = 0.69). +- Si succès < échecs, le rapport < 1 et le logit a une valeur négative (log(.5) = -0.69). +] + +-- + +> Ici, l'estimé est positif ce qui signifie, qu'en moyenne, il y a plus de succès que d'échecs. + +--- +# GAM avec d'autres distributions + +```{r, echo = FALSE} +prop_summary$s.table +``` + +.comment[Qu'est ce que le terme de lissage indique?] + +Cela représente comment le ratio de succès vs échecs change sur l'échelle de $x1$. + +-- + +> Puisque les **EDF** > 1, la proportion des succès augmente plus rapidement avec $x1$ + +```{r, fig.height=4.3, fig.width=4.3, eval = FALSE} +plot(prop_model) +``` + +```{r, fig.height=4.3, fig.width=4.3, echo = FALSE} +par(mar = c(4,4,0,0)) +plot(prop_model) +``` + +--- +# Visualiser la tendance au fil du temps + +Il y a différente façon de représenter cette relation graphiquement : + +- **Contribution/effet partiel** correspond aux effets isolés d'une interaction ou prédiction particulière. Si vous visualisez votre modèle GAM avec `plot()`, vous obtenez les effets partiels. +- **effets additionnés** correspond aux mesures réponse prédites pour une valeur ou niveau donné de prédicteurs. Si vous visualisez votre GAM avec `itsadug::plot_smooth()`, vous obtenez les effets additionnés. + + +--- +# Visualiser la tendance au fil du temps + +Que nous disent ces graphes sur les succès et échecs? + +```{r, echo = FALSE, fig.height = 4.5, fig.width = 9} +library(itsadug) +par(mfrow=c(1,2), mar = c(4,4,0,0)) +plot(prop_model, select = 1, scale = 0, shade=TRUE) +abline(h=0) + +out <- plot_smooth(prop_model, view = "x1",main = "", print.summary = F) +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, "succès > échecs", col = 'red', font = 2) +``` + + + +.pull-left[ +**Contribution / effets partiels** + +La valeur logit augmente, donc les succès augmentent et les échecs diminuent.] + +.pull-right[ +**Valeurs ajustées, effets additionnés, intercepte inclu** + +Quantités égales de succès et d'échecs jusqu'à $x1 = 400$. +] + + +--- +# Visualiser la tendance au fil du temps + +Enfin, pour nous aider à interpréter les résultats, nous pouvons re-transformer l'effet sur une échelle de proportions avec la fonction `itsadug::plot_smooth()` : + +```{r, echo = -1, fig.height=3.5, fig.width==4} +par(mar=c(3.8,4,0,0)) +plot_smooth(prop_model, view = "x1", main = "", + transform = plogis, ylim = c(0,1), print.summary = F) +abline(h = 0.5, v = diff$start, col = 'red', lty = 2) +``` + +Comme précédemment, la proportion de succès augmente au-dessus de 0.5 à $x1 = 400$. + + + From 0ea20b83fb91a78531df5682347d0d1ed8e6ebff Mon Sep 17 00:00:00 2001 From: Daniel Schoenig Date: Wed, 10 Mar 2021 19:05:32 -0500 Subject: [PATCH 2/4] switch to method='REML' --- pres-en/workshop08-pres-en.Rmd | 74 +++++++++++++++++++++++----------- pres-fr/workshop08-pres-fr.Rmd | 74 +++++++++++++++++++++++----------- 2 files changed, 100 insertions(+), 48 deletions(-) diff --git a/pres-en/workshop08-pres-en.Rmd b/pres-en/workshop08-pres-en.Rmd index 35bd929..35a1f51 100644 --- a/pres-en/workshop08-pres-en.Rmd +++ b/pres-en/workshop08-pres-en.Rmd @@ -82,6 +82,7 @@ install.packages(c('ggplot2', 'itsadug', 'mgcv')) 7. Quick intro to GAMM --- + # Learning objectives 1. Use the mgcv package to fit non-linear relationships, @@ -291,7 +292,7 @@ data_plot <- data_plot + data_plot ``` -Note: as opposed to one fixed coefficient, $\beta$ in linear regression, the smooth function can continually change over the range of the predictor $x$ +Note: as opposed to one fixed coefficient, $\beta$ in linear regression, the smooth function can continually change over the range of the predictor $x$. --- exclude:true @@ -563,9 +564,10 @@ Ideally, it would be good to choose $\lambda$ so that the predicted $\hat{f}$ is $$M = 1/n \times \sum_{i=1}^n (\hat{f_i} - f_i)^2$$ -Since $f$ is unknown, $M$ is estimated using a generalized cross validation technique that leaves out each datum from the data in turn and considers the average ability of models fitted to the remaining data to predict the left out datum. +Since $f$ is unknown, $M$ must be estimated. The recommend methods for this are maximum likelihood (*ML*) or restricted maximum likelihood estimation (*REML*). Generalized cross validation (*GCV*) is another possibility. --- +exclude:true # Principle behind cross validation .center[ @@ -574,17 +576,18 @@ Since $f$ is unknown, $M$ is estimated using a generalized cross validation tech 1. fits many of the data poorly and does no better with the missing point. --- + 2. fits the underlying signal quite well, smoothing through the noise and the missing datum is reasonably well predicted. --- + 3. fits the noise as well as the signal and the extra variability induced causes it to predict the missing datum rather poorly. --- +exclude:true # Principle behind cross validation .center[ @@ -597,6 +600,7 @@ Since $f$ is unknown, $M$ is estimated using a generalized cross validation tech class: inverse, center, middle ## 3. GAM with multiple smooth terms + --- # GAM with multiple variables @@ -623,7 +627,7 @@ isit$Season <- as.factor(isit$Season) ```{r} isit$Season <- as.factor(isit$Season) -basic_model <- gam(Sources ~ Season + s(SampleDepth), data = isit) +basic_model <- gam(Sources ~ Season + s(SampleDepth), data = isit, method = "REML") basic_summary <- summary(basic_model) ``` @@ -727,7 +731,8 @@ We will talk about choosing $k$ in section 5. We can add a second term, `RelativeDepth`, but specify a linear relationship with `Sources` ```{r} -two_term_model <- gam(Sources ~ Season + s(SampleDepth) + RelativeDepth, data = isit) +two_term_model <- gam(Sources ~ Season + s(SampleDepth) + RelativeDepth, + data = isit, method = "REML") two_term_summary <- summary(two_term_model) ``` @@ -759,7 +764,8 @@ plot(two_term_model, page = 1, all.terms = TRUE) We can also explore whether the relationship between `Sources` and `RelativeDepth` is non-linear ```{r} -two_smooth_model <- gam(Sources ~ Season + s(SampleDepth) + s(RelativeDepth), data = isit) +two_smooth_model <- gam(Sources ~ Season + s(SampleDepth) + s(RelativeDepth), + data = isit, method = "REML") two_smooth_summary <- summary(two_smooth_model) ``` @@ -875,8 +881,16 @@ exclude:true
```{r,eval=FALSE} -three_term_model <- gam(Sources ~ Season + s(SampleDepth) + s(RelativeDepth) + Latitude, data = isit) -three_smooth_model <- gam(Sources ~ Season + s(SampleDepth) + s(RelativeDepth) + s(Latitude), data = isit) +three_term_model <- gam(Sources ~ + Season + s(SampleDepth) + s(RelativeDepth) + + Latitude, + data = isit, method = "REML") + +three_smooth_model <- gam(Sources ~ + Season + s(SampleDepth) + s(RelativeDepth) + + s(Latitude), + data = isit, method = "REML") + three_smooth_summary <- summary(three_smooth_model) three_smooth_summary ``` @@ -884,8 +898,16 @@ three_smooth_summary # Challenge 2 - Solution ![:cube]() ```{r,eval=TRUE,echo=FALSE} -three_term_model <- gam(Sources ~ Season + s(SampleDepth) + s(RelativeDepth) + Latitude, data = isit) -three_smooth_model <- gam(Sources ~ Season + s(SampleDepth) + s(RelativeDepth) + s(Latitude), data = isit) +three_term_model <- gam(Sources ~ + Season + s(SampleDepth) + s(RelativeDepth) + + Latitude, + data = isit, method = "REML") + +three_smooth_model <- gam(Sources ~ + Season + s(SampleDepth) + s(RelativeDepth) + + s(Latitude), + data = isit, method = "REML") + three_smooth_summary <- summary(three_smooth_model) three_smooth_summary ``` @@ -1000,7 +1022,7 @@ We will examine interaction effect using our categorical variable `Season` and a factor_interact <- gam(Sources ~ Season + s(SampleDepth, by=Season) + s(RelativeDepth), - data = isit) + data = isit, method = "REML") summary(factor_interact)$s.table ``` @@ -1051,7 +1073,8 @@ From the plots, we saw that the shape of the smooth terms were comparable among Finally we'll look at the interactions between 2 smoothed terms, `SampleDepth` and `RelativeDepth`. ```{r} -smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth), data = isit) +smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth), + data = isit, method = "REML") summary(smooth_interact)$s.table ``` @@ -1142,7 +1165,8 @@ Just like generalized linear models (GLM), we can formulate **generalized** addi Lets recall the interaction model for the bioluminescence data: ```{r} -smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth), data = isit) +smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth), + data = isit, method = "REML") summary(smooth_interact)$p.table @@ -1195,7 +1219,7 @@ Let's refit the model with a larger `k`: ```{r} smooth_interact_k30 <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), - data = isit) + data = isit, method = "REML") summary(smooth_interact_k30)$p.table summary(smooth_interact_k30)$s.table ``` @@ -1297,11 +1321,11 @@ To get an overview of families available in `mgcv`: ```{r} # Normal distribution smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), - data = isit) + data = isit, method = "REML") # Tweedie family with log link smooth_interact_tw <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), family = tw(link = "log"), - data = isit) + data = isit, method = "REML") ``` @@ -1314,7 +1338,7 @@ Fit the model: ```{r} smooth_interact_tw <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), family = tw(link = "log"), - data = isit) + data = isit, method = "REML") summary(smooth_interact_tw)$p.table summary(smooth_interact_tw)$s.table ``` @@ -1442,7 +1466,7 @@ qplot(nottem_month, nottem, colour = factor(nottem_year), geom = "line") + We can model both the cyclic change of temperature across months and the non-linear trend through years, using a cyclical cubic spline, or `cc`, for the month variable and a regular smooth for the year variable. ```{r, fig.height=4, fig.width=8} -year_gam <- gam(nottem ~ s(nottem_year) + s(nottem_month, bs = "cc")) +year_gam <- gam(nottem ~ s(nottem_year) + s(nottem_month, bs = "cc"), method = "REML") summary(year_gam)$s.table ``` @@ -1560,7 +1584,7 @@ str(gam_data2) ```{r, fig.height=5, fig.width=5.5, echo = -1} par(mar=c(4,4,1,1)) -gamm_intercept <- gam(y ~ s(x0) + s(fac, bs = "re"), data = gam_data2) +gamm_intercept <- gam(y ~ s(x0) + s(fac, bs = "re"), data = gam_data2, method = "REML") summary(gamm_intercept)$s.table plot(gamm_intercept, select = 2) ``` @@ -1611,7 +1635,7 @@ plot_smooth(gamm_intercept, view ="x0", cond = list(fac = "4"), add = T, col = ' # GAMM with a random slope ```{r} -gamm_slope <- gam(y ~ s(x0) + s(x0, fac, bs = "re"), data = gam_data2) +gamm_slope <- gam(y ~ s(x0) + s(x0, fac, bs = "re"), data = gam_data2, method = "REML") summary(gamm_slope)$s.table ``` @@ -1659,7 +1683,7 @@ plot_smooth(gamm_slope, view = "x0", cond = list(fac = "4"), add = T, col = 'tur ```{r} gamm_int_slope <- gam(y ~ s(x0) + s(fac, bs = "re") + s(fac, x0, bs = "re"), - data = gam_data2) + data = gam_data2, method = "REML") summary(gamm_int_slope)$s.table ``` @@ -1713,7 +1737,8 @@ plot(gamm_int_slope, select = 3) # GAMM with a random smooth ```{r} -gamm_smooth <- gam(y ~ s(x0) + s(x0, fac, bs = "fs", m = 1), data = gam_data2) +gamm_smooth <- gam(y ~ s(x0) + s(x0, fac, bs = "fs", m = 1), + data = gam_data2, method = "REML") summary(gamm_smooth)$s.table ``` @@ -1843,7 +1868,8 @@ lines(avg$x1, avg$prop, col = "orange", lwd = 2) We will test if this trend is linear or not using a logistic GAM (we use a binomial family distribution given that our response is proportion data). ```{r, warning=F} -prop_model <- gam(prop ~ s(x1), data = gam_data3, weights = total, family = "binomial") +prop_model <- gam(prop ~ s(x1), data = gam_data3, weights = total, + family = "binomial", method = "REML") prop_summary <- summary(prop_model) ``` diff --git a/pres-fr/workshop08-pres-fr.Rmd b/pres-fr/workshop08-pres-fr.Rmd index 2af16ac..81dc0ff 100644 --- a/pres-fr/workshop08-pres-fr.Rmd +++ b/pres-fr/workshop08-pres-fr.Rmd @@ -294,7 +294,7 @@ data_plot <- data_plot + data_plot ``` -.comment[Note: contrairement à un coefficient fixe $\beta$, la fonction de lissage peut changer tout au long du gradient $x$] +Note: contrairement à un coefficient fixe $\beta$, la fonction de lissage peut changer tout au long du gradient $x$. --- exclude:true @@ -571,9 +571,10 @@ Si $λ$ est trop élevé, les données seront trop lissées et si elle est trop $$M = 1/n \times \sum_{i=1}^n (\hat{f_i} - f_i)^2$$ -Étant donné que $f$ est inconnue, $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é. +Étant donné que $f$ est inconnue, $M$ doit être estimé. Les méthodes recommandées pour ce faire sont le maximum de vraisemblance (maximum likelihood, *ML*) ou l'estimation par maximum de vraisemblance restreint (restricted maximum likelihood, *REML*). La validation croisée généralisée (*GCV*) est une autre possibilité. --- +exclude:true # Principe de validation croisée .center[ @@ -582,15 +583,16 @@ $$M = 1/n \times \sum_{i=1}^n (\hat{f_i} - f_i)^2$$ 1. ajustement faible par rapport aux données et ne fait pas mieux avec le point manquant. --- + 2. très bon ajustement de la courbe du signal sous-jacent, le lissage passe à travers le bruit et la donnée manquante est plutôt bien prédite. --- + 3. la courbe ajuste le bruit aussi bien que le signal, la variabilité supplémentaire amène à prédire la donnée manquante plutôt mal. --- +exclude:true # Principe de validation croisée .center[ @@ -628,7 +630,7 @@ isit$Season <- as.factor(isit$Season) ```{r} isit$Season <- as.factor(isit$Season) -basic_model <- gam(Sources ~ Season + s(SampleDepth), data = isit) +basic_model <- gam(Sources ~ Season + s(SampleDepth), data = isit, method = "REML") basic_summary <- summary(basic_model) ``` @@ -722,6 +724,8 @@ La limite supérieure d'**EDF** est déterminée par les dimensions de base $k$ En pratique, le choix exact de $k$ est arbitraire, mais il devrait être **suffisamment grand** pour permettre une fonction lisse suffisamment complexe. +Nous discuterons du choix de $k$ dans la section 5. + --- # GAM à plusieurs variables @@ -729,7 +733,8 @@ En pratique, le choix exact de $k$ est arbitraire, mais il devrait être **suffi Nous pouvons ajouter un second terme, `RelativeDepth`, mais spécifier une relation linéaire avec `Sources` ```{r} -two_term_model <- gam(Sources ~ Season + s(SampleDepth) + RelativeDepth, data = isit) +two_term_model <- gam(Sources ~ Season + s(SampleDepth) + RelativeDepth, + data = isit, method = "REML") two_term_summary <- summary(two_term_model) ``` @@ -762,7 +767,8 @@ plot(two_term_model, page=1, all.terms = T) Nous pouvons aussi vérifier que la relation entre `Sources` et `RelativeDepth` est non-linéaire. ```{r} -two_smooth_model <- gam(Sources ~ Season + s(SampleDepth) + s(RelativeDepth), data = isit) +two_smooth_model <- gam(Sources ~ Season + s(SampleDepth) + s(RelativeDepth), + data = isit, method = "REML") two_smooth_summary <- summary(two_smooth_model) ``` @@ -881,8 +887,16 @@ exclude:true --- # Défi 2 - Solution ![:cube]() ```{r,eval=FALSE} -three_term_model <- gam(Sources ~ Season + s(SampleDepth) + s(RelativeDepth) + Latitude, data = isit) -three_smooth_model <- gam(Sources ~ Season + s(SampleDepth) + s(RelativeDepth) + s(Latitude), data = isit) +three_term_model <- gam(Sources ~ + Season + s(SampleDepth) + s(RelativeDepth) + + Latitude, + data = isit, method = "REML") + +three_smooth_model <- gam(Sources ~ + Season + s(SampleDepth) + s(RelativeDepth) + + s(Latitude), + data = isit, method = "REML") + three_smooth_summary <- summary(three_smooth_model) three_smooth_summary ``` @@ -890,8 +904,16 @@ three_smooth_summary --- # Défi 2 - Solution ![:cube]() ```{r,eval=TRUE,echo=FALSE} -three_term_model <- gam(Sources ~ Season + s(SampleDepth) + s(RelativeDepth) + Latitude, data = isit) -three_smooth_model <- gam(Sources ~ Season + s(SampleDepth) + s(RelativeDepth) + s(Latitude), data = isit) +three_term_model <- gam(Sources ~ + Season + s(SampleDepth) + s(RelativeDepth) + + Latitude, + data = isit, method = "REML") + +three_smooth_model <- gam(Sources ~ + Season + s(SampleDepth) + s(RelativeDepth) + + s(Latitude), + data = isit, method = "REML") + three_smooth_summary <- summary(three_smooth_model) three_smooth_summary ``` @@ -956,7 +978,7 @@ Nous allons examiner l'effet de l'interaction en utilisant notre variable qualit factor_interact <- gam(Sources ~ Season + s(SampleDepth,by=Season) + s(RelativeDepth), - data = isit) + data = isit, method = "REML") summary(factor_interact)$s.table ``` @@ -1007,7 +1029,8 @@ AIC(two_smooth_model, factor_interact) Finalement, nous regardons les interactions entre 2 termes non linéaires, `SampleDepth` et `RelativeDepth`. ```{r} -smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth), data = isit) +smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth), + data = isit, method = "REML") summary(smooth_interact)$s.table ``` @@ -1081,7 +1104,8 @@ Tout comme les modèles linéaires généralisés (GLM), nous pouvons formuler d Rappelons le modèle d'interaction pour les données de bioluminescence : ```{r} -smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth), data = isit) +smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth), + data = isit, method = "REML") summary(smooth_interact)$p.table @@ -1135,7 +1159,7 @@ Refaisons le modèle avec un `k` plus grand : ```{r} smooth_interact_k30 <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), - data = isit) + data = isit, method = "REML") summary(smooth_interact_k30)$p.table summary(smooth_interact_k30)$s.table ``` @@ -1236,11 +1260,11 @@ Pour en savoir plus sur les familles disponibles dans `mgcv` : ```{r} # Distribution normale smooth_interact <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), - data = isit) + data = isit, method = "REML") # Tweedie avec log comme fonction de liaison smooth_interact_tw <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), family = tw(link = "log"), - data = isit) + data = isit, method = "REML") ``` @@ -1253,7 +1277,7 @@ Ajuster le modèle : ```{r} smooth_interact_tw <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), family = tw(link = "log"), - data = isit) + data = isit, method = "REML") summary(smooth_interact_tw)$p.table summary(smooth_interact_tw)$s.table ``` @@ -1380,7 +1404,7 @@ qplot(nottem_month, nottem, colour = factor(nottem_year), geom = "line") + Nous pouvons modéliser le changement cyclique de température à travers les mois et la tendance non-linéaire à travers les années, en utilisant une spline cubique, ou `cc` pour modéliser les effets de mois ainsi qu'un terme non-linéaire pour la variable année. ```{r, fig.height=4, fig.width=8} -year_gam <- gam(nottem ~ s(nottem_year) + s(nottem_month, bs = "cc")) +year_gam <- gam(nottem ~ s(nottem_year) + s(nottem_month, bs = "cc"), method = "REML") summary(year_gam)$s.table ``` @@ -1497,7 +1521,7 @@ str(gam_data2) ```{r, fig.height=4.5, fig.width=5.5, echo = -1} par(mar=c(4,4,1,1)) -gamm_intercept <- gam(y ~ s(x0) + s(fac, bs = "re"), data = gam_data2) +gamm_intercept <- gam(y ~ s(x0) + s(fac, bs = "re"), data = gam_data2, method = "REML") summary(gamm_intercept)$s.table plot(gamm_intercept, select = 2) ``` @@ -1548,7 +1572,7 @@ plot_smooth(gamm_intercept, view = "x0", cond = list(fac = "4"), add = T, col = # GAMM avec une pente aléatoire ```{r} -gamm_slope <- gam(y ~ s(x0) + s(x0, fac, bs = "re"), data = gam_data2) +gamm_slope <- gam(y ~ s(x0) + s(x0, fac, bs = "re"), data = gam_data2, method = "REML") summary(gamm_slope)$s.table ``` @@ -1596,7 +1620,7 @@ plot_smooth(gamm_slope, view = "x0", cond = list(fac = "4"), add = T, col = 'tur ```{r} gamm_int_slope <- gam(y ~ s(x0) + s(fac, bs = "re") + s(fac, x0, bs = "re"), - data = gam_data2) + data = gam_data2, method = "REML") summary(gamm_int_slope)$s.table ``` @@ -1650,7 +1674,8 @@ plot(gamm_int_slope, select = 3) # GAMM avec une surface lisse aléatoire ```{r} -gamm_smooth <- gam(y ~ s(x0) + s(x0, fac, bs = "fs", m = 1), data = gam_data2) +gamm_smooth <- gam(y ~ s(x0) + s(x0, fac, bs = "fs", m = 1), + data = gam_data2, method = "REML") summary(gamm_smooth)$s.table ``` @@ -1788,7 +1813,8 @@ lines(avg$x1, avg$prop, col = "orange", lwd = 2) Nous allons tester si cette tendance est linéaire ou non avec un GAM logistique (nous utilisons une famille de distributions binomiales parce que nous avons des données de proportion). ```{r, warning=F} -prop_model <- gam(prop ~ s(x1), data = gam_data3, weights = total, family = "binomial") +prop_model <- gam(prop ~ s(x1), data = gam_data3, weights = total, + family = "binomial", method = "REML") prop_summary <- summary(prop_model) ``` From 5ca945c55606efdf2494f5ecc528f6727380d5e6 Mon Sep 17 00:00:00 2001 From: Daniel Schoenig Date: Wed, 10 Mar 2021 19:36:04 -0500 Subject: [PATCH 3/4] typos --- pres-en/workshop08-pres-en.Rmd | 4 ++-- pres-fr/workshop08-pres-fr.Rmd | 13 +++++++------ 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/pres-en/workshop08-pres-en.Rmd b/pres-en/workshop08-pres-en.Rmd index 35a1f51..91dff4c 100644 --- a/pres-en/workshop08-pres-en.Rmd +++ b/pres-en/workshop08-pres-en.Rmd @@ -1118,7 +1118,7 @@ vis.gam(smooth_interact, view = c("SampleDepth", "RelativeDepth"), AIC(two_smooth_model, smooth_interact) ``` -The model with the interaction between `s(SampleDepth)` and `s(RelativeDepth)` have a lower AIC and the 2D plot nicely illustrates this non-linear interactions, where `Sources` is lower at high values of `SampleDepth` but mid to high values of `RelativeDepth`. +The model with the interaction between `s(SampleDepth)` and `s(RelativeDepth)` has a lower AIC and the plots nicely illustrate the non-linear interaction, where `Sources` is lower at high values of `SampleDepth` but increases with `RelativeDepth`. @@ -1155,7 +1155,7 @@ But what can we do if: .comment[These cases are very common!] -Just like generalized linear models (GLM), we can formulate **generalized** additive models to deal with these issues +Just like generalized linear models (GLM), we can formulate **generalized** additive models to deal with these issues. --- diff --git a/pres-fr/workshop08-pres-fr.Rmd b/pres-fr/workshop08-pres-fr.Rmd index 81dc0ff..e357458 100644 --- a/pres-fr/workshop08-pres-fr.Rmd +++ b/pres-fr/workshop08-pres-fr.Rmd @@ -1071,9 +1071,10 @@ class: inverse, center, middle Le modèle additif de base peut être étendu de plusieurs façons : -1. Utilisation de différents types d'effets aléatoires pour ajuster des modèles à effets mixtes. +1. Utiliser d'autres distributions pour la variable de réponse avec l'argument `family` (comme dans un GLM), 2. Utiliser de différents types de fonctions de base, -3. Utiliser d'autres distributions pour la variable de réponse avec l'argument `famille` (comme dans un GLM), +3. Utilisation de différents types d'effets aléatoires pour ajuster des modèles à effets mixtes. + Nous allons maintenant examiner ces aspects. @@ -1141,7 +1142,7 @@ Fonctions incluses dans `mgcv` : Avons-nous choisi `k` assez grand ? -.comment [Le défaut pour les interactions lisses est `k = 30`]. +.comment[Le défaut pour les interactions lisses est `k = 30`]. ```{r} k.check(smooth_interact) @@ -1193,7 +1194,7 @@ par(mfrow = c(2,2)) # Afficher les 4 tracés à la fois gam.check(smooth_interact) ``` -.comment [En plus des tracés, `gam.check()` fournit également la sortie de `k.check()`]. +.comment[En plus des graphiques, `gam.check()` fournit également la sortie de `k.check()`]. --- @@ -1291,9 +1292,9 @@ k.check(smooth_interact_tw) --- -# Challenge 3 - Solution ![:cube]() +# Défi 3 - Solution ![:cube]() -Tracé des résidus : +Tracés des résidus : ```{r, eval = FALSE} par(mfrow=c(2,2)) From 9c29b5d8eea9341e66d6f56aea6dcc216d68905a Mon Sep 17 00:00:00 2001 From: Laurie Maynard Date: Thu, 11 Mar 2021 08:37:41 -0400 Subject: [PATCH 4/4] Last check Fixed a mistake in the script, with the wrong object being called. Clarified challenge 2. Fixed a typo. changed an object name to better fit the material presented. --- pres-en/workshop08-pres-en.Rmd | 18 +++++++++--------- pres-fr/workshop08-pres-fr.Rmd | 16 ++++++++-------- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/pres-en/workshop08-pres-en.Rmd b/pres-en/workshop08-pres-en.Rmd index 91dff4c..1156cc8 100644 --- a/pres-en/workshop08-pres-en.Rmd +++ b/pres-en/workshop08-pres-en.Rmd @@ -772,13 +772,13 @@ two_smooth_summary <- summary(two_smooth_model) Information on parametric effects (linear terms): ```{r} -two_term_summary$p.table +two_smooth_summary$p.table ``` Information on additive effects (non-linear terms): ```{r} -two_term_summary$s.table +two_smooth_summary$s.table ``` --- @@ -865,7 +865,7 @@ anova(basic_model, two_term_model, two_smooth_model,test='Chisq')
-1. Create 2 new models, with `Latitude` as a linear and smoothed term. +1. Create 2 new models, by adding `Latitude` to `two_smooth_model` as a linear and smoothed term. 2. Determine if `Latitude` is an important term to include using plots, coefficient tables and the AIC function. --- @@ -1179,7 +1179,7 @@ summary(smooth_interact)$s.table As for a GLM, it is essential to check whether we correctly specifified the model, especially the *distribution* for the response. -**We have two verify:** +**We have to verify:** 1. The choice of basis dimension `k`. 2. The residuals plots (just as for a GLM). @@ -1218,10 +1218,10 @@ The **EDF are very close to** `k`, this could be problematic. Let's refit the model with a larger `k`: ```{r} -smooth_interact_k30 <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), +smooth_interact_k60 <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), data = isit, method = "REML") -summary(smooth_interact_k30)$p.table -summary(smooth_interact_k30)$s.table +summary(smooth_interact_k60)$p.table +summary(smooth_interact_k60)$s.table ``` -- @@ -1229,7 +1229,7 @@ summary(smooth_interact_k30)$s.table Is `k` large enough this time? ```{r} -k.check(smooth_interact_k30) +k.check(smooth_interact_k60) ``` -- @@ -1237,7 +1237,7 @@ k.check(smooth_interact_k30) Looks better, let's replace the old model: ```{r} -smooth_interact <- smooth_interact_k30 +smooth_interact <- smooth_interact_k60 ``` --- diff --git a/pres-fr/workshop08-pres-fr.Rmd b/pres-fr/workshop08-pres-fr.Rmd index e357458..b593d42 100644 --- a/pres-fr/workshop08-pres-fr.Rmd +++ b/pres-fr/workshop08-pres-fr.Rmd @@ -775,13 +775,13 @@ two_smooth_summary <- summary(two_smooth_model) Informations sur les effets paramétriques (termes linéaires) : ```{r} -two_term_summary$p.table +two_smooth_summary$p.table ``` Informations sur les effets additifs (termes non linéaires) : ```{r} -two_term_summary$s.table +two_smooth_summary$s.table ``` @@ -872,7 +872,7 @@ AIC(basic_model, two_term_model, two_smooth_model)
-1. Créez deux nouveaux modèles avec la variable `Latitude` comme paramètre linéaire et non linéaire. +1. Créez deux nouveaux modèles en ajoutant la variable `Latitude` comme paramètre linéaire et non linéaire au modèle précédent. 2. Utilisez des graphiques, les tables des coefficients et la fonction `AIC()` afin de déterminer s'il est nécessaire d'inclure `Latitude` dans le modèle. --- @@ -1159,10 +1159,10 @@ Les **EDF se rapprochent beaucoup de** `k`, cela pourrait être problématique. Refaisons le modèle avec un `k` plus grand : ```{r} -smooth_interact_k30 <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), +smooth_interact_k60 <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60), data = isit, method = "REML") -summary(smooth_interact_k30)$p.table -summary(smooth_interact_k30)$s.table +summary(smooth_interact_k60)$p.table +summary(smooth_interact_k60)$s.table ``` -- @@ -1170,7 +1170,7 @@ summary(smooth_interact_k30)$s.table Est-ce que `k` est assez grand ? ```{r} -k.check(smooth_interact_k30) +k.check(smooth_interact_k60) ``` -- @@ -1178,7 +1178,7 @@ k.check(smooth_interact_k30) C'est mieux, remplaçons l'ancien modèle : ```{r} -smooth_interact <- smooth_interact_k30 +smooth_interact <- smooth_interact_k60 ``` ---