Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Updates for Regression and Other Stories #432

Merged
merged 27 commits into from
Jul 8, 2020
Merged

Updates for Regression and Other Stories #432

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
27 commits
Select commit Hold shift + click to select a range
8d43c7f
prior scale for intercept now 2.5 instead of 10
jgabry Apr 27, 2020
8bfe1d2
change autoscale usage
jgabry Apr 27, 2020
0973436
add default_prior_intercept, default_prior_coef
jgabry Apr 27, 2020
dd800ab
Update priors.R
jgabry Apr 27, 2020
76a3806
Update NAMESPACE
jgabry Apr 27, 2020
fa19d13
introduce posterior_epred
bgoodri May 3, 2020
7add53f
predict.R: remove error message for RAOS
bgoodri May 3, 2020
90d49bc
add (currently unused) family argument to default_prior_*
jgabry May 14, 2020
b469f1e
update posterior_linpred/epred doc
jgabry May 14, 2020
92d10a6
import posterior_epred generic from rstantools
jgabry May 14, 2020
fef4102
use posterior_epred instead of posterior_linpred in *_R2() functions
jgabry May 15, 2020
dfd7d33
use posterior_epred instead of linpred with transform
jgabry May 15, 2020
6ebfa90
logit and invlogit convenience functions
jgabry May 18, 2020
cc050a0
Update DESCRIPTION
jgabry Jun 3, 2020
64126fc
Merge branch 'master' into new-default-prior
jgabry Jun 10, 2020
c5b67a4
Update posterior_linpred.R
jgabry Jun 21, 2020
f97753c
always divide by sd(x) (unless x is constant)
jgabry Jun 26, 2020
eb5cbd4
Update priors.Rmd
jgabry Jun 26, 2020
4c3b8e2
Merge branch 'new-default-prior' of https://github.com/stan-dev/rstan…
jgabry Jun 26, 2020
dcf587a
more doc updates after updating default priors
jgabry Jun 26, 2020
5a43cad
Add note at top of priors vignette about changes
jgabry Jun 26, 2020
27e9988
remove autoscale=FALSE from vignettes
jgabry Jun 26, 2020
efc6d16
Update zzz.R
jgabry Jun 27, 2020
781a131
mention version number in vignette
jgabry Jun 29, 2020
3a412a9
averages -> expectations in posterior_epred doc
jgabry Jun 29, 2020
8708832
update message and doc for posterior_linpred(transform=TRUE)
jgabry Jul 1, 2020
5e8f452
Merge branch 'master' into new-default-prior
bgoodri Jul 8, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions NAMESPACE
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,6 +95,8 @@ export(bayes_R2)
export(cauchy)
export(compare_models)
export(decov)
export(default_prior_coef)
export(default_prior_intercept)
export(dirichlet)
export(exponential)
export(fixef)
Expand Down
40 changes: 27 additions & 13 deletions R/priors.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,7 @@
#' or equal to two, the mode of this Beta distribution does not exist
#' and an error will prompt the user to specify another choice for
#' \code{what}.
#' @param autoscale A logical scalar, defaulting to \code{TRUE}. If \code{TRUE}
#' @param autoscale If \code{TRUE}
#' then the scales of the priors on the intercept and regression coefficients
#' may be additionally modified internally by \pkg{rstanarm} in the following
#' cases. First, for Gaussian models only, the prior scales for the intercept,
Expand Down Expand Up @@ -104,12 +104,11 @@
#' approaches the normal distribution and if the degrees of freedom are one,
#' then the Student t distribution is the Cauchy distribution.
#'
#' If \code{scale} is not specified it will default to \eqn{10} for the
#' intercept and \eqn{2.5} for the other coefficients, unless the probit link
#' function is used, in which case these defaults are scaled by a factor of
#' \code{dnorm(0)/dlogis(0)}, which is roughly \eqn{1.6}.
#' If \code{scale} is not specified it will default to \eqn{2.5}, unless the
#' probit link function is used, in which case these defaults are scaled by a
#' factor of \code{dnorm(0)/dlogis(0)}, which is roughly \eqn{1.6}.
#'
#' If the \code{autoscale} argument is \code{TRUE} (the default), then the
#' If the \code{autoscale} argument is \code{TRUE}, then the
#' scales will be further adjusted as described above in the documentation of
#' the \code{autoscale} argument in the \strong{Arguments} section.
#' }
Expand Down Expand Up @@ -174,7 +173,7 @@
#'
#' It is also common in supervised learning to standardize the predictors
#' before training the model. We do not recommend doing so. Instead, it is
#' better to specify \code{autoscale = TRUE} (the default value), which
#' better to specify \code{autoscale = TRUE}, which
#' will adjust the scales of the priors according to the dispersion in the
#' variables. See the documentation of the \code{autoscale} argument above
#' and also the \code{\link{prior_summary}} page for more information.
Expand Down Expand Up @@ -436,22 +435,22 @@ NULL

