Skip to content

Commit

Permalink
Merge pull request #301 from fweber144/rstanarm_issue542_fixed
Browse files Browse the repository at this point in the history 
rstanarm issue 542 fixed; handling of offsets
  • Loading branch information
@fweber144
fweber144 authored Apr 20, 2022
2 parents c0d21f1 + 387cd30 commit dd4023a
Show file tree
Hide file tree
Showing 5 changed files with 78 additions and 120 deletions.
6 changes: 5 additions & 1 deletion NEWS.md
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,9 @@
# projpred 2.1.1.9000

## Minor changes

* Account for changes concerning the handling of offsets in **rstanarm** version 2.21.3. In particular, issue stan-dev/rstanarm#542 was fixed in **rstanarm** 2.21.3.

## Bug fixes

* Throw a more informative error message in case of special group-level terms which are currently not supported (in particular, nested ones).
Expand All @@ -10,7 +14,7 @@

* Fix the order of the package authors.
* Fix failing CRAN checks.
* Add an input check for argument `solution_terms` of `project()` to fix a test failure on R versions > 4.1.
* Add an input check for argument `solution_terms` of `project()` to fix a test failure in R versions >= 4.2.

# projpred 2.1.0

Expand Down
62 changes: 26 additions & 36 deletions R/refmodel.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -455,23 +455,14 @@ get_refmodel.stanreg <- function(object, ...) {
}

