diff --git a/pres-en/workshop08-pres-en.Rmd b/pres-en/workshop08-pres-en.Rmd
index fad4499..1156cc8 100644
--- a/pres-en/workshop08-pres-en.Rmd
+++ b/pres-en/workshop08-pres-en.Rmd
@@ -77,11 +77,12 @@ 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
---
+
# Learning objectives
1. Use the mgcv package to fit non-linear relationships,
@@ -291,7 +292,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
@@ -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)
```
@@ -695,9 +699,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 +723,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
@@ -726,7 +731,8 @@ In practice, the exact choice of $k$ is arbitrary, but it should be **large enou
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)
```
@@ -758,20 +764,21 @@ 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)
```
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
```
---
@@ -858,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.
---
@@ -874,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
```
@@ -883,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
```
@@ -974,6 +997,7 @@ AIC(two_smooth_model, three_term_model, test = "Chisq")
---
+
class: inverse, center, middle
## 4. Interactions
@@ -998,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
```
@@ -1049,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
```
@@ -1093,274 +1118,367 @@ 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`.
+
+
---
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()`  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()`  for modeling 2- or n-dimensional interaction surfaces of variables that are not on the same scale. Includes main effects.
+--
-`ti()`  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`  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, method = "REML")
+
+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`  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`  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 to 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_k60 <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60),
+ data = isit, method = "REML")
+summary(smooth_interact_k60)$p.table
+summary(smooth_interact_k60)$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?
-## 6. Other distributions
+```{r}
+k.check(smooth_interact_k60)
+```
+
+--
+Looks better, let's replace the old model:
+
+```{r}
+smooth_interact <- smooth_interact_k60
+```
---
-# 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, method = "REML")
+# Tweedie family with log link
+smooth_interact_tw <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60),
+ family = tw(link = "log"),
+ data = isit, method = "REML")
```
+---
+
+# 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, method = "REML")
+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()`  for modeling a 1-dimensional smooth or for modeling interactions among variables measured using the same unit and the same scale
+`te()`  for modeling 2- or n-dimensional interaction surfaces of variables that are not on the same scale. Includes main effects.
+
+`ti()`  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`  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`  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`  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"), method = "REML")
+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.
-
---
@@ -1466,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)
```
@@ -1517,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
```
@@ -1565,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
```
@@ -1619,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
```
@@ -1709,3 +1828,172 @@ class: inverse, center, bottom
# Thank you for attending this workshop!
+
+---
+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", method = "REML")
+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..b593d42 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
@@ -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,20 +767,21 @@ 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)
```
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
```
@@ -866,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.
---
@@ -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
```
@@ -1039,7 +1062,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 +1071,351 @@ 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),
+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. Utilisation de différents types d'effets aléatoires pour ajuster des modèles à effets mixtes.
+
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()`  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()`  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*.
+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`  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`  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`  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, method = "REML")
-
+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_k60 <- gam(Sources ~ Season + s(SampleDepth, RelativeDepth, k = 60),
+ data = isit, method = "REML")
+summary(smooth_interact_k60)$p.table
+summary(smooth_interact_k60)$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_k60)
```
-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
+--
-## 6. Autres distributions
+C'est mieux, remplaçons l'ancien modèle :
+
+```{r}
+smooth_interact <- smooth_interact_k60
+```
---
-# 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 graphiques, `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, 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, method = "REML")
```
-.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, method = "REML")
+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
+# Défi 3 - Solution ![:cube]()
+
+Tracés 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()`  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()`  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*.
---
-# 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`  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 :
-.pull-left[
-**Contribution / effets partiels**
+`d`  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.
-La valeur logit augmente, donc les succès augmentent et les échecs diminuent.]
+`bs`  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*).
-.pull-right[
-**Valeurs ajustées, effets additionnés, intercepte inclu**
+
-Quantités égales de succès et d'échecs jusqu'à $x1 = 400$.
-]
+---
+# 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.
+
+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"), method = "REML")
+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.
---
@@ -1404,7 +1522,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)
```
@@ -1455,7 +1573,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
```
@@ -1503,7 +1621,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
```
@@ -1557,7 +1675,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
```
@@ -1656,3 +1775,168 @@ class: inverse, center, bottom
+
+---
+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", method = "REML")
+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$.
+
+
+