---
title: "Getting started with gratia"
output: 
  rmarkdown::html_vignette:
    fig_width: 8
    fig_height: 5.3333
    dev: "png"
bibliography: getting-started.bib
vignette: >
  %\VignetteIndexEntry{Getting started with gratia}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup-knitr, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library("gratia")
library("mgcv")
```

gratia is a package to make working with generalized additive models (GAMs) in R
easier, including producing plots of estimated smooths using the ggplot2
📦.

This introduction will cover some of the basic functionality of gratia to get
you started. We'll work with some classic simulated data often used to
illustrate properties of GAMs

```{r data-sim}
df <- data_sim("eg1", seed = 42)
df
```

and the following GAM

```{r fit-gam}
m <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = df, method = "REML")
summary(m)
```

## Plotting

gratia provides the `draw()` function to produce plots using the ggplot2
&#128230;. To plot the estimated smooths from the GAM we fitted above, use

```{r draw-gam}
draw(m)
```

The plots produced are *partial effect* plots, which show the component
contributions, on the link scale, of each model term to the linear predictor.
The y axis on these plots is typically centred around 0 due to most smooths
having a sum-to-zero identifiability constraint applied to them. This
constraint is what allows the model to include multiple smooths and remain
identifiable. These plots allow you to read off the contributions of each
smooth to the fitted response (on the link scale); they show link-scale
predictions of the response for each smooth, conditional upon all other terms
in the model, including any parametric effects and the intercept, having zero 
contribution. In the parlance of the marginaleffects package
[@Arel-BundockGreiferHeiss:2024], these plots show adjusted predictions, just
where the adjustment includes setting the contribution of all other model terms
to the predicted value to zero. For partial derivatives (what *marginaleffects*
would call a marginal effect or slope), gratia provides `derivatives()`.

The resulting plot is intended as reasonable overview of the estimated model,
but it offers limited option to modify the resulting plot. If you want full
control, you can obtain the data used to create the plot above with 
`smooth_estimates()`
```{r smooth-estimates}
sm <- smooth_estimates(m)
sm
```
which will evaluate all smooths at values that are evenly spaced over the range
of the covariate(s). If you want to evaluate only selected smooths, you can
specify which via the `smooth` argument. This takes the *smooth labels* which
are the names of the smooths as they are known to mgcv. To list the labels for
the smooths in use
```{r smooths}
smooths(m)
```
To evaluate only $f(x_2)$ use
```{r smooth-estimates-x2}
sm <- smooth_estimates(m, smooth = "s(x2)")
sm
```
Then you can generate your own plot using the ggplot2 package, for example
```{r ggplot-smooth}
library("ggplot2")
library("dplyr")
sm |>
  add_confint() |>
  ggplot(aes(y = .estimate, x = x2)) +
  geom_ribbon(aes(ymin = .lower_ci, ymax = .upper_ci),
    alpha = 0.2, fill = "forestgreen"
  ) +
  geom_line(colour = "forestgreen", linewidth = 1.5) +
  labs(
    y = "Partial effect",
    title = expression("Partial effect of" ~ f(x[2])),
    x = expression(x[2])
  )
```

## Model diagnostics

The `appraise()` function provides standard diagnostic plots for GAMs

```{r appraise}
appraise(m)
```

The plots produced are (from left-to-right, top-to-bottom),

* a quantile-quantile (QQ) plot of deviance residuals,
* a scatterplot of deviance residuals against the linear predictor,
* a histogram of deviance residuals, and
* a scatterplot of observed vs fitted values.

Adding partial residuals to the partial effect plots produced by `draw()`
can also help diagnose problems with the model, such as oversmoothing
```{r draw-partial-residuals}
draw(m, residuals = TRUE)
```

## Want to learn more?

*gratia* is in very active development and an area of development that is
currently lacking is documentation. To find out more about the package, look at
the [help pages for the package](https://gavinsimpson.github.io/gratia/reference/index.html)
and look at the examples for more code to help you get going.

# References