Skip to content

Commit

Permalink
Fix issue #156 (#269)
Browse files Browse the repository at this point in the history
* Use `t2()` instead of `s()` in the demonstration of

a smooth term with multiple predictors. The reason is that `nitro` and `bv` measure different things (unlike spatial coordinates, for example) and `s()` uses a thin plate regression spline by default (which is *isotropic*).

* Fix issue #156.
  • Loading branch information
fweber144 authored Jan 24, 2022
1 parent a9ba38b commit ce3b2fd
Show file tree
Hide file tree
Showing 4 changed files with 33 additions and 17 deletions.
10 changes: 2 additions & 8 deletions R/formula.R
Original file line number Diff line number Diff line change
Expand Up @@ -84,14 +84,8 @@ extract_response <- function(response) {
return(response_name_ch)
}

## Parse additive terms from a list of individual terms. These may include
## individual smoothers s(x), interaction smoothers s(x, z), smooth intercepts
## s(factor, bs = "re"), group level independent smoothers s(x, by = g), or
## shared smoothers for low dimensional factors s(x, g, bs = "fs"), interaction
## group level smoothers with same shape for high dimensional factors t2(x, z,
## g, bs = c(., ., "re")), interaction group level independent smoothers te(x,
## z, by = g), population level interaction smoothers te(x, z, bs = c(., .)),
## s(x, z).
## Parse additive terms (smooth terms) from a list of individual terms. See
## `?init_refmodel` for allowed smooth terms.
## @param terms list of terms to parse
## @return a vector of smooth terms
parse_additive_terms <- function(terms) {
Expand Down
17 changes: 14 additions & 3 deletions R/refmodel.R
Original file line number Diff line number Diff line change
Expand Up @@ -21,9 +21,10 @@
#' @param data Data used for fitting the reference model.
#' @param formula Reference model's formula. For general information on formulas
#' in \R, see [`formula`]. For multilevel formulas, see also package
#' \pkg{lme4}, in particular [lme4::lmer()] and [lme4::glmer()]. For additive
#' formulas, see also packages \pkg{mgcv}, in particular [mgcv::gam()], and
#' \pkg{gamm4}, in particular [gamm4::gamm4()].
#' \pkg{lme4} (in particular, functions [lme4::lmer()] and [lme4::glmer()]).
#' For additive formulas, see also packages \pkg{mgcv} (in particular,
#' function [mgcv::gam()]) and \pkg{gamm4} (in particular, function
#' [gamm4::gamm4()]) as well as the notes in section "Formula terms" below.
#' @param ref_predfun Prediction function for the linear predictor of the
#' reference model, including offsets (if existing). See also section
#' "Arguments `ref_predfun`, `proj_predfun`, and `div_minimizer`" below. If
Expand Down Expand Up @@ -67,6 +68,16 @@
#'
#' @details
#'
#' # Formula terms
#'
#' For additive models (still an experimental feature), only [mgcv::s()] and
#' [mgcv::t2()] are currently supported as smooth terms. Furthermore, these need
#' to be called without any arguments apart from the predictor names (symbols).
#' For example, for smoothing the effect of a predictor `x`, only `s(x)` or
#' `t2(x)` are allowed. As another example, for smoothing the joint effect of
#' two predictors `x` and `z`, only `s(x, z)` or `t2(x, z)` are allowed (and
#' analogously for higher-order joint effects, e.g., of three predictors).
#'
#' # Arguments `ref_predfun`, `proj_predfun`, and `div_minimizer`
#'
#' Arguments `ref_predfun`, `proj_predfun`, and `div_minimizer` may be `NULL`
Expand Down
17 changes: 14 additions & 3 deletions man/refmodel-init-get.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

6 changes: 3 additions & 3 deletions vignettes/projpred.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -252,7 +252,7 @@ This PPPC shows that our final projection is able to generate predictions simila

## Supported types of reference models {#refmodtypes}

Apart from the `gaussian()` response family used in this vignette, **projpred** also supports the `binomial()`^[The `binomial()` family includes the `brms::bernoulli()` family as a special case which may only be used via the **brms** package.] and the `poisson()` family. On the side of the predictors, **projpred** not only supports linear main effects as shown in this vignette, but also interactions, multilevel^[Multilevel models are also known as *hierarchical* models or models with *partially pooled*, *group-level*, or---in frequentist terms---*random* effects.], and---as an experimental feature---additive^[Additive terms are also known as *smoothing* terms.] terms.
Apart from the `gaussian()` response family used in this vignette, **projpred** also supports the `binomial()`^[The `binomial()` family includes the `brms::bernoulli()` family as a special case which may only be used via the **brms** package.] and the `poisson()` family. On the side of the predictors, **projpred** not only supports linear main effects as shown in this vignette, but also interactions, multilevel^[Multilevel models are also known as *hierarchical* models or models with *partially pooled*, *group-level*, or---in frequentist terms---*random* effects.], and---as an experimental feature---additive^[Additive terms are also known as *smooth* terms.] terms.

Transferring this vignette to such more complex reference models is straightforward: Basically, only the code for fitting the reference model via **rstanarm** or **brms** needs to be adapted. The **projpred** code stays almost the same. Only note that in case of multilevel or additive reference models,
<!-- the parallelization of the projection is not recommended and that -->
Expand All @@ -274,10 +274,10 @@ As an example for an additive (non-multilevel) reference model, consider the `la
data("lasrosas.corn", package = "agridat")
# Convert `year` to a `factor` (this could also be solved by using
# `factor(year)` in the formula, but we avoid that here to put more emphasis on
# the demonstration of the smoothing term):
# the demonstration of the smooth term):
lasrosas.corn$year <- as.factor(lasrosas.corn$year)
refm_fit <- stan_gamm4(
yield ~ year + topo + s(nitro, bv),
yield ~ year + topo + t2(nitro, bv),
family = gaussian(),
data = lasrosas.corn,
seed = 4919670, QR = TRUE, refresh = 0
Expand Down

0 comments on commit ce3b2fd

Please sign in to comment.