diff --git a/DESCRIPTION b/DESCRIPTION index a35cdab7e..fe741e520 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Package: brms Encoding: UTF-8 Type: Package Title: Bayesian Regression Models using Stan -Version: 1.8.0.1 +Version: 1.8.0.2 Date: 2017-07-19 Authors@R: person("Paul-Christian", "Bürkner", email = "paul.buerkner@gmail.com", role = c("aut", "cre")) diff --git a/NAMESPACE b/NAMESPACE index eb41f593f..cee59af9f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,7 @@ # Generated by roxygen2: do not edit by hand +S3method("+",brmsformula) +S3method("+",brmsprior) S3method(LOO,brmsfit) S3method(VarCorr,brmsfit) S3method(WAIC,brmsfit) @@ -11,8 +13,6 @@ S3method(as.data.frame,brmsVarCorr) S3method(as.data.frame,brmsfit) S3method(as.matrix,brmsfit) S3method(as.mcmc,brmsfit) -S3method(auxpar_family,default) -S3method(auxpar_family,mixfamily) S3method(bayes_R2,brmsfit) S3method(bayes_factor,brmsfit) S3method(bridge_sampler,brmsfit) @@ -25,11 +25,15 @@ S3method(coef,brmsfit) S3method(control_params,brmsfit) S3method(data_effects,btl) S3method(data_effects,btnl) +S3method(dpar_family,default) +S3method(dpar_family,mixfamily) S3method(expose_functions,brmsfit) S3method(extract_draws,brmsfit) S3method(extract_draws,btl) S3method(extract_draws,btnl) S3method(family,brmsfit) +S3method(family_names,brmsformula) +S3method(family_names,brmsterms) S3method(family_names,default) S3method(family_names,family) S3method(family_names,mixfamily) @@ -56,6 +60,9 @@ S3method(loo,brmsfit) S3method(loo_linpred,brmsfit) S3method(loo_predict,brmsfit) S3method(loo_predictive_interval,brmsfit) +S3method(make_Jmo_list,brmsterms) +S3method(make_Jmo_list,btl) +S3method(make_Jmo_list,btnl) S3method(make_gp_list,brmsterms) S3method(make_gp_list,btl) S3method(make_gp_list,btnl) @@ -122,8 +129,8 @@ S3method(summary,family) S3method(summary,mixfamily) S3method(update,brmsfit) S3method(update,brmsformula) -S3method(valid_auxpars,default) -S3method(valid_auxpars,mixfamily) +S3method(valid_dpars,default) +S3method(valid_dpars,mixfamily) S3method(vcov,brmsfit) S3method(waic,brmsfit) export("add_ic<-") @@ -211,6 +218,7 @@ export(kfold) export(lasso) export(launch_shiny) export(launch_shinystan) +export(lf) export(log_lik) export(log_posterior) export(logit_scaled) @@ -233,6 +241,7 @@ export(monotonic) export(neff_ratio) export(negbinomial) export(ngrps) +export(nlf) export(nsamples) export(nuts_params) export(parnames) @@ -284,6 +293,7 @@ export(rskew_normal) export(rstudent_t) export(rvon_mises) export(rwiener) +export(set_nl) export(set_prior) export(skew_normal) export(sratio) diff --git a/R/brm.R b/R/brm.R index ac82c6e07..9a18d2d93 100644 --- a/R/brm.R +++ b/R/brm.R @@ -462,7 +462,7 @@ brm <- function(formula, data, family = gaussian(), prior = NULL, if (future) { require_package("future") if (cores > 1L) { - warning("Argument 'cores' is ignored when using 'future'.") + warning2("Argument 'cores' is ignored when using 'future'.") } args$chains <- 1L futures <- fits <- vector("list", chains) diff --git a/R/brmsfit-helpers.R b/R/brmsfit-helpers.R index c1d046e5c..1515ac055 100644 --- a/R/brmsfit-helpers.R +++ b/R/brmsfit-helpers.R @@ -102,7 +102,7 @@ restructure <- function(x, rstr_summary = FALSE) { } } if (version <= "0.10.0.9000") { - if (length(bterms$auxpars$mu$nlpars)) { + if (length(bterms$dpars$mu$nlpars)) { # nlpar and group have changed positions change <- change_old_re(x$ranef, pars = parnames(x), dims = x$fit@sim$dims_oi) @@ -137,6 +137,9 @@ restructure <- function(x, rstr_summary = FALSE) { any(grepl("^Xme_", parnames(x))) } } + if (version <= "1.8.0.1") { + x$prior[, c("resp", "dpar")] <- "" + } stan_env <- attributes(x$fit)$.MISC if (rstr_summary && exists("summary", stan_env)) { stan_summary <- get("summary", stan_env) @@ -238,7 +241,7 @@ prepare_conditions <- function(x, conditions = NULL, effects = NULL, lapply(get_effect(bterms, "offset"), rhs), re$form, lapply(re$gcall, "[[", "weightvars"), bterms$adforms[c("se", "disp", "trials", "cat")], - bterms$auxpars$mu$covars + bterms$dpars$mu$covars ) req_vars <- unique(ulapply(req_vars, all.vars)) req_vars <- setdiff(req_vars, rsv_vars) @@ -662,8 +665,8 @@ get_cov_matrix_ident <- function(sigma, nrows, se2 = 0) { mat } -get_auxpar <- function(x, i = NULL) { - # get samples of an auxiliary parameter +get_dpar <- function(x, i = NULL) { + # get samples of an distributional parameter # Args: # x: object to extract postarior samples from # i: the current observation number @@ -673,7 +676,7 @@ get_auxpar <- function(x, i = NULL) { family <- x[["f"]] x <- get_eta(x, i = i) if (!nzchar(family$family)) { - # apply links for auxiliary parameters only + # apply links for distributional parameters only # the main family link is applied later on x <- ilink(x, family$link) } @@ -692,21 +695,21 @@ get_auxpar <- function(x, i = NULL) { get_sigma <- function(x, data, i = NULL, dim = NULL) { # get the residual standard devation of linear models # Args: - # see get_auxpar + # see get_dpar # dim: target dimension of output matrices (used in fitted) stopifnot(is.atomic(x) || is.list(x)) out <- get_se(data = data, i = i, dim = dim) if (!is.null(x)) { - out <- sqrt(out^2 + get_auxpar(x, i = i)^2) + out <- sqrt(out^2 + get_dpar(x, i = i)^2) } mult_disp(out, data = data, i = i, dim = dim) } get_shape <- function(x, data, i = NULL, dim = NULL) { # get the shape parameter of gamma, weibull and negbinomial models - # Args: see get_auxpar + # Args: see get_dpar stopifnot(is.atomic(x) || is.list(x)) - x <- get_auxpar(x, i = i) + x <- get_dpar(x, i = i) mult_disp(x, data = data, i = i, dim = dim) } @@ -715,14 +718,14 @@ get_zi_hu <- function(draws, i = NULL, par = c("zi", "hu")) { # also works with deprecated models fitted with brms < 1.0.0 # which were using multivariate syntax # Args: - # see get_auxpar + # see get_dpar # par: parameter to extract; either 'zi' or 'hu' par <- match.arg(par) if (!is.null(draws$data$N_trait)) { j <- if (!is.null(i)) i else seq_len(draws$data$N_trait) out <- ilink(get_eta(draws$mu, j + draws$data$N_trait), "logit") } else { - out <- get_auxpar(draws[[par]], i = i) + out <- get_dpar(draws[[par]], i = i) } out } @@ -736,7 +739,7 @@ get_theta <- function(draws, i = NULL) { families <- family_names(draws$f) theta <- vector("list", length(families)) for (j in seq_along(families)) { - theta[[j]] <- get_auxpar(draws[[paste0("theta", j)]], i = i) + theta[[j]] <- get_dpar(draws[[paste0("theta", j)]], i = i) } theta <- do.call(abind, c(theta, along = 3)) for (n in seq_len(dim(theta)[2])) { @@ -752,10 +755,10 @@ get_theta <- function(draws, i = NULL) { get_disc <- function(draws, i = NULL, ncat = NULL) { # convenience function to extract discrimination parameters # Args: - # see get_auxpar + # see get_dpar # ncat: number of response categories if (!is.null(draws[["disc"]])) { - disc <- get_auxpar(draws[["disc"]], i) + disc <- get_dpar(draws[["disc"]], i) if (!is.null(dim(disc))) { stopifnot(is.numeric(ncat)) disc <- array(disc, dim = c(dim(disc), ncat - 1)) @@ -768,7 +771,7 @@ get_disc <- function(draws, i = NULL, ncat = NULL) { get_se <- function(data, i = NULL, dim = NULL) { # extract user-defined standard errors - # Args: see get_auxpar + # Args: see get_dpar se <- data[["se"]] if (!is.null(se)) { if (!is.null(i)) { @@ -785,7 +788,7 @@ get_se <- function(data, i = NULL, dim = NULL) { mult_disp <- function(x, data, i = NULL, dim = NULL) { # multiply existing samples by 'disp' data - # Args: see get_auxpar + # Args: see get_dpar if (!is.null(data$disp)) { if (!is.null(i)) { x <- x * data$disp[i] @@ -884,7 +887,7 @@ fixef_pars <- function() { default_plot_pars <- function() { # list all parameter classes to be included in plots by default c(fixef_pars(), "^sd_", "^cor_", "^sigma_", "^rescor_", - paste0("^", auxpars(), "[[:digit:]]*$"), "^delta$", + paste0("^", dpars(), "[[:digit:]]*$"), "^delta$", "^theta", "^ar", "^ma", "^arr", "^lagsar", "^errorsar", "^car", "^sdcar", "^sigmaLL", "^sds_", "^sdgp_", "^lscale_") } diff --git a/R/brmsfit-methods.R b/R/brmsfit-methods.R index 6c7d18a56..085ee9eed 100644 --- a/R/brmsfit-methods.R +++ b/R/brmsfit-methods.R @@ -176,9 +176,8 @@ ranef.brmsfit <- function(object, summary = TRUE, robust = FALSE, groups <- unique(ranef$group) out <- named_list(groups) for (g in groups) { - take <- ranef$group == g - nlpars_usc <- usc(ranef[take, "nlpar"], "suffix") - coefs <- paste0(nlpars_usc, ranef[take, "coef"]) + r <- subset2(ranef, group = g) + coefs <- paste0(usc(combine_prefix(r), "suffix"), r$coef) levels <- attr(ranef, "levels")[[g]] rpars <- pars[grepl(paste0("^r_", g, "(__.+\\[|\\[)"), pars)] out[[g]] <- as.matrix(object, rpars, exact_match = TRUE) @@ -385,8 +384,8 @@ VarCorr.brmsfit <- function(x, sigma = 1, summary = TRUE, robust = FALSE, if (nrow(x$ranef)) { get_names <- function(group) { # get names of group-level parameters - r <- x$ranef[x$ranef$group == group, ] - rnames <- paste0(usc(r$nlpar, "suffix"), r$coef) + r <- subset2(x$ranef, group = group) + rnames <- paste0(usc(combine_prefix(r), "suffix"), r$coef) cor_type <- paste0("cor_", group) sd_pars <- paste0("sd_", group, "__", rnames) cor_pars <- get_cornames(rnames, type = cor_type, brackets = FALSE) @@ -430,7 +429,7 @@ posterior_samples.brmsfit <- function(x, pars = NA, parameters = NA, ...) { pars <- use_alias(pars, parameters, default = NA) add_chain <- use_alias(add_chain, add_chains, default = FALSE) - if (all(c(as.matrix, as.array))) { + if (as.matrix && as.array) { stop2("Cannot use 'as.matrix' and 'as.array' at the same time.") } if (add_chain && as.array) { @@ -518,22 +517,27 @@ as.mcmc.brmsfit <- function(x, pars = NA, exact_match = FALSE, combine_chains = FALSE, inc_warmup = FALSE, ...) { contains_samples(x) - pars <- extract_pars(pars, all_pars = parnames(x), - exact_match = exact_match, ...) + pars <- extract_pars( + pars, all_pars = parnames(x), + exact_match = exact_match, ... + ) if (combine_chains) { if (inc_warmup) { stop2("Cannot include warmup samples when 'combine_chains' is TRUE.") } out <- as.matrix(x$fit, pars) - mcpar <- c(x$fit@sim$warmup * x$fit@sim$chain + 1, - x$fit@sim$iter * x$fit@sim$chain, x$fit@sim$thin) + mcpar <- c( + x$fit@sim$warmup * x$fit@sim$chain + 1, + x$fit@sim$iter * x$fit@sim$chain, x$fit@sim$thin + ) attr(out, "mcpar") <- mcpar class(out) <- "mcmc" } else { - ps <- extract(x$fit, pars, permuted = FALSE, - inc_warmup = inc_warmup) - mcpar <- c(if (inc_warmup) 1 else x$fit@sim$warmup + 1, - x$fit@sim$iter, x$fit@sim$thin) + ps <- extract(x$fit, pars, permuted = FALSE, inc_warmup = inc_warmup) + mcpar <- c( + if (inc_warmup) 1 else x$fit@sim$warmup + 1, + x$fit@sim$iter, x$fit@sim$thin + ) out <- vector("list", length = dim(ps)[2]) for (i in seq_along(out)) { out[[i]] <- ps[, i, ] @@ -754,7 +758,7 @@ summary.brmsfit <- function(object, waic = FALSE, loo = FALSE, R2 = FALSE, rownames(out$fixed) <- gsub(fixef_pars(), "", fe_pars) # summary of family specific parameters - spec_pars <- c(auxpars(), "delta", "theta", "rescor") + spec_pars <- c(dpars(), "delta", "theta", "rescor") spec_pars <- paste0("^(", paste0(spec_pars, collapse = "|"), ")") spec_pars <- pars[grepl(spec_pars, pars)] spec_pars <- setdiff(spec_pars, "sigmaLL") @@ -775,8 +779,7 @@ summary.brmsfit <- function(object, waic = FALSE, loo = FALSE, R2 = FALSE, # summary of group-level effects for (g in out$group) { r <- object$ranef[object$ranef$group == g, ] - nlpar_usc <- usc(r$nlpar, "suffix") - rnames <- paste0(nlpar_usc, r$coef) + rnames <- paste0(usc(combine_prefix(r), "suffix"), r$coef) sd_pars <- paste0("sd_", g, "__", rnames) sd_names <- paste0("sd", "(", rnames ,")") # construct correlation names @@ -934,7 +937,7 @@ launch_shinystan.brmsfit <- function( object, rstudio = getOption("shinystan.rstudio"), ... ) { contains_samples(object) - shinystan::launch_shinystan(object$fit, rstudio = rstudio, ...) + launch_shinystan(object$fit, rstudio = rstudio, ...) } #' Trace and Density Plots for MCMC Samples @@ -1004,8 +1007,9 @@ plot.brmsfit <- function(x, pars = NA, parameters = NA, pars <- default_plot_pars() exact_match <- FALSE } - samples <- posterior_samples(x, pars = pars, add_chain = TRUE, - exact_match = exact_match) + samples <- posterior_samples( + x, pars = pars, add_chain = TRUE, exact_match = exact_match + ) pars <- names(samples)[!names(samples) %in% c("chain", "iter")] if (!length(pars)) { stop2("No valid parameters selected.") @@ -1021,8 +1025,9 @@ plot.brmsfit <- function(x, pars = NA, parameters = NA, for (i in seq_len(n_plots)) { sub_pars <- pars[((i - 1) * N + 1):min(i * N, length(pars))] sub_samples <- samples[, c(sub_pars, "chain"), drop = FALSE] - plots[[i]] <- bayesplot::mcmc_combo(sub_samples, combo = combo, - gg_theme = theme, ...) + plots[[i]] <- bayesplot::mcmc_combo( + sub_samples, combo = combo, gg_theme = theme, ... + ) if (plot) { plot(plots[[i]], newpage = newpage || i > 1) if (i == 1) { @@ -1345,9 +1350,11 @@ marginal_effects.brmsfit <- function(x, effects = NULL, conditions = NULL, } else if (is_categorical(x$family)) { stop2("Marginal plots are not yet implemented for categorical models.") } else if (is_ordinal(x$family)) { - warning2("Predictions are treated as continuous variables ", - "in marginal plots, \nwhich is likely an invalid ", - "assumption for family ", x$family$family, ".") + warning2( + "Predictions are treated as continuous variables ", + "in marginal_effects, which is likely an invalid ", + "assumption for family ", x$family$family, "." + ) } if (!is.null(transform) && method != "predict") { stop2("'transform' is only allowed when 'method' is set to 'predict'.") @@ -1379,15 +1386,19 @@ marginal_effects.brmsfit <- function(x, effects = NULL, conditions = NULL, if (sum(matches) > 0 && sum(matches > 0) < length(effects)) { invalid <- effects[setdiff(seq_along(effects), sort(matches))] invalid <- ulapply(invalid, paste, collapse = ":") - warning2("Some specified effects are invalid for this model: ", - collapse_comma(invalid), "\nValid effects are ", - "(combinations of): ", collapse_comma(ae_coll)) + warning2( + "Some specified effects are invalid for this model: ", + collapse_comma(invalid), "\nValid effects are ", + "(combinations of): ", collapse_comma(ae_coll) + ) } effects <- unique(effects[sort(matches)]) if (!length(effects)) { - stop2("All specified effects are invalid for this model.\n", - "Valid effects are (combinations of): ", - collapse_comma(ae_coll)) + stop2( + "All specified effects are invalid for this model.\n", + "Valid effects are (combinations of): ", + collapse_comma(ae_coll) + ) } } if (length(probs) != 2L) { @@ -1514,12 +1525,12 @@ marginal_smooths.brmsfit <- function(x, smooths = NULL, lee <- list() if (length(bterms$response) > 1L) { for (r in bterms$response) { - lee <- c(lee, setNames(bterms$auxpars["mu"], r)) + lee <- c(lee, setNames(bterms$dpars["mu"], r)) } - bterms$auxpars[["mu"]] <- NULL + bterms$dpars[["mu"]] <- NULL } - for (ap in names(bterms$auxpars)) { - bt <- bterms$auxpars[ap] + for (dp in names(bterms$dpars)) { + bt <- bterms$dpars[dp] if (is.btnl(bt[[1]])) { lee <- c(lee, bt[[1]]$nlpars) } else { @@ -1781,10 +1792,10 @@ predict.brmsfit <- function(object, newdata = NULL, re_formula = NULL, if (is.list(draws$mu[["mv"]])) { draws$mu <- get_eta(draws$mu) } - auxpars <- intersect(valid_auxpars(family(object)), names(draws)) - for (ap in auxpars) { - if (is.list(draws[[ap]])) { - draws[[ap]] <- get_auxpar(draws[[ap]]) + dpars <- intersect(valid_dpars(family(object)), names(draws)) + for (dp in dpars) { + if (is.list(draws[[dp]])) { + draws[[dp]] <- get_dpar(draws[[dp]]) } } # see predict.R @@ -1868,7 +1879,7 @@ posterior_predict.brmsfit <- function(object, newdata = NULL, re_formula = NULL, #' If \code{"response"} results are returned on the scale #' of the response variable. If \code{"linear"} #' fitted values are returned on the scale of the linear predictor. -#' @param auxpar Optional name of a predicted auxiliary parameter. +#' @param dpar Optional name of a predicted distributional parameter. #' If specified, fitted values of this parameters are returned. #' #' @return Fitted values extracted from \code{object}. @@ -1906,10 +1917,12 @@ fitted.brmsfit <- function(object, newdata = NULL, re_formula = NULL, allow_new_levels = FALSE, sample_new_levels = "uncertainty", new_objects = list(), incl_autocor = TRUE, - auxpar = NULL, subset = NULL, + dpar = NULL, subset = NULL, nsamples = NULL, sort = FALSE, nug = NULL, summary = TRUE, robust = FALSE, probs = c(0.025, 0.975), ...) { + dots <- list(...) + dpar <- use_alias(dpar, dots[["auxpar"]]) scale <- match.arg(scale) contains_samples(object) object <- restructure(object) @@ -1918,14 +1931,14 @@ fitted.brmsfit <- function(object, newdata = NULL, re_formula = NULL, sample_new_levels, new_objects, subset, nsamples, nug ) draws <- do.call(extract_draws, draws_args) - auxpars <- intersect(valid_auxpars(family(object)), names(draws)) - if (!length(auxpar)) { + dpars <- intersect(valid_dpars(family(object)), names(draws)) + if (!length(dpar)) { if (is.list(draws[["mu"]][["mv"]])) { draws$mu <- get_eta(draws$mu) } - for (ap in auxpars) { - if (is.list(draws[[ap]])) { - draws[[ap]] <- get_auxpar(draws[[ap]]) + for (dp in dpars) { + if (is.list(draws[[dp]])) { + draws[[dp]] <- get_dpar(draws[[dp]]) } } if (grepl("_mv$", draws$f$family) && length(dim(draws$mu)) == 3L) { @@ -1940,23 +1953,23 @@ fitted.brmsfit <- function(object, newdata = NULL, re_formula = NULL, draws$mu <- fitted_fun(draws) } } else { - auxpars <- auxpars[auxpar_class(auxpars) != "mu"] - if (length(auxpar) != 1L || !auxpar %in% auxpars) { - stop2("Invalid argument 'auxpar'. Valid auxiliary ", - "parameters are: ", collapse_comma(auxpars)) + dpars <- dpars[dpar_class(dpars) != "mu"] + if (length(dpar) != 1L || !dpar %in% dpars) { + stop2("Invalid argument 'dpar'. Valid distributional ", + "parameters are: ", collapse_comma(dpars)) } - if (!isTRUE(attr(draws[[auxpar]], "predicted"))) { - stop2("Auxiliary parameter '", auxpar, "' was not predicted.") + if (!isTRUE(attr(draws[[dpar]], "predicted"))) { + stop2("Distributional parameter '", dpar, "' was not predicted.") } - if (scale == "linear" && is.list(draws[[auxpar]])) { - draws[[auxpar]]$f$link <- "identity" + if (scale == "linear" && is.list(draws[[dpar]])) { + draws[[dpar]]$f$link <- "identity" } - if (auxpar_class(auxpar) == "theta" && scale == "response") { - ap_id <- as.numeric(auxpar_id(auxpar)) + if (dpar_class(dpar) == "theta" && scale == "response") { + ap_id <- as.numeric(dpar_id(dpar)) draws$mu <- get_theta(draws)[, , ap_id, drop = FALSE] dim(draws$mu) <- dim(draws$mu)[c(1, 2)] } else { - draws$mu <- get_auxpar(draws[[auxpar]]) + draws$mu <- get_dpar(draws[[dpar]]) } } if (is.null(dim(draws$mu))) { @@ -2252,8 +2265,10 @@ update.brmsfit <- function(object, formula., newdata = NULL, } } - arg_names <- c("prior", "autocor", "nonlinear", "threshold", - "cov_ranef", "sparse", "sample_prior") + arg_names <- c( + "prior", "autocor", "nonlinear", "threshold", + "cov_ranef", "sparse", "sample_prior" + ) new_args <- intersect(arg_names, names(dots)) old_args <- setdiff(arg_names, new_args) dots[old_args] <- object[old_args] @@ -2640,10 +2655,10 @@ log_lik.brmsfit <- function(object, newdata = NULL, re_formula = NULL, if (is.list(draws$mu[["mv"]])) { draws$mu <- get_eta(draws$mu) } - auxpars <- intersect(valid_auxpars(family(object)), names(draws)) - for (ap in auxpars) { - if (is.list(draws[[ap]])) { - draws[[ap]] <- get_auxpar(draws[[ap]]) + dpars <- intersect(valid_dpars(family(object)), names(draws)) + for (dp in dpars) { + if (is.list(draws[[dp]])) { + draws[[dp]] <- get_dpar(draws[[dp]]) } } loglik <- do.call(cbind, lapply(seq_len(N), loglik_fun, draws = draws)) @@ -2691,10 +2706,10 @@ pp_mixture.brmsfit <- function(x, newdata = NULL, re_formula = NULL, draws <- do.call(extract_draws, draws_args) draws$pp_mixture <- TRUE - auxpars <- intersect(valid_auxpars(family(x)), names(draws)) - for (ap in auxpars) { - if (is.list(draws[[ap]])) { - draws[[ap]] <- get_auxpar(draws[[ap]]) + dpars <- intersect(valid_dpars(family(x)), names(draws)) + for (dp in dpars) { + if (is.list(draws[[dp]])) { + draws[[dp]] <- get_dpar(draws[[dp]]) } } N <- choose_N(draws) diff --git a/R/brmsformula.R b/R/brmsformula.R index c4257c4ef..43826efce 100644 --- a/R/brmsformula.R +++ b/R/brmsformula.R @@ -11,10 +11,10 @@ #' a symbolic description of the model to be fitted. #' The details of model specification are given in 'Details'. #' @param ... Additional \code{formula} objects to specify -#' predictors of non-linear and auxiliary parameters. +#' predictors of non-linear and distributional parameters. #' Formulas can either be named directly or contain #' names on their left-hand side. -#' The following are auxiliary parameters of specific families +#' The following are distributional parameters of specific families #' (all other parameters are treated as non-linear parameters): #' \code{sigma} (residual standard deviation or scale of #' the \code{gaussian}, \code{student}, \code{lognormal} @@ -37,7 +37,7 @@ #' \code{bs}, \code{ndt}, and \code{bias} (boundary separation, #' non-decision time, and initial bias of the \code{wiener} #' diffusion model). -#' All auxiliary parameters are modeled +#' All distributional parameters are modeled #' on the log or logit scale to ensure correct definition #' intervals after transformation. #' See 'Details' for more explanation. @@ -55,6 +55,8 @@ #' is essentially a \code{list} containing all model #' formulas as well as some additional information. #' +#' @seealso \code{\link[brms:brmsformula-helpers]{brmsformula-helpers}} +#' #' @details #' #' \bold{General formula structure} @@ -299,7 +301,7 @@ #' #' As of \pkg{brms} 1.0.0, zero-inflated and hurdle models are specfied #' in the same way as as their non-inflated counterparts. -#' However, they have additional auxiliary parameters +#' However, they have additional distributional parameters #' (named \code{zi} and \code{hu} respectively) #' modeling the zero-inflation / hurdle probability depending on which #' model you choose. These parameters can also be affected by predictors @@ -365,15 +367,15 @@ #' For this reason it is mandatory to specify priors on the non-linear parameters. #' For instructions on how to do that, see \code{\link[brms:set_prior]{set_prior}}. #' -#' \bold{Formula syntax for predicting auxiliary parameters} +#' \bold{Formula syntax for predicting distributional parameters} #' -#' It is also possible to predict auxiliary parameters of the response +#' It is also possible to predict parameters of the response #' distribution such as the residual standard deviation \code{sigma} #' in gaussian models or the hurdle probability \code{hu} in hurdle models. #' The syntax closely resembles that of a non-linear #' parameter, for instance \code{sigma ~ x + s(z) + (1+x|g)}. #' -#' Alternatively, one may fix auxiliary parameters to certain values. +#' Alternatively, one may fix distributional parameters to certain values. #' However, this is mainly useful when models become too #' complicated and otherwise have convergence issues. #' We thus suggest to be generally careful when making use of this option. @@ -390,19 +392,19 @@ #' the negative binomial distribution with \code{shape = 1}. #' Furthermore, the parameter \code{disc} ('discrimination') in ordinal #' models is fixed to \code{1} by default and not estimated, -#' but may be modeled as any other auxiliary parameter if desired +#' but may be modeled as any other distributional parameter if desired #' (see examples). For reasons of identification, \code{'disc'} #' can only be positive, which is achieved by applying the log-link. #' -#' All auxiliary parameters currently supported by \code{brmsformula} +#' All distributional parameters currently supported by \code{brmsformula} #' have to positive (a negative standard deviation or precision parameter #' doesn't make any sense) or are bounded between 0 and 1 (for zero-inflated / #' hurdle proabilities, quantiles, or the intial bias parameter of #' drift-diffusion models). #' However, linear predictors can be positive or negative, and thus the log link #' (for positive parameters) or logit link (for probability parameters) are used -#' by default to ensure that auxiliary parameters are within their valid intervals. -#' This implies that, by default, effects for auxiliary parameters are estimated +#' by default to ensure that distributional parameters are within their valid intervals. +#' This implies that, by default, effects for distributional parameters are estimated #' on the log / logit scale and one has to apply the inverse link function to get #' to the effects on the original scale. #' Alternatively, it is possible to use the identity link to predict parameters @@ -421,17 +423,17 @@ #' terms allowed in non-mixture models are allowed in mixture models #' as well. #' -#' Auxiliary parameters of mixture distributions have the same +#' distributional parameters of mixture distributions have the same #' name as those of the corresponding ordinary distributions, but with #' a number at the end to indicate the mixture component. For instance, if -#' you use family \code{mixture(gaussian, gaussian)}, the auxiliary +#' you use family \code{mixture(gaussian, gaussian)}, the distributional #' parameters are \code{sigma1} and \code{sigma2}. -#' Auxiliary parameters of the same class can be fixed to the same value. +#' distributional parameters of the same class can be fixed to the same value. #' For the above example, we could write \code{sigma2 = "sigma1"} to make #' sure that both components have the same residual standard deviation, #' which is in turn estimated from the data. #' -#' In addition, there are two types of special auxiliary parameters. +#' In addition, there are two types of special distributional parameters. #' The first are named \code{mu}, that allow for modeling different #' predictors for the mean parameters of different mixture components. #' For instance, if you want to predict the mean of the first component @@ -510,6 +512,12 @@ #' # fix both residual standard deviations to the same value #' bf(y ~ x, sigma2 = "sigma1", family = mix) #' +#' # use the '+' operator to specify models +#' bf(y ~ 1) + +#' nlf(sigma ~ a * exp(b * x), a ~ x) + +#' lf(b ~ z + (1|g), dpar = "sigma") + +#' gaussian() +#' #' @export brmsformula <- function(formula, ..., flist = NULL, family = NULL, nl = NULL, nonlinear = NULL) { @@ -519,9 +527,11 @@ brmsformula <- function(formula, ..., flist = NULL, family = NULL, class(formula) <- "formula" } if (!is.null(nonlinear)) { - warning2("Argument 'nonlinear' is deprecated. ", - "See help(brmsformula) for the new way ", - "of specifying non-linear models.") + warning2( + "Argument 'nonlinear' is deprecated. ", + "See help(brmsformula) for the new way ", + "of specifying non-linear models." + ) } old_nonlinear <- attr(formula, "nonlinear") if (is.list(old_nonlinear)) { @@ -530,8 +540,8 @@ brmsformula <- function(formula, ..., flist = NULL, family = NULL, if (length(nonlinear)) { nl <- TRUE } - old_forms <- rmNULL(attributes(formula)[auxpars()]) - attributes(formula)[c(auxpars(), "nonlinear")] <- NULL + old_forms <- rmNULL(attributes(formula)[dpars()]) + attributes(formula)[c(dpars(), "nonlinear")] <- NULL if (is.brmsformula(formula)) { out <- formula @@ -541,73 +551,93 @@ brmsformula <- function(formula, ..., flist = NULL, family = NULL, out$pforms[names(old_forms)] <- old_forms # parse and validate dots arguments - dots <- c(list(...), flist, nonlinear) + dots <- c(out$pforms, out$pfix, list(...), flist, nonlinear) + dots <- lapply(dots, function(x) if (is.list(x)) x else list(x)) + dots <- unlist(dots, recursive = FALSE) forms <- list() for (i in seq_along(dots)) { forms <- c(forms, prepare_auxformula(dots[[i]], par = names(dots)[i])) } - dupl_pars <- names(forms)[duplicated(names(forms))] - if (length(dupl_pars)) { - dupl_pars <- collapse_comma(dupl_pars) - stop2("Duplicated specification of parameters ", dupl_pars) + is_dupl_pars <- duplicated(names(forms), fromLast = TRUE) + if (any(is_dupl_pars)) { + dupl_pars <- collapse_comma(names(forms)[is_dupl_pars]) + message("Replacing initial definitions of parameters ", dupl_pars) + forms[is_dupl_pars] <- NULL } not_form <- ulapply(forms, function(x) !is.formula(x)) fix <- forms[not_form] forms[names(fix)] <- NULL - fix_theta <- fix[auxpar_class(names(fix)) %in% "theta"] + out$pforms <- forms + # validate fixed distributional parameters + fix_theta <- fix[dpar_class(names(fix)) %in% "theta"] if (length(fix_theta)) { # normalize mixing proportions sum_theta <- sum(unlist(fix_theta)) fix_theta <- lapply(fix_theta, "/", sum_theta) fix[names(fix_theta)] <- fix_theta } - out$pforms[names(forms)] <- forms - out$pfix[names(fix)] <- fix - for (ap in names(out$pfix)) { - if (is.character(out$pfix[[ap]])) { - if (identical(ap, out$pfix[[ap]])) { - stop2("Equating '", ap, "' with itself is not meaningful.") + out$pfix <- fix + for (dp in names(out$pfix)) { + if (is.character(out$pfix[[dp]])) { + if (identical(dp, out$pfix[[dp]])) { + stop2("Equating '", dp, "' with itself is not meaningful.") } - ap_class <- auxpar_class(ap) + ap_class <- dpar_class(dp) if (ap_class == "mu") { stop2("Equating parameters of class 'mu' is not allowed.") } - if (!identical(ap_class, auxpar_class(out$pfix[[ap]]))) { + if (!identical(ap_class, dpar_class(out$pfix[[dp]]))) { stop2("Can only equate parameters of the same class.") } - if (out$pfix[[ap]] %in% names(out$pfix)) { + if (out$pfix[[dp]] %in% names(out$pfix)) { stop2("Cannot use fixed parameters on ", "the right-hand side of an equation.") } - if (out$pfix[[ap]] %in% names(out$pforms)) { + if (out$pfix[[dp]] %in% names(out$pforms)) { stop2("Cannot use predicted parameters on ", "the right-hand side of an equation.") } } } if (!is.null(nl)) { - nl <- as.logical(nl)[1] - if (!nl %in% c(TRUE, FALSE)) { - stop2("Argument 'nl' could not be coerced to logical.") - } - out[["nl"]] <- nl + attr(out$formula, "nl") <- as_one_logical(nl) + } else if (!is.null(out[["nl"]])) { + # for backwards compatibility with brms <= 1.8.0 + attr(out$formula, "nl") <- out[["nl"]] + out[["nl"]] <- NULL + } + if (is.null(attr(out$formula, "nl"))) { + attr(out$formula, "nl") <- FALSE } if (!is.null(family)) { out[["family"]] <- check_family(family) } if (!is.null(out[["family"]])) { # check for the presence of non-linear parameters - auxpars <- is_auxpar_name(names(out$pforms), out$family) - auxpars <- names(out$pforms)[auxpars] - for (ap in names(out$pforms)) { - if (!ap %in% auxpars) { - if (!isTRUE(out$nl)) { - stop2("The parameter '", ap, "' is not an auxiliary parameter.\n", - "Set 'nl = TRUE' if you want to specify non-linear models.") + dpars <- is_dpar_name(names(out$pforms), out$family) + dpars <- names(out$pforms)[dpars] + for (dp in names(out$pforms)) { + if (!dp %in% dpars) { + # indicate the correspondence to distributional parameter + if (is.null(attr(out$pforms[[dp]], "dpar"))) { + attr(out$pforms[[dp]], "dpar") <- "mu" + } + dpar <- attr(out$pforms[[dp]], "dpar") + if (!is.null(out$pforms[[dpar]])) { + nl_allowed <- isTRUE(attr(out$pforms[[dpar]], "nl")) + } else { + if (dpar_class(dpar) == "mu") { + nl_allowed <- isTRUE(attr(out$formula, "nl")) + } else { + nl_allowed <- FALSE + } } - # indicate the correspondence to an auxiliary parameter - if (is.null(attr(out$pforms[[ap]], "auxpar"))) { - attr(out$pforms[[ap]], "auxpar") <- "mu" + if (!nl_allowed) { + stop2( + "The parameter '", dp, "' is not a valid ", + "distributional or non-linear parameter. ", + "Did you forget to set 'nl = TRUE'?" + ) } } } @@ -615,7 +645,7 @@ brmsformula <- function(formula, ..., flist = NULL, family = NULL, # add default values for unspecified elements defs <- list( pforms = list(), pfix = list(), family = NULL, - nl = FALSE, response = NULL, old_mv = FALSE + response = NULL, old_mv = FALSE ) defs <- defs[setdiff(names(defs), names(rmNULL(out, FALSE)))] out[names(defs)] <- defs @@ -631,8 +661,117 @@ bf <- function(formula, ..., flist = NULL, family = NULL, nl = nl, nonlinear = nonlinear) } +#' Linear and Non-linear formulas in \pkg{brms} +#' +#' Helper functions to specify linear and non-linear +#' formulas for use with \code{\link[brms:brmsformula]{brmsformula}}. +#' +#' @name brmsformula-helpers +#' @aliases bf-helpers nlf lf set_nl +#' +#' @param formula Non-linear formula for a distributional parameter. +#' The name of the distributional parameter can either be specified +#' on the left-hand side of \code{formula} or via argument \code{dpar}. +#' @param dpar Optional character string specifying the distributional +#' parameter to which the formulas passed via \code{...} and +#' \code{flist} belong. +#' @inheritParams brmsformula +#' +#' @return For \code{lf} and \code{nlf} a \code{list} that can be +#' passed to \code{\link[brms:brmsformula]{brmsformula}} or added +#' to an existing \code{brmsformula} object. For \code{set_nl} +#' a \code{list} that can be added to an existing +#' \code{brmsformula} object. +#' +#' @seealso \code{\link[brms:brmsformula]{brmsformula}} +#' +#' @examples +#' # add more formulas to the model +#' bf(y ~ 1) + +#' nlf(sigma ~ a * exp(b * x), a ~ x) + +#' lf(b ~ z + (1|g), dpar = "sigma") + +#' gaussian() +#' +#' # specify 'nl' later on +#' bf(y ~ a * inv_logit(x * b)) + +#' lf(a + b ~ z) + +#' set_nl(TRUE) +#' +#' @export +nlf <- function(formula, ..., flist = NULL, dpar = NULL) { + formula <- as.formula(formula) + resp_pars <- all.vars(formula[[2]]) + if (length(resp_pars) == 0L) { + if (is.null(dpar)) { + stop2("No parameter name passed via the LHS of ", + "'formula' or argument 'dpar'.") + } + } else if (length(resp_pars) == 1L) { + dpar <- resp_pars + } else { + stop2("LHS of non-linear formula should contain only one variable.") + } + attr(formula, "nl") <- TRUE + c(setNames(list(formula), dpar), lf(..., flist = flist, dpar = dpar)) +} + +#' @rdname brmsformula-helpers +#' @export +lf <- function(..., flist = NULL, dpar = NULL) { + out <- c(list(...), flist) + if (!is.null(dpar)) { + dpar <- as.character(dpar) + if (length(dpar) != 1L) { + stop2("Argument 'dpar' should be of length 1.") + } + for (i in seq_along(out)) { + attr(out[[i]], "dpar") <- dpar + } + } + out +} + +#' @rdname brmsformula-helpers +#' @export +set_nl <- function(nl = TRUE, dpar = NULL) { + nl <- as_one_logical(nl) + if (!is.null(dpar)) { + dpar <- as.character(dpar) + if (length(dpar) != 1L) { + stop2("Argument 'dpar' should be of length 1.") + } + } + structure(nlist(nl, dpar), class = "setnl") +} + +#' @export +"+.brmsformula" <- function(e1, e2) { + if (is.function(e2)) { + e2 <- try(e2(), silent = TRUE) + if (!is.family(e2)) { + stop2("Don't know how to handle non-family functions.") + } + } + if (is.family(e2)) { + e1 <- bf(e1, family = e2) + } else if (inherits(e2, "setnl")) { + if (is.null(e2$dpar)) { + e1 <- bf(e1, nl = e2$nl) + } else { + if (is.null(e1$pforms[[e2$dpar]])) { + stop2("Parameter '", e2$dpar, "' has no formula.") + } + attr(e1$pforms[[e2$dpar]], "nl") <- e2$nl + e1 <- bf(e1) + } + } else { + e1 <- bf(e1, e2) + } + e1 +} + prepare_auxformula <- function(formula, par = NULL, rsv_pars = NULL) { - # validate and prepare a formula of an auxiliary parameter + # validate and prepare a formula of an distributional parameter # Args: # formula: an object of class formula # par: optional name of the parameter; if not specified @@ -642,7 +781,7 @@ prepare_auxformula <- function(formula, par = NULL, rsv_pars = NULL) { try_formula <- try(as.formula(formula), silent = TRUE) if (is(try_formula, "try-error")) { if (length(formula) != 1L) { - stop2("Expecting a single value when fixing auxiliary parameters.") + stop2("Expecting a single value when fixing parameters.") } scalar <- SW(as.numeric(formula)) if (!is.na(scalar)) { @@ -679,16 +818,16 @@ prepare_auxformula <- function(formula, par = NULL, rsv_pars = NULL) { out } -auxpars <- function() { - # names of auxiliary parameters +dpars <- function() { + # names of distributional parameters c("mu", "sigma", "shape", "nu", "phi", "kappa", "beta", "xi", "zi", "hu", "zoi", "coi", "disc", "bs", "ndt", "bias", "quantile", "alpha", "theta") } -links_auxpars <- function(ap) { - # link functions for auxiliary parameters - switch(ap, +links_dpars <- function(dp) { + # link functions for distributional parameters + switch(dp, mu = "identity", sigma = c("log", "identity"), shape = c("log", "identity"), @@ -708,13 +847,13 @@ links_auxpars <- function(ap) { xi = c("log1p", "identity"), alpha = c("identity", "log"), theta = c("identity"), - stop2("Parameter '", ap, "' is not supported.") + stop2("Parameter '", dp, "' is not supported.") ) } #' @export -valid_auxpars.default <- function(family, bterms = NULL, ...) { - # convenience function to find relevant auxiliary parameters +valid_dpars.default <- function(family, bterms = NULL, ...) { + # convenience function to find relevant distributional parameters x <- c( mu = TRUE, sigma = has_sigma(family, bterms = bterms), @@ -739,45 +878,45 @@ valid_auxpars.default <- function(family, bterms = NULL, ...) { } #' @export -valid_auxpars.mixfamily <- function(family, ...) { - out <- lapply(family$mix, valid_auxpars, ...) +valid_dpars.mixfamily <- function(family, ...) { + out <- lapply(family$mix, valid_dpars, ...) for (i in seq_along(out)) { out[[i]] <- paste0(out[[i]], i) } c(unlist(out), paste0("theta", seq_along(out))) } -is_auxpar_name <- function(auxpars, family = NULL, ...) { +is_dpar_name <- function(dpars, family = NULL, ...) { # check if provided names of auxiliar parameters are valid # Args: - # auxpars: character vector to be checked + # dpars: character vector to be checked # family: the model family - # ...: further arguments passed to valid_auxpars - auxpars <- as.character(auxpars) - if (!length(auxpars)) { + # ...: further arguments passed to valid_dpars + dpars <- as.character(dpars) + if (!length(dpars)) { return(logical(0)) } if (is.null(family)) { - patterns <- paste0("^", auxpars(), "[[:digit:]]*$") - .is_auxpar_name <- function(auxpar, ...) { - any(ulapply(patterns, grepl, x = auxpar)) + patterns <- paste0("^", dpars(), "[[:digit:]]*$") + .is_dpar_name <- function(dpar, ...) { + any(ulapply(patterns, grepl, x = dpar)) } - out <- ulapply(auxpars, .is_auxpar_name) + out <- ulapply(dpars, .is_dpar_name) } else { - out <- auxpars %in% valid_auxpars(family, ...) + out <- dpars %in% valid_dpars(family, ...) } as.logical(out) } -auxpar_class <- function(auxpar) { - # class of an auxiliary parameter - out <- get_matches("^[^[:digit:]]+", auxpar, simplify = FALSE) +dpar_class <- function(dpar) { + # class of an distributional parameter + out <- get_matches("^[^[:digit:]]+", dpar, simplify = FALSE) ulapply(out, function(x) ifelse(length(x), x, "")) } -auxpar_id <- function(auxpar) { - # id of an auxiliary parameter - out <- get_matches("[[:digit:]]+$", auxpar, simplify = FALSE) +dpar_id <- function(dpar) { + # id of an distributional parameter + out <- get_matches("[[:digit:]]+$", dpar, simplify = FALSE) ulapply(out, function(x) ifelse(length(x), x, "")) } @@ -805,18 +944,19 @@ amend_formula <- function(formula, data = NULL, family = gaussian(), if (is.null(out$family)) { out <- bf(out, family = family) } - # to allow the '.' symbol in formula - try_terms <- try(terms(out$formula, data = data), silent = TRUE) - if (!is(try_terms, "try-error")) { - out$formula <- formula(try_terms) + # allow the '.' symbol in the formulas + out$formula <- expand_dot_formula(out$formula) + for (i in seq_along(out$pforms)) { + out$pforms[[i]] <- expand_dot_formula(out$pforms[[i]]) } if (is_ordinal(out$family)) { # fix discrimination to 1 by default if (!"disc" %in% c(names(pforms(out)), names(pfix(out)))) { out <- bf(out, disc = 1) } - no_intercept <- isTRUE(attr(try_terms, "intercept", TRUE) == 0) - if (!is(try_terms, "try-error") && no_intercept) { + try_terms <- try(stats::terms(out$formula), silent = TRUE) + intercept <- attr(try_terms, "intercept", TRUE) + if (!is(try_terms, "try-error") && isTRUE(intercept == 0)) { stop2("Cannot remove the intercept in an ordinal model.") } } @@ -916,5 +1056,5 @@ is.brmsformula <- function(x) { is_nonlinear <- function(x) { stopifnot(is.brmsfit(x)) - bf(x$formula)[["nl"]] + isTRUE(attr(bf(x$formula)$formula, "nl", TRUE)) } diff --git a/R/data-helpers.R b/R/data-helpers.R index c2a6b2f66..c05db8bd4 100644 --- a/R/data-helpers.R +++ b/R/data-helpers.R @@ -7,7 +7,7 @@ update_data <- function(data, bterms, na.action = na.omit, # bterms: object of class brmsterms # na.action: function defining how to treat NAs # drop.unused.levels: indicates if unused factor levels - # should be removed + # should be removed # terms_attr: a list of attributes of the terms object of # the original model.frame; only used with newdata; # this ensures that (1) calls to 'poly' work correctly @@ -15,7 +15,7 @@ update_data <- function(data, bterms, na.action = na.omit, # of variable names; fixes issue #73 # knots: a list of knot values for GAMMS # Returns: - # model.frame in long format with combined grouping variables if present + # model.frame for use in brms functions if (missing(data)) { stop2("Argument 'data' is missing.") } @@ -45,8 +45,10 @@ update_data <- function(data, bterms, na.action = na.omit, stop2("The following variables are missing in 'data':\n", collapse_comma(missing_vars)) } - data <- model.frame(bterms$allvars, data, na.action = na.pass, - drop.unused.levels = drop.unused.levels) + data <- model.frame( + bterms$allvars, data, na.action = na.pass, + drop.unused.levels = drop.unused.levels + ) nrow_with_NA <- nrow(data) data <- na.action(data) if (nrow(data) != nrow_with_NA) { @@ -146,11 +148,13 @@ check_data_old_mv <- function(data, bterms) { } data_rsv_intercept <- function(data, bterms) { - # introduces the resevered variable "intercept" to the data + # add the resevered variable 'intercept' to the data # Args: # data: data.frame or list # bterms: object of class brmsterms - rsv_int <- ulapply(bterms$auxpars, function(x) attr(x$fe, "rsv_intercept")) + rsv_int <- ulapply(bterms$dpars, + function(x) attr(x$fe, "rsv_intercept") + ) if (any(rsv_int)) { if ("intercept" %in% names(data)) { stop2("Variable name 'intercept' is resevered in models ", @@ -172,7 +176,7 @@ combine_groups <- function(data, ...) { group <- c(...) for (i in seq_along(group)) { sgroup <- unlist(strsplit(group[[i]], ":")) - if (length(sgroup) > 1) { + if (length(sgroup) > 1L) { new.var <- get(sgroup[1], data) for (j in 2:length(sgroup)) { new.var <- paste0(new.var, "_", get(sgroup[j], data)) @@ -189,7 +193,7 @@ fix_factor_contrasts <- function(data, optdata = NULL, ignore = NULL) { # Args: # data: a data.frame # optdata: optional data.frame from which contrasts - # are taken if present + # are taken if present # ignore: names of variables for which not to fix contrasts # Returns: # a data.frame with amended contrasts attributes @@ -254,10 +258,13 @@ amend_newdata <- function(newdata, fit, re_formula = NULL, } # standata will be based on an updated formula if re_formula is specified new_formula <- update_re_terms(formula(fit), re_formula = re_formula) - bterms <- parse_bf(new_formula, family = family(fit), - autocor = fit$autocor, resp_rhs_all = FALSE) - resp_only_vars <- setdiff(all.vars(bterms$respform), - all.vars(rhs(bterms$allvars))) + bterms <- parse_bf( + new_formula, family = family(fit), + autocor = fit$autocor, resp_rhs_all = FALSE + ) + resp_only_vars <- setdiff( + all.vars(bterms$respform), all.vars(rhs(bterms$allvars)) + ) # fixes #162 resp_only_vars <- c(resp_only_vars, all.vars(bterms$adforms$dec)) missing_resp <- setdiff(resp_only_vars, names(newdata)) @@ -322,9 +329,11 @@ amend_newdata <- function(newdata, fit, re_formula = NULL, old_contrasts <- rbind(old_contrasts, zero__ = 0) } if (any(!new_levels %in% old_levels)) { - stop2("New factor levels are not allowed.", - "\nLevels found: ", paste(new_levels, collapse = ", ") , - "\nLevels allowed: ", paste(old_levels, collapse = ", ")) + stop2( + "New factor levels are not allowed.", + "\nLevels allowed: ", paste(old_levels, collapse = ", "), + "\nLevels found: ", paste(new_levels, collapse = ", ") + ) } newdata[[factor_names[i]]] <- factor(new_factor, old_levels) # don't use contrasts(.) here to avoid dimension checks @@ -342,11 +351,13 @@ amend_newdata <- function(newdata, fit, re_formula = NULL, new_values <- get(v, newdata) min_value <- min(list_data[[v]]) invalid <- new_values < min_value | - new_values > max(list_data[[v]]) | - !is_wholenumber(new_values) + new_values > max(list_data[[v]]) | + !is_wholenumber(new_values) if (sum(invalid)) { - stop2("Invalid values in variable '", v, "': ", - paste(new_values[invalid], collapse = ",")) + stop2( + "Invalid values in variable '", v, "': ", + paste0(new_values[invalid], collapse = ", ") + ) } attr(newdata[[v]], "min") <- min_value } @@ -368,9 +379,11 @@ amend_newdata <- function(newdata, fit, re_formula = NULL, unknown_levels <- setdiff(new_levels[[g]], old_levels[[g]]) if (!allow_new_levels && length(unknown_levels)) { unknown_levels <- paste0("'", unknown_levels, "'", collapse = ", ") - stop2("Levels ", unknown_levels, " of grouping factor '", g, "' ", - "cannot be found in the fitted model. ", - "Consider setting argument 'allow_new_levels' to TRUE.") + stop2( + "Levels ", unknown_levels, " of grouping factor '", g, "' ", + "cannot be found in the fitted model. ", + "Consider setting argument 'allow_new_levels' to TRUE." + ) } } if (return_standata) { @@ -386,16 +399,10 @@ amend_newdata <- function(newdata, fit, re_formula = NULL, old_terms <- attr(model.frame(fit), "terms") terms_attr <- c("variables", "predvars") control$terms_attr <- attributes(old_terms)[terms_attr] - # some components should not be computed based on newdata - has_mo <- length(get_effect(bterms, "mo")) > 0L - if (has_trials(fit$family) || has_cat(fit$family) || has_mo) { - pars <- c(names(bterms$nlpars), intersect(auxpars(), names(bterms))) - comp <- c("trials", "ncat", paste0("Jmo", c("", paste0("_", pars)))) - old_standata <- rmNULL(standata(fit)[comp]) + if (has_trials(fit$family) || has_cat(fit$family)) { + # trials and cat should not be computed based on newdata + old_standata <- standata(fit, control = list(only_response = TRUE)) control[c("trials", "ncat")] <- old_standata[c("trials", "ncat")] - Jmo <- old_standata[grepl("^Jmo", names(old_standata))] - names(Jmo) <- sub("^Jmo$", "mu", sub("^Jmo_", "", names(Jmo))) - control[["Jmo"]] <- Jmo } if (is.cor_car(fit$autocor)) { if (isTRUE(nzchar(bterms$time$group))) { @@ -405,6 +412,7 @@ amend_newdata <- function(newdata, fit, re_formula = NULL, } control$smooths <- make_smooth_list(bterms, model.frame(fit)) control$gps <- make_gp_list(bterms, model.frame(fit)) + control$Jmo <- make_Jmo_list(bterms, model.frame(fit)) knots <- attr(model.frame(fit), "knots") newdata <- make_standata( new_formula, data = newdata, family = fit$family, @@ -479,14 +487,12 @@ get_model_matrix <- function(formula, data = environment(formula), X } -prepare_mo_vars <- function(formula, data, check = TRUE) { - # prepare monotonic variables for use in Stan +mo_design_matrix <- function(formula, data, check = TRUE) { + # compute the design matrix of monotonic effects # Args: # formula: formula containing mononotic terms # data: a data.frame or named list # check: check the number of levels? - # Returns: - # 'data' with amended monotonic variables data <- model.frame(formula, data) vars <- names(data) for (i in seq_along(vars)) { @@ -501,13 +507,17 @@ prepare_mo_vars <- function(formula, data, check = TRUE) { } data[[vars[i]]] <- data[[vars[i]]] - min_value } else { - stop2("Monotonic predictors must be either integers or ", - "ordered factors. Error occured for variable '", - vars[i], "'.") + stop2( + "Monotonic predictors must be either integers or ", + "ordered factors. Error occured for variable '", + vars[i], "'." + ) } if (check && max(data[[vars[i]]]) < 2L) { - stop2("Monotonic predictors must have at least 3 different ", - "values. Error occured for variable '", vars[i], "'.") + stop2( + "Monotonic predictors must have at least 3 different ", + "values. Error occured for variable '", vars[i], "'." + ) } } out <- get_model_matrix(formula, data, cols2remove = "(Intercept)") @@ -517,13 +527,75 @@ prepare_mo_vars <- function(formula, data, check = TRUE) { out } +arr_design_matrix <- function(Y, r, group) { + # compute the design matrix for ARR effects + # Args: + # Y: a vector containing the response variable + # r: ARR order + # group: vector containing the grouping variable + # Notes: + # expects Y to be sorted after group already + # Returns: + # the design matrix for ARR effects + stopifnot(length(Y) == length(group)) + if (r > 0) { + U_group <- unique(group) + N_group <- length(U_group) + out <- matrix(0, nrow = length(Y), ncol = r) + ptsum <- rep(0, N_group + 1) + for (j in seq_len(N_group)) { + ptsum[j + 1] <- ptsum[j] + sum(group == U_group[j]) + for (i in seq_len(r)) { + if (ptsum[j] + i + 1 <= ptsum[j + 1]) { + out[(ptsum[j] + i + 1):ptsum[j + 1], i] <- + Y[(ptsum[j] + 1):(ptsum[j + 1] - i)] + } + } + } + } else { + out <- NULL + } + out +} + +#' @export +make_Jmo_list.btl <- function(x, data, ...) { + if (!is.null(x$mo)) { + Xmo <- mo_design_matrix(x$mo, data, ...) + out <- as.array(apply(Xmo, 2, max)) + } else { + out <- NULL + } + out +} + +#' @export +make_Jmo_list.btnl <- function(x, data, ...) { + out <- named_list(names(x$nlpars)) + for (i in seq_along(out)) { + out[[i]] <- make_Jmo_list(x$nlpars[[i]], data, ...) + } + out +} + +#' @export +make_Jmo_list.brmsterms <- function(x, data, ...) { + out <- named_list(names(x$dpars)) + for (i in seq_along(out)) { + out[[i]] <- make_Jmo_list(x$dpars[[i]], data, ...) + } + out +} + #' @export make_smooth_list.btl <- function(x, data, ...) { if (has_smooths(x)) { knots <- attr(data, "knots") data <- rm_attr(data, "terms") - gam_args <- list(data = data, knots = knots, - absorb.cons = TRUE, modCon = 3) + gam_args <- list( + data = data, knots = knots, + absorb.cons = TRUE, modCon = 3 + ) sm_labels <- get_sm_labels(x) out <- named_list(sm_labels) for (i in seq_along(sm_labels)) { @@ -547,9 +619,9 @@ make_smooth_list.btnl <- function(x, data, ...) { #' @export make_smooth_list.brmsterms <- function(x, data, ...) { - out <- named_list(names(x$auxpars)) + out <- named_list(names(x$dpars)) for (i in seq_along(out)) { - out[[i]] <- make_smooth_list(x$auxpars[[i]], data, ...) + out[[i]] <- make_smooth_list(x$dpars[[i]], data, ...) } out } @@ -577,40 +649,9 @@ make_gp_list.btnl <- function(x, data, ...) { #' @export make_gp_list.brmsterms <- function(x, data, ...) { - out <- named_list(names(x$auxpars)) + out <- named_list(names(x$dpars)) for (i in seq_along(out)) { - out[[i]] <- make_gp_list(x$auxpars[[i]], data, ...) + out[[i]] <- make_gp_list(x$dpars[[i]], data, ...) } out } - -arr_design_matrix <- function(Y, r, group) { - # compute the design matrix for ARR effects - # Args: - # Y: a vector containing the response variable - # r: ARR order - # group: vector containing the grouping variable - # Notes: - # expects Y to be sorted after group already - # Returns: - # the design matrix for ARR effects - stopifnot(length(Y) == length(group)) - if (r > 0) { - U_group <- unique(group) - N_group <- length(U_group) - out <- matrix(0, nrow = length(Y), ncol = r) - ptsum <- rep(0, N_group + 1) - for (j in seq_len(N_group)) { - ptsum[j + 1] <- ptsum[j] + sum(group == U_group[j]) - for (i in seq_len(r)) { - if (ptsum[j] + i + 1 <= ptsum[j + 1]) { - out[(ptsum[j] + i + 1):ptsum[j + 1], i] <- - Y[(ptsum[j] + 1):(ptsum[j + 1] - i)] - } - } - } - } else { - out <- NULL - } - out -} diff --git a/R/data-predictor.R b/R/data-predictor.R index e6b0bbb1f..bb9aa9bef 100644 --- a/R/data-predictor.R +++ b/R/data-predictor.R @@ -1,8 +1,8 @@ #' @export data_effects.btl <- function(x, data, ranef = empty_ranef(), prior = brmsprior(), knots = NULL, - nlpar = "", not4stan = FALSE, - smooths = NULL, gps = NULL, Jmo = NULL) { + not4stan = FALSE, smooths = NULL, + gps = NULL, Jmo = NULL) { # prepare data for all types of effects for use in Stan # Args: # data: the data passed by the user @@ -18,28 +18,23 @@ data_effects.btl <- function(x, data, ranef = empty_ranef(), # Jmo: optional precomputed values of Jmo for monotonic effects # Returns: # A named list of data to be passed to Stan - nlpar <- check_nlpar(nlpar) - data_fe <- data_fe(x, data = data, knots = knots, - nlpar = nlpar, not4stan = not4stan, - smooths = smooths) - data_mo <- data_mo(x, data = data, ranef = ranef, - prior = prior, Jmo = Jmo, nlpar = nlpar) - data_re <- data_re(ranef, data = data, nlpar = nlpar, - not4stan = not4stan) - data_me <- data_me(x, data = data, nlpar = nlpar) - data_cs <- data_cs(x, data = data, nlpar = nlpar) - data_gp <- data_gp(x, data = data, nlpar = nlpar, gps = gps) - data_offset <- data_offset(x, data = data, nlpar = nlpar) - data_prior <- data_prior(prior, data = data, nlpar = nlpar) - c(data_fe, data_mo, data_re, data_me, data_cs, - data_gp, data_offset, data_prior) + c( + data_fe(x, data, knots = knots, not4stan = not4stan, smooths = smooths), + data_mo(x, data, ranef = ranef, prior = prior, Jmo = Jmo), + data_re(x, data, ranef = ranef, not4stan = not4stan), + data_me(x, data), + data_cs(x, data), + data_gp(x, data, gps = gps), + data_offset(x, data), + data_prior(x, data, prior = prior) + ) } #' @export data_effects.btnl <- function(x, data, ranef = empty_ranef(), prior = brmsprior(), knots = NULL, - nlpar = "", not4stan = FALSE, - smooths = NULL, gps = NULL, Jmo = NULL) { + not4stan = FALSE, smooths = NULL, + gps = NULL, Jmo = NULL) { # prepare data for non-linear parameters for use in Stan # matrix of covariates appearing in the non-linear formula out <- list() @@ -49,14 +44,15 @@ data_effects.btnl <- function(x, data, ranef = empty_ranef(), } # fixes issue #127 occuring for factorial covariates colnames(C) <- all.vars(x$covars) + p <- usc(combine_prefix(x)) if (not4stan) { - out <- c(out, nlist(C)) + out[[paste0("C", p)]] <- C } else { # use vectors as indexing matrices in Stan is slow if (ncol(C)) { out <- c(out, setNames( as.list(as.data.frame(C)), - paste0("C_", seq_len(ncol(C))) + paste0("C", p, "_", seq_len(ncol(C))) )) } } @@ -64,22 +60,21 @@ data_effects.btnl <- function(x, data, ranef = empty_ranef(), out <- c(out, data_effects( x$nlpars[[nlp]], data, ranef = ranef, - prior = prior, knots = knots, nlpar = nlp, - not4stan = not4stan, smooths = smooths[[nlp]], - gps = gps, Jmo = Jmo + prior = prior, knots = knots, not4stan = not4stan, + smooths = smooths[[nlp]], gps = gps[[nlp]], + Jmo = Jmo[[nlp]] ) ) } out } -data_fe <- function(bterms, data, knots = NULL, nlpar = "", +data_fe <- function(bterms, data, knots = NULL, not4stan = FALSE, smooths = NULL) { # prepare data for fixed effects for use in Stan # Args: see data_effects - stopifnot(length(nlpar) == 1L) out <- list() - p <- usc(nlpar, "prefix") + p <- usc(combine_prefix(bterms)) is_ordinal <- is_ordinal(bterms$family) is_bsts <- inherits(bterms$autocor, "cor_bsts") # the intercept is removed inside the Stan code for ordinal models @@ -129,33 +124,36 @@ data_fe <- function(bterms, data, knots = NULL, nlpar = "", X <- structure(X, smooth_cols = scols, by_levels = by_levels) colnames(X) <- rename(colnames(X)) } - avoid_auxpars(colnames(X), bterms = bterms) + avoid_dpars(colnames(X), bterms = bterms) c(out, setNames(list(ncol(X), X), paste0(c("K", "X"), p))) } data_mo <- function(bterms, data, ranef = empty_ranef(), - prior = brmsprior(), nlpar = "", - Jmo = NULL) { + prior = brmsprior(), Jmo = NULL) { # prepare data for monotonic effects for use in Stan # Args: see data_effects - stopifnot(length(nlpar) == 1L) - p <- if (nchar(nlpar)) paste0("_", nlpar) else "" + px <- check_prefix(bterms) + p <- usc(combine_prefix(px)) out <- list() if (is.formula(bterms[["mo"]])) { - Xmo <- prepare_mo_vars(bterms$mo, data, check = is.null(Jmo)) - avoid_auxpars(colnames(Xmo), bterms = bterms) + Xmo <- mo_design_matrix(bterms$mo, data, check = is.null(Jmo)) + avoid_dpars(colnames(Xmo), bterms = bterms) if (is.null(Jmo)) { Jmo <- as.array(apply(Xmo, 2, max)) } out <- c(out, - setNames(list(ncol(Xmo), Xmo, Jmo), paste0(c("Kmo", "Xmo", "Jmo"), p)) + setNames( + list(ncol(Xmo), Xmo, Jmo), + paste0(c("Kmo", "Xmo", "Jmo"), p) + ) ) # validate and assign vectors for dirichlet prior monef <- colnames(Xmo) for (i in seq_along(monef)) { - take <- prior$class == "simplex" & prior$coef == monef[i] & - prior$nlpar == nlpar - simplex_prior <- prior$prior[take] + simplex_prior <- subset2(prior, + class = "simplex", coef = monef[i], ls = px + ) + simplex_prior <- simplex_prior$prior if (isTRUE(nzchar(simplex_prior))) { simplex_prior <- eval2(simplex_prior) if (length(simplex_prior) != Jmo[i]) { @@ -171,20 +169,21 @@ data_mo <- function(bterms, data, ranef = empty_ranef(), out } -data_re <- function(ranef, data, nlpar = "", not4stan = FALSE) { +data_re <- function(bterms, data, ranef, not4stan = FALSE) { # prepare data for group-level effects for use in Stan # Args: see data_effects - stopifnot(length(nlpar) == 1L) out <- list() - take <- ranef$nlpar == nlpar & !ranef$type %in% c("mo", "me") + px <- check_prefix(bterms) + take <- find_rows(ranef, ls = px) & + !find_rows(ranef, type = c("mo", "me")) ranef <- ranef[take, ] if (nrow(ranef)) { - Z <- lapply(ranef[!duplicated(ranef$gn), ]$form, - get_model_matrix, data = data) + unique_forms <- ranef[!duplicated(ranef$gn), "form"] + Z <- lapply(unique_forms, get_model_matrix, data = data) gn <- unique(ranef$gn) for (i in seq_along(gn)) { - r <- ranef[ranef$gn == gn[i], ] - idp <- paste0(r$id[1], usc(nlpar, "prefix")) + r <- subset2(ranef, gn = gn[i]) + idp <- paste0(r$id[1], usc(combine_prefix(px))) if (isTRUE(not4stan)) { # for internal use in S3 methods if (ncol(Z[[i]]) == 1L) { @@ -220,7 +219,7 @@ data_gr <- function(ranef, data, cov_ranef = NULL) { out <- list() ids <- unique(ranef$id) for (id in ids) { - id_ranef <- ranef[ranef$id == id, ] + id_ranef <- subset2(ranef, id = id) nranef <- nrow(id_ranef) group <- id_ranef$group[1] levels <- attr(ranef, "levels")[[group]] @@ -232,8 +231,10 @@ data_gr <- function(ranef, data, cov_ranef = NULL) { scale <- isTRUE(attr(weights, "scale")) weights <- as.matrix(eval_rhs(weights, data)) if (!identical(dim(weights), c(nrow(data), ngs))) { - stop2("Grouping structure 'mm' expects 'weights' to be a matrix ", - "with as many columns as grouping factors.") + stop2( + "Grouping structure 'mm' expects 'weights' to be ", + "a matrix with as many columns as grouping factors." + ) } if (scale) { if (any(weights < 0)) { @@ -291,32 +292,33 @@ data_gr <- function(ranef, data, cov_ranef = NULL) { out } -data_cs <- function(bterms, data, nlpar = "") { +data_cs <- function(bterms, data) { # prepare data for category specific effects for use in Stan # Args: # bterms: object of class brmsterms # data: the data passed by the user out <- list() + px <- check_prefix(bterms) if (length(all_terms(bterms[["cs"]]))) { - stopifnot(!nzchar(nlpar)) Xcs <- get_model_matrix(bterms$cs, data) - avoid_auxpars(colnames(Xcs), bterms = bterms) + avoid_dpars(colnames(Xcs), bterms = bterms) out <- c(out, list(Kcs = ncol(Xcs), Xcs = Xcs)) } out } -data_me <- function(bterms, data, nlpar = "") { +data_me <- function(bterms, data) { # prepare data of meausurement error variables for use in Stan # Args: # bterms: object of class brmsterms # data: the data passed by the user out <- list() + px <- check_prefix(bterms) meef <- get_me_labels(bterms, data) if (length(meef)) { - p <- usc(nlpar, "prefix") + p <- usc(combine_prefix(px)) Cme <- get_model_matrix(bterms$me, data) - avoid_auxpars(colnames(Cme), bterms = bterms) + avoid_dpars(colnames(Cme), bterms = bterms) Cme <- Cme[, attr(meef, "not_one"), drop = FALSE] Cme <- lapply(seq_len(ncol(Cme)), function(i) Cme[, i]) if (length(Cme)) { @@ -337,14 +339,15 @@ data_me <- function(bterms, data, nlpar = "") { out } -data_gp <- function(bterms, data, nlpar = "", gps = NULL) { +data_gp <- function(bterms, data, gps = NULL) { # prepare data for Gaussian process terms # Args: # see data_effects.btl out <- list() + px <- check_prefix(bterms) gpef <- get_gp_labels(bterms) if (length(gpef)) { - p <- usc(nlpar, "prefix") + p <- usc(combine_prefix(px)) for (i in seq_along(gpef)) { pi <- paste0(p, "_", i) gp <- eval2(gpef[i]) @@ -388,12 +391,14 @@ data_gp <- function(bterms, data, nlpar = "", gps = NULL) { out } -data_offset <- function(bterms, data, nlpar = "") { +data_offset <- function(bterms, data) { # prepare data of offsets for use in Stan out <- list() + px <- check_prefix(bterms) if (is.formula(bterms$offset)) { + p <- usc(combine_prefix(px)) mf <- model.frame(bterms$offset, rm_attr(data, "terms")) - out[[paste0("offset", usc(nlpar))]] <- model.offset(mf) + out[[paste0("offset", p)]] <- model.offset(mf) } out } @@ -404,10 +409,10 @@ data_mixture <- function(bterms, prior = brmsprior()) { out <- list() if (is.mixfamily(bterms$family)) { families <- family_names(bterms$family) - ap_classes <- auxpar_class( - names(c(bterms$auxpars, bterms$fauxpars)) + dp_classes <- dpar_class( + names(c(bterms$dpars, bterms$fdpars)) ) - if (!any(ap_classes %in% "theta")) { + if (!any(dp_classes %in% "theta")) { # estimate mixture probabilities directly take <- prior$class == "theta" theta_prior <- prior$prior[take] @@ -426,12 +431,13 @@ data_mixture <- function(bterms, prior = brmsprior()) { out } -data_prior <- function(prior, data, nlpar = "") { +data_prior <- function(bterms, data, prior) { # data for special priors such as horseshoe and lasso out <- list() - orig_nlpar <- ifelse(nzchar(nlpar), nlpar, "mu") - nlpar <- usc(nlpar) - special <- attr(prior, "special")[[orig_nlpar]] + px <- check_prefix(bterms) + p <- usc(combine_prefix(px)) + prefix <- combine_prefix(px, keep_mu = TRUE) + special <- attr(prior, "special")[[prefix]] if (!is.null(special[["hs_df"]])) { # data for the horseshoe prior hs_obj_names <- paste0("hs_", @@ -443,13 +449,13 @@ data_prior <- function(prior, data, nlpar = "") { } else { hs_data$hs_scale_global <- special$hs_par_ratio / sqrt(nrow(data)) } - names(hs_data) <- paste0(hs_obj_names, nlpar) + names(hs_data) <- paste0(hs_obj_names, p) out <- c(out, hs_data) } if (!is.null(special[["lasso_df"]])) { lasso_obj_names <- paste0("lasso_", c("df", "scale")) lasso_data <- special[lasso_obj_names] - names(lasso_data) <- paste0(lasso_obj_names, nlpar) + names(lasso_data) <- paste0(lasso_obj_names, p) out <- c(out, lasso_data) } out diff --git a/R/deprecated-methods.R b/R/deprecated-methods.R index da635a8fb..9b8b6c0d5 100644 --- a/R/deprecated-methods.R +++ b/R/deprecated-methods.R @@ -23,15 +23,15 @@ old_fixef_brmsfit <- function(object, estimate = "mean", ...) { old_ranef_brmsfit <- function(object, estimate = c("mean", "median"), var = FALSE, ...) { # ranef.brmsfit as implemented prior to brms 1.7.0 - .ranef <- function(group, nlpar = "") { + .ranef <- function(group, px = list()) { # get group-level effects of a grouping factor # Args: # group: name of a grouping factor # nlpar: name of a non-linear parameter - take <- object$ranef$group == group & object$ranef$nlpar == nlpar - rnames <- object$ranef[take, "coef"] - usc_nlpar <- usc(usc(nlpar)) - rpars <- pars[grepl(paste0("^r_", group, usc_nlpar, "\\["), pars)] + rnames <- subset2(object$ranef, group = group, ls = px)$coef + px <- check_prefix(px) + p <- usc(usc(combine_prefix(px))) + rpars <- pars[grepl(paste0("^r_", group, p, "\\["), pars)] if (!length(rpars)) { return(NULL) } @@ -71,21 +71,19 @@ old_ranef_brmsfit <- function(object, estimate = c("mean", "median"), attr(out, "var") <- Var } rownames(out) <- levels - if (nchar(nlpar)) { - attr(out, "nlpar") <- nlpar - } + attr(out, "prefix") <- px return(out) } stopifnot(is.brmsfit(object)) estimate <- match.arg(estimate) pars <- parnames(object) - group_nlpar <- unique(object$ranef[, c("group", "nlpar")]) - ranef <- named_list(group_nlpar$group) + sub_ranef <- unique(object$ranef[, c("group", vars_prefix())]) + ranef <- named_list(sub_ranef$group) for (i in seq_along(ranef)) { ranef[[i]] <- .ranef( - group = group_nlpar$group[i], - nlpar = group_nlpar$nlpar[i] + group = sub_ranef$group[i], + px = sub_ranef[i, vars_prefix()] ) } rmNULL(ranef) @@ -144,20 +142,22 @@ old_coef_brmsfit <- function(object, estimate = c("mean", "median"), ...) { estimate <- match.arg(estimate) fixef <- fixef(object, estimate = estimate, old = TRUE, ...) ranef <- ranef(object, estimate = estimate, old = TRUE, ...) - nlpars <- ulapply(ranef, attr, "nlpar") - if (length(nlpars)) { + px <- do.call(rbind, lapply(ranef, attr, "prefix")) + p_all <- combine_prefix(px) + if (length(px)) { # do not combine effects of different nlpars - unique_nlpars <- unique(nlpars) - coef <- named_list(unique_nlpars) - for (p in unique_nlpars) { - ranef_temp <- ranef[nlpars %in% p] - rx <- paste0("^", p, "_") + unique_px <- unique(px) + p <- combine_prefix(unique_px) + coef <- named_list(p) + for (j in seq_along(p)) { + ranef_temp <- ranef[p_all %in% p[j]] + rx <- paste0("^", usc(p[j], "suffix")) take_rows <- grepl(rx, rownames(fixef)) fixef_temp <- fixef[take_rows, , drop = FALSE] rownames(fixef_temp) <- sub(rx, "", rownames(fixef_temp)) - coef[[p]] <- .coef(ranef_temp, fixef_temp) - for (i in seq_along(coef[[p]])) { - attr(coef[[p]][[i]], "nlpar") <- p + coef[[p[j]]] <- .coef(ranef_temp, fixef_temp) + for (i in seq_along(coef[[p[j]]])) { + attr(coef[[p[j]]][[i]], "prefix") <- unique_px[j, ] } } coef <- unlist(unname(coef), recursive = FALSE) @@ -234,7 +234,7 @@ old_VarCorr_brmsfit <- function(x, estimate = "mean", ...) { get_names <- function(group) { # get names of group-level parameters r <- x$ranef[x$ranef$group == group, ] - rnames <- paste0(usc(r$nlpar, "suffix"), r$coef) + rnames <- paste0(usc(combine_prefix(r), "suffix"), r$coef) cor_type <- paste0("cor_", group) sd_pars <- paste0("sd_", group, "__", rnames) cor_pars <- get_cornames(rnames, type = cor_type, brackets = FALSE) @@ -247,7 +247,7 @@ old_VarCorr_brmsfit <- function(x, estimate = "mean", ...) { } # special treatment of residuals variances in linear models has_sigma <- has_sigma(family, bterms, incmv = TRUE) - if (has_sigma && !"sigma" %in% names(bterms$auxpars)) { + if (has_sigma && !"sigma" %in% names(bterms$dpars)) { cor_pars <- get_cornames( bterms$response, type = "rescor", brackets = FALSE ) diff --git a/R/extract_draws.R b/R/extract_draws.R index 891f40bf7..0994d49d2 100644 --- a/R/extract_draws.R +++ b/R/extract_draws.R @@ -25,9 +25,9 @@ extract_draws.brmsfit <- function(x, newdata = NULL, re_formula = NULL, data = do.call(amend_newdata, c(newd_args, list(...))) ) args <- c(newd_args, nlist( - sample_new_levels, subset, nsamples, nug, C = draws$data[["C"]] + sample_new_levels, subset, nsamples, nug )) - keep <- !grepl("^(X|Z|J|C)", names(draws$data)) + keep <- !grepl("^(X|Z|J)", names(draws$data)) draws$data <- subset_attr(draws$data, keep) # special treatment of multivariate models @@ -36,26 +36,30 @@ extract_draws.brmsfit <- function(x, newdata = NULL, re_formula = NULL, draws$mu[["mv"]] <- named_list(resp) for (j in seq_along(resp)) { r <- resp[j] - more_args <- list(x = bterms$auxpars[["mu"]], nlpar = r, mv = TRUE) + more_args <- list(x = bterms$dpars[["mu"]], mv = TRUE) + more_args$x$resp <- r draws$mu[["mv"]][[r]] <- do.call(extract_draws, c(args, more_args)) - draws$mu[["mv"]][[r]][["f"]] <- bterms$auxpars[["mu"]]$family + draws$mu[["mv"]][[r]][["f"]] <- bterms$dpars[["mu"]]$family if (isTRUE(ncol(draws$mu[["mv"]][[r]]$data$Y) > 1L)) { draws$mu[["mv"]][[r]]$data$Y <- draws$mu[["mv"]][[r]]$data$Y[, j] } } - bterms$auxpars[["mu"]] <- NULL + bterms$dpars[["mu"]] <- NULL } # extract draws of auxiliary parameters am_args <- nlist(x, subset) - valid_auxpars <- valid_auxpars(family(x), bterms = bterms) - for (ap in valid_auxpars) { + valid_dpars <- valid_dpars(family(x), bterms = bterms) + for (ap in valid_dpars) { ap_regex <- paste0("^", ap, "($|_)") - if (is.btl(bterms$auxpars[[ap]]) || is.btnl(bterms$auxpars[[ap]])) { - more_args <- nlist(x = bterms$auxpars[[ap]], nlpar = ap) + if (is.btl(bterms$dpars[[ap]]) || is.btnl(bterms$dpars[[ap]])) { + p <- usc(combine_prefix(bterms$dpars[[ap]])) + more_args <- nlist( + x = bterms$dpars[[ap]], C = draws$data[[paste0("C", p)]] + ) draws[[ap]] <- do.call(extract_draws, c(args, more_args)) - draws[[ap]][["f"]] <- bterms$auxpars[[ap]]$family - } else if (is.numeric(bterms$fauxpars[[ap]])) { - draws[[ap]] <- bterms$fauxpars[[ap]] + draws[[ap]][["f"]] <- bterms$dpars[[ap]]$family + } else if (is.numeric(bterms$fdpars[[ap]]$value)) { + draws[[ap]] <- bterms$fdpars[[ap]]$value } else if (any(grepl(ap_regex, parnames(x)))) { draws[[ap]] <- do.call(as.matrix, c(am_args, pars = ap_regex)) } @@ -100,16 +104,15 @@ extract_draws.brmsfit <- function(x, newdata = NULL, re_formula = NULL, } #' @export -extract_draws.btnl <- function(x, C, nlpar = "", ...) { +extract_draws.btnl <- function(x, C, ...) { # Args: # C: matrix containing covariates - # nlpar: unused and should NOT be passed further # ...: passed to extract_draws.btl nlpars <- names(x$nlpars) draws <- named_list(nlpars) for (nlp in nlpars) { rhs_formula <- x$nlpars[[nlp]]$formula - draws$nlpars[[nlp]] <- extract_draws(x$nlpars[[nlp]], nlpar = nlp, ...) + draws$nlpars[[nlp]] <- extract_draws(x$nlpars[[nlp]], ...) } nsamples <- draws$nlpars[[nlpars[1]]]$nsamples stopifnot(is.matrix(C)) @@ -127,12 +130,11 @@ extract_draws.btl <- function(x, fit, newdata = NULL, re_formula = NULL, allow_new_levels = FALSE, sample_new_levels = "uncertainty", new_objects = list(), incl_autocor = TRUE, - subset = NULL, nlpar = "", smooths_only = FALSE, + subset = NULL, smooths_only = FALSE, mv = FALSE, nug = NULL, ...) { # extract draws of all kinds of effects # Args: # fit: a brmsfit object - # nlpar: name of a non-linear parameter # mv: is the model multivariate? # ...: further elements to store in draws # other arguments: see extract_draws.brmsfit @@ -141,23 +143,22 @@ extract_draws.btl <- function(x, fit, newdata = NULL, re_formula = NULL, dots <- list(...) stopifnot(is.brmsfit(fit), is.brmsformula(fit$formula)) fit <- remove_autocor(fit, incl_autocor) - nlpar <- check_nlpar(nlpar) + px <- check_prefix(x) + p <- combine_prefix(px) nsamples <- nsamples(fit, subset = subset) fit$formula$formula <- update(fit$formula$formula, rhs(x$formula)) # ensure that auxiliary parameters are not included (fixes #154) fit$formula$pforms <- fit$formula$pfix <- NULL fit$formula$nl <- FALSE fit$ranef <- tidy_ranef(parse_bf(fit$formula), data = fit$data) - if (nzchar(nlpar)) { + if (nzchar(p)) { # make sure not to evaluate family specific stuff fit$formula[["response"]] <- NA - fit$family <- fit$formula$family <- par_family() + fit$family <- fit$formula$family <- .dpar_family() } new_formula <- update_re_terms(fit$formula, re_formula = re_formula) bterms <- parse_bf(new_formula, family = family(fit)) new_ranef <- tidy_ranef(bterms, model.frame(fit)) - nlpar_usc <- usc(nlpar, "suffix") - usc_nlpar <- usc(usc(nlpar)) newd_args <- nlist( fit, newdata, re_formula, allow_new_levels, new_objects, check_response = FALSE @@ -178,29 +179,30 @@ extract_draws.btl <- function(x, fit, newdata = NULL, re_formula = NULL, args <- nlist(x = fit, subset) new <- !is.null(newdata) + # deviate from the usual way of passing bterms and data + # as bterms and px do not coincide in extract_draws.btl fixef <- colnames(draws$data[["X"]]) monef <- colnames(draws$data[["Xmo"]]) csef <- colnames(draws$data[["Xcs"]]) - meef <- get_me_labels(bterms$auxpars$mu, fit$data) - smooths <- get_sm_labels(bterms$auxpars$mu, fit$data, covars = TRUE) - gpef <- get_gp_labels(bterms$auxpars$mu, fit$data, covars = TRUE) - sdata_old <- NULL - if (length(gpef) && new) { - oldd_args <- newd_args[!names(newd_args) %in% "newdata"] - sdata_old <- do.call(amend_newdata, c(oldd_args, list(newdata = NULL))) - } + meef <- get_me_labels(bterms$dpars$mu, fit$data) + smooths <- get_sm_labels(bterms$dpars$mu, fit$data, covars = TRUE) + gpef <- get_gp_labels(bterms$dpars$mu, fit$data, covars = TRUE) draws <- c(draws, - extract_draws_fe(fixef, args, nlpar = nlpar, old_cat = draws$old_cat), - extract_draws_mo(monef, args, sdata = draws$data, nlpar = nlpar), - extract_draws_cs(csef, args, nlpar = nlpar, old_cat = draws$old_cat), - extract_draws_me(meef, args, sdata = draws$data, nlpar = nlpar, new = new), - extract_draws_sm(smooths, args, sdata = draws$data, nlpar = nlpar), - extract_draws_gp(gpef, args, sdata = draws$data, sdata_old = sdata_old, - nlpar = nlpar, new = new, nug = nug), - extract_draws_re(new_ranef, args, sdata = draws$data, nlpar = nlpar, - sample_new_levels = sample_new_levels) + extract_draws_fe(fixef, args, px = px, old_cat = draws$old_cat), + extract_draws_mo(monef, args, sdata = draws$data, px = px), + extract_draws_cs(csef, args, px = px, old_cat = draws$old_cat), + extract_draws_me(meef, args, sdata = draws$data, px = px, new = new), + extract_draws_sm(smooths, args, sdata = draws$data, px = px), + extract_draws_gp( + gpef, args, sdata = draws$data, px = px, + new = new, nug = nug, newd_args = newd_args + ), + extract_draws_re( + new_ranef, args, sdata = draws$data, px = px, + sample_new_levels = sample_new_levels + ) ) - if (!use_cov(fit$autocor) && (!nzchar(nlpar) || mv)) { + if (!use_cov(fit$autocor) && (!nzchar(p) || mv)) { # only include autocorrelation parameters in draws for mu draws <- c(draws, extract_draws_autocor( @@ -211,7 +213,7 @@ extract_draws.btl <- function(x, fit, newdata = NULL, re_formula = NULL, structure(draws, predicted = TRUE) } -extract_draws_fe <- function(fixef, args, nlpar = "", old_cat = 0L) { +extract_draws_fe <- function(fixef, args, px = list(), old_cat = 0L) { # extract draws of ordinary population-level effects # Args: # fixef: names of the population-level effects @@ -221,10 +223,10 @@ extract_draws_fe <- function(fixef, args, nlpar = "", old_cat = 0L) { # Returns: # A named list to be interpreted by linear_predictor stopifnot("x" %in% names(args)) - nlpar_usc <- usc(nlpar, "suffix") + p <- usc(combine_prefix(px), "suffix") draws <- list() if (length(fixef) && old_cat != 1L) { - b_pars <- paste0("b_", nlpar_usc, fixef) + b_pars <- paste0("b_", p, fixef) draws[["b"]] <- do.call(as.matrix, c(args, list(pars = b_pars, exact = TRUE)) ) @@ -232,7 +234,7 @@ extract_draws_fe <- function(fixef, args, nlpar = "", old_cat = 0L) { draws } -extract_draws_mo <- function(monef, args, sdata, nlpar = "") { +extract_draws_mo <- function(monef, args, sdata, px = list()) { # extract draws of monotonic effects # Args: # monef: names of the monotonic effects @@ -242,7 +244,7 @@ extract_draws_mo <- function(monef, args, sdata, nlpar = "") { # Returns: # A named list to be interpreted by linear_predictor stopifnot("x" %in% names(args)) - nlpar_usc <- usc(nlpar, "suffix") + p <- usc(combine_prefix(px), "suffix") draws <- list() if (length(monef)) { draws[["bmo"]] <- draws$simplex <- named_list(monef) @@ -252,11 +254,11 @@ extract_draws_mo <- function(monef, args, sdata, nlpar = "") { ifelse(any(grepl("^bm_", parnames(args$x))), "bm_", "b_") ) for (i in seq_along(monef)) { - bmo_par <- paste0(bmo, nlpar_usc, monef[i]) + bmo_par <- paste0(bmo, p, monef[i]) draws[["bmo"]][[i]] <- do.call(as.matrix, c(args, list(pars = bmo_par, exact = TRUE)) ) - simplex_par <- paste0("simplex_", nlpar_usc, monef[i], + simplex_par <- paste0("simplex_", p, monef[i], "[", seq_len(sdata$Jm[i]), "]") draws[["simplex"]][[i]] <- do.call(as.matrix, c(args, list(pars = simplex_par, exact = TRUE)) @@ -266,7 +268,7 @@ extract_draws_mo <- function(monef, args, sdata, nlpar = "") { draws } -extract_draws_cs <- function(csef, args, nlpar = "", old_cat = 0L) { +extract_draws_cs <- function(csef, args, px = list(), old_cat = 0L) { # category specific effects # extract draws of category specific effects # Args: @@ -277,7 +279,7 @@ extract_draws_cs <- function(csef, args, nlpar = "", old_cat = 0L) { # Returns: # A named list to be interpreted by linear_predictor stopifnot("x" %in% names(args)) - nlpar_usc <- usc(nlpar, "suffix") + p <- usc(combine_prefix(px), "suffix") draws <- list() if (is_ordinal(family(args$x))) { draws[["Intercept"]] <- do.call(as.matrix, @@ -286,7 +288,7 @@ extract_draws_cs <- function(csef, args, nlpar = "", old_cat = 0L) { if (length(csef)) { # as of brms > 1.0.1 the original prefix 'bcs' is used bcs <- ifelse(any(grepl("^bcs_", parnames(args$x))), "^bcs_", "^b_") - cs_pars <- paste0(bcs, nlpar_usc, csef, "\\[") + cs_pars <- paste0(bcs, p, csef, "\\[") draws[["cs"]] <- do.call(as.matrix, c(args, list(pars = cs_pars))) } } else if (old_cat == 1L) { @@ -296,7 +298,7 @@ extract_draws_cs <- function(csef, args, nlpar = "", old_cat = 0L) { draws } -extract_draws_me <- function(meef, args, sdata, nlpar = "", +extract_draws_me <- function(meef, args, sdata, px = list(), new = FALSE) { # extract draws of noise-free effects # Args: @@ -308,14 +310,14 @@ extract_draws_me <- function(meef, args, sdata, nlpar = "", # Returns: # A named list to be interpreted by linear_predictor stopifnot("x" %in% names(args)) - nlpar_usc <- usc(nlpar, "suffix") + p <- usc(combine_prefix(px), "suffix") draws <- list() if (length(meef)) { if (new) { stop2("Predictions with noise-free variables are not yet ", "possible when passing new data.") } - if (!any(grepl(paste0("Xme_", nlpar_usc), parnames(args$x)))) { + if (!any(grepl(paste0("Xme_", p), parnames(args$x)))) { stop2("Noise-free variables were not saved. Please set ", "argument 'save_mevars' to TRUE when calling 'brm'.") } @@ -340,7 +342,7 @@ extract_draws_me <- function(meef, args, sdata, nlpar = "", draws[["bme"]] <- named_list(meef) attr(draws[["bme"]], "calls") <- parse(text = meef_terms) for (j in seq_along(draws[["bme"]])) { - bme_par <- paste0("bme_", nlpar_usc, meef[j]) + bme_par <- paste0("bme_", p, meef[j]) draws[["bme"]][[j]] <- do.call(as.matrix, c(args, list(pars = bme_par, exact = TRUE)) ) @@ -350,7 +352,7 @@ extract_draws_me <- function(meef, args, sdata, nlpar = "", draws[["Xme"]] <- named_list(uni_me) for (j in seq_along(draws[["Xme"]])) { Xme_pars <- paste0( - "Xme_", nlpar_usc, uni_me[j], "[", seq_len(nobs(args$x)), "]" + "Xme_", p, uni_me[j], "[", seq_len(nobs(args$x)), "]" ) draws[["Xme"]][[j]] <- do.call(as.matrix, c(args, list(pars = Xme_pars, exact = TRUE)) @@ -368,7 +370,7 @@ extract_draws_me <- function(meef, args, sdata, nlpar = "", draws } -extract_draws_sm <- function(smooths, args, sdata, nlpar = "") { +extract_draws_sm <- function(smooths, args, sdata, px = list()) { # extract draws of smooth terms # Args: # smooths: names of the smooth terms @@ -378,7 +380,7 @@ extract_draws_sm <- function(smooths, args, sdata, nlpar = "") { # Returns: # A named list to be interpreted by linear_predictor stopifnot("x" %in% names(args)) - nlpar_usc <- usc(nlpar, "suffix") + p <- usc(combine_prefix(px), "suffix") draws <- list() if (length(smooths)) { draws[["Zs"]] <- draws[["s"]] <- named_list(smooths) @@ -386,7 +388,7 @@ extract_draws_sm <- function(smooths, args, sdata, nlpar = "") { nb <- seq_len(attr(smooths, "nbases")[[i]]) for (j in nb) { draws[["Zs"]][[smooths[i]]][[j]] <- sdata[[paste0("Zs_", i, "_", j)]] - s_pars <- paste0("^s_", nlpar_usc, smooths[i], "_", j, "\\[") + s_pars <- paste0("^s_", p, smooths[i], "_", j, "\\[") draws[["s"]][[smooths[i]]][[j]] <- do.call(as.matrix, c(args, list(pars = s_pars))) } @@ -395,13 +397,19 @@ extract_draws_sm <- function(smooths, args, sdata, nlpar = "") { draws } -extract_draws_gp <- function(gpef, args, sdata, sdata_old = NULL, - nlpar = "", new = FALSE, nug = NULL) { +extract_draws_gp <- function(gpef, args, sdata, px = list(), + new = FALSE, nug = NULL, + newd_args = list()) { # extract draws for gaussian processes # Args: # gpef: names of the gaussian process terms stopifnot("x" %in% names(args)) - nlpar_usc <- usc(nlpar, "suffix") + sdata_old <- NULL + if (length(gpef) && new) { + newd_args["newdata"] <- list(NULL) + sdata_old <- do.call(amend_newdata, newd_args) + } + p <- usc(combine_prefix(px), "suffix") draws <- list() if (length(gpef)) { draws[["gp"]] <- named_list(gpef) @@ -412,11 +420,11 @@ extract_draws_gp <- function(gpef, args, sdata, sdata_old = NULL, gp <- list() by_levels <- attr(gpef, "by_levels")[[i]] gp_names <- paste0(gpef[i], usc(by_levels)) - sdgp <- paste0("^sdgp_", nlpar_usc, gp_names) + sdgp <- paste0("^sdgp_", p, gp_names) gp[["sdgp"]] <- do.call(as.matrix, c(args, list(pars = sdgp))) - lscale <- paste0("^lscale_", nlpar_usc, gp_names) + lscale <- paste0("^lscale_", p, gp_names) gp[["lscale"]] <- do.call(as.matrix, c(args, list(pars = lscale))) - zgp <- paste0("^zgp_", nlpar_usc, gpef[i], "\\[") + zgp <- paste0("^zgp_", p, gpef[i], "\\[") gp[["zgp"]] <- do.call(as.matrix, c(args, list(pars = zgp))) gp[["bynum"]] <- sdata[[paste0("Cgp_", i)]] Jgp_pos <- grepl(paste0("^Jgp_", i), names(sdata)) @@ -442,7 +450,7 @@ extract_draws_gp <- function(gpef, args, sdata, sdata_old = NULL, draws } -extract_draws_re <- function(ranef, args, sdata, nlpar = "", +extract_draws_re <- function(ranef, args, sdata, px = list(), sample_new_levels = "uncertainty") { # extract draws of group-level effects # Args: @@ -454,14 +462,14 @@ extract_draws_re <- function(ranef, args, sdata, nlpar = "", # Returns: # A named list to be interpreted by linear_predictor stopifnot("x" %in% names(args)) - usc_nlpar <- usc(usc(nlpar)) + p <- combine_prefix(px) draws <- list() groups <- unique(ranef$group) # assigning S4 objects requires initialisation of list elements draws[c("Z", "Zmo", "Zcs", "Zme")] <- list(named_list(groups)) for (g in groups) { - new_r <- ranef[ranef$group == g, ] - r_pars <- paste0("^r_", g, usc_nlpar, "\\[") + new_r <- subset2(ranef, group = g) + r_pars <- paste0("^r_", g, usc(usc(p)), "\\[") r <- do.call(as.matrix, c(args, list(pars = r_pars))) if (is.null(r)) { stop2("Group-level effects for each level of group ", @@ -469,7 +477,7 @@ extract_draws_re <- function(ranef, args, sdata, nlpar = "", "when calling brm.") } nlevels <- ngrps(args$x)[[g]] - old_r <- args$x$ranef[args$x$ranef$group == g, ] + old_r <- subset2(args$x$ranef, group = g) used_re <- match(new_r$coef, old_r$coef) used_re_pars <- outer(seq_len(nlevels), (used_re - 1) * nlevels, "+") used_re_pars <- as.vector(used_re_pars) @@ -506,12 +514,12 @@ extract_draws_re <- function(ranef, args, sdata, nlpar = "", } } else if (sample_new_levels == "gaussian") { # extract hyperparameters used to compute the covariance matrix - sd_pars <- paste0("sd_", g, usc_nlpar, "__", new_r$coef) + sd_pars <- paste0("sd_", g, usc(usc(p)), "__", new_r$coef) sd_samples <- do.call(as.matrix, c(args, list(pars = sd_pars, exact_match = TRUE)) ) cor_type <- paste0("cor_", g) - cor_pars <- paste0(usc(nlpar, "suffix"), new_r$coef) + cor_pars <- paste0(usc(p, "suffix"), new_r$coef) cor_pars <- get_cornames(cor_pars, type = cor_type, brackets = FALSE) cor_samples <- matrix( 0, nrow = nrow(sd_samples), ncol = length(cor_pars) @@ -559,7 +567,7 @@ extract_draws_re <- function(ranef, args, sdata, nlpar = "", levels <- unique(unlist(gf)) r <- subset_levels(r, levels, nranef) # monotonic group-level terms - new_r_mo <- new_r[new_r$type == "mo", ] + new_r_mo <- subset2(new_r, type = "mo") if (nrow(new_r_mo)) { Z <- matrix(1, length(gf[[1]])) draws[["Zmo"]][[g]] <- prepare_Z(Z, gf, max_level, weights) @@ -571,7 +579,7 @@ extract_draws_re <- function(ranef, args, sdata, nlpar = "", } } # noise-free group-level terms - new_r_me <- new_r[new_r$type == "me", ] + new_r_me <- subset2(new_r, type = "me") if (nrow(new_r_me)) { Z <- matrix(1, length(gf[[1]])) draws[["Zme"]][[g]] <- prepare_Z(Z, gf, max_level, weights) @@ -583,7 +591,7 @@ extract_draws_re <- function(ranef, args, sdata, nlpar = "", } } # category specific group-level terms - new_r_cs <- new_r[new_r$type == "cs", ] + new_r_cs <- subset2(new_r, type = "cs") if (nrow(new_r_cs)) { Z <- do.call(cbind, lapply(unique(new_r_cs$gn), function(k) sdata[[paste0("Z_", k)]] @@ -598,7 +606,7 @@ extract_draws_re <- function(ranef, args, sdata, nlpar = "", } } # basic group-level effects - new_r_basic <- new_r[!nzchar(new_r$type), ] + new_r_basic <- subset2(new_r, type = "") if (nrow(new_r_basic)) { Z <- do.call(cbind, lapply(unique(new_r_basic$gn), function(k) sdata[[paste0("Z_", k)]] @@ -668,7 +676,8 @@ subset_levels <- function(x, levels, nranef) { # levels: grouping factor levels to keep # nranef: number of group-level effects take_levels <- ulapply(levels, - function(l) ((l - 1) * nranef + 1):(l * nranef)) + function(l) ((l - 1) * nranef + 1):(l * nranef) + ) x[, take_levels, drop = FALSE] } @@ -682,8 +691,10 @@ prepare_Z <- function(Z, gf, max_level = max(unlist(gf)), # weights: optional list of weights of the same length as gf nranef <- ncol(Z) levels <- unique(unlist(gf)) - Z <- mapply(expand_matrix, x = gf, weights = weights, - MoreArgs = nlist(A = Z, max_level)) + Z <- mapply( + expand_matrix, x = gf, weights = weights, + MoreArgs = nlist(A = Z, max_level) + ) Z <- Reduce("+", Z) subset_levels(Z, levels, nranef) } diff --git a/R/families.R b/R/families.R index 9c8f51e5e..621ad0d11 100644 --- a/R/families.R +++ b/R/families.R @@ -276,19 +276,19 @@ brmsfamily <- function(family, link = NULL, link_sigma = "log", list(family = family, link = slink), class = c("brmsfamily", "family") ) - for (ap in valid_auxpars(out$family)) { - alink <- as.character(aux_links[[paste0("link_", ap)]]) + for (dp in valid_dpars(out$family)) { + alink <- as.character(aux_links[[paste0("link_", dp)]]) if (length(alink)) { if (length(alink) > 1L) { stop2("Link functions must be of length 1.") } - valid_links <- links_auxpars(ap) + valid_links <- links_dpars(dp) if (!alink %in% valid_links) { stop2("'", alink, "' is not a supported link ", - "for parameter '", ap, "'.\nSupported links are: ", + "for parameter '", dp, "'.\nSupported links are: ", collapse_comma(valid_links)) } - out[[paste0("link_", ap)]] <- alink + out[[paste0("link_", dp)]] <- alink } } if (is_ordinal(out$family)) { @@ -539,33 +539,6 @@ acat <- function(link = "logit", link_disc = "log", link_disc = link_disc, threshold = threshold) } -par_family <- function(par = NULL, link = NULL) { - # set up special family objects for parameters - # Args: - # par: name of the parameter - # link: optional link function of the parameter - ap_class <- auxpar_class(par) - if (!is.null(par) && !isTRUE(ap_class %in% auxpars())) { - stop2("Parameter '", par, "' is invalid.") - } - if (is.null(par)) { - link <- "identity" - } else { - links <- links_auxpars(ap_class) - if (is.null(link)) { - link <- links[1] - } else { - if (!isTRUE(link %in% links)) { - stop2("Link '", link, "' is invalid for parameter '", par, "'.") - } - } - } - structure( - list(family = "", link = link, par = par), - class = c("brmsfamily", "family") - ) -} - check_family <- function(family, link = NULL, threshold = NULL) { # checks and corrects validity of the model family # Args: @@ -765,22 +738,59 @@ family_names.mixfamily <- function(family, ...) { } #' @export -auxpar_family.default <- function(family, auxpar, ...) { - ap_class <- auxpar_class(auxpar) - if (!identical(ap_class, "mu")) { - link <- family[[paste0("link_", ap_class)]] - family <- par_family(auxpar, link) +family_names.brmsformula <- function(family, ...) { + family_names(family$family) +} + +#' @export +family_names.brmsterms <- function(family, ...) { + family_names(family$family) +} + +#' @export +dpar_family.default <- function(family, dpar, ...) { + dp_class <- dpar_class(dpar) + if (!identical(dp_class, "mu")) { + link <- family[[paste0("link_", dp_class)]] + family <- .dpar_family(dpar, link) } family } #' @export -auxpar_family.mixfamily <- function(family, auxpar, ...) { - ap_id <- as.numeric(auxpar_id(auxpar)) - if (!(length(ap_id) == 1L && is.numeric(ap_id))) { - stop2("Parameter '", auxpar, "' is not a valid mixture parameter.") +dpar_family.mixfamily <- function(family, dpar, ...) { + dp_id <- as.numeric(dpar_id(dpar)) + if (!(length(dp_id) == 1L && is.numeric(dp_id))) { + stop2("Parameter '", dpar, "' is not a valid mixture parameter.") + } + dpar_family(family$mix[[dp_id]], dpar, ...) +} + +.dpar_family <- function(dpar = NULL, link = NULL) { + # set up special family objects for distributional parameters + # Args: + # dpar: name of the distributional parameter + # link: optional link function of the parameter + dp_class <- dpar_class(dpar) + if (!is.null(dpar) && !isTRUE(dp_class %in% dpars())) { + stop2("Parameter '", dpar, "' is invalid.") + } + if (is.null(dpar)) { + link <- "identity" + } else { + links <- links_dpars(dp_class) + if (is.null(link)) { + link <- links[1] + } else { + if (!isTRUE(link %in% links)) { + stop2("Link '", link, "' is invalid for parameter '", dpar, "'.") + } + } } - auxpar_family(family$mix[[ap_id]], auxpar, ...) + structure( + nlist(family = "", link, dpar), + class = c("brmsfamily", "family") + ) } #' @export @@ -794,13 +804,13 @@ print.brmsfamily <- function(x, links = FALSE, newline = TRUE, ...) { cat("Threshold:", x$threshold, "\n") } if (isTRUE(links) || is.character(links)) { - ap_links <- x[grepl("^link_", names(x))] - names(ap_links) <- sub("^link_", "", names(ap_links)) + dp_links <- x[grepl("^link_", names(x))] + names(dp_links) <- sub("^link_", "", names(dp_links)) if (is.character(links)) { - ap_links <- rmNULL(ap_links[links]) + dp_links <- rmNULL(dp_links[links]) } - for (ap in names(ap_links)) { - cat(paste0("Link function of ", ap, ": ", ap_links[[ap]], " \n")) + for (dp in names(dp_links)) { + cat(paste0("Link function of ", dp, ": ", dp_links[[dp]], " \n")) } } if (newline) { diff --git a/R/fitted.R b/R/fitted.R index 436944a22..843ad6789 100644 --- a/R/fitted.R +++ b/R/fitted.R @@ -92,7 +92,7 @@ fitted_gen_extreme_value <- function(draws) { draws$sigma <- get_sigma( draws$sigma, data = draws$data, dim = dim_mu(draws) ) - draws$xi <- get_auxpar(draws$xi) + draws$xi <- get_dpar(draws$xi) draws$mu <- ilink(draws$mu, draws$f$link) if (!is_trunc(draws$data)) { draws$mu <- with(draws, mu + sigma * (gamma(1 - xi) - 1) / xi) @@ -107,16 +107,16 @@ fitted_inverse.gaussian <- function(draws) { } fitted_exgaussian <- function(draws) { - draws$mu <- ilink(draws$mu, draws$f$link) + get_auxpar(draws$beta) + draws$mu <- ilink(draws$mu, draws$f$link) + get_dpar(draws$beta) fitted_trunc(draws) } fitted_wiener <- function(draws) { # mu is the drift rate draws$mu <- ilink(draws$mu, draws$f$link) - draws$bs <- get_auxpar(draws$bs) - draws$ndt <- get_auxpar(draws$ndt) - draws$bias <- get_auxpar(draws$bias) + draws$bs <- get_dpar(draws$bs) + draws$ndt <- get_dpar(draws$ndt) + draws$bias <- get_dpar(draws$bias) with(draws, ndt - bias / mu + bs / mu * (exp(- 2 * mu * bias) - 1) / (exp(-2 * mu * bs) - 1) @@ -132,7 +132,7 @@ fitted_von_mises <- function(draws) { } fitted_asym_laplace <- function(draws) { - draws$quantile <- get_auxpar(draws$quantile) + draws$quantile <- get_dpar(draws$quantile) draws$sigma <- get_sigma( draws$sigma, data = draws$data, dim = dim_mu(draws) ) @@ -208,8 +208,8 @@ fitted_zero_inflated_beta <- function(draws) { } fitted_zero_one_inflated_beta <- function(draws) { - draws$zoi <- get_auxpar(draws$zoi) - draws$coi <- get_auxpar(draws$coi) + draws$zoi <- get_dpar(draws$zoi) + draws$coi <- get_dpar(draws$coi) draws$zoi * draws$coi + ilink(draws$mu, draws$f$link) * (1 - draws$zoi) } @@ -249,13 +249,13 @@ fitted_mixture <- function(draws) { for (j in seq_along(families)) { fitted_fun <- paste0("fitted_", families[j]) fitted_fun <- get(fitted_fun, asNamespace("brms")) - auxpars <- valid_auxpars(families[j]) + dpars <- valid_dpars(families[j]) tmp_draws <- list( f = draws$f$mix[[j]], nsamples = draws[["nsamples"]], data = draws[["data"]] ) - for (ap in auxpars) { + for (ap in dpars) { tmp_draws[[ap]] <- draws[[paste0(ap, j)]] } if (length(dim(draws$theta)) == 3L) { @@ -372,7 +372,7 @@ fitted_trunc_student <- function(draws, lb, ub) { draws$sigma <- get_sigma( draws$sigma, data = draws$data, dim = dim_mu(draws) ) - draws$nu <- get_auxpar(draws$nu) + draws$nu <- get_dpar(draws$nu) zlb <- (lb - draws$mu) / draws$sigma zub <- (ub - draws$mu) / draws$sigma # see Kim 2008: Moments of truncated Student-t distribution diff --git a/R/formula-helpers.R b/R/formula-helpers.R index f51b0e87f..786c0a800 100644 --- a/R/formula-helpers.R +++ b/R/formula-helpers.R @@ -611,3 +611,19 @@ formula2str <- function(formula, rm = c(0, 0), space = c("rm", "trim")) { is.formula <- function(x) { inherits(x, "formula") } + +expand_dot_formula <- function(formula, data = NULL) { + # expand the '.' variable in formula using stats::terms + if (isTRUE("." %in% all.vars(formula))) { + att <- attributes(formula) + try_terms <- try( + stats::terms(formula, data = data), + silent = TRUE + ) + if (!is(try_terms, "try-error")) { + formula <- formula(try_terms) + } + attributes(formula) <- att + } + formula +} diff --git a/R/generics.R b/R/generics.R index dcf7c591c..4a4f6046b 100644 --- a/R/generics.R +++ b/R/generics.R @@ -1204,8 +1204,14 @@ extract_draws <- function(x, ...) { UseMethod("extract_draws") } +make_Jmo_list <- function(x, data, ...) { + # compute Jmo values based on the original data + # as the basis for doing predictions with new data + UseMethod("make_Jmo_list") +} + make_smooth_list <- function(x, data, ...) { - # compute smoothing objects based on the original data + # compute smooth objects based on the original data # as the basis for doing predictions with new data UseMethod("make_smooth_list") } @@ -1221,9 +1227,9 @@ check_prior_special <- function(x, ...) { UseMethod("check_prior_special") } -auxpar_family <- function(family, auxpar, ...) { +dpar_family <- function(family, dpar, ...) { # generate a family object of an auxiliary parameter - UseMethod("auxpar_family") + UseMethod("dpar_family") } family_names <- function(family, ...) { @@ -1231,9 +1237,9 @@ family_names <- function(family, ...) { UseMethod("family_names") } -valid_auxpars <- function(family, ...) { +valid_dpars <- function(family, ...) { # get valid auxiliary parameters for a family - UseMethod("valid_auxpars") + UseMethod("valid_dpars") } stan_llh <- function(family, ...) { diff --git a/R/loglik.R b/R/loglik.R index 8454cb3a5..65a1784e5 100644 --- a/R/loglik.R +++ b/R/loglik.R @@ -20,7 +20,7 @@ loglik_gaussian <- function(i, draws, data = data.frame()) { } loglik_student <- function(i, draws, data = data.frame()) { - args <- list(df = get_auxpar(draws$nu, i = i), + args <- list(df = get_dpar(draws$nu, i = i), mu = ilink(get_eta(draws$mu, i), draws$f$link), sigma = get_sigma(draws$sigma, data = draws$data, i = i)) out <- loglik_censor(dist = "student_t", args = args, @@ -51,7 +51,7 @@ loglik_lognormal <- function(i, draws, data = data.frame()) { loglik_skew_normal <- function(i, draws, data = data.frame()) { sigma <- get_sigma(draws$sigma, data = draws$data, i = i) - alpha <- get_auxpar(draws$alpha, i = i) + alpha <- get_dpar(draws$alpha, i = i) mu <- ilink(get_eta(draws$mu, i), draws$f$link) args <- nlist(mu, sigma, alpha) out <- loglik_censor( @@ -84,7 +84,7 @@ loglik_student_mv <- function(i, draws, data = data.frame()) { } else { mu <- ilink(get_eta(draws$mu, i), draws$f$link)[, 1, ] } - nu <- get_auxpar(draws$nu, obs) + nu <- get_dpar(draws$nu, obs) out <- sapply(1:draws$nsamples, function(s) dmulti_student_t(draws$data$Y[i, ], df = nu[s, ], mu = mu[s, ], Sigma = draws$Sigma[s, , ], log = TRUE)) @@ -154,7 +154,7 @@ loglik_student_cov <- function(i, draws, data = data.frame()) { Sigma <- do.call(get_cov_matrix_arma1, args) } mu <- ilink(get_eta(draws$mu, obs), draws$f$link) - nu <- get_auxpar(draws$nu, obs) + nu <- get_dpar(draws$nu, obs) out <- sapply(1:draws$nsamples, function(s) dmulti_student_t( draws$data$Y[obs], df = nu[s, ], mu = mu[s, ], @@ -198,7 +198,7 @@ loglik_student_lagsar <- function(i, draws, data = data.frame()) { N <- draws$data$N mu <- ilink(get_eta(draws$mu), draws$f$link) sigma <- get_sigma(draws$sigma, data = draws$data) - nu <- get_auxpar(draws$nu) + nu <- get_dpar(draws$nu) # weights, truncation and censoring not yet allowed sapply(1:draws$nsamples, .loglik_student_lagsar) } @@ -228,7 +228,7 @@ loglik_student_errorsar <- function(i, draws, data = data.frame()) { } mu <- ilink(get_eta(draws$mu), draws$f$link) sigma <- get_sigma(draws$sigma, data = draws$data) - nu <- get_auxpar(draws$nu) + nu <- get_dpar(draws$nu) # weights, truncation and censoring not yet allowed sapply(1:draws$nsamples, .loglik_student_errorsar) } @@ -244,7 +244,7 @@ loglik_gaussian_fixed <- function(i, draws, data = data.frame()) { loglik_student_fixed <- function(i, draws, data = data.frame()) { stopifnot(i == 1) mu <- ilink(get_eta(draws$mu, 1:nrow(draws$data$V)), draws$f$link) - nu <- get_auxpar(draws$nu, 1:nrow(draws$data$V)) + nu <- get_dpar(draws$nu, 1:nrow(draws$data$V)) sapply(1:draws$nsamples, function(s) dmulti_student_t(draws$data$Y, df = nu[s, ], mu = mu[s, ], Sigma = draws$data$V, log = TRUE)) @@ -332,7 +332,7 @@ loglik_weibull <- function(i, draws, data = data.frame()) { } loglik_frechet <- function(i, draws, data = data.frame()) { - nu <- get_auxpar(draws$nu, i = i) + nu <- get_dpar(draws$nu, i = i) scale <- ilink(get_eta(draws$mu, i), draws$f$link) / gamma(1 - 1 / nu) args <- list(scale = scale, shape = nu) out <- loglik_censor(dist = "frechet", args = args, @@ -344,7 +344,7 @@ loglik_frechet <- function(i, draws, data = data.frame()) { loglik_gen_extreme_value <- function(i, draws, data = data.frame()) { sigma <- get_sigma(draws$sigma, data = draws$data, i = i) - xi <- get_auxpar(draws$xi, i = i) + xi <- get_dpar(draws$xi, i = i) mu <- ilink(get_eta(draws$mu, i), draws$f$link) args <- nlist(mu, sigma, xi) out <- loglik_censor(dist = "gen_extreme_value", args = args, @@ -367,7 +367,7 @@ loglik_inverse.gaussian <- function(i, draws, data = data.frame()) { loglik_exgaussian <- function(i, draws, data = data.frame()) { args <- list(mu = ilink(get_eta(draws$mu, i), draws$f$link), sigma = get_sigma(draws$sigma, data = draws$data, i = i), - beta = get_auxpar(draws$beta, i = i)) + beta = get_dpar(draws$beta, i = i)) out <- loglik_censor(dist = "exgaussian", args = args, i = i, data = draws$data) out <- loglik_truncate(out, cdf = pexgaussian, args = args, @@ -377,9 +377,9 @@ loglik_exgaussian <- function(i, draws, data = data.frame()) { loglik_wiener <- function(i, draws, data = data.frame()) { args <- list(delta = ilink(get_eta(draws$mu, i), draws$f$link), - alpha = get_auxpar(draws$bs, i = i), - tau = get_auxpar(draws$ndt, i = i), - beta = get_auxpar(draws$bias, i = i), + alpha = get_dpar(draws$bs, i = i), + tau = get_dpar(draws$ndt, i = i), + beta = get_dpar(draws$bias, i = i), resp = draws$data[["dec"]][i]) out <- do.call("dwiener", c(draws$data$Y[i], args, log = TRUE)) loglik_weight(out, i = i, data = draws$data) @@ -387,7 +387,7 @@ loglik_wiener <- function(i, draws, data = data.frame()) { loglik_beta <- function(i, draws, data = data.frame()) { mu <- ilink(get_eta(draws$mu, i), draws$f$link) - phi <- get_auxpar(draws$phi, i) + phi <- get_dpar(draws$phi, i) args <- list(shape1 = mu * phi, shape2 = (1 - mu) * phi) out <- loglik_censor(dist = "beta", args = args, i = i, data = draws$data) out <- loglik_truncate(out, cdf = pbeta, args = args, @@ -397,7 +397,7 @@ loglik_beta <- function(i, draws, data = data.frame()) { loglik_von_mises <- function(i, draws, data = data.frame()) { args <- list(mu = ilink(get_eta(draws$mu, i), draws$f$link), - kappa = get_auxpar(draws$kappa, i = i)) + kappa = get_dpar(draws$kappa, i = i)) out <- loglik_censor(dist = "von_mises", args = args, i = i, data = draws$data) out <- loglik_truncate(out, cdf = pvon_mises, args = args, @@ -408,7 +408,7 @@ loglik_von_mises <- function(i, draws, data = data.frame()) { loglik_asym_laplace <- function(i, draws, ...) { args <- list(mu = ilink(get_eta(draws$mu, i), draws$f$link), sigma = get_sigma(draws$sigma, data = draws$data, i = i), - quantile = get_auxpar(draws$quantile, i = i)) + quantile = get_dpar(draws$quantile, i = i)) out <- loglik_censor(dist = "asym_laplace", args = args, i = i, data = draws$data) out <- loglik_truncate(out, cdf = pvon_mises, args = args, @@ -481,7 +481,7 @@ loglik_zero_inflated_binomial <- function(i, draws, data = data.frame()) { loglik_zero_inflated_beta <- function(i, draws, data = data.frame()) { theta <- get_zi_hu(draws, i, par = "zi") mu <- ilink(get_eta(draws$mu, i), draws$f$link) - phi <- get_auxpar(draws$phi, i) + phi <- get_dpar(draws$phi, i) args <- list(shape1 = mu * phi, shape2 = (1 - mu) * phi) # zi_beta is technically a hurdle model out <- loglik_hurdle_continuous(pdf = dbeta, theta = theta, @@ -490,13 +490,13 @@ loglik_zero_inflated_beta <- function(i, draws, data = data.frame()) { } loglik_zero_one_inflated_beta <- function(i, draws, data = data.frame()) { - zoi <- get_auxpar(draws$zoi, i) - coi <- get_auxpar(draws$coi, i) + zoi <- get_dpar(draws$zoi, i) + coi <- get_dpar(draws$coi, i) if (draws$data$Y[i] %in% c(0, 1)) { out <- dbinom(1, size = 1, prob = zoi, log = TRUE) + dbinom(draws$data$Y[i], size = 1, prob = coi, log = TRUE) } else { - phi <- get_auxpar(draws$phi, i) + phi <- get_dpar(draws$phi, i) mu <- ilink(get_eta(draws$mu, i), draws$f$link) args <- list(shape1 = mu * phi, shape2 = (1 - mu) * phi) out <- dbinom(0, size = 1, prob = zoi, log = TRUE) + @@ -605,13 +605,13 @@ loglik_mixture <- function(i, draws, data = data.frame()) { for (j in seq_along(families)) { loglik_fun <- paste0("loglik_", families[j]) loglik_fun <- get(loglik_fun, asNamespace("brms")) - auxpars <- valid_auxpars(families[j]) + dpars <- valid_dpars(families[j]) tmp_draws <- list( f = draws$f$mix[[j]], nsamples = draws[["nsamples"]], data = draws[["data"]] ) - for (ap in auxpars) { + for (ap in dpars) { tmp_draws[[ap]] <- draws[[paste0(ap, j)]] } out[, j] <- exp(log(theta[, j]) + loglik_fun(i, tmp_draws)) diff --git a/R/make_stancode.R b/R/make_stancode.R index e2e8c8b32..de2e5f2b1 100644 --- a/R/make_stancode.R +++ b/R/make_stancode.R @@ -61,7 +61,7 @@ make_stancode <- function(formula, data, family = gaussian(), bterms, data = data, ranef = ranef, prior = prior, sparse = sparse ) - bterms$auxpars[["mu"]] <- NULL + bterms$dpars[["mu"]] <- NULL } else { text_effects_mv <- list() } @@ -85,8 +85,8 @@ make_stancode <- function(formula, data, family = gaussian(), autocor, bterms = bterms, family = family, prior = prior ) text_mv <- stan_mv(family, response = bterms$response, prior = prior) - disc <- "disc" %in% names(bterms$auxpars) || - isTRUE(bterms$fauxpars$disc != 1) + disc <- "disc" %in% names(bterms$dpars) || + isTRUE(bterms$fdpars$disc$value != 1) text_ordinal <- stan_ordinal( family, prior = prior, cs = has_cs(bterms), disc = disc ) @@ -215,7 +215,7 @@ make_stancode <- function(formula, data, family = gaussian(), ) # generate model block - # list auxpars before pred as part of fixing issue #124 + # list dpars before pred as part of fixing issue #124 text_model_loop <- paste0( text_effects$modelC2, text_autocor$modelC2, diff --git a/R/make_standata.R b/R/make_standata.R index 6c4b958e1..d7c4cc83a 100644 --- a/R/make_standata.R +++ b/R/make_standata.R @@ -101,7 +101,7 @@ make_standata <- function(formula, data, family = gaussian(), if (!factors_allowed && !is.numeric(out$Y)) { stop2("Family '", family4error, "' requires numeric responses.") } - # transform and check response variable for different families + # transform and check response variables for different families regex_pos_int <- "(^|_)(binomial|poisson|negbinomial|geometric)$" if (any(grepl(regex_pos_int, families))) { if (!all(is_wholenumber(out$Y)) || min(out$Y) < 0) { @@ -168,20 +168,21 @@ make_standata <- function(formula, data, family = gaussian(), # data for various kinds of effects only_response <- isTRUE(control$only_response) if (!only_response) { - ranef <- tidy_ranef(bterms, data, ncat = control$ncat, - old_levels = control$old_levels) + ranef <- tidy_ranef( + bterms, data, ncat = control$ncat, + old_levels = control$old_levels + ) args_eff <- nlist(data, ranef, prior, knots, not4stan) resp <- bterms$response if (length(resp) > 1L && !old_mv) { args_eff_spec <- list( - x = bterms$auxpars[["mu"]], smooths = control$smooths[["mu"]], + x = bterms$dpars[["mu"]], smooths = control$smooths[["mu"]], gps = control$gps[["mu"]], Jmo = control$Jmo[["mu"]] ) - bterms$auxpars[["mu"]] <- NULL + bterms$dpars[["mu"]] <- NULL for (r in resp) { - data_eff <- do.call( - data_effects, c(args_eff_spec, args_eff, nlpar = r) - ) + args_eff_spec$x$resp <- r + data_eff <- do.call(data_effects, c(args_eff_spec, args_eff)) out <- c(out, data_eff) } if (is_linear(family)) { @@ -190,18 +191,19 @@ make_standata <- function(formula, data, family = gaussian(), colnames(out$Y) <- resp } } - # data for predictors of auxiliary parameters - for (ap in names(bterms$auxpars)) { + # data for predictors of distributional parameters + for (dp in names(bterms$dpars)) { args_eff_spec <- list( - x = bterms$auxpars[[ap]], nlpar = ap, - smooths = control$smooths[[ap]], - gps = control$gps[[ap]], Jmo = control$Jmo[[ap]] + x = bterms$dpars[[dp]], + smooths = control$smooths[[dp]], + gps = control$gps[[dp]], + Jmo = control$Jmo[[dp]] ) data_aux_eff <- do.call(data_effects, c(args_eff_spec, args_eff)) out <- c(out, data_aux_eff) } - for (ap in names(bterms$fauxpars)) { - out[[ap]] <- bterms$fauxpars[[ap]] + for (dp in names(bterms$fdpars)) { + out[[dp]] <- bterms$fdpars[[dp]]$value } out <- c(out, data_gr(ranef, data, cov_ranef = cov_ranef), @@ -232,8 +234,8 @@ make_standata <- function(formula, data, family = gaussian(), "might be a more efficient choice.") } if (check_response && any(out$Y > out$trials)) { - stop2("Number of trials is smaller than the response ", - "variable would suggest.") + stop2("Number of trials is smaller than ", + "the number of events.") } out$trials <- as.array(out$trials) } diff --git a/R/misc.R b/R/misc.R index b96ff361e..ccccfc023 100644 --- a/R/misc.R +++ b/R/misc.R @@ -6,16 +6,51 @@ p <- function(x, i = NULL, row = TRUE) { # row: indicating if rows or cols should be indexed # only relevant if x has two dimensions if (!length(i)) { - x + out <- x } else if (length(dim(x)) == 2L) { if (row) { - x[i, , drop = FALSE] + out <- x[i, , drop = FALSE] } else { - x[, i, drop = FALSE] + out <- x[, i, drop = FALSE] } } else { - x[i] + out <- x[i] } + out +} + +match_rows <- function(x, y, ...) { + # match rows in x with rows in y + x <- as.data.frame(x) + y <- as.data.frame(y) + x <- do.call("paste", c(x, sep = "\r")) + y <- do.call("paste", c(y, sep = "\r")) + match(x, y, ...) +} + +find_rows <- function(x, ..., ls = list(), fun = '%in%') { + # finding rows matching columns passed via ls and ... + x <- as.data.frame(x) + if (!nrow(x)) { + return(logical(0)) + } + out <- rep(TRUE, nrow(x)) + ls <- c(ls, list(...)) + if (!length(ls)) { + return(out) + } + if (is.null(names(ls))) { + stop("Argument 'ls' must be named.") + } + for (name in names(ls)) { + out <- out & do.call(fun, list(x[[name]], ls[[name]])) + } + out +} + +subset2 <- function(x, ..., ls = list(), fun = '%in%') { + # subset x using arguments passed via ls and ... + x[find_rows(x, ..., ls = ls, fun = fun), ] } select_indices <- function(x, i) { diff --git a/R/predict.R b/R/predict.R index 4e8a4d9b6..4280474ef 100644 --- a/R/predict.R +++ b/R/predict.R @@ -17,7 +17,7 @@ predict_gaussian <- function(i, draws, ...) { } predict_student <- function(i, draws, ...) { - args <- list(df = get_auxpar(draws$nu, i = i), + args <- list(df = get_dpar(draws$nu, i = i), mu = ilink(get_eta(draws$mu, i), draws$f$link), sigma = get_sigma(draws$sigma, data = draws$data, i = i)) rng_continuous(nrng = draws$nsamples, dist = "student_t", args = args, @@ -40,7 +40,7 @@ predict_lognormal <- function(i, draws, ...) { predict_skew_normal <- function(i, draws, ...) { sigma <- get_sigma(draws$sigma, data = draws$data, i = i) - alpha <- get_auxpar(draws$alpha, i = i) + alpha <- get_dpar(draws$alpha, i = i) mu <- ilink(get_eta(draws$mu, i), draws$f$link) args <- nlist(mu, sigma, alpha) rng_continuous(nrng = draws$nsamples, dist = "skew_normal", args = args, @@ -69,7 +69,7 @@ predict_student_mv <- function(i, draws, ...) { } else { mu <- ilink(get_eta(draws$mu, i), draws$f$link)[, 1, ] } - nu <- get_auxpar(draws$nu, i = i) + nu <- get_dpar(draws$nu, i = i) .fun <- function(s) { rmulti_student_t(1, df = nu[s, ], mu = mu[s, ], Sigma = draws$Sigma[s, , ]) @@ -142,7 +142,7 @@ predict_student_cov <- function(i, draws, ...) { # identity matrix Sigma <- do.call(get_cov_matrix_ident, args) } - nu <- get_auxpar(draws$nu, i = i) + nu <- get_dpar(draws$nu, i = i) .fun <- function(s) { rmulti_student_t(1, df = nu[s, ], mu = mu[s, ], Sigma = Sigma[s, , ]) @@ -179,7 +179,7 @@ predict_student_lagsar <- function(i, draws, ...) { N <- draws$data$N mu <- ilink(get_eta(draws$mu), draws$f$link) sigma <- get_sigma(draws$sigma, data = draws$data) - nu <- get_auxpar(draws$nu) + nu <- get_dpar(draws$nu) do.call(rbind, lapply(1:draws$nsamples, .predict_student_lagsar)) } @@ -204,7 +204,7 @@ predict_student_errorsar <- function(i, draws, ...) { } mu <- ilink(get_eta(draws$mu), draws$f$link) sigma <- get_sigma(draws$sigma, data = draws$data) - nu <- get_auxpar(draws$nu) + nu <- get_dpar(draws$nu) do.call(rbind, lapply(1:draws$nsamples, .predict_student_errorsar)) } @@ -220,7 +220,7 @@ predict_gaussian_fixed <- function(i, draws, ...) { predict_student_fixed <- function(i, draws, ...) { stopifnot(i == 1) mu <- ilink(get_eta(draws$mu, 1:nrow(draws$data$V)), draws$f$link) - nu <- get_auxpar(draws$nu, 1:nrow(draws$data$V)) + nu <- get_dpar(draws$nu, 1:nrow(draws$data$V)) .fun <- function(s) { rmulti_student_t(1, df = nu[s, ], mu = mu[s, ], Sigma = draws$data$V) } @@ -292,7 +292,7 @@ predict_weibull <- function(i, draws, ...) { } predict_frechet <- function(i, draws, ...) { - nu <- get_auxpar(draws$nu, i = i) + nu <- get_dpar(draws$nu, i = i) scale <- ilink(get_eta(draws$mu, i), draws$f$link) / gamma(1 - 1 / nu) args <- list(scale = scale, shape = nu) rng_continuous(nrng = draws$nsamples, dist = "frechet", args = args, @@ -301,7 +301,7 @@ predict_frechet <- function(i, draws, ...) { predict_gen_extreme_value <- function(i, draws, ...) { sigma <- get_sigma(draws$sigma, data = draws$data, i = i) - xi <- get_auxpar(draws$xi, i = i) + xi <- get_dpar(draws$xi, i = i) mu <- ilink(get_eta(draws$mu, i), draws$f$link) args <- nlist(mu, sigma, xi) rng_continuous(nrng = draws$nsamples, dist = "gen_extreme_value", @@ -318,7 +318,7 @@ predict_inverse.gaussian <- function(i, draws, ...) { predict_exgaussian <- function(i, draws, ...) { args <- list(mu = ilink(get_eta(draws$mu, i), draws$f$link), sigma = get_sigma(draws$sigma, data = draws$data, i = i), - beta = get_auxpar(draws$beta, i = i)) + beta = get_dpar(draws$beta, i = i)) rng_continuous(nrng = draws$nsamples, dist = "exgaussian", args = args, lb = draws$data$lb[i], ub = draws$data$ub[i]) } @@ -326,9 +326,9 @@ predict_exgaussian <- function(i, draws, ...) { predict_wiener <- function(i, draws, negative_rt = FALSE, ...) { args <- list( delta = ilink(get_eta(draws$mu, i), draws$f$link), - alpha = get_auxpar(draws$bs, i = i), - tau = get_auxpar(draws$ndt, i = i), - beta = get_auxpar(draws$bias, i = i), + alpha = get_dpar(draws$bs, i = i), + tau = get_dpar(draws$ndt, i = i), + beta = get_dpar(draws$bias, i = i), types = if (negative_rt) c("q", "resp") else "q" ) out <- rng_continuous( @@ -344,7 +344,7 @@ predict_wiener <- function(i, draws, negative_rt = FALSE, ...) { predict_beta <- function(i, draws, ...) { mu <- ilink(get_eta(draws$mu, i), draws$f$link) - phi <- get_auxpar(draws$phi, i = i) + phi <- get_dpar(draws$phi, i = i) args <- list(shape1 = mu * phi, shape2 = (1 - mu) * phi) rng_continuous(nrng = draws$nsamples, dist = "beta", args = args, lb = draws$data$lb[i], ub = draws$data$ub[i]) @@ -352,7 +352,7 @@ predict_beta <- function(i, draws, ...) { predict_von_mises <- function(i, draws, ...) { args <- list(mu = ilink(get_eta(draws$mu, i), draws$f$link), - kappa = get_auxpar(draws$kappa, i = i)) + kappa = get_dpar(draws$kappa, i = i)) rng_continuous(nrng = draws$nsamples, dist = "von_mises", args = args, lb = draws$data$lb[i], ub = draws$data$ub[i]) } @@ -360,7 +360,7 @@ predict_von_mises <- function(i, draws, ...) { predict_asym_laplace <- function(i, draws, ...) { args <- list(mu = ilink(get_eta(draws$mu, i), draws$f$link), sigma = get_sigma(draws$sigma, data = draws$data, i = i), - quantile = get_auxpar(draws$quantile, i = i)) + quantile = get_dpar(draws$quantile, i = i)) rng_continuous(nrng = draws$nsamples, dist = "asym_laplace", args = args, lb = draws$data$lb[i], ub = draws$data$ub[i]) } @@ -418,7 +418,7 @@ predict_zero_inflated_beta <- function(i, draws, ...) { # theta is the bernoulli hurdle parameter theta <- get_zi_hu(draws, i, par = "zi") mu <- ilink(get_eta(draws$mu, i), draws$f$link) - phi <- get_auxpar(draws$phi, i = i) + phi <- get_dpar(draws$phi, i = i) # compare with theta to incorporate the hurdle process hu <- runif(draws$nsamples, 0, 1) ifelse(hu < theta, 0, rbeta(draws$nsamples, shape1 = mu * phi, @@ -426,10 +426,10 @@ predict_zero_inflated_beta <- function(i, draws, ...) { } predict_zero_one_inflated_beta <- function(i, draws, ...) { - zoi <- get_auxpar(draws$zoi, i) - coi <- get_auxpar(draws$coi, i) + zoi <- get_dpar(draws$zoi, i) + coi <- get_dpar(draws$coi, i) mu <- ilink(get_eta(draws$mu, i), draws$f$link) - phi <- get_auxpar(draws$phi, i = i) + phi <- get_dpar(draws$phi, i = i) hu <- runif(draws$nsamples, 0, 1) one_or_zero <- runif(draws$nsamples, 0, 1) ifelse(hu < zoi, @@ -512,13 +512,13 @@ predict_mixture <- function(i, draws, ...) { if (length(sample_ids)) { predict_fun <- paste0("predict_", families[j]) predict_fun <- get(predict_fun, asNamespace("brms")) - auxpars <- valid_auxpars(families[j]) + dpars <- valid_dpars(families[j]) tmp_draws <- list( f = draws$f$mix[[j]], nsamples = length(sample_ids), data = draws[["data"]] ) - for (ap in auxpars) { + for (ap in dpars) { tmp_draws[[ap]] <- p(draws[[paste0(ap, j)]], sample_ids, row = TRUE) } out[sample_ids] <- predict_fun(i, tmp_draws, ...) diff --git a/R/priors.R b/R/priors.R index 4dc8ca5c4..25e06c098 100644 --- a/R/priors.R +++ b/R/priors.R @@ -10,11 +10,12 @@ #' See 'Details' for other valid parameter classes. #' @param coef Name of the (population- or group-level) parameter. #' @param group Grouping factor of group-level parameters. -#' @param nlpar Name of a non-linear / auxiliary parameter. -#' Only used in non-linear / distributional models. #' @param resp Name of the response variable / category. #' Only used in multivariate and categorical models. -#' Is internally handled as an alias of \code{nlpar}. +#' @param dpar Name of a distributional parameter. +#' Only used in distributional models. +#' @param nlpar Name of a non-linear parameter. +#' Only used in non-linear models. #' @param lb Lower bound for parameter restriction. Currently only allowed #' for classes \code{"b"}, \code{"ar"}, \code{"ma"}, and \code{"arr"}. #' Defaults to \code{NULL}, that is no restriction. @@ -312,57 +313,43 @@ #' #' @export set_prior <- function(prior, class = "b", coef = "", group = "", - nlpar = "", resp = NULL, lb = NULL, ub = NULL, - check = TRUE) { + resp = "", dpar = "", nlpar = "", + lb = NULL, ub = NULL, check = TRUE) { prior <- as.character(prior) class <- as.character(class) group <- as.character(group) coef <- as.character(coef) - nlpar <- as.character(use_alias(nlpar, resp, warn = FALSE)) + resp <- as.character(resp) + dpar <- as.character(dpar) + nlpar <- as.character(nlpar) lb <- as.numeric(lb) ub <- as.numeric(ub) - check <- as.logical(check)[1] + check <- as_one_logical(check) if (length(prior) != 1 || length(class) != 1 || length(coef) != 1 || - length(group) != 1 || length(nlpar) != 1 || length(lb) > 1 || - length(ub) > 1) { + length(group) != 1 || length(resp) != 1 || length(dpar) != 1 || + length(nlpar) != 1 || length(lb) > 1 || length(ub) > 1) { stop2("All arguments of set_prior must be of length 1.") } - - valid_classes <- c( - "Intercept", "b", "sd", "sds", "sdgp", "lscale", "simplex", "cor", - "L", "ar", "ma", "arr", "lagsar", "errorsar", "car", "sdcar", - "sigmaLL", "rescor", "Lrescor", "delta", "theta", if (!check) "" - ) - if (!(class %in% valid_classes || auxpar_class(class) %in% auxpars())) { - stop2("'", class, "' is not a valid parameter class.") - } - if (nchar(group) && !class %in% c("sd", "cor", "L")) { - stop2("Argument 'group' is not meaningful for class '", class, "'.") - } - coef_classes <- c( - "Intercept", "b", "sd", "sds", "sigma", - "simplex", "sdgp", "lscale" - ) - if (nchar(coef) && !class %in% coef_classes) { - stop2("Argument 'coef' ia not meaningful for class '", class, "'.") - } - if (nchar(nlpar) && !class %in% valid_classes[1:7]) { - stop2("Argument 'nlpar' is not meaningful for class '", class, "'.") - } + # validate boundaries is_arma <- class %in% c("ar", "ma") if (length(lb) || length(ub) || is_arma) { - if (!(class %in% c("b", "arr") || is_arma)) - stop2("Currently boundaries are only allowed for ", - "population-level and autocorrelation parameters.") - if (nchar(coef)) { + if (!(class %in% c("b", "ar", "ma", "arr"))) { + stop2( + "Currently boundaries are only allowed for ", + "population-level and autocorrelation parameters." + ) + } + if (nzchar(coef)) { stop2("Argument 'coef' may not be specified when using boundaries.") } if (is_arma) { lb <- ifelse(length(lb), lb, -1) ub <- ifelse(length(ub), ub, 1) if (is.na(lb) || is.na(ub) || abs(lb) > 1 || abs(ub) > 1) { - warning2("Setting boundaries of autocorrelation parameters ", - "outside of [-1,1] may not be appropriate.") + warning2( + "Setting boundaries of autocorrelation parameters ", + "outside of [-1,1] may not be appropriate." + ) } } # don't put spaces in boundary declarations @@ -378,9 +365,11 @@ set_prior <- function(prior, class = "b", coef = "", group = "", } if (!check) { # prior will be added to the log-posterior as is - class <- coef <- group <- nlpar <- bound <- "" + class <- coef <- group <- resp <- dpar <- nlpar <- bound <- "" } - do.call(brmsprior, nlist(prior, class, coef, group, nlpar, bound)) + do.call(brmsprior, + nlist(prior, class, coef, group, resp, dpar, nlpar, bound) + ) } #' @describeIn set_prior Alias of \code{set_prior} allowing to @@ -401,7 +390,6 @@ prior_ <- function(prior, ...) { call <- nlist(prior, ...) seval <- rmNULL(call[prior_seval_args()]) call[prior_seval_args()] <- NULL - as_string <- function(x) { if (is.formula(x) && length(x) == 2) { deparse_no_string(x[[2]]) @@ -476,7 +464,8 @@ get_prior <- function(formula, data, family = gaussian(), data <- update_data(data, bterms = bterms) ranef <- tidy_ranef(bterms, data) - # ensure that RE and residual SDs only have a weakly informative prior by default + # ensure that RE and residual SDs only have + # a weakly informative prior by default Y <- unname(model.response(data)) prior_scale <- 10 if (is_lognormal(family)) { @@ -498,23 +487,24 @@ get_prior <- function(formula, data, family = gaussian(), if (length(bterms$response) > 1L) { # priors for effects in multivariate models for (r in bterms$response) { + bterms$dpars[["mu"]]$resp <- r prior_eff <- prior_effects( - bterms$auxpars[["mu"]], data = data, + bterms$dpars[["mu"]], data = data, def_scale_prior = def_scale_prior, - nlpar = r, internal = internal + internal = internal ) - prior <- rbind(prior, prior_eff) + prior <- prior + prior_eff } - bterms$auxpars[["mu"]] <- NULL + bterms$dpars[["mu"]] <- NULL # add "global" priors for population-level effects # in 1.8.0 as users keep asking about this for (cl in c("b", "Intercept")) { if (any(with(prior, class == cl & coef == ""))) { - prior <- rbind(prior, brmsprior(class = cl)) + prior <- prior + brmsprior(class = cl) } } } - # priors for auxiliary parameters + # priors for distributional parameters def_auxprior <- c( sigma = def_scale_prior, shape = "gamma(0.01, 0.01)", @@ -535,20 +525,20 @@ get_prior <- function(formula, data, family = gaussian(), disc = NA, mu = NA ) - valid_auxpars <- valid_auxpars(family, bterms = bterms) - for (ap in valid_auxpars) { - ap_class <- auxpar_class(ap) - if (!is.null(bterms$auxpars[[ap]])) { - auxprior <- prior_effects( - bterms$auxpars[[ap]], data = data, nlpar = ap, + valid_dpars <- valid_dpars(family, bterms = bterms) + for (dp in valid_dpars) { + dp_class <- dpar_class(dp) + if (!is.null(bterms$dpars[[dp]])) { + dp_prior <- prior_effects( + bterms$dpars[[dp]], data = data, def_scale_prior = def_scale_prior ) - } else if (!is.na(def_auxprior[ap_class])) { - auxprior <- brmsprior(class = ap, prior = def_auxprior[ap_class]) + } else if (!is.na(def_auxprior[dp_class])) { + dp_prior <- brmsprior(class = dp, prior = def_auxprior[dp_class]) } else { - auxprior <- empty_brmsprior() + dp_prior <- empty_brmsprior() } - prior <- rbind(prior, auxprior) + prior <- prior + dp_prior } # priors of group-level parameters prior_re <- prior_re( @@ -556,73 +546,72 @@ get_prior <- function(formula, data, family = gaussian(), global_sd = length(bterms$response) > 1L, internal = internal ) - prior <- rbind(prior, prior_re) + prior <- prior + prior_re # prior for the delta parameter for equidistant thresholds if (is_ordinal(family) && is_equal(family$threshold, "equidistant")) { bound <- ifelse(family$family == "cumulative", "", "") - prior <- rbind(prior, brmsprior(class = "delta", bound = bound)) + prior <- prior + brmsprior(class = "delta", bound = bound) } # priors for mixture models - ap_classes <- auxpar_class(names(c(bterms$auxpars, bterms$fauxpars))) + ap_classes <- dpar_class(names(c(bterms$dpars, bterms$fdpars))) if (is.mixfamily(family) && !any(ap_classes == "theta")) { - prior <- rbind(prior, brmsprior(class = "theta")) + prior <- prior + brmsprior(class = "theta") } - # priors for auxiliary parameters of multivariate models + # priors for distributional parameters of multivariate models if (is_linear(family) && length(bterms$response) > 1L) { sigma_coef <- c("", bterms$response) sigma_prior <- c(def_scale_prior, rep("", length(bterms$response))) sigma_prior <- brmsprior( class = "sigma", coef = sigma_coef, prior = sigma_prior ) - prior <- rbind(prior, sigma_prior) + prior <- prior + sigma_prior if (internal) { - prior <- rbind(prior, + prior <- prior + brmsprior(class = "Lrescor", prior = "lkj_corr_cholesky(1)") - ) } else { - prior <- rbind(prior, brmsprior(class = "rescor", prior = "lkj(1)")) + prior <- prior + brmsprior(class = "rescor", prior = "lkj(1)") } } # priors for autocor parameters cbound <- "" if (get_ar(autocor)) { - prior <- rbind(prior, brmsprior(class = "ar", bound = cbound)) + prior <- prior + brmsprior(class = "ar", bound = cbound) } if (get_ma(autocor)) { - prior <- rbind(prior, brmsprior(class = "ma", bound = cbound)) + prior <- prior + brmsprior(class = "ma", bound = cbound) } if (get_arr(autocor)) { - prior <- rbind(prior, brmsprior(class = "arr")) + prior <- prior + brmsprior(class = "arr") } if (is.cor_sar(autocor)) { if (identical(autocor$type, "lag")) { - prior <- rbind(prior, brmsprior(class = "lagsar")) + prior <- prior + brmsprior(class = "lagsar") } if (identical(autocor$type, "error")) { - prior <- rbind(prior, brmsprior(class = "errorsar")) + prior <- prior + brmsprior(class = "errorsar") } } if (is.cor_car(autocor)) { - prior <- rbind(prior, + prior <- prior + brmsprior(def_scale_prior, class = "sdcar") - ) if (identical(autocor$type, "escar")) { - prior <- rbind(prior, brmsprior(class = "car")) + prior <- prior + brmsprior(class = "car") } } if (is(autocor, "cor_bsts")) { - prior <- rbind(prior, + prior <- prior + brmsprior(class = "sigmaLL", prior = def_scale_prior) - ) } # do not remove unique(.) - prior <- unique(prior[with(prior, order(nlpar, class, group, coef)), ]) + prior <- unique(prior[with(prior, + order(resp, dpar, nlpar, class, group, coef) + ), ]) rownames(prior) <- NULL structure(prior, class = c("brmsprior", "data.frame")) } #' @export -prior_effects.btl <- function(x, data, nlpar = "", spec_intercept = TRUE, +prior_effects.btl <- function(x, data, spec_intercept = TRUE, def_scale_prior = "", ...) { # collect default priors for various kinds of effects # Args: @@ -630,21 +619,12 @@ prior_effects.btl <- function(x, data, nlpar = "", spec_intercept = TRUE, # def_scale_prior: default prior for SD parameters # Return: # An object of class brmsprior - nlpar <- check_nlpar(nlpar) - fixef <- colnames(data_fe(x, data)$X) - spec_intercept <- has_intercept(x$fe) && spec_intercept - prior_fe <- prior_fe( - fixef, spec_intercept = spec_intercept, nlpar = nlpar - ) - monef <- all_terms(x$mo) - prior_mo <- prior_mo(monef, fixef = fixef, nlpar = nlpar) - smooths <- get_sm_labels(x) - prior_sm <- prior_sm(smooths, def_scale_prior, nlpar = nlpar) - csef <- colnames(get_model_matrix(x$cs, data = data)) - prior_cs <- prior_cs(csef, fixef = fixef, nlpar = nlpar) - prior_me <- prior_me(get_me_labels(x, data), nlpar = nlpar) - prior_gp <- prior_gp(get_gp_labels(x), def_scale_prior, nlpar = nlpar) - rbind(prior_fe, prior_mo, prior_sm, prior_cs, prior_me, prior_gp) + prior_fe(x, data, spec_intercept = spec_intercept) + + prior_cs(x, data) + + prior_mo(x, data) + + prior_sm(x, data, def_scale_prior = def_scale_prior) + + prior_me(x, data) + + prior_gp(x, data, def_scale_prior = def_scale_prior) } #' @export @@ -656,110 +636,95 @@ prior_effects.btnl <- function(x, data, def_scale_prior = "", ...) { prior <- empty_brmsprior() for (i in seq_along(nlpars)) { prior_eff <- prior_effects( - x$nlpars[[i]], data = data, nlpar = nlpars[i], + x$nlpars[[i]], data = data, def_scale_prior = def_scale_prior, spec_intercept = FALSE ) - prior <- rbind(prior, prior_eff) + prior <- prior + prior_eff } prior } -prior_fe <- function(fixef, spec_intercept = TRUE, nlpar = "") { +prior_fe <- function(bterms, data, spec_intercept = TRUE) { # priors for population-level parameters # Args: - # fixef: names of the population-level effects # spec_intercept: special parameter class for the Intercept? # Returns: # an object of class brmsprior prior <- empty_brmsprior() - if (spec_intercept) { - prior <- rbind(prior, - brmsprior(class = "Intercept", coef = "", nlpar = nlpar), - brmsprior(class = "b", coef = "Intercept", nlpar = nlpar) - ) + fixef <- colnames(data_fe(bterms, data)$X) + px <- check_prefix(bterms) + if (has_intercept(bterms$fe) && spec_intercept) { + prior <- prior + + brmsprior(class = "Intercept", coef = "", ls = px) + + brmsprior(class = "b", coef = "Intercept", ls = px) fixef <- setdiff(fixef, "Intercept") } if (length(fixef)) { - prior <- rbind(prior, - brmsprior(class = "b", coef = c("", fixef), nlpar = nlpar) - ) + prior <- prior + + brmsprior(class = "b", coef = c("", fixef), ls = px) } prior } -prior_mo <- function(monef, fixef = NULL, nlpar = "") { +prior_mo <- function(bterms, data) { # priors for monotonic effects parameters - # Args: - # monef: names of the monotonic effects - # fixef: names of the population-level effects - # nlpar: optional name of a non-linear parameter # Returns: # an object of class brmsprior prior <- empty_brmsprior() + monef <- all_terms(bterms$mo) if (length(monef)) { - invalid <- intersect(fixef, monef) - if (length(invalid)) { - stop2("Variables cannot be modeled as fixed and ", - "monotonic effects at the same time. ", - "\nError occured for variables: ", - collapse_comma(invalid)) - } - prior <- rbind( - brmsprior(class = "b", coef = c("", monef), nlpar = nlpar), - brmsprior(class = "simplex", coef = monef, nlpar = nlpar) - ) + px <- check_prefix(bterms) + prior <- prior + + brmsprior(class = "b", coef = c("", monef), ls = px) + + brmsprior(class = "simplex", coef = monef, ls = px) } prior } -prior_cs <- function(csef, fixef = NULL, nlpar = "") { +prior_cs <- function(bterms, data) { # priors for category spcific effects parameters - # Args: - # csef: names of the category specific effects - # fixef: names of the population-level effects # Returns: # an object of class brmsprior prior <- empty_brmsprior() + csef <- colnames(get_model_matrix(bterms$cs, data = data)) if (length(csef)) { - stopifnot(!nzchar(nlpar)) - invalid <- intersect(fixef, csef) - if (length(invalid)) { - stop2("Variables cannot be modeled as fixed and ", - "category specific effects at the same time. ", - "\nError occured for variables: ", - collapse_comma(invalid)) - } - prior <- brmsprior(class = "b", coef = c("", csef)) + px <- check_prefix(bterms) + prior <- prior + + brmsprior(class = "b", coef = c("", csef), ls = px) } prior } -prior_me <- function(meef, nlpar = "") { +prior_me <- function(bterms, data) { # default priors of coefficients of noisy terms - # Args: - # meef: terms containing noisy variables + # Returns: + # an object of class brmsprior prior <- empty_brmsprior() + meef <- get_me_labels(bterms, data) if (length(meef)) { - prior <- brmsprior( - class = "b", coef = c("", rename(meef)), nlpar = nlpar - ) + px <- check_prefix(bterms) + prior <- prior + + brmsprior(class = "b", coef = c("", rename(meef)), ls = px) } prior } -prior_gp <- function(gpef, def_scale_prior, nlpar = "") { - # default priors of coefficients of noisy terms - # Args: - # meef: terms containing noisy variables +prior_gp <- function(bterms, data, def_scale_prior) { + # default priors of gaussian processes + # Returns: + # an object of class brmsprior + # def_scale_prior: a character string defining + # the default prior for random effects SDs prior <- empty_brmsprior() + gpef <- get_gp_labels(bterms) if (length(gpef)) { - prior <- rbind( - brmsprior(class = "sdgp", prior = def_scale_prior, nlpar = nlpar), - brmsprior(class = "sdgp", coef = gpef, nlpar = nlpar), - brmsprior(class = "lscale", prior = "normal(0, 0.5)", nlpar = nlpar), - brmsprior(class = "lscale", coef = gpef, nlpar = nlpar) - ) + px <- check_prefix(bterms) + prior <- prior + + brmsprior(class = "sdgp", prior = def_scale_prior, ls = px) + + brmsprior(class = "sdgp", coef = gpef, ls = px) + + brmsprior(class = "lscale", prior = "normal(0, 0.5)", ls = px) + + brmsprior(class = "lscale", coef = gpef, ls = px) } prior } @@ -777,68 +742,74 @@ prior_re <- function(ranef, def_scale_prior, global_sd = FALSE, # Returns: # an object of class brmsprior prior <- empty_brmsprior() - if (nrow(ranef)) { - # global sd class - nlpars <- unique(ranef$nlpar) - if (global_sd) { - global_sd_prior <- rep("", length(setdiff(nlpars, ""))) - global_sd_prior <- c(def_scale_prior, global_sd_prior) - global_sd_prior <- brmsprior( - class = "sd", prior = global_sd_prior, nlpar = union("", nlpars) - ) - } else { - global_sd_prior <- - brmsprior(class = "sd", prior = def_scale_prior, nlpar = nlpars) + if (!nrow(ranef)) { + return(prior) + } + # global sd class + px <- check_prefix(ranef) + upx <- unique(px) + if (global_sd) { + global_sd_prior <- rep("", nrow(upx)) + global_sd_prior <- c(def_scale_prior, global_sd_prior) + upx <- lapply(upx, function(x) union("", x)) + global_sd_prior <- brmsprior( + class = "sd", prior = global_sd_prior, ls = upx + ) + } else { + global_sd_prior <- brmsprior( + class = "sd", prior = def_scale_prior, ls = px + ) + } + prior <- prior + global_sd_prior + for (id in unique(ranef$id)) { + r <- subset2(ranef, id = id) + group <- r$group[1] + rpx <- check_prefix(r) + urpx <- unique(rpx) + # include group-level standard deviations + prior <- prior + + brmsprior(class = "sd", group = group, ls = urpx) + + brmsprior(class = "sd", coef = r$coef, group = group, ls = rpx) + # detect duplicated group-level effects + J <- with(prior, class == "sd" & nzchar(coef)) + dupli <- duplicated(prior[J, ]) + if (any(dupli)) { + stop2("Duplicated group-level effects detected for group ", group) } - prior <- rbind(prior, global_sd_prior) - for (id in unique(ranef$id)) { - r <- ranef[ranef$id == id, ] - group <- r$group[1] - # include group-level standard deviations - prior <- rbind(prior, - brmsprior(class = "sd", group = group, - nlpar = unique(r$nlpar)), - brmsprior(class = "sd", coef = r$coef, - group = group, nlpar = r$nlpar) - ) - # detect duplicated group-level effects - J <- with(prior, class == "sd" & nchar(coef)) - dupli <- duplicated(prior[J, ]) - if (any(dupli)) { - stop2("Duplicated group-level effects detected for group ", group) - } - # include correlation parameters - if (isTRUE(r$cor[1]) && nrow(r) > 1L) { - if (internal) { - prior <- rbind(prior, - brmsprior(class = "L", group = c("", group), - prior = c("lkj_corr_cholesky(1)", "")) + # include correlation parameters + if (isTRUE(r$cor[1]) && nrow(r) > 1L) { + if (internal) { + prior <- prior + + brmsprior( + class = "L", group = c("", group), + prior = c("lkj_corr_cholesky(1)", "") ) - } else { - prior <- rbind(prior, - brmsprior(class = "cor", group = c("", group), - prior = c("lkj(1)", "")) + } else { + prior <- prior + + brmsprior( + class = "cor", group = c("", group), + prior = c("lkj(1)", "") ) - } } } - } + } prior } -prior_sm <- function(smooths, def_scale_prior, nlpar = "") { +prior_sm <- function(bterms, data, def_scale_prior) { # priors for smooth terms # Args: - # smooths: names of the smooth terms - # def_scale_prior: a character string defining the default - # prior for smooth SDs - # nlpar: optional name of a non-linear parameter + # def_scale_prior: a character string defining + # the default prior for smooth SDs + prior <- empty_brmsprior() + smooths <- get_sm_labels(bterms) if (length(smooths)) { + px <- check_prefix(bterms) prior_strings <- c(def_scale_prior, rep("", length(smooths))) - prior <- brmsprior(class = "sds", coef = c("", smooths), - prior = prior_strings, nlpar = nlpar) - } else { - prior <- empty_brmsprior() + prior <- prior + brmsprior( + class = "sds", coef = c("", smooths), + prior = prior_strings, ls = px + ) } prior } @@ -870,7 +841,7 @@ check_prior <- function(prior, formula, data = NULL, family = NULL, if (is.null(prior)) { prior <- all_priors } - # temporarily exclude priors priors that should not be checked + # temporarily exclude priors that should not be checked no_checks <- !nzchar(prior$class) prior_no_checks <- prior[no_checks, ] prior <- prior[!no_checks, ] @@ -879,13 +850,34 @@ check_prior <- function(prior, formula, data = NULL, family = NULL, prior$class, symbols = c("^cor$", "^rescor$"), subs = c("L", "Lrescor"), fixed = FALSE ) - duplicated_input <- duplicated(prior[, 2:5]) + duplicated_input <- duplicated(prior[, 2:7]) if (any(duplicated_input)) { stop2("Duplicated prior specifications are not allowed.") } + valid_dpars <- valid_dpars(family, bterms) + nlpars_in_dpars <- prior$nlpar %in% valid_dpars + if (any(nlpars_in_dpars)) { + warning2( + "Specifying priors of distributional parameters via ", + "'nlpar' is deprecated. Please use 'dpar' instead." + ) + prior$dpar[nlpars_in_dpars] <- prior$nlpar[nlpars_in_dpars] + prior$nlpar[nlpars_in_dpars] <- "" + } + if (length(bterms$response) > 1L) { + nlpars_in_resp <- prior$nlpar %in% bterms$response + if (any(nlpars_in_resp)) { + warning2( + "Specifying priors in multivariate models via ", + "'nlpar' is deprecated. Please use 'resp' instead." + ) + prior$resp[nlpars_in_resp] <- prior$nlpar[nlpars_in_resp] + prior$nlpar[nlpars_in_resp] <- "" + } + } # check if parameters in prior are valid if (nrow(prior)) { - valid <- which(duplicated(rbind(all_priors[, 2:5], prior[, 2:5]))) + valid <- which(duplicated(rbind(all_priors[, 2:7], prior[, 2:7]))) invalid <- which(!1:nrow(prior) %in% (valid - nrow(all_priors))) if (length(invalid)) { msg_priors <- .print_prior(prior[invalid, ]) @@ -899,27 +891,32 @@ check_prior <- function(prior, formula, data = NULL, family = NULL, prior$prior <- sub("^(lkj|lkj_corr)\\(", "lkj_corr_cholesky(", prior$prior) check_prior_content(prior, family = family, warn = warn) # check if priors for non-linear parameters are defined - nlpars <- names(bterms$auxpars$mu$nlpars) - for (nlp in nlpars) { - nlp_prior <- prior$prior[with(prior, nlpar == nlp & class == "b")] - if (!any(nzchar(nlp_prior))) { - stop2("Priors on population-level effects are required in ", - "non-linear models,\nbut none were found for parameter ", - "'", nlp, "'. \nSee help(set_prior) for more details.") + for (dp in names(bterms$dpars)) { + nlpars <- names(bterms$dpars[[dp]]$nlpars) + dp <- ifelse(dp == "mu", "", dp) + for (nlp in nlpars) { + nlp_prior <- subset2(prior, dpar = dp, nlpar = nlp, class = "b") + if (!any(nzchar(nlp_prior$prior))) { + stop2( + "Priors on population-level effects are required in ", + "non-linear models,but none were found for parameter ", + "'", nlp, "'. See help(set_prior) for more details." + ) + } } } # prepare special priors for use in Stan prior <- check_prior_special(bterms, prior) # merge user-specified priors with default priors - prior <- rbind(prior, all_priors) - prior <- prior[!duplicated(prior[, 2:5]), ] + prior <- prior + all_priors + prior <- prior[!duplicated(prior[, 2:7]), ] rows2remove <- NULL # copy over the global population-level prior in MV models if (length(bterms$response) > 1L) { for (cl in c("b", "Intercept")) { - g_index <- with(prior, nlpar == "" & class == cl & coef == "") - for (resp in bterms$response) { - r_index <- with(prior, nlpar == resp & class == cl & coef == "") + g_index <- find_rows(prior, class = cl, coef = "", resp = "") + for (r in bterms$response) { + r_index <- find_rows(prior, class = cl, coef = "", resp = r) if (isTRUE(!nzchar(prior$prior[r_index]))) { prior$prior[r_index] <- prior$prior[g_index] } @@ -982,8 +979,8 @@ check_prior <- function(prior, formula, data = NULL, family = NULL, if (length(rows2remove)) { prior <- prior[-rows2remove, ] } - prior <- prior[with(prior, order(nlpar, class, group, coef)), ] - prior <- rbind(prior, prior_no_checks) + prior <- prior[with(prior, order(resp, dpar, nlpar, class, group, coef)), ] + prior <- prior + prior_no_checks rownames(prior) <- NULL attr(prior, "sample_prior") <- sample_prior attr(prior, "checked") <- TRUE @@ -1049,9 +1046,11 @@ check_prior_content <- function(prior, family = gaussian(), warn = TRUE) { } } else if (prior$class[i] %in% cor_pars) { if (nchar(prior$prior[i]) && !grepl("^lkj", prior$prior[i])) { - stop2("Currently 'lkj' is the only valid prior ", - "for group-level correlations. See help(set_prior) ", - "for more details.") + stop2( + "Currently 'lkj' is the only valid prior ", + "for group-level correlations. See help(set_prior) ", + "for more details." + ) } } else if (prior$class[i] %in% autocor_pars) { if (prior$bound[i] != "") { @@ -1059,29 +1058,37 @@ check_prior_content <- function(prior, family = gaussian(), warn = TRUE) { } } else if (prior$class[i] %in% c("simplex", "theta")) { if (nchar(prior$prior[i]) && !grepl("^dirichlet\\(", prior$prior[i])) { - stop2("Currently 'dirichlet' is the only valid prior ", - "for simplex parameters. See help(set_prior) ", - "for more details.") + stop2( + "Currently 'dirichlet' is the only valid prior ", + "for simplex parameters. See help(set_prior) ", + "for more details." + ) } } } if (nchar(lb_warning) && warn) { - warning2("It appears as if you have specified a lower bounded ", - "prior on a parameter that has no natural lower bound.", - "\nIf this is really what you want, please specify ", - "argument 'lb' of 'set_prior' appropriately.", - "\nWarning occurred for prior \n", lb_warning) + warning2( + "It appears as if you have specified a lower bounded ", + "prior on a parameter that has no natural lower bound.", + "\nIf this is really what you want, please specify ", + "argument 'lb' of 'set_prior' appropriately.", + "\nWarning occurred for prior \n", lb_warning + ) } if (nchar(ub_warning) && warn) { - warning2("It appears as if you have specified an upper bounded ", - "prior on a parameter that has no natural upper bound.", - "\nIf this is really what you want, please specify ", - "argument 'ub' of 'set_prior' appropriately.", - "\nWarning occurred for prior \n", ub_warning) + warning2( + "It appears as if you have specified an upper bounded ", + "prior on a parameter that has no natural upper bound.", + "\nIf this is really what you want, please specify ", + "argument 'ub' of 'set_prior' appropriately.", + "\nWarning occurred for prior \n", ub_warning + ) } if (autocor_warning && warn) { - warning2("Changing the boundaries of autocorrelation ", - "parameters is not recommended.") + warning2( + "Changing the boundaries of autocorrelation ", + "parameters is not recommended." + ) } } invisible(TRUE) @@ -1095,41 +1102,36 @@ check_prior_special.brmsterms <- function(x, prior = NULL, ...) { if (is.null(prior)) { prior <- empty_brmsprior() } - simple_sigma <- has_sigma(x$family, x) && is.null(x$auxpars$sigma) - for (ap in names(x$auxpars)) { - allow_autoscale <- simple_sigma && identical(ap, "mu") + simple_sigma <- has_sigma(x$family, x) && is.null(x$dpars$sigma) + for (dp in names(x$dpars)) { + allow_autoscale <- simple_sigma && identical(dp, "mu") prior <- check_prior_special( - x$auxpars[[ap]], prior, nlpar = ap, - allow_autoscale = allow_autoscale, ... + x$dpars[[dp]], prior, allow_autoscale = allow_autoscale, ... ) } prior } #' @export -check_prior_special.btnl <- function(x, prior, nlpar = "", ...) { +check_prior_special.btnl <- function(x, prior, ...) { stopifnot(is.brmsprior(prior)) for (nlp in names(x$nlpars)) { - prior <- check_prior_special(x$nlpars[[nlp]], prior, nlpar = nlp, ...) + prior <- check_prior_special(x$nlpars[[nlp]], prior, ...) } prior } #' @export -check_prior_special.btl <- function(x, prior, nlpar = "", - allow_autoscale = TRUE, ...) { +check_prior_special.btl <- function(x, prior, allow_autoscale = TRUE, ...) { # prepare special priors such as horseshoe or lasso # Args: # prior: an object of class brmsprior # allow_autoscale: allow autoscaling using sigma? # Returns: # a possibly amended brmsprior object with additional attributes - nlpar_original <- nlpar - nlpar <- check_nlpar(nlpar) + px <- check_prefix(x) prior_special <- list() - b_index <- which( - prior$class == "b" & !nchar(prior$coef) & prior$nlpar == nlpar - ) + b_index <- which(find_rows(prior, class = "b", coef = "", ls = px)) stopifnot(length(b_index) <= 1L) if (length(b_index)) { b_prior <- prior$prior[b_index] @@ -1144,13 +1146,15 @@ check_prior_special.btl <- function(x, prior, nlpar = "", "in models with special population-level effects.") } b_coef_indices <- which( - prior$class == "b" & nchar(prior$coef) & - prior$coef != "Intercept" & prior$nlpar == nlpar + find_rows(prior, class = "b", ls = px) & + !find_rows(prior, coef = c("", "Intercept")) ) if (any(nchar(prior$prior[b_coef_indices]))) { - stop2("Defining priors for single population-level parameters", - "is not allowed when using horseshoe or lasso priors", - "(except for the Intercept).") + stop2( + "Defining priors for single population-level parameters", + "is not allowed when using horseshoe or lasso priors", + "(except for the Intercept)." + ) } if (grepl("^horseshoe\\(", b_prior)) { hs <- eval2(b_prior) @@ -1169,9 +1173,9 @@ check_prior_special.btl <- function(x, prior, nlpar = "", # the parameterization via double_exponential appears to be more # efficient than an indirect parameterization via normal and # exponential distributions; tested on 2017-06-09 - usc_nlpar <- usc(nlpar) + p <- usc(combine_prefix(px)) lasso_scale <- paste0( - "lasso_scale", usc_nlpar, " * lasso_inv_lambda", usc_nlpar + "lasso_scale", p, " * lasso_inv_lambda", p ) lasso_prior <- paste0( "double_exponential(0, ", lasso_scale, ")" @@ -1183,7 +1187,8 @@ check_prior_special.btl <- function(x, prior, nlpar = "", } } } - attributes(prior)$special[[nlpar_original]] <- prior_special + prefix <- combine_prefix(px, keep_mu = TRUE) + attributes(prior)$special[[prefix]] <- prior_special prior } @@ -1198,7 +1203,7 @@ check_sample_prior <- function(sample_prior) { } get_bound <- function(prior, class = "b", coef = "", - group = "", nlpar = "") { + group = "", px = list()) { # extract the boundaries of a parameter described by class etc. # Args: # prior: object of class brmsprior @@ -1206,30 +1211,48 @@ get_bound <- function(prior, class = "b", coef = "", stopifnot(length(class) == 1L) if (!length(coef)) coef <- "" if (!length(group)) group <- "" - if (!length(nlpar)) nlpar <- "" - take <- prior$class == class & prior$coef == coef & - prior$group == group & prior$nlpar == nlpar - if (sum(take) > 1L) { + bound <- subset2(prior, ls = c(nlist(class, coef, group), px))$bound + if (length(bound) > 1L) { stop("extracted more than one boundary at once") } - prior$bound[take] + bound } brmsprior <- function(prior = "", class = "", coef = "", group = "", - nlpar = "", bound = "") { - # helper function to create data.frames containing prior information - out <- data.frame(prior = prior, class = class, coef = coef, - group = group, nlpar = nlpar, bound = bound, - stringsAsFactors = FALSE) + resp = "", dpar = "", nlpar = "", bound = "", + ls = list()) { + # helper function to create data.frames containing prior information + if (length(ls)) { + if (is.null(names(ls))) { + stop("Argument 'ls' must be named.") + } + names <- c( + "prior", "class", "coef", "group", + "resp", "dpar", "nlpar", "bound" + ) + if (!all(names(ls) %in% names)) { + stop("Names of 'ls' must some of ", collapse_comma(names)) + } + for (v in names(ls)) { + assign(v, ls[[v]]) + } + } + out <- data.frame( + prior, class, coef, group, resp, dpar, nlpar, bound, + stringsAsFactors = FALSE + ) class(out) <- c("brmsprior", "data.frame") out } empty_brmsprior <- function() { # define a brmsprior object with zero rows - brmsprior(prior = character(0), class = character(0), - coef = character(0), group = character(0), - nlpar = character(0), bound = character(0)) + char0 <- character(0) + brmsprior( + prior = char0, class = char0, coef = char0, + group = char0, resp = char0, dpar = char0, + nlpar = char0, bound = char0 + ) } prior_bounds <- function(prior) { @@ -1334,20 +1357,23 @@ print.brmsprior <- function(x, show_df, ...) { .print_prior <- function(x) { # prepare text for print.brmsprior group <- usc(x$group) + resp <- usc(x$resp) + dpar <- usc(x$dpar) nlpar <- usc(x$nlpar) coef <- usc(x$coef) - nlpar <- ifelse(nchar(group), usc(nlpar), nlpar) - coef <- ifelse(nchar(group) & !nchar(nlpar), usc(coef), coef) - bound <- ifelse(nchar(x$bound), paste0(x$bound, " "), "") - tilde <- ifelse(nchar(x$class) + nchar(group) + nchar(coef), " ~ ", "") - prior <- ifelse(nchar(x$prior), x$prior, "(no prior)") - paste0(bound, x$class, group, nlpar, coef, tilde, prior) + if (any(nzchar(c(resp, dpar, nlpar, coef)))) { + group <- usc(group, "suffix") + } + bound <- ifelse(nzchar(x$bound), paste0(x$bound, " "), "") + tilde <- ifelse(nzchar(x$class) | nzchar(group) | nzchar(coef), " ~ ", "") + prior <- ifelse(nzchar(x$prior), x$prior, "(no prior)") + paste0(bound, x$class, group, resp, dpar, nlpar, coef, tilde, prior) } #' @export c.brmsprior <- function(x, ...) { - # combines multiple brmsprior objects into one brmsprior - if (all(sapply(list(...), is, class2 = "brmsprior"))) { + # combine multiple brmsprior objects into one brmsprior + if (all(sapply(list(...), is.brmsprior))) { out <- do.call(rbind, list(x, ...)) } else { out <- c(as.data.frame(x), ...) @@ -1355,6 +1381,11 @@ c.brmsprior <- function(x, ...) { out } +#' @export +"+.brmsprior" <- function(e1, e2) { + c(e1, e2) +} + dirichlet <- function(...) { # dirichlet prior of simplex parameters out <- as.numeric(c(...)) @@ -1390,7 +1421,7 @@ dirichlet <- function(...) { #' @param autoscale Logical; indicating whether the horseshoe #' prior should be scaled using the residual standard deviation #' \code{sigma} if possible and sensible (defaults to \code{TRUE}). -#' Autoscaling is not applied for auxiliary parameters or +#' Autoscaling is not applied for distributional parameters or #' when the model does not contain the parameter \code{sigma}. #' #' @return A character string obtained by \code{match.call()} with diff --git a/R/rename.R b/R/rename.R index 1078eb9df..58d599942 100644 --- a/R/rename.R +++ b/R/rename.R @@ -17,30 +17,27 @@ rename_pars <- function(x) { change <- list() resp <- bterms$response if (length(resp) > 1L) { - # rename effects in multivaraite models + # rename effects in multivariate models for (r in resp) { + bterms$dpars[["mu"]]$resp <- r change_eff <- change_effects( - bterms$auxpars[["mu"]], model.frame(x), pars, - dims = x$fit@sim$dims_oi, nlpar = r, - stancode = stancode(x) + bterms$dpars[["mu"]], data = model.frame(x), + pars = pars, stancode = stancode(x) ) change <- c(change, change_eff) } - bterms$auxpars[["mu"]] <- NULL + bterms$dpars[["mu"]] <- NULL } - # rename effects of auxilliary parameters - for (ap in names(bterms$auxpars)) { + # rename effects of distributional parameters + for (ap in names(bterms$dpars)) { change_eff <- change_effects( - bterms$auxpars[[ap]], model.frame(x), pars, - dims = x$fit@sim$dims_oi, nlpar = ap, - stancode = stancode(x) + bterms$dpars[[ap]], data = model.frame(x), + pars = pars, stancode = stancode(x) ) change <- c(change, change_eff) } # rename group-level parameters separately - change <- c(change, - change_re(x$ranef, pars = pars, dims = x$fit@sim$dims_oi) - ) + change <- c(change, change_re(x$ranef, pars = pars)) # rename residual parameters of multivariate linear models if (is_linear(family) && length(bterms$response) > 1L) { corfnames <- paste0("sigma_", bterms$response) @@ -74,53 +71,46 @@ rename_pars <- function(x) { } #' @export -change_effects.btl <- function(x, data, pars, dims, nlpar = "", - stancode = "", ...) { +change_effects.btl <- function(x, data, pars, stancode = "", ...) { # helps in renaming various kinds of effects - nlpar <- check_nlpar(nlpar) - fixef <- colnames(data_fe(x, data)$X) - fixef <- rm_int_fe(fixef, stancode, nlpar = nlpar) - change_fe <- change_fe(fixef, pars, nlpar = nlpar) - smooths = get_sm_labels(x, data, covars = TRUE) - change_sm <- change_sm(smooths, pars, nlpar = nlpar) - csef <- colnames(data_cs(x, data)$Xcs) - change_cs <- change_cs(csef, pars, nlpar = nlpar) - monef <- colnames(data_mo(x, data)$Xmo) - change_mo <- change_mo(monef, pars, nlpar = nlpar) - meef <- get_me_labels(x, data) - change_me <- change_me(meef, pars, dims = dims, nlpar = nlpar) - gpef <- get_gp_labels(x, data = data, covars = TRUE) - change_gp <- change_gp(gpef, pars, dims = dims, nlpar = nlpar) - c(change_fe, change_sm, change_cs, change_mo, change_me, change_gp) + # Returns: + # a list whose elements can be interpreted by do_renaming + c( + change_fe(x, data, pars, stancode = stancode), + change_sm(x, data, pars), + change_cs(x, data, pars), + change_mo(x, data, pars), + change_me(x, data, pars), + change_gp(x, data, pars) + ) } -change_effects.btnl <- function(x, data, pars, dims, ...) { +change_effects.btnl <- function(x, data, pars, ...) { # helps in renaming effects for non-linear parameters + # Returns: + # a list whose elements can be interpreted by do_renaming change <- list() for (nlp in names(x$nlpars)) { change <- c(change, - change_effects(x$nlpars[[nlp]], data, pars, dims, nlpar = nlp) + change_effects(x$nlpars[[nlp]], data, pars, ...) ) } change } -change_fe <- function(fixef, pars, nlpar = "") { +change_fe <- function(bterms, data, pars, stancode = "") { # helps in renaming fixed effects parameters - # Args: - # fixef: names of the fixed effects - # pars: names of all model parameters - # nlpar: optional string naming a non-linear parameter # Returns: # a list whose elements can be interpreted by do_renaming change <- list() + px <- check_prefix(bterms) + fixef <- colnames(data_fe(bterms, data)$X) + fixef <- rm_int_fe(fixef, stancode, px = px) if (length(fixef)) { - b <- paste0("b", usc(nlpar, "prefix")) + b <- paste0("b", usc(combine_prefix(px), "prefix")) pos <- grepl(paste0("^", b, "\\["), pars) bnames <- paste0(b, "_", fixef) - change <- lc(change, - list(pos = pos, fnames = bnames) - ) + change <- lc(change, list(pos = pos, fnames = bnames)) change <- c(change, change_prior(class = b, pars = pars, names = fixef) ) @@ -128,17 +118,14 @@ change_fe <- function(fixef, pars, nlpar = "") { change } -change_mo <- function(monef, pars, nlpar = "") { +change_mo <- function(bterms, data, pars) { # helps in renaming monotonous effects parameters - # Args: - # monef: names of the monotonous effects - # pars: names of all model parameters - # nlpar: optional string naming a non-linear parameter # Returns: # a list whose elements can be interpreted by do_renaming change <- list() + monef <- colnames(data_mo(bterms, data)$Xmo) if (length(monef)) { - p <- usc(nlpar, "prefix") + p <- usc(combine_prefix(bterms), "prefix") bmo <- paste0("bmo", p) newnames <- paste0("bmo", p, "_", monef) change <- lc(change, @@ -168,16 +155,13 @@ change_mo <- function(monef, pars, nlpar = "") { change } -change_cs <- function(csef, pars, nlpar = "") { +change_cs <- function(bterms, data, pars) { # helps in renaming category specific effects parameters - # Args: - # csef: names of the category specific effects - # pars: names of all model parameters # Returns: # a list whose elements can be interpreted by do_renaming change <- list() + csef <- colnames(data_cs(bterms, data)$Xcs) if (length(csef)) { - stopifnot(!nzchar(nlpar)) ncse <- length(csef) thres <- sum(grepl("^b_Intercept\\[", pars)) csenames <- t(outer(csef, paste0("[", 1:thres, "]"), FUN = paste0)) @@ -193,11 +177,14 @@ change_cs <- function(csef, pars, nlpar = "") { change } -change_me <- function(meef, pars, dims, nlpar = "") { - # rename parameters of noisy variables +change_me <- function(bterms, data, pars) { + # helps in renaming parameters of noisy variables + # Returns: + # a list whose elements can be interpreted by do_renaming change <- list() + meef <- get_me_labels(bterms, data) if (length(meef)) { - p <- usc(nlpar, "prefix") + p <- usc(combine_prefix(bterms), "prefix") meef <- rename(meef) # rename coefficients of noise free terms bme <- paste0("bme", p) @@ -224,9 +211,13 @@ change_me <- function(meef, pars, dims, nlpar = "") { change } -change_gp <- function(gpef, pars, dims, nlpar = "") { +change_gp <- function(bterms, data, pars) { + # helps in renaming parameters of gaussian processes + # Returns: + # a list whose elements can be interpreted by do_renaming change <- list() - p <- usc(nlpar, "prefix") + gpef <- get_gp_labels(bterms, data = data, covars = TRUE) + p <- usc(combine_prefix(bterms), "prefix") for (i in seq_along(gpef)) { # rename GP hyperparameters by_levels = attr(gpef, "by_levels")[[i]] @@ -263,18 +254,15 @@ change_gp <- function(gpef, pars, dims, nlpar = "") { change } -change_sm <- function(smooths, pars, nlpar = "") { +change_sm <- function(bterms, data, pars) { # helps in renaming smoothing term parameters - # Args: - # smooths: smoothing terms - # pars: names of all model parameters - # nlpar: optional string naming a non-linear parameter # Returns: # a list whose elements can be interpreted by do_renaming change <- list() + smooths <- get_sm_labels(bterms, data, covars = TRUE) if (length(smooths)) { stopifnot(!is.null(attr(smooths, "nbases"))) - p <- usc(nlpar, "prefix") + p <- usc(combine_prefix(bterms), "prefix") sds <- paste0("sds", p) sds_names <- paste0(sds, "_", smooths) s <- paste0("s", p) @@ -304,21 +292,19 @@ change_sm <- function(smooths, pars, nlpar = "") { } -change_re <- function(ranef, pars, dims) { +change_re <- function(ranef, pars) { # helps in renaming group-level parameters # Args: # ranef: list returned by tidy_ranef # pars: names of all model parameters - # dims: named list containing parameter dimensions - # nlpar: optional string naming a non-linear parameter # Returns: # a list whose elements can be interpreted by do_renaming change <- list() if (nrow(ranef)) { for (id in unique(ranef$id)) { - r <- ranef[ranef$id == id, ] + r <- subset2(ranef, id = id) g <- r$group[1] - suffix <- paste0(usc(r$nlpar, "suffix"), r$coef) + suffix <- paste0(usc(combine_prefix(r), "suffix"), r$coef) rfnames <- paste0("sd_", g, "__", suffix) rpos <- grepl(paste0("^sd_", id, "(\\[|$)"), pars) change <- lc(change, list(pos = rpos, fnames = rfnames)) @@ -346,29 +332,27 @@ change_re <- function(ranef, pars, dims) { } } if (any(grepl("^r_", pars))) { - rc_args <- nlist(ranef, pars, dims) - change <- c(change, do.call(change_re_levels, rc_args)) + change <- c(change, change_re_levels(ranef, pars = pars)) } } change } -change_re_levels <- function(ranef, pars, dims) { +change_re_levels <- function(ranef, pars) { # helps in renaming random effects 'r_.' parameters # Args: # ranef: output of tidy_ranef # pars: names of all model parameters - # dims: named list containing parameter dimensions # Returns: # a list whose elements can be interpreted by do_renaming change <- list() for (i in seq_len(nrow(ranef))) { r <- ranef[i, ] - usc_nlpar <- usc(r$nlpar) - r_parnames <- paste0("r_", r$id, usc_nlpar, "_", r$cn) + p <- usc(combine_prefix(r)) + r_parnames <- paste0("r_", r$id, p, "_", r$cn) r_regex <- paste0("^", r_parnames, "(\\[|$)") change_rl <- list(pos = grepl(r_regex, pars)) - r_new_parname <- paste0("r_", r$group, usc(usc_nlpar)) + r_new_parname <- paste0("r_", r$group, usc(p)) # rstan doesn't like whitespaces in parameter names levels <- gsub("[ \t\r\n]", ".", attr(ranef, "levels")[[r$group]]) index_names <- make_index_names(levels, r$coef, dim = 2) @@ -429,7 +413,7 @@ change_old_re <- function(ranef, pars, dims) { # a list whose elements can be interpreted by do_renaming change <- list() for (id in unique(ranef$id)) { - r <- ranef[ranef$id == id, ] + r <- subset2(ranef, id = id) g <- r$group[1] nlpar <- r$nlpar[1] stopifnot(nzchar(nlpar)) @@ -483,7 +467,7 @@ change_old_re2 <- function(ranef, pars, dims) { # a list whose elements can be interpreted by do_renaming change <- list() for (id in unique(ranef$id)) { - r <- ranef[ranef$id == id, ] + r <- subset2(ranef, id = id) g <- r$group[1] nlpars_usc <- usc(r$nlpar, "suffix") # rename sd-parameters @@ -530,12 +514,11 @@ change_old_sm <- function(bterms, pars, dims) { # change names of spline parameters fitted with brms <= 1.0.1 # this became necessary after allowing smooths with multiple covariates stopifnot(is.brmsterms(bterms)) - .change_old_sm <- function(bt, nlpar = "") { - nlpar <- check_nlpar(nlpar) + .change_old_sm <- function(bt) { change <- list() sm_labels <- get_sm_labels(bt) if (length(sm_labels)) { - p <- usc(nlpar, "suffix") + p <- usc(combine_prefix(bt), "suffix") old_smooths <- rename(paste0(p, sm_labels)) new_smooths <- rename(paste0(p, get_sm_labels(bt, covars = TRUE))) old_sds_pars <- paste0("sds_", old_smooths) @@ -562,20 +545,19 @@ change_old_sm <- function(bterms, pars, dims) { change <- list() if (length(bterms$response) > 1L) { for (r in bterms$response) { - change <- c(change, .change_old_sm(bterms$auxpars$mu, nlpar = r)) + bterms$dpars$mu$resp <- r + change <- c(change, .change_old_sm(bterms$dpars$mu)) } - bterms$auxpars$mu <- NULL + bterms$dpars$mu <- NULL } - for (ap in names(bterms$auxpars)) { - bt <- bterms$auxpars[[ap]] + for (dp in names(bterms$dpars)) { + bt <- bterms$dpars[[dp]] if (length(bt$nlpars)) { for (nlp in names(bt$nlpars)) { - change <- c(change, - .change_old_sm(bt$nlpars[[nlp]], nlpar = nlp) - ) + change <- c(change, .change_old_sm(bt$nlpars[[nlp]])) } } else { - change <- c(change, .change_old_sm(bt, nlpar = ap)) + change <- c(change, .change_old_sm(bt)) } } change @@ -594,10 +576,10 @@ change_simple <- function(oldname, fnames, pars, dims, out } -rm_int_fe <- function(fixef, stancode, nlpar = "") { +rm_int_fe <- function(fixef, stancode, px = list()) { # identifies if the intercept has to be removed from fixef # and returns adjusted fixef names - p <- usc(nlpar, "suffix") + p <- usc(combine_prefix(px), "suffix") int <- paste0("b_", p, "Intercept = temp_", p, "Intercept") loclev <- "vector[N] loclev;" if (any(ulapply(c(int, loclev), grepl, stancode, fixed = TRUE))) { @@ -737,7 +719,7 @@ reorder_pars <- function(x) { all_classes <- c( "b", "bmo", "bcs", "bme", "ar", "ma", "arr", "lagsar", "errorsar", "car", "sdcar", "sigmaLL", "sd", "cor", "sds", - "sdgp", "lscale", auxpars(), "temp", "rescor", "delta", + "sdgp", "lscale", dpars(), "temp", "rescor", "delta", "lasso", "simplex", "r", "s", "zgp", "rcar", "loclev", "Xme", "prior", "lp" ) diff --git a/R/stan-helpers.R b/R/stan-helpers.R index 30fd8e71c..5f9335c38 100644 --- a/R/stan-helpers.R +++ b/R/stan-helpers.R @@ -53,7 +53,7 @@ stan_autocor <- function(autocor, bterms, family, prior) { if (is.formula(bterms$adforms$disp)) { stop2(err_msg, " when specifying 'disp'.") } - if (any(c("sigma", "nu") %in% names(bterms$auxpars))) { + if (any(c("sigma", "nu") %in% names(bterms$dpars))) { stop2(err_msg, " when predicting 'sigma' or 'nu'.") } str_add(out$data) <- paste0( @@ -108,7 +108,7 @@ stan_autocor <- function(autocor, bterms, family, prior) { if (is.formula(bterms$adforms$se)) { stop2(err_msg, " when specifying 'se'.") } - if (length(bterms$auxpars[["mu"]]$nlpars)) { + if (length(bterms$dpars[["mu"]]$nlpars)) { stop2(err_msg, " in non-linear models.") } if (!identical(family$link, "identity")) { @@ -148,7 +148,7 @@ stan_autocor <- function(autocor, bterms, family, prior) { if (Karr) { # autoregressive effects of the response err_msg <- "ARR models are not yet working" - if (length(bterms$auxpars[["mu"]]$nlpars)) { + if (length(bterms$dpars[["mu"]]$nlpars)) { stop2(err_msg, " in non-linear models.") } if (is_mv) { @@ -176,7 +176,7 @@ stan_autocor <- function(autocor, bterms, family, prior) { if (is.formula(bterms$adforms$disp)) { stop2(err_msg, " when specifying 'disp'.") } - if (any(c("sigma", "nu") %in% names(bterms$auxpars))) { + if (any(c("sigma", "nu") %in% names(bterms$dpars))) { stop2(err_msg, " when predicting 'sigma' or 'nu'.") } str_add(out$data) <- paste0( @@ -221,7 +221,7 @@ stan_autocor <- function(autocor, bterms, family, prior) { if (is_mv) { stop2(err_msg, " in multivariate models.") } - if (length(bterms$auxpars[["mu"]]$nlpars)) { + if (length(bterms$dpars[["mu"]]$nlpars)) { stop2(err_msg, " in non-linear models.") } str_add(out$data) <- paste0( @@ -288,7 +288,7 @@ stan_autocor <- function(autocor, bterms, family, prior) { if (is_mv) { stop2(err_msg, " in multivariate models.") } - if (length(bterms$auxpars[["mu"]]$nlpars)) { + if (length(bterms$dpars[["mu"]]$nlpars)) { stop2(err_msg, " in non-linear models.") } str_add(out$data) <- @@ -575,7 +575,7 @@ stan_families <- function(family, bterms) { # as suggested by Stephen Martin use sigma and mu of CP # but the skewness parameter alpha of DP str_add(out$tdataD) <- " real sqrt_2_div_pi = sqrt(2 / pi()); \n" - ap_names <- names(bterms$auxpars) + ap_names <- names(bterms$dpars) for (i in which(families %in% "skew_normal")) { id <- ifelse(length(families) == 1L, "", i) ns <- ifelse(paste0("sigma", id) %in% ap_names, "[n]", "") @@ -611,7 +611,7 @@ stan_families <- function(family, bterms) { " #include 'fun_gen_extreme_value.stan' \n", " #include 'fun_scale_xi.stan' \n" ) - ap_names <- c(names(bterms$auxpars), names(bterms$fauxpars)) + ap_names <- c(names(bterms$dpars), names(bterms$fdpars)) for (i in which(families %in% "gen_extreme_value")) { id <- ifelse(length(families) == 1L, "", i) xi <- paste0("xi", id) @@ -620,7 +620,7 @@ stan_families <- function(family, bterms) { " real ", xi, "; // scaled shape parameter \n" ) sigma <- paste0("sigma", id) - v <- ifelse(sigma %in% names(bterms$auxpars), "_vector", "") + v <- ifelse(sigma %in% names(bterms$dpars), "_vector", "") args <- sargs(paste0("temp_", xi), "Y", paste0("mu", id), sigma) str_add(out$modelC) <- paste0( " ", xi, " = scale_xi", v, "(", args, "); \n" @@ -636,10 +636,10 @@ stan_mixture <- function(bterms, prior) { out <- list() if (is.mixfamily(bterms$family)) { nmix <- length(bterms$family$mix) - theta_pred <- grepl("^theta", names(bterms$auxpars)) - theta_pred <- bterms$auxpars[theta_pred] - theta_fix <- grepl("^theta", names(bterms$fauxpars)) - theta_fix <- bterms$fauxpars[theta_fix] + theta_pred <- grepl("^theta", names(bterms$dpars)) + theta_pred <- bterms$dpars[theta_pred] + theta_fix <- grepl("^theta", names(bterms$fdpars)) + theta_fix <- bterms$fdpars[theta_fix] def_thetas <- collapse( " real theta", 1:nmix, ";", " // mixing proportion \n" @@ -648,7 +648,7 @@ stan_mixture <- function(bterms, prior) { if (length(theta_pred) != nmix - 1) { stop2("Can only predict all but one mixing proportion.") } - missing_id <- setdiff(1:nmix, auxpar_id(names(theta_pred))) + missing_id <- setdiff(1:nmix, dpar_id(names(theta_pred))) str_add(out$modelD) <- paste0( " vector[N] theta", missing_id, " = rep_vector(0, N); \n", " real log_sum_exp_theta; \n" @@ -795,7 +795,7 @@ stan_misc_functions <- function(family, prior, kronecker) { } stan_prior <- function(prior, class, coef = "", group = "", - nlpar = "", prefix = "", suffix = "", + px = list(), prefix = "", suffix = "", wsp = 2, matrix = FALSE) { # Define priors for parameters in Stan language # Args: @@ -816,32 +816,36 @@ stan_prior <- function(prior, class, coef = "", group = "", tp <- tp(wsp) wsp <- wsp(nsp = wsp) prior_only <- identical(attr(prior, "sample_prior"), "only") - keep <- prior$class == class & - prior$coef %in% c(coef, "") & prior$group %in% c(group, "") + prior <- subset2(prior, + class = class, coef = c(coef, ""), group = c(group, "") + ) if (class %in% c("sd", "cor")) { # only sd and cor parameters have global priors - keep <- keep & prior$nlpar %in% c(nlpar, "") + px_tmp <- lapply(px, function(x) c(x, "")) + prior <- subset2(prior, ls = px_tmp) } else { - keep <- keep & prior$nlpar %in% nlpar + prior <- subset2(prior, ls = px) } - prior <- prior[keep, ] if (!nchar(class) && nrow(prior)) { # unchecked prior statements are directly passed to Stan return(collapse(wsp, prior$prior, "; \n")) } - unlpar <- unique(nlpar) - if (length(unlpar) > 1L) { + px <- as.data.frame(px) + upx <- unique(px) + if (nrow(upx) > 1L) { # can only happen for SD parameters of the same ID - base_prior <- rep(NA, length(unlpar)) - for (i in seq_along(unlpar)) { - nlpar_prior <- prior[prior$nlpar %in% c("", unlpar[i]), ] - base_prior[i] <- stan_base_prior(nlpar_prior) + base_prior <- rep(NA, nrow(upx)) + for (i in seq_len(nrow(upx))) { + sub_upx <- lapply(upx[i, ], function(x) c(x, "")) + sub_prior <- subset2(prior, ls = sub_upx) + base_prior[i] <- stan_base_prior(sub_prior) } if (length(unique(base_prior)) > 1L) { # define prior for single coefficients manually # as there is not single base_prior anymore - take <- match(prior[nzchar(prior$coef), "nlpar"], unlpar) + prior_of_coefs <- prior[nzchar(prior$coef), vars_prefix()] + take <- match_rows(prior_of_coefs, upx) prior[nzchar(prior$coef), "prior"] <- base_prior[take] } base_prior <- base_prior[1] @@ -851,15 +855,15 @@ stan_prior <- function(prior, class, coef = "", group = "", bound <- prior[!nzchar(prior$coef), "bound"] } - individual_prior <- function(i, max_index) { + individual_prior <- function(i, prior, max_index) { # individual priors for each parameter of a class if (max_index > 1L || matrix) { index <- paste0("[", i, "]") } else { index <- "" } - if (length(nlpar) > 1L) { - prior <- prior[prior$nlpar == nlpar[i], ] + if (nrow(px) > 1L) { + prior <- subset2(prior, ls = px[i, ]) } uc_prior <- prior$prior[match(coef[i], prior$coef)] if (!is.na(uc_prior) & nchar(uc_prior)) { @@ -885,8 +889,10 @@ stan_prior <- function(prior, class, coef = "", group = "", class <- paste0(prefix, class, suffix) if (any(with(prior, nchar(coef) & nchar(prior)))) { # generate a prior for each coefficient - out <- sapply(seq_along(coef), individual_prior, - max_index = length(coef)) + out <- sapply( + seq_along(coef), individual_prior, + prior = prior, max_index = length(coef) + ) } else if (nchar(base_prior) > 0) { if (matrix) { class <- paste0("to_vector(", class, ")") @@ -899,7 +905,7 @@ stan_prior <- function(prior, class, coef = "", group = "", out <- "" } special_prior <- stan_special_prior( - class, prior, ncoef = length(coef), nlpar = nlpar + class, prior, ncoef = length(coef), px = px ) out <- collapse(c(out, special_prior)) if (prior_only && nzchar(class) && !nchar(out)) { @@ -915,21 +921,19 @@ stan_base_prior <- function(prior) { # Args: # prior: a prior.frame stopifnot(length(unique(prior$class)) <= 1L) - igroup <- which(with(prior, !nchar(coef) & nchar(group) & nchar(prior))) - inlpar <- which(with(prior, !nchar(coef) & nchar(nlpar) & nchar(prior))) - iclass <- which(with(prior, !nchar(coef) & !nchar(group) & nchar(prior))) - if (length(igroup)) { - # if there is a global prior for this group - base_prior <- prior[igroup, "prior"] - } else if (length(inlpar)) { - # if there is a global prior for this non-linear parameter - base_prior <- prior[inlpar, "prior"] - } else if (length(iclass)) { - # if there is a global prior for this class - base_prior <- prior[iclass, "prior"] - } else { - # no proper prior for this class - base_prior <- "" + prior <- prior[with(prior, !nzchar(coef) & nzchar(prior)), ] + vars <- c("group", "nlpar", "dpar", "resp", "class") + i <- 1 + found <- FALSE + base_prior <- "" + take <- rep(FALSE, nrow(prior)) + while (!found && i <= length(vars)) { + take <- nzchar(prior[[vars[i]]]) & !take + if (any(take)) { + base_prior <- prior[take, "prior"] + found <- TRUE + } + i <- i + 1 } stopifnot(length(base_prior) == 1L) base_prior @@ -979,16 +983,16 @@ stan_target_prior <- function(prior, par, ncoef = 1, bound = "") { out } -stan_special_prior <- function(class, prior, ncoef, nlpar = "") { +stan_special_prior <- function(class, prior, ncoef, px = list()) { # add special priors such as horseshoe and lasso out <- "" - p <- usc(nlpar) + p <- usc(combine_prefix(px)) if (all(class == paste0("b", p))) { - stopifnot(length(nlpar) == 1L) + stopifnot(length(p) == 1L) tp <- tp() # add horseshoe and lasso shrinkage priors - orig_nlpar <- ifelse(nzchar(nlpar), nlpar, "mu") - special <- attributes(prior)$special[[orig_nlpar]] + prefix <- combine_prefix(px, keep_mu = TRUE) + special <- attributes(prior)$special[[prefix]] if (!is.null(special$hs_df)) { local_args <- paste0("0.5 * hs_df", p) local_args <- sargs(local_args, local_args) @@ -1108,11 +1112,11 @@ stan_rngprior <- function(sample_prior, prior, par_declars, # use parameters sampled from priors for use in other priors spars <- NULL # cannot sample from the horseshoe prior anymore as of brms 1.5.0 - lasso_nlpars <- nzchar(ulapply(prior_special, "[[", "lasso_df")) - lasso_nlpars <- names(prior_special)[lasso_nlpars] - lasso_nlpars <- usc(ulapply(lasso_nlpars, check_nlpar)) - if (length(lasso_nlpars)) { - spars <- c(spars, paste0("lasso_inv_lambda", lasso_nlpars)) + lasso_prefix <- nzchar(ulapply(prior_special, "[[", "lasso_df")) + lasso_prefix <- names(prior_special)[lasso_prefix] + lasso_prefix <- usc(sub("^mu(_|$)", "", lasso_prefix)) + if (length(lasso_prefix)) { + spars <- c(spars, paste0("lasso_inv_lambda", lasso_prefix)) } if (length(spars)) { bpars <- grepl("^b(|mo|cs|me)(_|$)", pars) @@ -1196,7 +1200,7 @@ stan_has_built_in_fun <- function(family) { # family: a list with elements 'family' and 'link' stopifnot(all(c("family", "link") %in% names(family))) link <- family$link - par <- family$par + dpar <- family$dpar family <- family$family log_families <- c( "poisson", "negbinomial", "geometric", @@ -1210,7 +1214,7 @@ stan_has_built_in_fun <- function(family) { isTRUE( family %in% log_families && link == "log" || family %in% logit_families && link == "logit" || - isTRUE(par %in% c("zi", "hu")) && link == "logit" + isTRUE(dpar %in% c("zi", "hu")) && link == "logit" ) } diff --git a/R/stan-likelihood.R b/R/stan-likelihood.R index 0bc80245f..6e2eaf3a6 100644 --- a/R/stan-likelihood.R +++ b/R/stan-likelihood.R @@ -32,30 +32,30 @@ stan_llh.default <- function(family, bterms, data, autocor, has_trunc <- any(bounds$lb > -Inf) || any(bounds$ub < Inf) llh_adj <- stan_llh_adj(bterms$adforms) - auxpars <- names(bterms$auxpars) + dpars <- names(bterms$dpars) reqn <- llh_adj || is_categorical || is_ordinal || is_hurdle || is_zero_inflated || is_mix || is_zero_one_inflated(family) || is_wiener(family) || is_exgaussian(family) || is_asym_laplace(family) || is_gev(family) || has_sigma && has_se && !use_cov(autocor) || - any(c("phi", "kappa") %in% auxpars) + any(c("phi", "kappa") %in% dpars) n <- ifelse(reqn, "[n]", "") - # prepare auxiliary parameters - p <- named_list(auxpars()) + # prepare distributional parameters + p <- named_list(dpars()) p$mu <- paste0(ifelse(is_mv || is_categorical, "Mu", "mu"), mix, n) p$sigma <- stan_llh_sigma(family, bterms, mix) p$shape <- stan_llh_shape(family, bterms, mix) - for (ap in setdiff(auxpars(), c("mu", "sigma", "shape"))) { - p[[ap]] <- paste0(ap, mix, if (reqn && ap %in% auxpars) "[n]") + for (ap in setdiff(dpars(), c("mu", "sigma", "shape"))) { + p[[ap]] <- paste0(ap, mix, if (reqn && ap %in% dpars) "[n]") } if (family == "skew_normal") { # required because of CP parameterization of mu and sigma - nomega <- if (reqn && any(c("sigma", "alpha") %in% auxpars)) "[n]" + nomega <- if (reqn && any(c("sigma", "alpha") %in% dpars)) "[n]" p$omega <- paste0("omega", mix, nomega) } ord_args <- sargs(p$mu, if (has_cs) "mucs[n]", "temp_Intercept", p$disc) - usc_logit <- stan_llh_auxpar_usc_logit(c("zi", "hu"), bterms) + usc_logit <- stan_llh_dpar_usc_logit(c("zi", "hu"), bterms) trials <- ifelse(llh_adj || is_zero_inflated, "trials[n]", "trials") if (is_mv) { @@ -84,7 +84,7 @@ stan_llh.default <- function(family, bterms, data, autocor, } simplify <- stan_has_built_in_fun(nlist(family, link)) && - !has_trunc && !has_cens && !"disc" %in% auxpars + !has_trunc && !has_cens && !"disc" %in% dpars if (simplify) { llh_pre <- switch(family, poisson = c( @@ -337,14 +337,14 @@ stan_llh.default <- function(family, bterms, data, autocor, #' @export stan_llh.mixfamily <- function(family, bterms, ...) { - ap_ids <- auxpar_id(names(bterms$auxpars)) - fap_ids <- auxpar_id(names(bterms$fauxpars)) - ptheta <- any(auxpar_class(names(bterms$auxpars)) %in% "theta") + ap_ids <- dpar_id(names(bterms$dpars)) + fap_ids <- dpar_id(names(bterms$fdpars)) + ptheta <- any(dpar_class(names(bterms$dpars)) %in% "theta") llh <- rep(NA, length(family$mix)) for (i in seq_along(family$mix)) { sbterms <- bterms - sbterms$auxpars <- sbterms$auxpars[ap_ids == i] - sbterms$fauxpars <- sbterms$fauxpars[fap_ids == i] + sbterms$dpars <- sbterms$dpars[ap_ids == i] + sbterms$fdpars <- sbterms$fdpars[fap_ids == i] llh[i] <- stan_llh( family$mix[[i]], sbterms, mix = i, ptheta = ptheta, ... ) @@ -470,10 +470,10 @@ stan_llh_sigma <- function(family, bterms, mix = "") { has_se <- is.formula(bterms$adforms$se) has_disp <- is.formula(bterms$adforms$disp) llh_adj <- stan_llh_adj(bterms$adforms) - auxpars <- names(bterms$auxpars) + dpars <- names(bterms$dpars) nsigma <- llh_adj || has_se || nzchar(mix) || family %in% c("exgaussian", "gen_extreme_value", "asym_laplace") - nsigma <- nsigma && (has_disp || paste0("sigma", mix) %in% auxpars) + nsigma <- nsigma && (has_disp || paste0("sigma", mix) %in% dpars) nsigma <- if (nsigma) "[n]" nse <- if (llh_adj) "[n]" if (has_sigma) { @@ -496,19 +496,19 @@ stan_llh_shape <- function(family, bterms, mix = "") { # prepare the code for 'shape' in the likelihood statement has_disp <- is.formula(bterms$adforms$disp) llh_adj <- stan_llh_adj(bterms$adforms) - auxpars <- names(bterms$auxpars) + dpars <- names(bterms$dpars) nshape <- (llh_adj || is_forked(family) || nzchar(mix)) && - (has_disp || paste0("shape", mix) %in% auxpars) + (has_disp || paste0("shape", mix) %in% dpars) nshape <- if (nshape) "[n]" paste0(if (has_disp) "disp_", "shape", mix, nshape) } -stan_llh_auxpar_usc_logit <- function(auxpars, bterms) { - # prepare _logit suffix for auxiliary parameters +stan_llh_dpar_usc_logit <- function(dpars, bterms) { + # prepare _logit suffix for distributional parameters # currently only used in zero-inflated and hurdle models stopifnot(is.brmsterms(bterms)) - use_logit <- any(ulapply(auxpars, function(ap) - isTRUE(bterms$auxpars[[ap]]$family$link == "logit") + use_logit <- any(ulapply(dpars, function(dp) + isTRUE(bterms$dpars[[dp]]$family$link == "logit") )) ifelse(use_logit, "_logit", "") } diff --git a/R/stan-predictor.R b/R/stan-predictor.R index 80d63abfb..4e1485f85 100644 --- a/R/stan-predictor.R +++ b/R/stan-predictor.R @@ -1,52 +1,31 @@ stan_effects.btl <- function(x, data, ranef, prior, center_X = TRUE, - sparse = FALSE, nlpar = "", eta = "mu", - ilink = rep("", 2), order_mixture = 'none', - ...) { + sparse = FALSE, ilink = rep("", 2), + order_mixture = 'none', ...) { # combine effects for the predictors of a single (non-linear) parameter # Args: # center_X: center population-level design matrix if possible? - # eta: prefix of the linear predictor variable - nlpar <- check_nlpar(nlpar) - if (nzchar(eta) && nzchar(nlpar)) { - eta <- usc(eta, "suffix") - } - eta <- paste0(eta, nlpar) - stopifnot(nzchar(eta)) + # sparse: should the population-level design matrix be treated as sparse? + # ilink: character vector of lenght 2 defining the link to be applied + # order_mixture: indicates how to identify mixture models via ordering stopifnot(length(ilink) == 2L) - - ranef <- ranef[ranef$nlpar == nlpar, ] - out <- list() - # include population-level effects + px <- check_prefix(x) + eta <- combine_prefix(px, keep_mu = TRUE) + stopifnot(nzchar(eta)) + ranef <- subset2(ranef, ls = px) center_X <- center_X && has_intercept(x$fe) && !is(x$autocor, "cor_bsts") && !sparse - rm_int <- center_X || is(x$autocor, "cor_bsts") || is_ordinal(x$family) - cols2remove <- if (rm_int) "Intercept" - fixef <- setdiff(colnames(data_fe(x, data)$X), cols2remove) - text_fe <- stan_fe( - fixef, center_X = center_X, family = x$family, - prior = prior, nlpar = nlpar, sparse = sparse, - order_mixture = order_mixture - ) - # include smooth terms - smooths <- get_sm_labels(x, data = data) - text_sm <- stan_sm(smooths, prior = prior, nlpar = nlpar) - # include category specific effects - csef <- colnames(get_model_matrix(x$cs, data)) - text_cs <- stan_cs(csef, ranef, prior = prior) - # include monotonic effects - monef <- all_terms(x$mo) - text_mo <- stan_mo(monef, ranef, prior = prior, nlpar = nlpar) - # include measurement error variables - meef <- get_me_labels(x, data = data) - text_me <- stan_me(meef, ranef, prior = prior, nlpar = nlpar) - # include gaussian processes - gpef <- get_gp_labels(x, data = data) - text_gp <- stan_gp(gpef, prior = prior, nlpar = nlpar) - out <- collapse_lists( - out, text_fe, text_cs, text_mo, text_me, text_sm, text_gp + text_fe <- stan_fe( + x, data, prior = prior, center_X = center_X, + sparse = sparse, order_mixture = order_mixture + ), + text_cs <- stan_cs(x, data, ranef = ranef, prior = prior), + text_mo <- stan_mo(x, data, ranef = ranef, prior = prior), + text_me <- stan_me(x, data, ranef = ranef, prior = prior), + text_sm <- stan_sm(x, data, prior = prior), + text_gp <- stan_gp(x, data, prior = prior) ) - p <- usc(nlpar, "prefix") + p <- usc(combine_prefix(px)) if (is.formula(x$offset)) { str_add(out$data) <- paste0( " vector[N] offset", p, "; \n" @@ -67,9 +46,9 @@ stan_effects.btl <- function(x, data, ranef, prior, center_X = TRUE, # repare loop over eta eta_loop <- paste0( - stan_eta_re(ranef, nlpar = nlpar), + stan_eta_re(ranef, px = px), text_mo$eta, text_me$eta, - stan_eta_ma(x$autocor, nlpar = p), + stan_eta_ma(x$autocor, px = px), stan_eta_car(x$autocor), stan_eta_bsts(x$autocor) ) @@ -88,7 +67,7 @@ stan_effects.btl <- function(x, data, ranef, prior, center_X = TRUE, # possibly transform eta before it is passed to the likelihood if (sum(nzchar(ilink))) { # make sure mu comes last as it might depend on other parameters - not_mu <- nzchar(nlpar) && auxpar_class(nlpar) != "mu" + not_mu <- nzchar(x$dpar) && dpar_class(x$dpar) != "mu" position <- ifelse(not_mu, "modelC3", "modelC4") str_add(out[[position]]) <- paste0( " ", eta, "[n] = ", ilink[1], eta, "[n]", ilink[2], "; \n" @@ -97,8 +76,8 @@ stan_effects.btl <- function(x, data, ranef, prior, center_X = TRUE, out } -stan_effects.btnl <- function(x, data, ranef, prior, eta = "mu", - nlpar = "", ilink = rep("", 2), ...) { +stan_effects.btnl <- function(x, data, ranef, prior, + ilink = rep("", 2), ...) { # prepare Stan code for non-linear models # Args: # data: data.frame supplied by the user @@ -106,95 +85,97 @@ stan_effects.btnl <- function(x, data, ranef, prior, eta = "mu", # prior: a brmsprior object # nlpar: currently unused but should not be part of ... # ilink: character vector of length 2 containing - # Stan code for the link function + # Stan code for the link function # ...: passed to stan_effects.btl stopifnot(length(ilink) == 2L) out <- list() - if (length(x$nlpars)) { - nlpars <- names(x$nlpars) - for (nlp in nlpars) { - nl_text <- stan_effects( - x = x$nlpars[[nlp]], data = data, - ranef = ranef, prior = prior, - nlpar = nlp, center_X = FALSE, - ... - ) - out <- collapse_lists(out, nl_text) - } - # prepare non-linear model - new_nlpars <- paste0(" ", eta, "_", nlpars, "[n] ") - # covariates in the non-linear model - covars <- wsp(setdiff(all.vars(rhs(x$formula)), nlpars)) - if (length(covars)) { - # use vectors as indexing matrices in Stan is slow - str_add(out$data) <- paste0( - " // covariate vectors \n", - collapse(" vector[N] C_", seq_along(covars), ";\n") - ) - new_covars <- paste0(" C_", seq_along(covars), "[n] ") - } else { - new_covars <- NULL - } - # add whitespaces to be able to replace parameters and covariates - meta_sym <- c("+", "-", "*", "/", "^", ")", "(", ",") - nlmodel <- gsub(" ", "", collapse(deparse(x$formula[[2]]))) - nlmodel <- wsp(rename(nlmodel, meta_sym, wsp(meta_sym))) - nlmodel <- rename(nlmodel, - c(wsp(nlpars), covars, " ( ", " ) "), - c(new_nlpars, new_covars, "(", ")") + if (!length(x$nlpars)) { + return(out) + } + nlpars <- names(x$nlpars) + for (nlp in nlpars) { + nl_text <- stan_effects( + x = x$nlpars[[nlp]], data = data, + ranef = ranef, prior = prior, + center_X = FALSE, ... ) - # possibly transform eta in the transformed params block - str_add(out$modelD) <- paste0(" vector[N] ", eta, "; \n") - str_add(out$modelC4) <- paste0( - " // compute non-linear predictor \n", - " ", eta, "[n] = ", ilink[1], trimws(nlmodel), ilink[2], "; \n" + out <- collapse_lists(out, nl_text) + } + # prepare non-linear model + par <- combine_prefix(x, keep_mu = TRUE) + new_nlpars <- paste0(" ", par, "_", nlpars, "[n] ") + # covariates in the non-linear model + covars <- wsp(setdiff(all.vars(rhs(x$formula)), nlpars)) + if (length(covars)) { + # use vectors as indexing matrices in Stan is slow + p <- usc(combine_prefix(x), "suffix") + str_add(out$data) <- paste0( + " // covariate vectors \n", + collapse(" vector[N] C_", p, seq_along(covars), ";\n") ) + new_covars <- paste0(" C_", p, seq_along(covars), "[n] ") + } else { + new_covars <- NULL } + # add whitespaces to be able to replace parameters and covariates + meta_sym <- c("+", "-", "*", "/", "^", ")", "(", ",") + nlmodel <- gsub(" ", "", collapse(deparse(x$formula[[2]]))) + nlmodel <- wsp(rename(nlmodel, meta_sym, wsp(meta_sym))) + nlmodel <- rename(nlmodel, + c(wsp(nlpars), covars, " ( ", " ) "), + c(new_nlpars, new_covars, "(", ")") + ) + # possibly transform eta in the transformed params block + str_add(out$modelD) <- paste0(" vector[N] ", par, "; \n") + str_add(out$modelC4) <- paste0( + " // compute non-linear predictor \n", + " ", par, "[n] = ", ilink[1], trimws(nlmodel), ilink[2], "; \n" + ) out } stan_effects.brmsterms <- function(x, data, ranef, prior, sparse = FALSE, ...) { - # Stan code for auxiliary parameters + # Stan code for distributional parameters # Args: # bterms: object of class brmsterms # other arguments: same as make_stancode out <- list() - valid_auxpars <- valid_auxpars(x$family, bterms = x) + valid_dpars <- valid_dpars(x$family, bterms = x) args <- nlist(data, ranef, prior) - for (ap in valid_auxpars) { - ap_terms <- x$auxpars[[ap]] + for (dp in valid_dpars) { + ap_terms <- x$dpars[[dp]] if (is.btl(ap_terms) || is.btnl(ap_terms)) { ilink <- stan_eta_ilink( - ap_terms$family, auxpars = names(x$auxpars), - adforms = x$adforms, mix = auxpar_id(ap) + ap_terms$family, dpars = names(x$dpars), + adforms = x$adforms, mix = dpar_id(dp) ) - eta <- ifelse(ap == "mu", "mu", "") + eta <- ifelse(dp == "mu", "mu", "") ap_args <- list( - ap_terms, nlpar = ap, eta = eta, ilink = ilink, + ap_terms, eta = eta, ilink = ilink, sparse = sparse, order_mixture = x$family$order ) - out[[ap]] <- do.call(stan_effects, c(ap_args, args)) - } else if (is.numeric(x$fauxpars[[ap]])) { - out[[ap]] <- list(data = stan_auxpar_defs(ap)) - } else if (is.character(x$fauxpars[[ap]])) { - if (!x$fauxpars[[ap]] %in% valid_auxpars) { - stop2("Parameter '", x$fauxpars[[ap]], "' cannot be found.") + out[[dp]] <- do.call(stan_effects, c(ap_args, args)) + } else if (is.numeric(x$fdpars[[dp]]$value)) { + out[[dp]] <- list(data = stan_dpar_defs(dp)) + } else if (is.character(x$fdpars[[dp]]$value)) { + if (!x$fdpars[[dp]]$value %in% valid_dpars) { + stop2("Parameter '", x$fdpars[[dp]]$value, "' cannot be found.") } - out[[ap]] <- list( - tparD = stan_auxpar_defs(ap), - tparC1 = paste0(" ", ap, " = ", x$fauxpars[[ap]], "; \n") + out[[dp]] <- list( + tparD = stan_dpar_defs(dp), + tparC1 = paste0(" ", dp, " = ", x$fdpars[[dp]]$value, "; \n") ) } else { - def_temp <- stan_auxpar_defs_temp(ap) - def <- stan_auxpar_defs(ap) + def_temp <- stan_dpar_defs_temp(dp) + def <- stan_dpar_defs(dp) if (nzchar(def_temp)) { - out[[ap]] <- list(par = def_temp, - prior = stan_prior(prior, class = ap, prefix = "temp_") + out[[dp]] <- list(par = def_temp, + prior = stan_prior(prior, class = dp, prefix = "temp_") ) } else if (nzchar(def)) { - out[[ap]] <- list(par = def, - prior = stan_prior(prior, class = ap) + out[[dp]] <- list(par = def, + prior = stan_prior(prior, class = dp) ) } } @@ -211,13 +192,14 @@ stan_effects_mv <- function(bterms, data, ranef, prior, sparse = FALSE) { resp <- bterms$response if (length(resp) > 1L) { args <- nlist( - x = bterms$auxpars[["mu"]], data, ranef, prior, sparse, + x = bterms$dpars[["mu"]], data, ranef, prior, sparse, ilink = stan_eta_ilink(bterms$family, adforms = bterms$adforms) ) resp <- bterms$response tmp_list <- named_list(resp) for (r in resp) { - tmp_list[[r]] <- do.call(stan_effects, c(args, nlpar = r)) + args$x$resp <- r + tmp_list[[r]] <- do.call(stan_effects, args) } out <- collapse_lists(ls = tmp_list) if (is_linear(bterms$family)) { @@ -239,22 +221,24 @@ stan_effects_mv <- function(bterms, data, ranef, prior, sparse = FALSE) { out } -stan_fe <- function(fixef, prior, family = gaussian(), - center_X = TRUE, nlpar = "", sparse = FALSE, - order_mixture = 'none') { +stan_fe <- function(bterms, data, prior, center_X = TRUE, + sparse = FALSE, order_mixture = 'none') { # Stan code for population-level effects # Args: - # fixef: names of the population-level effects # center_X: center the design matrix? - # family: the model family - # prior: a data.frame containing user defined priors - # as returned by check_prior + # sparse: should the design matrix be treated as sparse? # order_mixture: order intercepts to identify mixture models? # Returns: # a list containing Stan code related to population-level effects - p <- usc(nlpar, "prefix") - ct <- ifelse(center_X, "c", "") out <- list() + fixef <- colnames(data_fe(bterms, data)$X) + rm_intercept <- center_X || is.cor_bsts(bterms$autocor) || + is_ordinal(bterms$family) + if (rm_intercept) { + fixef <- setdiff(fixef, "Intercept") + } + px <- check_prefix(bterms) + p <- usc(combine_prefix(px)) if (length(fixef)) { str_add(out$data) <- paste0( " int K", p, ";", @@ -275,8 +259,9 @@ stan_fe <- function(fixef, prior, family = gaussian(), ) } # prepare population-level coefficients - orig_nlpar <- ifelse(nzchar(nlpar), nlpar, "mu") - special <- attr(prior, "special")[[orig_nlpar]] + ct <- ifelse(center_X, "c", "") + prefix <- combine_prefix(px, keep_mu = TRUE) + special <- attr(prior, "special")[[prefix]] if (!is.null(special[["hs_df"]])) { str_add(out$data) <- paste0( " real hs_df", p, "; \n", @@ -307,7 +292,7 @@ stan_fe <- function(fixef, prior, family = gaussian(), " = horseshoe(", hs_args, "); \n" ) } else { - bound <- get_bound(prior, class = "b", nlpar = nlpar) + bound <- get_bound(prior, class = "b", px = px) str_add(out$par) <- paste0( " vector", bound, "[K", ct, p, "] b", p, ";", " // population-level effects \n" @@ -324,8 +309,7 @@ stan_fe <- function(fixef, prior, family = gaussian(), ) } str_add(out$prior) <- stan_prior( - prior, class = "b", coef = fixef, - nlpar = nlpar, suffix = p + prior, class = "b", coef = fixef, px = px, suffix = p ) } if (center_X) { @@ -344,7 +328,7 @@ stan_fe <- function(fixef, prior, family = gaussian(), " } \n" ) # cumulative and sratio models are parameterized as thres - eta - use_plus <- family$family %in% c("cumulative", "sratio") + use_plus <- bterms$family$family %in% c("cumulative", "sratio") sub_X_means <- paste0( ifelse(use_plus, " + ", " - "), "dot_product(means_X", p, ", b", p, ")" @@ -352,7 +336,7 @@ stan_fe <- function(fixef, prior, family = gaussian(), } else { sub_X_means <- "" } - if (is_ordinal(family)) { + if (is_ordinal(bterms$family)) { # temp intercepts for ordinal models are defined in stan_ordinal str_add(out$genD) <- paste0( " // compute actual thresholds \n", @@ -360,9 +344,9 @@ stan_fe <- function(fixef, prior, family = gaussian(), " = temp_Intercept", sub_X_means, "; \n" ) } else { - if (identical(auxpar_class(nlpar), order_mixture)) { + if (identical(dpar_class(px$dpar), order_mixture)) { # identify mixtures via ordering of the intercepts - ap_id <- auxpar_id(nlpar) + ap_id <- dpar_id(px$dpar) str_add(out$tparD) <- paste0( " // identify mixtures via ordering of the intercepts \n", " real temp", p, "_Intercept", @@ -381,19 +365,21 @@ stan_fe <- function(fixef, prior, family = gaussian(), } # for equidistant thresholds only temp_Intercept1 is a parameter prefix <- paste0("temp", p, "_") - suffix <- ifelse(is_equal(family$threshold, "equidistant"), "1", "") + suffix <- ifelse( + is_equal(bterms$family$threshold, "equidistant"), "1", "" + ) str_add(out$prior) <- stan_prior( - prior, class = "Intercept", nlpar = nlpar, + prior, class = "Intercept", px = px, prefix = prefix, suffix = suffix ) } else { - if (identical(auxpar_class(nlpar), order_mixture)) { + if (identical(dpar_class(px$dpar), order_mixture)) { stop2("Identifying mixture components via ordering requires ", "population-level intercepts to be present.\n", "Try setting order = 'none' in function 'mixture'.") } } - out$eta <- stan_eta_fe(fixef, center_X, sparse, nlpar) + out$eta <- stan_eta_fe(fixef, center_X, sparse, px = px) out } @@ -406,11 +392,12 @@ stan_re <- function(id, ranef, prior, cov_ranef = NULL) { # cov_ranef: a list of custom covariance matrices # Returns: # A list of strings containing Stan code - r <- ranef[ranef$id == id, ] + out <- list() + r <- subset2(ranef, id = id) ccov <- r$group[1] %in% names(cov_ranef) ng <- seq_along(r$gcall[[1]]$groups) - idp <- paste0(r$id, usc(r$nlpar, "prefix")) - out <- list() + px <- check_prefix(r) + idp <- paste0(r$id, usc(combine_prefix(px))) str_add(out$data) <- paste0( " // data for group-level effects of ID ", id, " \n", if (r$gtype[1] == "mm") { @@ -429,8 +416,8 @@ stan_re <- function(id, ranef, prior, cov_ranef = NULL) { ) ) str_add(out$prior) <- stan_prior( - prior, class = "sd", group = r$group[1], coef = r$coef, - nlpar = r$nlpar, suffix = paste0("_", id) + prior, class = "sd", group = r$group[1], coef = r$coef, + px = px, suffix = paste0("_", id) ) J <- seq_len(nrow(r)) has_def_type <- !r$type %in% c("mo", "me") @@ -518,16 +505,12 @@ stan_re <- function(id, ranef, prior, cov_ranef = NULL) { out } -stan_sm <- function(smooths, prior, nlpar = "") { +stan_sm <- function(bterms, data, prior) { # Stan code of smooth terms - # Args: - # smooths: names of the smooth terms - # prior: object of class brmsprior - # nlpar: optional name of a non-linear parameter - # Returns: - # A list of strings containing Stan code out <- list() - p <- usc(nlpar) + smooths <- get_sm_labels(bterms, data = data) + px <- check_prefix(bterms) + p <- usc(combine_prefix(px)) if (length(smooths)) { stopifnot(!is.null(attr(smooths, "nbases"))) for (i in seq_along(smooths)) { @@ -558,22 +541,20 @@ stan_sm <- function(smooths, prior, nlpar = "") { " target += normal_lpdf(zs", pi, "_", nb, " | 0, 1); \n" ), stan_prior(prior, class = "sds", coef = smooths[i], - nlpar = nlpar, suffix = paste0(pi, "_", nb)) + px = px, suffix = paste0(pi, "_", nb)) ) } - out$eta <- stan_eta_sm(smooths, nlpar = nlpar) + out$eta <- stan_eta_sm(smooths, px = px) } out } -stan_mo <- function(monef, ranef, prior, nlpar = "") { +stan_mo <- function(bterms, data, ranef, prior) { # Stan code for monotonic effects - # Args: - # monef: names of the monotonic effects - # prior: a data.frame containing user defined priors - # as returned by check_prior - p <- usc(nlpar) out <- list() + monef <- all_terms(bterms$mo) + px <- check_prefix(bterms) + p <- usc(combine_prefix(px)) if (length(monef)) { I <- seq_along(monef) str_add(out$data) <- paste0( @@ -585,7 +566,7 @@ stan_mo <- function(monef, ranef, prior, nlpar = "") { " con_simplex", p, "_", I, "; \n" ) ) - bound <- get_bound(prior, class = "b", nlpar = nlpar) + bound <- get_bound(prior, class = "b", px = px) str_add(out$par) <- paste0( " // monotonic effects \n", " vector", bound, "[Kmo", p, "] bmo", p, "; \n", @@ -595,28 +576,27 @@ stan_mo <- function(monef, ranef, prior, nlpar = "") { ) ) str_add(out$prior) <- paste0( - stan_prior(prior, class = "b", coef = monef, - nlpar = nlpar, suffix = paste0("mo", p)), + stan_prior(prior, class = "b", coef = monef, + px = px, suffix = paste0("mo", p)), collapse( " target += dirichlet_lpdf(", "simplex", p, "_", I, " | con_simplex", p, "_", I, "); \n" ) ) - str_add(out$eta) <- stan_eta_mo(monef, ranef = ranef, nlpar = nlpar) + str_add(out$eta) <- stan_eta_mo( + monef, ranef = ranef, px = px + ) } out } -stan_cs <- function(csef, ranef, prior, nlpar = "") { +stan_cs <- function(bterms, data, ranef, prior) { # Stan code for category specific effects - # Args: - # csef: names of the category specific effects - # prior: a data.frame containing user defined priors - # as returned by check_prior - # (!) Not yet implemented for non-linear models - stopifnot(!nzchar(nlpar)) - ranef <- ranef[ranef$nlpar == nlpar & ranef$type == "cs", ] + # (!) Not implemented for non-linear models out <- list() + csef <- colnames(get_model_matrix(bterms$cs, data)) + px <- check_prefix(bterms) + ranef <- subset2(ranef, type = "cs", ls = px) if (length(csef)) { str_add(out$data) <- paste0( " int Kcs; // number of category specific effects \n", @@ -642,7 +622,7 @@ stan_cs <- function(csef, ranef, prior, nlpar = "") { str_add(out$modelD) <- paste0( " // linear predictor for category specific effects \n", " matrix[N, ncat - 1] mucs = rep_matrix(0, N, ncat - 1); \n" - ) + ) } cats <- get_matches("\\[[[:digit:]]+\\]$", ranef$coef) ncatM1 <- max(as.numeric(substr(cats, 2, nchar(cats) - 1))) @@ -653,7 +633,8 @@ stan_cs <- function(csef, ranef, prior, nlpar = "") { ) for (id in unique(r_cat$id)) { r <- r_cat[r_cat$id == id, ] - idp <- paste0(r$id, usc(r$nlpar, "prefix")) + rpx <- check_prefix(r) + idp <- paste0(r$id, usc(combine_prefix(rpx))) str_add(out$modelC2) <- collapse( " + r_", idp, "_", r$cn, "[J_", r$id, "[n]]", " * Z_", idp, "_", r$cn, "[n]" @@ -665,15 +646,15 @@ stan_cs <- function(csef, ranef, prior, nlpar = "") { out } -stan_me <- function(meef, ranef, prior, nlpar = "") { - # stan code for measurement error effects - # Args: - # meef: vector of terms containing noisy predictors +stan_me <- function(bterms, data, ranef, prior) { + # Stan code for measurement error effects out <- list() + meef <- get_me_labels(bterms, data = data) if (length(meef)) { not_one <- attr(meef, "not_one") uni_me <- attr(meef, "uni_me") - p <- usc(nlpar) + px <- check_prefix(bterms) + p <- usc(combine_prefix(px)) pK <- paste0(p, "_", seq_along(uni_me)) me_sp <- strsplit(gsub("[[:space:]]", "", meef), ":") @@ -697,14 +678,14 @@ stan_me <- function(meef, ranef, prior, nlpar = "") { # prepare linear predictor component meef <- rename(meef) meef_terms <- gsub(":", " .* ", meef_terms) - ranef <- ranef[ranef$nlpar == nlpar & ranef$type == "me", ] + ranef <- subset2(ranef, type = "me", ls = px) invalid_coef <- setdiff(ranef$coef, meef) if (length(invalid_coef)) { stop2("Noisy group-level terms require ", "corresponding population-level terms.") } for (i in seq_along(meef)) { - r <- ranef[ranef$coef == meef[i], ] + r <- subset2(ranef, coef = meef[i]) if (nrow(r)) { rpars <- paste0(" + ", stan_eta_r(r)) } else { @@ -737,7 +718,7 @@ stan_me <- function(meef, ranef, prior, nlpar = "") { ) str_add(out$prior) <- paste0( stan_prior(prior, class = "b", coef = meef, - nlpar = nlpar, suffix = paste0("me", p)), + px = px, suffix = paste0("me", p)), collapse( " target += normal_lpdf(Xme", pK, " | Xn", pK, ", noise", pK,"); \n" ) @@ -746,12 +727,12 @@ stan_me <- function(meef, ranef, prior, nlpar = "") { out } -stan_gp <- function(gpef, prior, nlpar = "") { +stan_gp <- function(bterms, data, prior) { # Stan code for latent gaussian processes - # Args: - # gpef: names of the gaussian process terms - p <- usc(nlpar) out <- list() + gpef <- get_gp_labels(bterms, data = data) + px <- check_prefix(bterms) + p <- usc(combine_prefix(px)) for (i in seq_along(gpef)) { pi <- paste0(p, "_", i) byvar <- attr(gpef, "byvars")[[i]] @@ -784,14 +765,14 @@ stan_gp <- function(gpef, prior, nlpar = "") { rgpef <- rename(gpef[i]) str_add(out$prior) <- paste0( stan_prior(prior, class = "sdgp", coef = rgpef, - nlpar = nlpar, suffix = pi), + px = px, suffix = pi), stan_prior(prior, class = "lscale", coef = rgpef, - nlpar = nlpar, suffix = pi), + px = px, suffix = pi), collapse(tp(), "normal_lpdf(zgp", pi, " | 0, 1); \n") ) if (byfac) { Jgp <- paste0("Jgp", pi, "_", J) - eta <- paste0(ifelse(nzchar(nlpar), nlpar, "mu"), "[", Jgp, "]") + eta <- paste0(combine_prefix(px, keep_mu = TRUE), "[", Jgp, "]") gp_args <- paste0( "Xgp", pi, "[", Jgp, "], sdgp", pi, "[", J, "], ", "lscale", pi, "[", J, "], zgp", pi, "[", Jgp, "]" @@ -810,15 +791,15 @@ stan_gp <- function(gpef, prior, nlpar = "") { out } -stan_eta_fe <- function(fixef, center_X = TRUE, - sparse = FALSE, nlpar = "") { +stan_eta_fe <- function(fixef, center_X = TRUE, sparse = FALSE, + px = list()) { # define Stan code to compute the fixef part of eta # Args: # fixef: names of the population-level effects # center_X: use the centered design matrix? # sparse: use sparse matrix multiplication? # nlpar: optional name of a non-linear parameter - p <- usc(nlpar) + p <- usc(combine_prefix(px)) if (length(fixef)) { if (sparse) { stopifnot(!center_X) @@ -836,16 +817,17 @@ stan_eta_fe <- function(fixef, center_X = TRUE, eta_fe } -stan_eta_re <- function(ranef, nlpar = "") { +stan_eta_re <- function(ranef, px = list()) { # write the group-level part of the linear predictor # Args: # ranef: a named list returned by tidy_ranef # nlpar: optional name of a non-linear parameter eta_re <- "" - ranef <- ranef[ranef$nlpar == nlpar & !nzchar(ranef$type), ] + ranef <- subset2(ranef, type = "", ls = px) for (id in unique(ranef$id)) { - r <- ranef[ranef$id == id, ] - idp <- paste0(r$id, usc(r$nlpar, "prefix")) + r <- subset2(ranef, id = id) + rpx <- check_prefix(r) + idp <- paste0(r$id, usc(combine_prefix(rpx))) str_add(eta_re) <- collapse( " + (", stan_eta_r(r), ") * Z_", idp, "_", r$cn, "[n]" ) @@ -860,7 +842,8 @@ stan_eta_r <- function(r) { # Returns: # A character vector, one element per row of 'r' stopifnot(nrow(r) > 0L, length(unique(r$gtype)) == 1L) - idp <- paste0(r$id, usc(r$nlpar, "prefix")) + rpx <- check_prefix(r) + idp <- paste0(r$id, usc(combine_prefix(rpx))) if (r$gtype[1] == "mm") { ng <- seq_along(r$gcall[[1]]$groups) out <- rep("", nrow(r)) @@ -876,22 +859,22 @@ stan_eta_r <- function(r) { out } -stan_eta_mo <- function(monef, ranef, nlpar = "") { +stan_eta_mo <- function(monef, ranef, px = list()) { # write the linear predictor for monotonic effects # Args: # monef: names of the monotonic effects # nlpar: an optional character string to add to the varnames # (used for non-linear models) - p <- usc(nlpar) + p <- usc(combine_prefix(px)) eta_mo <- "" - ranef <- ranef[ranef$nlpar == nlpar & ranef$type == "mo", ] + ranef <- subset2(ranef, type = "mo", ls = px) invalid_coef <- setdiff(ranef$coef, monef) if (length(invalid_coef)) { stop2("Monotonic group-level terms require ", "corresponding population-level terms.") } for (i in seq_along(monef)) { - r <- ranef[ranef$coef == monef[i], ] + r <- subset2(ranef, coef = monef[i]) if (nrow(r)) { rpars <- paste0(" + ", stan_eta_r(r)) } else { @@ -905,12 +888,12 @@ stan_eta_mo <- function(monef, ranef, nlpar = "") { eta_mo } -stan_eta_sm <- function(smooths, nlpar = "") { +stan_eta_sm <- function(smooths, px = list()) { # write the linear predictor for smooth terms # Args: # smooths: names of the smooth terms # nlpar: optional character string to add to the varnames - p <- usc(nlpar) + p <- usc(combine_prefix(px)) eta_smooths <- "" if (length(smooths)) { stopifnot(!is.null(attr(smooths, "nbases"))) @@ -925,14 +908,15 @@ stan_eta_sm <- function(smooths, nlpar = "") { eta_smooths } -stan_eta_ma <- function(autocor, nlpar = "") { +stan_eta_ma <- function(autocor, px = list()) { # write the linear predictor for MA terms # Args: # autocor: object of class cor_brms # nlpar: optional character string to add to the varnames out <- "" if (get_ma(autocor) && !use_cov(autocor)) { - str_add(out) <- paste0(" + head(E", nlpar, "[n], Kma) * ma") + px <- combine_prefix(px) + str_add(out) <- paste0(" + head(E", px, "[n], Kma) * ma") } out } @@ -972,12 +956,12 @@ stan_eta_transform <- function(family, llh_adj = FALSE) { (llh_adj || !stan_has_built_in_fun(family)) } -stan_eta_ilink <- function(family, auxpars = NULL, +stan_eta_ilink <- function(family, dpars = NULL, adforms = NULL, mix = "") { # correctly apply inverse link to eta # Args: # family: a list with elements 'family' and 'link - # auxpars: names of auxiliary parameters + # dpars: names of distributional parameters # adforms: list of formulas containing addition terms stopifnot(all(c("family", "link") %in% names(family))) llh_adj <- stan_llh_adj(adforms, c("cens", "trunc")) @@ -987,10 +971,10 @@ stan_eta_ilink <- function(family, auxpars = NULL, shape <- paste0("shape", mix) shape <- ifelse( is.formula(adforms$disp), paste0("disp_", shape, "[n]"), - ifelse(shape %in% auxpars, paste0(shape, "[n]"), shape) + ifelse(shape %in% dpars, paste0(shape, "[n]"), shape) ) nu <- paste0("nu", mix) - nu <- ifelse(nu %in% auxpars, paste0(nu, "[n]"), nu) + nu <- ifelse(nu %in% dpars, paste0(nu, "[n]"), nu) fl <- ifelse( family %in% c("gamma", "hurdle_gamma", "exponential"), paste0(family, "_", link), family @@ -1022,8 +1006,8 @@ stan_eta_ilink <- function(family, auxpars = NULL, out } -stan_auxpar_defs <- function(auxpar) { - # default Stan definitions for auxiliary parameters +stan_dpar_defs <- function(dpar) { + # default Stan definitions for distributional parameters default_defs <- list( sigma = c( " real ", @@ -1095,26 +1079,26 @@ stan_auxpar_defs <- function(auxpar) { ) # theta is handled in stan_mixture ) - def <- default_defs[[auxpar_class(auxpar)]] + def <- default_defs[[dpar_class(dpar)]] if (!is.null(def)) { - def <- paste0(def[1], auxpar, def[2]) + def <- paste0(def[1], dpar, def[2]) } else { def <- "" } def } -stan_auxpar_defs_temp <- function(auxpar) { - # default Stan definitions for temporary auxiliary parameters +stan_dpar_defs_temp <- function(dpar) { + # default Stan definitions for temporary distributional parameters default_defs <- list( xi = c( " real temp_", "; // unscaled shape parameter \n" ) ) - def <- default_defs[[auxpar_class(auxpar)]] + def <- default_defs[[dpar_class(dpar)]] if (!is.null(def)) { - def <- paste0(def[1], auxpar, def[2]) + def <- paste0(def[1], dpar, def[2]) } else { def <- "" } diff --git a/R/sysdata.rda b/R/sysdata.rda index 7da1593ed..be81b879e 100644 Binary files a/R/sysdata.rda and b/R/sysdata.rda differ diff --git a/R/validate.R b/R/validate.R index 50e806bfd..b15945857 100644 --- a/R/validate.R +++ b/R/validate.R @@ -59,6 +59,7 @@ parse_bf <- function(formula, family = NULL, autocor = NULL, advars <- str2formula(ulapply(adforms, all.vars)) y$adforms[names(adforms)] <- adforms + # copy stuff from the formula to parameter 'mu' str_rhs_form <- formula2str(rhs(formula)) terms <- try(terms(rhs(formula)), silent = TRUE) has_terms <- is(terms, "try-error") || @@ -70,12 +71,14 @@ parse_bf <- function(formula, family = NULL, autocor = NULL, mui <- paste0("mu", i) if (!is.formula(x$pforms[[mui]])) { x$pforms[[mui]] <- eval2(paste0(mui, str_rhs_form)) + attr(x$pforms[[mui]], "nl") <- attr(formula, "nl") rhs_needed <- TRUE } } } else { if (!is.formula(x$pforms[["mu"]])) { x$pforms[["mu"]] <- eval2(paste0("mu", str_rhs_form)) + attr(x$pforms[["mu"]], "nl") <- attr(formula, "nl") rhs_needed <- TRUE } x$pforms <- x$pforms[c("mu", setdiff(names(x$pforms), "mu"))] @@ -85,51 +88,60 @@ parse_bf <- function(formula, family = NULL, autocor = NULL, "the right-hand side of 'formula' is unused.") } - auxpars <- is_auxpar_name(names(x$pforms), family, bterms = y) - auxpars <- names(x$pforms)[auxpars] - # amend when generalizing non-linear models to auxiliary parameters - nlpars <- setdiff(names(x$pforms), auxpars) - if (isTRUE(x[["nl"]])) { - if (is.mixfamily(family) || is_ordinal(family) || is_categorical(family)) { - stop2("Non-linear formulas are not yet allowed for this family.") - } - y$auxpars[["mu"]] <- parse_nlf(x$pforms[["mu"]], x$pforms[nlpars]) - y$auxpars[["mu"]]$family <- auxpar_family(family, "mu") - auxpars <- setdiff(auxpars, "mu") - } else { - if (length(nlpars)) { - stop2("Parameter '", nlpars[1], "' is not part of the model.") + # predicted distributional parameters + dpars <- is_dpar_name(names(x$pforms), family, bterms = y) + dpars <- names(x$pforms)[dpars] + dpar_forms <- x$pforms[dpars] + nlpars <- setdiff(names(x$pforms), dpars) + nlpar_forms <- x$pforms[nlpars] + dpar_of_nlpars <- rep(NA, length(nlpars)) + for (i in seq_along(nlpar_forms)) { + dpar <- attr(nlpar_forms[[i]], "dpar", TRUE) + if (length(dpar) != 1L) { + stop2("Parameter '", nlpars[i], "' is not part of the model") } + dpar_of_nlpars[i] <- dpar } - - # predicted auxiliary parameters - for (ap in auxpars) { - y$auxpars[[ap]] <- parse_lf(x$pforms[[ap]], family = family) - y$auxpars[[ap]]$family <- auxpar_family(family, ap) + for (dp in dpars) { + if (isTRUE(attr(dpar_forms[[dp]], "nl"))) { + if (is.mixfamily(family) || is_ordinal(family) || is_categorical(family)) { + stop2("Non-linear formulas are not yet allowed for this family.") + } + dpar_nlpar_forms <- nlpar_forms[dpar_of_nlpars %in% dp] + y$dpars[[dp]] <- parse_nlf(dpar_forms[[dp]], dpar_nlpar_forms, dp) + } else { + y$dpars[[dp]] <- parse_lf(dpar_forms[[dp]], family = family) + } + y$dpars[[dp]]$family <- dpar_family(family, dp) + y$dpars[[dp]]$dpar <- dp + # temporary until brms 2.0.0 + y$dpars[[dp]]$resp <- "" } - if (!is.null(y$auxpars[["mu"]])) { - y$auxpars[["mu"]][["autocor"]] <- autocor + # fixed distributional parameters + inv_fixed_dpars <- setdiff(names(x$pfix), valid_dpars(family, y)) + if (length(inv_fixed_dpars)) { + stop2("Invalid fixed parameters: ", collapse_comma(inv_fixed_dpars)) } - # fixed auxiliary parameters - inv_fauxpars <- setdiff(names(x$pfix), valid_auxpars(family, y)) - if (length(inv_fauxpars)) { - stop2("Invalid auxiliary parameters: ", collapse_comma(inv_fauxpars)) + for (dp in names(x$pfix)) { + y$fdpars[[dp]] <- list(value = x$pfix[[dp]], dpar = dp) } - y$fauxpars <- x$pfix - check_fauxpars(y$fauxpars) + check_fdpars(y$fdpars) # check for illegal use of cs terms if (has_cs(y) && !(is.null(family) || allows_cs(family))) { stop2("Category specific effects are only meaningful for ", "families 'sratio', 'cratio', and 'acat'.") } # parse autocor formula + if (!is.null(y$dpars[["mu"]])) { + y$dpars[["mu"]][["autocor"]] <- autocor + } y$time <- parse_time(autocor) # make a formula containing all required variables lhsvars <- if (resp_rhs_all) all.vars(y$respform) lformula <- c( lhsvars, advars, - lapply(y$auxpars, "[[", "allvars"), + lapply(y$dpars, "[[", "allvars"), y$time$allvars ) y$allvars <- allvars_formula(lformula) @@ -143,8 +155,8 @@ parse_bf <- function(formula, family = NULL, autocor = NULL, stop2("Multivariate models currently allow only ", "addition argument 'weights'.") } - if (length(y$auxpars) > 1L) { - stop2("Auxiliary parameters cannot yet be ", + if (length(y$dpars) > 1L) { + stop2("Distributional parameters cannot yet be ", "predicted in multivariate models.") } } @@ -226,7 +238,7 @@ parse_lf <- function(formula, family = NULL) { structure(y, class = "btl") } -parse_nlf <- function(formula, nlpar_forms) { +parse_nlf <- function(formula, nlpar_forms, dpar) { # parse non-linear formulas # Args: # formula: formula of the non-linear model @@ -235,6 +247,7 @@ parse_nlf <- function(formula, nlpar_forms) { # Returns: # object of class 'btnl' stopifnot(is.list(nlpar_forms)) + stopifnot(length(dpar) == 1L) formula <- rhs(as.formula(formula)) y <- list() if (length(nlpar_forms)) { @@ -243,6 +256,10 @@ parse_nlf <- function(formula, nlpar_forms) { y$nlpars <- named_list(nlpars) for (nlp in nlpars) { y$nlpars[[nlp]] <- parse_lf(nlpar_forms[[nlp]]) + y$nlpars[[nlp]]$nlpar <- nlp + y$nlpars[[nlp]]$dpar <- dpar + # temporary until brms 2.0.0 + y$nlpars[[nlp]]$resp <- "" } model_vars <- all.vars(formula) missing_pars <- setdiff(names(y$nlpars), model_vars) @@ -567,46 +584,78 @@ is.btnl <- function(x) { inherits(x, "btnl") } -avoid_auxpars <- function(names, bterms) { +avoid_dpars <- function(names, bterms) { # avoid ambiguous parameter names # Args: # names: names to check for ambiguity # bterms: object of class brmsterms #stopifnot(is.brmsterms(bterms)) - auxpars <- c(names(bterms$auxpars), "mo", "cs", "me") - if (length(auxpars)) { - auxpars_prefix <- paste0("^", auxpars, "_") - invalid <- any(ulapply(auxpars_prefix, grepl, names)) + dpars <- c(names(bterms$dpars), "mo", "cs", "me") + if (length(dpars)) { + dpars_prefix <- paste0("^", dpars, "_") + invalid <- any(ulapply(dpars_prefix, grepl, names)) if (invalid) { - auxpars <- paste0("'", auxpars, "_'", collapse = ", ") - stop2("Variable names starting with ", auxpars, + dpars <- paste0("'", dpars, "_'", collapse = ", ") + stop2("Variable names starting with ", dpars, " are not allowed for this model.") } } invisible(NULL) } -check_nlpar <- function(nlpar) { - stopifnot(length(nlpar) == 1L) - nlpar <- as.character(nlpar) - ifelse(nlpar == "mu", "", nlpar) +vars_prefix <- function() { + c("dpar", "resp", "nlpar") } -check_fauxpars <- function(x) { - # check validity of fixed auxiliary parameters +check_prefix <- function(x, keep_mu = FALSE) { + vpx <- vars_prefix() + if (is.data.frame(x) && nrow(x) == 0) { + # avoids a bug in data.frames with zero rows + x <- list() + } + x[setdiff(vpx, names(x))] <- "" + x <- x[vpx] + for (i in seq_along(x)) { + x[[i]] <- as.character(x[[i]]) + if (!length(x[[i]])) { + x[[i]] <- "" + } + x[[i]] <- ifelse( + !keep_mu & names(x)[i] == "dpar" & x[[i]] %in% "mu", + yes = "", no = x[[i]] + ) + x[[i]] <- ifelse( + keep_mu & names(x)[i] == "dpar" & x[[i]] %in% "", + yes = "mu", no = x[[i]] + ) + } + x +} + +combine_prefix <- function(prefix, keep_mu = FALSE) { + prefix <- check_prefix(prefix, keep_mu = keep_mu) + prefix <- lapply(prefix, function(x) + ifelse(nzchar(x), paste0("_", x), "") + ) + sub("^_", "", do.call(paste0, prefix)) +} + +check_fdpars <- function(x) { + # check validity of fixed distributional parameters stopifnot(is.null(x) || is.list(x)) pos_pars <- c( "sigma", "shape", "nu", "phi", "kappa", "beta", "disc", "bs", "ndt", "theta" ) prob_pars <- c("zi", "hu", "bias", "quantile") - for (ap in names(x)) { - apc <- auxpar_class(ap) - if (apc %in% pos_pars && x[[ap]] < 0) { - stop2("Parameter '", ap, "' must be positive.") + for (dp in names(x)) { + apc <- dpar_class(dp) + value <- x[[dp]]$value + if (apc %in% pos_pars && value < 0) { + stop2("Parameter '", dp, "' must be positive.") } - if (apc %in% prob_pars && (x[[ap]] < 0 || x[[ap]] > 1)) { - stop2("Parameter '", ap, "' must be between 0 and 1.") + if (apc %in% prob_pars && (value < 0 || value > 1)) { + stop2("Parameter '", dp, "' must be between 0 and 1.") } } invisible(TRUE) @@ -714,8 +763,8 @@ check_re_formula <- function(re_formula, formula) { re_formula <- ~ 1 } else { re_formula <- get_re_terms(as.formula(re_formula), formula = TRUE) - new <- parse_bf(re_formula, check_response = FALSE)$auxpars$mu$re - old <- parse_bf(old_re_formula, check_response = FALSE)$auxpars$mu$re + new <- parse_bf(re_formula, check_response = FALSE)$dpars$mu$re + old <- parse_bf(old_re_formula, check_response = FALSE)$dpars$mu$re if (nrow(new)) { new_terms <- lapply(new$form, terms) found <- rep(FALSE, nrow(new)) @@ -759,7 +808,10 @@ update_re_terms <- function(x, re_formula = NULL) { # add valid group-level terms of re_formula # Args: # formula: object of class 'formula' - formula <- as.formula(formula) + if (isTRUE(attr(formula, "nl"))) { + # non-linear formulas contain no group-level effects + return(formula) + } re_formula <- check_re_formula(re_formula, formula) new_formula <- formula2str(formula) old_re_terms <- get_re_terms(formula) @@ -775,16 +827,14 @@ update_re_terms <- function(x, re_formula = NULL) { new_re_terms <- get_re_terms(re_formula) new_formula <- paste(c(new_formula, new_re_terms), collapse = "+") new_formula <- formula(new_formula) - environment(new_formula) <- environment(formula) + attributes(new_formula) <- attributes(formula) return(new_formula) } if (is.formula(x)) { x <- .update_re_terms(x, re_formula) } else if (is.brmsformula(x)) { - if (!x[["nl"]]) { - x$formula <- .update_re_terms(x$formula, re_formula) - } + x$formula <- .update_re_terms(x$formula, re_formula) x$pforms <- lapply(pforms(x), .update_re_terms, re_formula) } else { stop("Don't know how to handle objects of class ", @@ -854,37 +904,37 @@ get_re.brmsterms <- function(x, all = TRUE, ...) { # all: logical; include ranefs of nl and aux parameters? old_mv <- isTRUE(attr(x$formula, "old_mv")) if (all) { - re <- named_list(names(x$auxpars)) - for (ap in names(re)) { - re[[ap]] <- get_re( - x$auxpars[[ap]], nlpar = ap, - response = x$response, old_mv = old_mv + re <- named_list(names(x$dpars)) + for (dp in names(re)) { + re[[dp]] <- get_re( + x$dpars[[dp]], response = x$response, old_mv = old_mv ) } re <- do.call(rbind, re) } else { - x$auxpars[["mu"]]$nlpars <- NULL - re <- get_re(x$auxpars[["mu"]]) + x$dpars[["mu"]]$nlpars <- NULL + re <- get_re(x$dpars[["mu"]]) } re } #' @export -get_re.btl <- function(x, nlpar = "", response = "", old_mv = FALSE, ...) { +get_re.btl <- function(x, response = "", old_mv = FALSE, ...) { stopifnot(is.data.frame(x$re)) - nlpar <- check_nlpar(nlpar) + px <- check_prefix(x) re <- x$re nresp <- length(response) if (!old_mv && nresp > 1L && nrow(re)) { - # new MV models are also using the 'nlpar' argument re <- replicate(nresp, re, simplify = FALSE) for (i in seq_len(nresp)) { - re[[i]]$nlpar <- rep(response[i], nrow(re[[i]])) + re[[i]]$resp <- rep(response[i], nrow(re[[i]])) } re <- do.call(rbind, re) } else { - re$nlpar <- rep(nlpar, nrow(re)) + re$resp <- rep(px$resp, nrow(re)) } + re$dpar <- rep(px$dpar, nrow(re)) + re$nlpar <- rep(px$nlpar, nrow(re)) re } @@ -892,7 +942,7 @@ get_re.btl <- function(x, nlpar = "", response = "", old_mv = FALSE, ...) { get_re.btnl <- function(x, ...) { re <- named_list(names(x$nlpars)) for (nlp in names(re)) { - re[[nlp]] <- get_re(x$nlpars[[nlp]], nlpar = nlp) + re[[nlp]] <- get_re(x$nlpars[[nlp]]) } do.call(rbind, re) } @@ -902,15 +952,15 @@ get_effect.brmsterms <- function(x, target = "fe", all = TRUE, ...) { # get formulas of certain effects in a list # Args: # target: type of effects to return - # all: logical; include effects of nlpars and auxpars? + # all: logical; include effects of nlpars and dpars? if (all) { - out <- named_list(names(x$auxpars)) - for (ap in names(out)) { - out[[ap]] <- get_effect(x$auxpars[[ap]], target = target) + out <- named_list(names(x$dpars)) + for (dp in names(out)) { + out[[dp]] <- get_effect(x$dpars[[dp]], target = target) } } else { - x$auxpars[["mu"]]$nlpars <- NULL - out <- get_effect(x$auxpars[["mu"]], target = target) + x$dpars[["mu"]]$nlpars <- NULL + out <- get_effect(x$dpars[["mu"]], target = target) } unlist(out, recursive = FALSE) } @@ -958,8 +1008,8 @@ get_all_effects.brmsterms <- function(x, rsv_vars = NULL, #stopifnot(is.brmsterms(bterms)) stopifnot(is.atomic(rsv_vars)) out <- list() - for (ap in names(x$auxpars)) { - out <- c(out, get_all_effects(x$auxpars[[ap]])) + for (dp in names(x$dpars)) { + out <- c(out, get_all_effects(x$dpars[[dp]])) } out <- rmNULL(lapply(out, setdiff, y = rsv_vars)) if (comb_all) { @@ -1011,7 +1061,7 @@ get_sm_labels <- function(x, data = NULL, covars = FALSE, # or just return the covariate names (FALSE)? if (is.formula(x)) { x <- parse_bf(x, check_response = FALSE) - sm_form <- x$auxpars$mu[["sm"]] + sm_form <- x$dpars$mu[["sm"]] } else { sm_form <- x[["sm"]] } @@ -1040,7 +1090,7 @@ get_sm_labels <- function(x, data = NULL, covars = FALSE, stop("Invalid combination of arguments. Please report a bug.") } sdata_fe <- data_fe(x, data, knots = attr(data, "knots")) - by_levels <- attr(sdata_fe[["X"]], "by_levels") + by_levels <- attr(sdata_fe$X, "by_levels") nby <- lengths(by_levels) sms <- as.list(sms) for (i in seq_along(sms)) { @@ -1066,7 +1116,7 @@ get_me_labels <- function(x, data) { # data: data frame containing the noisy variables if (is.formula(x)) { x <- parse_bf(x, check_response = FALSE) - me_form <- x$auxpars$mu[["me"]] + me_form <- x$dpars$mu[["me"]] } else { me_form <- x[["me"]] } @@ -1088,7 +1138,7 @@ get_gp_labels <- function(x, data = NULL, covars = FALSE) { # returned instead of the full term names? if (is.formula(x)) { x <- parse_bf(x, check_response = FALSE) - gp_form <- x$auxpars$mu[["gp"]] + gp_form <- x$dpars$mu[["gp"]] } else { gp_form <- x[["gp"]] } @@ -1241,7 +1291,7 @@ tidy_ranef <- function(bterms, data = NULL, all = TRUE, # Args: # bterms: object of class brmsterms # data: data passed to brm after updating - # all: include REs of non-linear and auxiliary parameters? + # all: include REs of non-linear and distributional parameters? # ncat: optional number of response categories # only used for category specific group-level effects # Returns: @@ -1263,7 +1313,7 @@ tidy_ranef <- function(bterms, data = NULL, all = TRUE, j <- 1 for (i in seq_len(nrow(re))) { if (re$type[[i]] == "mo") { - coef <- prepare_mo_vars(re$form[[i]], data, check = FALSE) + coef <- mo_design_matrix(re$form[[i]], data, check = FALSE) coef <- colnames(coef) } else if (re$type[[i]] == "cs") { coef <- colnames(get_model_matrix(re$form[[i]], data = data)) @@ -1279,13 +1329,15 @@ tidy_ranef <- function(bterms, data = NULL, all = TRUE, } else { coef <- colnames(get_model_matrix(re$form[[i]], data = data)) } - avoid_auxpars(coef, bterms = bterms) + avoid_dpars(coef, bterms = bterms) rdat <- data.frame( id = re$id[[i]], group = re$group[[i]], gn = re$gn[[i]], gtype = re$gtype[[i]], coef = coef, cn = NA, + resp = re$resp[[i]], + dpar = re$dpar[[i]], nlpar = re$nlpar[[i]], cor = re$cor[[i]], type = re$type[[i]], @@ -1326,7 +1378,7 @@ tidy_ranef <- function(bterms, data = NULL, all = TRUE, stop2("Grouping factor names ", inv_groups, " are resevered.") } # check for duplicated and thus not identified effects - dup <- duplicated(ranef[, c("group", "coef", "nlpar")]) + dup <- duplicated(ranef[, c("group", "coef", vars_prefix())]) if (any(dup)) { stop2("Duplicated group-level effects are not allowed.") } @@ -1356,8 +1408,9 @@ tidy_ranef <- function(bterms, data = NULL, all = TRUE, empty_ranef <- function() { data.frame( id = numeric(0), group = character(0), gn = numeric(0), - coef = character(0), cn = numeric(0), nlpar = character(0), - cor = logical(0), type = character(0), form = character(0), + coef = character(0), cn = numeric(0), resp = character(0), + dpar = character(0), nlpar = character(0), cor = logical(0), + type = character(0), form = character(0), stringsAsFactors = FALSE ) } @@ -1371,7 +1424,7 @@ rsv_vars <- function(bterms, incl_intercept = TRUE) { nresp <- length(bterms$response) family <- bterms$family old_mv <- isTRUE(attr(bterms$formula, "old_mv")) - rsv_intercept <- any(ulapply(bterms$auxpars, has_rsv_intercept)) + rsv_intercept <- any(ulapply(bterms$dpars, has_rsv_intercept)) if (old_mv) { if (is_linear(family) && nresp > 1L || is_categorical(family)) { rsv <- c("trait", "response") @@ -1428,23 +1481,23 @@ exclude_pars <- function(bterms, data = NULL, ranef = empty_ranef(), # save_mevars: should samples of noise-free variables be saved? # Returns: # a vector of parameters to be excluded - .exclude_pars <- function(bt, nlpar = "") { + .exclude_pars <- function(bt) { stopifnot(is.btl(bt)) - nlpar <- usc(check_nlpar(nlpar)) + p <- usc(combine_prefix(bt)) out <- c( - paste0("temp", nlpar, "_Intercept"), - paste0(c("hs_local", "hs_global", "zb"), nlpar) + paste0("temp", p, "_Intercept"), + paste0(c("hs_local", "hs_global", "zb"), p) ) sms <- get_sm_labels(bt, data) if (length(sms) && !is.null(data)) { for (i in seq_along(sms)) { nb <- seq_len(attr(sms, "nbases")[[i]]) - out <- c(out, paste0("zs", nlpar, "_", i, "_", nb)) + out <- c(out, paste0("zs", p, "_", i, "_", nb)) } } meef <- get_me_labels(bt, data) if (!save_mevars && length(meef)) { - out <- c(out, paste0("Xme", nlpar, "_", seq_along(meef))) + out <- c(out, paste0("Xme", p, "_", seq_along(meef))) } return(out) } @@ -1455,7 +1508,7 @@ exclude_pars <- function(bterms, data = NULL, ranef = empty_ranef(), save_all_pars <- as_one_logical(save_all_pars) out <- c( "Rescor", "Sigma", "res_cov_matrix", - intersect(auxpars(), names(bterms$auxpars)) + intersect(dpars(), names(bterms$dpars)) ) if (!save_all_pars) { out <- c(out, @@ -1464,18 +1517,19 @@ exclude_pars <- function(bterms, data = NULL, ranef = empty_ranef(), ) if (length(bterms$response) > 1L) { for (r in bterms$response) { - out <- c(out, .exclude_pars(bterms$auxpars$mu, nlpar = r)) + bterms$dpars$mu$resp <- r + out <- c(out, .exclude_pars(bterms$dpars$mu)) } - bterms$auxpars$mu <- NULL + bterms$dpars$mu <- NULL } - for (ap in names(bterms$auxpars)) { - bt <- bterms$auxpars[[ap]] + for (dp in names(bterms$dpars)) { + bt <- bterms$dpars[[dp]] if (length(bt$nlpars)) { for (nlp in names(bt$nlpars)) { - out <- c(out, .exclude_pars(bt$nlpars[[nlp]], nlpar = nlp)) + out <- c(out, .exclude_pars(bt$nlpars[[nlp]])) } } else { - out <- c(out, .exclude_pars(bt, nlpar = ap)) + out <- c(out, .exclude_pars(bt)) } } } @@ -1486,8 +1540,8 @@ exclude_pars <- function(bterms, data = NULL, ranef = empty_ranef(), out <- c(out, paste0(rm_re_pars, "_", id)) } if (!save_ranef) { - usc_nlpar <- usc(ranef$nlpar, "prefix") - out <- c(out, paste0("r_", ranef$id, usc_nlpar, "_", ranef$cn)) + p <- usc(combine_prefix(ranef)) + out <- c(out, paste0("r_", ranef$id, p, "_", ranef$cn)) } } att <- nlist(save_ranef, save_mevars, save_all_pars) diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd index 9c7eaa771..35a4fdc49 100644 --- a/inst/NEWS.Rd +++ b/inst/NEWS.Rd @@ -10,6 +10,26 @@ \code{bridge_sampler}, \code{bayes_factor}, and \code{post_prob} all powered by the \pkg{bridgesampling} package. + \item Specify non-linear models for all + distributional parameters. + \item Combine multiple model formulas using + the \code{+} operator and the helper functions + \code{lf}, \code{nlf}, and \code{set_nl}. + \item Combine multiple priors using the + \code{+} operator. + \item Split the \code{nlpar} argument of + \code{set_prior} into the three arguments + \code{resp}, \code{dpar}, and \code{nlpar} + to allow for more flexible prior specifications. + } + } + \subsection{OTHER CHANGES}{ + \itemize{ + \item Refactor parts of the package to prepare + for the implementation of more flexible + multivariate models in future updates. + \item Rename argument \code{auxpar} of + \code{fitted.brmsfit} to \code{dpar}. } } } diff --git a/man/brmsformula-helpers.Rd b/man/brmsformula-helpers.Rd new file mode 100644 index 000000000..f9ec0a5cd --- /dev/null +++ b/man/brmsformula-helpers.Rd @@ -0,0 +1,93 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/brmsformula.R +\name{brmsformula-helpers} +\alias{brmsformula-helpers} +\alias{nlf} +\alias{bf-helpers} +\alias{lf} +\alias{set_nl} +\alias{lf} +\alias{set_nl} +\title{Linear and Non-linear formulas in \pkg{brms}} +\usage{ +nlf(formula, ..., flist = NULL, dpar = NULL) + +lf(..., flist = NULL, dpar = NULL) + +set_nl(nl = TRUE, dpar = NULL) +} +\arguments{ +\item{formula}{Non-linear formula for a distributional parameter. +The name of the distributional parameter can either be specified +on the left-hand side of \code{formula} or via argument \code{dpar}.} + +\item{...}{Additional \code{formula} objects to specify +predictors of non-linear and distributional parameters. +Formulas can either be named directly or contain +names on their left-hand side. +The following are distributional parameters of specific families +(all other parameters are treated as non-linear parameters): +\code{sigma} (residual standard deviation or scale of +the \code{gaussian}, \code{student}, \code{lognormal} +\code{exgaussian}, and \code{asym_laplace} families); +\code{shape} (shape parameter of the \code{Gamma}, +\code{weibull}, \code{negbinomial}, and related +zero-inflated / hurdle families); \code{nu} +(degrees of freedom parameter of the \code{student} family); +\code{phi} (precision parameter of the \code{beta} +and \code{zero_inflated_beta} families); +\code{kappa} (precision parameter of the \code{von_mises} family); +\code{beta} (mean parameter of the exponential componenent +of the \code{exgaussian} family); +\code{quantile} (quantile parameter of the \code{asym_laplace} family); +\code{zi} (zero-inflation probability); +\code{hu} (hurdle probability); +\code{zoi} (zero-one-inflation probability); +\code{coi} (conditional one-inflation probability); +\code{disc} (discrimination) for ordinal models; +\code{bs}, \code{ndt}, and \code{bias} (boundary separation, +non-decision time, and initial bias of the \code{wiener} +diffusion model). +All distributional parameters are modeled +on the log or logit scale to ensure correct definition +intervals after transformation. +See 'Details' for more explanation.} + +\item{flist}{Optional list of formulas, which are treated in the +same way as formulas passed via the \code{...} argument.} + +\item{dpar}{Optional character string specifying the distributional +parameter to which the formulas passed via \code{...} and +\code{flist} belong.} + +\item{nl}{Logical; Indicates whether \code{formula} should be +treated as specifying a non-linear model. By default, \code{formula} +is treated as an ordinary linear model formula.} +} +\value{ +For \code{lf} and \code{nlf} a \code{list} that can be + passed to \code{\link[brms:brmsformula]{brmsformula}} or added + to an existing \code{brmsformula} object. For \code{set_nl} + a \code{list} that can be added to an existing + \code{brmsformula} object. +} +\description{ +Helper functions to specify linear and non-linear +formulas for use with \code{\link[brms:brmsformula]{brmsformula}}. +} +\examples{ +# add more formulas to the model +bf(y ~ 1) + + nlf(sigma ~ a * exp(b * x), a ~ x) + + lf(b ~ z + (1|g), dpar = "sigma") + + gaussian() + +# specify 'nl' later on +bf(y ~ a * inv_logit(x * b)) + + lf(a + b ~ z) + + set_nl(TRUE) + +} +\seealso{ +\code{\link[brms:brmsformula]{brmsformula}} +} diff --git a/man/brmsformula.Rd b/man/brmsformula.Rd index ce87612cf..7e0fc997b 100644 --- a/man/brmsformula.Rd +++ b/man/brmsformula.Rd @@ -15,10 +15,10 @@ a symbolic description of the model to be fitted. The details of model specification are given in 'Details'.} \item{...}{Additional \code{formula} objects to specify -predictors of non-linear and auxiliary parameters. +predictors of non-linear and distributional parameters. Formulas can either be named directly or contain names on their left-hand side. -The following are auxiliary parameters of specific families +The following are distributional parameters of specific families (all other parameters are treated as non-linear parameters): \code{sigma} (residual standard deviation or scale of the \code{gaussian}, \code{student}, \code{lognormal} @@ -41,7 +41,7 @@ of the \code{exgaussian} family); \code{bs}, \code{ndt}, and \code{bias} (boundary separation, non-decision time, and initial bias of the \code{wiener} diffusion model). -All auxiliary parameters are modeled +All distributional parameters are modeled on the log or logit scale to ensure correct definition intervals after transformation. See 'Details' for more explanation.} @@ -323,7 +323,7 @@ models for all parameters of the assumed response distribution. As of \pkg{brms} 1.0.0, zero-inflated and hurdle models are specfied in the same way as as their non-inflated counterparts. - However, they have additional auxiliary parameters + However, they have additional distributional parameters (named \code{zi} and \code{hu} respectively) modeling the zero-inflation / hurdle probability depending on which model you choose. These parameters can also be affected by predictors @@ -389,15 +389,15 @@ models for all parameters of the assumed response distribution. For this reason it is mandatory to specify priors on the non-linear parameters. For instructions on how to do that, see \code{\link[brms:set_prior]{set_prior}}. - \bold{Formula syntax for predicting auxiliary parameters} + \bold{Formula syntax for predicting distributional parameters} - It is also possible to predict auxiliary parameters of the response + It is also possible to predict parameters of the response distribution such as the residual standard deviation \code{sigma} in gaussian models or the hurdle probability \code{hu} in hurdle models. The syntax closely resembles that of a non-linear parameter, for instance \code{sigma ~ x + s(z) + (1+x|g)}. - Alternatively, one may fix auxiliary parameters to certain values. + Alternatively, one may fix distributional parameters to certain values. However, this is mainly useful when models become too complicated and otherwise have convergence issues. We thus suggest to be generally careful when making use of this option. @@ -414,19 +414,19 @@ models for all parameters of the assumed response distribution. the negative binomial distribution with \code{shape = 1}. Furthermore, the parameter \code{disc} ('discrimination') in ordinal models is fixed to \code{1} by default and not estimated, - but may be modeled as any other auxiliary parameter if desired + but may be modeled as any other distributional parameter if desired (see examples). For reasons of identification, \code{'disc'} can only be positive, which is achieved by applying the log-link. - All auxiliary parameters currently supported by \code{brmsformula} + All distributional parameters currently supported by \code{brmsformula} have to positive (a negative standard deviation or precision parameter doesn't make any sense) or are bounded between 0 and 1 (for zero-inflated / hurdle proabilities, quantiles, or the intial bias parameter of drift-diffusion models). However, linear predictors can be positive or negative, and thus the log link (for positive parameters) or logit link (for probability parameters) are used - by default to ensure that auxiliary parameters are within their valid intervals. - This implies that, by default, effects for auxiliary parameters are estimated + by default to ensure that distributional parameters are within their valid intervals. + This implies that, by default, effects for distributional parameters are estimated on the log / logit scale and one has to apply the inverse link function to get to the effects on the original scale. Alternatively, it is possible to use the identity link to predict parameters @@ -445,17 +445,17 @@ models for all parameters of the assumed response distribution. terms allowed in non-mixture models are allowed in mixture models as well. - Auxiliary parameters of mixture distributions have the same + distributional parameters of mixture distributions have the same name as those of the corresponding ordinary distributions, but with a number at the end to indicate the mixture component. For instance, if - you use family \code{mixture(gaussian, gaussian)}, the auxiliary + you use family \code{mixture(gaussian, gaussian)}, the distributional parameters are \code{sigma1} and \code{sigma2}. - Auxiliary parameters of the same class can be fixed to the same value. + distributional parameters of the same class can be fixed to the same value. For the above example, we could write \code{sigma2 = "sigma1"} to make sure that both components have the same residual standard deviation, which is in turn estimated from the data. - In addition, there are two types of special auxiliary parameters. + In addition, there are two types of special distributional parameters. The first are named \code{mu}, that allow for modeling different predictors for the mean parameters of different mixture components. For instance, if you want to predict the mean of the first component @@ -534,4 +534,13 @@ bf(y ~ 1, mu1 ~ x, mu2 ~ z, family = mix) # fix both residual standard deviations to the same value bf(y ~ x, sigma2 = "sigma1", family = mix) +# use the '+' operator to specify models +bf(y ~ 1) + + nlf(sigma ~ a * exp(b * x), a ~ x) + + lf(b ~ z + (1|g), dpar = "sigma") + + gaussian() + +} +\seealso{ +\code{\link[brms:brmsformula-helpers]{brmsformula-helpers}} } diff --git a/man/fitted.brmsfit.Rd b/man/fitted.brmsfit.Rd index 5dbf10f36..ffe446cbd 100644 --- a/man/fitted.brmsfit.Rd +++ b/man/fitted.brmsfit.Rd @@ -7,7 +7,7 @@ \method{fitted}{brmsfit}(object, newdata = NULL, re_formula = NULL, scale = c("response", "linear"), allow_new_levels = FALSE, sample_new_levels = "uncertainty", new_objects = list(), - incl_autocor = TRUE, auxpar = NULL, subset = NULL, nsamples = NULL, + incl_autocor = TRUE, dpar = NULL, subset = NULL, nsamples = NULL, sort = FALSE, nug = NULL, summary = TRUE, robust = FALSE, probs = c(0.025, 0.975), ...) } @@ -57,7 +57,7 @@ parameters should be included in the predictions. Defaults to correlation structures such as \code{\link[brms:cor_bsts]{cor_bsts}}, or \code{\link[brms:cor_fixed]{cor_fixed}}.} -\item{auxpar}{Optional name of a predicted auxiliary parameter. +\item{dpar}{Optional name of a predicted distributional parameter. If specified, fitted values of this parameters are returned.} \item{subset}{A numeric vector specifying diff --git a/man/horseshoe.Rd b/man/horseshoe.Rd index 97190828e..5a93df652 100644 --- a/man/horseshoe.Rd +++ b/man/horseshoe.Rd @@ -34,7 +34,7 @@ of observations in the data.} \item{autoscale}{Logical; indicating whether the horseshoe prior should be scaled using the residual standard deviation \code{sigma} if possible and sensible (defaults to \code{TRUE}). -Autoscaling is not applied for auxiliary parameters or +Autoscaling is not applied for distributional parameters or when the model does not contain the parameter \code{sigma}.} } \value{ diff --git a/man/set_prior.Rd b/man/set_prior.Rd index f5e3e8119..35c9cf0ad 100644 --- a/man/set_prior.Rd +++ b/man/set_prior.Rd @@ -9,8 +9,8 @@ \alias{prior_string} \title{Prior Definitions for \pkg{brms} Models} \usage{ -set_prior(prior, class = "b", coef = "", group = "", nlpar = "", - resp = NULL, lb = NULL, ub = NULL, check = TRUE) +set_prior(prior, class = "b", coef = "", group = "", resp = "", + dpar = "", nlpar = "", lb = NULL, ub = NULL, check = TRUE) prior(prior, ...) @@ -29,12 +29,14 @@ See 'Details' for other valid parameter classes.} \item{group}{Grouping factor of group-level parameters.} -\item{nlpar}{Name of a non-linear / auxiliary parameter. -Only used in non-linear / distributional models.} - \item{resp}{Name of the response variable / category. -Only used in multivariate and categorical models. -Is internally handled as an alias of \code{nlpar}.} +Only used in multivariate and categorical models.} + +\item{dpar}{Name of a distributional parameter. +Only used in distributional models.} + +\item{nlpar}{Name of a non-linear parameter. +Only used in non-linear models.} \item{lb}{Lower bound for parameter restriction. Currently only allowed for classes \code{"b"}, \code{"ar"}, \code{"ma"}, and \code{"arr"}. diff --git a/tests/brmsfit_examples.R b/tests/brmsfit_examples.R index 0f2ee1fe1..17238048e 100644 --- a/tests/brmsfit_examples.R +++ b/tests/brmsfit_examples.R @@ -24,7 +24,7 @@ brmsfit_example1 <- brm( autocor = cor_arma(~visit|patient, 1, 1), prior = c(set_prior("normal(0,5)", class = "b"), set_prior("cauchy(0,2)", class = "sd"), - set_prior("normal(0,3)", nlpar = "sigma")), + set_prior("normal(0,3)", dpar = "sigma")), sample_prior = TRUE, warmup = 150, iter = 200, chains = 2, save_dso = FALSE, testmode = TRUE ) @@ -56,9 +56,9 @@ brmsfit_example4 <- brm( brmsfit_example5 <- brm( bf(count ~ Age + (1|visit), mu2 ~ Age), dat, family = mixture(gaussian, exponential), - prior = c(prior(normal(0, 10), Intercept, nlpar = mu1), - prior(normal(0, 1), Intercept, nlpar = mu2), - prior(normal(0, 1), nlpar = mu2)), + prior = c(prior(normal(0, 10), Intercept, dpar = mu1), + prior(normal(0, 1), Intercept, dpar = mu2), + prior(normal(0, 1), dpar = mu2)), warmup = 150, iter = 200, chains = 2, save_dso = FALSE, testmode = TRUE ) diff --git a/tests/local/tests.models_new.R b/tests/local/tests.models_new.R index 8767edf48..5fb4d23d3 100644 --- a/tests/local/tests.models_new.R +++ b/tests/local/tests.models_new.R @@ -68,9 +68,23 @@ test_that("Binomial model from brm doc works correctly", { family = binomial("probit"), cores = 2) print(fit4) - expect_message(me4 <- marginal_effects(fit4), - "median number of trials") - expect_ggplot(plot(me4, ask = FALSE)[[1]]) + expect_message( + me <- marginal_effects(fit4), "median number of trials" + ) + expect_ggplot(plot(me, ask = FALSE)[[1]]) +}) + +test_that("Non-linear model from brm doc works correctly", { + x <- rnorm(100) + y <- rnorm(100, mean = 2 - 1.5^x, sd = 1) + data5 <- data.frame(x, y) + fit5 <- brm(bf(y ~ a1 - a2^x, a1 + a2 ~ 1, nl = TRUE), + data = data5, + prior = c(prior(normal(0, 2), nlpar = a1), + prior(normal(0, 2), nlpar = a2))) + print(fit5) + me <- marginal_effects(fit5) + expect_ggplot(plot(me, ask = FALSE)[[1]]) }) test_that("Models from hypothesis doc work correctly", { @@ -106,7 +120,7 @@ test_that("bridgesampling methods work correctly", { # model with the treatment effect fit1 <- brm( count ~ log_Age_c + log_Base4_c + Trt_c, - data = epilepsy, family = negbinomial(), + data = epilepsy, family = negbinomial(), prior = prior(normal(0, 1), class = b), save_all_pars = TRUE ) @@ -114,12 +128,12 @@ test_that("bridgesampling methods work correctly", { # model without the treatment effect fit2 <- brm( count ~ log_Age_c + log_Base4_c, - data = epilepsy, family = negbinomial(), + data = epilepsy, family = negbinomial(), prior = prior(normal(0, 1), class = b), save_all_pars = TRUE ) print(fit2) - + # compute the bayes factor expect_gt(bayes_factor(fit1, fit2)$bf, 1) # compute the posterior model probabilities @@ -141,7 +155,7 @@ test_that("varying slopes without a fixed effect work", { # test reloo loo1 <- LOO(fit1) reloo1 <- reloo(loo1, fit1, chains = 1, iter = 100) - expect_range(reloo1$looic, 1620, 1700) + expect_range(reloo1$looic, 1600, 1700) conditions <- data.frame(log_Age_c = 0, log_Base4_c = 0, Trt_c = 0) me <- marginal_effects(fit1, conditions = conditions) @@ -215,7 +229,7 @@ test_that("multivariate normal models work correctly", { tim <- sample(1:N, 100) data <- data.frame(y1, y2, x, month, id, tim) - fit_mv1 <- brm(cbind(y1,y2) ~ s(x) + poly(month, 3) + (1|x|id), + fit_mv1 <- brm(cbind(y1, y2) ~ s(x) + poly(month, 3) + (1|x|id), data = data, autocor = cor_arma(~tim|id, p = 1), prior = c(prior_(~normal(0,5), resp = "y1"), prior_(~normal(0,5), resp = "y2"), @@ -296,6 +310,21 @@ test_that("Non-linear models work correctly", { expect_range(LOO(fit_loss)$looic, 700, 720) }) +test_that("Non-linear models of auxiliary parameters work correctly", { + x <- rnorm(100) + y <- rnorm(100, 1, exp(x)) + dat <- data.frame(y, x, g = rep(1:10, each = 10)) + bform <- bf(y ~ x + (1|V|g)) + + nlf(sigma ~ a, a ~ x + (1|V|g)) + + gaussian() + bprior <- prior(normal(0, 3), dpar = sigma, nlpar = a) + fit <- brm(bform, dat, prior = bprior, chains = 2, cores = 2) + print(fit) + expect_ggplot(plot(marginal_effects(fit, method = "predict"))[[1]]) + expect_equal(dim(fitted(fit, dat[1:10, ])), c(10, 4)) + expect_range(LOO(fit)$looic, 240, 350) +}) + test_that("Multivariate GAMMs work correctly", { n <- 200 sig <- 2 @@ -373,7 +402,7 @@ test_that("generalized extreme value models work correctly", { expect_ggplot(plot(me, points = TRUE, ask = FALSE)[[1]]) }) -test_that("update works correclty for some special cases", { +test_that("update works correctly for some special cases", { # models are recompiled when changing number of FEs from 0 to 1 fit1 <- brm(count ~ 1, data = epilepsy, cores = 2) fit2 <- update(fit1, ~ . + Trt_c, newdata = epilepsy, cores = 2) @@ -392,7 +421,7 @@ test_that("update works correclty for some special cases", { test_that("Wiener diffusion models work correctly", { x <- rnorm(100, mean = 1) - dat <- brms:::rwiener(n=1, alpha=2, tau=.3, beta=.5, delta=.5+x) + dat <- rwiener(n=1, alpha=2, tau=.3, beta=.5, delta=.5+x) dat$x <- x fit_d1 <- brm(bf(q | dec(resp) ~ x), dat, diff --git a/tests/testthat/tests.brm.R b/tests/testthat/tests.brm.R index b51fac9c0..fb9ce499d 100644 --- a/tests/testthat/tests.brm.R +++ b/tests/testthat/tests.brm.R @@ -28,7 +28,7 @@ test_that("brm produces expected errors", { expect_error(brm(cbind(y, x) | se(z) ~ x, dat, family = gaussian()), "allow only addition argument 'weights'") expect_error(brm(bf(y ~ x, shape ~ x), family = gaussian()), - "The parameter 'shape' is not an auxiliary parameter") + "The parameter 'shape' is not a valid distributional") expect_error(brm(y ~ x + (1|abc|g/x), dat), "Can only combine group-level terms") expect_error(brm(y ~ x + (1|g) + (x|g), dat), @@ -45,12 +45,6 @@ test_that("brm produces expected errors", { "not yet implemented for family 'poisson'") # ordinal models - expect_error(brm(rating ~ treat + period + carry + cse(treat) + (1|subject), - data = inhaler, family = cratio("logit")), - "Error occured for variables: 'treat'") - expect_error(brm(rating ~ treat + period + carry + monotonic(carry), - data = inhaler, family = cratio("logit")), - "Error occured for variables: 'carry'") expect_error(brm(rating ~ treat + (cs(period)|subject), data = inhaler, family = categorical()), "Category specific effects are only meaningful") @@ -62,7 +56,7 @@ test_that("brm produces expected errors", { "'inverse' is not a supported link for family 'poisson'") expect_error(brm(y ~ x, dat, family = c("weibull", "sqrt")), "'sqrt' is not a supported link for family 'weibull'") - expect_error(brm(y ~ x, dat, family = c("categorical","probit")), + expect_error(brm(y ~ x, dat, family = c("categorical", "probit")), "'probit' is not a supported link for family 'categorical'") expect_error(brm(y ~ x, dat, family = "ordinal"), "ordinal is not a supported family") diff --git a/tests/testthat/tests.brmsfit-methods.R b/tests/testthat/tests.brmsfit-methods.R index f6c1c1784..e1ec35c60 100644 --- a/tests/testthat/tests.brmsfit-methods.R +++ b/tests/testthat/tests.brmsfit-methods.R @@ -55,7 +55,7 @@ test_that("all S3 methods have reasonable ouputs", { expect_equal(dim(coef1$visit), c(4, 8)) coef2 <- coef(fit2, old = TRUE) expect_equal(dim(coef2[[2]]), c(59, 2)) - expect_equal(attr(coef2[[1]], "nlpar"), "a") + expect_equal(attr(coef2[[1]], "prefix")$nlpar, "a") coef4 <- coef(fit4, old = TRUE) expect_equal(dim(coef4$subject), c(10, 8)) @@ -104,16 +104,16 @@ test_that("all S3 methods have reasonable ouputs", { Age = 0, visit = c("a", "b"), Trt = 0, count = 20, patient = 1, Exp = 2 ) - fi <- fitted(fit1, auxpar = "sigma") + fi <- fitted(fit1, dpar = "sigma") expect_equal(dim(fi), c(nobs(fit1), 4)) expect_true(all(fi > 0)) - fi_lin <- fitted(fit1, auxpar = "sigma", scale = "linear") + fi_lin <- fitted(fit1, dpar = "sigma", scale = "linear") expect_equal(dim(fi_lin), c(nobs(fit1), 4)) expect_true(!isTRUE(all.equal(fi, fi_lin))) - expect_error(fitted(fit1, auxpar = "inv"), - "Invalid argument 'auxpar'") - expect_error(fitted(fit1, auxpar = "nu"), - "Auxiliary parameter 'nu' was not predicted") + expect_error(fitted(fit1, dpar = "inv"), + "Invalid argument 'dpar'") + expect_error(fitted(fit1, dpar = "nu"), + "Distributional parameter 'nu' was not predicted") fi <- fitted(fit2) expect_equal(dim(fi), c(nobs(fit2), 4)) @@ -570,11 +570,11 @@ test_that("all S3 methods have reasonable ouputs", { # do not actually refit the model as is causes CRAN checks to fail up <- update(fit1, testmode = TRUE) expect_true(is(up, "brmsfit")) + new_data <- data.frame(Age = rnorm(18), visit = rep(c(3, 2, 4), 6), Trt = rep(c(0, 0.5, -0.5), 6), count = rep(c(5, 17, 28), 6), patient = 1, Exp = 4) - up <- update(fit1, newdata = new_data, ranef = FALSE, testmode = TRUE) expect_true(is(up, "brmsfit")) expect_equal(up$data.name, "new_data") @@ -604,7 +604,7 @@ test_that("all S3 methods have reasonable ouputs", { up <- update(fit2, formula. = bf(. ~ ., a + b ~ 1, nl = TRUE), testmode = TRUE) expect_true(is(up, "brmsfit")) - up <- update(fit2, formula. = count ~ a + b, testmode = TRUE) + up <- update(fit2, formula. = bf(count ~ a + b, nl = TRUE), testmode = TRUE) expect_true(is(up, "brmsfit")) up <- update(fit3, family = acat(), testmode = TRUE) expect_true(is(up, "brmsfit")) diff --git a/tests/testthat/tests.brmsformula.R b/tests/testthat/tests.brmsformula.R index 7ae8fc356..969fdef20 100644 --- a/tests/testthat/tests.brmsformula.R +++ b/tests/testthat/tests.brmsformula.R @@ -11,7 +11,7 @@ test_that("brmsformula validates formulas of auxiliary parameters", { expect_error(bf(y ~ a, ~ 1, sigma ~ 1), "Additional formulas must be named") expect_error(bf(y ~ a^x, a ~ 1, family = gaussian()), - "The parameter 'a' is not an auxiliary parameter") + "The parameter 'a' is not a valid distributional") }) test_that("brmsformula does not change a 'brmsformula' object", { @@ -26,20 +26,20 @@ test_that("brmsformula is backwards compatible", { nonlinear = a + b ~ 1), "Argument 'nonlinear' is deprecated") expect_equivalent(pforms(form), list(a ~ 1, b ~ 1)) - expect_true(form[["nl"]]) + expect_true(attr(form$formula, "nl")) expect_warning(form <- bf(y ~ a * exp(-b * x), nonlinear = list(a ~ x, b ~ 1)), "Argument 'nonlinear' is deprecated") expect_equivalent(pforms(form), list(a ~ x, b ~ 1)) - expect_true(form[["nl"]]) + expect_true(attr(form$formula, "nl")) form <- structure(y ~ x + z, sigma = sigma ~ x) class(form) <- c("brmsformula", "formula") form <- bf(form) expect_equal(form$formula, y ~ x + z) expect_equal(pforms(form), list(sigma = sigma ~ x)) - expect_true(!form[["nl"]]) + expect_true(!attr(form$formula, "nl")) form <- structure(y ~ a * exp(-b * x), nonlinear = list(a = a ~ x, b = b ~ 1)) @@ -47,7 +47,7 @@ test_that("brmsformula is backwards compatible", { form <- bf(form) expect_equal(form$formula, y ~ a * exp(-b * x)) expect_equal(pforms(form), list(a = a ~ x, b = b ~ 1)) - expect_true(form[["nl"]]) + expect_true(attr(form$formula, "nl")) }) test_that("brmsformula detects auxiliary parameter equations", { diff --git a/tests/testthat/tests.make_stancode.R b/tests/testthat/tests.make_stancode.R index d953cd71b..408588594 100644 --- a/tests/testthat/tests.make_stancode.R +++ b/tests/testthat/tests.make_stancode.R @@ -48,7 +48,7 @@ test_that("specified priors appear in the Stan code", { prior(lkj(2), class = cor, group = g)) scode <- make_stancode( bf(y ~ a * exp(-b * x1), a + b ~ (1|ID|g), nl = TRUE), - data = dat, prior = prior,sample_prior = TRUE + data = dat, prior = prior, sample_prior = TRUE ) expect_match2(scode, "target += normal_lpdf(b_a | 0, 5)") expect_match2(scode, "target += normal_lpdf(b_b | 0, 10)") @@ -97,7 +97,7 @@ test_that("specified priors appear in the Stan code", { # test for problem described in #213 prior <- c(prior(normal(0, 1), coef = x1), - prior(normal(0, 2), coef = x1, nlpar = sigma)) + prior(normal(0, 2), coef = x1, dpar = sigma)) scode <- make_stancode(bf(y ~ x1, sigma ~ x1), dat, prior = prior) expect_match2(scode, "target += normal_lpdf(b | 0, 1);") expect_match2(scode, "target += normal_lpdf(b_sigma | 0, 2);") @@ -116,8 +116,8 @@ test_that("specified priors appear in the Stan code", { "no natural upper bound") }) -test_that("deprecated priors for the population-level intercept are used", { - dat <- data.frame(y = 1:10) +test_that("deprecated priors are used", { + dat <- data.frame(y = 1:10, x = 1:10) prior <- prior(normal(0, 5), coef = Intercept) expect_warning(scode <- make_stancode(y~1, dat, prior = prior), "Setting a prior on the population-level intercept") @@ -126,6 +126,20 @@ test_that("deprecated priors for the population-level intercept are used", { prior <- c(prior, prior(normal(0, 10), class = Intercept)) expect_error(make_stancode(y~1, dat, prior = prior), "Duplicated prior definitions detected") + + prior <- prior(normal(0, 10), nlpar = sigma) + expect_warning( + scode <- make_stancode(bf(y~1, sigma~x), dat, prior = prior), + "Specifying priors of distributional parameters via 'nlpar' is deprecated" + ) + expect_match2(scode, "normal_lpdf(b_sigma | 0, 10)") + + prior <- prior(normal(0, 5), Intercept, nlpar = y) + expect_warning( + scode <- make_stancode(cbind(y, x) ~ 1, dat, prior = prior), + "Specifying priors in multivariate models via 'nlpar' is deprecated" + ) + expect_match2(scode, "normal_lpdf(temp_y_Intercept | 0, 5)") }) test_that("special shrinkage priors appear in the Stan code", { @@ -535,7 +549,7 @@ test_that("ordinal disc parameters appear in the Stan code", { scode <- make_stancode( bf(rating ~ period + carry + treat, disc ~ period), data = inhaler, family = cumulative(), - prior = prior(normal(0,5), nlpar = disc) + prior = prior(normal(0,5), dpar = disc) ) expect_match2(scode, "target += cumulative_lpmf(Y[n] | mu[n], temp_Intercept, disc[n])" @@ -583,12 +597,32 @@ test_that("Stan code for non-linear models is correct", { flist <- list(a1 ~ 1, a2 ~ z + (x|g)) prior <- c(set_prior("beta(1,1)", nlpar = "a1", lb = 0, ub = 1), set_prior("normal(0,1)", nlpar = "a2")) - scode <- make_stancode(bf(y ~ a1 * exp(-x/(a2 + z)), - flist = flist, nl = TRUE), - data = data, family = Gamma("log"), prior = prior) + scode <- make_stancode( + bf(y ~ a1 * exp(-x/(a2 + z)), + flist = flist, nl = TRUE), + data = data, family = Gamma("log"), + prior = prior + ) expect_match2(scode, paste("mu[n] = shape * exp(-(mu_a1[n] *", "exp( - C_1[n] / (mu_a2[n] + C_2[n]))));")) + + bform <- bf(y ~ x) + + nlf(sigma ~ a1 * exp(-x/(a2 + z)), + a1 ~ 1, a2 ~ z + (x|g)) + + lf(alpha ~ x) + scode <- make_stancode( + bform, data, family = skew_normal(), + prior = c( + prior(normal(0, 1), dpar = sigma, nlpar = a1), + prior(normal(0, 5), dpar = sigma, nlpar = a2) + ) + ) + expect_match2(scode, "sigma_a1 = X_sigma_a1 * b_sigma_a1") + expect_match2(scode, + "sigma[n] = exp(sigma_a1[n] * exp( - C_sigma_1[n] / (sigma_a2[n] + C_sigma_2[n])))" + ) + expect_match2(scode, "target += normal_lpdf(b_sigma_a2 | 0, 5)") }) test_that("make_stancode accepts very long non-linear formulas", { @@ -863,8 +897,8 @@ test_that("priors on intercepts appear in the Stan code", { ) scode <- make_stancode(cbind(count, count) ~ log_Age_c + log_Base4_c * Trt_c, data = epilepsy, family = gaussian(), - prior = c(prior(student_t(5,0,10), class = b, nlpar = count), - prior(normal(0,5), class = Intercept, nlpar = count)), + prior = c(prior(student_t(5,0,10), class = b, resp = count), + prior(normal(0,5), class = Intercept, resp = count)), sample_prior = TRUE) expect_match2(scode, paste0("prior_b_count_Intercept = prior_temp_count_Intercept ", @@ -1114,14 +1148,20 @@ test_that("Stan code for Gaussian processes is correct", { "sdgp_2[2], lscale_2[2], zgp_2[Jgp_2_2]);" )) - prior <- c(prior(normal(0, 10), lscale, nlpar = eta), - prior(gamma(0.1, 0.1), sdgp, nlpar = eta), - prior(normal(0, 1), b, nlpar = eta)) - scode <- make_stancode(bf(y ~ eta, eta ~ gp(x1), nl = TRUE), + prior <- c(prior(normal(0, 10), lscale, nlpar = a), + prior(gamma(0.1, 0.1), sdgp, nlpar = a), + prior(normal(0, 1), b, nlpar = a)) + scode <- make_stancode(bf(y ~ a, a ~ gp(x1), nl = TRUE), data = dat, prior = prior) - expect_match2(scode, "target += normal_lpdf(lscale_eta_1 | 0, 10)") - expect_match2(scode, "target += gamma_lpdf(sdgp_eta_1 | 0.1, 0.1)") - expect_match2(scode, "gp(Xgp_eta_1, sdgp_eta_1[1], lscale_eta_1[1], zgp_eta_1)") + expect_match2(scode, "target += normal_lpdf(lscale_a_1 | 0, 10)") + expect_match2(scode, "target += gamma_lpdf(sdgp_a_1 | 0.1, 0.1)") + expect_match2(scode, "gp(Xgp_a_1, sdgp_a_1[1], lscale_a_1[1], zgp_a_1)") + + scode <- make_stancode(bf(y ~ a, a ~ gp(x1, by = z), nl = TRUE), + data = dat, prior = prior, silent = TRUE) + expect_match2(scode, + "mu_a[Jgp_a_1_1] = mu_a[Jgp_a_1_1] + gp(Xgp_a_1[Jgp_a_1_1]," + ) }) test_that("Stan code for SAR models is correct", { diff --git a/tests/testthat/tests.make_standata.R b/tests/testthat/tests.make_standata.R index ed00c015e..39bd96938 100644 --- a/tests/testthat/tests.make_standata.R +++ b/tests/testthat/tests.make_standata.R @@ -181,7 +181,7 @@ test_that("make_standata rejects incorrect addition terms", { expect_error(make_standata(y | cens(c) ~ 1, data = dat)) expect_error(make_standata(z | trials(t) ~ 1, data = dat, family = "binomial"), - "Number of trials is smaller than the response variable") + "Number of trials is smaller than the number of events") }) test_that("make_standata handles multivariate models", { @@ -296,6 +296,16 @@ test_that("make_standata correctly prepares data for non-linear models", { sdata <- make_standata(bf(y ~ a - b^z, flist = flist, nl = TRUE), data = data, control = list(not4stan = TRUE)) expect_equal(colnames(sdata$C), "z") + + bform <- bf(y ~ x) + + nlf(sigma ~ a1 * exp(-x/(a2 + z)), + a1 ~ 1, a2 ~ z + (x|g)) + + lf(alpha ~ x) + sdata <- make_standata(bform, data, family = skew_normal()) + sdata_names <- c("C_sigma_1", "C_sigma_2", "X_sigma_a2", "Z_1_sigma_a2_1") + expect_true( + all(sdata_names %in% names(sdata)) + ) }) test_that("make_standata correctly prepares data for monotonic effects", { @@ -309,7 +319,7 @@ test_that("make_standata correctly prepares data for monotonic effects", { expect_equal(sdata$con_simplex_1, rep(1, 3)) prior <- set_prior("dirichlet(1:3)", coef = "x1", - class = "simplex", nlpar = "sigma") + class = "simplex", dpar = "sigma") sdata <- make_standata(bf(y ~ 1, sigma ~ mono(x1)), data = data, prior = prior) expect_equal(sdata$con_simplex_sigma_1, 1:3) @@ -463,7 +473,7 @@ test_that("make_standata allows fixed auxiliary parameters", { expect_equal(make_standata(bf(y ~ 1, nu = 3), dat, student())$nu, 3) expect_equal(make_standata(y ~ 1, dat, acat())$disc, 1) expect_error(make_standata(bf(y ~ 1, bias = 0.5), dat), - "Invalid auxiliary parameters: 'bias'") + "Invalid fixed parameters: 'bias'") }) test_that("make_standata correctly includes offsets", { diff --git a/tests/testthat/tests.priors.R b/tests/testthat/tests.priors.R index 07a4ec4f6..d1c41a236 100644 --- a/tests/testthat/tests.priors.R +++ b/tests/testthat/tests.priors.R @@ -56,7 +56,7 @@ test_that("get_prior returns correct prior names for auxiliary parameters", { dat <- data.frame(y = rnorm(10), x = rnorm(10), z = rnorm(10), g = rep(1:2, 5)) prior <- get_prior(bf(y ~ 1, phi ~ z + (1|g)), data = dat, family = Beta()) - prior <- prior[prior$nlpar == "phi", ] + prior <- prior[prior$dpar == "phi", ] pdata <- data.frame(class = c("b", "b", "b", "Intercept", rep("sd", 3)), coef = c("", "Intercept", "z", "", "", "", "Intercept"), group = c(rep("", 5), "g", "g"), @@ -71,16 +71,16 @@ test_that("get_prior returns global priors in multivariate models", { # MV normal prior <- get_prior(cbind(y1, y2) ~ x + (x|ID1|g), data = dat, family = gaussian()) - expect_equal(prior[prior$nlpar == "y1" & prior$class == "b", "coef"], + expect_equal(prior[prior$resp == "y1" & prior$class == "b", "coef"], c("", "Intercept", "x")) - expect_equal(prior[prior$nlpar == "" & prior$class == "sd", "prior"], + expect_equal(prior[prior$resp == "" & prior$class == "sd", "prior"], c("student_t(3, 0, 10)")) # categorical prior <- get_prior(y2 ~ x + (x|ID1|g), data = dat, family = categorical()) - expect_equal(prior[prior$nlpar == "X2" & prior$class == "b", "coef"], + expect_equal(prior[prior$resp == "X2" & prior$class == "b", "coef"], c("", "Intercept", "x")) - expect_equal(prior[prior$nlpar == "" & prior$class == "sd", "prior"], + expect_equal(prior[prior$resp == "" & prior$class == "sd", "prior"], c("student_t(3, 0, 10)")) }) diff --git a/tests/testthat/tests.rename.R b/tests/testthat/tests.rename.R index d0692e7b7..a7eb797d0 100644 --- a/tests/testthat/tests.rename.R +++ b/tests/testthat/tests.rename.R @@ -18,7 +18,7 @@ test_that("rm_int_fe works as expected", { code <- make_stancode(y ~ 1, data = dat) expect_equal(rm_int_fe("Intercept", code), character(0)) code <- make_stancode(cbind(y, y) ~ x, data = dat) - expect_equal(rm_int_fe(c("Intercept", "x"), code, nlpar = "y"), "x") + expect_equal(rm_int_fe(c("Intercept", "x"), code, px = list(resp = "y")), "x") code <- make_stancode(y ~ x, data = dat, family = sratio()) expect_equal(rm_int_fe(c("Intercept", "x"), code), "x") code <- make_stancode(y ~ 0 + intercept + x, data = dat) @@ -42,7 +42,8 @@ test_that("change_prior returns expected lists", { test_that("change_old_re and change_old_re2 return expected lists", { data <- data.frame(y = rnorm(10), x = rnorm(10), g = 1:10) - bterms <- parse_bf(bf(y ~ a, a ~ x + (1+x|g), nl = TRUE)) + bterms <- parse_bf(bf(y ~ a, a ~ x + (1+x|g), + family = gaussian(), nl = TRUE)) ranef <- brms:::tidy_ranef(bterms, data = data) target <- list( list(pos = c(rep(FALSE, 2), TRUE, rep(FALSE, 22)), diff --git a/tests/testthat/tests.validate.R b/tests/testthat/tests.validate.R index 8b0dde30a..df20c7696 100644 --- a/tests/testthat/tests.validate.R +++ b/tests/testthat/tests.validate.R @@ -11,13 +11,13 @@ test_that("parse_bf handles very long RE terms", { formula <- paste(sprintf("y ~ 0 + trait + trait:(%s)", covariate_vector), sprintf("(1+%s|id)", covariate_vector), sep = " + ") bterms <- parse_bf(as.formula(formula)) - expect_equal(bterms$auxpars$mu$re$group, "id") + expect_equal(bterms$dpars$mu$re$group, "id") }) test_that("parse_bf correctly handles auxiliary parameter 'mu'", { bterms1 <- parse_bf(y ~ x + (x|g)) bterms2 <- parse_bf(bf(y~1, mu ~ x + (x|g))) - expect_equal(bterms1$auxpars$mu, bterms2$auxpars$mu) + expect_equal(bterms1$dpars$mu, bterms2$dpars$mu) expect_error(parse_bf(bf(y ~ z, mu ~ x + (x|g))), "All 'mu' parameters are specified") }) @@ -94,6 +94,6 @@ test_that("exclude_pars returns expected parameter names", { bterms <- parse_bf(y ~ x + s(z)) data <- data.frame(y = rnorm(20), x = rnorm(20), z = rnorm(20)) expect_true("zs_1_1" %in% exclude_pars(bterms, data)) - bterms <- parse_bf(bf(y ~ eta, eta ~ x + s(z), nl = TRUE)) + bterms <- parse_bf(bf(y ~ eta, eta ~ x + s(z), family = gaussian(), nl = TRUE)) expect_true("zs_eta_1_1" %in% exclude_pars(bterms, data)) })