DESCRIPTION

@@ -1,7 +1,7 @@
 Package: tidymv
 Type: Package
 Title: Tidy Model Visualisation
-Version: 1.2.0
+Version: 1.3.0
 Authors@R: person("Stefano", "Coretta", email = paste0("stefano.coretta", "@", "gmail.com"),
                   role = c("aut", "cre"))
 Description: This package provides functions for model visualisation using tidy
@@ -10,13 +10,15 @@ License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
 RoxygenNote: 6.0.1
-Depends: itsadug
-Imports: tidyverse,
-    cowplot,
+Imports: cowplot,
     dplyr,
-    itsadug,
     ggplot2,
-    magrittr
+    grDevices,
+    itsadug,
+    magrittr,
+    rlang,
+    tidyr,
+    tidyverse
 Suggests: knitr,
     rmarkdown
 VignetteBuilder: knitr

NAMESPACE

@@ -1,5 +1,9 @@
 # Generated by roxygen2: do not edit by hand
 
 export(create_event_start)
+export(get_gam_predictions)
 export(plot_gamsd)
+export(plot_smooths)
 importFrom(magrittr,"%>%")
+importFrom(rlang,":=")
+importFrom(stats,"predict")

NEWS.md

@@ -1,5 +1,22 @@
 # Change Log
 
+## [1.3.0] - 2018-02-26
+### Added
+- `get_gam_predictions()` for predicting with a `gam` object
+- `plot_smooths()` for plotting smooths
+- `plot-smooths.Rmd` vignette
+- support for models with smooths with terms not in `time_series`
+
+### Fixed
+- typo in documentation of `create_event_start()`
+- double call of `itsadug` in DESCRIPTION
+- note about `lag` on check
+- note about non-imported `itsadug`
+- incompatibility with `rlang@v0.2.0`
+
+### Deprecated
+- `plot_gamsd()` will be deprecated: use `plot_smooths` instead. Plotting both smooths and difference smooth is not supported yet.
+
 ## 1.2.0 - 2017-12-07
 ### Added
 - `rm_re` parameter for removing random effects in `plot_gamsd`

R/functions.R

@@ -5,22 +5,193 @@
 #' with tibbles.
 #'
 #' @param tibble A tibble arranged according to the time series.
-#' @param event.col A string with the name of the column that defines the time series.
+#' @param event_col A string with the name of the column that defines the time series.
 #'
 #' @export
 create_event_start <- function(tibble, event_col) {
     dplyr::mutate(
         tibble,
         start.event = ifelse(
-            as.character(tibble[[event_col]]) == lag(as.character(tibble[[event_col]]), default = FALSE),
+            as.character(tibble[[event_col]]) == dplyr::lag(as.character(tibble[[event_col]]), default = FALSE),
             FALSE,
             TRUE
             )
     )
 }
 