# Offsets:
# Element `stan_function` (needed for handling the offsets) is not documented
# in `?rstanarm::`stanreg-objects``, so check at least its length and type:
if (length(object$stan_function) != 1 ||
!is.vector(object$stan_function, mode = "character")) {
stop("Unexpected value of `object$stan_function`. Please notify the ",
"package maintainer.")
}
if (length(object$offset) > 0) {
if (is.null(attr(terms(formula(object), data = data), "offset"))) {
# In this case, we would have to use argument `offset` of
# posterior_linpred.stanreg() to allow for new offsets, requiring changes
# in all ref_predfun() calls. Furthermore, there is rstanarm issue #541.
# Thus, throw an error:
stop("It looks like `object` was fitted with offsets specified via ",
"argument `offset`. Currently, projpred does not support offsets ",
"specified this way. Please use an `offset()` term in the model ",
"formula instead.")
# Element `stan_function` (needed here for handling rstanarm issue #546) is
# not documented in `?rstanarm::`stanreg-objects``, so check at least its
# length and type:
if (length(object$stan_function) != 1 ||
!is.vector(object$stan_function, mode = "character")) {
stop("Unexpected value of `object$stan_function`. Please notify the ",
"package maintainer.")
}
if (object$stan_function == "stan_gamm4") {
stop("Because of rstanarm issue #546 (see GitHub), projpred cannot ",
Expand Down Expand Up @@ -531,27 +522,26 @@ get_refmodel.stanreg <- function(object, ...) {
}

ref_predfun <- function(fit, newdata = NULL) {
linpred_out <- t(
posterior_linpred(fit, newdata = newdata)
)
# Use a workaround for rstanarm issue #541 and rstanarm issue #542. This
# workaround consists of using `cond_no_offs` which indicates whether
# posterior_linpred() excluded (`TRUE`) or included (`FALSE`) the offsets:
cond_no_offs <- (
fit$stan_function %in% c("stan_lmer", "stan_glmer") &&
!is.null(attr(terms(formula), "offset"))
) || (
fit$stan_function %in% c("stan_lm", "stan_glm") &&
!is.null(newdata) && length(fit$offset) > 0
)
if (cond_no_offs) {
# Observation weights are not needed here, so use `wrhs = NULL` to avoid
# potential conflicts for a non-`NULL` default `wrhs`:
offs <- extract_model_data(fit, newdata = newdata, wrhs = NULL)$offset
stopifnot(length(offs) %in% c(1L, nrow(linpred_out)))
linpred_out <- linpred_out + offs
# The easiest way to deal with rstanarm issue #541 and rstanarm issue #542,
# changes between rstanarm versions 2.21.2 and 2.21.3 with respect to these
# issues, and the fact that offsets may be specified via argument `offset`
# of the respective model-fitting function (e.g., rstanarm::stan_glm()) is
# to include offsets explicitly in the call to
# rstanarm:::posterior_linpred.stanreg().

# Observation weights are not needed here, so use `wrhs = NULL` to avoid
# potential conflicts for a non-`NULL` default `wrhs`:
offs <- extract_model_data(fit, newdata = newdata, wrhs = NULL)$offset
n_obs <- nrow(newdata %||% data)
if (length(offs) == 0) {
offs <- rep(0, n_obs)
} else if (length(offs) == 1) {
offs <- rep(offs, n_obs)
} else if (length(offs) != n_obs) {
stop("Unexpected length of element `offset` returned by ",
"extract_model_data() (see `?init_refmodel`).")
}
return(linpred_out)
return(t(posterior_linpred(fit, newdata = newdata, offset = offs)))
}

cvfun <- function(folds) {
Expand Down
7 changes: 3 additions & 4 deletions tests/testthat/setup.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -360,11 +360,10 @@ seed_fit <- 74345
### Formula ---------------------------------------------------------------

# Notes:
# * Argument `offset` has an issue for rstanarm::stan_glmer() (see rstanarm
# issue #541). Instead, use offset() in the formula.
# * Argument `offset` is not supported by rstanarm::stan_gamm4(). Instead, use
# offset() in the formula. However, because of rstanarm issue #546 and
# rstanarm issue #253, omit the offsets in GAMs and GAMMs.
# offset() in the formula (like for all other models). However, because of
# rstanarm issue #546 and rstanarm issue #253, omit the offsets in GAMs and
# GAMMs.
# * In rstanarm::stan_gamm4(), multilevel terms are specified via argument
# `random`.

Expand Down
80 changes: 36 additions & 44 deletions tests/testthat/test_refmodel.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,18 +62,37 @@ test_that("`formula` as a character string fails", {
"^inherits\\(formula, \"formula\"\\) is not TRUE$")
})

test_that("offsets specified via argument `offset` fail", {
fit_offs_arg <- suppressWarnings(rstanarm::stan_glm(
y_glm_gauss ~ xco.1 + xco.2 + xco.3 + xca.1 + xca.2,
family = f_gauss, data = dat,
weights = wobs_tst, offset = offs_tst,
chains = chains_tst, seed = seed_fit, iter = iter_tst, QR = TRUE,
refresh = 0
test_that("offsets specified via argument `offset` work", {
args_fit_i <- args_fit$rstanarm.glm.gauss.stdformul.with_wobs.with_offs
skip_if_not(!is.null(args_fit_i))
fit_fun_nm <- switch(args_fit_i$pkg_nm,
"rstanarm" = switch(args_fit_i$mod_nm,
"glm" = "stan_glm",
"glmm" = "stan_glmer",
"stan_gamm4"),
"brms" = "brm",
stop("Unknown `pkg_nm`."))
upd_no_offs <- paste(". ~", sub(" \\+ offset\\(offs_col\\)", "",
as.character(args_fit_i$formula[3])))
fit_offs_arg <- suppressWarnings(do.call(
get(fit_fun_nm, asNamespace(args_fit_i$pkg_nm)),
c(list(formula = update(args_fit_i$formula, upd_no_offs),
offset = offs_tst),
excl_nonargs(args_fit_i, nms_excl_add = "formula"))
))
expect_error(
get_refmodel(fit_offs_arg),
paste("^It looks like `object` was fitted with offsets specified via",
"argument `offset`\\.")
refmod_offs_arg <- get_refmodel(fit_offs_arg)
expect_equal(
as.matrix(fit_offs_arg),
as.matrix(fits$rstanarm.glm.gauss.stdformul.with_wobs.with_offs),
tolerance = 1e-12,
info = "rstanarm.glm.gauss.stdformul.with_wobs.with_offs"
)
nms_compare <- c("mu", "dis", "y", "loglik", "wobs", "wsample", "offset")
expect_equal(
refmod_offs_arg[nms_compare],
refmods$rstanarm.glm.gauss.stdformul.with_wobs.with_offs[nms_compare],
tolerance = 1e-12,
info = "rstanarm.glm.gauss.stdformul.with_wobs.with_offs"
)
})

Expand Down Expand Up @@ -131,44 +150,17 @@ test_that(paste(
mod_crr <- args_ref[[tstsetup]]$mod_nm
fam_crr <- args_ref[[tstsetup]]$fam_nm

# We expect a warning which in fact should be suppressed, though (see
# issue #162):
warn_expected <- if (pkg_crr == "rstanarm" &&
mod_crr == "glm" &&
grepl("\\.with_offs", tstsetup)) {
paste("^'offset' argument is NULL but it looks like you",
"estimated the model using an offset term\\.$")
} else {
NA
}

y_crr <- dat[, paste("y", mod_crr, fam_crr, sep = "_")]

# Without `ynew`:
expect_warning(
predref_resp <- predict(refmods[[tstsetup]], dat, type = "response"),
warn_expected,
info = tstsetup
)
expect_warning(
predref_link <- predict(refmods[[tstsetup]], dat, type = "link"),
warn_expected,
info = tstsetup
)
predref_resp <- predict(refmods[[tstsetup]], dat, type = "response")
predref_link <- predict(refmods[[tstsetup]], dat, type = "link")

# With `ynew`:
expect_warning(
predref_ynew_resp <- predict(refmods[[tstsetup]], dat, ynew = y_crr,
type = "response"),
warn_expected,
info = tstsetup
)
expect_warning(
predref_ynew_link <- predict(refmods[[tstsetup]], dat, ynew = y_crr,
type = "link"),
warn_expected,
info = tstsetup
)
predref_ynew_resp <- predict(refmods[[tstsetup]], dat, ynew = y_crr,
type = "response")
predref_ynew_link <- predict(refmods[[tstsetup]], dat, ynew = y_crr,
type = "link")

# Checks without `ynew`:
expect_true(is.vector(predref_resp, "double"), info = tstsetup)
Expand Down
43 changes: 8 additions & 35 deletions tests/testthat/test_varsel.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,24 +91,10 @@ test_that("`d_test` works", {
type = "test",
offset = refmod_crr$offset
)
# We expect a warning which in fact should be suppressed, though (see
# issue #162):
warn_expected <- if (pkg_crr == "rstanarm" &&
mod_crr == "glm" &&
grepl("\\.with_offs", tstsetup)) {
paste("^'offset' argument is NULL but it looks like you estimated the",
"model using an offset term\\.$")
} else {
NA
}
expect_warning(
vs_repr <- do.call(varsel, c(
list(object = refmod_crr, d_test = d_test_crr),
excl_nonargs(args_vs_i)
)),
warn_expected,
info = tstsetup
)
vs_repr <- do.call(varsel, c(
list(object = refmod_crr, d_test = d_test_crr),
excl_nonargs(args_vs_i)
))
meth_exp_crr <- args_vs_i$method
if (is.null(meth_exp_crr)) {
meth_exp_crr <- ifelse(mod_crr == "glm", "L1", "forward")
Expand Down Expand Up @@ -756,23 +742,10 @@ test_that(paste(
refmod_crr <- get_refmodel(fit_crr, cvfits = kfold_obj)

# Run cv_varsel():
# We expect a warning which in fact should be suppressed, though (see
# issue #162):
warn_expected <- if (mod_crr == "glm" &&
grepl("\\.with_offs", tstsetup)) {
paste("^'offset' argument is NULL but it looks like you estimated the",
"model using an offset term\\.$")
} else {
NA
}
expect_warning(
cvvs_cvfits <- do.call(cv_varsel, c(
list(object = refmod_crr),
excl_nonargs(args_cvvs_i, nms_excl_add = "K")
)),
warn_expected,
info = tstsetup
)
cvvs_cvfits <- do.call(cv_varsel, c(
list(object = refmod_crr),
excl_nonargs(args_cvvs_i, nms_excl_add = "K")
))

# Checks:
vsel_tester(
Expand Down

0 comments on commit dd4023a

Please sign in to comment.