#' @rdname priors
#' @export
normal <- function(location = 0, scale = NULL, autoscale = TRUE) {
normal <- function(location = 0, scale = NULL, autoscale = FALSE) {
validate_parameter_value(scale)
nlist(dist = "normal", df = NA, location, scale, autoscale)
}

#' @rdname priors
#' @export
student_t <- function(df = 1, location = 0, scale = NULL, autoscale = TRUE) {
student_t <- function(df = 1, location = 0, scale = NULL, autoscale = FALSE) {
validate_parameter_value(scale)
validate_parameter_value(df)
nlist(dist = "t", df, location, scale, autoscale)
}

#' @rdname priors
#' @export
cauchy <- function(location = 0, scale = NULL, autoscale = TRUE) {
cauchy <- function(location = 0, scale = NULL, autoscale = FALSE) {
student_t(df = 1, location = location, scale = scale, autoscale)
}

Expand Down Expand Up @@ -485,13 +484,13 @@ hs_plus <- function(df1 = 1, df2 = 1, global_df = 1, global_scale = 0.01,

#' @rdname priors
#' @export
laplace <- function(location = 0, scale = NULL, autoscale = TRUE) {
laplace <- function(location = 0, scale = NULL, autoscale = FALSE) {
nlist(dist = "laplace", df = NA, location, scale, autoscale)
}

#' @rdname priors
#' @export
lasso <- function(df = 1, location = 0, scale = NULL, autoscale = TRUE) {
lasso <- function(df = 1, location = 0, scale = NULL, autoscale = FALSE) {
nlist(dist = "lasso", df, location, scale, autoscale)
}

Expand All @@ -510,7 +509,7 @@ product_normal <- function(df = 2, location = 0, scale = 1) {
#' \code{1}. For the exponential distribution, the rate parameter is the
#' \emph{reciprocal} of the mean.
#'
exponential <- function(rate = 1, autoscale = TRUE) {
exponential <- function(rate = 1, autoscale = FALSE) {
stopifnot(length(rate) == 1)
validate_parameter_value(rate)
nlist(dist = "exponential",
Expand Down Expand Up @@ -562,7 +561,22 @@ R2 <- function(location = NULL, what = c("mode", "mean", "median", "log")) {
list(dist = "R2", location = location, what = what, df = 0, scale = 0)
}

#' @rdname priors
#' @export
default_prior_intercept = function() {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think default_prior_* should take a family argument and then stan_glm, etc. should default to default_prior(family). This will give us more flexibility to change things in the future.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok. What about for functions like stan_lmer(), stan_glm.nb(), stan_glmer.nb(), where there's no family argument?

That is, for stan_glmer() we can have

stan_glmer(..., prior_intercept = default_prior_intercept(family))

but for stan_lmer() and others should we just fill in the family like this?

stan_lmer(..., prior_intercept = default_prior_intercept(family=gaussian()))
stan_glmer.nb(..., prior_intercept = default_prior_intercept(family=neg_binomial_2()))
stan_glm.nb(..., prior_intercept = default_prior_intercept(family=neg_binomial_2()))

Or some other option?

out <- normal(0, 2.5, autoscale=TRUE)
out$location <- NULL # not determined yet
out$default <- TRUE
out
}

#' @rdname priors
#' @export
default_prior_coef = function() {
out <- normal(0, 2.5, autoscale=TRUE)
out$default <- TRUE
out
}


# internal ----------------------------------------------------------------
Expand Down
10 changes: 5 additions & 5 deletions R/stan_betareg.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -117,11 +117,11 @@ stan_betareg <-
y = TRUE,
x = FALSE,
...,
prior = normal(),
prior_intercept = normal(),
prior_z = normal(),
prior_intercept_z = normal(),
prior_phi = exponential(),
prior = normal(autoscale=TRUE),
prior_intercept = normal(autoscale=TRUE),
prior_z = normal(autoscale=TRUE),
prior_intercept_z = normal(autoscale=TRUE),
prior_phi = exponential(autoscale=TRUE),
prior_PD = FALSE,
algorithm = c("sampling", "optimizing", "meanfield", "fullrank"),
adapt_delta = NULL,
Expand Down
14 changes: 7 additions & 7 deletions R/stan_betareg.fit.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,11 +26,11 @@ stan_betareg.fit <-
offset = rep(0, NROW(x)),
link = c("logit", "probit", "cloglog", "cauchit", "log", "loglog"),
link.phi = NULL, ...,
prior = normal(),
prior_intercept = normal(),
prior_z = normal(),
prior_intercept_z = normal(),
prior_phi = exponential(),
prior = normal(autoscale=TRUE),
prior_intercept = normal(autoscale=TRUE),
prior_z = normal(autoscale=TRUE),
prior_intercept_z = normal(autoscale=TRUE),
prior_phi = exponential(autoscale=TRUE),
prior_PD = FALSE,
algorithm = c("sampling", "optimizing", "meanfield", "fullrank"),
adapt_delta = NULL,
Expand Down Expand Up @@ -101,7 +101,7 @@ stan_betareg.fit <-
assign(i, prior_stuff[[i]])

prior_intercept_stuff <- handle_glm_prior(prior_intercept, nvars = 1,
default_scale = 10, link = link,
default_scale = 2.5, link = link,
ok_dists = ok_intercept_dists)
names(prior_intercept_stuff) <- paste0(names(prior_intercept_stuff),
"_for_intercept")
Expand All @@ -115,7 +115,7 @@ stan_betareg.fit <-
assign(paste0(i,"_z"), prior_stuff_z[[i]])

prior_intercept_stuff_z <- handle_glm_prior(prior_intercept_z, nvars = 1,
link = link.phi, default_scale = 10,
link = link.phi, default_scale = 2.5,
ok_dists = ok_intercept_dists)
names(prior_intercept_stuff_z) <- paste0(names(prior_intercept_stuff_z),
"_for_intercept")
Expand Down
2 changes: 1 addition & 1 deletion R/stan_clogit.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,7 @@
#'
#' @importFrom lme4 findbars
stan_clogit <- function(formula, data, subset, na.action = NULL, ...,
strata, prior = normal(),
strata, prior = normal(autoscale=TRUE),
prior_covariance = decov(), prior_PD = FALSE,
algorithm = c("sampling", "optimizing",
"meanfield", "fullrank"),
Expand Down
6 changes: 3 additions & 3 deletions R/stan_gamm4.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,10 +127,10 @@ stan_gamm4 <-
knots = NULL,
drop.unused.levels = TRUE,
...,
prior = normal(),
prior_intercept = normal(),
prior = default_prior_coef(),
prior_intercept = default_prior_intercept(),
prior_smooth = exponential(autoscale = FALSE),
prior_aux = exponential(),
prior_aux = exponential(autoscale=TRUE),
prior_covariance = decov(),
prior_PD = FALSE,
algorithm = c("sampling", "meanfield", "fullrank"),
Expand Down
14 changes: 7 additions & 7 deletions R/stan_glm.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -138,7 +138,7 @@
#' counts ~ outcome + treatment,
#' data = count_data,
#' family = poisson(link="log"),
#' prior = normal(0, 2, autoscale = FALSE),
#' prior = normal(0, 2),
#' refresh = 0,
#' # for speed of example only
#' chains = 2, iter = 250
Expand Down Expand Up @@ -207,9 +207,9 @@ stan_glm <-
y = TRUE,
contrasts = NULL,
...,
prior = normal(),
prior_intercept = normal(),
prior_aux = exponential(),
prior = default_prior_coef(),
prior_intercept = default_prior_intercept(),
prior_aux = exponential(autoscale=TRUE),
prior_PD = FALSE,
algorithm = c("sampling", "optimizing", "meanfield", "fullrank"),
mean_PPD = algorithm != "optimizing",
Expand Down Expand Up @@ -319,9 +319,9 @@ stan_glm.nb <-
contrasts = NULL,
link = "log",
...,
prior = normal(),
prior_intercept = normal(),
prior_aux = exponential(),
prior = default_prior_coef(),
prior_intercept = default_prior_intercept(),
prior_aux = exponential(autoscale=TRUE),
prior_PD = FALSE,
algorithm = c("sampling", "optimizing", "meanfield", "fullrank"),
mean_PPD = algorithm != "optimizing",
Expand Down
16 changes: 12 additions & 4 deletions R/stan_glm.fit.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,9 +41,9 @@ stan_glm.fit <-
offset = rep(0, NROW(y)),
family = gaussian(),
...,
prior = normal(),
prior_intercept = normal(),
prior_aux = exponential(),
prior = default_prior_coef(),
prior_intercept = default_prior_intercept(),
prior_aux = exponential(autoscale=TRUE),
prior_smooth = exponential(autoscale = FALSE),
prior_ops = NULL,
group = list(),
Expand Down Expand Up @@ -132,6 +132,14 @@ stan_glm.fit <-
ok_aux_dists <- c(ok_dists[1:3], exponential = "exponential")

# prior distributions
if (isTRUE(prior_intercept$default)) {
m_y <- 0
if (family$family == "gaussian" && family$link == "identity") {
if (!is.null(y)) m_y <- mean(y) # y can be NULL if prior_PD=TRUE
}
prior_intercept$location <- m_y
}

prior_stuff <- handle_glm_prior(
prior,
nvars,
Expand All @@ -147,7 +155,7 @@ stan_glm.fit <-
prior_intercept_stuff <- handle_glm_prior(
prior_intercept,
nvars = 1,
default_scale = 10,
default_scale = 2.5,
link = family$link,
ok_dists = ok_intercept_dists
)
Expand Down
18 changes: 9 additions & 9 deletions R/stan_glmer.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,9 +103,9 @@ stan_glmer <-
offset,
contrasts = NULL,
...,
prior = normal(),
prior_intercept = normal(),
prior_aux = exponential(),
prior = default_prior_coef(),
prior_intercept = default_prior_intercept(),
prior_aux = exponential(autoscale=TRUE),
prior_covariance = decov(),
prior_PD = FALSE,
algorithm = c("sampling", "meanfield", "fullrank"),
Expand Down Expand Up @@ -202,9 +202,9 @@ stan_lmer <-
offset,
contrasts = NULL,
...,
prior = normal(),
prior_intercept = normal(),
prior_aux = exponential(),
prior = default_prior_coef(),
prior_intercept = default_prior_intercept(),
prior_aux = exponential(autoscale=TRUE),
prior_covariance = decov(),
prior_PD = FALSE,
algorithm = c("sampling", "meanfield", "fullrank"),
Expand Down Expand Up @@ -244,9 +244,9 @@ stan_glmer.nb <-
contrasts = NULL,
link = "log",
...,
prior = normal(),
prior_intercept = normal(),
prior_aux = exponential(),
prior = default_prior_coef(),
prior_intercept = default_prior_intercept(),
prior_aux = exponential(autoscale=TRUE),
prior_covariance = decov(),
prior_PD = FALSE,
algorithm = c("sampling", "meanfield", "fullrank"),
Expand Down
8 changes: 4 additions & 4 deletions R/stan_jm.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -523,10 +523,10 @@ stan_jm <- function(formulaLong, dataLong, formulaEvent, dataEvent, time_var,
lag_assoc = 0, grp_assoc, epsilon = 1E-5,
basehaz = c("bs", "weibull", "piecewise"), basehaz_ops,
qnodes = 15, init = "prefit", weights,
priorLong = normal(), priorLong_intercept = normal(),
priorLong_aux = cauchy(0, 5), priorEvent = normal(),
priorEvent_intercept = normal(), priorEvent_aux = cauchy(),
priorEvent_assoc = normal(), prior_covariance = lkj(),
priorLong = normal(autoscale=TRUE), priorLong_intercept = normal(autoscale=TRUE),
priorLong_aux = cauchy(0, 5, autoscale=TRUE), priorEvent = normal(autoscale=TRUE),
priorEvent_intercept = normal(autoscale=TRUE), priorEvent_aux = cauchy(autoscale=TRUE),
priorEvent_assoc = normal(autoscale=TRUE), prior_covariance = lkj(autoscale=TRUE),
prior_PD = FALSE, algorithm = c("sampling", "meanfield", "fullrank"),
adapt_delta = NULL, max_treedepth = 10L, QR = FALSE,
sparse = FALSE, ...) {
Expand Down
8 changes: 4 additions & 4 deletions R/stan_jm.fit.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,10 +27,10 @@ stan_jm.fit <- function(formulaLong = NULL, dataLong = NULL, formulaEvent = NULL
assoc = "etavalue", lag_assoc = 0, grp_assoc,
epsilon = 1E-5, basehaz = c("bs", "weibull", "piecewise"),
basehaz_ops, qnodes = 15, init = "prefit", weights,
priorLong = normal(), priorLong_intercept = normal(),
priorLong_aux = cauchy(0, 5), priorEvent = normal(),
priorEvent_intercept = normal(), priorEvent_aux = cauchy(),
priorEvent_assoc = normal(), prior_covariance = lkj(), prior_PD = FALSE,
priorLong = normal(autoscale=TRUE), priorLong_intercept = normal(autoscale=TRUE),
priorLong_aux = cauchy(0, 5, autoscale=TRUE), priorEvent = normal(autoscale=TRUE),
priorEvent_intercept = normal(autoscale=TRUE), priorEvent_aux = cauchy(autoscale=TRUE),
priorEvent_assoc = normal(autoscale=TRUE), prior_covariance = lkj(autoscale=TRUE), prior_PD = FALSE,
algorithm = c("sampling", "meanfield", "fullrank"),
adapt_delta = NULL, max_treedepth = 10L,
QR = FALSE, sparse = FALSE, ...) {
Expand Down
6 changes: 3 additions & 3 deletions R/stan_mvmer.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -129,9 +129,9 @@
#' }
#'
stan_mvmer <- function(formula, data, family = gaussian, weights,
prior = normal(), prior_intercept = normal(),
prior_aux = cauchy(0, 5),
prior_covariance = lkj(), prior_PD = FALSE,
prior = normal(autoscale=TRUE), prior_intercept = normal(autoscale=TRUE),
prior_aux = cauchy(0, 5, autoscale=TRUE),
prior_covariance = lkj(autoscale=TRUE), prior_PD = FALSE,
algorithm = c("sampling", "meanfield", "fullrank"),
adapt_delta = NULL, max_treedepth = 10L,
init = "random", QR = FALSE, sparse = FALSE, ...) {
Expand Down
4 changes: 2 additions & 2 deletions R/stan_nlmer.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,8 +95,8 @@ stan_nlmer <-
offset,
contrasts = NULL,
...,
prior = normal(),
prior_aux = exponential(),
prior = normal(autoscale=TRUE),
prior_aux = exponential(autoscale=TRUE),
prior_covariance = decov(),
prior_PD = FALSE,
algorithm = c("sampling", "meanfield", "fullrank"),
Expand Down
Loading