+#' Get predictions from a GAM model.
+#'
+#' It returns a tibble with the predictions from a a \link[mgcv]{gam} or \link[mgcv]{bam} object.
+#'
+#' @param model A \code{gam} or \code{bam} model object.
+#' @param time_series An unquoted expression indicating the model term that defines the time series.
+#' @param series_length An integer indicating how many values along the time series to use for predicting the outcome term.
+#' @param conditions A list of quosures with \link[rlang]{quos} specifying the levels to plot from the model terms.
+#' @param exclude_random Whether to exclude random smooths (the default is \code{TRUE}).
+#'
+#' @export
+get_gam_predictions <- function(model, time_series, series_length = 25, conditions = NULL, exclude_random = TRUE) {
+    time_series_q <- dplyr::enquo(time_series)
+    time_series_name <- dplyr::quo_name(time_series_q)
+    outcome_q <- model$formula[[2]]
+
+    fitted <- model$model
+
+    random_effects <- list()
+    random_effects_terms <- NULL
+
+    if (exclude_random == TRUE) {
+        for (i in 1:length(model[["smooth"]])) {
+            smooth_class <- attr(model$smooth[[i]],"class")[1]
+            if (smooth_class %in% c("random.effect", "fs.interaction")) {
+                random_effects <- c(
+                    random_effects,
+                    list(model$smooth[[i]]$label)
+                )
+                random_effects_terms <- c(
+                    random_effects_terms,
+                    model$smooth[[i]]$fterm
+                )
+            }
+        }
+    }
+
+    time_series_min <- dplyr::select(fitted, !!time_series_q) %>% min()
+    time_series_max <- dplyr::select(fitted, !!time_series_q) %>% max()
+
+    fitted <- fitted %>%
+        dplyr::select(-!!time_series_q, -!!outcome_q)
+
+    if ("(AR.start)" %in% colnames(fitted)) {
+        fitted$`(AR.start)` <- NULL
+    }
+
+    fitted_series <- fitted %>%
+        unique()
+
+    fitted_series <- fitted_series %>%
+        dplyr::mutate(
+            !!dplyr::quo_name(time_series_q) := rep(
+                list(seq(time_series_min, time_series_max, length.out = series_length)),
+                nrow(fitted_series)
+            )
+        ) %>%
+        tidyr::unnest(!!time_series_q)
+
+    if (exclude_random) {
+        if (rlang::is_empty(random_effects)) {
+            exclude_random_effects <- as.null()
+        } else {
+            exclude_random_effects <- random_effects
+        }
+    } else {
+        exclude_random_effects <- as.null()
+    }
+
+    # Exclude smooth terms which are not the time series to be plotted
+    exclude_smooths <- as.null()
+    exclude_terms <- as.null()
+    for (smooth in 1:length(model[["smooth"]])) {
+        smooth_term <- model[["smooth"]][[smooth]][["term"]]
+        if (smooth_term != time_series_name) {
+            exclude_terms <- c(exclude_terms, smooth_term)
+            smooth_label <- model[["smooth"]][[smooth]][["label"]]
+            exclude_smooths <- c(exclude_smooths, smooth_label)
+        }
+    }
+
+    exclude_these <- c(exclude_random_effects, exclude_smooths)
+
+    predicted <- stats::predict(
+        model,
+        fitted_series,
+        se.fit = TRUE,
+        exclude = exclude_these
+    )
+
+    predicted_tbl <- cbind(fitted_series, predicted) %>%
+        dplyr::mutate(
+            CI_upper = fit + 1.96 * se.fit,
+            CI_lower = fit - 1.96 * se.fit
+        ) %>%
+        dplyr::rename(
+            !!outcome_q := fit,
+            SE = se.fit
+        )
+
+    if (!is.null(exclude_random_effects)) {
+        predicted_tbl <- predicted_tbl %>%
+            dplyr::select(-(!!!rlang::syms(random_effects_terms))) %>%
+            unique()
+    }
 
-#' Plot GLM estimate smooths and difference curve.
+    if (!is.null(exclude_smooths)) {
+        predicted_tbl <- predicted_tbl %>%
+            dplyr::select(-(!!!rlang::syms(exclude_terms))) %>%
+            unique()
+    }
+
+    if (!is.null(conditions)) {
+        predicted_tbl <- predicted_tbl %>%
+            dplyr::filter(!!!conditions)
+    }
+
+    return(predicted_tbl)
+}
+
+#' Plot GAM smooths.
+#'
+#' It plots the smooths from the estimates of a \link[mgcv]{gam} or \link[mgcv]{bam} object.
+#'
+#' @inheritParams get_gam_predictions
+#' @param comparison An unquoted expression indicating the model term for which the comparison will be plotted.
+#' @param facet_terms An unquoted formula with the terms used for faceting.
+#' @param conditions A list of quosures with \link[rlang]{quos} specifying the levels to plot from the model terms not among \code{time_series}, \code{comparison}, or \code{facet_terms}.
+#'
+#' @importFrom magrittr "%>%"
+#' @importFrom rlang ":="
+#' @importFrom stats "predict"
+#' @export
+plot_smooths <- function(model, time_series, comparison, facet_terms = NULL, conditions = NULL, exclude_random = TRUE, series_length = 25) {
+    time_series_q <- dplyr::enquo(time_series)
+    comparison_q <- dplyr::enquo(comparison)
+    facet_terms_q <- dplyr::enquo(facet_terms)
+    if (facet_terms_q == dplyr::quo(NULL)) {
+        facet_terms_q <- NULL
+    }
+    outcome_q <- model$formula[[2]]
+
+    predicted_tbl <- get_gam_predictions(model, !!time_series_q, conditions, exclude_random = exclude_random, series_length = series_length)
+
+    smooths_plot <- predicted_tbl %>%
+        ggplot2::ggplot(
+            ggplot2::aes_string(
+                dplyr::quo_name(time_series_q), dplyr::quo_name(outcome_q)
+            )
+        ) +
+        ggplot2::geom_ribbon(
+            ggplot2::aes_string(
+                ymin = "CI_lower",
+                ymax = "CI_upper",
+                fill = dplyr::quo_name(comparison_q)
+            ),
+            alpha = 0.2
+        ) +
+        ggplot2::geom_path(
+            ggplot2::aes_string(
+                colour = dplyr::quo_name(comparison_q),
+                linetype = dplyr::quo_name(comparison_q)
+            )
+        ) +
+        {if (!is.null(facet_terms_q)) {
+            ggplot2::facet_wrap(facet_terms_q)
+        }}
+
+    return(smooths_plot)
+}
+
+
+#' Plot GAM estimate smooths and difference curve.
 #'
 #' It plots comparison smooths from the estimates of a \link[mgcv]{gam} or \link[mgcv]{bam}
 #' and the difference curve. Significant differences are marked with red areas.
