diff --git a/README.md b/README.md index 87e639d..ebb4ca8 100644 --- a/README.md +++ b/README.md @@ -223,6 +223,41 @@ In the figure "Negative binomial example" above, suppose we had decided that `ca For a more in-depth discussion of the relationship between sensitivity and robustness, see Appendix C of our paper [3]. +# Development + +The following commands are to be run from the root of a clone of the +repository. + +To install a local version, in `R`, run + +``` +> library(devtools) +> install_local("rstansensitivity", force=TRUE) +``` + +To run the R unit tests, execute the following in your shell: + +``` +$ cd rstansensitivity/tests +$ ./run_tests.R +``` + +To run the python tests, execute the following in your shell: + +``` +$ rstansensitivity/inst/generate_models_unittest.py +``` + +To build the documentation (and exports), run the following commands in `R` from +the `rstansensitivity` directory: + +``` +> library(roxygen2) +> library(pkgdown) +> roxygen2::roxygenize() +> pkgdown::build_site() +``` + # References [1]: Local posterior robustness with parametric priors: Maximum and average sensitivity, Basu, Jammalamadaka, Rao, Liu (1996) diff --git a/rstansensitivity/DESCRIPTION b/rstansensitivity/DESCRIPTION index 824c290..5aef0dd 100644 --- a/rstansensitivity/DESCRIPTION +++ b/rstansensitivity/DESCRIPTION @@ -1,12 +1,12 @@ Package: rstansensitivity Title: Tools for calculating hyperparameter sensitivity in Stan. -Version: 0.0.0.9000 -Authors@R: person("Ryan", "Giordano", email = "rgiordano@berkeley.edu", role = c("aut", "cre")) -Maintainer@R: person("Ryan", "Giordano", email = "rgiordano@berkeley.edu", role = c("aut", "cre")) +Version: 0.0.1.0000 +Authors@R: person("Ryan", "Giordano", email = "rgiordan@mit.edu", role = c("aut", "cre")) +Maintainer@R: person("Ryan", "Giordano", email = "rgiordan@mit.edu", role = c("aut", "cre")) Description: Use Stan's autodifferentiation capabilities to automatically calculate Monte Carlo estimates of local sensitivity of posterior means to hyperparameters. -Depends: R (>= 3.2.3), rstan, dplyr, reshape2 +Depends: R (>= 3.2.3), rstan, dplyr, reshape2, tidybayes License: ALv2 Encoding: UTF-8 LazyData: true -Suggests: testthat -RoxygenNote: 6.0.1 +Suggests: testthat, pkgdown, roxygen2 +RoxygenNote: 6.1.1 diff --git a/rstansensitivity/NAMESPACE b/rstansensitivity/NAMESPACE index 5bb2983..d87bd01 100644 --- a/rstansensitivity/NAMESPACE +++ b/rstansensitivity/NAMESPACE @@ -7,14 +7,19 @@ export(GetHyperparameterDataFrame) export(GetImportanceSamplingFromModelFit) export(GetMCMCDataFrame) export(GetSamplingModelFilename) +export(GetSensitivityFromGrads) export(GetSensitivityStandardErrors) export(GetStanSensitivityFromModelFit) +export(GetStanSensitivityMatricesFromModelFit) export(GetStanSensitivityModel) export(GetTidyResult) +export(GetTidySensitivityResults) +export(GetTidyStanSummary) +export(GetWeightMatrixSecondDerivative) +export(LoadStanData) export(PlotSensitivities) export(PredictSensitivityFromStanData) export(SetSensitivityParameterList) export(SummarizeSensitivityMatrices) +export(TidyColumn) export(TransformSensitivityResult) -export(GetCovarianceSE) -export(GetNormalizedCovarianceSE) diff --git a/rstansensitivity/R/covariance_standard_errors_lib.R b/rstansensitivity/R/covariance_standard_errors_lib.R index 1400da0..b5277e5 100644 --- a/rstansensitivity/R/covariance_standard_errors_lib.R +++ b/rstansensitivity/R/covariance_standard_errors_lib.R @@ -1,3 +1,5 @@ +# Utilities using `mcmcse` to evaluate standard errors for correlations. + library(mcmcse) @@ -11,7 +13,10 @@ GetCovarianceSE <- function(x_draws, y_draws, fix_mean=FALSE) { arg_draws <- cbind(x_draws * y_draws, x_draws, y_draws) arg_cov_mat <- mcmcse::mcse.multi(arg_draws)$cov grad_g <- c(1, -1 * y_mean, -1 * x_mean) - g_se <- as.numeric(sqrt(t(grad_g) %*% arg_cov_mat %*% grad_g / nrow(arg_draws))) + g_se <- as.numeric( + sqrt(t(grad_g) %*% + arg_cov_mat %*% + grad_g / nrow(arg_draws))) } return(g_se) } @@ -27,6 +32,7 @@ PackNormalizedCovariancePar <- function(x_draws, y_draws) { } GetNormalizedCovarianceGradient <- function(par) { + # We use this to apply the delta method to normalizes sensitivities. mean_xy <- par[1] mean_x2 <- par[2] mean_x <- par[3] @@ -76,6 +82,7 @@ GetNormalizedCovarianceSE <- function(x_draws, y_draws) { GetSensitivityStandardErrors <- function( draws_mat, grad_mat, fix_mean=FALSE, normalized=FALSE) { + # TODO: fix the order of these arguments to match everything else. if (normalized && fix_mean) { stop("You cannot specify both fix_mean and normalized.") } diff --git a/rstansensitivity/R/result_processing_lib.R b/rstansensitivity/R/result_processing_lib.R index 88ff24d..d9bba80 100644 --- a/rstansensitivity/R/result_processing_lib.R +++ b/rstansensitivity/R/result_processing_lib.R @@ -1,8 +1,25 @@ -library(ggplot2) -library(dplyr) -library(reshape2) + +CheckDimensions <- function(grad_mat, draws_mat) { + if (ncol(grad_mat) != nrow(draws_mat)) { + stop(paste0( + "The dimensions of `grad_mat` must be ", + "# hyper parameters x # MCMC draws, and the dimensions ", + "of draws_mat must be # draws x # parameters.")) + } +} +#' Get a sensitivity matrix from a matrix of partial derivatives and +#' a matrix of draws. +#' @param grad_mat A matrix of gradients of the log probability with respect +#' to hyperparameters. The rows are hyperparameters, and the columns are +#' draws. +#' @param draws_mat A matrix of draws of parameters, the means of which +#' are the quantities of interest. The rows are draws, and the columns +#' are parameters. +#' @return A matrix of estimated derivatives, dE[parameter | hyper] / dhyper. +#' The rows are hyperparameters and the columns are parameters. +#' @export GetSensitivityFromGrads <- function(grad_mat, draws_mat) { # This should in fact match cov() but without having to transpose, # which gives a speedup. @@ -10,6 +27,8 @@ GetSensitivityFromGrads <- function(grad_mat, draws_mat) { # Should match: # sens_mat <- cov(t(grad_mat), draws_mat) + CheckDimensions(grad_mat, draws_mat) + grad_means <- rowMeans(grad_mat) draw_means <- colMeans(draws_mat) @@ -53,183 +72,3 @@ TransformSensitivityResult <- function(sens_result, transform_matrix) { transform_matrix %*% transform_sens_result$grad_mat return(transform_sens_result) } - - -SensitivityMatrixToDataframe <- function( - sens_mat, hyperparameter_names, parameter_names) { - colnames(sens_mat) <- parameter_names - return(data.frame(sens_mat) %>% - mutate(hyperparameter=hyperparameter_names) %>% - melt(id.var="hyperparameter") %>% - rename(parameter=variable)) -} - - -#' Make a tidy dataframe out of a sensitivity matrix and its standard errors. -#' @param sens_mat A matrix of sensitivities. -#' @param sens_se A matrix of standard errors of \code{sens_mat}. -#' @param measure What to call these sensitivites. -#' @param num_se The number of standard errors for the upper and lower bounds. -#' @return A tidy dataframe with columns for the parameters, hyperparameters, -#' sensitivities, and their standard errors. -#' @export -SummarizeSensitivityMatrices <- function( - sens_mat, sens_se, measure, num_se=2) { - - sens_df <- rbind( - SensitivityMatrixToDataframe( - sens_mat, - hyperparameter_names=rownames(sens_mat), - parameter_names=colnames(sens_mat)) %>% - mutate(measure=measure), - SensitivityMatrixToDataframe( - sens_se, - hyperparameter_names=rownames(sens_se), - parameter_names=colnames(sens_se)) %>% - mutate(measure=paste(measure, "se", sep="_")), - SensitivityMatrixToDataframe( - sens_mat + num_se * sens_se, - hyperparameter_names=rownames(sens_mat), - parameter_names=colnames(sens_mat)) %>% - mutate(measure=paste(measure, "upper", sep="_")), - SensitivityMatrixToDataframe( - sens_mat - num_se * sens_se, - hyperparameter_names=rownames(sens_mat), - parameter_names=colnames(sens_mat)) %>% - mutate(measure=paste(measure, "lower", sep="_")) - ) %>% - dcast(hyperparameter + parameter ~ measure) %>% - filter(parameter != "lp__") - return(sens_df) -} - -#' Process the results of GetStanSensitivityFromModelFit into a -#' tidy dataframe with standard errors. -#' -#' @param sens_result The output of \code{GetStanSensitivityFromModelFit}. -#' @param num_se The number of standard errors for the upper and lower bounds. -#' @return A dataframe summarizing the sensitivity of the model posterior means -#' to the hyperparameters. The reported standard errors are based on a -#' multivariate normal and delta method approximation which may not be -#' accurate for highly variable or non-normal draws. The standard errors should -#' be interpreted with caution. -#' \itemize{ -#' \item{hyperparameter: }{The name of the hyperparameter.} -#' \item{parameter: }{The name of the parameter.} -#' \item{sensitivity:} -#' {The MCMC estimate of \code{dE[parameter | X, hyperparameter] / d hyperparameter}.} -#' \item{sensitivity_se: } -#' {The estimated Monte Carlo error of \code{sensitivity}.} -#' \item{normalized_sensitivity: } -#' {The MCMC estimate of \code{sensitivity / sd(parameter}).} -#' \item{normalized_sensitivity_se: } -#' {The estimated Monte Carlo error of \code{normalized_sensitivity}.} -#' } -#' @export -GetTidyResult <- function(sens_result, num_se=2) { - sens_se <- GetSensitivityStandardErrors( - sens_result$draws_mat, sens_result$grad_mat, - fix_mean=FALSE, normalized=FALSE) - norm_sens_se <- GetSensitivityStandardErrors( - sens_result$draws_mat, sens_result$grad_mat, - fix_mean=FALSE, normalized=TRUE) - sens_df <- SummarizeSensitivityMatrices( - sens_result$sens_mat, sens_se, - measure="sensitivity", num_se=num_se) - sens_norm_df <- SummarizeSensitivityMatrices( - sens_result$sens_mat_normalized, norm_sens_se, - measure="normalized_sensitivity", num_se=num_se) - - result <- inner_join( - sens_df,sens_norm_df, by=c("hyperparameter", "parameter")) - - return(result) -} - - -#' Make a barchart from the output of \code{GetTidyResult}. -#' -#' @param sens_df The output of \code{GetTidyResult}. -#' @param normalized Whether to display the normalized sensitivities. -#' @param se_num The number of standard errors to plot with the errorbars. -#' @return A \code{ggplot} object showing the sensitivites of each parameter -#' to each hyperparameter. -#' @export -PlotSensitivities <- function(sens_df, normalized=TRUE, se_num=2) { - if (normalized) { - y_axis_label <- "Normalized sensitivity" - sens_df_plot <- - mutate(sens_df, s=normalized_sensitivity, - s_se=normalized_sensitivity_se) - } else { - y_axis_label <- "Unnormalized sensitivity" - sens_df_plot <- - mutate(sens_df, s=sensitivity, s_se=sensitivity_se) - } - return( - ggplot(sens_df_plot) + - geom_bar(aes(x=parameter, y=s, fill=hyperparameter), - stat="identity", position="dodge") + - geom_errorbar(aes(x=parameter, - ymin=s - se_num * s_se, - ymax=s + se_num * s_se, - group=hyperparameter), - position=position_dodge(0.9), width=0.2) + - theme(axis.text.x = element_text(angle = 90, hjust = 1)) + - scale_fill_discrete(name="Hyperparameter") + - ylab(y_axis_label) + xlab("Parameter") - ) -} - - -#' Convert a Stan sampling result to a form that can be joined with tidy -#' sensitivity results. -#' -#' @param sampling_result The output of \code{stan::sampling} -#' @param cols The columns of the summary to keep. -#' @return A data frame with the MCMC reuslts that can be joined with tidy -#' sensitivity results. -#' @export -GetMCMCDataFrame <- function( - sampling_result, cols=c("mean", "se_mean", "sd", "n_eff", "Rhat")) { - mcmc_result <- as.data.frame(rstan::summary(sampling_result)$summary) - mcmc_result$parameter <- make.names(rownames(mcmc_result)) - rownames(mcmc_result) <- NULL - mcmc_result <- - select(mcmc_result, "parameter", cols) %>% - filter(parameter != "lp__") - return(mcmc_result) -} - - -#' Use the linear approximation to predict the sensitivity to a new -#' stan data list. -#' -#' @param stan_sensitivity_list The output of \code{GetStanSensitivityModel} -#' @param stan_result The output of \code{GetStanSensitivityFromModelFit} -#' @param stan_data The original stan data at which -#' \code{stan_sensitivity_list} was calculated. -#' @param stan_data_perturb A new stan data file with different hyperparameters. -#' @param description A hyperparameter name to describe this perturbation. -#' @return A tidy sensitivity dataframe where the sensitivity is in the -#' direction of the difference between the hyperparameters in the two stan -#' data lists. -#' @export -PredictSensitivityFromStanData <- function( - stan_sensitivity_list, sens_result, stan_data, stan_data_perturb, - description="perturbation") { - - hyperparameter_df <- - inner_join( - GetHyperparameterDataFrame(stan_sensitivity_list, stan_data) %>% - rename(hyperparameter_val_orig=hyperparameter_val), - GetHyperparameterDataFrame(stan_sensitivity_list, stan_data_perturb), - by="hyperparameter") %>% - mutate(hyperparameter_diff=hyperparameter_val - hyperparameter_val_orig) - linear_comb <- matrix(hyperparameter_df$hyperparameter_diff, nrow=1) - colnames(linear_comb) <- hyperparameter_df$hyperparameter - linear_comb <- linear_comb[, rownames(sens_result$sens_mat), drop=FALSE] - rownames(linear_comb) <- description - sens_result_pert <- TransformSensitivityResult(sens_result, linear_comb) - return(GetTidyResult(sens_result_pert)) -} diff --git a/rstansensitivity/R/reweighting_lib.R b/rstansensitivity/R/reweighting_lib.R new file mode 100644 index 0000000..802b648 --- /dev/null +++ b/rstansensitivity/R/reweighting_lib.R @@ -0,0 +1,87 @@ +# These are experimental ways of doing a different design. i'm not going to +# export any for the moment, but I'll keep them here as a placeholder. + +library(tidybayes) +library(dplyr) + + +GetWeightMatrixSecondDerivative <- function(stanfit, log_lik_name="log_lik") { + log_lik <- t(as.matrix(stanfit, log_lik_name)) + log_lik <- log_lik - rowMeans(log_lik) + + SecondDerivative <- function(w1, w2, draws, base_w=1) { + log_lik1 <- t(w1 - base_w) %*% log_lik + log_lik2 <- t(w2 - base_w) %*% log_lik + d2dw_mat <- log_lik1 * log_lik2 + d2dw_mat <- d2dw_mat - mean(d2dw_mat) + d2dw_mat %*% draws / nrow(draws) + } + + return(SecondDerivative) +} + + +GetWeightMatrixDerivative <- function(stanfit, log_lik_name="log_lik") { + log_lik <- t(as.matrix(stanfit, log_lik_name)) + log_lik <- log_lik - rowMeans(log_lik) + + Derivative <- function(draws) { + log_lik %*% draws / nrow(draws) + } + + return(Derivative) +} + + +GetWeightMatrixPredictor <- function(stanfit, log_lik_name="log_lik") { + Derivative <- GetWeightMatrixDerivative( + stanfit=stanfit, log_lik_name=log_lik_name) + PredictDiff <- function(w, draws, base_w=1) { + t(w - base_w) %*% Derivative(draws) + } + + return(PredictDiff) +} + + +GetTidyWeightPredictor <- function(..., + stanfit=NULL, + PredictDiff=NULL, + draws=NULL, + base_w=1) { + pars <- enquos(...) + if (!xor(is.null(stanfit), is.null(PredictDiff))) { + stop("You must specify `stanfit` or `PredictDiff` but not both.") + } + if (!xor(is.null(stanfit), is.null(draws))) { + stop("You must specify `stanfit` or `draws` but not both.") + } + if (!is.null(stanfit)) { + PredictDiff <- GetWeightMatrixPredictor(stanfit) + # Take the quoted variables, convert to strings, and strip the index + # for passing to Stan directly, + par_names <- + lapply(pars, as_label) %>% + sapply(function(x) { sub("\\[.*$", "", x) }) + draws <- as.matrix(stanfit, par_names) + } + par_means <- + gather_draws(t(colMeans(draws)), !!!pars) %>% + RemoveExtraTidyColumns() %>% + rename(base_mean=.value) + join_vars <- setdiff(names(par_means), "base_mean") + + Predict <- function(w) { + pred_mat <- PredictDiff(w=w, draws=draws, base_w=base_w) + pred_df <- + inner_join( + gather_draws(pred_mat, !!!pars) %>% + rename(pred_diff=.value) %>% + RemoveExtraTidyColumns(), + par_means, + by=join_vars) %>% + mutate(pred_mean=base_mean + pred_diff) + return(pred_df) + } + return(Predict) +} diff --git a/rstansensitivity/R/stan_sensitivity_lib.R b/rstansensitivity/R/stan_sensitivity_lib.R index a84f9b0..25a2367 100644 --- a/rstansensitivity/R/stan_sensitivity_lib.R +++ b/rstansensitivity/R/stan_sensitivity_lib.R @@ -1,6 +1,7 @@ +# Tools for getting `grad_mat` (the matrix of partial derivatives +# of the log posterior with respect to the hyperparameters) out of Stan. + library(rstan) -library(dplyr) -library(reshape2) # Just a more readable shortcut for the Stan attribute. @@ -17,10 +18,7 @@ GetParamNames <- function(model_fit) { #' #' @return The full path of the generated file to be used for sampling. #' @export -#' @examples -#' model_name <- GenerateSensitivityFromModel("models/my_model.stan") -#' model <- stan_model(GetSamplingModelFilename(model_name)) -#' sampling_result <- sampling(model, data=stan_data) +#' @example inst/examples/GenerateSensitivityFromModel.R GetSamplingModelFilename <- function(model_name) { return(paste(model_name, "_generated.stan", sep="")) } @@ -36,15 +34,16 @@ GetSamplingModelFilename <- function(model_name) { #' functions in this library. The function also generates the sensitivity model #' files in the same location as the original base model. #' @export -#' @examples -#' GenerateSensitivityFromModel("models/my_model.stan") +#' @example inst/examples/GenerateSensitivityFromModel.R GenerateSensitivityFromModel <- function( base_model_name, python_script=system.file("generate_models.py", package="rstansensitivity")) { model_suffix <- - substr(base_model_name, nchar(base_model_name) - 4, nchar(base_model_name)) + substr(base_model_name, + nchar(base_model_name) - 4, + nchar(base_model_name)) stopifnot(model_suffix == ".stan") system(paste("python ", python_script, " --base_model=", base_model_name, sep="")) @@ -59,14 +58,15 @@ GenerateSensitivityFromModel <- function( #' @param model_sens_params A stanfit object with the sensitivity model #' parameters. In order to guarantee legal initial values, it may be #' preferable to use an empty model block. -#' @param model_sens_params A stanfit object for the original model. +#' @param model_params A stanfit object for the original model. #' @param stan_data The stan data list for the original model with #' hyperparameters specified for reading in the data block. Each hyperparameter -#' in the model_sens_params stanfit model must be specified in \code{stan_data}. +#' in the \code{model_sens_params} stanfit model must be specified in +#' \code{stan_data}. #' @return A list of valid model parameters than can be passed to the -#' sensivitiy model, in which the sampled parameters are taken from +#' sensivity model, in which the sampled parameters are taken from #' \code{model_params} and the hyperparameters are taken from \code{stan_data}. -#' @export +##' @export SetSensitivityParameterList <- function( model_sens_params, model_params, stan_data) { @@ -118,9 +118,10 @@ GetStanSensitivityModel <- function(model_name, stan_data) { # This is the stan fit object that we will use to evaluate the gradient # of the log probability at the MCMC draws. - model_sens_fit <- stan(paste(model_name, "_sensitivity.stan", sep=""), - data=stan_data, algorithm="Fixed_param", - iter=1, chains=1, init=list(sens_par_list)) + model_sens_fit <- rstan::stan( + paste(model_name, "_sensitivity.stan", sep=""), + data=stan_data, algorithm="Fixed_param", + iter=1, chains=1, init=list(sens_par_list)) # These names help sort through the vectors of sensitivity. param_names <- GetParamNames(model_params) @@ -136,67 +137,94 @@ GetStanSensitivityModel <- function(model_name, stan_data) { #' Evaluate a model (represented as a stanfit object) at MCMC draws, possibly -#' from a different model. +#' from a different model. The model \code{model_stanfit} is evaluated +#' at the draws found in \code{samples_stanfit}. If there are parameters in +#' \code{model_stanfit} that are not contained in \code{samples_stanfit}, +#' values from \code{model_par_list} are used. +#' +#' This was designed to evaluate gradients with respect to hyperparameters +#' which are Stan parameters in \code{model_stanfit} but not in +#' \code{samples_stanfit}. The values of those hyperparameters would be +#' specified in \code{model_par_list} and reused for each evaluation of +#' \code{model_stanfit}. Since this is the intended use case, the function +#' requires that \code{model_stanfit} contains a (non-strict) superset of the +#' parameters in \code{samples_stanfit}. #' -#' @param sampling_result The output of Stan sampling (e.g. rstan::sampling()) -#' @param model_fit A stanfit object from a model, possibly containing a -#' superset of the parameters in sampling_result. -#' @param A list of parameters for model_fit which can be passed to -#' \code{unconstrain_pars} for the \code{model_fit} model. Every parameter -#' in \code{get_inits} applied to \code{sampling_result} must be also be found +#' @param samples_stanfit (sampling_result) The output of Stan sampling +#' containing the draws at which you want to evaluate the model. +#' (e.g. \code{rstan::sampling()}) +#' @param model_stanfit (model_fit) A stanfit object from a model, possibly containing a +#' superset of the parameters in samples_stanfit. +#' @param model_par_list A list of parameters for model_stanfit which can be +#' passed to +#' \code{unconstrain_pars} for the \code{model_stanfit} model. Every parameter +#' in \code{get_inits} applied to \code{samples_stanfit} must be also be found #' in \code{model_par_list}. #' @param compute_grads If FALSE, only return the log probability. +#' @param max_chains How many chains to evaluate. The default is to +#' evaluate the number of chains in \code{samples_stanfit}. +#' @param max_num_samples How many samples to evaluate. The default is to +#' evaluate the number of non-warmup samples in \code{samples_stanfit}. +#' @param verbose If true, display progress messages. #' @return #' A list with \code{lp_vec} containing a vector of log probabilities and, #' if \code{compute_grads} is \code{TRUE}, gradients of the log probability, -#' all with resepct to the parameters in \code{model_fit} at the draws in -#' \code{sampling_result}. If the sampler contains multiple chains the +#' all with resepct to the parameters in \code{model_stanfit} at the draws in +#' \code{samples_stanfit}. If the sampler contains multiple chains the #' chains are concatenated in order. -#' @export +##' @export EvaluateModelAtDraws <- function( - sampling_result, model_fit, model_par_list, - compute_grads=FALSE, max_chains=Inf, max_num_samples=Inf) { + samples_stanfit, model_stanfit, model_par_list, + compute_grads=FALSE, max_chains=Inf, max_num_samples=Inf, + verbose=TRUE) { - num_warmup_samples <- sampling_result@sim$warmup - num_samples <- min(sampling_result@sim$iter - num_warmup_samples, + num_warmup_samples <- samples_stanfit@sim$warmup + num_samples <- min(samples_stanfit@sim$iter - num_warmup_samples, max_num_samples) - num_chains <- min(sampling_result@sim$chains, max_chains) + num_chains <- min(samples_stanfit@sim$chains, max_chains) - # Check that every parameter in the sampling results is also in the model - # parameter list. - par_list_check <- get_inits(sampling_result, iter=1)[[1]] - missing_pars <- setdiff(names(par_list_check), names(model_par_list)) + # Check that every parameter in the samples_stanfit is also in model_stanfit. + samples_stanfit_pars <- rstan::get_inits(samples_stanfit, iter=1)[[1]] + missing_pars <- setdiff(names(samples_stanfit_pars), names(model_par_list)) if (length(missing_pars) > 0) { - stop(sprintf( - "Parameters from sampling result not in new model: %s", - paste(missing_pars, collapse=", "))) + err_msg <- paste0( + "Every parameter in `samples_stanfit` must also be found in ", + "`model_stanfit`. ", + sprintf( + "Parameters from sampling result not in new model: %s", + paste(missing_pars, collapse=", "))) + stop(err_msg) } # Get the model gradients with respect to the hyperparameters (and parameters). lp_vec <- rep(NA, num_samples) if (compute_grads) { - param_names <- GetParamNames(model_fit) - grad_mat <- matrix(NA, length(param_names), num_samples * num_chains) + param_names <- GetParamNames(model_stanfit) + grad_mat <- matrix( + NA, nrow=length(param_names), ncol=num_samples * num_chains) rownames(grad_mat) <- param_names } else { grad_mat <- matrix() } for (chain in 1:num_chains) { - cat("Evaluating model at the MCMC draws for chain ", - chain, ".\n") - prog_bar <- txtProgressBar(min=1, max=num_samples, style=3) + if (verbose) { + cat("Evaluating model at the MCMC draws for chain ", chain, ".\n") + prog_bar <- txtProgressBar(min=1, max=num_samples, style=3) + } for (n in 1:num_samples) { - setTxtProgressBar(prog_bar, value=n) + if (verbose) { + setTxtProgressBar(prog_bar, value=n) + } # We rely on get_inits to return the draws at iteration n in a form # that is easy to parse. par_list <- get_inits( - sampling_result, iter=n + num_warmup_samples)[[chain]] + samples_stanfit, iter=n + num_warmup_samples)[[chain]] for (par in ls(par_list)) { model_par_list[[par]] <- par_list[[par]] } - pars_free <- unconstrain_pars(model_fit, model_par_list) + pars_free <- rstan::unconstrain_pars(model_stanfit, model_par_list) # The index in the matrix stacked by chain. # This needs to match the stacking done to the Stan samples in @@ -205,14 +233,16 @@ EvaluateModelAtDraws <- function( # put them into compatible shapes... ix <- (chain - 1) * num_samples + n if (compute_grads) { - glp <- grad_log_prob(model_fit, pars_free) + glp <- rstan::grad_log_prob(model_stanfit, pars_free) grad_mat[, ix] <- glp lp_vec[ix] <- attr(glp, "log_prob") } else { - lp_vec[ix] <- log_prob(model_fit, pars_free) + lp_vec[ix] <- rstan::log_prob(model_stanfit, pars_free) } } - close(prog_bar) + if (verbose) { + close(prog_bar) + } } return(list(lp_vec=lp_vec, grad_mat=grad_mat)) } @@ -239,12 +269,12 @@ CheckModelEquivalence <- function( check_grads=TRUE, tol=1e-8, max_chains=1, max_num_samples=100) { - model_fit_1 <- sampling( + model_fit_1 <- rstan::sampling( object=model_1, data=stan_data_1, algorithm="Fixed_param", iter=1, chains=1) model_par_list_1 <- get_inits(model_fit_1, 1)[[1]] - model_fit_2 <- sampling( + model_fit_2 <- rstan::sampling( object=model_2, data=stan_data_2, algorithm="Fixed_param", iter=1, chains=1) model_par_list_2 <- get_inits(model_fit_2, 1)[[1]] @@ -327,15 +357,12 @@ StackChainArray <- function(draws_array) { #' @param stan_sensitivity_list The output of \code{GetStanSensitivityModel} #' @return A list of matrices. The elements of the list are #' \itemize{ -#' \item{sens_mat: }{The local sensitivity of each posterior parameter to each hyperparameter.} -#' \item{sens_mat_normalized: }{The same quantities as \code{sens_mat}, but normalized by the posterior standard deviation.} #' \item{grad_mat: }{The gradients of the log posterior evaluatated at the draws.} #' \item{lp_vec: }{The log probability of the model at each draw.} #' \item{draws_mat: }{The parameter draws in the same order as that of \code{grad_mat}}. #' } -#' The result can be used directly, or passed to \code{GetTidyResult}. #' @export -GetStanSensitivityFromModelFit <- function( +GetStanSensitivityMatricesFromModelFit <- function( sampling_result, stan_sensitivity_list) { draws_mat <- StackChainArray(extract(sampling_result, permute=FALSE)) @@ -353,15 +380,43 @@ GetStanSensitivityFromModelFit <- function( grad_mat <- model_at_draws$grad_mat[ setdiff(sens_param_names, param_names),, drop=FALSE] + return(list(grad_mat=grad_mat, + lp_vec=model_at_draws$lp_vec, + draws_mat=draws_mat)) +} + + +#' Process the results of Stan samples and GetStanSensitivityModel into a +#' sensitivity matrix. Note that currently, only the first chain is supported. +#' +#' @param sampling_result The output of \code{stan::sampling} +#' @param stan_sensitivity_list The output of \code{GetStanSensitivityModel} +#' @return A list of matrices. The elements of the list are +#' \itemize{ +#' \item{sens_mat: }{The local sensitivity of each posterior parameter to each hyperparameter.} +#' \item{sens_mat_normalized: }{\code{sens_mat}, but normalized by the posterior standard deviation.} +#' \item{grad_mat: }{The gradients of the log posterior evaluatated at the draws.} +#' \item{lp_vec: }{The log probability of the model at each draw.} +#' \item{draws_mat: }{The parameter draws in the same order as that of \code{grad_mat}}. +#' } +#' The result can be used directly, or passed to \code{GetTidyResult}. +#' @export +GetStanSensitivityFromModelFit <- function( + sampling_result, stan_sensitivity_list) { + + # TODO: deprecate this for something more general. + sens_mats <- GetStanSensitivityMatricesFromModelFit( + sampling_result, stan_sensitivity_list) + # Calculate the sensitivity. - sens_mat <- GetSensitivityFromGrads(grad_mat, draws_mat) + sens_mats$sens_mat <- GetSensitivityFromGrads( + sens_mats$grad_mat, sens_mats$draws_mat) # Normalize by the marginal standard deviation. - sens_mat_normalized <- NormalizeSensitivityMatrix(sens_mat, draws_mat) + sens_mats$sens_mat_normalized <- NormalizeSensitivityMatrix( + sens_mats$sens_mat, sens_mats$draws_mat) - return(list(sens_mat=sens_mat, sens_mat_normalized=sens_mat_normalized, - grad_mat=grad_mat, lp_vec=model_at_draws$lp_vec, - draws_mat=draws_mat)) + return(sens_mats) } @@ -389,30 +444,9 @@ GetImportanceSamplingFromModelFit <- function( imp_means <- colSums(imp_weights * draws_mat) - return(list(imp_weights=imp_weights, imp_lp_vec=imp_lp_vec, lp_vec=lp_vec, + return(list(imp_weights=imp_weights, + imp_lp_vec=imp_lp_vec, + lp_vec=lp_vec, eff_num_imp_samples=eff_num_imp_samples, imp_means=imp_means)) } - - -#' Get a dataframe with the hyperparameter values contained in a stan data list. -#' -#' @param stan_sensitivity_list The output of \code{GetStanSensitivityModel} -#' @param stan_data The Stan data for the original model. -#' @return A tidy dataframe with the hyperparameter names and values. -#' @export -GetHyperparameterDataFrame <- function(stan_sensitivity_list, stan_data) { - sens_par_list <- SetSensitivityParameterList( - stan_sensitivity_list$model_sens_params, - stan_sensitivity_list$model_params, stan_data) - pars_free <- unconstrain_pars( - stan_sensitivity_list$model_sens_fit, sens_par_list) - names(pars_free) <- stan_sensitivity_list$sens_param_names - hyperparam_names <- setdiff(stan_sensitivity_list$sens_param_names, - stan_sensitivity_list$param_names) - hyperparameter_df <- - data.frame(hyperparameter=hyperparam_names, - hyperparameter_val=pars_free[hyperparam_names]) - rownames(hyperparameter_df) <- NULL - return(hyperparameter_df) -} diff --git a/rstansensitivity/R/tidy_helpers_lib.R b/rstansensitivity/R/tidy_helpers_lib.R new file mode 100644 index 0000000..aec3f8a --- /dev/null +++ b/rstansensitivity/R/tidy_helpers_lib.R @@ -0,0 +1,333 @@ +# Functions to process stan and sensitivity results into tidy dataframes. +# +# TODO: move most of this over to tidybayes + +library(ggplot2) +library(dplyr) +library(reshape2) + + +#' Get a dataframe with the hyperparameter values contained in a stan data list. +#' +#' @param stan_sensitivity_list The output of \code{GetStanSensitivityModel} +#' @param stan_data The Stan data for the original model. +#' @return A tidy dataframe with the hyperparameter names and values. +#' @export +GetHyperparameterDataFrame <- function(stan_sensitivity_list, stan_data) { + sens_par_list <- SetSensitivityParameterList( + stan_sensitivity_list$model_sens_params, + stan_sensitivity_list$model_params, + stan_data) + pars_free <- rstan::unconstrain_pars( + stan_sensitivity_list$model_sens_fit, sens_par_list) + names(pars_free) <- stan_sensitivity_list$sens_param_names + hyperparam_names <- setdiff(stan_sensitivity_list$sens_param_names, + stan_sensitivity_list$param_names) + hyperparameter_df <- + data.frame(hyperparameter=hyperparam_names, + hyperparameter_val=pars_free[hyperparam_names]) + rownames(hyperparameter_df) <- NULL + return(hyperparameter_df) +} + + +#' Execute an R file and store its environment variables in a list. +#' @param data_file The filename of a code containing R code.. +#' @return A list with the variables defined when running `data_file`. +#' @export +LoadStanData <- function(data_file) { + stan_data <- new.env() + source(data_file, local=stan_data) + stan_data <- as.list(stan_data) + return(stan_data) +} + +SensitivityMatrixToDataframe <- function( + sens_mat, hyperparameter_names, parameter_names) { + colnames(sens_mat) <- parameter_names + return(data.frame(sens_mat) %>% + mutate(hyperparameter=hyperparameter_names) %>% + melt(id.var="hyperparameter") %>% + rename(parameter=variable)) +} + + +#' Make a tidy dataframe out of a sensitivity matrix and its standard errors. +#' @param sens_mat A matrix of sensitivities. Hyperparameters should be in +#' the rows and parameters in the columns. +#' @param sens_se A matrix of standard errors of \code{sens_mat}. +#' @param measure What to call these sensitivites. +#' @param num_se The number of standard errors for the upper and lower bounds. +#' @return A tidy dataframe with columns for the parameters, hyperparameters, +#' sensitivities, and their standard errors. +#' @export +SummarizeSensitivityMatrices <- function( + sens_mat, sens_se, measure, num_se=qnorm(0.975)) { + + sens_df <- rbind( + SensitivityMatrixToDataframe( + sens_mat, + hyperparameter_names=rownames(sens_mat), + parameter_names=colnames(sens_mat)) %>% + mutate(measure=measure), + SensitivityMatrixToDataframe( + sens_se, + hyperparameter_names=rownames(sens_se), + parameter_names=colnames(sens_se)) %>% + mutate(measure=paste(measure, "se", sep="_")), + SensitivityMatrixToDataframe( + sens_mat + num_se * sens_se, + hyperparameter_names=rownames(sens_mat), + parameter_names=colnames(sens_mat)) %>% + mutate(measure=paste(measure, "upper", sep="_")), + SensitivityMatrixToDataframe( + sens_mat - num_se * sens_se, + hyperparameter_names=rownames(sens_mat), + parameter_names=colnames(sens_mat)) %>% + mutate(measure=paste(measure, "lower", sep="_")) + ) %>% + dcast(hyperparameter + parameter ~ measure) %>% + filter(parameter != "lp__") + return(sens_df) +} + +#' Process the results of GetStanSensitivityFromModelFit into a +#' tidy dataframe with standard errors. +#' +#' @param sens_result The output of \code{GetStanSensitivityFromModelFit}. +#' @param num_se The number of standard errors for the upper and lower bounds. +#' @return A dataframe summarizing the sensitivity of the model posterior means +#' to the hyperparameters. The reported standard errors are based on a +#' multivariate normal and delta method approximation which may not be +#' accurate for highly variable or non-normal draws. The standard errors should +#' be interpreted with caution. +#' \itemize{ +#' \item{hyperparameter: }{The name of the hyperparameter.} +#' \item{parameter: }{The name of the parameter.} +#' \item{sensitivity:} +#' {The MCMC estimate of \code{dE[parameter | X, hyperparameter] / d hyperparameter}.} +#' \item{sensitivity_se: } +#' {The estimated Monte Carlo error of \code{sensitivity}.} +#' \item{normalized_sensitivity: } +#' {The MCMC estimate of \code{sensitivity / sd(parameter}).} +#' \item{normalized_sensitivity_se: } +#' {The estimated Monte Carlo error of \code{normalized_sensitivity}.} +#' } +#' @export +GetTidyResult <- function(sens_result, num_se=2) { + sens_se <- GetSensitivityStandardErrors( + sens_result$draws_mat, sens_result$grad_mat, + fix_mean=FALSE, normalized=FALSE) + norm_sens_se <- GetSensitivityStandardErrors( + sens_result$draws_mat, sens_result$grad_mat, + fix_mean=FALSE, normalized=TRUE) + sens_df <- SummarizeSensitivityMatrices( + sens_result$sens_mat, sens_se, + measure="sensitivity", num_se=num_se) + sens_norm_df <- SummarizeSensitivityMatrices( + sens_result$sens_mat_normalized, norm_sens_se, + measure="normalized_sensitivity", num_se=num_se) + + result <- inner_join( + sens_df, sens_norm_df, by=c("hyperparameter", "parameter")) + + return(result) +} + + +#' Make a barchart from the output of \code{GetTidyResult}. +#' +#' @param sens_df The output of \code{GetTidyResult}. +#' @param normalized Whether to display the normalized sensitivities. +#' @param se_num The number of standard errors to plot with the errorbars. +#' @return A \code{ggplot} object showing the sensitivites of each parameter +#' to each hyperparameter. +#' @export +PlotSensitivities <- function(sens_df, normalized=TRUE, se_num=2) { + if (normalized) { + y_axis_label <- "Normalized sensitivity" + sens_df_plot <- + mutate(sens_df, s=normalized_sensitivity, + s_se=normalized_sensitivity_se) + } else { + y_axis_label <- "Unnormalized sensitivity" + sens_df_plot <- + mutate(sens_df, s=sensitivity, s_se=sensitivity_se) + } + return( + ggplot(sens_df_plot) + + geom_bar(aes(x=parameter, y=s, fill=hyperparameter), + stat="identity", position="dodge") + + geom_errorbar(aes(x=parameter, + ymin=s - se_num * s_se, + ymax=s + se_num * s_se, + group=hyperparameter), + position=position_dodge(0.9), width=0.2) + + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + + scale_fill_discrete(name="Hyperparameter") + + ylab(y_axis_label) + xlab("Parameter") + ) +} + + +#' Convert a Stan sampling result to a form that can be joined with tidy +#' sensitivity results. +#' +#' @param sampling_result The output of \code{stan::sampling} +#' @param cols The columns of the summary to keep. +#' @return A data frame with the MCMC reuslts that can be joined with tidy +#' sensitivity results. +#' @export +GetMCMCDataFrame <- function( + sampling_result, cols=c("mean", "se_mean", "sd", "n_eff", "Rhat")) { + mcmc_result <- as.data.frame(rstan::summary(sampling_result)$summary) + mcmc_result$parameter <- make.names(rownames(mcmc_result)) + rownames(mcmc_result) <- NULL + mcmc_result <- + select(mcmc_result, "parameter", cols) %>% + filter(parameter != "lp__") + return(mcmc_result) +} + + +#' Use the linear approximation to predict the sensitivity to a new +#' stan data list. +#' +#' @param stan_sensitivity_list The output of \code{GetStanSensitivityModel} +#' @param stan_result The output of \code{GetStanSensitivityFromModelFit} +#' @param stan_data The original stan data at which +#' \code{stan_sensitivity_list} was calculated. +#' @param stan_data_perturb A new stan data file with different hyperparameters. +#' @param description A hyperparameter name to describe this perturbation. +#' @return A tidy sensitivity dataframe where the sensitivity is in the +#' direction of the difference between the hyperparameters in the two stan +#' data lists. +#' @export +PredictSensitivityFromStanData <- function( + stan_sensitivity_list, sens_result, stan_data, stan_data_perturb, + description="perturbation") { + + hyperparameter_df <- + inner_join( + GetHyperparameterDataFrame(stan_sensitivity_list, stan_data) %>% + rename(hyperparameter_val_orig=hyperparameter_val), + GetHyperparameterDataFrame(stan_sensitivity_list, stan_data_perturb), + by="hyperparameter") %>% + mutate(hyperparameter_diff=hyperparameter_val - hyperparameter_val_orig) + linear_comb <- matrix(hyperparameter_df$hyperparameter_diff, nrow=1) + colnames(linear_comb) <- hyperparameter_df$hyperparameter + linear_comb <- linear_comb[, rownames(sens_result$sens_mat), drop=FALSE] + rownames(linear_comb) <- description + sens_result_pert <- TransformSensitivityResult(sens_result, linear_comb) + return(GetTidyResult(sens_result_pert)) +} + + +# Get a tidy version of a Stan summary with tidybayes. +# +#' @param stanfit A \code{stanfit} object (e.g., from \code{sampling}). +#' @param ... Parameter names in the style of \code{tidybayes}. +#' @param spread Optional. If \code{TRUE}, return a wide dataframe in which +#' each summary metric is a column. Otherwise, return a tall dataframe +#' in which the summary metric is in a column called \code{metric}. +#' @return A tidy dataframe containing the \code{stan} summary object. +#' @export +GetTidyStanSummary <- function(stanfit, ..., spread=FALSE) { + pars <- enquos(...) + summary_mat <- t(rstan::summary(stanfit)$summary) + tidy_summary <- + gather_draws(summary_mat, !!!pars) %>% + inner_join(data.frame(.draw=1:nrow(summary_mat), + metric=rownames(summary_mat)), + by=".draw") %>% + select(-.chain, -.draw, -.iteration) + + if (spread) { + key_string <- paste(setdiff(names(tidy_summary), + c("metric", ".value")), collapse=" + ") + tidy_summary <- dcast( + tidy_summary, + formula(sprintf("%s ~ metric", key_string)), + value.var=".value") + } + return(tidy_summary) +} + + +RemoveExtraTidyColumns <- function(df) { + select(df, -.chain, -.iteration, -.draw) +} + + +#' Convert a column of variable names to a tidy format. +#' @param df A dataframe with a parameter column to be tidied. +#' @param col The string name of the column to be tidied. +#' @param ... A parameter listing in \code{tidybayes} for the parameters in the +#' column \code{df[[col]]}. +#' @param variable_name Optional. A string for the new column's variable name. +#' The default is ".variable", as is standard for \code{tidybayes}. +#' @return A new dataframe with extra \code{tidybayes}-style columns. +#' @export +TidyColumn <- function(df, col, ..., variable_name=".variable") { + pars <- enquos(...) + par_names <- unique(df[[col]]) + tidy_pars_df <- + matrix(par_names, nrow=1, dimnames=list("", par_names)) %>% + tidybayes::gather_draws(!!!pars) %>% + rename(!!sym(col):=.value, !!sym(variable_name):=.variable) %>% + select(-.chain, -.iteration, -.draw) + return(inner_join(df, tidy_pars_df, by=col)) +} + + +#' Gather a matrix as with \code{gather_draws}, adding a column for the +#' rownames. +GatherRowNamedMatrix <- function(mat, ...) { + pars <- enquos(...) + df <- + tidybayes::gather_draws(mat, !!!pars) %>% + inner_join(data.frame(.draw=1:nrow(mat), .rowname=rownames(mat)), + by=".draw") %>% + RemoveExtraTidyColumns() + return(df) +} + + +#' @export +GetTidySensitivityResults <- function(grad_mat, + draws_mat, + ..., + normalize=TRUE, + calculate_se=TRUE) { + pars <- enquos(...) + GetSensitivityDf <- function(mat, measure_name) { + GatherRowNamedMatrix(mat, !!!pars) %>% + mutate(measure=measure_name) %>% + rename(hyperparameter=.rowname) + } + sens_mat <- GetSensitivityFromGrads(grad_mat, draws_mat) + results <- GetSensitivityDf(sens_mat, "sensitivity") + if (calculate_se) { + sens_mat_se <- GetSensitivityStandardErrors( + draws_mat, grad_mat, normalized=FALSE) + results <- bind_rows( + results, + GetSensitivityDf(sens_mat_se, "sensitivity_se")) + } + + if (normalize) { + # Normalize by the marginal standard deviation. + sens_mat_norm <- NormalizeSensitivityMatrix(sens_mat, draws_mat) + results <- bind_rows( + results, + GetSensitivityDf(sens_mat_norm, "normalized_sensitivity")) + if (calculate_se) { + sens_mat_norm_se <- GetSensitivityStandardErrors( + draws_mat, grad_mat, normalized=TRUE) + results <- bind_rows( + results, + GetSensitivityDf(sens_mat_norm_se, "normalized_sensitivity_se")) + } + } + return(results) +} diff --git a/rstansensitivity/docs/404.html b/rstansensitivity/docs/404.html new file mode 100644 index 0000000..56645c3 --- /dev/null +++ b/rstansensitivity/docs/404.html @@ -0,0 +1,127 @@ + + + + + + + + +Page not found (404) • rstansensitivity + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ + + + +
+ +
+
+
+

Page not found (404)

+
+ +Content not found. Please use links in the navbar. + +
+ +
+ + + + +
+ + + + + + + + diff --git a/rstansensitivity/docs/authors.html b/rstansensitivity/docs/authors.html index 76e050a..ce9a0d1 100644 --- a/rstansensitivity/docs/authors.html +++ b/rstansensitivity/docs/authors.html @@ -1,6 +1,6 @@ - + @@ -9,20 +9,22 @@ Authors • rstansensitivity - + - - + + - + + - + - - + + + @@ -35,7 +37,8 @@ - + + + @@ -51,14 +55,15 @@