diff --git a/DESCRIPTION b/DESCRIPTION index 96c51232b..0e224ff5f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -7,23 +7,23 @@ Encoding: UTF-8 Authors@R: c(person("Jonah", "Gabry", email = "jsg2201@columbia.edu", role = "aut"), person("Imad", "Ali", role = "ctb"), person("Sam", "Brilleman", role = "ctb"), - person(given = "Jacqueline Buros", family = "Novik", + person(given = "Jacqueline Buros", family = "Novik", role = "ctb", comment = "R/stan_jm.R"), person("AstraZeneca", role = "ctb", comment = "R/stan_jm.R"), person("Trustees of", "Columbia University", role = "cph"), person("Simon", "Wood", role = "cph", comment = "R/stan_gamm4.R"), - person("R Core", "Deveopment Team", role = "cph", + person("R Core", "Deveopment Team", role = "cph", comment = "R/stan_aov.R"), person("Douglas", "Bates", role = "cph", comment = "R/pp_data.R"), person("Martin", "Maechler", role = "cph", comment = "R/pp_data.R"), person("Ben", "Bolker", role = "cph", comment = "R/pp_data.R"), person("Steve", "Walker", role = "cph", comment = "R/pp_data.R"), - person("Brian", "Ripley", role = "cph", + person("Brian", "Ripley", role = "cph", comment = "R/stan_aov.R, R/stan_polr.R"), person("William", "Venables", role = "cph", comment = "R/stan_polr.R"), person("Paul-Christian", "Burkner", email = "paul.buerkner@gmail.com", role = "cph", comment = "R/misc.R"), - person("Ben", "Goodrich", email = "benjamin.goodrich@columbia.edu", + person("Ben", "Goodrich", email = "benjamin.goodrich@columbia.edu", role = c("cre", "aut"))) Description: Estimates previously compiled regression models using the 'rstan' package, which provides the R interface to the Stan C++ library for Bayesian @@ -54,6 +54,7 @@ Suggests: digest, gridExtra, HSAUR3, + INLA, knitr (>= 1.15.1), MASS, mgcv (>= 1.8-13), diff --git a/NAMESPACE b/NAMESPACE index 5e865a6b9..a0509c702 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -89,10 +89,12 @@ S3method(update,stanmvreg) S3method(update,stanreg) S3method(vcov,stanreg) S3method(waic,stanreg) +export() export(R2) export(Surv) export(VarCorr) export(bayes_R2) +export(beta) export(cauchy) export(compare_models) export(decov) @@ -148,6 +150,7 @@ export(ranef) export(se) export(sigma) export(stan_aov) +export(stan_besag) export(stan_betareg) export(stan_betareg.fit) export(stan_biglm) @@ -168,6 +171,7 @@ export(stan_mvmer) export(stan_nlmer) export(stan_polr) export(stan_polr.fit) +export(stan_spatial.fit) export(stanjm_list) export(stanmvreg_list) export(stanreg_list) diff --git a/R/doc-datasets.R b/R/doc-datasets.R index 8daa1d214..103ac5ab9 100644 --- a/R/doc-datasets.R +++ b/R/doc-datasets.R @@ -1,34 +1,34 @@ # 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. # #' Datasets for rstanarm examples -#' +#' #' Small datasets for use in \pkg{rstanarm} examples and vignettes. #' #' @name rstanarm-datasets #' @aliases kidiq roaches wells bball1970 bball2006 mortality tumors radon pbcLong pbcSurv -#' @format +#' @format #' \describe{ #' \item{\code{bball1970}}{ #' Data on hits and at-bats from the 1970 Major League Baseball season for 18 #' players. -#' +#' #' Source: Efron and Morris (1975). -#' +#' #' 18 obs. of 5 variables #' \itemize{ #' \item \code{Player} Player's last name @@ -41,9 +41,9 @@ #' \item{\code{bball2006}}{ #' Hits and at-bats for the entire 2006 American League season of Major League #' Baseball. -#' +#' #' Source: Carpenter (2009) -#' +#' #' 302 obs. of 2 variables #' \itemize{ #' \item \code{y} Number of hits @@ -51,11 +51,11 @@ #' } #' } #' \item{\code{kidiq}}{ -#' Data from a survey of adult American women and their children +#' Data from a survey of adult American women and their children #' (a subsample from the National Longitudinal Survey of Youth). -#' +#' #' Source: Gelman and Hill (2007) -#' +#' #' 434 obs. of 4 variables #' \itemize{ #' \item \code{kid_score} Child's IQ score @@ -67,7 +67,7 @@ #' \item{\code{mortality}}{ #' Surgical mortality rates in 12 hospitals performing cardiac surgery #' in babies. -#' +#' #' Source: Spiegelhalter et al. (1996). #' #' 12 obs. of 2 variables @@ -77,38 +77,38 @@ #' } #' } #' \item{\code{pbcLong,pbcSurv}}{ -#' Longitudinal biomarker and time-to-event survival data for 40 patients -#' with primary biliary cirrhosis who participated in a randomised +#' Longitudinal biomarker and time-to-event survival data for 40 patients +#' with primary biliary cirrhosis who participated in a randomised #' placebo controlled trial of D-penicillamine conducted at the Mayo #' Clinic between 1974 and 1984. -#' +#' #' Source: Therneau and Grambsch (2000) -#' +#' #' 304 obs. of 8 variables (\code{pbcLong}) and 40 obs. of 7 variables (\code{pbcSurv}) #' \itemize{ #' \item \code{age} {in years} #' \item \code{albumin} {serum albumin (g/dl)} #' \item \code{logBili} {logarithm of serum bilirubin} -#' \item \code{death} {indicator of death at endpoint} -#' \item \code{futimeYears} {time (in years) between baseline and +#' \item \code{death} {indicator of death at endpoint} +#' \item \code{futimeYears} {time (in years) between baseline and #' the earliest of death, transplantion or censoring} #' \item \code{id} {numeric ID unique to each individual} #' \item \code{platelet} {platelet count} #' \item \code{sex} {gender (m = male, f = female)} -#' \item \code{status} {status at endpoint (0 = censored, +#' \item \code{status} {status at endpoint (0 = censored, #' 1 = transplant, 2 = dead)} -#' \item \code{trt} {binary treatment code (0 = placebo, 1 = +#' \item \code{trt} {binary treatment code (0 = placebo, 1 = #' D-penicillamine)} #' \item \code{year} {time (in years) of the longitudinal measurements, #' taken as time since baseline)} #' } #' } -#' +#' #' \item{\code{radon}}{ -#' Data on radon levels in houses in the state of Minnesota. -#' +#' Data on radon levels in houses in the state of Minnesota. +#' #' Source: Gelman and Hill (2007) -#' +#' #' 919 obs. of 4 variables #' \itemize{ #' \item \code{log_radon} Radon measurement from the house (log scale) @@ -121,9 +121,9 @@ #' \item{\code{roaches}}{ #' Data on the efficacy of a pest management system at reducing the number of #' roaches in urban apartments. -#' +#' #' Source: Gelman and Hill (2007) -#' +#' #' 262 obs. of 6 variables #' \itemize{ #' \item \code{y} Number of roaches caught @@ -136,10 +136,10 @@ #' \item{\code{tumors}}{ #' Tarone (1982) provides a data set of tumor incidence in historical #' control groups of rats; specifically endometrial stromal polyps in -#' female lab rats of type F344. -#' +#' female lab rats of type F344. +#' #' Source: Gelman and Hill (2007) -#' +#' #' 71 obs. of 2 variables #' \itemize{ #' \item \code{y} Number of rats with tumors @@ -153,9 +153,9 @@ #' safe public or private well in the nearby area and the survey was conducted #' several years later to learn which of the affected residents had switched #' wells. -#' +#' #' Souce: Gelman and Hill (2007) -#' +#' #' 3020 obs. of 5 variables #' \itemize{ #' \item \code{switch} Indicator for well-switching @@ -167,27 +167,46 @@ #' \item \code{educ} Years of education (head of household) #' } #' } +#' \item{\code{lattice}}{ The \code{grid_sim15} and \code{grid_sim30} +#' SpatialPolygonsDataFrame objects are a simulated lattice of 225 and 900 +#' spatial units, respectively. These objects contain spatially dependent +#' observations simulated as a Gaussian Markov Random Field (GMRF). The data +#' slot in \code{grid_sim15} contains simulated outcomes using the GMRF and +#' various supported likelihoods. +#' +#' Source: Simulated data +#' +#' \itemize{ +#' \item \code{gmrf} GMRF data simulated using Cholesky decomposition sampling +#' from the precision form of the Multivariate Normal distribution centered at +#' zero. +#' \item \code{r} X-position of a location on the grid (origin is at the +#' bottom-left). +#' \item \code{c} Y-position of a location on the grid (origin is at the +#' bottom-left). #' } -#' -#' @references +#' } +#' } +#' +#' @references #' Carpenter, B. (2009) Bayesian estimators for the beta-binomial model of #' batting ability. \url{http://lingpipe-blog.com/2009/09/23/} -#' +#' #' Efron, B. and Morris, C. (1975) Data analysis using Stein's estimator and its #' generalizations. \emph{Journal of the American Statistical Association} #' \strong{70}(350), 311--319. -#' +#' #' @templateVar armRef \url{http://stat.columbia.edu/~gelman/arm/} #' @template reference-gelman-hill -#' +#' #' @references -#' Spiegelhalter, D., Thomas, A., Best, N., & Gilks, W. (1996) BUGS 0.5 +#' Spiegelhalter, D., Thomas, A., Best, N., & Gilks, W. (1996) BUGS 0.5 #' Examples. MRC Biostatistics Unit, Institute of Public health, Cambridge, UK. -#' +#' #' Tarone, R. E. (1982) The use of historical control information in testing for #' a trend in proportions. \emph{Biometrics} \strong{38}(1):215--220. -#' -#' Therneau, T. and Grambsch, P. (2000) \emph{Modeling Survival Data: Extending +#' +#' Therneau, T. and Grambsch, P. (2000) \emph{Modeling Survival Data: Extending #' the Cox Model}. Springer-Verlag, New York, US. #' #' @examples @@ -200,7 +219,7 @@ #' pp_check(fit, nreps = 20) #' \donttest{ #' bayesplot::color_scheme_set("brightblue") -#' pp_check(fit, plotfun = "stat_grouped", stat = "median", +#' pp_check(fit, plotfun = "stat_grouped", stat = "median", #' group = factor(kidiq$mom_hs, labels = c("No HS", "HS"))) #' } #' } diff --git a/R/log_lik.R b/R/log_lik.R index d83e6bea4..15523da2b 100644 --- a/R/log_lik.R +++ b/R/log_lik.R @@ -21,7 +21,7 @@ #' \eqn{S} by \eqn{N} pointwise log-likelihood matrix, where \eqn{S} is the size #' of the posterior sample and \eqn{N} is the number of data points, or in the #' case of the \code{stanmvreg} method (when called on \code{\link{stan_jm}} -#' model objects) an \eqn{S} by \eqn{Npat} matrix where \eqn{Npat} is the number +#' model objects) an \eqn{S} by \eqn{Npat} matrix where \eqn{Npat} is the number #' of individuals. #' #' @aliases log_lik @@ -36,12 +36,11 @@ #' @param offset A vector of offsets. Only required if \code{newdata} is #' specified and an \code{offset} was specified when fitting the model. #' -#' @return For the \code{stanreg} and \code{stanmvreg} methods an \eqn{S} by -#' \eqn{N} matrix, where \eqn{S} is the size of the posterior sample and -#' \eqn{N} is the number of data points. For the \code{stanjm} method +#' @return For the \code{stanreg} and \code{stanmvreg} methods an \eqn{S} by +#' \eqn{N} matrix, where \eqn{S} is the size of the posterior sample and +#' \eqn{N} is the number of data points. For the \code{stanjm} method #' an \eqn{S} by \eqn{Npat} matrix where \eqn{Npat} is the number of individuals. -#' -#' +#' #' @examples #' if (.Platform$OS.type != "windows" || .Platform$r_arch != "i386") { #' \donttest{ @@ -75,7 +74,7 @@ log_lik.stanreg <- function(object, newdata = NULL, offset = NULL, ...) { dots <- list(...) if (is.stanmvreg(object)) { m <- dots[["m"]] - if (is.null(m)) + if (is.null(m)) STOP_arg_required_for_stanmvreg(m) if (!is.null(offset)) stop2("'offset' cannot be specified for stanmvreg objects.") @@ -123,7 +122,7 @@ log_lik.stanreg <- function(object, newdata = NULL, offset = NULL, ...) { #' @export #' @templateVar mArg m #' @template args-m -#' +#' log_lik.stanmvreg <- function(object, m = 1, newdata = NULL, ...) { validate_stanmvreg_object(object) out <- log_lik.stanreg(object, newdata = newdata, m = m, ...) @@ -132,16 +131,16 @@ log_lik.stanmvreg <- function(object, m = 1, newdata = NULL, ...) { #' @rdname log_lik.stanreg #' @export -#' @param newdataLong,newdataEvent Optional data frames containing new data -#' (e.g. holdout data) to use when evaluating the log-likelihood for a -#' model estimated using \code{\link{stan_jm}}. If the fitted model +#' @param newdataLong,newdataEvent Optional data frames containing new data +#' (e.g. holdout data) to use when evaluating the log-likelihood for a +#' model estimated using \code{\link{stan_jm}}. If the fitted model #' was a multivariate joint model (i.e. more than one longitudinal outcome), -#' then \code{newdataLong} is allowed to be a list of data frames. If supplying +#' then \code{newdataLong} is allowed to be a list of data frames. If supplying #' new data, then \code{newdataEvent} should also include variables corresponding #' to the event time and event indicator as these are required for evaluating the -#' log likelihood for the event submodel. For more details, see the description +#' log likelihood for the event submodel. For more details, see the description #' of \code{newdataLong} and \code{newdataEvent} for \code{\link{posterior_survfit}}. -#' +#' log_lik.stanjm <- function(object, newdataLong = NULL, newdataEvent = NULL, ...) { if (!used.sampling(object)) STOP_sampling_only("Pointwise log-likelihood matrix") @@ -154,11 +153,11 @@ log_lik.stanjm <- function(object, newdataLong = NULL, newdataEvent = NULL, ...) stop("Both newdataLong and newdataEvent must be supplied together.") if (!is.null(newdataLong)) { newdatas <- validate_newdatas(object, newdataLong, newdataEvent) - newdataLong <- newdatas[1:M] + newdataLong <- newdatas[1:M] newdataEvent <- newdatas[["Event"]] } pars <- extract_pars(object) # full array of draws - data <- .pp_data_jm(object, newdataLong, newdataEvent) + data <- .pp_data_jm(object, newdataLong, newdataEvent) calling_fun <- as.character(sys.call(-1))[1] reloo_or_kfold <- calling_fun %in% c("kfold", "reloo") val <- .ll_jm(object, data, pars, reloo_or_kfold = reloo_or_kfold, ...) @@ -175,11 +174,11 @@ ll_fun <- function(x, m = NULL) { f <- family(x, m = m) if (!is(f, "family") || is_scobit(x)) return(.ll_polr_i) - else if (is_clogit(x)) + else if (is_clogit(x)) return(.ll_clogit_i) - else if (is.nlmer(x)) + else if (is.nlmer(x)) return(.ll_nlmer_i) - + fun <- paste0(".ll_", family(x, m = m)$family, "_i") get(fun, mode = "function") } @@ -191,8 +190,8 @@ ll_fun <- function(x, m = NULL) { # newdata is specified) # @param m Integer specifying which submodel for stanmvreg objects # @param reloo_or_kfold logical. TRUE if ll_args is for reloo or kfold -# @param ... For models without group-specific terms (i.e., not stan_[g]lmer), -# if reloo_or_kfold is TRUE and 'newdata' is specified then ... is used to +# @param ... For models without group-specific terms (i.e., not stan_[g]lmer), +# if reloo_or_kfold is TRUE and 'newdata' is specified then ... is used to # pass 'newx' and 'stanmat' from reloo or kfold (bypassing pp_data). This is a # workaround in case there are issues with newdata containing factors with # only a single level. Or for stanmvreg objects, then ... can be used to pass @@ -201,15 +200,15 @@ ll_fun <- function(x, m = NULL) { # @return a named list with elements data, draws, S (posterior sample size) and # N = number of observations ll_args <- function(object, ...) UseMethod("ll_args") -ll_args.stanreg <- function(object, newdata = NULL, offset = NULL, m = NULL, +ll_args.stanreg <- function(object, newdata = NULL, offset = NULL, m = NULL, reloo_or_kfold = FALSE, ...) { validate_stanreg_object(object) f <- family(object, m = m) draws <- nlist(f) has_newdata <- !is.null(newdata) - + dots <- list(...) - + z_betareg <- NULL if (has_newdata && reloo_or_kfold && !is.mer(object)) { x <- dots$newx @@ -236,8 +235,8 @@ ll_args.stanreg <- function(object, newdata = NULL, offset = NULL, m = NULL, } if (is.stanmvreg(object) && !is.null(dots$stanmat)) { stanmat <- dots$stanmat # potentially use a stanmat with a single draw - } - + } + if (!is_polr(object)) { # not polr or scobit model fname <- f$family if (is.nlmer(object)) { @@ -250,7 +249,7 @@ ll_args.stanreg <- function(object, newdata = NULL, offset = NULL, m = NULL, } data$i_ <- seq_len(nrow(data)) # for nlmer need access to i inside .ll_nlmer_i return(nlist(data, draws, S = NROW(draws$mu), N = nrow(data))) - + } else if (!is.binomial(fname)) { data <- data.frame(y, x) if (!is.null(z_betareg)) { @@ -270,26 +269,29 @@ ll_args.stanreg <- function(object, newdata = NULL, offset = NULL, m = NULL, trials <- 1L } else { trials <- 1 - if (is.factor(y)) + if (is.factor(y)) y <- fac2bin(y) - stopifnot(all(y %in% c(0, 1))) + if (!is_car(object)) + stopifnot(all(y %in% c(0, 1))) + else + trials <- object$trials } data <- data.frame(y, trials, x) } - nms <- if (is.stanmvreg(object)) + nms <- if (is.stanmvreg(object)) collect_nms(colnames(stanmat), - M = get_M(object), - stub = get_stub(object)) else NULL + M = get_M(object), + stub = get_stub(object)) else NULL beta_sel <- if (is.null(nms)) seq_len(ncol(x)) else nms$y[[m]] draws$beta <- stanmat[, beta_sel, drop = FALSE] m_stub <- get_m_stub(m, stub = get_stub(object)) - if (is.gaussian(fname)) + if (is.gaussian(fname)) draws$sigma <- stanmat[, paste0(m_stub, "sigma")] - if (is.gamma(fname)) + if (is.gamma(fname)) draws$shape <- stanmat[, paste0(m_stub, "shape")] - if (is.ig(fname)) + if (is.ig(fname)) draws$lambda <- stanmat[, paste0(m_stub, "lambda")] - if (is.nb(fname)) + if (is.nb(fname)) draws$size <- stanmat[, paste0(m_stub, "reciprocal_dispersion")] if (is.beta(fname)) { draws$f_phi <- object$family_phi @@ -299,7 +301,7 @@ ll_args.stanreg <- function(object, newdata = NULL, offset = NULL, m = NULL, } else { if (has_newdata) { if (!is.null(z_betareg)) { - colnames(data) <- c("y", colnames(get_x(object)), + colnames(data) <- c("y", colnames(get_x(object)), paste0("(phi)_", colnames(z_betareg))) } } else { @@ -323,21 +325,21 @@ ll_args.stanreg <- function(object, newdata = NULL, offset = NULL, m = NULL, zetas <- grep("|", colnames(stanmat), fixed = TRUE, value = TRUE) draws$zeta <- stanmat[, zetas, drop = FALSE] draws$max_y <- max(y) - if ("alpha" %in% colnames(stanmat)) { + if ("alpha" %in% colnames(stanmat)) { stopifnot(is_scobit(object)) # scobit draws$alpha <- stanmat[, "alpha"] draws$f <- object$method } } - + data$offset <- if (has_newdata) offset else object$offset if (model_has_weights(object)) { - if (is.stanmvreg(object)) + if (is.stanmvreg(object)) STOP_if_stanmvreg("posterior_survfit with weights") data$weights <- object$weights } - + if (is.mer(object)) { b_sel <- if (is.null(nms)) b_names(colnames(stanmat)) else nms$y_b[[m]] b <- stanmat[, b_sel, drop = FALSE] @@ -356,25 +358,30 @@ ll_args.stanreg <- function(object, newdata = NULL, offset = NULL, m = NULL, data <- cbind(data, as.matrix(z)[1:NROW(x),, drop = FALSE]) draws$beta <- cbind(draws$beta, b) } - + if (is_car(object)) { + psi_indx <- grep("^psi\\[[[:digit:]]+\\]", colnames(stanmat)) + psi <- stanmat[, psi_indx, drop = FALSE] + data$psi <- t(psi) + out <- nlist(data, draws, S = NROW(draws$beta), N = nrow(data)) + } if (is_clogit(object)) { data$strata <- strata out <- nlist(data, draws, S = NROW(draws$beta), N = nlevels(strata)) } else { - out <- nlist(data, draws, S = NROW(draws$beta), N = nrow(data)) + out <- nlist(data, draws, S = NROW(draws$beta), N = nrow(data)) } return(out) } # check intercept for polr models ----------------------------------------- -# Check if a model fit with stan_polr has an intercept (i.e. if it's actually a +# Check if a model fit with stan_polr has an intercept (i.e. if it's actually a # bernoulli model). If it doesn't have an intercept then the intercept column in # x is dropped. This is only necessary if newdata is specified because otherwise # the correct x is taken from the fitted model object. .validate_polr_x <- function(object, x) { x0 <- get_x(object) - has_intercept <- colnames(x0)[1L] == "(Intercept)" + has_intercept <- colnames(x0)[1L] == "(Intercept)" if (!has_intercept && colnames(x)[1L] == "(Intercept)") x <- x[, -1L, drop = FALSE] x @@ -387,20 +394,25 @@ ll_args.stanreg <- function(object, newdata = NULL, offset = NULL, m = NULL, val } else { val * w - } + } } .xdata <- function(data) { - sel <- c("y", "weights","offset", "trials","strata") + if (!is.null(data$psi)) + sel <- c("y", "weights","offset", "trials", "psi") + else + sel <- c("y", "weights","offset", "trials","strata") data[, -which(colnames(data) %in% sel)] } .mu <- function(data, draws) { eta <- as.vector(linear_predictor(draws$beta, .xdata(data), data$offset)) + if (!is.null(data$psi)) # for CAR models (adding spatial effect) + eta <- eta + as.vector(data$psi) draws$f$linkinv(eta) } # for stan_betareg only -.xdata_beta <- function(data) { +.xdata_beta <- function(data) { sel <- c("y", "weights","offset", "trials") data[, -c(which(colnames(data) %in% sel), grep("(phi)_", colnames(data), fixed = TRUE))] } @@ -440,15 +452,15 @@ ll_args.stanreg <- function(object, newdata = NULL, offset = NULL, m = NULL, .weighted(val, data_i$weights) } .ll_Gamma_i <- function(data_i, draws) { - val <- dgamma(data_i$y, shape = draws$shape, + val <- dgamma(data_i$y, shape = draws$shape, rate = draws$shape / .mu(data_i,draws), log = TRUE) .weighted(val, data_i$weights) } .ll_inverse.gaussian_i <- function(data_i, draws) { mu <- .mu(data_i, draws) - val <- 0.5 * log(draws$lambda / (2 * pi)) - + val <- 0.5 * log(draws$lambda / (2 * pi)) - 1.5 * log(data_i$y) - - 0.5 * draws$lambda * (data_i$y - mu)^2 / + 0.5 * draws$lambda * (data_i$y - mu)^2 / (data_i$y * mu^2) .weighted(val, data_i$weights) } @@ -464,7 +476,7 @@ ll_args.stanreg <- function(object, newdata = NULL, offset = NULL, m = NULL, } else if (y_i == J) { val <- log1p(-linkinv(draws$zeta[, J-1] - eta)) } else { - val <- log(linkinv(draws$zeta[, y_i] - eta) - + val <- log(linkinv(draws$zeta[, y_i] - eta) - linkinv(draws$zeta[, y_i - 1L] - eta)) } } else { @@ -505,7 +517,7 @@ ll_args.stanreg <- function(object, newdata = NULL, offset = NULL, m = NULL, # @param pars Output from extract_pars # @param m Integer specifying which submodel # @param reloo_or_kfold logical. TRUE if ll_args is for reloo or kfold -ll_args.stanjm <- function(object, data, pars, m = 1, +ll_args.stanjm <- function(object, data, pars, m = 1, reloo_or_kfold = FALSE, ...) { validate_stanjm_object(object) if (model_has_weights(object)) @@ -521,8 +533,8 @@ ll_args.stanjm <- function(object, data, pars, m = 1, y <- data$y[[m]] x <- data$yX[[m]] z <- t(data$yZt[[m]]) - Z_names <- data$yZ_names[[m]] - } else { + Z_names <- data$yZ_names[[m]] + } else { # for stan_mvmer models, log_lik is only ever called for # one submodel at a time, so data is for one submodel y <- data$y @@ -538,22 +550,22 @@ ll_args.stanjm <- function(object, data, pars, m = 1, y <- y[, 1L] } else { trials <- 1 - if (is.factor(y)) + if (is.factor(y)) y <- fac2bin(y) stopifnot(all(y %in% c(0, 1))) } dat <- data.frame(y, trials, x) - } + } dat <- cbind(dat, as.matrix(z)) draws$beta <- stanmat[, nms$y[[m]], drop = FALSE] m_stub <- get_m_stub(m) - if (is.gaussian(fname)) + if (is.gaussian(fname)) draws$sigma <- stanmat[, paste0(m_stub, "sigma")] - if (is.gamma(fname)) + if (is.gamma(fname)) draws$shape <- stanmat[, paste0(m_stub, "shape")] - if (is.ig(fname)) + if (is.ig(fname)) draws$lambda <- stanmat[, paste0(m_stub, "lambda")] - if (is.nb(fname)) + if (is.nb(fname)) draws$size <- stanmat[, paste0(m_stub, "reciprocal_dispersion")] b <- stanmat[, nms$y_b[[m]], drop = FALSE] b <- pp_b_ord(b, Z_names) @@ -567,9 +579,9 @@ ll_args.stanjm <- function(object, data, pars, m = 1, # with elements $basehaz, $family, $assoc # @param data Output from .pp_data_jm # @param pars Output from extract_pars -# @param include_long A logical, if TRUE then the log likelihood for the +# @param include_long A logical, if TRUE then the log likelihood for the # longitudinal submodels are included in the log likelihood calculation. -# @param include_b A logical, if TRUE then the log likelihood for the random +# @param include_b A logical, if TRUE then the log likelihood for the random # effects distribution is also included in the log likelihood calculation. # @param sum A logical. If TRUE then the log likelihood is summed across all # individuals. If FALSE then the log likelihood is returned for each @@ -579,20 +591,20 @@ ll_args.stanjm <- function(object, data, pars, m = 1, # a logical specifying whether the function calling ll_jm was reloo or kfold. # @return Either a matrix, a vector or a scalar, depending on the input types # and whether sum is set to TRUE. -.ll_jm <- function(object, data, pars, include_long = TRUE, +.ll_jm <- function(object, data, pars, include_long = TRUE, include_b = FALSE, sum = FALSE, ...) { - + M <- get_M(object) - + # Log likelihood for event submodel ll_event <- .ll_survival(object, data, pars) - + # Log likelihoods for longitudinal submodels if (include_long) { - ll_long <- lapply(1:M, function(m) + ll_long <- lapply(1:M, function(m) .ll_long(object, data, pars, m = m, ...)) } - + # Log likelihood for random effects submodel # NB this is only used in the Metropolis algorithm in 'posterior_survfit' # when drawing random effects for new individuals. But it is not used @@ -614,7 +626,7 @@ ll_args.stanjm <- function(object, data, pars, m = 1, b <- as.vector(pp_b_ord(b, Z_names)) Sigma_id <- VarCorr(object, stanmat = pars$stanmat)[[id_var]] if (length(cnms) > 1L) { - b2_var <- grep(utils::glob2rx(id_var), names(cnms), + b2_var <- grep(utils::glob2rx(id_var), names(cnms), value = TRUE, invert = TRUE) Sigma_b2 <- VarCorr(object, stanmat = pars$stanmat)[[b2_var]] Sigma_list <- rep(list(Sigma_b2), data$Ni) @@ -628,17 +640,17 @@ ll_args.stanjm <- function(object, data, pars, m = 1, } else { Sigma <- Sigma_id } - ll_b <- -0.5 * (c(determinant(Sigma, logarithm = TRUE)$modulus) + + ll_b <- -0.5 * (c(determinant(Sigma, logarithm = TRUE)$modulus) + (b %*% chol2inv(chol(Sigma)) %*% b)[1] + length(b) * log(2 * pi)) } else { ll_b <- NULL } - + # Check the dimensions of the various components if (is.matrix(ll_event)) { # S * Npat matrices if (include_long) { mats <- unique(sapply(c(ll_long, list(ll_event)), is.matrix)) - dims <- unique(lapply(c(ll_long, list(ll_event)), dim)) + dims <- unique(lapply(c(ll_long, list(ll_event)), dim)) if ((length(dims) > 1L) || (length(mats) > 1L)) stop("Bug found: elements of 'll_long' should be same class and ", "dimension as 'll_event'.") @@ -654,8 +666,8 @@ ll_args.stanjm <- function(object, data, pars, m = 1, } if (include_b && !identical(length(ll_b), length(ll_event))) stop("Bug found: length of 'll_b' should be equal to length of 'll_event'.") - } - + } + # Sum the various components (long + event + random effects) if (include_long) { val <- Reduce('+', c(ll_long, list(ll_event))) @@ -663,18 +675,18 @@ ll_args.stanjm <- function(object, data, pars, m = 1, val <- ll_event } if (include_b && is.matrix(val)) { - val <- sweep(val, 2L, ll_b, `+`) + val <- sweep(val, 2L, ll_b, `+`) } else if (include_b && is.vector(val)) { val <- val + ll_b } - + # Return log likelihood for joint model if (!sum) { return(val) # S * Npat matrix or length Npat vector } else if (is.matrix(val)) { return(rowSums(val)) # length S vector } else { - return(sum(val)) # scalar + return(sum(val)) # scalar } } @@ -684,11 +696,11 @@ ll_args.stanjm <- function(object, data, pars, m = 1, # @param data Output from .pp_data_jm. # @param pars Output from extract_pars. # @param m Integer specifying the longitudinal submodel. -# @param reloo_or_kfold Logical specifying whether the call came from +# @param reloo_or_kfold Logical specifying whether the call came from # reloo or kfold. # @return An S*Npat matrix. .ll_long <- function(object, data, pars, m = 1, reloo_or_kfold = FALSE) { - args <- ll_args.stanjm(object, data, pars, m = m, + args <- ll_args.stanjm(object, data, pars, m = m, reloo_or_kfold = reloo_or_kfold) fun <- ll_fun(object, m = m) ll <- lapply(seq_len(args$N), function(j) as.vector( @@ -697,7 +709,7 @@ ll_args.stanjm <- function(object, data, pars, m = 1, # return S*Npat matrix by summing log-lik for y within each individual res <- apply(ll, 1L, function(row) tapply(row, data$flist[[m]], sum)) res <- if (is.vector(res) & (args$S > 1L)) cbind(res) else t(res) - return(res) + return(res) } # Return survival probability or log-likelihood for event submodel @@ -705,34 +717,34 @@ ll_args.stanjm <- function(object, data, pars, m = 1, # @param object A stanjm object. # @param data Output from .pp_data_jm. # @param pars Output from extract_pars. -# @param one_draw A logical specifying whether the parameters provided in the +# @param one_draw A logical specifying whether the parameters provided in the # pars argument are vectors for a single realisation of the parameter (e.g. # a single MCMC draw, or a posterior mean) (TRUE) or a stanmat array (FALSE). -# @param survprob A logical specifying whether to return the survival probability +# @param survprob A logical specifying whether to return the survival probability # (TRUE) or the log likelihood for the event submodel (FALSE). # @param An S by Npat matrix, or a length Npat vector, depending on the inputs -# (where S is the size of the posterior sample and Npat is the number of +# (where S is the size of the posterior sample and Npat is the number of # individuals). .ll_survival <- function(object, data, pars, one_draw = FALSE, survprob = FALSE) { basehaz <- object$basehaz family <- object$family - assoc <- object$assoc + assoc <- object$assoc etimes <- attr(data$assoc_parts, "etimes") estatus <- attr(data$assoc_parts, "estatus") qnodes <- attr(data$assoc_parts, "qnodes") qtimes <- attr(data$assoc_parts, "qtimes") - qwts <- attr(data$assoc_parts, "qwts") + qwts <- attr(data$assoc_parts, "qwts") times <- c(etimes, qtimes) - - # To avoid an error in log(times) replace times equal to zero with a small + + # To avoid an error in log(times) replace times equal to zero with a small # non-zero value. Note that these times correspond to individuals where the, - # event time (etimes) was zero, and therefore the cumhaz (at baseline) will - # be forced to zero for these individuals further down in the code anyhow. - times[times == 0] <- 1E-10 - + # event time (etimes) was zero, and therefore the cumhaz (at baseline) will + # be forced to zero for these individuals further down in the code anyhow. + times[times == 0] <- 1E-10 + # Linear predictor for the survival submodel - e_eta <- linear_predictor(pars$ebeta, data$eXq) - + e_eta <- linear_predictor(pars$ebeta, data$eXq) + # Add on contribution from assoc structure if (length(pars$abeta)) { M <- get_M(object) @@ -747,26 +759,26 @@ ll_args.stanjm <- function(object, data, pars, m = 1, pp_b_ord(if (is.matrix(b_m)) b_m else t(b_m), Z_names_m) }) if (one_draw) { - aXq <- make_assoc_terms(parts = data$assoc_parts, assoc = assoc, + aXq <- make_assoc_terms(parts = data$assoc_parts, assoc = assoc, family = family, beta = pars$beta, b = pars$b) e_eta <- e_eta + linear_predictor.default(pars$abeta, aXq) } else { - aXq <- make_assoc_terms(parts = data$assoc_parts, assoc = assoc, + aXq <- make_assoc_terms(parts = data$assoc_parts, assoc = assoc, family = family, beta = pars$beta, b = pars$b) for (k in 1:length(aXq)) { e_eta <- e_eta + sweep(aXq[[k]], 1L, pars$abeta[,k], `*`) } - } + } } - + # Log baseline hazard at etimes (if not NULL) and qtimes - log_basehaz <- evaluate_log_basehaz(times = times, - basehaz = basehaz, + log_basehaz <- evaluate_log_basehaz(times = times, + basehaz = basehaz, coefs = pars$bhcoef) - + # Log hazard at etimes (if not NULL) and qtimes - log_haz <- log_basehaz + e_eta - + log_haz <- log_basehaz + e_eta + # Extract log hazard at qtimes only if (is.vector(log_haz)) { q_log_haz <- tail(log_haz, length(qtimes)) @@ -774,19 +786,19 @@ ll_args.stanjm <- function(object, data, pars, m = 1, sel_cols <- tail(1:ncol(log_haz), length(qtimes)) q_log_haz <- log_haz[, sel_cols, drop = FALSE] } - + # Evaluate log survival log_surv <- evaluate_log_survival(log_haz = q_log_haz, qnodes = qnodes, qwts = qwts) - + # Force surv prob to 1 (ie. log surv prob to 0) if evaluating # at time t = 0; this avoids possible numerical errors log_surv[etimes == 0] <- 0 - + # Possibly return surv prob at time t (upper limit of integral) if (survprob) - return(exp(log_surv)) - + return(exp(log_surv)) + # Otherwise return log likelihood at time t if (is.null(etimes) || is.null(estatus)) stop("'etimes' and 'estatus' cannot be NULL if 'survprob = FALSE'.") @@ -802,7 +814,7 @@ ll_args.stanjm <- function(object, data, pars, m = 1, e_log_haz <- log_haz[, 1:length(etimes), drop = FALSE] return(sweep(e_log_haz, 2L, estatus, `*`) + log_surv) } -} +} # Evaluate the log baseline hazard at the specified times # given the vector or matrix of MCMC draws for the baseline @@ -814,12 +826,12 @@ ll_args.stanjm <- function(object, data, pars, m = 1, # @return A vector or matrix, depending on the input type of coefs. evaluate_log_basehaz <- function(times, basehaz, coefs) { type <- basehaz$type_name - if (type == "weibull") { + if (type == "weibull") { X <- log(times) # log times B1 <- log(coefs) # log shape B2 <- coefs - 1 # shape - 1 log_basehaz <- as.vector(B1) + linear_predictor(B2,X) - } else if (type == "bs") { + } else if (type == "bs") { X <- predict(basehaz$bs_basis, times) # b-spline basis B <- coefs # b-spline coefs log_basehaz <- linear_predictor(B,X) @@ -837,7 +849,7 @@ evaluate_log_basehaz <- function(times, basehaz, coefs) { # individual, evaluated at each of the quadrature points. The # vector should be ordered such that the first N elements contain # the log_haz evaluated for each individual at quadrature point 1, -# then the next N elements are the log_haz evaluated for each +# then the next N elements are the log_haz evaluated for each # individual at quadrature point 2, and so on. # @param qnodes Integer specifying the number of quadrature nodes # at which the log hazard was evaluated for each individual. @@ -851,7 +863,7 @@ evaluate_log_survival.default <- function(log_haz, qnodes, qwts) { # convert log hazard to hazard haz <- exp(log_haz) # apply GK quadrature weights - weighted_haz <- qwts * haz + weighted_haz <- qwts * haz # sum quadrature points for each individual to get cum_haz splitting_vec <- rep(1:qnodes, each = length(haz) / qnodes) cum_haz <- Reduce('+', split(weighted_haz, splitting_vec)) @@ -868,4 +880,4 @@ evaluate_log_survival.matrix <- function(log_haz, qnodes, qwts) { cum_haz <- Reduce('+', array2list(weighted_haz, nsplits = qnodes)) # return: -cum_haz == log survival probability -cum_haz -} +} diff --git a/R/loo.R b/R/loo.R index 9b1f4e0d7..0e839c8e1 100644 --- a/R/loo.R +++ b/R/loo.R @@ -406,9 +406,6 @@ loo_compare.stanreg <- .loo_comparison(fits, criterion = criterion, detail = detail) } - -#' @rdname loo.stanreg -#' @export loo_compare.stanreg_list <- function(x, ..., @@ -722,7 +719,7 @@ kfold_and_reloo_data <- function(x) { } d <- d[, all_vars, drop=FALSE] } - + if (!is.null(sub)) { d <- d[keep,, drop=FALSE] } @@ -733,7 +730,7 @@ kfold_and_reloo_data <- function(x) { strata_var <- as.character(getCall(x)$strata) d[[strata_var]] <- model.weights(model.frame(x)) } - + return(d) } diff --git a/R/misc.R b/R/misc.R index 71a26ed7f..1687239e7 100644 --- a/R/misc.R +++ b/R/misc.R @@ -135,6 +135,11 @@ is_polr <- function(object) { is(object, "polr") } +# test if a stanreg object has class car (conditional autoregressive) +is_car <- function(object) { + is(object, "car") +} + # test if a stanreg object is a scobit model is_scobit <- function(object) { validate_stanreg_object(object) diff --git a/R/posterior_predict.R b/R/posterior_predict.R index 0b4e4c2b2..fe33a772a 100644 --- a/R/posterior_predict.R +++ b/R/posterior_predict.R @@ -62,10 +62,10 @@ #' model. #' @param ... For \code{stanmvreg} objects, argument \code{m} can be specified #' indicating the submodel for which you wish to obtain predictions. -#' +#' #' @return A \code{draws} by \code{nrow(newdata)} matrix of simulations from the -#' posterior predictive distribution. Each row of the matrix is a vector of -#' predictions generated using a single draw of the model parameters from the +#' posterior predictive distribution. Each row of the matrix is a vector of +#' predictions generated using a single draw of the model parameters from the #' posterior distribution. The returned matrix will also have class #' \code{"ppd"} to indicate it contains draws from the posterior predictive #' distribution. @@ -83,13 +83,13 @@ #' probably with \code{successes} set to \code{0} and \code{trials} specifying #' the number of trials. See the Examples section below and the #' \emph{How to Use the rstanarm Package} for examples. -#' @note For models estimated with \code{\link{stan_clogit}}, the number of +#' @note For models estimated with \code{\link{stan_clogit}}, the number of #' successes per stratum is ostensibly fixed by the research design. Thus, when #' doing posterior prediction with new data, the \code{data.frame} passed to #' the \code{newdata} argument must contain an outcome variable and a stratifying -#' factor, both with the same name as in the original \code{data.frame}. Then, the +#' factor, both with the same name as in the original \code{data.frame}. Then, the #' posterior predictions will condition on this outcome in the new data. -#' +#' #' @seealso \code{\link{pp_check}} for graphical posterior predictive checks. #' Examples of posterior predictive checking can also be found in the #' \pkg{rstanarm} vignettes and demos. @@ -161,15 +161,14 @@ posterior_predict.stanreg <- function(object, newdata = NULL, draws = NULL, if (is.stanmvreg(object)) { m <- dots[["m"]] # submodel to predict for stanmat <- dots[["stanmat"]] # possibly incl. new b pars (dynamic preds) - if (is.null(m)) + if (is.null(m)) STOP_arg_required_for_stanmvreg(m) if (!is.null(offset)) stop2("'offset' cannot be specified for stanmvreg objects.") } else { m <- NULL stanmat <- NULL - } - + } newdata <- validate_newdata(object, newdata = newdata, m = m) pp_data_args <- c(list(object, newdata = newdata, @@ -196,17 +195,17 @@ posterior_predict.stanreg <- function(object, newdata = NULL, draws = NULL, } else if (is.stanjm(object)) { ppargs <- pp_args(object, data = pp_eta(object, dat, draws, m = m, stanmat = stanmat), m = m) - + } else { if (!is.null(newdata) && is_clogit(object)) { y <- eval(formula(object)[[2L]], newdata) strata <- as.factor(eval(object$call$strata, newdata)) formals(object$family$linkinv)$g <- strata - formals(object$family$linkinv)$successes <- + formals(object$family$linkinv)$successes <- aggregate(y, by = list(strata), FUN = sum)$x } ppargs <- pp_args(object, data = pp_eta(object, dat, draws, m = m), m = m) - } + } if (is_clogit(object)) { if (is.null(newdata)) ppargs$strata <- model.frame(object)[,"(weights)"] @@ -215,7 +214,8 @@ posterior_predict.stanreg <- function(object, newdata = NULL, draws = NULL, } else if (!is_polr(object) && is.binomial(family(object, m = m)$family)) { ppargs$trials <- pp_binomial_trials(object, newdata, m = m) } - + if (is_car(object) && is.binomial(family(object)$family)) + ppargs$trials <- object$trials ppfun <- pp_fun(object, m = m) ytilde <- do.call(ppfun, ppargs) @@ -227,15 +227,22 @@ posterior_predict.stanreg <- function(object, newdata = NULL, draws = NULL, ytilde <- do.call(fun, list(ytilde)) if (is_polr(object) && !is_scobit(object)) ytilde <- matrix(levels(get_y(object))[ytilde], nrow(ytilde), ncol(ytilde)) - - if (is.null(newdata)) colnames(ytilde) <- rownames(model.frame(object, m = m)) - else colnames(ytilde) <- rownames(newdata) - + + # fix me! + if (is_car(object)) { + if (is.null(newdata)) colnames(ytilde) <- rownames(model.frame(object)) + else colnames(ytilde) <- rownames(newdata) + } + else { + if (is.null(newdata)) colnames(ytilde) <- rownames(model.frame(object, m = m)) + else colnames(ytilde) <- rownames(newdata) + } + # if function is called from posterior_traj then add mu as attribute fn <- tryCatch(sys.call(-3)[[1]], error = function(e) NULL) if (!is.null(fn) && grepl("posterior_traj", deparse(fn), fixed = TRUE)) return(structure(ytilde, mu = ppargs$mu, class = c("ppd", class(ytilde)))) - + structure(ytilde, class = c("ppd", class(ytilde))) } @@ -243,7 +250,7 @@ posterior_predict.stanreg <- function(object, newdata = NULL, draws = NULL, #' @export #' @templateVar mArg m #' @template args-m -#' +#' posterior_predict.stanmvreg <- function(object, m = 1, newdata = NULL, draws = NULL, re.form = NULL, fun = NULL, seed = NULL, ...) { validate_stanmvreg_object(object) @@ -256,15 +263,15 @@ posterior_predict.stanmvreg <- function(object, m = 1, newdata = NULL, draws = N re.form = re.form, fun = fun, seed = seed, offset = NULL, m = m, ...) out -} - - +} + + # internal ---------------------------------------------------------------- # functions to draw from the various posterior predictive distributions pp_fun <- function(object, m = NULL) { - suffix <- if (is_polr(object)) "polr" else - if (is_clogit(object)) "clogit" else + suffix <- if (is_polr(object)) "polr" else + if (is_clogit(object)) "clogit" else family(object, m = m)$family get(paste0(".pp_", suffix), mode = "function") @@ -394,7 +401,7 @@ pp_args <- function(object, data, m = NULL) { # @param m Optional integer specifying the submodel for stanmvreg objects # @param stanmat Optionally pass a stanmat that has been amended to include # new b parameters for individuals in the prediction data but who were not -# included in the model estimation; relevant for dynamic predictions for +# included in the model estimation; relevant for dynamic predictions for # stan_jm objects only # @return Linear predictor "eta" and matrix of posterior draws "stanmat". For # some stan_betareg models "" is also included in the list. @@ -416,11 +423,11 @@ pp_eta <- function(object, data, draws = NULL, m = NULL, stanmat = NULL) { M <- get_M(object) } if (is.null(stanmat)) { - stanmat <- if (is.null(data$Zt)) + stanmat <- if (is.null(data$Zt)) as.matrix.stanreg(object) else as.matrix(object$stanfit) } - nms <- if (is.stanmvreg(object)) - collect_nms(colnames(stanmat), M, stub = get_stub(object)) else NULL + nms <- if (is.stanmvreg(object)) + collect_nms(colnames(stanmat), M, stub = get_stub(object)) else NULL beta_sel <- if (is.null(nms)) seq_len(ncol(x)) else nms$y[[m]] beta <- stanmat[, beta_sel, drop = FALSE] if (some_draws) @@ -443,19 +450,28 @@ pp_eta <- function(object, data, draws = NULL, m = NULL, stanmat = NULL) { else eta <- linkinv(object)(eta, data$arg1, data$arg2) eta <- t(eta) } - + + if (is_car(object)) { + psi_indx <- grep("psi", colnames(stanmat)) + psi <- stanmat[, psi_indx, drop = FALSE] + if (some_draws) + psi <- psi[samp, , drop = FALSE] + eta <- eta + psi + } +# nlist(eta, stanmat) + out <- nlist(eta, stanmat) - + if (inherits(object, "betareg")) { z_vars <- colnames(stanmat)[grepl("(phi)", colnames(stanmat))] omega <- stanmat[, z_vars] if (length(z_vars) == 1 && z_vars == "(phi)") { - out$phi <- stanmat[, "(phi)"] + out$phi <- stanmat[, "(phi)"] } else { out$phi_linpred <- linear_predictor(as.matrix(omega), as.matrix(data$z_betareg), data$offset) } } - + return(out) } @@ -501,15 +517,15 @@ pp_binomial_trials <- function(object, newdata = NULL, m = NULL) { if (is.stanmvreg(object) && is.null(m)) { STOP_arg_required_for_stanmvreg(m) } - - y <- get_y(object, m) + + y <- get_y(object, m) is_bernoulli <- NCOL(y) == 1L - + if (is_bernoulli) { - trials <- if (is.null(newdata)) + trials <- if (is.null(newdata)) rep(1, NROW(y)) else rep(1, NROW(newdata)) } else { - trials <- if (is.null(newdata)) + trials <- if (is.null(newdata)) rowSums(y) else rowSums(eval(formula(object, m = m)[[2L]], newdata)) } return(trials) diff --git a/R/pp_data.R b/R/pp_data.R index a2b2e5bd6..25e4b0f35 100644 --- a/R/pp_data.R +++ b/R/pp_data.R @@ -1,17 +1,17 @@ # Part of the rstanarm package for estimating model parameters # Copyright 2015 Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker # 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. @@ -41,7 +41,7 @@ pp_data <- requireNamespace("mgcv", quietly = TRUE) if (is.null(newdata)) x <- predict(object$jam, type = "lpmatrix") else x <- predict(object$jam, newdata = newdata, type = "lpmatrix") - if (is.null(offset)) + if (is.null(offset)) offset <- object$offset %ORifNULL% rep(0, nrow(x)) return(nlist(x, offset)) } @@ -53,27 +53,29 @@ pp_data <- if (inherits(object, "betareg")) { return(nlist(x, offset, z_betareg = object$z)) } - + return(nlist(x, offset)) } - + #if (!is.null(newdata) & is(object, "car")) { + # return(nlist(x)) + #} offset <- .pp_data_offset(object, newdata, offset) Terms <- delete.response(terms(object)) m <- model.frame(Terms, newdata, xlev = object$xlevels) - if (!is.null(cl <- attr(Terms, "dataClasses"))) + if (!is.null(cl <- attr(Terms, "dataClasses"))) .checkMFClasses(cl, m) x <- model.matrix(Terms, m, contrasts.arg = object$contrasts) - if (is(object, "polr") && !is_scobit(object)) + if (is(object, "polr") && !is_scobit(object)) x <- x[,colnames(x) != "(Intercept)", drop = FALSE] - + if (inherits(object, "betareg")) { - mf <- model.frame(delete.response(object$terms$precision), - data = newdata, na.action = object$na.action, + mf <- model.frame(delete.response(object$terms$precision), + data = newdata, na.action = object$na.action, xlev = object$levels$precision) z_betareg <- model.matrix(object$terms$precision, mf, contrasts = object$contrasts$precision) return(nlist(x, offset, z_betareg)) } - + return(nlist(x, offset)) } @@ -134,19 +136,19 @@ pp_data <- offset <- .pp_data_offset(object, newdata, offset) group <- with(nlf$reTrms, pad_reTrms(Ztlist, cnms, flist)) - if (!is.null(re.form) && !is(re.form, "formula") && is.na(re.form)) + if (!is.null(re.form) && !is(re.form, "formula") && is.na(re.form)) group$Z@x <- 0 return(nlist(x = nlf$X, offset = offset, Z = group$Z, Z_names = make_b_nms(group), arg1, arg2)) } -# the functions below are heavily based on a combination of -# lme4:::predict.merMod and lme4:::mkNewReTrms, although they do also have +# the functions below are heavily based on a combination of +# lme4:::predict.merMod and lme4:::mkNewReTrms, although they do also have # substantial modifications .pp_data_mer_x <- function(object, newdata, m = NULL, ...) { x <- get_x(object, m = m) if (is.null(newdata)) return(x) - form <- if (is.null(m)) attr(object$glmod$fr, "formula") else + form <- if (is.null(m)) attr(object$glmod$fr, "formula") else formula(object, m = m) L <- length(form) form[[L]] <- lme4::nobars(form[[L]]) @@ -158,7 +160,7 @@ pp_data <- mf <- mf[vars] isFac <- vapply(mf, is.factor, FUN.VALUE = TRUE) isFac[attr(Terms, "response")] <- FALSE - orig_levs <- if (length(isFac) == 0) + orig_levs <- if (length(isFac) == 0) NULL else lapply(mf[isFac], levels) mfnew <- model.frame(delete.response(Terms), newdata, xlev = orig_levs) x <- model.matrix(RHS, data = mfnew, contrasts.arg = attr(x, "contrasts")) @@ -166,22 +168,22 @@ pp_data <- } .pp_data_mer_z <- function(object, newdata, re.form = NULL, - allow.new.levels = TRUE, na.action = na.pass, + allow.new.levels = TRUE, na.action = na.pass, m = NULL, ...) { NAcheck <- !is.null(re.form) && !is(re.form, "formula") && is.na(re.form) - fmla0check <- (is(re.form, "formula") && - length(re.form) == 2 && + fmla0check <- (is(re.form, "formula") && + length(re.form) == 2 && identical(re.form[[2]], 0)) if (NAcheck || fmla0check) return(list()) if (is.null(newdata) && is.null(re.form)) { - Z <- get_z(object, m = m) + Z <- get_z(object, m = m) if (!is.stanmvreg(object)) { # Z_names not needed for stanreg with no newdata return(list(Zt = t(Z))) } else { # must supply Z_names for stanmvreg since b pars # might be for multiple submodels and Zt will only - # be for one submodel, so their elements may not + # be for one submodel, so their elements may not # correspond exactly ReTrms <- object$glmod[[m]]$reTrms Z_names <- make_b_nms(ReTrms, m = m, stub = get_stub(object)) @@ -190,7 +192,7 @@ pp_data <- } else if (is.null(newdata)) { rfd <- mfnew <- model.frame(object, m = m) - } + } else if (inherits(object, "gamm4")) { requireNamespace("mgcv", quietly = TRUE) if (is.null(newdata)) x <- predict(object$jam, type = "lpmatrix") @@ -210,7 +212,7 @@ pp_data <- if (!is.null(fixed.na.action)) attr(rfd,"na.action") <- fixed.na.action } - if (is.null(re.form)) + if (is.null(re.form)) re.form <- justRE(formula(object, m = m)) if (!inherits(re.form, "formula")) stop("'re.form' must be NULL, NA, or a formula.") @@ -242,7 +244,7 @@ null_or_zero <- function(x) { .pp_data_offset <- function(object, newdata = NULL, offset = NULL) { if (is.null(newdata)) { # get offset from model object (should be null if no offset) - if (is.null(offset)) + if (is.null(offset)) offset <- object$offset %ORifNULL% model.offset(model.frame(object)) } else { if (!is.null(offset)) @@ -250,12 +252,12 @@ null_or_zero <- function(x) { else { # if newdata specified but not offset then confirm that model wasn't fit # with an offset (warning, not error) - if (!is.null(object$call$offset) || - !null_or_zero(object$offset) || + if (!is.null(object$call$offset) || + !null_or_zero(object$offset) || !null_or_zero(model.offset(model.frame(object)))) { warning( - "'offset' argument is NULL but it looks like you estimated ", - "the model using an offset term.", + "'offset' argument is NULL but it looks like you estimated ", + "the model using an offset term.", call. = FALSE ) } @@ -272,82 +274,82 @@ null_or_zero <- function(x) { # log-likelihood in post-estimation functions for a \code{stan_jm} model # # @param object A stanmvreg object -# @param newdataLong A data frame or list of data frames with the new +# @param newdataLong A data frame or list of data frames with the new # covariate data for the longitudinal submodel # @param newdataEvent A data frame with the new covariate data for the # event submodel # @param ids An optional vector of subject IDs specifying which individuals # should be included in the returned design matrices. # @param etimes An optional vector of times at which the event submodel -# design matrices should be evaluated (also used to determine the +# design matrices should be evaluated (also used to determine the # quadrature times). If NULL then times are taken to be the eventimes in # the fitted object (if newdataEvent is NULL) or in newdataEvent. # @param long_parts,event_parts A logical specifying whether to return the # design matrices for the longitudinal and/or event submodels. -# @return A named list (with components M, Npat, ndL, ndE, yX, tZt, -# yZnames, eXq, assoc_parts) -.pp_data_jm <- function(object, newdataLong = NULL, newdataEvent = NULL, - ids = NULL, etimes = NULL, long_parts = TRUE, +# @return A named list (with components M, Npat, ndL, ndE, yX, tZt, +# yZnames, eXq, assoc_parts) +.pp_data_jm <- function(object, newdataLong = NULL, newdataEvent = NULL, + ids = NULL, etimes = NULL, long_parts = TRUE, event_parts = TRUE) { M <- get_M(object) id_var <- object$id_var time_var <- object$time_var - + if (!is.null(newdataLong) || !is.null(newdataEvent)) newdatas <- validate_newdatas(object, newdataLong, newdataEvent) - + # prediction data for longitudinal submodels - ndL <- if (is.null(newdataLong)) + ndL <- if (is.null(newdataLong)) get_model_data(object)[1:M] else newdatas[1:M] - + # prediction data for event submodel - ndE <- if (is.null(newdataEvent)) - get_model_data(object)[["Event"]] else newdatas[["Event"]] - + ndE <- if (is.null(newdataEvent)) + get_model_data(object)[["Event"]] else newdatas[["Event"]] + # possibly subset if (!is.null(ids)) { ndL <- subset_ids(object, ndL, ids) ndE <- subset_ids(object, ndE, ids) } id_list <- unique(ndE[[id_var]]) # unique subject id list - + # evaluate the last known survival time and status if (!is.null(newdataEvent) && is.null(etimes)) { - # prediction data for the event submodel was provided but + # prediction data for the event submodel was provided but # no event times were explicitly specified by the user, so # they must be evaluated using the data frame surv <- eval(formula(object, m = "Event")[[2L]], ndE) etimes <- unclass(surv)[,"time"] - estatus <- unclass(surv)[,"status"] + estatus <- unclass(surv)[,"status"] } else if (is.null(etimes)) { - # if no prediction data was provided then event times are + # if no prediction data was provided then event times are # taken from the fitted model etimes <- object$eventtime[as.character(id_list)] estatus <- object$status[as.character(id_list)] - } else { - # otherwise, event times ('etimes') are only directly specified for dynamic - # predictions via posterior_survfit in which case the 'etimes' correspond + } else { + # otherwise, event times ('etimes') are only directly specified for dynamic + # predictions via posterior_survfit in which case the 'etimes' correspond # to the last known survival time and therefore we assume everyone has survived - # up to that point (ie, set estatus = 0 for all individuals), this is true + # up to that point (ie, set estatus = 0 for all individuals), this is true # even if there is an event indicated in the data supplied by the user. estatus <- rep(0, length(etimes)) } res <- nlist(M, Npat = length(id_list), ndL, ndE) - + if (long_parts && event_parts) lapply(ndL, function(x) { - if (!time_var %in% colnames(x)) + if (!time_var %in% colnames(x)) STOP_no_var(time_var) - if (!id_var %in% colnames(x)) + if (!id_var %in% colnames(x)) STOP_no_var(id_var) if (any(x[[time_var]] < 0)) stop2("Values for the time variable (", time_var, ") should not be negative.") mt <- tapply(x[[time_var]], factor(x[[id_var]]), max) if (any(mt > etimes)) stop2("There appears to be observation times in the longitudinal data that ", - "are later than the event time specified in the 'etimes' argument.") - }) - + "are later than the event time specified in the 'etimes' argument.") + }) + # response and design matrices for longitudinal submodels if (long_parts) { y <- lapply(1:M, function(m) eval(formula(object, m = m)[[2L]], ndL[[m]])) @@ -358,7 +360,7 @@ null_or_zero <- function(x) { flist <- lapply(ndL, function(x) factor(x[[id_var]])) res <- c(res, nlist(y, yX, yZt, yZ_names, flist)) } - + # design matrices for event submodel and association structure if (event_parts) { qnodes <- object$qnodes @@ -376,13 +378,13 @@ null_or_zero <- function(x) { grp_stuff <- object$grp_stuff[[m]] if (grp_stuff$has_grp) { grp_stuff <- get_extra_grp_info( # update grp_info with new data - grp_stuff, flist = ymf, id_var = id_var, + grp_stuff, flist = ymf, id_var = id_var, grp_assoc = grp_stuff$grp_assoc) } ymf <- prepare_data_table(ymf, id_var = id_var, time_var = time_var, grp_var = grp_stuff$grp_var) make_assoc_parts( - ymf, assoc = object$assoc[,m], id_var = id_var, time_var = time_var, + ymf, assoc = object$assoc[,m], id_var = id_var, time_var = time_var, ids = id_rep, times = times, grp_stuff = grp_stuff, use_function = pp_data, object = object, m = m) }) @@ -390,7 +392,7 @@ null_or_zero <- function(x) { assoc_parts <- do.call("structure", assoc_attr) res <- c(res, nlist(eXq, assoc_parts)) } - + return(res) } @@ -398,13 +400,13 @@ null_or_zero <- function(x) { # (1) only includes variables used in the model formula # (2) only includes rows contained in the glmod/coxmod model frames # (3) ensures that additional variables that are required -# such as the ID variable or variables used in the +# such as the ID variable or variables used in the # interaction-type association structures, are included. # -# It is necessary to drop unneeded variables though so that -# errors are not encountered if the original data contained +# It is necessary to drop unneeded variables though so that +# errors are not encountered if the original data contained # NA values for variables unrelated to the model formula. -# We generate a data frame here for in-sample predictions +# We generate a data frame here for in-sample predictions # rather than using a model frame, since some quantities will # need to be recalculated at quadrature points etc, for example # in posterior_survfit. @@ -418,25 +420,25 @@ get_model_data <- function(object, m = NULL) { validate_stanmvreg_object(object) M <- get_M(object) terms <- terms(object, fixed.only = FALSE) - + # identify variables to add to the terms objects if (is.jm(object)) { extra_vars <- lapply(1:M, function(m) { - # for each submodel loop over the four possible assoc + # for each submodel loop over the four possible assoc # interaction formulas and collect any variables used forms_m <- object$assoc["which_formulas",][[m]] uapply(forms_m, function(x) { if (length(x)) { - rownames(attr(terms.formula(x), "factors")) + rownames(attr(terms.formula(x), "factors")) } else NULL }) }) # also ensure that id_var is in the event data extra_vars$Event <- object$id_var - + if (!identical(length(terms), length(extra_vars))) stop2("Bug found: terms and extra_vars should be same length.") - + # add the extra variables to the terms formula for each submodel terms <- xapply(terms, extra_vars, FUN = function(x, y) { lhs <- x[[2L]] @@ -445,20 +447,20 @@ get_model_data <- function(object, m = NULL) { rhs <- c(rhs, y) reformulate(rhs, response = lhs) }) - + datas <- c(object$dataLong, list(object$dataEvent)) } else { datas <- object$data } - + # identify rows that were in the model frame row_nms <- lapply(model.frame(object), rownames) - + # drop rows and variables not required for predictions mfs <- xapply(w = terms, x = datas, y = row_nms, - FUN = function(w, x, y) + FUN = function(w, x, y) get_all_vars(w, x)[y, , drop = FALSE]) - + mfs <- list_nms(mfs, M, stub = get_stub(object)) if (is.null(m)) mfs else mfs[[m]] } diff --git a/R/pp_validate.R b/R/pp_validate.R index ab3fe09a5..640f50df7 100644 --- a/R/pp_validate.R +++ b/R/pp_validate.R @@ -135,6 +135,7 @@ pp_validate <- function(object, nreps = 20, seed = 12345, ...) { post_mat <- as.matrix(post) data_mat <- posterior_predict(post) constant <- apply(data_mat, 1, FUN = function(x) all(duplicated(x)[-1L])) + if (any(constant)) stop("'pp_validate' cannot proceed because some simulated outcomes are constant. ", "Try again with better priors on the parameters.") diff --git a/R/predict.R b/R/predict.R index ae5386d6f..6f0002e07 100644 --- a/R/predict.R +++ b/R/predict.R @@ -72,6 +72,14 @@ predict.stanreg <- function(object, stanmat <- as.matrix.stanreg(object) beta <- stanmat[, seq_len(ncol(dat$x))] eta <- linear_predictor(beta, dat$x, dat$offset) + if (is(object, "car")) { + if (nrow(object$x) > nrow(dat$x)) + stop("'newdata' is less than the number of spatial regions.") + psi_indx <- grep("psi", colnames(as.matrix(object$stanfit))) + psi <- as.matrix(object$stanfit)[,psi_indx] + psi <- unname(colMeans(psi)) + eta <- eta + psi + } if (type == "response") { inverse_link <- linkinv(object) eta <- inverse_link(eta) diff --git a/R/print-and-summary.R b/R/print-and-summary.R index bfd6f9e82..827200e39 100644 --- a/R/print-and-summary.R +++ b/R/print-and-summary.R @@ -1,27 +1,27 @@ # 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. #' Print method for stanreg objects -#' +#' #' The \code{print} method for stanreg objects displays a compact summary of the #' fitted model. See the \strong{Details} section below for descriptions of the #' different components of the printed output. For additional summary statistics #' and diagnostics use the \code{\link[=summary.stanreg]{summary}} method. -#' +#' #' @export #' @method print stanreg #' @templateVar stanregArg x @@ -31,9 +31,9 @@ #' @param digits Number of digits to use for formatting numbers. #' @param ... Ignored. #' @return Returns \code{x}, invisibly. -#' @details +#' @details #' \subsection{Point estimates}{ -#' Regardless of the estimation algorithm, point estimates are medians computed +#' Regardless of the estimation algorithm, point estimates are medians computed #' from simulations. For models fit using MCMC (\code{"sampling"}) the posterior #' sample is used. For optimization (\code{"optimizing"}), the simulations are #' generated from the asymptotic Gaussian sampling distribution of the @@ -56,15 +56,15 @@ #' output also shows point estimates of the standard deviations of the group #' effects (and correlations if there are both intercept and slopes that vary by #' group). -#' +#' #' \item For analysis of variance models (see \code{\link{stan_aov}}) models, an #' ANOVA-like table is also displayed. -#' +#' #' \item For joint longitudinal and time-to-event (see \code{\link{stan_jm}}) models -#' the estimates are presented separately for each of the distinct submodels. +#' the estimates are presented separately for each of the distinct submodels. #' } #' } -#' +#' #' @seealso \code{\link{summary.stanreg}}, \code{\link{stanreg-methods}} #' print.stanreg <- function(x, digits = 1, detail = TRUE, ...) { @@ -86,15 +86,15 @@ print.stanreg <- function(x, digits = 1, detail = TRUE, ...) { ord <- is_polr(x) && !("(Intercept)" %in% rownames(x$stan_summary)) aux_nms <- .aux_name(x) - + if (!used.optimizing(x)) { - + if (isTRUE(x$stan_function %in% c("stan_lm", "stan_aov"))) { aux_nms <- c("R2", "log-fit_ratio", aux_nms) } mat <- as.matrix(x$stanfit) # don't used as.matrix.stanreg method b/c want access to mean_PPD nms <- setdiff(rownames(x$stan_summary), c("log-posterior", aux_nms)) - + if (gamm) { smooth_sd_nms <- grep("^smooth_sd\\[", nms, value = TRUE) nms <- setdiff(nms, smooth_sd_nms) @@ -110,17 +110,18 @@ print.stanreg <- function(x, digits = 1, detail = TRUE, ...) { cut_mat <- mat[, cut_nms, drop = FALSE] cut_estimates <- .median_and_madsd(cut_mat) } - + if ("car" %in% class(x)) + nms <- setdiff(rownames(x$stan_summary), c("log-posterior", paste0("psi[", 1:length(x$y), "]"))) ppd_nms <- grep("^mean_PPD", nms, value = TRUE) nms <- setdiff(nms, ppd_nms) coef_mat <- mat[, nms, drop = FALSE] estimates <- .median_and_madsd(coef_mat) - + if (mer) { estimates <- estimates[!grepl("^Sigma\\[", rownames(estimates)),, drop=FALSE] } .printfr(estimates, digits, ...) - + if (length(aux_nms)) { aux_estimates <- .median_and_madsd(mat[, aux_nms, drop=FALSE]) cat("\nAuxiliary parameter(s):\n") @@ -137,27 +138,25 @@ print.stanreg <- function(x, digits = 1, detail = TRUE, ...) { if (mer) { cat("\nError terms:\n") print(VarCorr(x), digits = digits + 1, ...) - cat("Num. levels:", + cat("Num. levels:", paste(names(ngrps(x)), unname(ngrps(x)), collapse = ", "), "\n") } if (is(x, "aov")) { print_anova_table(x, digits, ...) - } - + } } else { # used optimization nms <- names(x$coefficients) ppd_nms <- grep("^mean_PPD", rownames(x$stan_summary), value = TRUE) - + estimates <- x$stan_summary[nms, 1:2, drop=FALSE] .printfr(estimates, digits, ...) - + if (length(aux_nms)) { cat("\nAuxiliary parameter(s):\n") .printfr(x$stan_summary[aux_nms, 1:2, drop=FALSE], digits, ...) } - } - + } if (detail) { cat("\n------\n") cat("* For help interpreting the printed output see ?print.stanreg\n") @@ -180,16 +179,16 @@ print.stanmvreg <- function(x, digits = 3, ...) { for (m in 1:M) { cat("\n formula", stubs[m], formula_string(formula(x, m = m))) cat("\n family ", stubs[m], family_plus_link(x, m = m)) - } + } } if (surv) { cat("\n formula (Event):", formula_string(formula(x, m = "Event"))) - cat("\n baseline hazard:", x$basehaz$type_name) + cat("\n baseline hazard:", x$basehaz$type_name) } if (jm) { sel <- grep("^which", rownames(x$assoc), invert = TRUE, value = TRUE) assoc <- lapply(1:M, function(m) { - vals <- sel[which(x$assoc[sel,m] == TRUE)] + vals <- sel[which(x$assoc[sel,m] == TRUE)] paste0(vals, " (Long", m, ")") }) cat("\n assoc: ", paste(unlist(assoc), collapse = ", ")) @@ -197,9 +196,9 @@ print.stanmvreg <- function(x, digits = 3, ...) { cat("\n------\n") mat <- as.matrix(x$stanfit) - nms <- collect_nms(rownames(x$stan_summary), M, + nms <- collect_nms(rownames(x$stan_summary), M, stub = get_stub(x), value = TRUE) - + # Estimates table for longitudinal submodel(s) if (mvmer) { link <- sapply(1:M, function(m) x$family[[m]]$link) @@ -208,54 +207,54 @@ print.stanmvreg <- function(x, digits = 3, ...) { sel <- attr(terms_m, "response") yvar <- rownames(attr(terms_m, "factors"))[sel] if (is.jm(x)) { - cat(paste0("\nLongitudinal submodel", if (M > 1) paste0(" ", m), + cat(paste0("\nLongitudinal submodel", if (M > 1) paste0(" ", m), ": ", yvar,"\n")) } else { cat(paste0("\nSubmodel for y", m, ": ", yvar,"\n")) } coef_mat <- mat[, c(nms$y[[m]], nms$y_extra[[m]]), drop = FALSE] - + # Calculate median and MAD estimates <- .median_and_madsd(coef_mat) - + # Add column with eform - if (link[m] %in% c("log", "logit")) - estimates <- cbind(estimates, - "exp(Median)" = c(exp(estimates[nms$y[[m]], "Median"]), + if (link[m] %in% c("log", "logit")) + estimates <- cbind(estimates, + "exp(Median)" = c(exp(estimates[nms$y[[m]], "Median"]), rep(NA, length(nms$y_extra[[m]])))) - + # Print estimates - rownames(estimates) <- - gsub(paste0("^", get_stub(x), m, "\\|"), "", rownames(estimates)) + rownames(estimates) <- + gsub(paste0("^", get_stub(x), m, "\\|"), "", rownames(estimates)) .printfr(estimates, digits, ...) - } + } } - + # Estimates table for event submodel if (surv) { - cat("\nEvent submodel:\n") + cat("\nEvent submodel:\n") coef_mat <- mat[, c(nms$e, nms$a, nms$e_extra), drop = FALSE] - + # Calculate median and MAD estimates <- .median_and_madsd(coef_mat) - + # Add column with eform - estimates <- cbind(estimates, - "exp(Median)" = c(exp(estimates[c(nms$e, nms$a), "Median"]), + estimates <- cbind(estimates, + "exp(Median)" = c(exp(estimates[c(nms$e, nms$a), "Median"]), rep(NA, length(nms$e_extra)))) - - rownames(estimates) <- gsub("^Event\\|", "", rownames(estimates)) - rownames(estimates) <- gsub("^Assoc\\|", "", rownames(estimates)) + + rownames(estimates) <- gsub("^Event\\|", "", rownames(estimates)) + rownames(estimates) <- gsub("^Assoc\\|", "", rownames(estimates)) .printfr(estimates, digits, ...) } - + # Estimates table for group-level random effects if (mvmer) { - cat("\nGroup-level error terms:\n") + cat("\nGroup-level error terms:\n") print(VarCorr(x), digits = digits + 1, ...) - cat("Num. levels:", paste(names(ngrps(x)), unname(ngrps(x)), - collapse = ", "), "\n") - + cat("Num. levels:", paste(names(ngrps(x)), unname(ngrps(x)), + collapse = ", "), "\n") + # Sample average of the PPD ppd_mat <- mat[, nms$ppd, drop = FALSE] ppd_estimates <- .median_and_madsd(ppd_mat) @@ -263,71 +262,71 @@ print.stanmvreg <- function(x, digits = 3, ...) { if (is.jm(x)) "longitudinal outcomes:\n" else "y:\n") .printfr(ppd_estimates, digits, ...) } - + cat("\n------\n") cat("For info on the priors used see help('prior_summary.stanreg').") - + invisible(x) } #' Summary method for stanreg objects -#' -#' Summaries of parameter estimates and MCMC convergence diagnostics +#' +#' Summaries of parameter estimates and MCMC convergence diagnostics #' (Monte Carlo error, effective sample size, Rhat). -#' +#' #' @export #' @method summary stanreg -#' +#' #' @templateVar stanregArg object #' @template args-stanreg-object #' @template args-regex-pars -#' +#' #' @param ... Currently ignored. #' @param pars An optional character vector specifying a subset of parameters to -#' display. Parameters can be specified by name or several shortcuts can be -#' used. Using \code{pars="beta"} will restrict the displayed parameters to -#' only the regression coefficients (without the intercept). \code{"alpha"} -#' can also be used as a shortcut for \code{"(Intercept)"}. If the model has -#' varying intercepts and/or slopes they can be selected using \code{pars = +#' display. Parameters can be specified by name or several shortcuts can be +#' used. Using \code{pars="beta"} will restrict the displayed parameters to +#' only the regression coefficients (without the intercept). \code{"alpha"} +#' can also be used as a shortcut for \code{"(Intercept)"}. If the model has +#' varying intercepts and/or slopes they can be selected using \code{pars = #' "varying"}. -#' -#' In addition, for \code{stanmvreg} objects there are some additional shortcuts -#' available. Using \code{pars = "long"} will display the +#' +#' In addition, for \code{stanmvreg} objects there are some additional shortcuts +#' available. Using \code{pars = "long"} will display the #' parameter estimates for the longitudinal submodels only (excluding group-specific #' pparameters, but including auxiliary parameters). -#' Using \code{pars = "event"} will display the +#' Using \code{pars = "event"} will display the #' parameter estimates for the event submodel only, including any association -#' parameters. -#' Using \code{pars = "assoc"} will display only the -#' association parameters. +#' parameters. +#' Using \code{pars = "assoc"} will display only the +#' association parameters. #' Using \code{pars = "fixef"} will display all fixed effects, but not -#' the random effects or the auxiliary parameters. -#' \code{pars} and \code{regex_pars} are set to \code{NULL} then all -#' fixed effect regression coefficients are selected, as well as any -#' auxiliary parameters and the log posterior. -#' +#' the random effects or the auxiliary parameters. +#' \code{pars} and \code{regex_pars} are set to \code{NULL} then all +#' fixed effect regression coefficients are selected, as well as any +#' auxiliary parameters and the log posterior. +#' #' If \code{pars} is \code{NULL} all parameters are selected for a \code{stanreg} -#' object, while for a \code{stanmvreg} object all -#' fixed effect regression coefficients are selected as well as any -#' auxiliary parameters and the log posterior. See +#' object, while for a \code{stanmvreg} object all +#' fixed effect regression coefficients are selected as well as any +#' auxiliary parameters and the log posterior. See #' \strong{Examples}. -#' @param probs For models fit using MCMC or one of the variational algorithms, -#' an optional numeric vector of probabilities passed to +#' @param probs For models fit using MCMC or one of the variational algorithms, +#' an optional numeric vector of probabilities passed to #' \code{\link[stats]{quantile}}. -#' @param digits Number of digits to use for formatting numbers when printing. -#' When calling \code{summary}, the value of digits is stored as the +#' @param digits Number of digits to use for formatting numbers when printing. +#' When calling \code{summary}, the value of digits is stored as the #' \code{"print.digits"} attribute of the returned object. -#' -#' @return The \code{summary} method returns an object of class -#' \code{"summary.stanreg"} (or \code{"summary.stanmvreg"}, inheriting -#' \code{"summary.stanreg"}), which is a matrix of -#' summary statistics and +#' +#' @return The \code{summary} method returns an object of class +#' \code{"summary.stanreg"} (or \code{"summary.stanmvreg"}, inheriting +#' \code{"summary.stanreg"}), which is a matrix of +#' summary statistics and #' diagnostics, with attributes storing information for use by the #' \code{print} method. The \code{print} method for \code{summary.stanreg} or -#' \code{summary.stanmvreg} objects is called for its side effect and just returns -#' its input. The \code{as.data.frame} method for \code{summary.stanreg} -#' objects converts the matrix to a data.frame, preserving row and column +#' \code{summary.stanmvreg} objects is called for its side effect and just returns +#' its input. The \code{as.data.frame} method for \code{summary.stanreg} +#' objects converts the matrix to a data.frame, preserving row and column #' names but dropping the \code{print}-related attributes. #' #' @details @@ -344,22 +343,22 @@ print.stanmvreg <- function(x, digits = 3, ...) { #' #' @seealso \code{\link{prior_summary}} to extract or print a summary of the #' priors used for a particular model. -#' +#' #' @examples #' if (.Platform$OS.type != "windows" || .Platform$r_arch != "i386") { #' if (!exists("example_model")) example(example_model) #' summary(example_model, probs = c(0.1, 0.9)) -#' -#' # These produce the same output for this example, +#' +#' # These produce the same output for this example, #' # but the second method can be used for any model -#' summary(example_model, pars = c("(Intercept)", "size", +#' summary(example_model, pars = c("(Intercept)", "size", #' paste0("period", 2:4))) #' summary(example_model, pars = c("alpha", "beta")) -#' +#' #' # Only show parameters varying by group #' summary(example_model, pars = "varying") #' as.data.frame(summary(example_model, pars = "varying")) -#' } +#' #' @importMethodsFrom rstan summary summary.stanreg <- function(object, pars = NULL, @@ -369,7 +368,6 @@ summary.stanreg <- function(object, digits = 1) { mer <- is.mer(object) pars <- collect_pars(object, pars, regex_pars) - if (!used.optimizing(object)) { args <- list(object = object$stanfit, probs = probs) out <- do.call("summary", args)$summary @@ -381,7 +379,7 @@ summary.stanreg <- function(object, pars <- allow_special_parnames(object, pars) out <- out[rownames(out) %in% pars, , drop = FALSE] } - + out <- out[!grepl(":_NEW_", rownames(out), fixed = TRUE), , drop = FALSE] stats <- colnames(out) if ("n_eff" %in% stats) { @@ -401,22 +399,22 @@ summary.stanreg <- function(object, if (is.null(pars)) { famname <- family(object)$family mark <- names(object$coefficients) - if (is.gaussian(famname)) + if (is.gaussian(famname)) mark <- c(mark, "sigma") - if (is.nb(famname)) + if (is.nb(famname)) mark <- c(mark, "reciprocal_dispersion") } else { mark <- NA - if ("alpha" %in% pars) + if ("alpha" %in% pars) mark <- c(mark, "(Intercept)") - if ("beta" %in% pars) + if ("beta" %in% pars) mark <- c(mark, setdiff(names(object$coefficients), "(Intercept)")) mark <- c(mark, setdiff(pars, c("alpha", "beta"))) mark <- mark[!is.na(mark)] } out <- object$stan_summary[mark, , drop=FALSE] } - + structure( out, call = object$call, @@ -534,14 +532,14 @@ as.data.frame.summary.stanreg <- function(x, ...) { #' @rdname summary.stanreg #' @export #' @method summary stanmvreg -summary.stanmvreg <- function(object, pars = NULL, regex_pars = NULL, +summary.stanmvreg <- function(object, pars = NULL, regex_pars = NULL, probs = NULL, ..., digits = 3) { pars <- collect_pars(object, pars, regex_pars) M <- object$n_markers mvmer <- is.mvmer(object) surv <- is.surv(object) - jm <- is.jm(object) - + jm <- is.jm(object) + if (mvmer) { # Outcome variable for each longitudinal submodel y_vars <- sapply(1:M, function(m, object) { @@ -549,41 +547,41 @@ summary.stanmvreg <- function(object, pars = NULL, regex_pars = NULL, sel <- attr(terms_m, "response") ret <- rownames(attr(terms_m, "factors"))[sel] }, object = object) - + # Family and link for each longitudinal submodel fam <- lapply(1:M, function(m) family_plus_link(object, m = m)) } if (jm) { # Association structure sel <- grep("^which", rownames(object$assoc), invert = TRUE, value = TRUE) - assoc <- list_nms(lapply(1:M, function(m) - sel[which(object$assoc[sel,m] == TRUE)]), M) + assoc <- list_nms(lapply(1:M, function(m) + sel[which(object$assoc[sel,m] == TRUE)]), M) } - - # Construct summary table + + # Construct summary table args <- list(object = object$stanfit) - if (!is.null(probs)) + if (!is.null(probs)) args$probs <- probs out <- do.call("summary", args)$summary - nms <- collect_nms(rownames(object$stan_summary), M, + nms <- collect_nms(rownames(object$stan_summary), M, stub = get_stub(object), value = TRUE) if (!is.null(pars)) { - pars2 <- NA + pars2 <- NA if ("alpha" %in% pars) pars2 <- c(pars2, nms$alpha) if ("beta" %in% pars) pars2 <- c(pars2, nms$beta) if ("long" %in% pars) pars2 <- c(pars2, unlist(nms$y), unlist(nms$y_extra)) if ("event" %in% pars) pars2 <- c(pars2, nms$e, nms$a, nms$e_extra) - if ("assoc" %in% pars) pars2 <- c(pars2, nms$a) + if ("assoc" %in% pars) pars2 <- c(pars2, nms$a) if ("fixef" %in% pars) pars2 <- c(pars2, unlist(nms$y), nms$e, nms$a) if ("b" %in% pars) pars2 <- c(pars2, nms$b) - pars2 <- c(pars2, setdiff(pars, + pars2 <- c(pars2, setdiff(pars, c("alpha", "beta", "varying", "b", "long", "event", "assoc", "fixef"))) pars <- pars2[!is.na(pars2)] } else { pars <- rownames(object$stan_summary) pars <- setdiff(pars, b_names(pars, value = TRUE)) - if (used.variational(object)) + if (used.variational(object)) pars <- setdiff(pars, "log-posterior") } out <- out[rownames(out) %in% pars, , drop = FALSE] @@ -593,28 +591,28 @@ summary.stanmvreg <- function(object, pars = NULL, regex_pars = NULL, out[, "n_eff"] <- round(out[, "n_eff"]) if ("se_mean" %in% stats) # So people don't confuse se_mean and sd colnames(out)[stats %in% "se_mean"] <- "mcse" - + # Reorder rows of output table - nms_tmp <- rownames(out) - nms_tmp_y <- lapply(1:M, function(m) + nms_tmp <- rownames(out) + nms_tmp_y <- lapply(1:M, function(m) grep(paste0("^", get_stub(object), m, "\\|"), nms_tmp, value = TRUE)) nms_tmp_e <- grep("^Event\\|", nms_tmp, value = TRUE) nms_tmp_a <- grep("^Assoc\\|", nms_tmp, value = TRUE) nms_tmp_b <- b_names(nms_tmp, value = TRUE) nms_tmp_Sigma <- grep("^Sigma", nms_tmp, value = TRUE) nms_tmp_lp <- grep("^log-posterior$", nms_tmp, value = TRUE) - out <- out[c(unlist(nms_tmp_y), nms_tmp_e, nms_tmp_a, nms_tmp_b, + out <- out[c(unlist(nms_tmp_y), nms_tmp_e, nms_tmp_a, nms_tmp_b, nms_tmp_Sigma, nms_tmp_lp), , drop = FALSE] - + # Output object if (mvmer) out <- structure( - out, y_vars = y_vars, family = fam, n_markers = object$n_markers, + out, y_vars = y_vars, family = fam, n_markers = object$n_markers, n_yobs = object$n_yobs, n_grps = object$n_grps) if (surv) out <- structure( out, n_subjects = object$n_subjects, n_events = object$n_events, - basehaz = object$basehaz) + basehaz = object$basehaz) if (jm) out <- structure( out, id_var = object$id_var, time_var = object$time_var, assoc = assoc) @@ -629,7 +627,7 @@ summary.stanmvreg <- function(object, pars = NULL, regex_pars = NULL, #' @rdname summary.stanreg #' @export #' @method print summary.stanmvreg -print.summary.stanmvreg <- function(x, digits = max(1, attr(x, "print.digits")), +print.summary.stanmvreg <- function(x, digits = max(1, attr(x, "print.digits")), ...) { atts <- attributes(x) mvmer <- atts$stan_function %in% c("stan_mvmer", "stan_jm") @@ -639,11 +637,11 @@ print.summary.stanmvreg <- function(x, digits = max(1, attr(x, "print.digits")), cat("\n function: ", tab, atts$stan_function) if (mvmer) { M <- atts$n_markers - stubs <- paste0("(", if (jm) "Long" else "y", 1:M, "):") + stubs <- paste0("(", if (jm) "Long" else "y", 1:M, "):") for (m in 1:M) { cat("\n formula", stubs[m], formula_string(atts$formula[[m]])) cat("\n family ", stubs[m], atts$family[[m]]) - } + } } if (jm) { cat("\n formula (Event):", formula_string(atts$formula[["Event"]])) @@ -662,18 +660,18 @@ print.summary.stanmvreg <- function(x, digits = max(1, attr(x, "print.digits")), } if (jm) { cat("\n num subjects: ", atts$n_subjects) - cat(paste0("\n num events: ", atts$n_events, " (", + cat(paste0("\n num events: ", atts$n_events, " (", round(100 * atts$n_events/atts$n_subjects, 1), "%)")) - } + } if (!is.null(atts$n_grps)) - cat("\n groups: ", tab, - paste0(names(atts$n_grps), " (", unname(atts$n_grps), ")", collapse = ", ")) + cat("\n groups: ", tab, + paste0(names(atts$n_grps), " (", unname(atts$n_grps), ")", collapse = ", ")) if (atts$algorithm == "sampling") { maxtime <- max(atts$runtime[, "total"]) if (maxtime == 0) maxtime <- "<0.1" cat("\n runtime: ", tab, maxtime, "mins") - } - + } + cat("\n\nEstimates:\n") sel <- which(colnames(x) %in% c("mcse", "n_eff", "Rhat")) if (!length(sel)) { @@ -683,13 +681,13 @@ print.summary.stanmvreg <- function(x, digits = max(1, attr(x, "print.digits")), colnames(xtemp) <- paste(" ", colnames(xtemp)) .printfr(xtemp, digits) cat("\nDiagnostics:\n") - mcse_rhat <- format(round(x[, c("mcse", "Rhat"), drop = FALSE], digits), + mcse_rhat <- format(round(x[, c("mcse", "Rhat"), drop = FALSE], digits), nsmall = digits) n_eff <- format(x[, "n_eff", drop = FALSE], drop0trailing = TRUE) print(cbind(mcse_rhat, n_eff), quote = FALSE) - cat("\nFor each parameter, mcse is Monte Carlo standard error, ", - "n_eff is a crude measure of effective sample size, ", - "and Rhat is the potential scale reduction factor on split chains", + cat("\nFor each parameter, mcse is Monte Carlo standard error, ", + "n_eff is a crude measure of effective sample size, ", + "and Rhat is the potential scale reduction factor on split chains", " (at convergence Rhat=1).\n", sep = '') } invisible(x) @@ -737,7 +735,7 @@ allow_special_parnames <- function(object, pars) { pars2[!is.na(pars2)] } -# Family name with link in parenthesis +# Family name with link in parenthesis # @param x stanreg object # @param ... Optionally include m to specify which submodel for stanmvreg models family_plus_link <- function(x, ...) { @@ -784,7 +782,7 @@ formula_string <- function(formula, break_and_indent = TRUE) { print_anova_table <- function(x, digits, ...) { labels <- attributes(x$terms)$term.labels patterns <- gsub(":", ".*:", labels) - dnms <- dimnames(extract(x$stanfit, pars = "beta", + dnms <- dimnames(extract(x$stanfit, pars = "beta", permuted = FALSE))$parameters groups <- sapply(patterns, simplify = FALSE, FUN = grep, x = dnms) names(groups) <- gsub(".*", "", names(groups), fixed = TRUE) diff --git a/R/prior_summary.R b/R/prior_summary.R index e13a8ad9f..2a3ea0d11 100644 --- a/R/prior_summary.R +++ b/R/prior_summary.R @@ -1,5 +1,5 @@ #' Summarize the priors used for an rstanarm model -#' +#' #' The \code{prior_summary} method provides a summary of the prior distributions #' used for the parameters in a given model. In some cases the user-specified #' prior does not correspond exactly to the prior used internally by @@ -9,48 +9,48 @@ #' or alternatively by fitting the model with the \code{prior_PD} argument set #' to \code{TRUE} (to draw from the prior predictive distribution instead of #' conditioning on the outcome) and then plotting the parameters. -#' +#' #' @aliases prior_summary #' @export #' @templateVar stanregArg object #' @template args-stanreg-object #' @param digits Number of digits to use for rounding. #' @param ... Currently ignored by the method for stanreg objects. -#' -#' @section Intercept (after predictors centered): -#' For \pkg{rstanarm} modeling functions that accept a \code{prior_intercept} -#' argument, the specified prior for the intercept term applies to the -#' intercept after \pkg{rstanarm} internally centers the predictors so they -#' each have mean zero. The estimate of the intercept returned to the user -#' correspond to the intercept with the predictors as specified by the user -#' (unmodified by \pkg{rstanarm}), but when \emph{specifying} the prior the +#' +#' @section Intercept (after predictors centered): +#' For \pkg{rstanarm} modeling functions that accept a \code{prior_intercept} +#' argument, the specified prior for the intercept term applies to the +#' intercept after \pkg{rstanarm} internally centers the predictors so they +#' each have mean zero. The estimate of the intercept returned to the user +#' correspond to the intercept with the predictors as specified by the user +#' (unmodified by \pkg{rstanarm}), but when \emph{specifying} the prior the #' intercept can be thought of as the expected outcome when the predictors are -#' set to their means. The only exception to this is for models fit with the +#' set to their means. The only exception to this is for models fit with the #' \code{sparse} argument set to \code{TRUE} (which is only possible with a #' subset of the modeling functions and never the default). -#' +#' #' @section Adjusted scales: For some models you may see "\code{adjusted scale}" -#' in the printed output and adjusted scales included in the object returned -#' by \code{prior_summary}. These adjusted scale values are the prior scales -#' actually used by \pkg{rstanarm} and are computed by adjusting the prior -#' scales specified by the user to account for the scales of the predictors +#' in the printed output and adjusted scales included in the object returned +#' by \code{prior_summary}. These adjusted scale values are the prior scales +#' actually used by \pkg{rstanarm} and are computed by adjusting the prior +#' scales specified by the user to account for the scales of the predictors #' (as described in the documentation for the \code{\link[=priors]{autoscale}} -#' argument). To disable internal prior scale adjustments set the +#' argument). To disable internal prior scale adjustments set the #' \code{autoscale} argument to \code{FALSE} when setting a prior using one of #' the distributions that accepts an \code{autoscale} argument. For example, #' \code{normal(0, 5, autoscale=FALSE)} instead of just \code{normal(0, 5)}. -#' +#' #' @section Coefficients in Q-space: -#' For the models fit with an \pkg{rstanarm} modeling function that supports -#' the \code{QR} argument (see e.g, \code{\link{stan_glm}}), if \code{QR} is +#' For the models fit with an \pkg{rstanarm} modeling function that supports +#' the \code{QR} argument (see e.g, \code{\link{stan_glm}}), if \code{QR} is #' set to \code{TRUE} then the prior distributions for the regression #' coefficients specified using the \code{prior} argument are not relative to #' the original predictor variables \eqn{X} but rather to the variables in the -#' matrix \eqn{Q} obtained from the \eqn{QR} decomposition of \eqn{X}. -#' +#' matrix \eqn{Q} obtained from the \eqn{QR} decomposition of \eqn{X}. +#' #' In particular, if \code{prior = normal(location,scale)}, then this prior on -#' the coefficients in \eqn{Q}-space can be easily translated into a joint -#' multivariate normal (MVN) prior on the coefficients on the original +#' the coefficients in \eqn{Q}-space can be easily translated into a joint +#' multivariate normal (MVN) prior on the coefficients on the original #' predictors in \eqn{X}. Letting \eqn{\theta} denote the coefficients on #' \eqn{Q} and \eqn{\beta} the coefficients on \eqn{X} then if \eqn{\theta #' \sim N(\mu, \sigma)}{\theta ~ N(\mu, \sigma)} the corresponding prior on @@ -63,42 +63,42 @@ #' default), in which case the matrices actually used are #' \eqn{Q^\ast = Q \sqrt{n-1}}{Q* = Q (n-1)^0.5} and \eqn{R^\ast = #' \frac{1}{\sqrt{n-1}} R}{R* = (n-1)^(-0.5) R}. If \code{autoscale = FALSE} -#' we instead scale such that the lower-right element of \eqn{R^\ast}{R*} is -#' \eqn{1}, which is useful if you want to specify a prior on the coefficient -#' of the last predictor in its original units (see the documentation for the +#' we instead scale such that the lower-right element of \eqn{R^\ast}{R*} is +#' \eqn{1}, which is useful if you want to specify a prior on the coefficient +#' of the last predictor in its original units (see the documentation for the #' \code{\link[=stan_glm]{QR}} argument). -#' +#' #' If you are interested in the prior on \eqn{\beta} implied by the prior on #' \eqn{\theta}, we strongly recommend visualizing it as described above in #' the \strong{Description} section, which is simpler than working it out #' analytically. -#' +#' #' @return A list of class "prior_summary.stanreg", which has its own print #' method. -#' +#' #' @seealso The \link[=priors]{priors help page} and the \emph{Prior #' Distributions} vignette. -#' +#' #' @examples #' if (.Platform$OS.type != "windows" || .Platform$r_arch != "i386") { #' if (!exists("example_model")) example(example_model) #' prior_summary(example_model) -#' +#' #' priors <- prior_summary(example_model) #' names(priors) #' priors$prior$scale #' priors$prior$adjusted_scale -#' -#' # for a glm with adjusted scales (see Details, above), compare -#' # the default (rstanarm adjusting the scales) to setting +#' +#' # for a glm with adjusted scales (see Details, above), compare +#' # the default (rstanarm adjusting the scales) to setting #' # autoscale=FALSE for prior on coefficients -#' fit <- stan_glm(mpg ~ wt + am, data = mtcars, -#' prior = normal(0, c(2.5, 4)), -#' prior_intercept = normal(0, 5), -#' iter = 10, chains = 1) # only for demonstration +#' fit <- stan_glm(mpg ~ wt + am, data = mtcars, +#' prior = normal(0, c(2.5, 4)), +#' prior_intercept = normal(0, 5), +#' iter = 10, chains = 1) # only for demonstration #' prior_summary(fit) -#' -#' fit2 <- update(fit, prior = normal(0, c(2.5, 4), autoscale=FALSE), +#' +#' fit2 <- update(fit, prior = normal(0, c(2.5, 4), autoscale=FALSE), #' prior_intercept = normal(0, 5, autoscale=FALSE)) #' prior_summary(fit2) #' } @@ -107,18 +107,18 @@ prior_summary.stanreg <- function(object, digits = 2,...) { if (is.null(x)) { message("Priors not found in stanreg object.") return(invisible(NULL)) - } + } if (is.stanmvreg(object)) { M <- get_M(object) - x <- structure(x, M = M) + x <- structure(x, M = M) } structure(x, class = "prior_summary.stanreg", QR = used.QR(object), sparse = used.sparse(object), - model_name = deparse(substitute(object)), + model_name = deparse(substitute(object)), stan_function = object$stan_function, print_digits = digits) - + } #' @export @@ -134,21 +134,21 @@ print.prior_summary.stanreg <- function(x, digits, ...) { sparse <- attr(x, "sparse") model_name <- attr(x, "model_name") stan_function <- attr(x, "stan_function") - + msg <- paste0("Priors for model '", model_name, "'") cat(msg, "\n------") - + if (!stan_function == "stan_mvmer") { if (!is.null(x[["prior_intercept"]])) .print_scalar_prior( - x[["prior_intercept"]], - txt = paste0("Intercept", if (!sparse) " (after predictors centered)"), + x[["prior_intercept"]], + txt = paste0("Intercept", if (!sparse) " (after predictors centered)"), formatters ) if (!is.null(x[["prior"]])) .print_vector_prior( - x[["prior"]], - txt = paste0("\nCoefficients", if (QR) " (in Q-space)"), + x[["prior"]], + txt = paste0("\nCoefficients", if (QR) " (in Q-space)"), formatters = formatters ) if (!is.null(x[["prior_aux"]])) { @@ -157,25 +157,25 @@ print.prior_summary.stanreg <- function(x, digits, ...) { if (aux_dist %in% c("normal", "student_t", "cauchy")) x[["prior_aux"]][["dist"]] <- paste0("half-", aux_dist) .print_scalar_prior( - x[["prior_aux"]], - txt = paste0("\nAuxiliary (", aux_name, ")"), + x[["prior_aux"]], + txt = paste0("\nAuxiliary (", aux_name, ")"), formatters ) - } + } } else { # unique to stan_mvmer M <- attr(x, "M") for (m in 1:M) { if (!is.null(x[["prior_intercept"]][[m]])) .print_scalar_prior( - x[["prior_intercept"]][[m]], - txt = paste0(if (m > 1) "\n", "y", m, "|Intercept", if (!sparse) - " (after predictors centered)"), + x[["prior_intercept"]][[m]], + txt = paste0(if (m > 1) "\n", "y", m, "|Intercept", if (!sparse) + " (after predictors centered)"), formatters ) if (!is.null(x[["prior"]][[m]])) .print_vector_prior( - x[["prior"]][[m]], - txt = paste0("\ny", m, "|Coefficients", if (QR) " (in Q-space)"), + x[["prior"]][[m]], + txt = paste0("\ny", m, "|Coefficients", if (QR) " (in Q-space)"), formatters = formatters ) if (!is.null(x[["prior_aux"]][[m]])) { @@ -184,39 +184,39 @@ print.prior_summary.stanreg <- function(x, digits, ...) { if (aux_dist %in% c("normal", "student_t", "cauchy")) x[["prior_aux"]][[m]][["dist"]] <- paste0("half-", aux_dist) .print_scalar_prior( - x[["prior_aux"]][[m]], - txt = paste0("\ny", m, "|Auxiliary (", aux_name, ")"), + x[["prior_aux"]][[m]], + txt = paste0("\ny", m, "|Auxiliary (", aux_name, ")"), formatters ) } - } + } } - + # unique to stan_betareg if (!is.null(x[["prior_intercept_z"]])) .print_scalar_prior( - x[["prior_intercept_z"]], - txt = paste0("\nIntercept_z", if (!sparse) " (after predictors centered)"), + x[["prior_intercept_z"]], + txt = paste0("\nIntercept_z", if (!sparse) " (after predictors centered)"), formatters ) if (!is.null(x[["prior_z"]])) .print_vector_prior(x[["prior_z"]], txt = "\nCoefficients_z", formatters) - + # unique to stan_jm if (stan_function == "stan_jm") { M <- attr(x, "M") for (m in 1:M) { if (!is.null(x[["priorLong_intercept"]][[m]])) .print_scalar_prior( - x[["priorLong_intercept"]][[m]], - txt = paste0(if (m > 1) "\n", "Long", m, "|Intercept", if (!sparse) - " (after predictors centered)"), + x[["priorLong_intercept"]][[m]], + txt = paste0(if (m > 1) "\n", "Long", m, "|Intercept", if (!sparse) + " (after predictors centered)"), formatters ) if (!is.null(x[["priorLong"]][[m]])) .print_vector_prior( - x[["priorLong"]][[m]], - txt = paste0("\nLong", m, "|Coefficients", if (QR) " (in Q-space)"), + x[["priorLong"]][[m]], + txt = paste0("\nLong", m, "|Coefficients", if (QR) " (in Q-space)"), formatters = formatters ) if (!is.null(x[["priorLong_aux"]][[m]])) { @@ -225,22 +225,22 @@ print.prior_summary.stanreg <- function(x, digits, ...) { if (aux_dist %in% c("normal", "student_t", "cauchy")) x[["priorLong_aux"]][[m]][["dist"]] <- paste0("half-", aux_dist) .print_scalar_prior( - x[["priorLong_aux"]][[m]], - txt = paste0("\nLong", m, "|Auxiliary (", aux_name, ")"), + x[["priorLong_aux"]][[m]], + txt = paste0("\nLong", m, "|Auxiliary (", aux_name, ")"), formatters ) } } if (!is.null(x[["priorEvent_intercept"]])) .print_scalar_prior( - x[["priorEvent_intercept"]], - txt = paste0("\nEvent|Intercept", if (!sparse) " (after predictors centered)"), + x[["priorEvent_intercept"]], + txt = paste0("\nEvent|Intercept", if (!sparse) " (after predictors centered)"), formatters ) if (!is.null(x[["priorEvent"]])) .print_vector_prior( - x[["priorEvent"]], - txt = "\nEvent|Coefficients", + x[["priorEvent"]], + txt = "\nEvent|Coefficients", formatters = formatters ) if (!is.null(x[["priorEvent_aux"]])) { @@ -250,30 +250,30 @@ print.prior_summary.stanreg <- function(x, digits, ...) { (aux_dist %in% c("normal", "student_t", "cauchy"))) { # weibull x[["priorEvent_aux"]][["dist"]] <- paste0("half-", aux_dist) .print_scalar_prior( - x[["priorEvent_aux"]], - txt = paste0("\nEvent|Auxiliary (", aux_name, ")"), + x[["priorEvent_aux"]], + txt = paste0("\nEvent|Auxiliary (", aux_name, ")"), formatters - ) + ) } else { # bs or piecewise .print_vector_prior( - x[["priorEvent_aux"]], - txt = paste0("\nEvent|Auxiliary (", aux_name, ")"), + x[["priorEvent_aux"]], + txt = paste0("\nEvent|Auxiliary (", aux_name, ")"), formatters ) } } if (!is.null(x[["priorEvent_assoc"]])) .print_vector_prior( - x[["priorEvent_assoc"]], - txt = "\nAssociation parameters", + x[["priorEvent_assoc"]], + txt = "\nAssociation parameters", formatters = formatters ) - } - + } + # unique to stan_(g)lmer, stan_gamm4, stan_mvmer, or stan_jm if (!is.null(x[["prior_covariance"]])) .print_covariance_prior(x[["prior_covariance"]], txt = "\nCovariance", formatters) - + # unique to stan_polr if (!is.null(x[["prior_counts"]])) { p <- x[["prior_counts"]] @@ -284,10 +284,30 @@ print.prior_summary.stanreg <- function(x, digits, ...) { if (!is.null(x[["scobit_exponent"]])) { p <- x[["scobit_exponent"]] cat("\n\nScobit Exponent\n ~", - paste0(p$dist, "(shape = ", .fr2(p$shape), + paste0(p$dist, "(shape = ", .fr2(p$shape), ", rate = ", .fr2(p$rate), ")")) } + # unique to stan_besag/stan_bym + if (!is.null(x[["prior_sigma"]])) + .print_scalar_prior( + x[["prior_sigma"]], + txt = paste0("\nSigma", if (!sparse) " (after predictors centered)"), + formatters + ) + if (!is.null(x[["prior_rho"]])) + .print_scalar_prior( + x[["prior_rho"]], + txt = paste0("\nRho", if (!sparse) " (after predictors centered)"), + formatters + ) + if (!is.null(x[["prior_tau"]])) + .print_scalar_prior( + x[["prior_tau"]], + txt = paste0("\nTau", if (!sparse) " (after predictors centered)"), + formatters + ) + cat("\n------\n") cat("See help('prior_summary.stanreg') for more details\n") invisible(x) @@ -306,7 +326,7 @@ used.sparse <- function(x) { isTRUE(getCall(x)[["sparse"]]) } -# +# # @param x numeric vector # @param formatter a formatting function to apply (see .fr2, .fr3 above) # @param N the maximum number of values to include before replacing the rest @@ -316,9 +336,9 @@ used.sparse <- function(x) { if (K < 2) return(x) paste0( - "[", - paste(c(formatter(x[1:min(N, K)]), if (N < K) "..."), - collapse = ","), + "[", + paste(c(formatter(x[1:min(N, K)]), if (N < K) "..."), + collapse = ","), "]" ) } @@ -331,13 +351,10 @@ used.sparse <- function(x) { # in prior_summary.stanreg). The first is used for format all numbers except # for adjusted scales, for which the second function is used. This is kind of # hacky and should be replaced at some point. -# - .print_scalar_prior <- function(p, txt = "Intercept", formatters = list()) { stopifnot(length(formatters) == 2) .f1 <- formatters[[1]] - .f2 <- formatters[[2]] - + .f2 <- formatters[[2]] .cat_scalar_prior <- function(p, adjusted = FALSE, prepend_chars = "\n ~") { if (adjusted) { p$scale <- p$adjusted_scale @@ -369,7 +386,6 @@ used.sparse <- function(x) { cat("\n Adjusted prior:") .cat_scalar_prior(p, adjusted = TRUE, prepend_chars = "\n ~") } - } .print_covariance_prior <- function(p, txt = "Covariance", formatters = list()) { @@ -380,12 +396,12 @@ used.sparse <- function(x) { p$shape <- .format_pars(p$shape, .f1) p$scale <- .format_pars(p$scale, .f1) cat(paste0("\n", txt, "\n ~"), - paste0(p$dist, "(", + paste0(p$dist, "(", "reg. = ", .f1(p$regularization), - ", conc. = ", .f1(p$concentration), + ", conc. = ", .f1(p$concentration), ", shape = ", .f1(p$shape), ", scale = ", .f1(p$scale), ")") - ) + ) } else if (p$dist == "lkj") { .f1 <- formatters[[1]] .f2 <- formatters[[2]] @@ -395,11 +411,11 @@ used.sparse <- function(x) { if (!is.null(p$adjusted_scale)) p$adjusted_scale <- .format_pars(p$adjusted_scale, .f2) cat(paste0("\n", txt, "\n ~"), - paste0(p$dist, "(", + paste0(p$dist, "(", "reg. = ", .f1(p$regularization), - ", df = ", .f1(p$df), + ", df = ", .f1(p$df), ", scale = ", .f1(p$scale), ")") - ) + ) if (!is.null(p$adjusted_scale)) cat("\n **adjusted scale =", .f2(p$adjusted_scale)) } diff --git a/R/priors.R b/R/priors.R index 241661517..d2ac45038 100644 --- a/R/priors.R +++ b/R/priors.R @@ -565,6 +565,12 @@ R2 <- function(location = NULL, what = c("mode", "mean", "median", "log")) { #' @rdname priors #' @export +beta <- function(alpha = 2, beta = 2) { + validate_parameter_value(alpha) + validate_parameter_value(beta) + nlist(dist = "beta", alpha, beta) +} + #' @param family Not currently used. default_prior_intercept = function(family) { # family arg not used, but we can use in the future to do different things diff --git a/R/stan_besag.R b/R/stan_besag.R new file mode 100644 index 000000000..1219ba337 --- /dev/null +++ b/R/stan_besag.R @@ -0,0 +1,157 @@ +# Part of the rstanarm package for estimating model parameters +# Copyright (C) 2013, 2014, 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. + +#' Bayesian spatial ICAR models via Stan +#' +#' Spatial regression modeling with an intrinsic conditional autoregressive (ICAR) prior. +#' +#' @rdname stan_besag +#' @export +#' +#' @templateVar fun stan_besag +#' @templateVar fitfun stan_spatial.fit +#' @templateVar pkg rstanarm +#' @templateVar pkgfun stan_glm +#' @templateVar sameargs family +#' @template return-stanreg-object +#' @template return-stanfit-object +#' @template args-formula-data-subset +#' @template args-same-as +#' @template args-x-y +#' @template args-dots +#' @template args-prior_intercept +#' @template args-priors +#' @template args-prior_PD +#' @template args-algorithm +#' @template args-adapt_delta +#' @template args-QR +#' +#' @param trials If \code{family = binomial()} then a vector of trials (equal in +#' length to the outcome) must be declared. +#' @param W An N-by-N spatial weight matrix. +#' @param prior_structured The prior distribution on the variance of the non-centered +#' structured (spatial) effect. +#' @param order Order of the spatial random walk. Specifying \code{order = 2} +#' will smooth the spatial variation. The default is \code{order = 1}. +#' +#' @details The \code{stan_besag} model is similar to the Besag model in R-INLA. +#' However, instead of using the integrated nested Laplace approximation +#' (INLA) method, full Bayesian estimation is performed (if \code{algorithm} +#' is \code{"sampling"}) via MCMC. The model includes priors on the intercept, +#' regression coefficients, overall spatial variation, and any applicable +#' auxiliary parameters. The \code{stan_besag} function calls the workhorse +#' \code{stan_spatial.fit} function, but it is also possible to call the +#' latter directly. +#' +#' @seealso The vignette for \code{stan_besag}. +#' +#' @references Riebler, A., Sorbye, S.H., Simpson, D., Rue, H. (2016). An +#' intuitive Bayesian spatial model for disease mapping that accounts for +#' scaling. arXiv preprint arXiv:1601.01180. +#' +#' Besag, J. (1974). Spatial Interaction and the Statistical Analysis of +#' Lattice Systems. Journal of the Royal Statistical Society. Vol 36, No 2, +#' p192-236. +#' +#' @examples +#' \dontrun{ +#' ### Simulated Data on a Lattice +#' +#' data("lattice", package = "rstanarm") +#' +#' # plot GMRF +#' grid_sim <- grid_sim15 +#' var_range_gmrf <- seq(min(grid_sim@data$gmrf), max(grid_sim@data$gmrf), length = 50) +#' spplot(grid_sim, "gmrf", at = var_range_gmrf, main = expression(paste(phi, " (GMRF)")), +#' col.regions = colorRampPalette(c("#ef8a62", "#f7f7f7", "#67a9cf"))(50)) +#' +#' # Convert a spatial polygon to an N-by-N weight matrix +#' sp2weightmatrix <- function(spatialpolygon) { +#' spdep::nb2mat(spdep::poly2nb(spatialpolygon, queen = TRUE), style = "B", zero.policy = TRUE) +#' } +#' +#' # convert spatial object to neighborhood matrix +#' W <- sp2weightmatrix(grid_sim) +#' # W_sparse <- Matrix(W, sparse = TRUE, dimnames = list(NULL,NULL)) +#' +#' # simulate predictor/outcome +#' x <- rnorm(nrow(W), 3, 1) +#' phi <- grid_sim@data$gmrf +#' tau <- 2 +#' y <- rnorm(nrow(W), 0.5 + 3*x + phi*tau, 1) +#' +#' # fit the model +#' fit_besag <- stan_besag(y ~ 1 + x, data = data.frame(y=y,x=x), family = gaussian(), +#' W = W, iter = 500, chains = 4) +#' fit_besag +#' pp_check(fit_besag) +#'} +#' + +stan_besag <- function(formula, + family = NULL, + data, + trials = NULL, + W, + order = 1, + ..., + prior = normal(), prior_intercept = normal(), + prior_structured = normal(), prior_aux = NULL, + prior_PD = FALSE, + algorithm = c("sampling", "meanfield", "fullrank"), + adapt_delta = NULL, + QR = FALSE) { + stan_function <- "stan_besag" + mc <- match.call(expand.dots = FALSE) + algorithm <- match.arg(algorithm) + family <- validate_family(family) + mf <- model.frame(mc, data) + mt <- terms(formula, data = data) + Y <- array1D_check(model.response(mf, type = "any")) + X <- model.matrix(formula, data) + + stanfit <- stan_spatial.fit(x = X, y = Y, w = W, + trials = trials, + family = family, + stan_function = stan_function, + order = order, + ..., + prior = prior, + prior_intercept = prior_intercept, + prior_aux = prior_aux, + prior_tau = prior_structured, + prior_PD = prior_PD, + algorithm = algorithm, adapt_delta = adapt_delta, + QR = QR) + fit <- nlist(stanfit, + algorithm, + data, + x = X, y = Y, + family, + formula, + model = mf, + terms = mt, + call = match.call(), + stan_function = stan_function) + + if (family$family == "binomial") { + # fit$family <- binomial(link = "logit") + fit$trials <- trials + } + out <- stanreg(fit) + structure(out, class = c("stanreg", "car")) +} diff --git a/R/stan_bym.R b/R/stan_bym.R new file mode 100644 index 000000000..475de20ed --- /dev/null +++ b/R/stan_bym.R @@ -0,0 +1,161 @@ +# Part of the rstanarm package for estimating model parameters +# Copyright (C) 2013, 2014, 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. + +#' Bayesian spatial CAR BYM models via Stan +#' +#' Spatial regression modeling with a variant of the Besag, York, Mollie (BYM) +#' conditional autoregressive (CAR) prior. +#' +#' @export +#' +#' @templateVar fun stan_bym +#' @templateVar fitfun stan_spatial.fit +#' @templateVar pkg rstanarm +#' @templateVar pkgfun stan_glm +#' @templateVar sameargs family +#' @template return-stanreg-object +#' @template return-stanfit-object +#' @template args-formula-data-subset +#' @template args-same-as +#' @template args-x-y +#' @template args-dots +#' @template args-prior_intercept +#' @template args-priors +#' @template args-prior_PD +#' @template args-algorithm +#' @template args-adapt_delta +#' @template args-QR +#' +#' @param trials If \code{family = binomial()} then a vector of trials (equal in +#' length to the outcome) must be declared. +#' @param W An N-by-N spatial weight matrix. +#' @param prior_unstructured The prior on the marginal variance contribution of the +#' unstructured (random) effect. +#' @param prior_structured The prior on the marginal variance contribution of the +#' structured (spatial) effect. +#' @param order Order of the spatial random walk. Specifying \code{order = 2} +#' will smooth the spatial variation. The default is \code{order = 1}. +#' +#' @details The \code{stan_bym} model is similar to the BYM model in R-INLA. +#' However, instead of using the integrated nested Laplace approximation +#' (INLA) method, full Bayesian estimation is performed (if \code{algorithm} +#' is \code{"sampling"}) via MCMC. The model includes priors on the intercept, +#' regression coefficients, spatial mixing parameter, overall spatial +#' variation, and any applicable auxiliary parameters. The \code{stan_bym} +#' function calls the workhorse \code{stan_spatial.fit} function, but it is +#' also possible to call the latter directly. +#' +#' @seealso The vignette for \code{stan_bym}. +#' +#' @references Riebler, A., Sorbye, S.H., Simpson, D., Rue, H. (2016). An +#' intuitive Bayesian spatial model for disease mapping that accounts for +#' scaling. arXiv preprint arXiv:1601.01180. +#' +#' Besag, J., York, J. and Mollié, A. (1991). Bayesian image restoration, with +#' two applications in spatial statistics. Annals of the Institute of +#' Statistical Mathematics. Vol. 43, No. 01, p1-20. +#' +#' @examples +#' \dontrun{ +#' ### Simulated Data on a Lattice +#' +#' data("lattice", package = "rstanarm") +#' +#' # plot GMRF +#' grid_sim <- grid_sim15 +#' var_range_gmrf <- seq(min(grid_sim@data$gmrf), max(grid_sim@data$gmrf), length = 50) +#' spplot(grid_sim, "gmrf", at = var_range_gmrf, main = expression(paste(phi, " (GMRF)")), +#' col.regions = colorRampPalette(c("#ef8a62", "#f7f7f7", "#67a9cf"))(50)) +#' +#' # Convert a spatial polygon to an N-by-N weight matrix +#' sp2weightmatrix <- function(spatialpolygon) { +#' spdep::nb2mat(spdep::poly2nb(spatialpolygon, queen = TRUE), style = "B", zero.policy = TRUE) +#' } +#' +#' # convert spatial object to neighborhood matrix +#' W <- sp2weightmatrix(grid_sim) +#' # W_sparse <- Matrix(W, sparse = TRUE, dimnames = list(NULL,NULL)) +#' +#' # simulate predictor/outcome +#' x <- rnorm(nrow(W), 3, 1) +#' phi <- grid_sim@data$gmrf +#' theta <- rnorm(nrow(W), 0, 1) +#' tau <- 1 +#' rho <- 2 +#' spatial_data <- data.frame(x) +#' y <- rnorm(nrow(W), 0.5 + 3*x + phi*rho + theta*tau, 1) +#' +#' # fit the model +#' fit_bym <- stan_bym(y ~ 1 + x, data = data.frame(y=y,x=x), +#' W = W, iter = 500) +#' fit_bym +#' pp_check(fit_bym) +#' } +#' + +stan_bym <- function(formula, + family = gaussian(), + data, + trials = NULL, + W, + order = 1, + ..., + prior = normal(), prior_intercept = normal(), + prior_unstructured = normal(), prior_structured = normal(), prior_aux = NULL, + prior_PD = FALSE, + algorithm = c("sampling", "meanfield", "fullrank"), + adapt_delta = NULL, + QR = FALSE) { + stan_function <- "stan_bym" + mc <- match.call(expand.dots = FALSE) + algorithm <- match.arg(algorithm) + family <- validate_family(family) + mf <- model.frame(mc, data) + mt <- terms(formula, data = data) + Y <- array1D_check(model.response(mf, type = "any")) + X <- model.matrix(formula, data) + + stanfit <- stan_spatial.fit(x = X, y = Y, w = W, + trials = trials, + family = family, + stan_function = stan_function, + order = order, + ..., + prior = prior, + prior_intercept = prior_intercept, + prior_aux = prior_aux, prior_rho = prior_structured, prior_tau = prior_unstructured, + prior_PD = prior_PD, + algorithm = algorithm, adapt_delta = adapt_delta, + QR = QR) + fit <- nlist(stanfit, + algorithm, + data, + x = X, y = Y, + family, + formula, + model = mf, + terms = mt, + call = match.call(), + stan_function = stan_function) + + if (family$family == "binomial") { + fit$family <- binomial(link = "logit") + fit$trials <- trials + } + out <- stanreg(fit) + structure(out, class = c("stanreg", "car")) +} diff --git a/R/stan_bym2.R b/R/stan_bym2.R new file mode 100644 index 000000000..a1d2933ad --- /dev/null +++ b/R/stan_bym2.R @@ -0,0 +1,174 @@ +# Part of the rstanarm package for estimating model parameters +# Copyright (C) 2013, 2014, 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. + +#' Bayesian spatial CAR BYM variant models via Stan +#' +#' Spatial regression modeling with a variant of the Besag, York, Mollie (BYM) +#' conditional autoregressive (CAR) prior that accounts for scaling. +#' +#' @export +#' +#' @templateVar fun stan_bym2 +#' @templateVar fitfun stan_spatial.fit +#' @templateVar pkg rstanarm +#' @templateVar pkgfun stan_glm +#' @templateVar sameargs family +#' @template return-stanreg-object +#' @template return-stanfit-object +#' @template args-formula-data-subset +#' @template args-same-as +#' @template args-x-y +#' @template args-dots +#' @template args-prior_intercept +#' @template args-priors +#' @template args-prior_PD +#' @template args-algorithm +#' @template args-adapt_delta +#' @template args-QR +#' +#' @param trials If \code{family = binomial()} then a vector of trials (equal in +#' length to the outcome) must be declared. +#' @param W An N-by-N spatial weight matrix. +#' @param prior_structured The prior on the marginal variance contribution of the +#' structured (spatial) and unstructured (random) effect. +#' @param prior_mixing The prior on the proportion of the marginal variance that is +#' explained by the structured (spatial) effect. The hyperparameter \code{rho} +#' is on the unit interval so users have the option of declaring a Beta prior +#' distribution or a flat prior. A prior distribution with most of the mass +#' around 1 is analogous to the prior belief that there exists a strong +#' spatial relationship on the graph. +#' @param order Order of the spatial random walk. Specifying \code{order = 2} +#' will smooth the spatial variation. The default is \code{order = 1}. +#' +#' @details The \code{stan_bym2} model is similar to the BYM2 model in R-INLA. +#' However, instead of using the integrated nested Laplace approximation +#' (INLA) method, full Bayesian estimation is performed (if \code{algorithm} +#' is \code{"sampling"}) via MCMC. The model includes priors on the intercept, +#' regression coefficients, spatial mixing parameter, overall spatial +#' variation, and any applicable auxiliary parameters. The \code{stan_bym2} +#' function calls the workhorse \code{stan_spatial.fit} function, but it is +#' also possible to call the latter directly. +#' +#' @seealso The vignette for \code{stan_bym2}. +#' +#' @references Riebler, A., Sorbye, S.H., Simpson, D., Rue, H. (2016). An +#' intuitive Bayesian spatial model for disease mapping that accounts for +#' scaling. arXiv preprint arXiv:1601.01180. +#' +#' Besag, J., York, J. and Mollié, A. (1991). Bayesian image restoration, with +#' two applications in spatial statistics. Annals of the Institute of +#' Statistical Mathematics. Vol. 43, No. 01, p1-20. +#' +#' Simpson, D., Rue, H., Martins, T.G., Riebler, A. and Sorbye, S.H. (2015). +#' Penalising model component complexity: A principled, practical approach to +#' constructing priors. arXiv preprint arXiv:1403.4630. +#' +#' @examples +#' \dontrun{ +#' ### Simulated Data on a Lattice +#' +#' data("lattice", package = "rstanarm") +#' +#' # plot GMRF +#' grid_sim <- grid_sim15 +#' var_range_gmrf <- seq(min(grid_sim@data$gmrf), max(grid_sim@data$gmrf), length = 50) +#' spplot(grid_sim, "gmrf", at = var_range_gmrf, main = expression(paste(phi, " (GMRF)")), +#' col.regions = colorRampPalette(c("#ef8a62", "#f7f7f7", "#67a9cf"))(50)) +#' +#' # Convert a spatial polygon to an N-by-N weight matrix +#' sp2weightmatrix <- function(spatialpolygon) { +#' spdep::nb2mat(spdep::poly2nb(spatialpolygon, queen = TRUE), style = "B", zero.policy = TRUE) +#' } +#' +#' # convert spatial object to neighborhood matrix +#' W <- sp2weightmatrix(grid_sim) +#' # W_sparse <- Matrix(W, sparse = TRUE, dimnames = list(NULL,NULL)) +#' +#' # simulate predictor/outcome +#' x <- rnorm(nrow(W), 3, 1) +#' phi <- grid_sim@data$gmrf +#' theta <- rnorm(nrow(W), 0, 1) +#' tau <- 1 +#' rho <- 0.7 +#' psi <- (1/tau) * (sqrt(1-rho)*theta + sqrt(rho)*phi) +#' y <- rnorm(nrow(W), 0 + 0.4 * x + psi, 1) +#' +#' # fit the model +#' fit_bym2 <- stan_bym2(y ~ 1 + x, data = data.frame(y=y,x=x), +#' W = W, iter = 1e3, chains = 4, cores = 2) +#' fit_bym2 +#' pp_check(fit_bym2) +#' } +#' + +stan_bym2 <- function(formula, + family = gaussian(), + data, + trials = NULL, + W, + order = 1, + ..., + prior = normal(), prior_intercept = normal(), + prior_structured = normal(), prior_mixing = beta(0.5,0.5), prior_aux = NULL, + prior_PD = FALSE, + algorithm = c("sampling", "meanfield", "fullrank"), + adapt_delta = NULL, + QR = FALSE) { + stan_function <- "stan_bym2" + mc <- match.call(expand.dots = FALSE) + algorithm <- match.arg(algorithm) + family <- validate_family(family) + mf <- model.frame(mc, data) + mt <- terms(formula, data = data) + Y <- array1D_check(model.response(mf, type = "any")) + X <- model.matrix(formula, data) + + if (!is.null(prior_mixing)) { + if (prior_mixing$dist != "beta") + stop("'prior_mixing' must be either beta() or NULL.") + } + + stanfit <- stan_spatial.fit(x = X, y = Y, w = W, + trials = trials, + family = family, + stan_function = stan_function, + order = order, + ..., + prior = prior, + prior_intercept = prior_intercept, + prior_aux = prior_aux, prior_rho = prior_mixing, prior_tau = prior_structured, + prior_PD = prior_PD, + algorithm = algorithm, adapt_delta = adapt_delta, + QR = QR) + fit <- nlist(stanfit, + algorithm, + data, + x = X, y = Y, + family, + formula, + model = mf, + terms = mt, + call = match.call(), + stan_function = stan_function) + + if (family$family == "binomial") { + fit$family <- binomial(link = "logit") + fit$trials <- trials + } + out <- stanreg(fit) + structure(out, class = c("stanreg", "car")) +} diff --git a/R/stan_glm.fit.R b/R/stan_glm.fit.R index db79062e1..6fbe834e7 100644 --- a/R/stan_glm.fit.R +++ b/R/stan_glm.fit.R @@ -289,7 +289,7 @@ stan_glm.fit <- if (length(weights) > 0 && all(weights == 1)) weights <- double() if (length(offset) > 0 && all(offset == 0)) offset <- double() - + # create entries in the data block of the .stan file standata <- nlist( N = nrow(xtemp), diff --git a/R/stan_spatial.fit.R b/R/stan_spatial.fit.R new file mode 100644 index 000000000..55d4a2a85 --- /dev/null +++ b/R/stan_spatial.fit.R @@ -0,0 +1,575 @@ +# Part of the rstanarm package for estimating model parameters +# Copyright (C) 2013, 2014, 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. + +#' Workhorse function for CAR models. +#' +#' Both \code{stan_besag} and \code{stan_bym2} call \code{stan_spatial.fit} to +#' fit the appropriate spatial model. See the documentation for either modeling +#' function for further details on the arguments of \code{stan_spatial.fit}. +#' +#' @export +#' + +stan_spatial.fit <- function(x, y, w, + stan_function = c("stan_besag", "stan_bym2"), + family = NULL, + trials = NULL, + order = c(1,2), + ..., + prior = normal(), prior_intercept = normal(), + prior_tau = normal(), prior_aux = NULL, prior_rho = NULL, + prior_PD = FALSE, + algorithm = c("sampling", "meanfield", "fullrank"), + adapt_delta = NULL, + QR = FALSE) { + + w[upper.tri(w)] <- 0 + + # convert W to a sparse matrix if not already sparse. + if(!is(w, "sparseMatrix")) + w <- Matrix(w, sparse = TRUE) + + # pull out adjacency pairs from W + edges <- summary(w) # analagous to `which(w == 1, arr.ind = TRUE)` on dense matrix + edges <- edges[,grep("^i$|^j$", colnames(edges))] + + algorithm <- match.arg(algorithm) + + family <- validate_family(family) + supported_families <- c("binomial", "gaussian", "poisson", "Gamma", "neg_binomial_2") + fam <- which(pmatch(supported_families, family$family, nomatch = 0L) == 1L) + if (!length(fam)) + stop("'family' must be one of ", paste(supported_families, collapse = ", ")) + + supported_links <- supported_glm_links(supported_families[fam]) + link <- which(supported_links == family$link) + if (!length(link)) + stop("'link' must be one of ", paste(supported_links, collapse = ", ")) + + family_num <- switch(family$family, + gaussian = 1, + poisson = 6, + neg_binomial_2 = 7, + binomial = 5, + Gamma = 2) + + # for when consistent-family-numbers gets merged + # family_num <- switch(family$family, + # gaussian = 1, + # Gamma = 2, + # inv_gaussian = 3, + # beta = 4, + # binomial = 5, + # poisson = 6, + # neg_binomial_2 = 7) + + if (family$family %in% c("gaussian", "Gamma")) { + is_continuous <- TRUE + y_real <- y + y_int <- array(0, dim = c(0)) + } + else { + is_continuous <- FALSE + y_real <- array(0, dim = c(0)) + y_int <- y + } + + if (family$family %in% c("binomial", "poisson")) + has_aux <- FALSE + else + has_aux <- TRUE + + if (family$family != "binomial") + trials <- array(0, dim = c(0)) + + if (family$family %in% c("binomial", "poisson", "neg_binomial_2")) { + if(!is.integer(y_int)) + stop("Outcome must be an integer for count likelihoods.") + if (family$family == "binomial" & (is.null(trials) | any(y > trials))) + stop("Outcome values must be less than or equal to the corresponding value in `trials`.") + } + + if (stan_function == "stan_besag") + model_type <- 1 + else if (stan_function == "stan_bym") + model_type <- 2 + else if (stan_function == "stan_bym2") + model_type <- 3 + + if (!(order %in% c(1,2))) + stop("Argument 'order' must be 1 or 2.") + + x_stuff <- center_x(x, sparse = FALSE) + for (i in names(x_stuff)) # xtemp, xbar, has_intercept + assign(i, x_stuff[[i]]) + nvars <- ncol(xtemp) + + ok_dists <- nlist("normal", student_t = "t", "cauchy", "hs", "hs_plus", + "laplace", "lasso", "product_normal") + ok_intercept_dists <- ok_dists + ok_scale_dists <- nlist("normal", student_t = "t", "cauchy", "exponential") + + # Deal with prior_intercept + prior_intercept_stuff <- handle_glm_prior(prior_intercept, nvars = 1, + default_scale = 10, link = family$link, + ok_dists = ok_intercept_dists) + names(prior_intercept_stuff) <- paste0(names(prior_intercept_stuff), + "_for_intercept") + for (i in names(prior_intercept_stuff)) # prior_{dist, mean, scale, df, autoscale}_for_intercept + assign(i, prior_intercept_stuff[[i]]) + + # Deal with prior + prior_stuff <- handle_glm_prior(prior, nvars, family$link, default_scale = 2.5, + ok_dists = ok_dists) + for (i in names(prior_stuff)) # prior_{dist, mean, scale, df, autoscale} + assign(i, prior_stuff[[i]]) + + # Deal with prior_tau + prior_tau_stuff <- handle_glm_prior(prior_tau, nvars = 1, family$link, default_scale = 1, + ok_dists = ok_scale_dists) + names(prior_tau_stuff) <- paste0(names(prior_tau_stuff), + "_for_tau") + for (i in names(prior_tau_stuff)) # prior_{dist, mean, scale, df, autoscale} + assign(i, prior_tau_stuff[[i]]) + + # Deal with prior_rho + if (stan_function == "stan_bym2") { + has_rho <- 1 + if (is.null(prior_rho)) { + prior_dist_for_rho <- 0 + prior_rho$alpha <- 0 + prior_rho$beta <- 0 + prior_dist_name_for_rho <- NA + } + else { + prior_dist_for_rho <- 1 + prior_dist_name_for_rho <- "beta" + } + prior_rho_stuff <- handle_glm_prior(NULL, nvars = 1, family$link, default_scale = 1, + ok_dists = ok_scale_dists) + names(prior_rho_stuff) <- paste0(names(prior_rho_stuff), + "_for_rho") + for (i in names(prior_rho_stuff)) # prior_{dist, mean, scale, df, autoscale} + assign(i, prior_rho_stuff[[i]]) + prior_rho_stuff$shape1 <- prior_rho$alpha + prior_rho_stuff$shape2 <- prior_rho$beta + } + else if (stan_function == "stan_besag") { + has_rho <- 0 + prior_dist_for_rho <- 0 + # prior_rho_stuff <- list(prior_dist_name_for_rho = NA) + # prior_scale_for_rho <- 0 + # prior_rho_stuff$shape1 <- 0 + # prior_rho_stuff$shape2 <- 0 + prior_rho_stuff <- handle_glm_prior(NULL, nvars = 1, family$link, default_scale = 1, + ok_dists = ok_scale_dists) + names(prior_rho_stuff) <- paste0(names(prior_rho_stuff), + "_for_rho") + for (i in names(prior_rho_stuff)) # prior_{dist, mean, scale, df, autoscale} + assign(i, prior_rho_stuff[[i]]) + prior_rho_stuff$shape1 <- 0 + prior_rho_stuff$shape2 <- 0 + } + else if (stan_function == "stan_bym") { + has_rho <- 1 + prior_rho_stuff <- handle_glm_prior(prior_rho, nvars = 1, family$link, default_scale = 1, + ok_dists = ok_scale_dists) + names(prior_rho_stuff) <- paste0(names(prior_rho_stuff), + "_for_rho") + for (i in names(prior_rho_stuff)) # prior_{dist, mean, scale, df, autoscale} + assign(i, prior_rho_stuff[[i]]) + prior_rho_stuff$shape1 <- 0 + prior_rho_stuff$shape2 <- 0 + } + + # deal with auxiliary parameter + if (has_aux) { + prior_aux_stuff <- handle_glm_prior(prior_aux, nvars = 1, family$link, default_scale = 1, + ok_dists = ok_dists) + names(prior_aux_stuff) <- paste0(names(prior_aux_stuff), + "_for_aux") + for (i in names(prior_aux_stuff)) # prior_{dist, mean, scale, df, autoscale} + assign(i, prior_aux_stuff[[i]]) + } + else { + prior_dist_for_aux <- 0 + prior_mean_for_aux <- 0 + prior_scale_for_aux <- 1 + prior_df_for_aux <- 1 + } + + # QR decomposition for x + if (QR) { + if (nvars <= 1) + stop("'QR' can only be specified when there are multiple predictors.") + else { + cn <- colnames(xtemp) + decomposition <- qr(xtemp) + sqrt_nm1 <- sqrt(nrow(xtemp) - 1L) + Q <- qr.Q(decomposition) + R_inv <- qr.solve(decomposition, Q) * sqrt_nm1 + xtemp <- Q * sqrt_nm1 + colnames(xtemp) <- cn + xbar <- c(xbar %*% R_inv) + } + } + + # need to use uncentered version + standata <- nlist(N = nrow(xtemp), + K = ncol(xtemp), + edges = edges, + E_n = nrow(edges), + family = family_num, + link = link, + is_continuous = is_continuous, + has_aux = has_aux, + X = xtemp, + xbar = as.array(xbar), + y_real = y_real, + y_int = y_int, + trials = trials, + shape1_rho = c(prior_rho_stuff$shape1), + shape2_rho = c(prior_rho_stuff$shape2), + prior_dist_for_intercept = prior_dist_for_intercept, + prior_dist = prior_dist, + prior_dist_rho = prior_dist_for_rho, + prior_dist_tau = prior_dist_for_tau, + prior_dist_for_aux = prior_dist_for_aux, + prior_mean_for_intercept = c(prior_mean_for_intercept), + prior_scale_for_intercept = c(prior_scale_for_intercept), + prior_df_for_intercept = c(prior_df_for_intercept), + prior_mean = as.array(prior_mean), + prior_scale = as.array(prior_scale), + prior_df = as.array(prior_df), + prior_mean_rho = c(prior_mean_for_rho), + prior_scale_rho = c(prior_scale_for_rho), + prior_df_rho = c(prior_df_for_rho), + prior_mean_tau = c(prior_mean_for_tau), + prior_scale_tau = c(prior_scale_for_tau), + prior_df_tau = c(prior_df_for_tau), + prior_mean_for_aux = c(prior_mean_for_aux), + prior_scale_for_aux = c(prior_scale_for_aux), + prior_df_for_aux = c(prior_df_for_aux), + has_intercept = has_intercept, + model_type = model_type, + global_prior_df, + global_prior_df_for_intercept, + global_prior_scale, + global_prior_scale_for_intercept, + num_normals = if(prior_dist == 7) as.integer(prior_df) else integer(0)) + + if (stan_function == "stan_bym2") + standata$scaling_factor <- create_scaling_factor(standata) + else + standata$scaling_factor <- 0 + + standata$order <- order + if (order == 2) { + Q <- Matrix::diag(Matrix::rowSums(w)) - w + Q <- Q %*% Q + sparse_stuff <- rstan::extract_sparse_parts(Q) + standata$Q_n <- as.array(length(sparse_stuff$w), dim = 1) + standata$w <- sparse_stuff$w + standata$v <- sparse_stuff$v + standata$u <- sparse_stuff$u + } + if (order == 1) { + standata$Q_n <- array(0, dim = c(0)) + standata$w <- array(0, dim = c(0)) + standata$v <- array(0, dim = c(0)) + standata$u <- array(0, dim = c(0)) + } + + pars <- c(if (has_intercept) "alpha", "beta", if(model_type != 1) "rho", "tau", if(has_aux) "aux", + "mean_PPD", "psi") + + switch_aux <- switch(family$family, + gaussian = "sigma", + poisson = NA, + neg_binomial_2 = "reciprocal_dispersion", + binomial = NA, + Gamma = "shape") + + prior_info <- summarize_spatial_prior( + user_prior = prior_stuff, + user_prior_intercept = prior_intercept_stuff, + user_prior_rho = prior_rho_stuff, + user_prior_aux = if (has_aux == 1) {prior_aux_stuff} else {NULL}, + user_prior_tau = prior_tau_stuff, + has_intercept = has_intercept, + has_predictors = nvars > 0, + has_aux = has_aux, + has_rho = has_rho, + has_tau = 1, + adjusted_prior_scale = prior_scale, + adjusted_prior_intercept_scale = prior_scale_for_intercept, + adjusted_prior_aux_scale = prior_scale_for_aux, + adjusted_prior_scale_tau = prior_scale_for_tau, + family = family) + + stanfit <- stanmodels$spatial + + # n.b. optimizing is not supported + if (algorithm == "optimizing") { + out <- + optimizing(stanfit, + data = standata, + draws = 1000, + constrained = TRUE, + hessian = TRUE, + ...) + check_stanfit(out) + out$par <- out$par[!grepl("(phi_raw|theta_raw)", names(out$par))] + new_names <- names(out$par) + mark <- grepl("^beta\\[[[:digit:]]+\\]$", new_names) + if (QR && ncol(xtemp) > 1) { + out$par[mark] <- R_inv %*% out$par[mark] + out$theta_tilde[,mark] <- out$theta_tilde[, mark] %*% t(R_inv) + } + new_names[mark] <- colnames(xtemp) + new_names[new_names == "alpha[1]"] <- "(Intercept)" + names(out$par) <- new_names + out$stanfit <- suppressMessages(sampling(stanfit, data = standata, chains = 0)) + return(structure(out, prior.info = prior_info)) + } + else { + if (algorithm == "sampling") { + sampling_args <- set_sampling_args( + object = stanfit, + prior = prior, + user_dots = list(...), + user_adapt_delta = adapt_delta, + data = standata, + pars = pars, + show_messages = FALSE) + stanfit <- do.call(sampling, sampling_args) + } + else { # algorithm either "meanfield" or "fullrank" + stanfit <- rstan::vb(stanfit, pars = pars, data = standata, + algorithm = algorithm, init = 0.001, ...) + if (!QR) + recommend_QR_for_vb() + } + check_stanfit(stanfit) + + if (QR) { + thetas <- extract(stanfit, pars = "beta", inc_warmup = TRUE, + permuted = FALSE) + betas <- apply(thetas, 1:2, FUN = function(theta) R_inv %*% theta) + end <- tail(dim(betas), 1L) + for (chain in 1:end) for (param in 1:nrow(betas)) { + stanfit@sim$samples[[chain]][[has_intercept + param]] <- + if (ncol(xtemp) > 1) betas[param, , chain] else betas[param, chain] + } + } + new_names <- c(if (has_intercept) "(Intercept)", + colnames(xtemp), + if(model_type == 1) {"structured"}, # tau + if(model_type == 2) {c("structured", "unstructured")}, # rho, tau + if(model_type == 3) {c("mixing", "structured")}, # rho, tau + if(has_aux) switch_aux, "mean_PPD", paste0("psi[", 1:standata$N, "]"), "log-posterior") + stanfit@sim$fnames_oi <- new_names + return(structure(stanfit, prior.info = prior_info)) + } +} + +# create scaling_factor a la Dan Simpson +create_scaling_factor <- function(dat) { + edges <- dat$edges + # Build the adjacency matrix + adj.matrix <- Matrix::sparseMatrix(i=edges[,1],j=edges[,2],x=1,symmetric=TRUE) + # The ICAR precision matrix (note! This is singular) + Q <- Matrix::Diagonal(dat$N, Matrix::rowSums(adj.matrix)) - adj.matrix + # Add a small jitter to the diagonal for numerical stability (optional but recommended) + Q_pert <- Q + Matrix::Diagonal(dat$N) * max(Matrix::diag(Q)) * sqrt(.Machine$double.eps) + # Compute the diagonal elements of the covariance matrix subject to the + # constraint that the entries of the ICAR sum to zero. + # See the function help for further details. + # Q_inv <- INLA::inla.qinv(Q_pert, constr=list(A = matrix(1,1,dat$N),e=0)) + Q_inv <- qinv(Q_pert, A = matrix(1,1,dat$N)) + + # Compute the geometric mean of the variances, which are on the diagonal of Q.inv + scaling_factor <- exp(mean(log(Matrix::diag(Q_inv)))) + return(scaling_factor) +} + +# qinv function (analagous to inla.qinv) + +qinv <- function(Q, A = NULL) { + # need to replace the line below with the sparse version, using recursions + Sigma <- Matrix::solve(Q) + if (is.null(A)) + return(Sigma) + else { + A <- matrix(1,1, nrow(Sigma)) + W <- Sigma %*% t(A) + Sigma_const <- Sigma - W %*% solve(A %*% W) %*% t(W) + return(Sigma_const) + } +} + +# Summarize spatial prior + +summarize_spatial_prior <- function(user_prior, + user_prior_intercept, + user_prior_aux, + user_prior_rho, + user_prior_tau, + has_intercept, + has_predictors, + has_aux, + has_rho, + has_tau, + adjusted_prior_scale, + adjusted_prior_intercept_scale, + adjusted_prior_scale_rho, + adjusted_prior_scale_tau, + adjusted_prior_aux_scale, + family) { + rescaled_coef <- + user_prior$prior_autoscale && has_predictors && + !is.na(user_prior$prior_dist_name) && + !all(user_prior$prior_scale == adjusted_prior_scale) + rescaled_int <- + user_prior_intercept$prior_autoscale_for_intercept && has_intercept && + !is.na(user_prior_intercept$prior_dist_name_for_intercept) && + (user_prior_intercept$prior_scale != adjusted_prior_intercept_scale) + if (has_aux) { + rescaled_aux <- user_prior_aux$prior_autoscale_for_aux && + !is.na(user_prior_aux$prior_dist_name_for_aux) && + (user_prior_aux$prior_scale_for_aux != adjusted_prior_aux_scale) + } + + if (has_predictors && user_prior$prior_dist_name %in% "t") { + if (all(user_prior$prior_df == 1)) { + user_prior$prior_dist_name <- "cauchy" + } else { + user_prior$prior_dist_name <- "student_t" + } + } + if (has_intercept && + user_prior_intercept$prior_dist_name_for_intercept %in% "t") { + if (all(user_prior_intercept$prior_df_for_intercept == 1)) { + user_prior_intercept$prior_dist_name_for_intercept <- "cauchy" + } else { + user_prior_intercept$prior_dist_name_for_intercept <- "student_t" + } + } + if (has_aux && + user_prior_aux$prior_dist_name_for_aux %in% "t") { + if (all(user_prior_aux$prior_df_for_aux == 1)) { + user_prior_aux$prior_dist_name_for_aux <- "cauchy" + } else { + user_prior_aux$prior_dist_name_for_aux <- "student_t" + } + } + + if (has_rho && + user_prior_rho$prior_dist_name_for_rho %in% "t") { + if (all(user_prior_rho$prior_df_for_rho == 1)) { + user_prior_rho$prior_dist_name_for_rho <- "cauchy" + } else if (has_rho && user_prior_rho$prior_dist_name_for_rho == "beta") { + user_prior_rho$prior_dist_name_for_rho <- "beta" + } else { + user_prior_rho$prior_dist_name_for_rho <- "student_t" + } + } + if (has_tau && + user_prior_tau$prior_dist_name_for_tau %in% "t") { + if (all(user_prior_tau$prior_df_for_tau == 1)) { + user_prior_tau$prior_dist_name_for_tau <- "cauchy" + } else { + user_prior_tau$prior_dist_name_for_tau <- "student_t" + } + } + prior_list <- list( + prior = + if (!has_predictors) NULL else with(user_prior, list( + dist = prior_dist_name, + location = prior_mean, + scale = prior_scale, + adjusted_scale = if (rescaled_coef) + adjusted_prior_scale else NULL, + df = if (prior_dist_name %in% c("student_t", "hs", "hs_plus", + "lasso", "product_normal")) + prior_df else NULL + )), + prior_intercept = + if (!has_intercept) NULL else with(user_prior_intercept, list( + dist = prior_dist_name_for_intercept, + location = prior_mean_for_intercept, + scale = prior_scale_for_intercept, + adjusted_scale = if (rescaled_int) + adjusted_prior_intercept_scale else NULL, + df = if (prior_dist_name_for_intercept %in% "student_t") + prior_df_for_intercept else NULL + )), + prior_aux = + if (!has_aux) NULL else with(user_prior_aux, list( + dist = prior_dist_name_for_aux, + location = prior_mean_for_aux, + scale = prior_scale_for_aux, + adjusted_scale = if (rescaled_int) + adjusted_prior_aux_scale else NULL, + df = if (prior_dist_name_for_aux %in% "student_t") + prior_df_for_aux else NULL + )), + prior_rho = + if (!has_rho) NULL else with(user_prior_rho, list( + dist = prior_dist_name_for_rho, + location = prior_mean_for_rho, + scale = prior_scale_for_rho, + shape1 = shape1, + shape2 = shape2, + adjusted_scale = if (rescaled_int) + adjusted_prior_rho_scale else NULL, + df = if (prior_dist_name_for_rho %in% "student_t") + prior_df_for_rho else NULL + )), + prior_tau = + if (!has_tau) NULL else with(user_prior_tau, list( + dist = prior_dist_name_for_tau, + location = prior_mean_for_tau, + scale = prior_scale_for_tau, + adjusted_scale = if (rescaled_int) + adjusted_prior_tau_scale else NULL, + df = if (prior_dist_name_for_tau %in% "student_t") + prior_df_for_tau else NULL + )) + ) + aux_name <- .rename_aux(family) + prior_list$prior_aux <- if (is.na(aux_name)) + NULL else with(user_prior_aux, list( + dist = prior_dist_name_for_aux, + location = if (!is.na(prior_dist_name_for_aux) && + prior_dist_name_for_aux != "exponential") + prior_mean_for_aux else NULL, + scale = if (!is.na(prior_dist_name_for_aux) && + prior_dist_name_for_aux != "exponential") + prior_scale_for_aux else NULL, + adjusted_scale = if (rescaled_aux) + adjusted_prior_aux_scale else NULL, + df = if (!is.na(prior_dist_name_for_aux) && + prior_dist_name_for_aux %in% "student_t") + prior_df_for_aux else NULL, + rate = if (!is.na(prior_dist_name_for_aux) && + prior_dist_name_for_aux %in% "exponential") + 1 / prior_scale_for_aux else NULL, + aux_name = aux_name + )) + return(prior_list) +} diff --git a/R/stanreg-methods.R b/R/stanreg-methods.R index fa1a9271a..6e747db67 100644 --- a/R/stanreg-methods.R +++ b/R/stanreg-methods.R @@ -452,7 +452,6 @@ terms.stanreg <- function(x, ..., fixed.only = TRUE, random.only = FALSE) { fixed.only <- FALSE if (fixed.only && random.only) stop("'fixed.only' and 'random.only' can't both be TRUE.", call. = FALSE) - Terms <- attr(fr, "terms") if (fixed.only) { Terms <- terms.formula(formula(x, fixed.only = TRUE)) diff --git a/R/stanreg.R b/R/stanreg.R index b4bce2638..2bdc41796 100644 --- a/R/stanreg.R +++ b/R/stanreg.R @@ -31,6 +31,7 @@ stanreg <- function(object) { nobs <- NROW(y) ynames <- if (is.matrix(y)) rownames(y) else names(y) + is_car <- object$stan_function %in% c("stan_besag", "stan_bym", "stan_bym2") is_betareg <- is.beta(family$family) if (is_betareg) { family_phi <- object$family_phi # pull out phi family/link @@ -44,6 +45,7 @@ stanreg <- function(object) { if (opt) { stanmat <- stanfit$theta_tilde probs <- c(0.025, .975) + stan_summary <- cbind(Median = apply(stanmat, 2L, median), MAD_SD = apply(stanmat, 2L, mad), t(apply(stanmat, 2L, quantile, probs))) @@ -90,8 +92,13 @@ stanreg <- function(object) { # linear predictor, fitted values eta <- linear_predictor(coefs, x, object$offset) + if (is_car) { + psi_indx <- grep("psi", colnames(as.matrix(object$stanfit))) + psi <- as.matrix(object$stanfit)[,psi_indx] + psi <- unname(colMeans(psi)) + eta <- eta + psi + } mu <- family$linkinv(eta) - if (NCOL(y) == 2L) { # residuals of type 'response', (glm which does 'deviance' residuals by default) residuals <- y[, 1L] / rowSums(y) - mu @@ -104,7 +111,7 @@ stanreg <- function(object) { eta_z <- linear_predictor(coefs_z, z, object$offset) phi <- family_phi$linkinv(eta_z) } - + out <- nlist( coefficients = unpad_reTrms(coefs), ses = unpad_reTrms(ses), @@ -149,6 +156,10 @@ stanreg <- function(object) { out$eta_z <- eta_z out$phi <- phi } + if (is_car) { + out$psi <- psi + out$trials <- object$trials + } structure(out, class = c("stanreg", "glm", "lm")) } diff --git a/cleanup b/cleanup index 3b59622ab..b1e2a1e77 100755 --- a/cleanup +++ b/cleanup @@ -1,5 +1,5 @@ #!/bin/sh -e # Note to Windows users: This is not actually platform specific. -# "${R_HOME}/bin/R" --vanilla --slave -e 'roxygen2::roxygenize(clean = FALSE)' +# "${R_HOME}/bin/R" --vanilla --slave -e 'roxygen2::roxygenize(clean = TRUE)' exit $? diff --git a/data/lattice.rda b/data/lattice.rda new file mode 100644 index 000000000..fecb830cf Binary files /dev/null and b/data/lattice.rda differ diff --git a/src/Makevars.win b/src/Makevars.win index 7750ec0ff..7a6e362a7 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -23,5 +23,5 @@ clean: %.cc: %.stan "$(R_HOME)/bin$(R_ARCH_BIN)/Rscript" -e "source(file.path('..', 'tools', 'make_cc.R')); make_cc(commandArgs(TRUE))" $< - + .phony: clean diff --git a/src/stan_files/continuous.stan b/src/stan_files/continuous.stan index 360c2b2e9..37f6eca19 100644 --- a/src/stan_files/continuous.stan +++ b/src/stan_files/continuous.stan @@ -35,9 +35,9 @@ functions { * @param b Vector that is multiplied from the left by the CSR matrix * @return A vector that is the product of the CSR matrix and b */ - vector test_csr_matrix_times_vector(int m, int n, vector w, + vector test_csr_matrix_times_vector(int m, int n, vector w, int[] v, int[] u, vector b) { - return csr_matrix_times_vector(m, n, w, v, u, b); + return csr_matrix_times_vector(m, n, w, v, u, b); } } data { @@ -115,7 +115,7 @@ parameters { } transformed parameters { // aux has to be defined first in the hs case - real aux = prior_dist_for_aux == 0 ? aux_unscaled : (prior_dist_for_aux <= 2 ? + real aux = prior_dist_for_aux == 0 ? aux_unscaled : (prior_dist_for_aux <= 2 ? prior_scale_for_aux * aux_unscaled + prior_mean_for_aux : prior_scale_for_aux * aux_unscaled); @@ -123,7 +123,7 @@ transformed parameters { // defines beta, b, theta_L #include /tparameters/tparameters_glm.stan #include /tparameters/tparameters_betareg.stan - + if (prior_dist_for_aux == 0) // none aux = aux_unscaled; else { @@ -144,7 +144,7 @@ transformed parameters { } } else { - theta_L = make_theta_L(len_theta_L, p, + theta_L = make_theta_L(len_theta_L, p, aux, tau, scale, zeta, rho, z_T); b = make_b(z_b, theta_L, p, l); } @@ -243,15 +243,6 @@ model { rows_dot_product((1 - mu) , mu_z)); } } - else { // weighted log-likelihoods - vector[N] summands; - if (family == 1) summands = pw_gauss(y, eta, aux, link); - else if (family == 2) summands = pw_gamma(y, eta, aux, link); - else if (family == 3) summands = pw_inv_gaussian(y, eta, aux, link, log_y, sqrt_y); - else if (family == 4 && link_phi == 0) summands = pw_beta(y, eta, aux, link); - else if (family == 4 && link_phi > 0) summands = pw_beta_z(y, eta, eta_z, link, link_phi); - target += dot_product(weights, summands); - } } // Log-priors @@ -261,14 +252,14 @@ model { target += normal_lpdf(aux_unscaled | 0, 1) - log_half; else if (prior_dist_for_aux == 2) target += student_t_lpdf(aux_unscaled | prior_df_for_aux, 0, 1) - log_half; - else + else target += exponential_lpdf(aux_unscaled | 1); } - + #include /model/priors_glm.stan #include /model/priors_betareg.stan if (t > 0) { - real dummy = decov_lp(z_b, z_T, rho, zeta, tau, + real dummy = decov_lp(z_b, z_T, rho, zeta, tau, regularization, delta, shape, t, p); } } @@ -276,15 +267,15 @@ generated quantities { real mean_PPD = compute_mean_PPD ? 0 : negative_infinity(); real alpha[has_intercept]; real omega_int[has_intercept_z]; - + if (has_intercept == 1) { if (dense_X) alpha[1] = gamma[1] - dot_product(xbar, beta); else alpha[1] = gamma[1]; } - if (has_intercept_z == 1) { - omega_int[1] = gamma_z[1] - dot_product(zbar, omega); // adjust betareg intercept + if (has_intercept_z == 1) { + omega_int[1] = gamma_z[1] - dot_product(zbar, omega); // adjust betareg intercept } - + if (compute_mean_PPD) { vector[N] eta_z; #include /model/make_eta.stan @@ -322,10 +313,10 @@ generated quantities { else { // has_intercept_z == 0 #include /model/eta_z_no_intercept.stan } - + if (SSfun > 0) { // nlmer vector[len_y] eta_nlmer; - matrix[len_y, K] P; + matrix[len_y, K] P; P = reshape_vec(eta, len_y, K); if (SSfun < 5) { if (SSfun <= 2) { @@ -361,7 +352,7 @@ generated quantities { vector[N] mu = link > 1 ? linkinv_inv_gaussian(eta, link) : eta; for (n in 1:len_y) mean_PPD += inv_gaussian_rng(mu[n], aux); } - else if (family == 4 && link_phi == 0) { + else if (family == 4 && link_phi == 0) { vector[N] mu = linkinv_beta(eta, link); for (n in 1:N) { real mu_n = mu[n]; diff --git a/src/stan_files/functions/common_functions.stan b/src/stan_files/functions/common_functions.stan index b288c1902..2e38a19d6 100644 --- a/src/stan_files/functions/common_functions.stan +++ b/src/stan_files/functions/common_functions.stan @@ -1,13 +1,13 @@ /* for multiple .stan files */ - - /** + + /** * Create group-specific block-diagonal Cholesky factor, see section 2 of * https://cran.r-project.org/web/packages/lme4/vignettes/lmer.pdf - * @param len_theta_L An integer indicating the length of returned vector, + * @param len_theta_L An integer indicating the length of returned vector, * which lme4 denotes as m * @param p An integer array with the number variables on the LHS of each | * @param dispersion Scalar standard deviation of the errors, calles sigma by lme4 - * @param tau Vector of scale parameters whose squares are proportional to the + * @param tau Vector of scale parameters whose squares are proportional to the * traces of the relative covariance matrices of the group-specific terms * @param scale Vector of prior scales that are multiplied by elements of tau * @param zeta Vector of positive parameters that are normalized into simplexes @@ -26,33 +26,33 @@ int theta_L_mark = 1; // each of these is a diagonal block of the implicit Cholesky factor - for (i in 1:size(p)) { + for (i in 1:size(p)) { int nc = p[i]; if (nc == 1) { // "block" is just a standard deviation theta_L[theta_L_mark] = tau[i] * scale[i] * dispersion; // unlike lme4, theta[theta_L_mark] includes the dispersion term in it theta_L_mark += 1; } - else { // block is lower-triangular - matrix[nc,nc] T_i; + else { // block is lower-triangular + matrix[nc,nc] T_i; real std_dev; real T21; real trace_T_i = square(tau[i] * scale[i] * dispersion) * nc; vector[nc] pi = segment(zeta, zeta_mark, nc); // gamma(zeta | shape, 1) pi /= sum(pi); // thus dirichlet(pi | shape) - + // unlike lme4, T_i includes the dispersion term in it zeta_mark += nc; std_dev = sqrt(pi[1] * trace_T_i); T_i[1,1] = std_dev; - + // Put a correlation into T_i[2,1] and scale by std_dev std_dev = sqrt(pi[2] * trace_T_i); T21 = 2.0 * rho[rho_mark] - 1.0; rho_mark += 1; T_i[2,2] = std_dev * sqrt(1.0 - square(T21)); T_i[2,1] = std_dev * T21; - + for (r in 2:(nc - 1)) { // scaled onion method to fill T_i int rp1 = r + 1; vector[r] T_row = segment(z_T, z_T_mark, r); @@ -63,7 +63,7 @@ T_i[rp1,rp1] = sqrt(1.0 - rho[rho_mark]) * std_dev; rho_mark += 1; } - + // now vech T_i for (c in 1:nc) for (r in c:nc) { theta_L[theta_L_mark] = T_i[r,c]; @@ -73,15 +73,15 @@ } return theta_L; } - - /** + + /** * Create group-specific coefficients, see section 2 of * https://cran.r-project.org/web/packages/lme4/vignettes/lmer.pdf * * @param z_b Vector whose elements are iid normal(0,sigma) a priori * @param theta Vector with covariance parameters as defined in lme4 * @param p An integer array with the number variables on the LHS of each | - * @param l An integer array with the number of levels for the factor(s) on + * @param l An integer array with the number of levels for the factor(s) on * the RHS of each | * @return A vector of group-specific coefficients */ @@ -93,7 +93,7 @@ int nc = p[i]; if (nc == 1) { real theta_L_start = theta_L[theta_L_mark]; - for (s in b_mark:(b_mark + l[i] - 1)) + for (s in b_mark:(b_mark + l[i] - 1)) b[s] = theta_L_start * z_b[s]; b_mark += l[i]; theta_L_mark += 1; @@ -119,7 +119,7 @@ return b; } - /** + /** * Prior on group-specific parameters * * @param z_b A vector of primitive coefficients @@ -160,7 +160,7 @@ target += gamma_lpdf(tau | shape, 1); return target(); } - + /** * Hierarchical shrinkage parameterization * @@ -172,7 +172,7 @@ * @param c2 A positive real number * @return A vector of coefficientes */ - vector hs_prior(vector z_beta, real[] global, vector[] local, + vector hs_prior(vector z_beta, real[] global, vector[] local, real global_prior_scale, real error_scale, real c2) { int K = rows(z_beta); vector[K] lambda = local[1] .* sqrt(local[2]); @@ -182,7 +182,7 @@ return z_beta .* lambda_tilde * tau; } - /** + /** * Hierarchical shrinkage plus parameterization * * @param z_beta A vector of primitive coefficients @@ -193,19 +193,19 @@ * @param c2 A positive real number * @return A vector of coefficientes */ - vector hsplus_prior(vector z_beta, real[] global, vector[] local, + vector hsplus_prior(vector z_beta, real[] global, vector[] local, real global_prior_scale, real error_scale, real c2) { int K = rows(z_beta); vector[K] lambda = local[1] .* sqrt(local[2]); vector[K] eta = local[3] .* sqrt(local[4]); real tau = global[1] * sqrt(global[2]) * global_prior_scale * error_scale; vector[K] lambda_eta2 = square(lambda .* eta); - vector[K] lambda_tilde = sqrt( c2 * lambda_eta2 ./ + vector[K] lambda_tilde = sqrt( c2 * lambda_eta2 ./ ( c2 + square(tau) * lambda_eta2) ); return z_beta .* lambda_tilde * tau; } - - /** + + /** * Cornish-Fisher expansion for standard normal to Student t * * See result 26.7.5 of @@ -229,7 +229,7 @@ + (79 * z9 + 776 * z7 + 1482 * z5 - 1920 * z3 - 945 * z) / (92160 * df4); } - /** + /** * Return two-dimensional array of group membership * * @param N An integer indicating the number of observations @@ -247,7 +247,7 @@ return V; } - /** + /** * faster version of csr_matrix_times_vector * declared here and defined in C++ * @@ -259,7 +259,7 @@ * @param b Vector that is multiplied from the left by the CSR matrix * @return A vector that is the product of the CSR matrix and b */ - vector csr_matrix_times_vector2(int m, int n, vector w, + vector csr_matrix_times_vector2(int m, int n, vector w, int[] v, int[] u, vector b); /** @@ -283,6 +283,9 @@ if (link == 2) return negative_infinity(); // log return 0; } + if (family == 6 && link != 1) return 0; // poisson + if (family == 7 && link != 1) return 0; // neg_binomial_2 + if (family == 2 && (link == 1 || link == 3)) return 0; // gamma(link="idenity|inverse") return negative_infinity(); } @@ -295,5 +298,6 @@ */ real make_upper(int family, int link) { if (family == 4 && link == 5) return 0; + if (family == 5 && link == 4) return 0; // binomial(link="log") return positive_infinity(); } diff --git a/src/stan_files/model/priors_glm.stan b/src/stan_files/model/priors_glm.stan index ffa775582..4e92d5ba5 100644 --- a/src/stan_files/model/priors_glm.stan +++ b/src/stan_files/model/priors_glm.stan @@ -35,26 +35,16 @@ target += normal_lpdf(z_beta | 0, 1); } /* else prior_dist is 0 and nothing is added */ - - // Log-prior for intercept - if (has_intercept == 1) { - if (prior_dist_for_intercept == 1) // normal - target += normal_lpdf(gamma | prior_mean_for_intercept, prior_scale_for_intercept); - else if (prior_dist_for_intercept == 2) // student_t - target += student_t_lpdf(gamma | prior_df_for_intercept, prior_mean_for_intercept, - prior_scale_for_intercept); - /* else prior_dist is 0 and nothing is added */ - } if (K_smooth) { target += normal_lpdf(z_beta_smooth | 0, 1); if (prior_dist_for_smooth > 0) { real log_half = -0.693147180559945286; - if (prior_dist_for_smooth == 1) + if (prior_dist_for_smooth == 1) target += normal_lpdf(smooth_sd_raw | 0, 1) - log_half; else if (prior_dist_for_smooth == 2) target += student_t_lpdf(smooth_sd_raw | prior_df_for_smooth, 0, 1) - log_half; - else if (prior_dist_for_smooth == 3) + else if (prior_dist_for_smooth == 3) target += exponential_lpdf(smooth_sd_raw | 1); } } diff --git a/src/stan_files/spatial.stan b/src/stan_files/spatial.stan new file mode 100644 index 000000000..24e9c8d8f --- /dev/null +++ b/src/stan_files/spatial.stan @@ -0,0 +1,325 @@ +#include /pre/Columbia_copyright.stan +#include /pre/license.stan +// CAR SPATIAL MODELS +functions { +#include /functions/continuous_likelihoods.stan +#include /functions/binomial_likelihoods.stan +#include /functions/count_likelihoods.stan +#include /functions/common_functions.stan +} +data { + int N; // number of regions + int K; // number of predictors (inc intercept) + matrix[N,K] X; // model matrix + int family; // 1 gaussian; 2 gamma; 5 binomial; 6 poisson; 7 neg_binomial_2 + int link; + int is_continuous; + int has_aux; + int model_type; // Besag = 1; BYM = 2; BYM2 = 3 + int has_intercept; + vector[K] xbar; + int trials[family == 5 ? N : 0]; + int y_int[is_continuous == 1 ? 0 : N]; + vector[is_continuous == 1 ? N : 0] y_real; + // pairwise difference version of CAR + int E_n; // number of adjacency pairs + int edges[E_n, 2]; // adjacency pairs + // matrix version of CAR + int order; + int Q_n[order == 2]; + vector[order == 2 ? Q_n[1] : 0] w; + int v[order == 2 ? Q_n[1] : 0]; + int u[order == 2 ? N+1 : 0]; + // prior stuff + int prior_dist_rho; + real prior_mean_rho; + real prior_scale_rho; + real prior_df_rho; + real shape1_rho; + real shape2_rho; + real scaling_factor; + int prior_dist_for_intercept; + int prior_dist; + int prior_dist_tau; + real prior_mean_for_intercept; + real prior_scale_for_intercept; + real prior_df_for_intercept; + vector[K] prior_mean; + vector[K] prior_scale; + vector[K] prior_df; + real prior_mean_tau; + real prior_scale_tau; + real prior_df_tau; + real global_prior_df; + real global_prior_scale; + int num_normals[prior_dist == 7 ? K : 0]; + int prior_dist_for_aux; + real prior_mean_for_aux; + real prior_scale_for_aux; + real prior_df_for_aux; + real slab_scale; // for hs prior only +} +transformed data { + real poisson_max = pow(2.0, 30.0); + int hs; + real sum_log_y = family == 1 ? not_a_number() : sum(log(y_real)); + if (prior_dist <= 2) hs = 0; + else if (prior_dist == 3) hs = 2; + else if (prior_dist == 4) hs = 4; + else hs = 0; +} +parameters { + real gamma[has_intercept]; + // real gamma[has_intercept]; + vector[K] z_beta; + vector[model_type == 1? 0 : N] theta_raw; + vector[N-1] phi_raw; + // interpretation of rho and tau depends on model_type! + real rho[model_type != 1]; + real tau; + real global[hs]; + vector[K] local[hs]; + vector[K] mix[prior_dist == 5 || prior_dist == 6]; + real aux_unscaled[has_aux]; // interpretation depends on family! + real one_over_lambda[prior_dist == 6]; + real caux[hs > 0]; +} +transformed parameters { + vector[K] beta; // predictors on covariates (including intercept) + // aux has to be defined first in the hs case + real aux = has_aux == 0 ? 0 : (prior_dist_for_aux == 0 ? aux_unscaled[1] : (prior_dist_for_aux <= 2 ? + prior_scale_for_aux * aux_unscaled[1] + prior_mean_for_aux : + prior_scale_for_aux * aux_unscaled[1])); + vector[N] phi; // non-centered random effect (spatial) + vector[N] psi; + phi[1:(N - 1)] = phi_raw; + phi[N] = -sum(phi_raw); + /* precision form + if (model_type == 1) + psi = phi * sqrt(inv(tau)); + else if (model_type == 2) + psi = phi * sqrt(inv(rho[1])) + theta_raw * sqrt(inv(tau)); + else if (model_type == 3) + psi = tau*(sqrt(1-rho[1])*theta_raw + sqrt(rho[1]/scaling_factor)*phi); + */ + if (model_type == 1) + psi = phi * tau; + else if (model_type == 2) + psi = phi * rho[1] + theta_raw * tau; + else if (model_type == 3) + psi = tau*(sqrt(1-rho[1])*theta_raw + sqrt(rho[1]/scaling_factor)*phi); + // for regression coefficients + // "tparameters.stan" + if (prior_dist == 0) beta = z_beta; + else if (prior_dist == 1) beta = z_beta .* prior_scale + prior_mean; + else if (prior_dist == 2) for (k in 1:K) { + beta[k] = CFt(z_beta[k], prior_df[k]) * prior_scale[k] + prior_mean[k]; + } + else if (prior_dist == 3) { + real c2 = square(slab_scale) * caux[1]; + if (is_continuous == 1 && family == 1) + beta = hs_prior(z_beta, global, local, global_prior_scale, aux, c2); + else beta = hs_prior(z_beta, global, local, global_prior_scale, 1, c2); + } + else if (prior_dist == 4) { + real c2 = square(slab_scale) * caux[1]; + if (is_continuous == 1 && family == 1) + beta = hsplus_prior(z_beta, global, local, global_prior_scale, aux, c2); + else beta = hsplus_prior(z_beta, global, local, global_prior_scale, 1, c2); + } + else if (prior_dist == 5) // laplace + beta = prior_mean + prior_scale .* sqrt(2 * mix[1]) .* z_beta; + else if (prior_dist == 6) // lasso + beta = prior_mean + one_over_lambda[1] * prior_scale .* sqrt(2 * mix[1]) .* z_beta; + else if (prior_dist == 7) { // product_normal + int z_pos = 1; + for (k in 1:K) { + beta[k] = z_beta[z_pos]; + z_pos = z_pos + 1; + for (n in 2:num_normals[k]) { + beta[k] = beta[k] * z_beta[z_pos]; + z_pos = z_pos + 1; + } + beta[k] = beta[k] * prior_scale[k] ^ num_normals[k] + prior_mean[k]; + } + } +} +model { + vector[N] eta; // linear predictor + spatial random effects + // deal with intercept + if (has_intercept == 1) { + eta = X * beta + psi; + if ((family == 5 && link == 4)) // binomial + eta -= max(eta); + else if ((family == 6 && link != 1) || + (family == 7 && link != 1) || + (family == 2 && (link == 1 || link == 3))) // poisson, neg_binomial_2, and gamma + eta -= min(eta); + eta += gamma[1]; + } + else + eta = X * beta + psi; + // likelihoods + if (family == 1) { + target+= normal_lpdf(y_real | linkinv_gauss(eta, link), aux); + } + else if (family == 6) { + target+= poisson_lpmf(y_int | linkinv_count(eta, link)); + } + else if (family == 7) { + target += neg_binomial_2_lpmf(y_int | linkinv_count(eta, link), aux); + } + else if (family == 5) { + target+= binomial_lpmf(y_int | trials, linkinv_binom(eta, link)); + } + else if (family == 2) { + target += GammaReg(y_real, eta, aux, link, sum_log_y); + } + // prior on spatial parameter vector (GMRF) + if (order == 1) + target += -0.5 * dot_self(phi[edges[,1]] - phi[edges[,2]]); + else if (order == 2) + target+= -0.5 * dot_product(phi, csr_matrix_times_vector(N, N, w, v, u, phi)); + // priors on coefficients + // "priors.stan" + // Log-priors for coefficients + if (prior_dist == 1) target += normal_lpdf(z_beta | 0, 1); + else if (prior_dist == 2) target += normal_lpdf(z_beta | 0, 1); // Student t + else if (prior_dist == 3) { // hs + real log_half = -0.693147180559945286; + target += normal_lpdf(z_beta | 0, 1); + target += normal_lpdf(local[1] | 0, 1) - log_half; + target += inv_gamma_lpdf(local[2] | 0.5 * prior_df, 0.5 * prior_df); + target += normal_lpdf(global[1] | 0, 1) - log_half; + target += inv_gamma_lpdf(global[2] | 0.5 * global_prior_df, 0.5 * global_prior_df); + } + else if (prior_dist == 4) { // hs+ + real log_half = -0.693147180559945286; + target += normal_lpdf(z_beta | 0, 1); + target += normal_lpdf(local[1] | 0, 1) - log_half; + target += inv_gamma_lpdf(local[2] | 0.5 * prior_df, 0.5 * prior_df); + target += normal_lpdf(local[3] | 0, 1) - log_half; + // unorthodox useage of prior_scale as another df hyperparameter + target += inv_gamma_lpdf(local[4] | 0.5 * prior_scale, 0.5 * prior_scale); + target += normal_lpdf(global[1] | 0, 1) - log_half; + target += inv_gamma_lpdf(global[2] | 0.5 * global_prior_df, 0.5 * global_prior_df); + } + else if (prior_dist == 5) { // laplace + target += normal_lpdf(z_beta | 0, 1); + target += exponential_lpdf(mix[1] | 1); + } + else if (prior_dist == 6) { // lasso + target += normal_lpdf(z_beta | 0, 1); + target += exponential_lpdf(mix[1] | 1); + target += chi_square_lpdf(one_over_lambda[1] | prior_df[1]); + } + else if (prior_dist == 7) { // product_normal + target += normal_lpdf(z_beta | 0, 1); + } + /* else prior_dist is 0 and nothing is added */ + + // Log-prior for intercept + if (has_intercept == 1) { + if (prior_dist_for_intercept == 1) // normal + target += normal_lpdf(gamma | prior_mean_for_intercept, prior_scale_for_intercept); + else if (prior_dist_for_intercept == 2) // student_t + target += student_t_lpdf(gamma | prior_df_for_intercept, prior_mean_for_intercept, + prior_scale_for_intercept); + /* else prior_dist is 0 and nothing is added */ + } + // model specific priors + if (model_type == 2) { + target+= normal_lpdf(theta_raw | 0, 1); // unstructured (random) effect + // prior on overall spatial variation + if (prior_dist_rho == 1) + target+= normal_lpdf(rho[1] | prior_mean_rho, prior_scale_rho); + else if (prior_dist_rho == 2) + target+= student_t_lpdf(rho[1] | prior_df_rho, prior_mean_rho, prior_scale_rho); + else if (prior_dist_rho == 3) + target+= exponential_lpdf(rho[1] | prior_scale_rho); + /* else prior_dist_tau is 0 and nothing is added */ + } + else if (model_type == 3) { // BYM + target+= normal_lpdf(theta_raw | 0, 1); // unstructured (random) effect + if (prior_dist_rho == 1) + target+= beta_lpdf(rho[1] | shape1_rho, shape2_rho); + /* else prior_dist_tau is 0 and nothing is added */ + } + // prior on overall spatial variation + if (prior_dist_tau == 1) + target+= normal_lpdf(tau | prior_mean_tau, prior_scale_tau); + else if (prior_dist_tau == 2) + target+= student_t_lpdf(tau | prior_df_tau, prior_mean_tau, prior_scale_tau); + else if (prior_dist_tau == 3) + target+= exponential_lpdf(tau | prior_scale_tau); + /* else prior_dist_tau is 0 and nothing is added */ + // priors on auxilliary parameters (Log-priors) + if (has_aux == 1) { + // "priors_aux.stan" + // Log-priors + if (prior_dist_for_aux > 0 && prior_scale_for_aux > 0) { + real log_half = -0.693147180559945286; + if (prior_dist_for_aux == 1) + target += normal_lpdf(aux_unscaled | 0, 1) - log_half; + else if (prior_dist_for_aux == 2) + target += student_t_lpdf(aux_unscaled | prior_df_for_aux, 0, 1) - log_half; + else + target += exponential_lpdf(aux_unscaled | 1); + } + } +} +generated quantities { + real mean_PPD = 0; + real alpha[has_intercept]; + { + vector[N] eta = X * beta + psi; + if (has_intercept == 1) { + alpha[1] = gamma[1] - dot_product(beta, xbar); + if ((family == 5 && link == 4)) { // binomial + eta -= max(eta); + alpha[1] -= max(eta); + } + else if ((family == 6 && link != 1) || + (family == 7 && link != 1) || + (family == 2 && (link == 1 || link == 3))) {// poisson, neg_binomial_2, and gamma + eta -= min(eta); + alpha[1] -= min(eta); + } + eta += gamma[1]; + } + if (family == 1) { + eta = linkinv_gauss(eta, link); + for (n in 1:N) mean_PPD += normal_rng(eta[n], aux); + } + else if (family == 6) { + eta = linkinv_count(eta, link); + if (link == 2) + for (n in 1:N) { + if (eta[n] < poisson_max) + mean_PPD += poisson_rng(eta[n]); + else + mean_PPD += normal_rng(eta[n], sqrt(eta[n])); + } + } + else if (family == 7) { + eta = linkinv_count(eta, link); + for (n in 1:N) { + real gamma_temp; + if (is_inf(aux)) gamma_temp = eta[n]; + else gamma_temp = gamma_rng(aux, aux / eta[n]); + if (gamma_temp < poisson_max) + mean_PPD += poisson_rng(gamma_temp); + else mean_PPD += normal_rng(gamma_temp, sqrt(gamma_temp)); + } + } + else if (family == 5) { + vector[N] lp_tf = linkinv_binom(eta, link); + for (n in 1:N) mean_PPD += binomial_rng(trials[n], lp_tf[n]); + } + else if (family == 2) { + if (link > 1) eta = linkinv_gamma(eta, link); + for (n in 1:N) mean_PPD += gamma_rng(aux, aux / eta[n]); + } + } + mean_PPD /= N; +} diff --git a/tests/testthat/test_stan_besag.R b/tests/testthat/test_stan_besag.R new file mode 100644 index 000000000..a8d10d969 --- /dev/null +++ b/tests/testthat/test_stan_besag.R @@ -0,0 +1,116 @@ +# 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. + +# tests can be run using devtools::test() or manually by loading testthat +# package and then running the code below possibly with options(mc.cores = 4). + +library(rstanarm) +SEED <- 12345 +set.seed(SEED) +ITER <- 25 +CHAINS <- 2 + +context("stan_besag") +# for line-byline testing only +# source(paste0("tests/testthat/",(file.path("helpers", "expect_stanreg.R")))) +# source(paste0("tests/testthat/",(file.path("helpers", "SW.R")))) +# for full package testing +source(file.path("helpers", "expect_stanreg.R")) +source(file.path("helpers", "SW.R")) + +data("lattice", package = "rstanarm") + +# Convert a spatial polygon to an N-by-N weight matrix +sp2weightmatrix <- function(spatialpolygon) { + spdep::nb2mat(spdep::poly2nb(spatialpolygon, queen = TRUE), style = "B", zero.policy = TRUE) +} +adj <- sp2weightmatrix(grid_sim15) +spatial_data <- grid_sim15@data + +SW(fit_gauss <- stan_besag(y_gauss ~ 1 + x, data = spatial_data, family = gaussian(link = "identity"), + prior_intercept = normal(0,1), prior = normal(0,1), prior_aux = normal(0,1), + W = adj, iter = ITER, chains = CHAINS)) +SW(fit_binom <- stan_besag(y_binom ~ 1 + x, trials = spatial_data$trials, data = spatial_data, + prior_intercept = normal(0,1), prior = normal(0,1), + family = binomial(link = "logit"), + W = adj, iter = ITER, chains = CHAINS)) +SW(fit_pois <- stan_besag(y_pois ~ 1 + x, data = spatial_data, + prior_intercept = normal(0,1), prior = normal(0,1), + family = poisson(link = "log"), + W = adj, iter = ITER, chains = CHAINS)) +SW(fit_nb2 <- stan_besag(y_pois ~ 1 + x, data = spatial_data, + prior_intercept = normal(0,1), prior = normal(0,1), prior_aux = normal(0,1), + family = neg_binomial_2(link = "log"), + W = adj, iter = ITER, chains = CHAINS)) +SW(fit_gamma <- stan_besag(y_gamma ~ 1 + x, data = spatial_data, family = Gamma(link = "log"), + prior_intercept = normal(0,1), prior = normal(0,1), prior_aux = normal(0,1), + W = adj, iter = ITER, chains = CHAINS)) + +# test QR +test_that("QR errors when number of predictors is <= 1", { + expect_error( + stan_besag(y_gauss ~ x, data = spatial_data, W = adj, family = gaussian(), seed = SEED, QR = TRUE), + "'QR' can only be specified when there are multiple predictors" + ) +}) + +test_that("QR works when number of x predictors is > 1", { + SW(fit_besag <- stan_besag(y_gauss ~ 1 + x + I(x^2), data = spatial_data, family = gaussian(), + W = adj, iter = ITER, chains = CHAINS, QR = TRUE)) + expect_stanreg(fit_besag) +}) + +test_that("stan_besag errors with algorithm = 'optimizing'", { + expect_error(stan_besag(y_gamma ~ 1 + x, data = spatial_data, family = gaussian(), + W = adj, iter = ITER, chains = CHAINS, algorithm = "optimizing"), + "'arg' should be one of “sampling”, “meanfield”, “fullrank”") +}) + +test_that("loo/waic for stan_besag works", { + loo(fit_gauss) + loo(fit_binom) + loo(fit_pois) + loo(fit_nb2) + loo(fit_gamma) +}) + +test_that("posterior_predict works for stan_besag", { + preds_gauss <- posterior_predict(fit_gauss) + preds_binom <- posterior_predict(fit_binom) + preds_pois <- posterior_predict(fit_pois) + preds_nb2 <- posterior_predict(fit_nb2) + preds_gamma <- posterior_predict(fit_gamma) +}) + +test_that("predict works for stan_besag", { + new_dat <- data.frame(x = rnorm(nrow(adj), 2, 1)) + predict(fit_gauss) + predict(fit_binom) + predict(fit_pois) + predict(fit_nb2) + predict(fit_gamma) + predict(fit_gauss, newdata = new_dat) + predict(fit_binom, newdata = new_dat) + predict(fit_pois, newdata = new_dat) + predict(fit_nb2, newdata = new_dat) + predict(fit_gamma, newdata = new_dat) +}) + +test_that("predict errors if nrow(newdata) < number of spatial regions", { + expect_error(predict(fit_gauss, newdata = data.frame(x = rnorm(10, 2, 1))), + "'newdata' is less than the number of spatial regions.") +}) diff --git a/tests/testthat/test_stan_bym.R b/tests/testthat/test_stan_bym.R new file mode 100644 index 000000000..0cd8b3639 --- /dev/null +++ b/tests/testthat/test_stan_bym.R @@ -0,0 +1,119 @@ +# 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. + +# tests can be run using devtools::test() or manually by loading testthat +# package and then running the code below possibly with options(mc.cores = 4). + +library(rstanarm) +SEED <- 12345 +set.seed(SEED) +ITER <- 25 +CHAINS <- 2 + +context("stan_bym") +# source(paste0("tests/testthat/",(file.path("helpers", "expect_stanreg.R")))) +# source(paste0("tests/testthat/",(file.path("helpers", "SW.R")))) +source(file.path("helpers", "expect_stanreg.R")) +source(file.path("helpers", "SW.R")) + +data("lattice", package = "rstanarm") + +# Convert a spatial polygon to an N-by-N weight matrix +sp2weightmatrix <- function(spatialpolygon) { + spdep::nb2mat(spdep::poly2nb(spatialpolygon, queen = TRUE), style = "B", zero.policy = TRUE) +} +adj <- sp2weightmatrix(grid_sim15) +spatial_data <- grid_sim15@data + +SW(fit_gauss <- stan_bym(y_gauss ~ 1 + x, data = spatial_data, family = gaussian(link = "identity"), + prior_intercept = normal(0,1), prior = normal(0,1), prior_aux = normal(0,1), + prior_unstructured = normal(0,1), prior_structured = normal(0,1), + W = adj, iter = ITER, chains = CHAINS)) +SW(fit_binom <- stan_bym(y_binom ~ 1 + x, trials = spatial_data$trials, data = spatial_data, + prior_intercept = normal(0,1), prior = normal(0,1), + prior_unstructured = normal(0,1), prior_structured = normal(0,1), + family = binomial(link = "logit"), + W = adj, iter = ITER, chains = CHAINS)) +SW(fit_pois <- stan_bym(y_pois ~ 1 + x, data = spatial_data, + prior_intercept = normal(0,1), prior = normal(0,1), + prior_unstructured = normal(0,1), prior_structured = normal(0,1), + family = poisson(link = "log"), + W = adj, iter = ITER, chains = CHAINS)) +SW(fit_nb2 <- stan_bym(y_pois ~ 1 + x, data = spatial_data, + prior_intercept = normal(0,1), prior = normal(0,1), prior_aux = normal(0,1), + prior_unstructured = normal(0,1), prior_structured = normal(0,1), + family = neg_binomial_2(link = "log"), + W = adj, iter = ITER, chains = CHAINS)) +SW(fit_gamma <- stan_bym(y_gamma ~ 1 + x, data = spatial_data, family = Gamma(link = "log"), + prior_intercept = normal(0,1), prior = normal(0,1), prior_aux = normal(0,1), + prior_unstructured = normal(0,1), prior_structured = normal(0,1), + W = adj, iter = ITER, chains = CHAINS)) + +# test QR +test_that("QR errors when number of predictors is <= 1", { + expect_error( + stan_bym(y_gauss ~ x, data = spatial_data, family = gaussian(), W = adj, seed = SEED, QR = TRUE), + "'QR' can only be specified when there are multiple predictors" + ) +}) + +test_that("QR works when number of x predictors is > 1", { + SW(fit_bym <- stan_bym(y_gauss ~ 1 + x + I(x^2), data = spatial_data, family = gaussian(), + W = adj, iter = ITER, chains = CHAINS, QR = TRUE)) + expect_stanreg(fit_bym) +}) + +test_that("stan_bym errors with algorithm = 'optimizing'", { + expect_error(stan_bym(y_gamma ~ 1 + x, data = spatial_data, family = gaussian(), + W = adj, iter = ITER, chains = CHAINS, algorithm = "optimizing"), + "'arg' should be one of “sampling”, “meanfield”, “fullrank”") +}) + +test_that("loo/waic for stan_bym works", { + loo(fit_gauss) + loo(fit_binom) + loo(fit_pois) + loo(fit_nb2) + loo(fit_gamma) +}) + +test_that("posterior_predict works for stan_bym", { + preds_gauss <- posterior_predict(fit_gauss) + preds_binom <- posterior_predict(fit_binom) + preds_pois <- posterior_predict(fit_pois) + preds_nb2 <- posterior_predict(fit_nb2) + preds_gamma <- posterior_predict(fit_gamma) +}) + +test_that("predict works for stan_bym", { + new_dat <- data.frame(x = rnorm(nrow(adj), 2, 1)) + predict(fit_gauss) + predict(fit_binom) + predict(fit_pois) + predict(fit_nb2) + predict(fit_gamma) + predict(fit_gauss, newdata = new_dat) + predict(fit_binom, newdata = new_dat) + predict(fit_pois, newdata = new_dat) + predict(fit_nb2, newdata = new_dat) + predict(fit_gamma, newdata = new_dat) +}) + +test_that("predict errors if nrow(newdata) < number of spatial regions", { + expect_error(predict(fit_gauss, newdata = data.frame(x = rnorm(10, 2, 1))), + "'newdata' is less than the number of spatial regions.") +}) diff --git a/tests/testthat/test_stan_bym2.R b/tests/testthat/test_stan_bym2.R new file mode 100644 index 000000000..a1fc3a902 --- /dev/null +++ b/tests/testthat/test_stan_bym2.R @@ -0,0 +1,119 @@ +# 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. + +# tests can be run using devtools::test() or manually by loading testthat +# package and then running the code below possibly with options(mc.cores = 4). + +library(rstanarm) +SEED <- 12345 +set.seed(SEED) +ITER <- 25 +CHAINS <- 2 + +context("stan_bym2") +# source(paste0("tests/testthat/",(file.path("helpers", "expect_stanreg.R")))) +# source(paste0("tests/testthat/",(file.path("helpers", "SW.R")))) +source(file.path("helpers", "expect_stanreg.R")) +source(file.path("helpers", "SW.R")) + +data("lattice", package = "rstanarm") + +# Convert a spatial polygon to an N-by-N weight matrix +sp2weightmatrix <- function(spatialpolygon) { + spdep::nb2mat(spdep::poly2nb(spatialpolygon, queen = TRUE), style = "B", zero.policy = TRUE) +} +adj <- sp2weightmatrix(grid_sim15) +spatial_data <- grid_sim15@data + +SW(fit_gauss <- stan_bym2(y_gauss ~ 1 + x, data = spatial_data, family = gaussian(link = "identity"), + prior_intercept = normal(0,1), prior = normal(0,1), prior_aux = normal(0,1), + prior_structured = normal(0,1), prior_mixing = beta(0.5,0.5), + W = adj, iter = ITER, chains = CHAINS)) +SW(fit_binom <- stan_bym2(y_binom ~ 1 + x, trials = spatial_data$trials, data = spatial_data, + prior_intercept = normal(0,1), prior = normal(0,1), + prior_structured = normal(0,1), prior_mixing = beta(0.5,0.5), + family = binomial(link = "logit"), + W = adj, iter = ITER, chains = CHAINS)) +SW(fit_pois <- stan_bym2(y_pois ~ 1 + x, data = spatial_data, + prior_intercept = normal(0,1), prior = normal(0,1), + prior_structured = normal(0,1), prior_mixing = beta(0.5,0.5), + family = poisson(link = "log"), + W = adj, iter = ITER, chains = CHAINS)) +SW(fit_nb2 <- stan_bym2(y_pois ~ 1 + x, data = spatial_data, + prior_intercept = normal(0,1), prior = normal(0,1), prior_aux = normal(0,1), + prior_structured = normal(0,1), prior_mixing = beta(0.5,0.5), + family = neg_binomial_2(link = "log"), + W = adj, iter = ITER, chains = CHAINS)) +SW(fit_gamma <- stan_bym2(y_gamma ~ 1 + x, data = spatial_data, family = Gamma(link = "log"), + prior_intercept = normal(0,1), prior = normal(0,1), prior_aux = normal(0,1), + prior_structured = normal(0,1), prior_mixing = beta(0.5,0.5), + W = adj, iter = ITER, chains = CHAINS)) + +# test QR +test_that("QR errors when number of predictors is <= 1", { + expect_error( + stan_bym2(y_gauss ~ x, data = spatial_data, family = gaussian(), W = adj, seed = SEED, QR = TRUE), + "'QR' can only be specified when there are multiple predictors" + ) +}) + +test_that("QR works when number of x predictors is > 1", { + SW(stan_bym <- stan_bym2(y_gauss ~ 1 + x + I(x^2), data = spatial_data, family = gaussian(), + W = adj, iter = ITER, chains = CHAINS, QR = TRUE)) + expect_stanreg(stan_bym) +}) + +test_that("stan_bym2 errors with algorithm = 'optimizing'", { + expect_error(stan_bym2(y_gamma ~ 1 + x, data = spatial_data, family = gaussian(), + W = adj, iter = ITER, chains = CHAINS, algorithm = "optimizing"), + "'arg' should be one of “sampling”, “meanfield”, “fullrank”") +}) + +test_that("loo/waic for stan_bym2 works", { + loo(fit_gauss) + loo(fit_binom) + loo(fit_pois) + loo(fit_nb2) + loo(fit_gamma) +}) + +test_that("posterior_predict works for stan_bym2", { + preds_gauss <- posterior_predict(fit_gauss) + preds_binom <- posterior_predict(fit_binom) + preds_pois <- posterior_predict(fit_pois) + preds_nb2 <- posterior_predict(fit_nb2) + preds_gamma <- posterior_predict(fit_gamma) +}) + +test_that("predict works for stan_bym2", { + new_dat <- data.frame(x = rnorm(nrow(adj), 2, 1)) + predict(fit_gauss) + predict(fit_binom) + predict(fit_pois) + predict(fit_nb2) + predict(fit_gamma) + predict(fit_gauss, newdata = new_dat) + predict(fit_binom, newdata = new_dat) + predict(fit_pois, newdata = new_dat) + predict(fit_nb2, newdata = new_dat) + predict(fit_gamma, newdata = new_dat) +}) + +test_that("predict errors if nrow(newdata) < number of spatial regions", { + expect_error(predict(fit_gauss, newdata = data.frame(x = rnorm(10, 2, 1))), + "'newdata' is less than the number of spatial regions.") +}) diff --git a/vignettes/spatial.Rmd b/vignettes/spatial.Rmd new file mode 100644 index 000000000..231afb5b3 --- /dev/null +++ b/vignettes/spatial.Rmd @@ -0,0 +1,257 @@ +--- +title: "CAR Spatial Models with rstanarm" +author: "Imad Ali, Jonah Gabry and Ben Goodrich" +date: "`r Sys.Date()`" +output: + html_vignette: + toc: yes +params: + EVAL: !r identical(Sys.getenv("NOT_CRAN"), "true") +--- + +```{r, child="children/SETTINGS-knitr.txt"} +``` +```{r, child="children/SETTINGS-gg.txt"} +``` +```{r, child="children/SETTINGS-rstan.txt"} +``` +```{r, child="children/SETTINGS-loo.txt"} +``` + +## Introduction + +This vignette explains how to account for spatial variation with conditional autoregressive (CAR) models when modeling discrete or continuous outcomes using the `stan_besag`, `stan_bym`, and `stan_bym2` functions in the __rstanarm__ package. The `stan_besag` function allows you to model space as an intrinsic conditional autoregressive model (ICAR), the `stan_bym` function allows you to fit the Besag, York, Mollie (1991) CAR model, and the `stan_bym2` is a variant of the BYM model where the priors are more interpretable (see Riebler et al (2016) for details). + +```{r, child="children/four_steps.txt"} +``` + +Steps 3 and 4 are covered in more depth by the vignette entitled ["How to Use +the __rstanarm__ Package"](rstanarm.html). + +## Likelihood + +The likelihoods supported in each CAR modeling function include Gaussian, Binomial, Poisson, Negative Binomial, and Gamma. For details on the link functions available and the contribution to the posterior of each of these likelihoods consult the appropriate [vignette](http://mc-stan.org/rstanarm/articles/index.html). + +The linear predictor takes the following form, +$$ +\boldsymbol{\eta} = \alpha + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\psi} +$$ +where $\alpha$ is the intercept, $\mathbf{X}$ is an $N$-by-$K$ matrix of predictors ($N$ being the number of observations and $K$ being the number of predictors), $\boldsymbol{\beta}$ is a $K$-dimensional vector of regression coefficients, and $\boldsymbol{\psi}$ is a $N$-dimensional vector representing the spatial effect. The construction of $\boldsymbol{\psi}$ depends on the model, which is discussed in the relevant sections below. + +Depending on the choice of likelihood there may or may not be an additional auxiliary parameter $\gamma$ in the model (e.g. in a Gaussian likelihood this would be the variation of the data). With all these components, for some probability density/mass function $f$, we can state the general form of the likelihood as, + +$$ +\mathcal{L}(\alpha, \boldsymbol{\beta}, \gamma | \mathbf{y}) = \prod_{i=1}^N f(y_i | \alpha, \boldsymbol{\beta}, \gamma ) +$$ + +## GMRF Hierarchical Component + +CAR models require that you define the spatial component as a Gaussian Markov Random Field (GMRF). The random vector $\boldsymbol{\phi}$ is a GMRF with respect to the graph $\mathcal{G} = (\mathcal{V} = \{1,\ldots,n\},\mathcal{E})$ with mean vector $\boldsymbol{\mu}$ and precision matrix $\mathbf{W}$ if its probability density function takes the precision form of the multivariate normal distribution, + +$$ +f(\boldsymbol{\phi} | \boldsymbol{\mu}) = (2\pi)^{-n/2} |\mathbf{W}|^{1/2} +\exp\bigg( -\frac{1}{2}(\boldsymbol{\phi-\mu})^{\top}\mathbf{W}(\boldsymbol{\phi-\mu}) \bigg) +$$ + +and $W_{i,j} \neq 0$ when $\{i,j\}\in\mathcal{E}$ for all $i \neq j$. Here, $\mathcal{V}$ refers to the verticies on the graph (i.e. the spatial units) and $\mathcal{E}$ refers to the edges on the graph (i.e. the spatial units that are neighbors). In other words, $\mathbf{W}$ is an $N$-by-$N$ matrix (with zeros on the diagonal) which describes spatial adjacency. + +Unfortunately there is no guarantee that $\mathbf{W}$ is positive definite so $\mathbf{Q} = \mbox{diag}(\mathbf{W1}) - W$ (which is guaranteed to be positive semi-definite) is used as the precision matrix. + +## Besag (ICAR) Spatial Prior + +The `stan_besag` modeling function fits the data to an ICAR model. This means that the spatial effect enters the linear predictor as + +$$ +\begin{align} +\boldsymbol{\psi} = \boldsymbol{\phi} \\ +f(\boldsymbol{\phi} | \boldsymbol{\mu}) &= (2\pi)^{-n/2} |\tau\mathbf{Q}|^{1/2} +\exp\bigg( -\frac{1}{2}\boldsymbol{\phi}^{\top}\tau\mathbf{Q}\boldsymbol{\phi} \bigg) +\end{align} +$$ + +where $\tau$ is a scalar that controls the overall spatial variation and has an appropriate hyperprior distribution. + +## BYM + +The ICAR model is limited in that it only accounts for spatial variation among the spatial units. Thus, the random variation is picked up by the spatial variation which may result in misleading parameter estimates and invalid inferences. The `stan_bym` model explains spatial variation as the sum of a structured (spatial) component $\boldsymbol{\phi}$ and an unstructured (random) component $\boldsymbol{\theta}$. Therefore, the spatial effect takes the following form, + +$$ +\begin{align} +\boldsymbol{\psi} &= \rho\boldsymbol{\phi} + \tau\boldsymbol{\theta} \\ +f(\boldsymbol{\phi} | \boldsymbol{\mu}) &= (2\pi)^{-n/2} |\mathbf{Q}|^{1/2} +\exp\bigg( -\frac{1}{2}\boldsymbol{\phi}^{\top}\mathbf{Q}\boldsymbol{\phi} \bigg) \\ +f(\theta_i) &= (2\pi)^{-1/2}\exp{\bigg( -\frac{\theta_i^2}{2} \bigg)} +\end{align} +$$ + +Note that the unstructured effect $\boldsymbol{\theta}$ is distributed standard normal. + +## BYM2 (Variant of the BYM Spatial Prior) + +The `stan_bym2` modeling function fits the data to a variant of the BYM model where the spatial effect enters as a convolution of the structured (spatial) effect and the unstructured (random) effect, + +$$ +\begin{align} +\boldsymbol{\psi} &= \tau(\boldsymbol{\theta}\sqrt{1-\rho} + \boldsymbol{\phi}\sqrt{\rho}) \\ +f(\boldsymbol{\phi} | \boldsymbol{\mu}) &= (2\pi)^{-n/2} |\mathbf{Q}|^{1/2} +\exp\bigg( -\frac{1}{2}\boldsymbol{\phi}^{\top}\mathbf{Q}\boldsymbol{\phi} \bigg) \\ +f(\theta_i) &= (2\pi)^{-1/2}\exp{\bigg( -\frac{\theta_i^2}{2} \bigg)} +\end{align} +$$ + +As in the BYM model $\boldsymbol{\theta}$ is distributed standard normal. However, now the parameter $\rho$ is on the unit interval and is interpreted as the proportion of spatial variation that is contributed to overall variation, and $\tau$ explains the overall (convolved) variation. + +Priors on $\rho$ should be chosen carefully as $\rho=0$ reduces to a model that not account for spatial variation and $\rho=1$ reduces to the ICAR model, which does not account for random variation among the spatial units. + +## Posterior + +Combining these components yields the following posteriors. For `stan_besag` we have, +$$ +f(\alpha, \boldsymbol{\beta},\boldsymbol{\phi}, \rho, \gamma | \mathbf{y},\mathbf{X}, \mathbf{Q}) \propto +\prod_{i=1}^N f(y_i | \gamma) \times +\prod_{k=1}^K f(\beta_k) \times +f(\boldsymbol{\phi} | \mathbf{Q}, \rho) \times +f(\rho) \times +f(\lambda) +$$ + + +and for both `stan_bym` and `stan_bym2` we have, +$$ +f(\alpha, \boldsymbol{\beta},\boldsymbol{\phi}, \boldsymbol{\theta}, \rho, \tau, \gamma | \mathbf{y},\mathbf{X}, \mathbf{Q}) \propto +\prod_{i=1}^N f(y_i | \gamma) \times +\prod_{k=1}^K f(\beta_k) \times +f(\boldsymbol{\phi} | \mathbf{Q}, \rho) \times +\prod_{i=1}^K f(\theta_i) \times +f(\rho) \times +f(\tau) \times +f(\lambda) +$$ + +## An Example Using Simulated Data on a Lattice + +As an example we use spatial units defined on a lattice. Below we plot a GMRF of 900 spatial units available in the rstanarm package. + +```{r load-lattice, fig.align='center', fig.height=8} +library(spdep) +# Load the preconstruced lattice/GMRF +data("lattice", package = "rstanarm") +grid_sim <- grid_sim30 +# plot the GMRF +var_range_gmrf <- seq(min(grid_sim@data$gmrf), max(grid_sim@data$gmrf), length = 50) +spplot(grid_sim, "gmrf", at = var_range_gmrf, main = expression(paste(phi, " (GMRF)")), + col.regions = colorRampPalette(c("#ef8a62", "#f7f7f7", "#67a9cf"))(50)) +``` + +Now we simulate predictors and a binomial outcome for the ICAR model. + +```{r simulate-data} +# Convert a spatial polygon to an N-by-N weight matrix +sp2weightmatrix <- function(spatialpolygon) { + spdep::nb2mat(spdep::poly2nb(spatialpolygon, queen = TRUE), style = "B", zero.policy = TRUE) +} +# Neighborhood matrix +W <- sp2weightmatrix(grid_sim) +# Simulate predictors and outcome +x <- rnorm(nrow(W), 2, 1) +z <- rnorm(nrow(W), -1, 1) +trials <- rep(10, nrow(W)) +spatial_data <- data.frame(x, z, phi = grid_sim@data$gmrf, trials) +eta <- binomial(link="logit")$linkinv(0.5 + 0.4 * x + 0.8 * z + spatial_data$phi) +spatial_data$y <- rbinom(nrow(W), trials, eta) +``` + +Plotting the outcome data on the lattice gives us a better understanding of its spatial variation. + +```{r plot-outcome} +grid_sim@data$y <- spatial_data$y +var_range_y <- seq(min(grid_sim@data$y), max(grid_sim@data$y) + 1, by = 1) +spplot(grid_sim, "y", at = var_range_y, main = expression(y), + col.regions = colorRampPalette(c("#ef8a62", "#f7f7f7", "#67a9cf"))(max(var_range_y) + 1)) +``` + +Now we can fit the model with `stan_besag` (we also fit a model that does not match that data generating process for model selection purposes). + +```{r fit-besag, results = "hide"} +library(rstanarm) +fit_besag <- stan_besag(y ~ 1 + x + z, data = spatial_data, W, + prior_intercept = normal(0,1), prior = normal(0,1), prior_rho = normal(0,1), + family = binomial(link="logit"), trials = spatial_data$trials, + chains = CHAINS, cores = CORES, seed = SEED, iter = ITER) +fit_besag_bad <- stan_besag(y ~ 1 + x + I(x^2), data = spatial_data, W, + prior_intercept = normal(0,1), prior = normal(0,1), prior_rho = normal(0,1), + family = binomial(link="logit"), trials = spatial_data$trials, + chains = CHAINS, cores = CORES, seed = SEED, iter = ITER) +print(fit_besag) +print(fit_besag_bad) +``` + +As a two-dimensional posterior predictive check we can plot the posterior predictions on the lattice using the `posterior_predict` function. +```{r ppcheck-2d, fig.align='center', fig.height=8} +grid_sim@data$y_pred <- colMeans(posterior_predict(fit_besag)) +var_range_y_pred <- seq(min(grid_sim@data$y_pred), max(grid_sim@data$y_pred) + 1, by = 1) +spplot(grid_sim, "y_pred", at = var_range_y_pred, main = expression(y[pred]), + col.regions = colorRampPalette(c("#ef8a62", "#f7f7f7", "#67a9cf"))(max(var_range_y_pred) + 1)) +``` + +Alternatively we can look at the conventional one-dimensional posterior predictive check with the `pp_check` function. + +```{r ppcheck-1d, fig.align='center', fig.height=8} +pp_check(fit_besag) +``` + +In order to compare the predictive performance between the models we need to use the [loo](http://mc-stan.org/loo) package. Looking at the results below we can confirm that the correct model outperforms the incorrect model in terms of predictive accuracy. + +```{r loo} +library(loo) +loo_correct <- loo(fit_besag) +loo_incorrect <- loo(fit_besag_bad) +compare(loo_correct, loo_incorrect) +``` + +## Smoothing the Spatial Random Walk + +In some cases modeling the GMRF spatial component with the precision matrix $\mathbf{Q}$ leads to rough spatial varation. This often occurs when dealing with spatial units on a fine lattice. Using the square of the precision matrix $\mathbf{Q}\mathbf{Q}$ allows us to smooth out the spatial variation across the spatial units using a spatial random walk of order 2. Below we use the `order = 2` argument to fit the model with smoothing of order 2. + +```{r smoothing, results="hide"} +fit_besag_smooth <- stan_besag(y ~ 1 + x + z, data = spatial_data, W, order = 2, + prior_intercept = normal(0,1), prior = normal(0,1), + prior_rho = normal(0,1), + family = binomial(link="logit"), trials = spatial_data$trials, + chains = CHAINS, cores = CORES, seed = SEED, iter = ITER) +``` + +Now we plot the results alongside the first fit with `order = 1`. + +```{r plot-smoothing, fig.align='center', fig.height=12} +grid_sim@data$y_pred_smooth <- colMeans(posterior_predict(fit_besag_smooth)) + +var_range_y_pred <- seq(min(grid_sim@data$y_pred_smooth), max(grid_sim@data$y_pred_smooth) + 1, by = 1) +gridExtra::grid.arrange(layout_matrix = matrix(1:2, nrow = 1), + spplot(grid_sim, "y_pred_smooth", at = var_range_y_pred, main = expression(y[pred_smooth]), + col.regions = colorRampPalette(c("#ef8a62", "#f7f7f7", "#67a9cf"))(max(var_range_y_pred) + 1)), + spplot(grid_sim, "y_pred", at = var_range_y_pred, main = expression(y[pred]), + col.regions = colorRampPalette(c("#ef8a62", "#f7f7f7", "#67a9cf"))(max(var_range_y_pred) + 1)) +) +``` + +The `stan_bym` and `stan_bym2` functions can be used in a similar way. + +## References + +Besag, J. (1974). Spatial Interaction and the Statistical Analysis of Lattice Systems. _Journal of the Royal Statistical Society._ Vol. 36, No. 2, p192-236. + +Besag, J., York, J. and Mollié, A. (1991). Bayesian image restoration, with two applications in spatial statistics. _Annals of the Institute of Statistical Mathematics._ Vol. 43, No. 01, p1-20. + +Riebler, A., Sorbye, S.H., Simpson, D. and Rue, H. (2016). An intuitive Bayesian spatial model for disease mapping that accounts for scaling. _arXiv_ preprint arXiv:1601.01180. + +Rue, H. and Held, L. (2005) _Gaussian Markov Random Fields: Theory And Applications (Monographs on Statistics and Applied Probability)._ Chapman & Hall/CRC, Boca Raton, FL. + +Rue, H. and Martino, S. (2007) Approximate Bayesian inference for hierarchical Gaussian Markov random fields models. _Journal of Statistical Planning and Inference._ Vol. 137, p3177-3192. + +Simpson, D., Rue, H., Martins, T.G., Riebler, A. and Sorbye, S.H. (2015). Penalising model component complexity: A principled, practical approach to constructing priors. _arXiv_ preprint arXiv:1403.4630. +