Skip to content

Commit

Permalink
Suppress compile warnings in tests, skip examples for Rwin
Browse files Browse the repository at this point in the history
  • Loading branch information
andrjohns committed Jun 6, 2024
1 parent 9c3f646 commit db0d7bc
Show file tree
Hide file tree
Showing 70 changed files with 2,950 additions and 2,950 deletions.
4 changes: 4 additions & 0 deletions .github/workflows/R-CMD-check.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -53,4 +53,8 @@ jobs:
- uses: r-lib/actions/setup-r-dependencies@v2-branch
with:
extra-packages: any::rcmdcheck any::betareg any::HSAUR3 any::biglm any::gamm4 any::V8
- name: Add flags to Makevars
run: |
mkdir -p ~/.R
echo "CXXFLAGS= -w -ftemplate-backtrace-limit=0" >> ~/.R/Makevars
- uses: r-lib/actions/check-r-package@v2-branch
65 changes: 32 additions & 33 deletions R/as.matrix.stanreg.R
Original file line number Diff line number Diff line change
@@ -1,71 +1,71 @@
# Part of the rstanarm package for estimating model parameters
# Copyright (C) 2015, 2016, 2017 Trustees of Columbia University
#
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 3
# of the License, or (at your option) any later version.
#
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

#' Extract the posterior sample
#'
#'
#' For models fit using MCMC (\code{algorithm="sampling"}), the posterior sample
#' ---the post-warmup draws from the posterior distribution--- can be extracted
#' from a fitted model object as a matrix, data frame, or array. The
#' \code{as.matrix} and \code{as.data.frame} methods merge all chains together,
#' whereas the \code{as.array} method keeps the chains separate. For models fit
#' using optimization (\code{"optimizing"}) or variational inference
#' (\code{"meanfield"} or \code{"fullrank"}), there is no posterior sample but
#' ---the post-warmup draws from the posterior distribution--- can be extracted
#' from a fitted model object as a matrix, data frame, or array. The
#' \code{as.matrix} and \code{as.data.frame} methods merge all chains together,
#' whereas the \code{as.array} method keeps the chains separate. For models fit
#' using optimization (\code{"optimizing"}) or variational inference
#' (\code{"meanfield"} or \code{"fullrank"}), there is no posterior sample but
#' rather a matrix (or data frame) of 1000 draws from either the asymptotic
#' multivariate Gaussian sampling distribution of the parameters or the
#' variational approximation to the posterior distribution.
#'
#'
#' @method as.matrix stanreg
#' @export
#' @templateVar stanregArg x
#' @template args-stanreg-object
#' @template args-pars
#' @template args-regex-pars
#' @param ... Ignored.
#'
#'
#' @return A matrix, data.frame, or array, the dimensions of which depend on
#' \code{pars} and \code{regex_pars}, as well as the model and estimation
#' algorithm (see the Description section above).
#'
#'
#' @seealso \code{\link{stanreg-draws-formats}}, \code{\link{stanreg-methods}}
#'
#'
#' @examples
#' if (.Platform$OS.type != "windows" || .Platform$r_arch != "i386") {
#' if (.Platform$OS.type != "windows") {
#' \donttest{
#' if (!exists("example_model")) example(example_model)
#' # Extract posterior sample after MCMC
#' draws <- as.matrix(example_model)
#' print(dim(draws))
#'
#' # For example, we can see that the median of the draws for the intercept
#'
#' # For example, we can see that the median of the draws for the intercept
#' # is the same as the point estimate rstanarm uses
#' print(median(draws[, "(Intercept)"]))
#' print(example_model$coefficients[["(Intercept)"]])
#'
#'
#' # The as.array method keeps the chains separate
#' draws_array <- as.array(example_model)
#' print(dim(draws_array)) # iterations x chains x parameters
#'
#' # Extract draws from asymptotic Gaussian sampling distribution
#'
#' # Extract draws from asymptotic Gaussian sampling distribution
#' # after optimization
#' fit <- stan_glm(mpg ~ wt, data = mtcars, algorithm = "optimizing")
#' draws <- as.data.frame(fit)
#' print(colnames(draws))
#' print(nrow(draws)) # 1000 draws are taken
#'
#'
#' # Extract draws from variational approximation to the posterior distribution
#' fit2 <- update(fit, algorithm = "meanfield")
#' draws <- as.data.frame(fit2, pars = "wt")
Expand All @@ -76,10 +76,10 @@
as.matrix.stanreg <- function(x, ..., pars = NULL, regex_pars = NULL) {
pars <- collect_pars(x, pars, regex_pars)
user_pars <- !is.null(pars)

if (used.optimizing(x)) {
mat <- x$asymptotic_sampling_dist
if (is.null(mat))
if (is.null(mat))
STOP_no_draws()
if (!user_pars) {
aux <- c("sigma", "scale", "shape", "lambda", "reciprocal_dispersion")
Expand Down Expand Up @@ -107,21 +107,21 @@ as.array.stanreg <- function(x, ..., pars = NULL, regex_pars = NULL) {
pars <- collect_pars(x, pars, regex_pars)
if (!used.sampling(x))
stop(
"For models not fit using MCMC ",
"For models not fit using MCMC ",
"use 'as.matrix' instead of 'as.array'"
)

arr <- as.array(x$stanfit)
if (identical(arr, numeric(0)))
STOP_no_draws()

if (!is.null(pars)) {
check_missing_pars(arr, pars)
} else {
pars <- exclude_lp_and_ppd(last_dimnames(arr))
}
arr <- arr[, , pars, drop = FALSE]

if (!is.mer(x))
return(arr)
unpad_reTrms(arr)
Expand All @@ -143,20 +143,19 @@ STOP_no_draws <- function() stop("No draws found.", call. = FALSE)

check_missing_pars <- function(x, pars) {
notfound <- which(!pars %in% last_dimnames(x))
if (length(notfound))
if (length(notfound))
stop(
"No parameter(s) ",
paste(pars[notfound], collapse = ", "),
"No parameter(s) ",
paste(pars[notfound], collapse = ", "),
call. = FALSE
)
}

exclude_lp_and_ppd <- function(pars) {
grep(
pattern = "mean_PPD|log-posterior",
x = pars,
invert = TRUE,
pattern = "mean_PPD|log-posterior",
x = pars,
invert = TRUE,
value = TRUE
)
}

50 changes: 25 additions & 25 deletions R/bayes_R2.R
Original file line number Diff line number Diff line change
Expand Up @@ -8,57 +8,57 @@
#' @param re.form For models with group-level terms, \code{re.form} is
#' passed to \code{\link{posterior_epred}} if specified.
#' @param ... Currently ignored.
#'
#'
#' @return A vector of R-squared values with length equal to the posterior
#' sample size (the posterior distribution of R-squared).
#'
#'
#' @references
#' Andrew Gelman, Ben Goodrich, Jonah Gabry, and Aki Vehtari (2018). R-squared
#' for Bayesian regression models. \emph{The American Statistician}, to appear.
#' \doi{10.1080/00031305.2018.1549100}
#' (\href{http://www.stat.columbia.edu/~gelman/research/published/bayes_R2_v3.pdf}{Preprint},
#' \href{https://avehtari.github.io/bayes_R2/bayes_R2.html}{Notebook})
#'
#'
#' @examples
#' if (.Platform$OS.type != "windows" || .Platform$r_arch != "i386") {
#' if (.Platform$OS.type != "windows") {
#' fit <- stan_glm(
#' mpg ~ wt + cyl,
#' data = mtcars,
#' QR = TRUE,
#' chains = 2,
#' mpg ~ wt + cyl,
#' data = mtcars,
#' QR = TRUE,
#' chains = 2,
#' refresh = 0
#' )
#' rsq <- bayes_R2(fit)
#' print(median(rsq))
#' hist(rsq)
#'
#'
#' loo_rsq <- loo_R2(fit)
#' print(median(loo_rsq))
#'
#'
#' # multilevel binomial model
#' if (!exists("example_model")) example(example_model)
#' print(example_model)
#' median(bayes_R2(example_model))
#' median(bayes_R2(example_model, re.form = NA)) # exclude group-level
#' }
bayes_R2.stanreg <- function(object, ..., re.form = NULL) {

if (!used.sampling(object))
STOP_sampling_only("bayes_R2")
if (is_polr(object))
stop("bayes_R2 is not available for stan_polr models.")

fam <- family(object)$family
if (!fam %in% c("gaussian", "binomial")) {
stop("bayes_R2 is only available for Gaussian and binomial models.")
}

mu_pred <- posterior_epred(object, re.form = re.form)
if (is.binomial(fam)) {
y <- get_y(object)
if (NCOL(y) == 2) {
trials <- rowSums(y)
trials_mat <- matrix(trials, nrow = nrow(mu_pred), ncol = ncol(mu_pred),
trials_mat <- matrix(trials, nrow = nrow(mu_pred), ncol = ncol(mu_pred),
byrow = TRUE)
tmp <- mu_pred * trials_mat
sigma2 <- rowMeans(tmp * (1 - mu_pred))
Expand All @@ -69,7 +69,7 @@ bayes_R2.stanreg <- function(object, ..., re.form = NULL) {
} else {
sigma2 <- drop(as.matrix(object, pars = "sigma"))^2
}

var_mu_pred <- apply(mu_pred, 1, var)
r_squared <- var_mu_pred / (var_mu_pred + sigma2)
return(r_squared)
Expand All @@ -80,40 +80,40 @@ bayes_R2.stanreg <- function(object, ..., re.form = NULL) {
#' @aliases loo_R2
#' @importFrom rstantools loo_R2
#' @export
#'
#'
loo_R2.stanreg <- function(object, ...) {
if (!used.sampling(object))
STOP_sampling_only("loo_R2")
if (is_polr(object))
stop("loo_R2 is not available for stan_polr models.")

fam <- family(object)$family
if (!fam %in% c("gaussian", "binomial")) {
stop("loo_R2 is only available for Gaussian and binomial models.")
}

y <- get_y(object)
log_ratios <- -log_lik(object)
psis_object <- object[["loo"]][["psis_object"]]
if (is.null(psis_object)) {
psis_object <- loo::psis(log_ratios, r_eff = NA)
}

mu_pred <- posterior_epred(object)
if (is.binomial(fam)) {
if (is.factor(y)) {
y <- fac2bin(y)
} else if (NCOL(y) == 2) {
trials <- rowSums(y)
y <- y[, 1]
trials_mat <- matrix(trials, nrow = nrow(mu_pred), ncol = ncol(mu_pred),
trials_mat <- matrix(trials, nrow = nrow(mu_pred), ncol = ncol(mu_pred),
byrow = TRUE)
mu_pred <- mu_pred * trials_mat
}
}
mu_pred_loo <- loo::E_loo(mu_pred, psis_object, log_ratios = log_ratios)$value
err_loo <- mu_pred_loo - y

S <- nrow(mu_pred)
N <- ncol(mu_pred)

Expand All @@ -123,16 +123,16 @@ loo_R2.stanreg <- function(object, ...) {
on.exit(assign(".Random.seed", rng_state_old, envir = .GlobalEnv))
set.seed(object$stanfit@stan_args[[1]]$seed)

# dirichlet weights
# dirichlet weights
exp_draws <- matrix(rexp(S * N, rate = 1), nrow = S, ncol = N)
wts <- exp_draws / rowSums(exp_draws)

var_y <- (rowSums(sweep(wts, 2, y^2, FUN = "*")) -
rowSums(sweep(wts, 2, y, FUN = "*"))^2) * (N/(N-1))

var_err_loo <- (rowSums(sweep(wts, 2, err_loo^2, FUN = "*")) -
rowSums(sweep(wts, 2, err_loo, FUN = "*")^2)) * (N/(N-1))

loo_r_squared <- 1 - var_err_loo / var_y
loo_r_squared[loo_r_squared < -1] <- -1
loo_r_squared[loo_r_squared > 1] <- 1
Expand Down
Loading

0 comments on commit db0d7bc

Please sign in to comment.