@@ -35,6 +206,8 @@ create_event_start <- function(tibble, event_col) {
 #' @importFrom magrittr "%>%"
 #' @export
 plot_gamsd <- function(model, view, comparison, conditions = NULL, rm_re = FALSE, bw = FALSE) {
+    .Deprecated("plot_smooth", msg = "'plot_gamsd' will be deprecated, use 'plot_smooths'. (Plotting with the difference smooth is not supported yet.)\n")
+
     diff.df <- itsadug::plot_diff(
         model,
         view = view,

R/zzz.R

@@ -0,0 +1 @@
+utils::globalVariables(c("fit", "se.fit"))

README.md

@@ -4,4 +4,6 @@ This is the repository of the `R` package `tidymv`. This package provides functi
 
 ## Installation
 
-To install the package, use `devtools::install_github("stefanocoretta/tidymv@v1.1.0", build_vignettes = TRUE)`. To learn how to use the package, do `vignette("plot-gamms", package = "tidymv")` after the installation.
+To install the package, use `devtools::install_github("stefanocoretta/tidymv@v1.3.0", build_vignettes = TRUE)`. To learn how to use the package, do `vignette("plot-smooths", package = "tidymv")` after the installation.
+
+If you wish to install the development version, use `devtools::install_github("stefanocoretta/tidymv", build_vignettes = TRUE)`.

vignettes/.gitignore

@@ -0,0 +1,2 @@
+*.html
+*.R

vignettes/plot-gamms.Rmd

@@ -17,7 +17,10 @@ library(itsadug)
 library(tidymv)
 ```
 
-Reference and difference smooths of factor predictos from a `gam` model can be plotted with `plot_gamsd()`.
+The estimated smooths and the difference smooth of factor predictors from a `gam` model can be plotted with `plot_gamsd()`.
+This function has now been deprecated and will be removed in future versions.
+This vignette is kept for historical reasons.
+
 To illustrate how to use `plot_gamsd()`, let's first prepare some dummy data with a factor variable and run `gam()` on this data. The `gam` model includes a reference smooth `s(x2)`, a by-factor difference smooth `s(x2, by = fac)`, and a smooth `s(x0)`.
 
 ```{r gam}