diff --git a/.gitignore b/.gitignore index 5b6a065..3fb4054 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,6 @@ .Rhistory .RData .Ruserdata + +# local testing files +tempfiles/ diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION index 9d146bf..1a4f00c 100644 --- a/CRAN-SUBMISSION +++ b/CRAN-SUBMISSION @@ -1,3 +1,3 @@ Version: 0.1.0 -Date: 2024-02-20 15:36:22 UTC -SHA: 40df86c9c02e76044b2499df231d4ddc43ce4940 +Date: 2024-02-27 07:01:18 UTC +SHA: 15721ee1ca7a92c51b130bafde6eb79d19e235da diff --git a/DESCRIPTION b/DESCRIPTION index ca102f2..7455f14 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: tsnet Type: Package Title: Fitting, Comparing, and Visualizing Networks Based on Time Series Data -Version: 0.1.0 +Version: 0.2.0 Authors@R: c(person("Björn S.", "Siepe", email = "bjoernsiepe@gmail.com", @@ -22,13 +22,14 @@ Encoding: UTF-8 LazyData: true URL: https://github.com/bsiepe/tsnet BugReports: https://github.com/bsiepe/tsnet/issues -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 Imports: cowplot, dplyr, ggdist, ggokabeito, ggplot2, + loo, methods, posterior, Rcpp (>= 0.12.0), diff --git a/NAMESPACE b/NAMESPACE index e9da964..cbe5956 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -26,6 +26,12 @@ importFrom(ggdist,stat_slab) importFrom(ggdist,theme_ggdist) importFrom(ggokabeito,palette_okabe_ito) importFrom(ggokabeito,scale_fill_okabe_ito) +importFrom(loo,loo) +importFrom(loo,loo_model_weights) +importFrom(loo,pareto_k_values) +importFrom(loo,psis) +importFrom(loo,relative_eff) +importFrom(loo,weights.importance_sampling) importFrom(posterior,as_draws_matrix) importFrom(rlang,.data) importFrom(rstan,extract) diff --git a/NEWS.md b/NEWS.md index dd2cb36..21c458f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,6 @@ +# tsnet 0.2.0 +* Fixed documentation issues (e.g., wrong display of LaTex equations) + # tsnet 0.1.0 * Initial CRAN submission. diff --git a/R/centrality.R b/R/centrality.R index 2a4bc92..bc10727 100644 --- a/R/centrality.R +++ b/R/centrality.R @@ -5,17 +5,18 @@ #' a network, while density describes the networks' overall connectedness. #' Specifically, it computes the in-strength, out-strength, contemporaneous #' strength, temporal network density, and contemporaneous network density. The -#' result can then be visualized using [plot_centrality()]. +#' result can then be visualized using \code{\link{plot_centrality}}. #' #' @param fitobj Fitted model object for a Bayesian GVAR model. This can be -#' `tsnet_fit` object (obtained from [stan_gvar()]), a BGGM object (obtained -#' from [BGGM::var_estimate()]), or extracted posterior samples (obtained from -#' [stan_fit_convert()). +#' `tsnet_fit` object (obtained from \code{\link{stan_gvar}}, +#' a BGGM object (obtained from \code{\link[BGGM]{var_estimate}} in \code{BGGM}), +#' or extracted posterior samples (obtained from \code{\link{stan_fit_convert}}). #' @param burnin An integer specifying the number of initial samples to discard -#' as burn-in. Default is 0. +#' as burn-in. Default is \code{0}. #' @param remove_ar A logical value specifying whether to remove the -#' autoregressive effects for centrality calculation. Default is TRUE. This is -#' only relevant for the calculation of temporal centrality/density measures. +#' autoregressive effects for centrality calculation. Default is \code{TRUE}. +#' This is only relevant for the calculation of temporal centrality/density +#' measures. #' #' @return A list containing the following centrality measures: #' \itemize{ @@ -106,11 +107,11 @@ get_centrality <- function(fitobj, #' centrality measures. #' #' @param obj An object containing the centrality measures obtained from -#' [get_centrality()]. +#' \code{\link{get_centrality}}. #' @param plot_type A character string specifying the type of plot. Accepts -#' "tiefighter" or "density". Default is "tiefighter". +#' "tiefighter" or "density". Default is \code{"tiefighter"}. #' @param cis A numeric value specifying the credible interval. Must be between -#' 0 and 1 (exclusive). Default is 0.95. +#' 0 and 1 (exclusive). Default is \code{0.95}. #' #' @return A ggplot object visualizing the centrality measures. For a #' "tiefighter" plot, each point represents the mean centrality measure for a diff --git a/R/check_eigen.R b/R/check_eigen.R index ef80385..5ffa20f 100644 --- a/R/check_eigen.R +++ b/R/check_eigen.R @@ -2,17 +2,17 @@ #' #' This function checks the eigenvalues of the Beta matrix (containing the #' temporal coefficients) to assure that the model is stationary. It uses the -#' same check as the `graphicalVAR` package. The function calculates the +#' same check as the \code{graphicalVAR} package. The function calculates the #' eigenvalues of the Beta matrix and checks if the sum of the squares of the #' real and imaginary parts of the eigenvalues is less than 1. If it is, the VAR #' model is considered stable. #' #' @param fitobj A fitted Bayesian GVAR object. This can be a tsnet_fit object -#' (obtained from [stan_gvar()]), a BGGM object (obtained from -#' [BGGM::var_estimate()]), or extracted posterior samples (obtained from -#' [stan_fit_convert()). +#' (obtained from \code{\link{stan_gvar}}, a \code{BGGM} object (obtained from +#' \code{\link[BGGM]{var_estimate}}), or extracted posterior samples +#' (obtained from \code{\link{stan_fit_convert}}. #' @param verbose Logical. If TRUE, a verbal summary of the results is printed. -#' Default is TRUE. +#' Default is \code{TRUE}. #' @examples #' data(fit_data) #' fitobj <- fit_data[[1]] diff --git a/R/compare_gvar.R b/R/compare_gvar.R index 2ebdde1..61a906b 100644 --- a/R/compare_gvar.R +++ b/R/compare_gvar.R @@ -10,39 +10,39 @@ #' . #' #' @param fit_a Fitted model object for Model A. This can be a tsnet_fit object -#' (obtained from [stan_gvar()]), a BGGM object (obtained from -#' [BGGM::var_estimate()]), or extracted posterior samples (obtained from -#' [stan_fit_convert()). +#' (obtained from \code{\link{stan_gvar}}, a \code{BGGM} object (obtained from +#' \code{\link[BGGM]{var_estimate}}, or extracted posterior samples (obtained from +#' \code{\link{stan_fit_convert}}). #' @param fit_b Fitted model object for Model B. This can be a tsnet_fit object -#' (obtained from [stan_gvar()]), a BGGM object (obtained from -#' [BGGM::var_estimate()]), or extracted posterior samples (obtained from -#' [stan_fit_convert()). -#' @param cutoff The percentage level of the test (default: 5\%) as integer. -#' @param dec_rule The decision rule to be used. Currently supports default "or" +#' (obtained from \code{\link{stan_gvar}}, a \code{BGGM} object (obtained from +#' \code{\link[BGGM]{var_estimate}}, or extracted posterior samples (obtained from +#' \code{\link{stan_fit_convert}}). +#' @param cutoff The percentage level of the test (default: \code{5}) as integer. +#' @param dec_rule The decision rule to be used. Currently supports default \code{"or"} #' (comparing against two reference distributions) and "comb" (combining the #' reference distributions). The use of "or" is recommended, as "comb" is less #' stable. #' @param n_draws The number of draws to use for reference distributions -#' (default: 1000). +#' (default: \code{1000}). #' @param comp The distance metric to use. Should be one of "frob" (Frobenius #' norm), "maxdiff" (maximum difference), or "l1" (L1 norm) (default: -#' "frob"). The use of the Frobenius norm is recommended. +#' \code{"frob"}). The use of the Frobenius norm is recommended. #' @param return_all Logical indicating whether to return all distributions -#' (default: FALSE). Has to be set to TRUE for plotting the results. +#' (default: \code{FALSE}). Has to be set to TRUE for plotting the results. #' @param sampling_method Draw sequential pairs of samples from the posterior, #' with certain distance between them ("sequential") or randomly from two #' halves of the posterior ("random"). The "random" method is preferred to #' account for potential autocorrelation between subsequent samples. Default: -#' "random". +#' \code{"random"}. #' @param indices A list of "beta" and "pcor" indices specifying which elements -#' of the matrices to consider when calculating distances. If NULL (default), +#' of the matrices to consider when calculating distances. If \code{NULL} (default), #' all elements of both matrices are considered. If provided, only the #' elements at these indices are considered. If only one of the matrices #' should have indices, the other one should be NULL. This can be useful if #' you want to calculate distances based on a subset of the elements in the #' matrices. -#' @param burnin The number of burn-in iterations to discard (default: 0). -#' @return A list (of class "compare_gvar") containing the results of the +#' @param burnin The number of burn-in iterations to discard (default: \code{0}). +#' @return A list (of class \code{compare_gvar}) containing the results of the #' comparison. The list includes: #' #' \item{sig_beta}{Binary decision on whether there is a significant difference between the temporal networks of A and B} diff --git a/R/data.R b/R/data.R index 35394a1..599b5e1 100644 --- a/R/data.R +++ b/R/data.R @@ -1,7 +1,7 @@ #' Simulated Time Series Dataset #' #' This dataset contains a simulated time series dataset for two individuals -#' generated using the `graphicalVAR` package. The dataset is useful for testing +#' generated using the \code{graphicalVAR} package. The dataset is useful for testing #' and demonstrating the functionality of the package. #' #' @format ## `ts_data` A data frame with 500 rows and 7 columns. @@ -9,7 +9,7 @@ #' \item{id}{A character string identifier for the individual. There are two unique ids, representing two individuals.} #' \item{V1-V6}{These columns represent six different variables in the time series data.} #' } -#' @source Simulated using the [graphicalVAR::graphicalVARsim()] function. +#' @source Simulated using the \code{\link[graphicalVAR]{graphicalVARsim}} function in the \code{graphicalVAR} package. #' #' @usage data(ts_data) #' @@ -28,19 +28,19 @@ #' Example Posterior Samples #' #' This dataset contains posterior samples of beta coefficients and partial correlations for two individuals. -#' It was generated by fitting a GVAR model using [stan_gvar()] with three variables from the [ts_data] dataset. +#' It was generated by fitting a GVAR model using \code{\link{stan_gvar}} with three variables from the \code{\link{ts_data}} dataset. #' #' @format ## `fit_data` #' A list with two elements, each containing posterior samples for one individual. #' #' @usage data(fit_data) #' -#' @source The data is generated using the [stan_gvar()] function on subsets -#' of the [ts_data] time series data. +#' @source The data is generated using the \code{\link{stan_gvar}} function on subsets +#' of the \code{\link{ts_data}} time series data. #' #' @details #' The list contains two elements, each containing posterior samples for one individual. -#' The samples were extracted using the [stan_fit_convert()] function. +#' The samples were extracted using the \code{\link{stan_fit_convert}} function. #' For each individual, the list elements contain the posterior means of the beta coefficients #' ("beta_mu") and the posterior means of the partial correlations ("pcor_mu"). #' The "fit" element contains all 1000 posterior samples of the beta coefficients and partial correlations. diff --git a/R/forecast.stan_gvar.R b/R/forecast.stan_gvar.R new file mode 100644 index 0000000..9913371 --- /dev/null +++ b/R/forecast.stan_gvar.R @@ -0,0 +1,323 @@ +#' Forecast Using Fitted GVAR Models +#' +#' @description +#' This function generates forecasts using one or more fitted GVAR models. +#' It supports forecasting with model combinations, optionally weighted by model performance, +#' and handles scenarios with or without new data for prediction. +#' +#' @param fitobj A \code{tsnet_fit} object or a \code{tsnet_list} object containing one or more fitted GVAR models. +#' @param ahead An integer specifying the forecast horizon. Default is \code{1}. +#' @param new_data Optional. A matrix or array of new data for forecasting. If not provided, forecasts are generated for an empty array. +#' @param interval A numeric vector specifying the prediction interval (e.g., \code{c(0.025, 0.975)}). Values must lie between 0 and 1. +#' @param weights An object of class \code{model_weights}, typically obtained from \code{\link{model_weights}}. +#' If not provided, the function computes model weights internally. +#' @param start Integer indicating the starting point for cross-validation. Inherited from \code{\link{lfo.stan_gvar}}. Default is \code{1}. +#' @param method A character string specifying the method for model combination. Default is \code{"stacking"}. +#' @param inc_samples Logical. If \code{TRUE}, includes all posterior samples in the output. Default is \code{FALSE}. +#' @param ... Additional arguments (currently unused). +#' +#' @details +#' The function processes the input models (either a single \code{tsnet_fit} or a list of them), +#' validates that all models are based on the same time series, and prepares data for forecasting. +#' If no model weights are provided, the function computes them using \code{\link{model_weights}}. +#' Forecasts are computed for the specified horizon (\code{ahead}) and optionally weighted by model performance. +#' +#' The forecasts include posterior predictive intervals based on the specified \code{interval} argument. +#' If \code{new_data} is not provided, forecasts assume no additional data and generate predictions accordingly. +#' +#' @return A \code{tsnet_forecast} object with the following components: +#' \describe{ +#' \item{forecast}{A list containing the forecasted means and optionally the full posterior samples.} +#' \item{mean}{The mean forecast values.} +#' \item{samples}{(Optional) A list of forecasted posterior samples if \code{inc_samples} is \code{TRUE}.} +#' } +#' +#' +#' @examples +#' \dontrun{ +#' # Assuming `fit1` and `fit2` are tsnet_fit objects +#' fits <- tsnet_list(fit1, fit2) +#' forecasts <- forecast.stan_gvar(fits, ahead = 5, method = "stacking") +#' print(forecasts) +#' } +#' @keywords internal +#' @noRd + + +forecast.stan_gvar <- function(fitobj, + ahead = 1, + new_data = NULL, + interval = c(0.025, 0.975), + weights = NULL, + start = 1, + method = "stacking", + inc_samples = FALSE, + ...) { + + ## Input checks + # Input can either be a single model or a list of models + # check if the input is a tsnet_fit or a tsnet_list + if (!inherits(fitobj, "tsnet_fit") & + !inherits(fitobj, "tsnet_list")) { + stop("Input must be a tsnet_fit or a tsnet_list object") + } + + # Check if the interval is between 0 and 1 + if (any(interval < 0) | any(interval > 1)) { + stop("Interval must be between 0 and 1") + } + # Check if weights are specified correctly + if(!is.null(weights) && !inherits(weights, "model_weights")){ + stop("Model weights must be obtained from the model_weights function.") + } + + + ## Prepare data + # number of models + n_mods <- 1 + if (inherits(fitobj, "tsnet_list")) { + n_mods <- length(fitobj) + } + else { + fitobj <- tsnet_list(fitobj) + } + + # check if they all used the same time series + col_names <- fitobj[[1]]$arguments$cnames + colnames_ident <- all(sapply(fitobj, function(x){ + all(x$names == col_names) + })) + if(!colnames_ident){ + stop("All models must use the same time series.") + } + + # newdata available? + if (is.null(new_data)) { + newdata <- array(0, dim = c(ahead, fitobj[[1]]$arguments$n_t)) + compute_log_lik <- 0 + } else { + compute_log_lik <- 1 + } + + # model weights available? + # Case 1: No model weights provided + # Case 2: model weights from a model_weight object + # Case 3: No model weights requested + if(n_mods > 1 & is.null( weights ) ) { + mw <- model_weights(fitobjs = fitobj, start = start) + weights <- mw$weights [] + } else if( n_mods > 1 & !is.null( weights ) ) { + weights <- weights$weights[] + } else if( n_mods == 1 ) { + weights <- 1 + } + + + ## Compute Forecasts + # Looping over all specified models + forecast_obj <- lapply(fitobj, function(m){ + + args <- m$arguments$fn_args + prior_delta <- (1 / m$arguments$priors$prior_Rho_marginal) - 1 + + standat <- list(Y = as.matrix(m$arguments$data), + K = m$arguments$p, + "T" = m$arguments$n_t, + beep = eval(args$beep), + prior_Beta_loc = m$arguments$priors$prior_Beta_loc, + prior_Beta_scale = m$arguments$priors$prior_Beta_scale, + prior_Rho_loc = m$arguments$priors$prior_Rho_loc, + prior_Rho_scale = m$arguments$priors$prior_Rho_scale, + prior_Eta = m$arguments$priors$prior_Eta, + prior_S = m$arguments$priors$prior_S, + prior_delta = prior_delta, + ahead = ahead, + Y_future = new_data, + compute_log_lik = compute_log_lik) + + # Obtain model name for sampling + gqs_model <- .get_gvar_name(m) + + # backcast <- max(m$mgarchP, m$mgarchQ) + # this is just lag-1 here + backcast <- 1 + + # TODO think really good about the indexing here + nt <- m$arguments$n_t + cast_start <- (nt - backcast + 1) + forecast_start <- (nt + 1) + forecast_end <- (nt + ahead) + + forecasted <- rstan::gqs(stanmodels[[gqs_model]], + draws = as.matrix(m$stan_fit), + data = standat) + return(forecasted) + }) + + ## Extract forecasts + # extract forecasted values with weighting + # f_mean <- + # restructure to correct format + + + + + ## Output + # Return a tsnet_forecast object + out <- list() + out$forecast <- list() + out$forecast$forecast_obj <- forecast_obj + # out$forecast$weighted_forecasts <- .weighted_samples(forecast_obj, + # params = "Y_forecast", + # weights = weights) + + + # should all samples be included? + # if(inc_samples){ + # f_samples <- .weighted_samples(forecast_obj, weights = weights, method = method) + # out$forecast$samples$mean <- f_samples + # } + if(compute_log_lik == 1 ) { + log_lik <- lapply(forecast_obj, function(x) { + rstan::extract(x, pars = "log_lik")$log_lik + }) + out$forecast$log_lik <- log_lik + } + + + # Store arguments + mc <- match.call() + fl <- formals() + missing_args <- setdiff(names(fl), names(mc)) + missing_args <- Map(function(arg, default) + if (!is.null(default)) mc[[arg]] <- default, missing_args, fl[missing_args]) + all_args <- c(as.list(mc), missing_args) + + args <- list( + ahead = ahead, + new_data = new_data, + interval = interval, + weights = weights, + start = start, + method = method, + inc_samples = inc_samples, + fn_args = all_args + ) + + out$arguments <- args + + class(out) <- c("tsnet_forecast", class(out)) + return(out) + +} + + + + +# Helpers + +#' Create a List of \code{tsnet_fit} Objects +#' +#' @description +#' This internal function creates a list of \code{tsnet_fit} objects, which are the +#' results of fitting a time series network model using the\code{\link{stan_gvar}} function. +#' The list is assigned the class \code{tsnet_list}, making it compatible with other +#' functions that expect such a list. +#' +#' @param ... Objects of class \code{tsnet_fit}, which are the results of fitting +#' time series network models. These objects are passed as arguments to the function. +#' +#' @details +#' The function takes any number of \code{tsnet_fit} objects (from the \code{\link{stan_gvar}} +#' model fitting process), combines them into a list, and assigns the class \code{tsnet_list}. +#' This allows the list to be used as input in other functions that expect a \code{tsnet_list} +#' object, providing compatibility with workflows that involve multiple time series models. +#' +#' @return A list of \code{tsnet_fit} objects, with the class \code{tsnet_list}. +#' +#' @seealso \code{\link{stan_gvar}}, \code{\link{forecast.stan_gvar}}, \code{\link{model_weights}} +#' +#' @keywords internal +#' @noRd + +tsnet_list <- function(...) { + out <- list(...) + class(out) <- c("tsnet_list", class(out)) + return(out) +} + +# Obtain model name from stan_gvar output + +#' Get gvar name +#' +#' Obtain the name of the \code{\link{stan_gvar}} model from a fitte model output. +#' +#' @param fitobj A \code{\link{stan_gvar}} object +#' +#' @return character name of the \code{\link{stan_gvar}} model +#' +#' @keywords internal +#' @noRd +.get_gvar_name <- function(fitobj){ + if(fitobj$arguments$fn_args$cov_prior == "IW"){ + if(isTRUE(fitobj$arguments$fn_args$rmv_overnight)){ + gqs_model <- "VAR_wishart_beep" + } + else{ + gqs_model <- "VAR_wishart" + } + } + if(fitobj$arguments$fn_args$cov_prior == "LKJ"){ + if(isTRUE(fitobj$arguments$fn_args$rmv_overnight)){ + gqs_model <- "VAR_LKJ_beep" + } + else{ + gqs_model <- "VAR_LKJ" + } + } + return(gqs_model) +} + + + + +#' +#' @description +#' This internal function applies weights to the posterior samples of different models +#' and sums them together. The function is used to combine the posterior samples +#' from multiple models based on their weights. +#' +#' @param fitobj A list of fitted model objects from which to extract posterior samples. +#' @param params A character vector specifying the parameters to extract from the models. +#' @param weights A numeric vector of weights to apply to the posterior samples of each model. +#' +#' @details +#' The function first extracts the posterior samples for the specified parameters from each model. +#' It then applies the weights to the samples and sums the weighted samples together. +#' The resulting combined samples have the same dimensionality as the original samples. +#' +#' @return A list of combined posterior samples for each specified parameter. +#' +#' @examples +#' \dontrun{ +#' # Assuming `fit1` and `fit2` are fitted model objects +#' fitobj <- tsnet_list(fit1, fit2) +#' params <- c("param1", "param2") +#' weights <- c(0.6, 0.4) +#' combined_samples <- .weighted_samples(fitobj, params, weights) +#' } +#' @keywords internal +#' @noRd +.weighted_samples <- function(fitobj, params, weights) { + samps <- lapply(fitobj, rstan::extract, pars = params) + for(i in seq_len(length(samps))) { # each model + samps[[i]] <- lapply(samps[[i]], function(p) { # each parameter + p * weights[i] + }) + } + samps_comb <- lapply(params, function(p) { # each parameter + Reduce("+", lapply(samps, function(m) {m[[p]]})) + }) + names(samps_comb) <- params + return(samps_comb) +} diff --git a/R/lfo.R b/R/lfo.R new file mode 100644 index 0000000..c4539e4 --- /dev/null +++ b/R/lfo.R @@ -0,0 +1,269 @@ +#' Leave-Future-Out Cross-Validation for Stan GVAR Models +#' +#' @description +#' Computes the Expected Log Predictive Density (ELPD) using Leave-One-Out (LOO) or Leave-Future-Out (LFO) approaches for models fitted using \code{\link{stan_gvar}}. +#' Currently, only the LOO method is implemented. This function is designed to work with objects of class \code{tsnet_fit}. +#' +#' @param x A \code{tsnet_fit} object obtained from \code{\link{stan_gvar}}. +#' @param type Character string specifying the type of cross-validation to perform. Options are: +#' \itemize{ +#' \item \code{"loo"}: Leave-One-Out Cross-Validation. +#' \item \code{"lfo"}: Leave-Future-Out Cross-Validation (not yet fully implemented). +#' } +#' Default is \code{"loo"}. +#' @param start Integer indicating the starting point for cross-validation. Default is 0. +#' @param ahead Integer indicating the number of steps ahead for predictions in Leave-Future-Out mode. Default is 1. +#' @param mode Character string specifying the cross-validation mode when \code{type = "lfo"}. Options are: +#' \itemize{ +#' \item \code{"backward"}: LFO with backward refitting. +#' \item \code{"forward"}: LFO with forward refitting (not yet implemented). +#' \item \code{"exact"}: Exact LFO refitting at each time point (not yet implemented). +#' } +#' Default is \code{"backward"}. +#' @param k_thres Numeric threshold for Pareto k values. Only used in \code{lfo} mode to determine when to refit the model. Default is 0.7. +#' @param ... Additional arguments passed to underlying methods. +#' +#' @details +#' This function evaluates the predictive performance of \code{stan_gvar} models using cross-validation. +#' When \code{type = "loo"}, the function performs Leave-One-Out Cross-Validation using efficient importance sampling via the \pkg{loo} package. +#' For \code{type = "lfo"}, Leave-Future-Out Cross-Validation is conceptually supported but not yet fully implemented. +#' +#' @importFrom loo relative_eff loo psis pareto_k_values weights.importance_sampling +#' +#' @return A list with the following components: +#' \describe{ +#' \item{elpd}{Expected Log Predictive Density.} +#' \item{mode}{Mode of cross-validation, as specified in the \code{mode} parameter.} +#' } +#' +#' @note +#' The LFO method is a work in progress and may not yet support all modes or configurations. The \code{backward} mode partially implements the LFO algorithm but does not include forecasting functionality. +#' This function borrows heavily from \code{\link[bmgarch:loo.bmgarch]{bmgarch::loo.bmgarch}}. +#' +#' @keywords internal +#' @noRd +#' + +lfo.stan_gvar <- function(x, + type = "loo", + start = 0, + ahead = 1, + mode = "backward", + k_thres = 0.7, + ...) { + obj <- x + + # Check if x is a tsnet_fit object + if (!inherits(obj, "tsnet_fit")) { + stop("x must be a tsnet_fit object.") + } + + # Check if type is valid + if (!type %in% c("loo", "lfo")) { + stop("type must be one of 'loo' or 'lfo'.") + } + + # Check if mode is valid + if (!mode %in% c("backward", "forward", "exact")) { + stop("mode must be one of 'backward', 'forward', or 'exact'.") + } + + + + # Setup output list + out <- list() + out$type <- type + out$mode <- mode + + # Obtain the number of observations in log likelihood (one less than full obs) + # TODO should I change this? + n <- obj$stan_fit@par_dims$log_lik + 1 + + if (type == "loo") { + + # Obtain log-likelihood + ll_full <- .log_lik_tsnet(obj) + ll <- ll_full[, (start + 1):n] + + ## obtain chain id vector for relative_eff + n_chains <- obj$stan_fit@sim$chains + n_samples <- obj$stan_fit@sim$iter - obj$stan_fit@sim$warmup + chain_id <- rep(seq_len(n_chains), each = n_samples) + r_eff <- loo::relative_eff(exp(ll), chain_id = chain_id) + backcast_loo <- loo::loo(ll, r_eff = r_eff) + out$elpd <- backcast_loo$estimates[, 'Estimate']['elpd_loo'] + } + + else if (type == "lfo") { + lll <- matrix(nrow = dim(.log_lik_tsnet(obj))[1], ncol = n) + approx_elpds_1sap <- rep(NA, n) + ks <- NULL + refits <- NULL + + if (mode == "backward") { + # Start with full model, go backwards until pareto k threshold exceeded + # then refit + fit_past <- x + i_refit <- n + + # loop across possible backfits + for (i in seq((n - ahead), start, by = -ahead)) { + # obtain log-likelihood + ll[, (i + 1):(i + ahead)] <- .log_lik_tsnet(fit_past)[, (i + 1):(i + ahead)] + # compute log of raw importance ratios + # TODO check this + logratio <- .sum_log_ratios(ll, (i + 1):i_refit) + # compute pareto-smoothed importance sampling + psis_obj <- suppressWarnings(loo::psis(logratio)) + # obtain k-vals + k <- loo::pareto_k_values(psis_obj) + ks <- c(ks, k) + + # check if k-threshold is exceeded + if (k > k_thres) { + # refit the model + i_refit <- i + refits <- c(refits, i) + past <- 1:i + oos <- (i + 1):(i + ahead) + + # contains data for the past + # TODO need to check data structure (array/matrix) + df_past <- obj$args$data[past, , drop = FALSE] + # contains data for the out-of-sample + # TODO why both combined? + df_oos <- obj$args$data[c(past, oos), , drop = FALSE] + # refitted model + fit_past <- .refit(obj, + data = df_past) + + # then forecast + # fc <- forecast(...) + # ll[, (i+1):(i+ahead) ] <- fc$forecast$ll[[1]] + # if(ahead == 1 ) { + # approx_elpds_1sap[ i+1 ] <- .log_mean_exp(ll, i+1 ]) + # } else { + # approx_elpds_1sap[(i+1):(i+ahead)] <- + # apply(ll[, (i+1):(i+ahead) ], MARGIN = 2, FUN = .log_mean_exp ) + } + # if k threshold not exceeded + else { + lw <- loo:weights.importance_sampling(psis_obj, normalize = TRUE)[, 1] + if (ahead == 1) { + approx_elpds_1sap[i + 1] <- .log_sum_exp(lw + ll[, i + 1]) + } else { + approx_elpds_1sap[(i + 1):(i + ahead)] <- + apply((lw + ll[, (i + 1):(i + ahead)]), 2, .log_sum_exp) + } + } + } # end for loop + out$elpd <- approx_elpds_1sap + + } + else if (mode == "forward") { + stop("Forward LFO not implemented yet.") + } + else if (mode == "exact") { + stop("Exact LFO not implemented yet.") + # refit at each time point + k_thres <- 0 + refits <- n - (start+1) + + # todo obtain data + df_dat <- obj$args$data + exact_elpds_1sap <- rep(NA, n) + + # for(i in start:( n-1 ) ) { + # fit_start <- .refit(obj, + # data = df_dat[1:i, , drop = FALSE]) + # fc <- tsnet::forecast(fit_start, + # ahead = ahead, + # newdata = df_dat[ i+1, ,drop = FALSE]) + # loglik[, i+1 ] <- fc$forecast$log_lik[[1]] + # } + # loglik_exact <- loglik[, (start+1):n] + # exact_elpds_1sap <- apply(loglik_exact, 2, .log_mean_exp ) + # out <- exact_elpds_1sap + + } + + } + + + class(out) <- c("lfo.stan_gvar", class(out)) + return(out) +} + + + +# Function for refitting a model with given specification +# needed in lfo function above +#' DOES NOT WORK YET WITH CUSTOM PRIORS +#' @keywords internal +#' @noRd +.refit <- function(x, data) { + + args <- x$arguments$fn_args + + fit <- stan_gvar(data, + beep = args$beep, + priors = x$args$priors, + method = args$method, + cov_prior = args$cov_prior, + rmv_overnight = args$rmv_overnight, + iter_sampling = args$iter_sampling, + iter_warmup = args$iter_warmup, + n_chains = args$n_chains, + n_cores = args$n_cores, + center_only = args$center_only, + ahead = args$ahead, + compute_log_lik = args$compute_log_lik) + return(fit) +} + + + + +# Following functions are obtained from bmgarch/lfo-cv tutorial +# Helper to obtain log-likelihood +#' @keywords internal +#' @noRd + +.log_lik_tsnet <- function(x) { + rstan::extract(x$stan_fit, pars = "log_lik")$log_lik +} + +# more stable than log(sum(exp(x))) +#' @keywords internal +#' @noRd +.log_sum_exp <- function(x) { + max_x <- max(x) + max_x + log(sum(exp(x - max_x))) +} + +## more stable than log(mean(exp(x))) +##' @keywords internal +##' @noRd +.log_mean_exp <- function(x) { + .log_sum_exp(x) - log(length(x)) +} + +## compute log of raw importance ratios +## sums over observations *not* over posterior samples +##' @keywords internal +##' @noRd +.sum_log_ratios <- function(ll, ids = NULL) { + if (!is.null(ids)) + ll <- ll[, ids , drop = FALSE] + - rowSums(ll) +} + +##' obtain relative efficiency +##' @keywords internal +##' @noRd +.rel_eff <- function(ll, x) { + warmup <- x$stan_fit@sim$warmup + iter <- x$stan_fit@sim$iter + n_chains <- x$stan_fit@sim$chains + loo::relative_eff(exp(ll), chain_id = rep(1:n_chains, each = iter - warmup)) +} diff --git a/R/model_weights.R b/R/model_weights.R new file mode 100644 index 0000000..7df9de8 --- /dev/null +++ b/R/model_weights.R @@ -0,0 +1,110 @@ +#' Compute Model Weights Based on Leave-Future-Out Cross-Validation +#' +#' @description +#' This function computes model weights based on Leave-Future-Out (LFO) cross-validation +#' for a list of fitted GVAR models. It leverages the LFO method to evaluate predictive performance +#' and calculates weights using the \code{\link[loo:loo_model_weights]{loo::loo_model_weights}} function. +#' +#' @param fits A list of \code{tsnet_fit} objects or an object of class \code{tsnet_list}. +#' @inheritParams lfo.stan_gvar +#' @param return_lfo Logical. If \code{TRUE}, the function also returns the LFO objects. Default is \code{FALSE}. +#' @param ... Additional arguments (currently not used). +#' +#' @details +#' The function processes the input list of fitted models, ensuring compatibility with the \code{tsnet_list} class. +#' It computes the LFO cross-validation results for each model in the list using \code{\link{lfo.stan_gvar}} +#' and extracts log-likelihoods and relative efficiency information. The final model weights are then calculated +#' based on the log-likelihoods and relative efficiency values using the \code{\link[loo:loo_model_weights]{loo::loo_model_weights}} function. +#' +#' The LFO cross-validation options (\code{start}, \code{ahead}, \code{mode}, \code{k_thres}) +#' are inherited from \code{\link{lfo.stan_gvar}}. +#' +#' @return A list with the following components: +#' \describe{ +#' \item{weights}{The computed model weights.} +#' \item{ll_list}{A list of log-likelihoods for each model.} +#' \item{r_eff_list}{A list of relative efficiency information for each model.} +#' \item{lfo_list}{(Optional) A list of LFO objects if \code{return_lfo} is \code{TRUE}.} +#' } +#' +#' @importFrom loo loo_model_weights +#' +#' @examples +#' \dontrun{ +#' # Assuming `fit1` and `fit2` are tsnet_fit objects +#' fits <- tsnet_list(fit1, fit2) +#' weights <- model_weights(fits) +#' print(weights) +#' } +#' @keywords internal +#' @noRd + +model_weights <- function(fits, + start = 0, + ahead = 1, + mode = "backward", + k_thres = 0.7, + return_lfo = FALSE, + ...) { + + # browser() + # Check if input is a list + if (!is.list(fits)) { + stop("Input must be a list of tsnet_fit objects") + } + + # Check if input is tsnet_list + # if not, convert to tsnet_list + if (!inherits(fits, "tsnet_list")) { + fits <- tsnet_list(fits) + } + + ## Compute model weights based on LFO + lfo_list <- lapply( + fits, + FUN = lfo.stan_gvar, + type = "lfo", + start = start, + ahead = ahead, + mode = mode + ) + + # TODO this is not correct yet + # we need to extract the log-likelihoods from the LFO objects, not the elpd + # see https://github.com/ph-rast/bmgarch/blob/master/R/model_weights.R#L3 + ll_list <- lapply( + lfo_list, + FUN = function(x) + x$loglik + ) + + # Obtain rel_eff info + r_eff_list <- list() + for (i in seq_along(lfo_list)) { + r_eff_list[[i]] <- .rel_eff(ll_list[[i]], fits[[i]]) + } + + + # Compute model weights + # adapted from bmgarch::model_weights + loo_method = "stacking" + + weights <- loo::loo_model_weights( + ll_list, + method = loo_method, + r_eff_list = r_eff_list, + optim_control = list(reltol = 1e-10) + ) + + out <- list() + out$weights <- weights + out$ll_list <- ll_list + out$r_eff_list <- r_eff_list + + if (isTRUE(return_lfo)) { + out$lfo_list <- lfo_list + } + class(out) <- c("model_weights", class(out)) + return(out) + +} diff --git a/R/plot_compare_gvar.R b/R/plot_compare_gvar.R index c20a273..316dc89 100644 --- a/R/plot_compare_gvar.R +++ b/R/plot_compare_gvar.R @@ -1,10 +1,10 @@ #' Plot compare_gvar #' #' This function is a plotting method for the class produced by -#' [compare_gvar()]. It generates a plot showing the density of posterior +#' \code{\link{compare_gvar}}. It generates a plot showing the density of posterior #' uncertainty distributions for distances and the empirical distance value for two GVAR models. #' -#' @param x An object of class "compare_gvar". +#' @param x An object of class \code{compare_gvar}. #' @param name_a Optional. The name for model A. If provided, it replaces #' "mod_a" in the plot. #' @param name_b Optional. The name for model B. If provided, it replaces @@ -12,21 +12,21 @@ #' @param ... Additional arguments to be passed to the plotting functions. #' #' @details The function first checks if the full reference distributions of -#' [compare_gvar()] are saved using the argument 'return_all' set to TRUE. If +#' \code{\link{compare_gvar}} are saved using the argument \code{return_all} set to \code{TRUE}. If #' not, an error is thrown. #' #' Using the "name_a" and "name_b" arguments allows for custom labeling of the #' two models in the plot. #' -#' The function generates two density plots using `ggplot2`, one for the +#' The function generates two density plots using \code{ggplot2}, one for the #' temporal network (beta) and another for the contemporaneous network (pcor). #' The density distributions are filled with different colors based on the #' corresponding models (mod_a and mod_b). The empirical distances between the #' networks are indicated by red vertical lines. #' -#' @return A ggplot object representing the density plots of the posterior +#' @return A \code{ggplot} object representing the density plots of the posterior #' uncertainty distributions for distances and the empirical distance for two GVAR models. -#' +#' #' @import ggplot2 #' @importFrom ggokabeito scale_fill_okabe_ito palette_okabe_ito #' @importFrom cowplot plot_grid get_legend diff --git a/R/post_distance_within.R b/R/post_distance_within.R index 19abf95..16fec17 100644 --- a/R/post_distance_within.R +++ b/R/post_distance_within.R @@ -7,15 +7,16 @@ #' two modes. Distances can be obtained either from posterior samples or #' posterior predictive draws. The distance between two models can currently #' be calculated based on three options: Frobenius norm, maximum difference, -#' or L1 norm. Used within [compare_gvar()]. The function is not intended to +#' or L1 norm. Used within \code{\link{compare_gvar}}. The function is not intended to #' be used directly by the user. #' #' @inheritParams compare_gvar -#' @param fitobj Fitted model object. This can be a tsnet_fit object (obtained -#' from [stan_gvar()]), a BGGM object (obtained from [BGGM::var_estimate()]), -#' or extracted posterior samples (obtained from [stan_fit_convert()). +#' @param fitobj Fitted model object. This can be a tsnet_fit object +#' (obtained from \code{\link{stan_gvar}}, a \code{BGGM} object (obtained from +#' \code{\link[BGGM]{var_estimate}}, or extracted posterior samples (obtained from +#' \code{\link{stan_fit_convert}}). #' @param pred A logical indicating whether the input is posterior predictive -#' draws (TRUE) or posterior samples (FALSE). Default: FALSE +#' draws (TRUE) or posterior samples (FALSE). Default: \code{FALSE} #' @return A list of distances between the specified pairs of fitted models. The #' list has length equal to the specified number of random pairs. Each list #' element contains two distance values, one for beta coefficients and one for diff --git a/R/posterior_plot.R b/R/posterior_plot.R index ab689b5..7157127 100644 --- a/R/posterior_plot.R +++ b/R/posterior_plot.R @@ -4,20 +4,21 @@ #' or the contemporaneous networks of a GVAR model. The posterior distributions #' are visualized as densities in a matrix layout. #' -#' @param fitobj Fitted model object. This can be a tsnet_fit object (obtained -#' from [stan_gvar()]) or a BGGM object (obtained from [BGGM::var_estimate()]). +#' @param fitobj Fitted model object. This can be a tsnet_fit object +#' (obtained from \code{\link{stan_gvar}}, a \code{BGGM} object (obtained from +#' \code{\link[BGGM]{var_estimate}}. #' @param mat A matrix to use for plotting. Possibilities include "beta" -#' (temporal network) and "pcor" (contemporaneous network). Default is "beta" +#' (temporal network) and "pcor" (contemporaneous network). Default is \code{"beta"} #' (temporal network). #' @param cis A numeric vector of credible intervals to use for plotting. -#' Default is c(0.8, 0.9, 0.95). +#' Default is \code{c(0.8, 0.9, 0.95)}. #' #' @details In the returned plot, posterior distributions for every parameter #' are shown. Lagged variables are displayed along the vertical line of the #' grid, and non-lagged variables along the horizontal line of the grids. #' -#' @return A ggplot object representing the posterior distributions of the parameters of the temporal -#' or the contemporaneous networks of a GVAR model. +#' @return A \code{ggplot} object representing the posterior distributions of the +#' parameters of the temporal or the contemporaneous networks of a GVAR model. #' #' @import ggplot2 #' @importFrom ggdist stat_pointinterval stat_slab diff --git a/R/posterior_samples.R b/R/posterior_samples.R index c3a0091..9c56832 100644 --- a/R/posterior_samples.R +++ b/R/posterior_samples.R @@ -2,14 +2,14 @@ #' #' @description This function extracts the posterior samples of partial #' correlations and beta coefficients from a Bayesian GVAR model that was fitted -#' with [BGGM::var_estimate()] or [stan_gvar()]. +#' with \code{\link[BGGM]{var_estimate}} in \code{BGGM} or \code{\link{stan_gvar}}. #' The function is not intended to be called directly by the user, #' but is used internally by the package. #' -#' @param fitobj A [BGGM::var_estimate()] or [stan_gvar()] +#' @param fitobj A \code{\link[BGGM]{var_estimate}} or \code{\link{stan_gvar}} #' fit object from which to extract the posterior samples. #' @param burnin An integer specifying the number of initial samples to discard -#' as burn-in. Default is 0. +#' as burn-in. Default is \code{0}. #' #' @return A matrix containing the posterior samples of partial correlations and #' beta coefficients. The rows represent the samples, and the columns @@ -178,17 +178,17 @@ prepare_samples_plot <- function(fitobj, #' #' This function converts a Stan fit object into an array of samples for the #' temporal coefficients and the innovation covariance or partial correlation -#' matrices. It supports rstan as a backend. It can be used to convert models -#' fit using [stan_gvar()] into 3D arrays, which is the standard data structure -#' used in `tsnet`. The function allows to select which parameters should be +#' matrices. It supports \code{rstan} as a backend. It can be used to convert models +#' fit using \code{\link{stan_gvar}} into 3D arrays, which is the standard data structure +#' used in \code{tsnet}. The function allows to select which parameters should be #' returned. #' #' @param stan_fit A Stan fit object obtained from rstan or a tsnet_fit object -#' from [stan_gvar()]. +#' from \code{\link{stan_gvar}}. #' @param return_params A character vector specifying which parameters to #' return. Options are "beta" (temporal network), "sigma" (innovation #' covariance), and "pcor" (partial correlations). Default is -#' c("beta","sigma", "pcor"). +#' \code{c("beta","sigma", "pcor")}. #' #' @return A list containing 3D arrays for the selected parameters. Each array #' represents the posterior samples for a parameter, and each slice of the @@ -306,3 +306,57 @@ stan_fit_convert <- function(stan_fit, } + + + + + + +#' Obtain Weighted Posterior Samples +#' +#' @description +#' This internal function computes weighted posterior samples from multiple model fits. +#' It takes a list of fitted models, extracts posterior samples for specified parameters, +#' and combines them using model weights. Created based on \code{bmgarch:::.weighted_samples}. +#' +#' @param model_fit A list of model fit objects, typically of class \code{tsnet_fit}. +#' @param params A character vector specifying the names of parameters to extract from the models. +#' @param weights A numeric vector of model weights, typically obtained from \code{\link{model_weights}}. +#' The length of \code{weights} must match the length of \code{model_fit}. +#' +#' @details +#' For each parameter, the function extracts posterior samples from all model fits +#' and multiplies them by the corresponding model weights. The weighted samples are then summed +#' across models to produce combined posterior samples for each parameter. +#' +#' This function is intended for internal use in workflows that involve model averaging. +#' +#' @return A named list of weighted posterior samples for the specified parameters. +#' +#' @seealso \code{\link{model_weights}}, \code{\link{forecast.stan_gvar}} +#' +#' @keywords internal +#' @noRd +.weighted_samples <- function(model_fit, + params, + weights) { + # Extract all samples + samps <- lapply(model_fit, rstan::extract, pars = params) + # Apply weights + for (i in seq_len(length(samps))) { + # Each model + samps[[i]] <- lapply(samps[[i]], function(p) { + # Each parameter + p * weights[i] # Weight them + }) + } + # Reduce + samps_comb <- lapply(params, function(p) { + # For each parameter + Reduce("+", lapply(samps, function(m) { + m[[p]] + })) # Sum samples together + }) + names(samps_comb) <- params + return(samps_comb) +} diff --git a/R/print.R b/R/print.R index e3e1a24..e220b7e 100644 --- a/R/print.R +++ b/R/print.R @@ -1,13 +1,13 @@ #' Print method for compare_gvar objects #' -#' This function prints a summary of the Norm-Based Comparison Test for a [compare_gvar()] object. +#' This function prints a summary of the Norm-Based Comparison Test for a \code{\link{compare_gvar}} object. #' -#' @param x A test object obtained from [compare_gvar()] +#' @param x A test object obtained from \code{\link{compare_gvar}} #' @param ... Additional arguments to be passed to the print method. (currently not used) #' #' @return Prints a summary of the Norm-Based Comparison Test to the console #' -#' @details This function prints a summary of the Norm-Based Comparison Test for a [compare_gvar()] object. +#' @details This function prints a summary of the Norm-Based Comparison Test for a \code{\link{compare_gvar}} object. # It displays the general summary and model-specific results, including the number of significant comparisons #' in the temporal and contemporaneous networks, as well as the number of reference distances that were larger #' than the empirical distance for each network. @@ -61,7 +61,7 @@ print.compare_gvar <- function(x, #' Print method for tsnet_fit objects #' -#' This method provides a summary of the Bayesian GVAR model fitted with [stan_gvar()]. +#' This method provides a summary of the Bayesian GVAR model fitted with \code{\link{stan_gvar}}. #' It prints general information about the model, including the estimation method and the number of chains and iterations #' It also prints the posterior mean of the temporal and contemporaneous coefficients. #' diff --git a/R/stan_gvar.R b/R/stan_gvar.R index f89210f..5de4729 100644 --- a/R/stan_gvar.R +++ b/R/stan_gvar.R @@ -4,16 +4,17 @@ #' using Stan. The estimation procedure is described further in Siepe, Kloft & #' Heck (2023) . The current implementation allows #' for a normal prior on the temporal effects and either an Inverse Wishart or -#' an LKJ prior on the contemporaneous effects. `rstan` is used as a backend +#' an LKJ prior on the contemporaneous effects. \code{rstan} is used as a backend #' for fitting the model in Stan. Data should be provided in long format, where #' the columns represent the variables and the rows represent the time points. -#' Data are automatically z-scaled for estimation. +#' Data are automatically z-scaled for estimation. The model currently does not +#' support missing data. #' #' @param data A data frame or matrix containing the time series data of a #' single subject. The data should be in long format, where the columns #' represent the variables and the rows represent the time points. See the -#' example data [ts_data] for the correct format. -#' @param beep A vector of beeps with length of `nrow(data)`. The beep indicator +#' example data \code{\link{ts_data}} for the correct format. +#' @param beep A vector of beeps with length of \code{nrow(data)}. The beep indicator #' can be used to remove overnight effects from the last beep of a day to the #' first beep of the next day. This should be a vector of positive integers. #' If left empty, the function will assume that there are no overnight @@ -23,34 +24,40 @@ #' values corresponding to the prior distributions. The following priors can #' be specified: #' \itemize{ -#' \item `prior_Beta_loc` A matrix of the same dimensions as the beta matrix +#' \item \code{prior_Beta_loc} A matrix of the same dimensions as the beta matrix #' `B` containing the mean of the prior distribution for the beta coefficients. -#' \item `prior_Beta_scale` A matrix of the same dimensions as the beta matrix +#' \item \code{prior_Beta_scale} A matrix of the same dimensions as the beta matrix #' `B` containing the standard deviation of the prior distribution for the beta #' coefficients.} #' #' @param method A string indicating the method to use for fitting the model. #' Options are "sampling" (for MCMC estimation) or "variational" (for #' variational inference). We currently recommend only using MCMC estimation. +#' Default: \code{"sampling"}. #' @param cov_prior A string indicating the prior distribution to use for the -#' covariance matrix. Options are "LKJ" or "IW" (Inverse-Wishart). +#' covariance matrix. Options are "LKJ" or "IW" (Inverse-Wishart). Default: \code{"LKJ"}. #' @param rmv_overnight A logical indicating whether to remove overnight -#' effects. Default is `FALSE`. If `TRUE`, the function will remove overnight +#' effects. Default is \code{FALSE}. If `TRUE`, the function will remove overnight #' effects from the last beep of a day to the first beep of the next day. -#' This requires the `beep` argument to be specified. +#' This requires the \code{beep} argument to be specified. #' @param iter_sampling An integer specifying the number of iterations for the -#' sampling method. Default is 500. +#' sampling method. Default is \code{500}. #' @param iter_warmup An integer specifying the number of warmup iterations for -#' the sampling method. Default is 500. +#' the sampling method. Default is \code{500}. #' @param n_chains An integer specifying the number of chains for the sampling #' method. Default is 4. If variational inference is used, the number of -#' iterations is calculated as `iter_sampling`*`n_chains`. +#' iterations is calculated as \code{iter_sampling}*\code{n_chains}. #' @param n_cores An integer specifying the number of cores to use for parallel -#' computation. Default is 1. [rstan] is used for parallel computation. +#' computation. Default is \code{1}. \code{rstan} is used for parallel computation. #' @param center_only A logical indicating whether to only center (and not -#' scale) the data. Default is `FALSE`. -#' @param ... Additional arguments passed to the `rstan::sampling` or -#' `rstan::vb` function. +#' scale) the data. Default is \code{FALSE}. +#' @param return_all A logical indicating whether to return all model inputs, including +#' the data and prior objects. Default is \code{TRUE}. +#' @param ahead An integer specifying the forecast horizon. Default is \code{0}. If `ahead` is +#' greater than 0, the function will return posterior predictive forecasts for the specified number of time points. +#' This functionality is experimental and may not work as expected. +#' @param ... Additional arguments passed to the \code{\link[rstan:sampling]{rstan::sampling}} or +#' \code{\link[rstan:vb]{rstan::vb}} function. #' #' @details \bold{General Information} #' @@ -64,33 +71,33 @@ #' lagged effects). This is typically represented in the partial correlation #' (pcor) matrix. If the model is represented and interpreted as a network, #' variables are called -#' *nodes*, *edges* represent the statistical association between the nodes, and -#' *edge weights* quantify the strength of these associations. +#' \emph{nodes}, \emph{edges} represent the statistical association between the nodes, and +#' \emph{edge weights} quantify the strength of these associations. #' #' \bold{Model} #' #' Let \eqn{Y} be a matrix with \eqn{n} rows and \eqn{p} columns, where \eqn{n_t} is the #' number of time points and \eqn{p} is the number of variables. The GVAR model is -#' given by the following equations: \deqn{Y_t = B* Y_{t-1} + zeta_t} -#' \deqn{zeta_t \sim N(0, \Sigma)} where \eqn{B} is a `p x p` matrix of VAR -#' coefficients between variables i and j (beta_ij), \eqn{\zeta_t} contains the -#' innovations at time point \eqn{t}, and \eqn{\Sigma} is a `p x p`covariance matrix. +#' given by the following equations: \deqn{Y_t = B* Y_{t-1} + \zeta_t} +#' \deqn{\zeta_t \sim N(0, \Sigma)} where \eqn{B} is a \eqn{p \times p} matrix of VAR +#' coefficients between variables \eqn{i} and \eqn{j} (\eqn{\beta_{ij}}), \eqn{\zeta_t} contains the +#' innovations at time point \eqn{t}, and \eqn{\Sigma} is a \eqn{p \times p} covariance matrix. #' The inverse of \eqn{\Sigma} is the precision matrix, which is used to obtain the -#' partial correlations between variables (rho_ij). The model setup is +#' partial correlations between variables (\eqn{\rho_{ij}}). The model setup is #' explained in more detail in Siepe, Kloft & Heck (2023) #' . #' #' #' \bold{Prior Setup} #' -#' For the p x p temporal matrix B (containing the beta coefficients), we use +#' For the \eqn{p \times p} temporal matrix B (containing the \eqn{\beta} coefficients), we use #' a normal prior distribution on each individual parameter: \deqn{\beta_{ij} -#' \sim N(PriorBetaLoc_{ij}, PriorBetaScale_{ij})} where `PriorBetaLoc` is the -#' mean of the prior distribution and `PriorBetaScale` is the standard +#' \sim N(PriorBetaLoc_{ij}, PriorBetaScale_{ij})} where \code{PriorBetaLoc} is the +#' mean of the prior distribution and \code{PriorBetaScale} is the standard #' deviation of the prior distribution. The default prior is a weakly #' informative normal distribution with mean 0 and standard deviation 0.5. The #' user can specify a different prior distribution by a matrix -#' `prior_Beta_loc` and a matrix `prior_Beta_scale` with the same dimensions +#' \code{prior_Beta_loc} and a matrix \code{prior_Beta_scale} with the same dimensions #' as \eqn{B}. #' #' @@ -109,54 +116,55 @@ #' an additional beta prior on each of the partial correlations. We first #' assigned an uninformed LKJ prior to the Cholesky factor decomposition of #' the correlation matrix of innovations: \deqn{\Omega_L \sim -#' LKJ-Cholesky(\eta)}. For \eqn{\eta = 1}, this implies a symmetric marginal +#' LKJ-Cholesky(\eta)} +#' For \eqn{\eta = 1}, this implies a symmetric marginal #' scaled beta distribution on the zero-order correlations \eqn{\omega_{ij}}. #' \deqn{(\omega_{ij}+1)/2 \sim Beta(p/2, p/2)} #' We can then obtain the covariance matrix and, -#' subsequently, the precision matrix (see Siepe, Kloft & Heck (2023)) -#' for details. +#' subsequently, the precision matrix (see Siepe, Kloft & Heck, 2023, +#' for details). #' The second part of the prior is a beta prior on each partial correlation #' \eqn{\rho_{ij}} (obtained from the off-diagonal elements of the precision matrix). #' This prior was assigned by transforming the partial correlations to the -#' interval of 0,1 and then assigning a proportional (mean-variance +#' interval of \eqn{[0,1]} and then assigning a proportional (mean-variance #' parameterized) beta prior: #' \deqn{(\rho_{ij}+1)/2 \sim Beta_{prop}(PriorRhoLoc, PriorRhoScale)} #' A beta location parameter of 0.5 translates to an expected correlation of 0. -#' The variance parameter of sqrt(0.5) implies a uniform distribution of +#' The variance parameter of \eqn{\sqrt(0.5)} implies a uniform distribution of #' partial correlations. #' The user can specify a different prior distribution by a matrix -#' `prior_Rho_loc` and a matrix `prior_Rho_scale` with the same dimensions as +#' \code{prior_Rho_loc} and a matrix \code{prior_Rho_scale} with the same dimensions as #' the partial correlation matrix. Additionally, the user can change \eqn{eta} -#' via the `prior_Eta` parameter. +#' via the \code{prior_Eta} parameter. #' #' The Inverse-Wishart prior is a distribution on the innovation covariance -#' matrix `Sigma`: +#' matrix \eqn{\Sigma}: #' \deqn{\Sigma \sim IW(\nu, S)} #' where \eqn{\nu} is the degrees of freedom and \eqn{S} is the scale matrix. We here #' use the default prior of \deqn{nu = delta + p - 1} for the degrees of freedom, #' where \eqn{\delta} is defined as \eqn{s_{\rho}^{-1}-1} and \eqn{s_{\rho}} is the #' standard deviation of the implied marginal beta distribution of the #' partial correlations. For the scale matrix \eqn{S}, we use the identity matrix -#' \eqn{I_p} of order p. +#' \eqn{I_p} of order \eqn{p}. #' The user can set a prior on the expected standard deviation of the partial -#' correlations by specifying a `prior_Rho_marginal` parameter. The default +#' correlations by specifying a \code{prior_Rho_marginal} parameter. The default #' value is 0.25, which has worked well in a simulation study. -#' Additionally, the user can specify a `prior_S` parameter to set a different +#' Additionally, the user can specify a \code{prior_S} parameter to set a different #' scale matrix. #' #' \bold{Sampling} #' The model can be fitted using either MCMC sampling or variational -#' inference via [rstan]. Per default, the model is fitted using the Stan +#' inference via \code{rstan}. Per default, the model is fitted using the Stan #' Hamiltonian Monte Carlo (HMC) No U-Turn (NUTS) sampler with 4 chains, #' 500 warmup iterations and 500 sampling iterations. We use a default -#' target average acceptance probability `adapt_delta` of 0.8. As the output -#' is returned as a standard `stanfit` object, the user can use the -#' `rstan` package to extract and analyze the results and obtain convergence +#' target average acceptance probability \code{adapt_delta} of 0.8. As the output +#' is returned as a standard \code{stanfit} object, the user can use the +#' \code{rstan} package to extract and analyze the results and obtain convergence #' diagnostics. #' #' #' -#' @return A `tsnet_fit` object in list format. The object contains the +#' @return A \code{tsnet_fit} object in list format. The object contains the #' following elements: #' \item{fit}{A stanfit object containing the fitted model.} #' \item{arguments}{The number of variables "p", the number of time points "n_t", the column names "cnames", and the arguments used in the function call.} @@ -178,166 +186,151 @@ #'} #' @export -stan_gvar <- - function(data, - beep = NULL, - priors = NULL, - method = "sampling", - cov_prior = "IW", - rmv_overnight = FALSE, - iter_sampling = 500, - iter_warmup = 500, - n_chains = 4, - n_cores = 1, - center_only = FALSE, - ...) { - if(isTRUE(center_only)){ - Y <- apply(data, MARGIN = 2, scale, center = TRUE, scale = FALSE) - } else{ - Y <- apply(data, MARGIN = 2, scale, center = TRUE, scale = TRUE) - } +stan_gvar <- function(data, + beep = NULL, + priors = NULL, + method = "sampling", + cov_prior = "IW", + rmv_overnight = FALSE, + iter_sampling = 500, + iter_warmup = 500, + n_chains = 4, + n_cores = 1, + center_only = FALSE, + return_all = TRUE, + ahead = 0, + ...) { + if (isTRUE(center_only)) { + Y <- apply(data, MARGIN = 2, scale, center = TRUE, scale = FALSE) + } else { + Y <- apply(data, MARGIN = 2, scale, center = TRUE, scale = TRUE) + } - K <- ncol(data) - n_t <- nrow(data) - cnames <- colnames(data) + K <- ncol(data) + n_t <- nrow(data) + cnames <- colnames(data) - # Store arguments - mc <- match.call() - fl <- formals() - missing_args <- setdiff(names(fl), names(mc)) - missing_args <- Map(function(arg, default) - if (!is.null(default)) mc[[arg]] <- default, missing_args, fl[missing_args]) - all_args <- c(as.list(mc), missing_args) + # Store arguments + mc <- match.call() + fl <- formals() + missing_args <- setdiff(names(fl), names(mc)) + missing_args <- Map(function(arg, default) + if (!is.null(default)) mc[[arg]] <- default, missing_args, fl[missing_args]) + all_args <- c(as.list(mc), missing_args) - if(is.null(beep)){ - beep <- seq(1, n_t) - } + if (is.null(beep)) { + beep <- seq(1, n_t) + } + # Initialize priors with default values + default_priors <- list( + prior_Beta_loc = matrix(rep(0, K * K), nrow = K, ncol = K), + prior_Beta_scale = matrix(rep(0.5, K * K), nrow = K, ncol = K), + prior_Rho_loc = matrix(rep(0.5, K * K), nrow = K, ncol = K), + prior_Rho_scale = matrix(rep(sqrt(0.5), K * K), nrow = K, ncol = K), + prior_Rho_marginal = 0.25, + prior_S = diag(K), + prior_Eta = 1 + ) - # Specify Priors - if(is.null(priors[["prior_Beta_loc"]])){ - prior_Beta_loc <- matrix(rep(0, K*K), nrow = K, ncol = K) - } else { - prior_Beta_loc <- priors[["prior_Beta_loc"]] - } - - if(is.null(priors[["prior_Beta_scale"]])){ - prior_Beta_scale <- matrix(rep(.5, K*K), nrow = K, ncol = K) - } else { - prior_Beta_scale <- priors[["prior_Beta_scale"]] - } + # Update default priors with any custom priors provided by the user + if (!is.null(priors)) { + priors <- modifyList(default_priors, priors) + } else { + priors <- default_priors + } - if(is.null(priors[["prior_Rho_loc"]])){ - prior_Rho_loc <- matrix(rep(.5, K*K), nrow = K, ncol = K) - } else { - prior_Rho_loc <- priors[["prior_Rho_loc"]] - } + # Convert SD to delta: SD = 1/(delta+1), delta = (1 / SD) - 1 + prior_delta <- (1 / priors$prior_Rho_marginal) - 1 - if(is.null(priors[["prior_Rho_scale"]])){ - prior_Rho_scale <- matrix(rep(sqrt(.5), K*K), nrow = K, ncol = K) - } else { - prior_Rho_scale <- priors[["prior_Rho_scale"]] - } + # Forecasting + if (ahead > 0){ + compute_log_lik <- 1 + } else { + compute_log_lik <- 0 + } - if(is.null(priors[["prior_Rho_marginal"]])){ - prior_Rho_marginal <- 0.25 - } else { - prior_Rho_marginal <- priors[["prior_Rho_marginal"]] - } + # Choose model to fit + if (cov_prior == "LKJ") { + # Stan Data + stan_data <- list( + K = K, + "T" = n_t, + Y = as.matrix(Y), + beep = beep, + prior_Rho_loc = priors$prior_Rho_loc, + prior_Rho_scale = priors$prior_Rho_scale, + prior_Beta_loc = priors$prior_Beta_loc, + prior_Beta_scale = priors$prior_Beta_scale, + prior_Eta = priors$prior_Eta, + ahead = ahead, + Y_future = matrix(rep(0, K), ncol = K), + compute_log_lik = compute_log_lik + ) - if(is.null(priors[["prior_S"]])){ - prior_S <- diag(K) + if (isTRUE(rmv_overnight)) { + # remove overnight effects + model_name <- "VAR_LKJ_beep" } else { - prior_S <- priors[["prior_S"]] + # standard model + model_name <- "VAR_LKJ" } + } + if (cov_prior == "IW") { + # Stan Data + stan_data <- list( + K = K, + "T" = n_t, + Y = as.matrix(Y), + beep = beep, + prior_Beta_loc = priors$prior_Beta_loc, + prior_Beta_scale = priors$prior_Beta_scale, + prior_S = priors$prior_S, + prior_delta = prior_delta, + ahead = ahead, + Y_future = matrix(rep(0, K), ncol = K), + compute_log_lik = compute_log_lik + ) - if(is.null(priors[["prior_Eta"]])){ - prior_Eta <- 1 + if (isTRUE(rmv_overnight)) { + # remove overnight effects + model_name <- "VAR_wishart_beep" } else { - prior_Eta <- priors[["prior_Eta"]] - } - - # Convert SD to delta: SD = 1/(delta+1), delta = (1 / SD) - 1 - prior_delta <- (1 / prior_Rho_marginal) - 1 - - - - # Choose model to fit - if (cov_prior == "LKJ") { - # Stan Data - stan_data <- list( - K = K, - "T" = n_t, - Y = as.matrix(Y), - beep = beep, - prior_Rho_loc = prior_Rho_loc, - prior_Rho_scale = prior_Rho_scale, - prior_Beta_loc = prior_Beta_loc, - prior_Beta_scale = prior_Beta_scale, - prior_Eta = prior_Eta - ) - - if (isTRUE(rmv_overnight)) { - # remove overnight effects - model_name <- "VAR_LKJ_beep" - } else{ - # standard model - model_name <- "VAR_LKJ" - } + # standard model + model_name <- "VAR_wishart" } - if (cov_prior == "IW") { - # Stan Data - stan_data <- list( - K = K, - "T" = n_t, - Y = as.matrix(Y), - beep = beep, - prior_Beta_loc = prior_Beta_loc, - prior_Beta_scale = prior_Beta_scale, - prior_S = prior_S, - prior_delta = prior_delta - ) + } - if (isTRUE(rmv_overnight)) { - # remove overnight effects - model_name <- "VAR_wishart_beep" - } else{ - # standard model - model_name <- "VAR_wishart" - } - } - - if (method == "sampling") { - default_args <- list( - object = stanmodels[[model_name]], - data = stan_data, - chains = n_chains, - cores = n_cores, - iter = iter_sampling + iter_warmup, - warmup = iter_warmup, - refresh = 500, - thin = 1, - init = .1, - control = list(adapt_delta = .8) - ) + if (method == "sampling") { + default_args <- list( + object = stanmodels[[model_name]], + data = stan_data, + chains = n_chains, + cores = n_cores, + iter = iter_sampling + iter_warmup, + warmup = iter_warmup, + refresh = 500, + thin = 1, + init = .1, + control = list(adapt_delta = .8) + ) - # Run sampler - stan_fit <- do.call(rstan::sampling, - utils::modifyList(default_args, list(...))) + stan_fit <- do.call(rstan::sampling, + utils::modifyList(default_args, list(...))) - } - if (method == "variational") { - default_args <- list( - object = stanmodels[[model_name]], - data = stan_data, - init = .1, - tol_rel_obj = .001, - output_samples = iter_sampling * n_chains - ) + } + if (method == "variational") { + default_args <- list( + object = stanmodels[[model_name]], + data = stan_data, + init = .1, + tol_rel_obj = .001, + output_samples = iter_sampling * n_chains + ) - stan_fit <- do.call(rstan::vb, - utils::modifyList(default_args, list(...))) - } + stan_fit <- do.call(rstan::vb, + utils::modifyList(default_args, list(...))) + } args <- list(p = K, @@ -345,6 +338,11 @@ stan_gvar <- cnames = cnames, fn_args = all_args ) + if(isTRUE(return_all)){ + args$data <- stan_data$Y + args$priors <- priors + } + ret_fit <- list( stan_fit = stan_fit, diff --git a/R/tsnet-package.R b/R/tsnet-package.R index 983ad25..a82571d 100644 --- a/R/tsnet-package.R +++ b/R/tsnet-package.R @@ -2,7 +2,6 @@ #' #' @description Time Series Network Analysis with R #' -#' @docType package #' @name tsnet-package #' @aliases tsnet #' @useDynLib tsnet, .registration = TRUE @@ -13,4 +12,4 @@ #' @importFrom RcppParallel RcppParallelLibs #' #' -NULL +"_PACKAGE" diff --git a/README.Rmd b/README.Rmd index a53ba41..76a023b 100644 --- a/README.Rmd +++ b/README.Rmd @@ -68,7 +68,7 @@ print(fit_stan) ### Comparing Network Models -This is an example of how to use the package to compare two network models. We use here BGGM to estimate the networks, but the `stan_gvar` function can be used as well. +This is an example of how to use the package to compare two network models. We here use `BGGM` to estimate the networks, but the `stan_gvar` function can be used as well. ```{r example-test, eval = FALSE} library(tsnet) diff --git a/README.md b/README.md index 8d69973..ef46aaf 100644 --- a/README.md +++ b/README.md @@ -69,7 +69,7 @@ print(fit_stan) ### Comparing Network Models This is an example of how to use the package to compare two network -models. We use here BGGM to estimate the networks, but the `stan_gvar` +models. We here use `BGGM` to estimate the networks, but the `stan_gvar` function can be used as well. ``` r diff --git a/cran-comments.md b/cran-comments.md index 15585ca..b42b59a 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -2,4 +2,4 @@ 0 errors | 0 warnings | 1 note -* This is a new release. Fixed previous documentation mistakes. +* This is a new release. Fixed missing roxygen documentation. diff --git a/inst/stan/VAR_LKJ_beep.stan b/inst/stan/VAR_LKJ_beep.stan index 94b6d62..4ba72d2 100644 --- a/inst/stan/VAR_LKJ_beep.stan +++ b/inst/stan/VAR_LKJ_beep.stan @@ -12,6 +12,11 @@ data { matrix[K,K] prior_Rho_loc; // locations for priors on partial correlations matrix[K,K] prior_Rho_scale; // scales for priors on partial correlations int prior_Eta; // prior for LKJ + // Forecast + int ahead; // forecasted time points + int compute_log_lik; // compute log likelihood? + array[ahead] vector[K] Y_future; // future responses for loglik computation + } //////////////////////////////////////////////////////////////////////////////// transformed data{ @@ -83,17 +88,38 @@ model { } } //////////////////////////////////////////////////////////////////////////////// -generated quantities{ +generated quantities { int min_beep = first_beep; - vector[T-1] log_lik; - { - // Cholesky decomposition of the covariance matrix - matrix[K, K] Sigma_chol = diag_pre_multiply(exp(sigma_theta), L_Theta); - for(t in 2:T){ - if(beep[t] > first_beep){ - vector[K] mu = Beta * Y[t-1,]; - log_lik[t-1] = multi_normal_cholesky_lpdf(Y[t, ] | mu, Sigma_chol); - } + vector[T + ahead] log_lik; + + // Initialize log_lik with zero for all time points + for (t in 1:(T + ahead)) { + log_lik[t] = 0; + } + + // Cholesky decomposition of the covariance matrix + matrix[K, K] Sigma_chol = diag_pre_multiply(exp(sigma_theta), L_Theta); + for (t in 2:T) { + if (beep[t] > first_beep) { + vector[K] mu = Beta * Y[t-1,]; + // use t here meaning that ll of timepoint two is stored in log_lik[2] + log_lik[t] = multi_normal_cholesky_lpdf(Y[t, ] | mu, Sigma_chol); } } + if(ahead > 0){ + // Forecasting ahead steps + array[ahead] vector[K] Y_forecast; // forecasted responses + vector[K] current_Y = Y[T]; // initialize current_Y to the last observed value + for (s in 1:ahead) { + vector[K] mu = Beta * current_Y; // mu for current step + Y_forecast[s] = multi_normal_cholesky_rng(mu, Sigma_chol); // generate forecast + print(Y_forecast[s]); + current_Y = Y_forecast[s]; // update current_Y to predicted value + // forecasted log-likelihood + if (compute_log_lik == 1) { + log_lik[T + s] = multi_normal_cholesky_lpdf(Y_future[s] | mu, Sigma_chol); + } + } + } + } diff --git a/man/check_eigen.Rd b/man/check_eigen.Rd index 974f295..0f6ff9b 100644 --- a/man/check_eigen.Rd +++ b/man/check_eigen.Rd @@ -8,12 +8,12 @@ check_eigen(fitobj, verbose = TRUE) } \arguments{ \item{fitobj}{A fitted Bayesian GVAR object. This can be a tsnet_fit object -(obtained from [stan_gvar()]), a BGGM object (obtained from -[BGGM::var_estimate()]), or extracted posterior samples (obtained from -[stan_fit_convert()).} +(obtained from \code{\link{stan_gvar}}, a \code{BGGM} object (obtained from +\code{\link[BGGM]{var_estimate}}), or extracted posterior samples +(obtained from \code{\link{stan_fit_convert}}.} \item{verbose}{Logical. If TRUE, a verbal summary of the results is printed. -Default is TRUE.} +Default is \code{TRUE}.} } \value{ A list containing the eigenvalues and a verbal summary of the @@ -22,7 +22,7 @@ A list containing the eigenvalues and a verbal summary of the \description{ This function checks the eigenvalues of the Beta matrix (containing the temporal coefficients) to assure that the model is stationary. It uses the -same check as the `graphicalVAR` package. The function calculates the +same check as the \code{graphicalVAR} package. The function calculates the eigenvalues of the Beta matrix and checks if the sum of the squares of the real and imaginary parts of the eigenvalues is less than 1. If it is, the VAR model is considered stable. diff --git a/man/compare_gvar.Rd b/man/compare_gvar.Rd index b0f587a..eef12bc 100644 --- a/man/compare_gvar.Rd +++ b/man/compare_gvar.Rd @@ -19,50 +19,50 @@ compare_gvar( } \arguments{ \item{fit_a}{Fitted model object for Model A. This can be a tsnet_fit object -(obtained from [stan_gvar()]), a BGGM object (obtained from -[BGGM::var_estimate()]), or extracted posterior samples (obtained from -[stan_fit_convert()).} +(obtained from \code{\link{stan_gvar}}, a \code{BGGM} object (obtained from +\code{\link[BGGM]{var_estimate}}, or extracted posterior samples (obtained from +\code{\link{stan_fit_convert}}).} \item{fit_b}{Fitted model object for Model B. This can be a tsnet_fit object -(obtained from [stan_gvar()]), a BGGM object (obtained from -[BGGM::var_estimate()]), or extracted posterior samples (obtained from -[stan_fit_convert()).} +(obtained from \code{\link{stan_gvar}}, a \code{BGGM} object (obtained from +\code{\link[BGGM]{var_estimate}}, or extracted posterior samples (obtained from +\code{\link{stan_fit_convert}}).} -\item{cutoff}{The percentage level of the test (default: 5\%) as integer.} +\item{cutoff}{The percentage level of the test (default: \code{5}) as integer.} -\item{dec_rule}{The decision rule to be used. Currently supports default "or" +\item{dec_rule}{The decision rule to be used. Currently supports default \code{"or"} (comparing against two reference distributions) and "comb" (combining the reference distributions). The use of "or" is recommended, as "comb" is less stable.} \item{n_draws}{The number of draws to use for reference distributions -(default: 1000).} +(default: \code{1000}).} \item{comp}{The distance metric to use. Should be one of "frob" (Frobenius norm), "maxdiff" (maximum difference), or "l1" (L1 norm) (default: -"frob"). The use of the Frobenius norm is recommended.} +\code{"frob"}). The use of the Frobenius norm is recommended.} \item{return_all}{Logical indicating whether to return all distributions -(default: FALSE). Has to be set to TRUE for plotting the results.} +(default: \code{FALSE}). Has to be set to TRUE for plotting the results.} \item{sampling_method}{Draw sequential pairs of samples from the posterior, with certain distance between them ("sequential") or randomly from two halves of the posterior ("random"). The "random" method is preferred to account for potential autocorrelation between subsequent samples. Default: -"random".} +\code{"random"}.} \item{indices}{A list of "beta" and "pcor" indices specifying which elements -of the matrices to consider when calculating distances. If NULL (default), +of the matrices to consider when calculating distances. If \code{NULL} (default), all elements of both matrices are considered. If provided, only the elements at these indices are considered. If only one of the matrices should have indices, the other one should be NULL. This can be useful if you want to calculate distances based on a subset of the elements in the matrices.} -\item{burnin}{The number of burn-in iterations to discard (default: 0).} +\item{burnin}{The number of burn-in iterations to discard (default: \code{0}).} } \value{ -A list (of class "compare_gvar") containing the results of the +A list (of class \code{compare_gvar}) containing the results of the comparison. The list includes: \item{sig_beta}{Binary decision on whether there is a significant difference between the temporal networks of A and B} diff --git a/man/fit_data.Rd b/man/fit_data.Rd index 6f898d1..bb7b474 100644 --- a/man/fit_data.Rd +++ b/man/fit_data.Rd @@ -9,19 +9,19 @@ A list with two elements, each containing posterior samples for one individual. } \source{ -The data is generated using the [stan_gvar()] function on subsets -of the [ts_data] time series data. +The data is generated using the \code{\link{stan_gvar}} function on subsets +of the \code{\link{ts_data}} time series data. } \usage{ data(fit_data) } \description{ This dataset contains posterior samples of beta coefficients and partial correlations for two individuals. -It was generated by fitting a GVAR model using [stan_gvar()] with three variables from the [ts_data] dataset. +It was generated by fitting a GVAR model using \code{\link{stan_gvar}} with three variables from the \code{\link{ts_data}} dataset. } \details{ The list contains two elements, each containing posterior samples for one individual. -The samples were extracted using the [stan_fit_convert()] function. +The samples were extracted using the \code{\link{stan_fit_convert}} function. For each individual, the list elements contain the posterior means of the beta coefficients ("beta_mu") and the posterior means of the partial correlations ("pcor_mu"). The "fit" element contains all 1000 posterior samples of the beta coefficients and partial correlations. diff --git a/man/get_centrality.Rd b/man/get_centrality.Rd index 533269a..d293707 100644 --- a/man/get_centrality.Rd +++ b/man/get_centrality.Rd @@ -8,16 +8,17 @@ get_centrality(fitobj, burnin = 0, remove_ar = TRUE) } \arguments{ \item{fitobj}{Fitted model object for a Bayesian GVAR model. This can be -`tsnet_fit` object (obtained from [stan_gvar()]), a BGGM object (obtained -from [BGGM::var_estimate()]), or extracted posterior samples (obtained from -[stan_fit_convert()).} +`tsnet_fit` object (obtained from \code{\link{stan_gvar}}, +a BGGM object (obtained from \code{\link[BGGM]{var_estimate}} in \code{BGGM}), +or extracted posterior samples (obtained from \code{\link{stan_fit_convert}}).} \item{burnin}{An integer specifying the number of initial samples to discard -as burn-in. Default is 0.} +as burn-in. Default is \code{0}.} \item{remove_ar}{A logical value specifying whether to remove the -autoregressive effects for centrality calculation. Default is TRUE. This is -only relevant for the calculation of temporal centrality/density measures.} +autoregressive effects for centrality calculation. Default is \code{TRUE}. +This is only relevant for the calculation of temporal centrality/density +measures.} } \value{ A list containing the following centrality measures: @@ -35,7 +36,7 @@ fit object. Centrality measures describe the "connectedness" of a variable in a network, while density describes the networks' overall connectedness. Specifically, it computes the in-strength, out-strength, contemporaneous strength, temporal network density, and contemporaneous network density. The -result can then be visualized using [plot_centrality()]. +result can then be visualized using \code{\link{plot_centrality}}. } \examples{ # Use first individual from example fit data from tsnet diff --git a/man/plot.compare_gvar.Rd b/man/plot.compare_gvar.Rd index d3efee3..b4cf6cf 100644 --- a/man/plot.compare_gvar.Rd +++ b/man/plot.compare_gvar.Rd @@ -7,7 +7,7 @@ \method{plot}{compare_gvar}(x, name_a = NULL, name_b = NULL, ...) } \arguments{ -\item{x}{An object of class "compare_gvar".} +\item{x}{An object of class \code{compare_gvar}.} \item{name_a}{Optional. The name for model A. If provided, it replaces "mod_a" in the plot.} @@ -18,23 +18,23 @@ \item{...}{Additional arguments to be passed to the plotting functions.} } \value{ -A ggplot object representing the density plots of the posterior +A \code{ggplot} object representing the density plots of the posterior uncertainty distributions for distances and the empirical distance for two GVAR models. } \description{ This function is a plotting method for the class produced by -[compare_gvar()]. It generates a plot showing the density of posterior +\code{\link{compare_gvar}}. It generates a plot showing the density of posterior uncertainty distributions for distances and the empirical distance value for two GVAR models. } \details{ The function first checks if the full reference distributions of -[compare_gvar()] are saved using the argument 'return_all' set to TRUE. If +\code{\link{compare_gvar}} are saved using the argument \code{return_all} set to \code{TRUE}. If not, an error is thrown. Using the "name_a" and "name_b" arguments allows for custom labeling of the two models in the plot. -The function generates two density plots using `ggplot2`, one for the +The function generates two density plots using \code{ggplot2}, one for the temporal network (beta) and another for the contemporaneous network (pcor). The density distributions are filled with different colors based on the corresponding models (mod_a and mod_b). The empirical distances between the diff --git a/man/plot_centrality.Rd b/man/plot_centrality.Rd index 92eac57..8839150 100644 --- a/man/plot_centrality.Rd +++ b/man/plot_centrality.Rd @@ -8,13 +8,13 @@ plot_centrality(obj, plot_type = "tiefighter", cis = 0.95) } \arguments{ \item{obj}{An object containing the centrality measures obtained from -[get_centrality()].} +\code{\link{get_centrality}}.} \item{plot_type}{A character string specifying the type of plot. Accepts -"tiefighter" or "density". Default is "tiefighter".} +"tiefighter" or "density". Default is \code{"tiefighter"}.} \item{cis}{A numeric value specifying the credible interval. Must be between -0 and 1 (exclusive). Default is 0.95.} +0 and 1 (exclusive). Default is \code{0.95}.} } \value{ A ggplot object visualizing the centrality measures. For a diff --git a/man/post_distance_within.Rd b/man/post_distance_within.Rd index 02d2e0e..d2b8dfb 100644 --- a/man/post_distance_within.Rd +++ b/man/post_distance_within.Rd @@ -16,35 +16,36 @@ post_distance_within( ) } \arguments{ -\item{fitobj}{Fitted model object. This can be a tsnet_fit object (obtained -from [stan_gvar()]), a BGGM object (obtained from [BGGM::var_estimate()]), -or extracted posterior samples (obtained from [stan_fit_convert()).} +\item{fitobj}{Fitted model object. This can be a tsnet_fit object +(obtained from \code{\link{stan_gvar}}, a \code{BGGM} object (obtained from +\code{\link[BGGM]{var_estimate}}, or extracted posterior samples (obtained from +\code{\link{stan_fit_convert}}).} \item{comp}{The distance metric to use. Should be one of "frob" (Frobenius norm), "maxdiff" (maximum difference), or "l1" (L1 norm) (default: -"frob"). The use of the Frobenius norm is recommended.} +\code{"frob"}). The use of the Frobenius norm is recommended.} \item{pred}{A logical indicating whether the input is posterior predictive -draws (TRUE) or posterior samples (FALSE). Default: FALSE} +draws (TRUE) or posterior samples (FALSE). Default: \code{FALSE}} \item{n_draws}{The number of draws to use for reference distributions -(default: 1000).} +(default: \code{1000}).} \item{sampling_method}{Draw sequential pairs of samples from the posterior, with certain distance between them ("sequential") or randomly from two halves of the posterior ("random"). The "random" method is preferred to account for potential autocorrelation between subsequent samples. Default: -"random".} +\code{"random"}.} \item{indices}{A list of "beta" and "pcor" indices specifying which elements -of the matrices to consider when calculating distances. If NULL (default), +of the matrices to consider when calculating distances. If \code{NULL} (default), all elements of both matrices are considered. If provided, only the elements at these indices are considered. If only one of the matrices should have indices, the other one should be NULL. This can be useful if you want to calculate distances based on a subset of the elements in the matrices.} -\item{burnin}{The number of burn-in iterations to discard (default: 0).} +\item{burnin}{The number of burn-in iterations to discard (default: \code{0}).} } \value{ A list of distances between the specified pairs of fitted models. The @@ -59,7 +60,7 @@ This function computes distances between posterior samples of a two modes. Distances can be obtained either from posterior samples or posterior predictive draws. The distance between two models can currently be calculated based on three options: Frobenius norm, maximum difference, - or L1 norm. Used within [compare_gvar()]. The function is not intended to + or L1 norm. Used within \code{\link{compare_gvar}}. The function is not intended to be used directly by the user. } \examples{ diff --git a/man/posterior_plot.Rd b/man/posterior_plot.Rd index 083751a..5ef8945 100644 --- a/man/posterior_plot.Rd +++ b/man/posterior_plot.Rd @@ -7,19 +7,20 @@ posterior_plot(fitobj, mat = "beta", cis = c(0.8, 0.9, 0.95)) } \arguments{ -\item{fitobj}{Fitted model object. This can be a tsnet_fit object (obtained -from [stan_gvar()]) or a BGGM object (obtained from [BGGM::var_estimate()]).} +\item{fitobj}{Fitted model object. This can be a tsnet_fit object +(obtained from \code{\link{stan_gvar}}, a \code{BGGM} object (obtained from +\code{\link[BGGM]{var_estimate}}.} \item{mat}{A matrix to use for plotting. Possibilities include "beta" -(temporal network) and "pcor" (contemporaneous network). Default is "beta" +(temporal network) and "pcor" (contemporaneous network). Default is \code{"beta"} (temporal network).} \item{cis}{A numeric vector of credible intervals to use for plotting. -Default is c(0.8, 0.9, 0.95).} +Default is \code{c(0.8, 0.9, 0.95)}.} } \value{ -A ggplot object representing the posterior distributions of the parameters of the temporal -or the contemporaneous networks of a GVAR model. +A \code{ggplot} object representing the posterior distributions of the +parameters of the temporal or the contemporaneous networks of a GVAR model. } \description{ Plots posterior distributions of the parameters of the temporal diff --git a/man/print.compare_gvar.Rd b/man/print.compare_gvar.Rd index e0863ed..b868aff 100644 --- a/man/print.compare_gvar.Rd +++ b/man/print.compare_gvar.Rd @@ -7,7 +7,7 @@ \method{print}{compare_gvar}(x, ...) } \arguments{ -\item{x}{A test object obtained from [compare_gvar()]} +\item{x}{A test object obtained from \code{\link{compare_gvar}}} \item{...}{Additional arguments to be passed to the print method. (currently not used)} } @@ -15,10 +15,10 @@ Prints a summary of the Norm-Based Comparison Test to the console } \description{ -This function prints a summary of the Norm-Based Comparison Test for a [compare_gvar()] object. +This function prints a summary of the Norm-Based Comparison Test for a \code{\link{compare_gvar}} object. } \details{ -This function prints a summary of the Norm-Based Comparison Test for a [compare_gvar()] object. +This function prints a summary of the Norm-Based Comparison Test for a \code{\link{compare_gvar}} object. in the temporal and contemporaneous networks, as well as the number of reference distances that were larger than the empirical distance for each network. } diff --git a/man/print.tsnet_fit.Rd b/man/print.tsnet_fit.Rd index 075f6eb..e432a75 100644 --- a/man/print.tsnet_fit.Rd +++ b/man/print.tsnet_fit.Rd @@ -15,7 +15,7 @@ Prints a summary to the console. } \description{ -This method provides a summary of the Bayesian GVAR model fitted with [stan_gvar()]. +This method provides a summary of the Bayesian GVAR model fitted with \code{\link{stan_gvar}}. It prints general information about the model, including the estimation method and the number of chains and iterations It also prints the posterior mean of the temporal and contemporaneous coefficients. } diff --git a/man/stan_fit_convert.Rd b/man/stan_fit_convert.Rd index c12950c..2c7fe38 100644 --- a/man/stan_fit_convert.Rd +++ b/man/stan_fit_convert.Rd @@ -8,12 +8,12 @@ stan_fit_convert(stan_fit, return_params = c("beta", "sigma", "pcor")) } \arguments{ \item{stan_fit}{A Stan fit object obtained from rstan or a tsnet_fit object -from [stan_gvar()].} +from \code{\link{stan_gvar}}.} \item{return_params}{A character vector specifying which parameters to return. Options are "beta" (temporal network), "sigma" (innovation covariance), and "pcor" (partial correlations). Default is -c("beta","sigma", "pcor").} +\code{c("beta","sigma", "pcor")}.} } \value{ A list containing 3D arrays for the selected parameters. Each array @@ -23,9 +23,9 @@ A list containing 3D arrays for the selected parameters. Each array \description{ This function converts a Stan fit object into an array of samples for the temporal coefficients and the innovation covariance or partial correlation -matrices. It supports rstan as a backend. It can be used to convert models -fit using [stan_gvar()] into 3D arrays, which is the standard data structure -used in `tsnet`. The function allows to select which parameters should be +matrices. It supports \code{rstan} as a backend. It can be used to convert models +fit using \code{\link{stan_gvar}} into 3D arrays, which is the standard data structure +used in \code{tsnet}. The function allows to select which parameters should be returned. } \examples{ diff --git a/man/stan_gvar.Rd b/man/stan_gvar.Rd index 4f2ccac..e18b65a 100644 --- a/man/stan_gvar.Rd +++ b/man/stan_gvar.Rd @@ -16,6 +16,8 @@ stan_gvar( n_chains = 4, n_cores = 1, center_only = FALSE, + return_all = TRUE, + ahead = 0, ... ) } @@ -23,9 +25,9 @@ stan_gvar( \item{data}{A data frame or matrix containing the time series data of a single subject. The data should be in long format, where the columns represent the variables and the rows represent the time points. See the -example data [ts_data] for the correct format.} +example data \code{\link{ts_data}} for the correct format.} -\item{beep}{A vector of beeps with length of `nrow(data)`. The beep indicator +\item{beep}{A vector of beeps with length of \code{nrow(data)}. The beep indicator can be used to remove overnight effects from the last beep of a day to the first beep of the next day. This should be a vector of positive integers. If left empty, the function will assume that there are no overnight @@ -36,45 +38,53 @@ should be a named list, with names corresponding to the parameter names and values corresponding to the prior distributions. The following priors can be specified: \itemize{ -\item `prior_Beta_loc` A matrix of the same dimensions as the beta matrix +\item \code{prior_Beta_loc} A matrix of the same dimensions as the beta matrix `B` containing the mean of the prior distribution for the beta coefficients. -\item `prior_Beta_scale` A matrix of the same dimensions as the beta matrix +\item \code{prior_Beta_scale} A matrix of the same dimensions as the beta matrix `B` containing the standard deviation of the prior distribution for the beta coefficients.}} \item{method}{A string indicating the method to use for fitting the model. Options are "sampling" (for MCMC estimation) or "variational" (for -variational inference). We currently recommend only using MCMC estimation.} +variational inference). We currently recommend only using MCMC estimation. +Default: \code{"sampling"}.} \item{cov_prior}{A string indicating the prior distribution to use for the -covariance matrix. Options are "LKJ" or "IW" (Inverse-Wishart).} +covariance matrix. Options are "LKJ" or "IW" (Inverse-Wishart). Default: \code{"LKJ"}.} \item{rmv_overnight}{A logical indicating whether to remove overnight -effects. Default is `FALSE`. If `TRUE`, the function will remove overnight +effects. Default is \code{FALSE}. If `TRUE`, the function will remove overnight effects from the last beep of a day to the first beep of the next day. -This requires the `beep` argument to be specified.} +This requires the \code{beep} argument to be specified.} \item{iter_sampling}{An integer specifying the number of iterations for the -sampling method. Default is 500.} +sampling method. Default is \code{500}.} \item{iter_warmup}{An integer specifying the number of warmup iterations for -the sampling method. Default is 500.} +the sampling method. Default is \code{500}.} \item{n_chains}{An integer specifying the number of chains for the sampling method. Default is 4. If variational inference is used, the number of -iterations is calculated as `iter_sampling`*`n_chains`.} +iterations is calculated as \code{iter_sampling}*\code{n_chains}.} \item{n_cores}{An integer specifying the number of cores to use for parallel -computation. Default is 1. [rstan] is used for parallel computation.} +computation. Default is \code{1}. \code{rstan} is used for parallel computation.} \item{center_only}{A logical indicating whether to only center (and not -scale) the data. Default is `FALSE`.} +scale) the data. Default is \code{FALSE}.} -\item{...}{Additional arguments passed to the `rstan::sampling` or -`rstan::vb` function.} +\item{return_all}{A logical indicating whether to return all model inputs, including +the data and prior objects. Default is \code{TRUE}.} + +\item{ahead}{An integer specifying the forecast horizon. Default is \code{0}. If `ahead` is +greater than 0, the function will return posterior predictive forecasts for the specified number of time points. +This functionality is experimental and may not work as expected.} + +\item{...}{Additional arguments passed to the \code{\link[rstan:sampling]{rstan::sampling}} or +\code{\link[rstan:vb]{rstan::vb}} function.} } \value{ -A `tsnet_fit` object in list format. The object contains the +A \code{tsnet_fit} object in list format. The object contains the following elements: \item{fit}{A stanfit object containing the fitted model.} \item{arguments}{The number of variables "p", the number of time points "n_t", the column names "cnames", and the arguments used in the function call.} @@ -84,10 +94,11 @@ This function fits a Bayesian GVAR model to the provided data using Stan. The estimation procedure is described further in Siepe, Kloft & Heck (2023) . The current implementation allows for a normal prior on the temporal effects and either an Inverse Wishart or - an LKJ prior on the contemporaneous effects. `rstan` is used as a backend + an LKJ prior on the contemporaneous effects. \code{rstan} is used as a backend for fitting the model in Stan. Data should be provided in long format, where the columns represent the variables and the rows represent the time points. - Data are automatically z-scaled for estimation. + Data are automatically z-scaled for estimation. The model currently does not + support missing data. } \details{ \bold{General Information} @@ -102,33 +113,33 @@ This function fits a Bayesian GVAR model to the provided data lagged effects). This is typically represented in the partial correlation (pcor) matrix. If the model is represented and interpreted as a network, variables are called -*nodes*, *edges* represent the statistical association between the nodes, and -*edge weights* quantify the strength of these associations. +\emph{nodes}, \emph{edges} represent the statistical association between the nodes, and +\emph{edge weights} quantify the strength of these associations. \bold{Model} Let \eqn{Y} be a matrix with \eqn{n} rows and \eqn{p} columns, where \eqn{n_t} is the number of time points and \eqn{p} is the number of variables. The GVAR model is - given by the following equations: \deqn{Y_t = B* Y_{t-1} + zeta_t} - \deqn{zeta_t \sim N(0, \Sigma)} where \eqn{B} is a `p x p` matrix of VAR - coefficients between variables i and j (beta_ij), \eqn{\zeta_t} contains the - innovations at time point \eqn{t}, and \eqn{\Sigma} is a `p x p`covariance matrix. + given by the following equations: \deqn{Y_t = B* Y_{t-1} + \zeta_t} + \deqn{\zeta_t \sim N(0, \Sigma)} where \eqn{B} is a \eqn{p \times p} matrix of VAR + coefficients between variables \eqn{i} and \eqn{j} (\eqn{\beta_{ij}}), \eqn{\zeta_t} contains the + innovations at time point \eqn{t}, and \eqn{\Sigma} is a \eqn{p \times p} covariance matrix. The inverse of \eqn{\Sigma} is the precision matrix, which is used to obtain the - partial correlations between variables (rho_ij). The model setup is + partial correlations between variables (\eqn{\rho_{ij}}). The model setup is explained in more detail in Siepe, Kloft & Heck (2023) . \bold{Prior Setup} - For the p x p temporal matrix B (containing the beta coefficients), we use + For the \eqn{p \times p} temporal matrix B (containing the \eqn{\beta} coefficients), we use a normal prior distribution on each individual parameter: \deqn{\beta_{ij} - \sim N(PriorBetaLoc_{ij}, PriorBetaScale_{ij})} where `PriorBetaLoc` is the - mean of the prior distribution and `PriorBetaScale` is the standard + \sim N(PriorBetaLoc_{ij}, PriorBetaScale_{ij})} where \code{PriorBetaLoc} is the + mean of the prior distribution and \code{PriorBetaScale} is the standard deviation of the prior distribution. The default prior is a weakly informative normal distribution with mean 0 and standard deviation 0.5. The user can specify a different prior distribution by a matrix - `prior_Beta_loc` and a matrix `prior_Beta_scale` with the same dimensions + \code{prior_Beta_loc} and a matrix \code{prior_Beta_scale} with the same dimensions as \eqn{B}. @@ -147,49 +158,50 @@ This function fits a Bayesian GVAR model to the provided data an additional beta prior on each of the partial correlations. We first assigned an uninformed LKJ prior to the Cholesky factor decomposition of the correlation matrix of innovations: \deqn{\Omega_L \sim - LKJ-Cholesky(\eta)}. For \eqn{\eta = 1}, this implies a symmetric marginal + LKJ-Cholesky(\eta)} + For \eqn{\eta = 1}, this implies a symmetric marginal scaled beta distribution on the zero-order correlations \eqn{\omega_{ij}}. \deqn{(\omega_{ij}+1)/2 \sim Beta(p/2, p/2)} We can then obtain the covariance matrix and, - subsequently, the precision matrix (see Siepe, Kloft & Heck (2023)) - for details. + subsequently, the precision matrix (see Siepe, Kloft & Heck, 2023, + for details). The second part of the prior is a beta prior on each partial correlation \eqn{\rho_{ij}} (obtained from the off-diagonal elements of the precision matrix). This prior was assigned by transforming the partial correlations to the - interval of 0,1 and then assigning a proportional (mean-variance + interval of \eqn{[0,1]} and then assigning a proportional (mean-variance parameterized) beta prior: \deqn{(\rho_{ij}+1)/2 \sim Beta_{prop}(PriorRhoLoc, PriorRhoScale)} A beta location parameter of 0.5 translates to an expected correlation of 0. - The variance parameter of sqrt(0.5) implies a uniform distribution of + The variance parameter of \eqn{\sqrt(0.5)} implies a uniform distribution of partial correlations. The user can specify a different prior distribution by a matrix - `prior_Rho_loc` and a matrix `prior_Rho_scale` with the same dimensions as + \code{prior_Rho_loc} and a matrix \code{prior_Rho_scale} with the same dimensions as the partial correlation matrix. Additionally, the user can change \eqn{eta} - via the `prior_Eta` parameter. + via the \code{prior_Eta} parameter. The Inverse-Wishart prior is a distribution on the innovation covariance - matrix `Sigma`: + matrix \eqn{\Sigma}: \deqn{\Sigma \sim IW(\nu, S)} where \eqn{\nu} is the degrees of freedom and \eqn{S} is the scale matrix. We here use the default prior of \deqn{nu = delta + p - 1} for the degrees of freedom, where \eqn{\delta} is defined as \eqn{s_{\rho}^{-1}-1} and \eqn{s_{\rho}} is the standard deviation of the implied marginal beta distribution of the partial correlations. For the scale matrix \eqn{S}, we use the identity matrix - \eqn{I_p} of order p. + \eqn{I_p} of order \eqn{p}. The user can set a prior on the expected standard deviation of the partial - correlations by specifying a `prior_Rho_marginal` parameter. The default + correlations by specifying a \code{prior_Rho_marginal} parameter. The default value is 0.25, which has worked well in a simulation study. - Additionally, the user can specify a `prior_S` parameter to set a different + Additionally, the user can specify a \code{prior_S} parameter to set a different scale matrix. \bold{Sampling} The model can be fitted using either MCMC sampling or variational - inference via [rstan]. Per default, the model is fitted using the Stan + inference via \code{rstan}. Per default, the model is fitted using the Stan Hamiltonian Monte Carlo (HMC) No U-Turn (NUTS) sampler with 4 chains, 500 warmup iterations and 500 sampling iterations. We use a default - target average acceptance probability `adapt_delta` of 0.8. As the output - is returned as a standard `stanfit` object, the user can use the - `rstan` package to extract and analyze the results and obtain convergence + target average acceptance probability \code{adapt_delta} of 0.8. As the output + is returned as a standard \code{stanfit} object, the user can use the + \code{rstan} package to extract and analyze the results and obtain convergence diagnostics. } \examples{ diff --git a/man/ts_data.Rd b/man/ts_data.Rd index 10a8d7d..be8d07c 100644 --- a/man/ts_data.Rd +++ b/man/ts_data.Rd @@ -12,14 +12,14 @@ } } \source{ -Simulated using the [graphicalVAR::graphicalVARsim()] function. +Simulated using the \code{\link[graphicalVAR]{graphicalVARsim}} function in the \code{graphicalVAR} package. } \usage{ data(ts_data) } \description{ This dataset contains a simulated time series dataset for two individuals -generated using the `graphicalVAR` package. The dataset is useful for testing +generated using the \code{graphicalVAR} package. The dataset is useful for testing and demonstrating the functionality of the package. } \details{ diff --git a/src/RcppExports.o b/src/RcppExports.o index a4dfcce..0908d5c 100644 Binary files a/src/RcppExports.o and b/src/RcppExports.o differ diff --git a/src/stanExports_VAR_LKJ.o b/src/stanExports_VAR_LKJ.o index 36eaeb3..ead4e79 100644 Binary files a/src/stanExports_VAR_LKJ.o and b/src/stanExports_VAR_LKJ.o differ diff --git a/src/stanExports_VAR_LKJ_beep.h b/src/stanExports_VAR_LKJ_beep.h index be78124..253cd1a 100644 --- a/src/stanExports_VAR_LKJ_beep.h +++ b/src/stanExports_VAR_LKJ_beep.h @@ -27,61 +27,78 @@ namespace model_VAR_LKJ_beep_namespace { using stan::model::model_base_crtp; using namespace stan::math; stan::math::profile_map profiles__; -static constexpr std::array locations_array__ = +static constexpr std::array locations_array__ = {" (found before start of program)", - " (in 'string', line 23, column 2 to column 23)", - " (in 'string', line 27, column 2 to column 34)", - " (in 'string', line 28, column 2 to column 24)", - " (in 'string', line 33, column 2 to column 67)", - " (in 'string', line 36, column 2 to line 37, column 68)", - " (in 'string', line 39, column 2 to column 18)", - " (in 'string', line 84, column 2 to column 28)", - " (in 'string', line 85, column 2 to column 22)", - " (in 'string', line 42, column 11 to column 12)", - " (in 'string', line 42, column 13 to column 14)", - " (in 'string', line 42, column 4 to column 43)", - " (in 'string', line 48, column 10 to column 23)", - " (in 'string', line 47, column 13 to line 49, column 9)", - " (in 'string', line 46, column 10 to column 65)", - " (in 'string', line 45, column 18 to line 47, column 9)", - " (in 'string', line 45, column 8 to line 49, column 9)", - " (in 'string', line 44, column 19 to line 50, column 7)", - " (in 'string', line 44, column 6 to line 50, column 7)", - " (in 'string', line 43, column 17 to line 51, column 5)", - " (in 'string', line 43, column 4 to line 51, column 5)", - " (in 'string', line 40, column 2 to line 52, column 3)", - " (in 'string', line 88, column 11 to column 12)", - " (in 'string', line 88, column 14 to column 15)", - " (in 'string', line 88, column 4 to column 75)", - " (in 'string', line 91, column 15 to column 16)", - " (in 'string', line 91, column 8 to column 38)", - " (in 'string', line 92, column 8 to column 75)", - " (in 'string', line 90, column 30 to line 93, column 7)", - " (in 'string', line 90, column 6 to line 93, column 7)", - " (in 'string', line 89, column 17 to line 94, column 5)", - " (in 'string', line 89, column 4 to line 94, column 5)", - " (in 'string', line 86, column 2 to line 95, column 3)", - " (in 'string', line 57, column 2 to column 50)", - " (in 'string', line 58, column 2 to column 55)", - " (in 'string', line 59, column 2 to column 49)", - " (in 'string', line 66, column 8 to line 67, column 73)", - " (in 'string', line 63, column 15 to line 68, column 9)", - " (in 'string', line 63, column 6 to line 68, column 9)", - " (in 'string', line 62, column 17 to line 69, column 7)", - " (in 'string', line 62, column 4 to line 69, column 7)", - " (in 'string', line 61, column 15 to line 70, column 5)", - " (in 'string', line 61, column 2 to line 70, column 5)", - " (in 'string', line 73, column 11 to column 12)", - " (in 'string', line 73, column 14 to column 15)", - " (in 'string', line 73, column 4 to column 75)", - " (in 'string', line 76, column 15 to column 16)", - " (in 'string', line 76, column 8 to column 38)", - " (in 'string', line 77, column 8 to column 69)", - " (in 'string', line 75, column 30 to line 78, column 7)", - " (in 'string', line 75, column 6 to line 78, column 7)", - " (in 'string', line 74, column 17 to line 79, column 5)", - " (in 'string', line 74, column 4 to line 79, column 5)", - " (in 'string', line 71, column 2 to line 80, column 3)", + " (in 'string', line 27, column 2 to column 23)", + " (in 'string', line 31, column 2 to column 34)", + " (in 'string', line 32, column 2 to column 24)", + " (in 'string', line 37, column 2 to column 67)", + " (in 'string', line 40, column 2 to line 41, column 68)", + " (in 'string', line 43, column 2 to column 18)", + " (in 'string', line 88, column 2 to column 28)", + " (in 'string', line 89, column 2 to column 28)", + " (in 'string', line 95, column 2 to column 73)", + " (in 'string', line 46, column 11 to column 12)", + " (in 'string', line 46, column 13 to column 14)", + " (in 'string', line 46, column 4 to column 43)", + " (in 'string', line 52, column 10 to column 23)", + " (in 'string', line 51, column 13 to line 53, column 9)", + " (in 'string', line 50, column 10 to column 65)", + " (in 'string', line 49, column 18 to line 51, column 9)", + " (in 'string', line 49, column 8 to line 53, column 9)", + " (in 'string', line 48, column 19 to line 54, column 7)", + " (in 'string', line 48, column 6 to line 54, column 7)", + " (in 'string', line 47, column 17 to line 55, column 5)", + " (in 'string', line 47, column 4 to line 55, column 5)", + " (in 'string', line 44, column 2 to line 56, column 3)", + " (in 'string', line 92, column 4 to column 19)", + " (in 'string', line 91, column 27 to line 93, column 3)", + " (in 'string', line 91, column 2 to line 93, column 3)", + " (in 'string', line 98, column 13 to column 14)", + " (in 'string', line 98, column 6 to column 36)", + " (in 'string', line 100, column 6 to column 71)", + " (in 'string', line 97, column 30 to line 101, column 5)", + " (in 'string', line 97, column 4 to line 101, column 5)", + " (in 'string', line 96, column 17 to line 102, column 3)", + " (in 'string', line 96, column 2 to line 102, column 3)", + " (in 'string', line 105, column 10 to column 15)", + " (in 'string', line 105, column 24 to column 25)", + " (in 'string', line 105, column 4 to column 38)", + " (in 'string', line 106, column 11 to column 12)", + " (in 'string', line 106, column 4 to column 31)", + " (in 'string', line 108, column 13 to column 14)", + " (in 'string', line 108, column 6 to column 38)", + " (in 'string', line 109, column 6 to column 64)", + " (in 'string', line 110, column 6 to column 27)", + " (in 'string', line 111, column 6 to column 32)", + " (in 'string', line 114, column 8 to column 82)", + " (in 'string', line 113, column 32 to line 115, column 5)", + " (in 'string', line 113, column 6 to line 115, column 5)", + " (in 'string', line 107, column 23 to line 116, column 3)", + " (in 'string', line 107, column 4 to line 116, column 3)", + " (in 'string', line 103, column 15 to line 117, column 3)", + " (in 'string', line 103, column 2 to line 117, column 3)", + " (in 'string', line 61, column 2 to column 50)", + " (in 'string', line 62, column 2 to column 55)", + " (in 'string', line 63, column 2 to column 49)", + " (in 'string', line 70, column 8 to line 71, column 73)", + " (in 'string', line 67, column 15 to line 72, column 9)", + " (in 'string', line 67, column 6 to line 72, column 9)", + " (in 'string', line 66, column 17 to line 73, column 7)", + " (in 'string', line 66, column 4 to line 73, column 7)", + " (in 'string', line 65, column 15 to line 74, column 5)", + " (in 'string', line 65, column 2 to line 74, column 5)", + " (in 'string', line 77, column 11 to column 12)", + " (in 'string', line 77, column 14 to column 15)", + " (in 'string', line 77, column 4 to column 75)", + " (in 'string', line 80, column 15 to column 16)", + " (in 'string', line 80, column 8 to column 38)", + " (in 'string', line 81, column 8 to column 69)", + " (in 'string', line 79, column 30 to line 82, column 7)", + " (in 'string', line 79, column 6 to line 82, column 7)", + " (in 'string', line 78, column 17 to line 83, column 5)", + " (in 'string', line 78, column 4 to line 83, column 5)", + " (in 'string', line 75, column 2 to line 84, column 3)", " (in 'string', line 5, column 2 to column 17)", " (in 'string', line 6, column 2 to column 17)", " (in 'string', line 7, column 8 to column 9)", @@ -102,18 +119,25 @@ static constexpr std::array locations_array__ = " (in 'string', line 13, column 11 to column 12)", " (in 'string', line 13, column 2 to column 30)", " (in 'string', line 14, column 2 to column 25)", - " (in 'string', line 18, column 2 to column 29)", - " (in 'string', line 23, column 9 to column 10)", - " (in 'string', line 23, column 11 to column 12)", - " (in 'string', line 27, column 23 to column 24)", - " (in 'string', line 28, column 9 to column 10)", - " (in 'string', line 33, column 9 to column 10)", - " (in 'string', line 33, column 11 to column 12)", - " (in 'string', line 36, column 9 to column 10)", - " (in 'string', line 36, column 11 to column 12)", - " (in 'string', line 39, column 9 to column 10)", - " (in 'string', line 39, column 11 to column 12)", - " (in 'string', line 85, column 9 to column 12)"}; + " (in 'string', line 16, column 2 to column 21)", + " (in 'string', line 17, column 2 to column 44)", + " (in 'string', line 18, column 8 to column 13)", + " (in 'string', line 18, column 22 to column 23)", + " (in 'string', line 18, column 2 to column 34)", + " (in 'string', line 22, column 2 to column 29)", + " (in 'string', line 27, column 9 to column 10)", + " (in 'string', line 27, column 11 to column 12)", + " (in 'string', line 31, column 23 to column 24)", + " (in 'string', line 32, column 9 to column 10)", + " (in 'string', line 37, column 9 to column 10)", + " (in 'string', line 37, column 11 to column 12)", + " (in 'string', line 40, column 9 to column 10)", + " (in 'string', line 40, column 11 to column 12)", + " (in 'string', line 43, column 9 to column 10)", + " (in 'string', line 43, column 11 to column 12)", + " (in 'string', line 89, column 9 to column 18)", + " (in 'string', line 95, column 9 to column 10)", + " (in 'string', line 95, column 12 to column 13)"}; #include class model_VAR_LKJ_beep final : public model_base_crtp { private: @@ -126,6 +150,9 @@ class model_VAR_LKJ_beep final : public model_base_crtp { Eigen::Matrix prior_Rho_loc_data__; Eigen::Matrix prior_Rho_scale_data__; int prior_Eta; + int ahead; + int compute_log_lik; + std::vector> Y_future; int first_beep; int log_lik_1dim__; Eigen::Map> prior_Beta_loc{nullptr, 0, 0}; @@ -153,35 +180,35 @@ class model_VAR_LKJ_beep final : public model_base_crtp { try { int pos__ = std::numeric_limits::min(); pos__ = 1; - current_statement__ = 54; + current_statement__ = 71; context__.validate_dims("data initialization", "K", "int", std::vector{}); K = std::numeric_limits::min(); - current_statement__ = 54; + current_statement__ = 71; K = context__.vals_i("K")[(1 - 1)]; - current_statement__ = 54; + current_statement__ = 71; stan::math::check_greater_or_equal(function__, "K", K, 0); - current_statement__ = 55; + current_statement__ = 72; context__.validate_dims("data initialization", "T", "int", std::vector{}); T = std::numeric_limits::min(); - current_statement__ = 55; + current_statement__ = 72; T = context__.vals_i("T")[(1 - 1)]; - current_statement__ = 55; + current_statement__ = 72; stan::math::check_greater_or_equal(function__, "T", T, 0); - current_statement__ = 56; + current_statement__ = 73; stan::math::validate_non_negative_index("beep", "T", T); - current_statement__ = 57; + current_statement__ = 74; context__.validate_dims("data initialization", "beep", "int", std::vector{static_cast(T)}); beep = std::vector(T, std::numeric_limits::min()); - current_statement__ = 57; + current_statement__ = 74; beep = context__.vals_i("beep"); - current_statement__ = 58; + current_statement__ = 75; stan::math::validate_non_negative_index("Y", "T", T); - current_statement__ = 59; + current_statement__ = 76; stan::math::validate_non_negative_index("Y", "K", K); - current_statement__ = 60; + current_statement__ = 77; context__.validate_dims("data initialization", "Y", "double", std::vector{static_cast(T), static_cast(K)}); Y = std::vector>(T, @@ -189,28 +216,28 @@ class model_VAR_LKJ_beep final : public model_base_crtp { std::numeric_limits::quiet_NaN())); { std::vector Y_flat__; - current_statement__ = 60; + current_statement__ = 77; Y_flat__ = context__.vals_r("Y"); - current_statement__ = 60; + current_statement__ = 77; pos__ = 1; - current_statement__ = 60; + current_statement__ = 77; for (int sym1__ = 1; sym1__ <= K; ++sym1__) { - current_statement__ = 60; + current_statement__ = 77; for (int sym2__ = 1; sym2__ <= T; ++sym2__) { - current_statement__ = 60; + current_statement__ = 77; stan::model::assign(Y, Y_flat__[(pos__ - 1)], "assigning variable Y", stan::model::index_uni(sym2__), stan::model::index_uni(sym1__)); - current_statement__ = 60; + current_statement__ = 77; pos__ = (pos__ + 1); } } } - current_statement__ = 61; + current_statement__ = 78; stan::math::validate_non_negative_index("prior_Beta_loc", "K", K); - current_statement__ = 62; + current_statement__ = 79; stan::math::validate_non_negative_index("prior_Beta_loc", "K", K); - current_statement__ = 63; + current_statement__ = 80; context__.validate_dims("data initialization", "prior_Beta_loc", "double", std::vector{static_cast(K), static_cast(K)}); @@ -221,28 +248,28 @@ class model_VAR_LKJ_beep final : public model_base_crtp { K, K); { std::vector prior_Beta_loc_flat__; - current_statement__ = 63; + current_statement__ = 80; prior_Beta_loc_flat__ = context__.vals_r("prior_Beta_loc"); - current_statement__ = 63; + current_statement__ = 80; pos__ = 1; - current_statement__ = 63; + current_statement__ = 80; for (int sym1__ = 1; sym1__ <= K; ++sym1__) { - current_statement__ = 63; + current_statement__ = 80; for (int sym2__ = 1; sym2__ <= K; ++sym2__) { - current_statement__ = 63; + current_statement__ = 80; stan::model::assign(prior_Beta_loc, prior_Beta_loc_flat__[(pos__ - 1)], "assigning variable prior_Beta_loc", stan::model::index_uni(sym2__), stan::model::index_uni(sym1__)); - current_statement__ = 63; + current_statement__ = 80; pos__ = (pos__ + 1); } } } - current_statement__ = 64; + current_statement__ = 81; stan::math::validate_non_negative_index("prior_Beta_scale", "K", K); - current_statement__ = 65; + current_statement__ = 82; stan::math::validate_non_negative_index("prior_Beta_scale", "K", K); - current_statement__ = 66; + current_statement__ = 83; context__.validate_dims("data initialization", "prior_Beta_scale", "double", std::vector{static_cast(K), static_cast(K)}); @@ -253,29 +280,29 @@ class model_VAR_LKJ_beep final : public model_base_crtp { K, K); { std::vector prior_Beta_scale_flat__; - current_statement__ = 66; + current_statement__ = 83; prior_Beta_scale_flat__ = context__.vals_r("prior_Beta_scale"); - current_statement__ = 66; + current_statement__ = 83; pos__ = 1; - current_statement__ = 66; + current_statement__ = 83; for (int sym1__ = 1; sym1__ <= K; ++sym1__) { - current_statement__ = 66; + current_statement__ = 83; for (int sym2__ = 1; sym2__ <= K; ++sym2__) { - current_statement__ = 66; + current_statement__ = 83; stan::model::assign(prior_Beta_scale, prior_Beta_scale_flat__[(pos__ - 1)], "assigning variable prior_Beta_scale", stan::model::index_uni(sym2__), stan::model::index_uni(sym1__)); - current_statement__ = 66; + current_statement__ = 83; pos__ = (pos__ + 1); } } } - current_statement__ = 67; + current_statement__ = 84; stan::math::validate_non_negative_index("prior_Rho_loc", "K", K); - current_statement__ = 68; + current_statement__ = 85; stan::math::validate_non_negative_index("prior_Rho_loc", "K", K); - current_statement__ = 69; + current_statement__ = 86; context__.validate_dims("data initialization", "prior_Rho_loc", "double", std::vector{static_cast(K), static_cast(K)}); @@ -286,28 +313,28 @@ class model_VAR_LKJ_beep final : public model_base_crtp { K, K); { std::vector prior_Rho_loc_flat__; - current_statement__ = 69; + current_statement__ = 86; prior_Rho_loc_flat__ = context__.vals_r("prior_Rho_loc"); - current_statement__ = 69; + current_statement__ = 86; pos__ = 1; - current_statement__ = 69; + current_statement__ = 86; for (int sym1__ = 1; sym1__ <= K; ++sym1__) { - current_statement__ = 69; + current_statement__ = 86; for (int sym2__ = 1; sym2__ <= K; ++sym2__) { - current_statement__ = 69; + current_statement__ = 86; stan::model::assign(prior_Rho_loc, prior_Rho_loc_flat__[(pos__ - 1)], "assigning variable prior_Rho_loc", stan::model::index_uni(sym2__), stan::model::index_uni(sym1__)); - current_statement__ = 69; + current_statement__ = 86; pos__ = (pos__ + 1); } } } - current_statement__ = 70; + current_statement__ = 87; stan::math::validate_non_negative_index("prior_Rho_scale", "K", K); - current_statement__ = 71; + current_statement__ = 88; stan::math::validate_non_negative_index("prior_Rho_scale", "K", K); - current_statement__ = 72; + current_statement__ = 89; context__.validate_dims("data initialization", "prior_Rho_scale", "double", std::vector{static_cast(K), static_cast(K)}); @@ -318,66 +345,120 @@ class model_VAR_LKJ_beep final : public model_base_crtp { K, K); { std::vector prior_Rho_scale_flat__; - current_statement__ = 72; + current_statement__ = 89; prior_Rho_scale_flat__ = context__.vals_r("prior_Rho_scale"); - current_statement__ = 72; + current_statement__ = 89; pos__ = 1; - current_statement__ = 72; + current_statement__ = 89; for (int sym1__ = 1; sym1__ <= K; ++sym1__) { - current_statement__ = 72; + current_statement__ = 89; for (int sym2__ = 1; sym2__ <= K; ++sym2__) { - current_statement__ = 72; + current_statement__ = 89; stan::model::assign(prior_Rho_scale, prior_Rho_scale_flat__[(pos__ - 1)], "assigning variable prior_Rho_scale", stan::model::index_uni(sym2__), stan::model::index_uni(sym1__)); - current_statement__ = 72; + current_statement__ = 89; pos__ = (pos__ + 1); } } } - current_statement__ = 73; + current_statement__ = 90; context__.validate_dims("data initialization", "prior_Eta", "int", std::vector{}); prior_Eta = std::numeric_limits::min(); - current_statement__ = 73; + current_statement__ = 90; prior_Eta = context__.vals_i("prior_Eta")[(1 - 1)]; - current_statement__ = 73; + current_statement__ = 90; stan::math::check_greater_or_equal(function__, "prior_Eta", prior_Eta, 1); - current_statement__ = 74; + current_statement__ = 91; + context__.validate_dims("data initialization", "ahead", "int", + std::vector{}); + ahead = std::numeric_limits::min(); + current_statement__ = 91; + ahead = context__.vals_i("ahead")[(1 - 1)]; + current_statement__ = 91; + stan::math::check_greater_or_equal(function__, "ahead", ahead, 0); + current_statement__ = 92; + context__.validate_dims("data initialization", "compute_log_lik", + "int", std::vector{}); + compute_log_lik = std::numeric_limits::min(); + current_statement__ = 92; + compute_log_lik = context__.vals_i("compute_log_lik")[(1 - 1)]; + current_statement__ = 92; + stan::math::check_greater_or_equal(function__, "compute_log_lik", + compute_log_lik, 0); + current_statement__ = 92; + stan::math::check_less_or_equal(function__, "compute_log_lik", + compute_log_lik, 1); + current_statement__ = 93; + stan::math::validate_non_negative_index("Y_future", "ahead", ahead); + current_statement__ = 94; + stan::math::validate_non_negative_index("Y_future", "K", K); + current_statement__ = 95; + context__.validate_dims("data initialization", "Y_future", "double", + std::vector{static_cast(ahead), + static_cast(K)}); + Y_future = std::vector>(ahead, + Eigen::Matrix::Constant(K, + std::numeric_limits::quiet_NaN())); + { + std::vector Y_future_flat__; + current_statement__ = 95; + Y_future_flat__ = context__.vals_r("Y_future"); + current_statement__ = 95; + pos__ = 1; + current_statement__ = 95; + for (int sym1__ = 1; sym1__ <= K; ++sym1__) { + current_statement__ = 95; + for (int sym2__ = 1; sym2__ <= ahead; ++sym2__) { + current_statement__ = 95; + stan::model::assign(Y_future, Y_future_flat__[(pos__ - 1)], + "assigning variable Y_future", stan::model::index_uni(sym2__), + stan::model::index_uni(sym1__)); + current_statement__ = 95; + pos__ = (pos__ + 1); + } + } + } + current_statement__ = 96; first_beep = std::numeric_limits::min(); - current_statement__ = 74; + current_statement__ = 96; first_beep = stan::math::min(beep); - current_statement__ = 75; + current_statement__ = 97; stan::math::validate_non_negative_index("Beta_raw", "K", K); - current_statement__ = 76; + current_statement__ = 98; stan::math::validate_non_negative_index("Beta_raw", "K", K); - current_statement__ = 77; + current_statement__ = 99; stan::math::validate_non_negative_index("L_Theta", "K", K); - current_statement__ = 77; + current_statement__ = 99; stan::math::validate_non_negative_index("L_Theta", "K", K); - current_statement__ = 78; + current_statement__ = 100; stan::math::validate_non_negative_index("sigma_theta", "K", K); - current_statement__ = 79; + current_statement__ = 101; stan::math::validate_non_negative_index("Beta", "K", K); - current_statement__ = 80; + current_statement__ = 102; stan::math::validate_non_negative_index("Beta", "K", K); - current_statement__ = 81; + current_statement__ = 103; stan::math::validate_non_negative_index("Sigma", "K", K); - current_statement__ = 82; + current_statement__ = 104; stan::math::validate_non_negative_index("Sigma", "K", K); - current_statement__ = 83; + current_statement__ = 105; stan::math::validate_non_negative_index("Rho", "K", K); - current_statement__ = 84; + current_statement__ = 106; stan::math::validate_non_negative_index("Rho", "K", K); - current_statement__ = 85; + current_statement__ = 107; log_lik_1dim__ = std::numeric_limits::min(); - current_statement__ = 85; - log_lik_1dim__ = (T - 1); - current_statement__ = 85; - stan::math::validate_non_negative_index("log_lik", "T - 1", + current_statement__ = 107; + log_lik_1dim__ = (T + ahead); + current_statement__ = 107; + stan::math::validate_non_negative_index("log_lik", "T + ahead", log_lik_1dim__); + current_statement__ = 108; + stan::math::validate_non_negative_index("Sigma_chol", "K", K); + current_statement__ = 109; + stan::math::validate_non_negative_index("Sigma_chol", "K", K); } catch (const std::exception& e) { stan::lang::rethrow_located(e, locations_array__[current_statement__]); } @@ -443,22 +524,22 @@ class model_VAR_LKJ_beep final : public model_base_crtp { Eigen::Matrix Rho = Eigen::Matrix::Constant(K, K, DUMMY_VAR__); { - current_statement__ = 9; - stan::math::validate_non_negative_index("Theta", "K", K); current_statement__ = 10; stan::math::validate_non_negative_index("Theta", "K", K); + current_statement__ = 11; + stan::math::validate_non_negative_index("Theta", "K", K); Eigen::Matrix Theta = Eigen::Matrix::Constant(K, K, DUMMY_VAR__); - current_statement__ = 11; + current_statement__ = 12; stan::model::assign(Theta, stan::math::inverse_spd(Sigma), "assigning variable Theta"); - current_statement__ = 20; + current_statement__ = 21; for (int i = 1; i <= K; ++i) { - current_statement__ = 18; + current_statement__ = 19; for (int j = 1; j <= K; ++j) { - current_statement__ = 16; + current_statement__ = 17; if (stan::math::logical_neq(i, j)) { - current_statement__ = 14; + current_statement__ = 15; stan::model::assign(Rho, (-stan::model::rvalue(Theta, "Theta", stan::model::index_uni(i), stan::model::index_uni(j)) / @@ -470,7 +551,7 @@ class model_VAR_LKJ_beep final : public model_base_crtp { "assigning variable Rho", stan::model::index_uni(i), stan::model::index_uni(j)); } else { - current_statement__ = 12; + current_statement__ = 13; stan::model::assign(Rho, 0, "assigning variable Rho", stan::model::index_uni(i), stan::model::index_uni(j)); } @@ -478,21 +559,21 @@ class model_VAR_LKJ_beep final : public model_base_crtp { } } { - current_statement__ = 33; + current_statement__ = 50; lp_accum__.add(stan::math::std_normal_lpdf( stan::math::to_vector(Beta_raw))); - current_statement__ = 34; + current_statement__ = 51; lp_accum__.add(stan::math::lkj_corr_cholesky_lpdf(L_Theta, prior_Eta)); - current_statement__ = 35; + current_statement__ = 52; lp_accum__.add(stan::math::student_t_lpdf(sigma_theta, 3, 0, 2)); - current_statement__ = 42; + current_statement__ = 59; for (int i = 1; i <= K; ++i) { - current_statement__ = 40; + current_statement__ = 57; for (int j = 1; j <= K; ++j) { - current_statement__ = 38; + current_statement__ = 55; if (stan::math::logical_lt(i, j)) { - current_statement__ = 36; + current_statement__ = 53; lp_accum__.add(stan::math::beta_proportion_lpdf( ((stan::model::rvalue(Rho, "Rho", stan::model::index_uni(i), @@ -508,35 +589,35 @@ class model_VAR_LKJ_beep final : public model_base_crtp { } } { - current_statement__ = 43; + current_statement__ = 60; stan::math::validate_non_negative_index("Sigma_chol", "K", K); - current_statement__ = 44; + current_statement__ = 61; stan::math::validate_non_negative_index("Sigma_chol", "K", K); Eigen::Matrix Sigma_chol = Eigen::Matrix::Constant(K, K, DUMMY_VAR__); - current_statement__ = 45; + current_statement__ = 62; stan::model::assign(Sigma_chol, stan::math::diag_pre_multiply(stan::math::exp(sigma_theta), L_Theta), "assigning variable Sigma_chol"); - current_statement__ = 52; + current_statement__ = 69; for (int t = 2; t <= T; ++t) { - current_statement__ = 50; + current_statement__ = 67; if (stan::math::logical_gt( stan::model::rvalue(beep, "beep", stan::model::index_uni(t)), first_beep)) { - current_statement__ = 46; + current_statement__ = 63; stan::math::validate_non_negative_index("mu", "K", K); Eigen::Matrix mu = Eigen::Matrix::Constant(K, DUMMY_VAR__); - current_statement__ = 47; + current_statement__ = 64; stan::model::assign(mu, stan::math::multiply(Beta, stan::model::rvalue(Y, "Y", stan::model::index_uni((t - 1)), stan::model::index_omni())), "assigning variable mu"); - current_statement__ = 48; + current_statement__ = 65; lp_accum__.add(stan::math::multi_normal_cholesky_lpdf( stan::model::rvalue(Y, "Y", stan::model::index_uni(t), @@ -629,23 +710,23 @@ class model_VAR_LKJ_beep final : public model_base_crtp { stan::math::diag_pre_multiply(stan::math::exp(sigma_theta), L_Theta))), "assigning variable Sigma"); { - current_statement__ = 9; - stan::math::validate_non_negative_index("Theta", "K", K); current_statement__ = 10; stan::math::validate_non_negative_index("Theta", "K", K); + current_statement__ = 11; + stan::math::validate_non_negative_index("Theta", "K", K); Eigen::Matrix Theta = Eigen::Matrix::Constant(K, K, std::numeric_limits::quiet_NaN()); - current_statement__ = 11; + current_statement__ = 12; stan::model::assign(Theta, stan::math::inverse_spd(Sigma), "assigning variable Theta"); - current_statement__ = 20; + current_statement__ = 21; for (int i = 1; i <= K; ++i) { - current_statement__ = 18; + current_statement__ = 19; for (int j = 1; j <= K; ++j) { - current_statement__ = 16; + current_statement__ = 17; if (stan::math::logical_neq(i, j)) { - current_statement__ = 14; + current_statement__ = 15; stan::model::assign(Rho, (-stan::model::rvalue(Theta, "Theta", stan::model::index_uni(i), stan::model::index_uni(j)) / @@ -657,7 +738,7 @@ class model_VAR_LKJ_beep final : public model_base_crtp { "assigning variable Rho", stan::model::index_uni(i), stan::model::index_uni(j)); } else { - current_statement__ = 12; + current_statement__ = 13; stan::model::assign(Rho, 0, "assigning variable Rho", stan::model::index_uni(i), stan::model::index_uni(j)); } @@ -678,45 +759,101 @@ class model_VAR_LKJ_beep final : public model_base_crtp { Eigen::Matrix log_lik = Eigen::Matrix::Constant(log_lik_1dim__, std::numeric_limits::quiet_NaN()); - { - current_statement__ = 22; - stan::math::validate_non_negative_index("Sigma_chol", "K", K); + current_statement__ = 25; + for (int t = 1; t <= (T + ahead); ++t) { current_statement__ = 23; - stan::math::validate_non_negative_index("Sigma_chol", "K", K); - Eigen::Matrix Sigma_chol = - Eigen::Matrix::Constant(K, K, + stan::model::assign(log_lik, 0, "assigning variable log_lik", + stan::model::index_uni(t)); + } + Eigen::Matrix Sigma_chol = + Eigen::Matrix::Constant(K, K, + std::numeric_limits::quiet_NaN()); + current_statement__ = 9; + stan::model::assign(Sigma_chol, + stan::math::diag_pre_multiply(stan::math::exp(sigma_theta), L_Theta), + "assigning variable Sigma_chol"); + current_statement__ = 32; + for (int t = 2; t <= T; ++t) { + current_statement__ = 30; + if (stan::math::logical_gt( + stan::model::rvalue(beep, "beep", stan::model::index_uni(t)), + first_beep)) { + current_statement__ = 26; + stan::math::validate_non_negative_index("mu", "K", K); + Eigen::Matrix mu = + Eigen::Matrix::Constant(K, + std::numeric_limits::quiet_NaN()); + current_statement__ = 27; + stan::model::assign(mu, + stan::math::multiply(Beta, + stan::model::rvalue(Y, "Y", stan::model::index_uni((t - 1)), + stan::model::index_omni())), "assigning variable mu"); + current_statement__ = 28; + stan::model::assign(log_lik, + stan::math::multi_normal_cholesky_lpdf( + stan::model::rvalue(Y, "Y", stan::model::index_uni(t), + stan::model::index_omni()), mu, Sigma_chol), + "assigning variable log_lik", stan::model::index_uni(t)); + } + } + current_statement__ = 49; + if (stan::math::logical_gt(ahead, 0)) { + current_statement__ = 33; + stan::math::validate_non_negative_index("Y_forecast", "ahead", ahead); + current_statement__ = 34; + stan::math::validate_non_negative_index("Y_forecast", "K", K); + std::vector> Y_forecast = + std::vector>(ahead, + Eigen::Matrix::Constant(K, + std::numeric_limits::quiet_NaN())); + current_statement__ = 36; + stan::math::validate_non_negative_index("current_Y", "K", K); + Eigen::Matrix current_Y = + Eigen::Matrix::Constant(K, std::numeric_limits::quiet_NaN()); - current_statement__ = 24; - stan::model::assign(Sigma_chol, - stan::math::diag_pre_multiply(stan::math::exp(sigma_theta), L_Theta), - "assigning variable Sigma_chol"); - current_statement__ = 31; - for (int t = 2; t <= T; ++t) { - current_statement__ = 29; - if (stan::math::logical_gt( - stan::model::rvalue(beep, "beep", stan::model::index_uni(t)), - first_beep)) { - current_statement__ = 25; - stan::math::validate_non_negative_index("mu", "K", K); - Eigen::Matrix mu = - Eigen::Matrix::Constant(K, - std::numeric_limits::quiet_NaN()); - current_statement__ = 26; - stan::model::assign(mu, - stan::math::multiply(Beta, - stan::model::rvalue(Y, "Y", stan::model::index_uni((t - 1)), - stan::model::index_omni())), "assigning variable mu"); - current_statement__ = 27; + current_statement__ = 37; + stan::model::assign(current_Y, + stan::model::rvalue(Y, "Y", stan::model::index_uni(T)), + "assigning variable current_Y"); + current_statement__ = 47; + for (int s = 1; s <= ahead; ++s) { + current_statement__ = 38; + stan::math::validate_non_negative_index("mu", "K", K); + Eigen::Matrix mu = + Eigen::Matrix::Constant(K, + std::numeric_limits::quiet_NaN()); + current_statement__ = 39; + stan::model::assign(mu, stan::math::multiply(Beta, current_Y), + "assigning variable mu"); + current_statement__ = 40; + stan::model::assign(Y_forecast, + stan::math::multi_normal_cholesky_rng(mu, Sigma_chol, base_rng__), + "assigning variable Y_forecast", stan::model::index_uni(s)); + current_statement__ = 41; + if (pstream__) { + stan::math::stan_print(pstream__, + stan::model::rvalue(Y_forecast, "Y_forecast", + stan::model::index_uni(s))); + *(pstream__) << std::endl; + } + current_statement__ = 42; + stan::model::assign(current_Y, + stan::model::rvalue(Y_forecast, "Y_forecast", + stan::model::index_uni(s)), "assigning variable current_Y"); + current_statement__ = 45; + if (stan::math::logical_eq(compute_log_lik, 1)) { + current_statement__ = 43; stan::model::assign(log_lik, stan::math::multi_normal_cholesky_lpdf( - stan::model::rvalue(Y, "Y", stan::model::index_uni(t), - stan::model::index_omni()), mu, Sigma_chol), - "assigning variable log_lik", stan::model::index_uni((t - 1))); + stan::model::rvalue(Y_future, "Y_future", + stan::model::index_uni(s)), mu, Sigma_chol), + "assigning variable log_lik", stan::model::index_uni((T + s))); } } } out__.write(min_beep); out__.write(log_lik); + out__.write(Sigma_chol); } catch (const std::exception& e) { stan::lang::rethrow_located(e, locations_array__[current_statement__]); } @@ -863,7 +1000,7 @@ class model_VAR_LKJ_beep final : public model_base_crtp { names__.insert(names__.end(), temp.begin(), temp.end()); } if (emit_generated_quantities__) { - std::vector temp{"min_beep", "log_lik"}; + std::vector temp{"min_beep", "log_lik", "Sigma_chol"}; names__.reserve(names__.size() + temp.size()); names__.insert(names__.end(), temp.begin(), temp.end()); } @@ -892,7 +1029,9 @@ class model_VAR_LKJ_beep final : public model_base_crtp { if (emit_generated_quantities__) { std::vector> temp{std::vector{}, - std::vector{static_cast(log_lik_1dim__)}}; + std::vector{static_cast(log_lik_1dim__)}, + std::vector{static_cast(K), + static_cast(K)}}; dimss__.reserve(dimss__.size() + temp.size()); dimss__.insert(dimss__.end(), temp.begin(), temp.end()); } @@ -943,6 +1082,12 @@ class model_VAR_LKJ_beep final : public model_base_crtp { param_names__.emplace_back(std::string() + "log_lik" + '.' + std::to_string(sym1__)); } + for (int sym1__ = 1; sym1__ <= K; ++sym1__) { + for (int sym2__ = 1; sym2__ <= K; ++sym2__) { + param_names__.emplace_back(std::string() + "Sigma_chol" + '.' + + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + } + } } } inline void @@ -989,13 +1134,19 @@ class model_VAR_LKJ_beep final : public model_base_crtp { param_names__.emplace_back(std::string() + "log_lik" + '.' + std::to_string(sym1__)); } + for (int sym1__ = 1; sym1__ <= K; ++sym1__) { + for (int sym2__ = 1; sym2__ <= K; ++sym2__) { + param_names__.emplace_back(std::string() + "Sigma_chol" + '.' + + std::to_string(sym2__) + '.' + std::to_string(sym1__)); + } + } } } inline std::string get_constrained_sizedtypes() const { - return std::string("[{\"name\":\"Beta_raw\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"parameters\"},{\"name\":\"L_Theta\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"parameters\"},{\"name\":\"sigma_theta\",\"type\":{\"name\":\"vector\",\"length\":" + std::to_string(K) + "},\"block\":\"parameters\"},{\"name\":\"Beta\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"transformed_parameters\"},{\"name\":\"Sigma\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"transformed_parameters\"},{\"name\":\"Rho\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"transformed_parameters\"},{\"name\":\"min_beep\",\"type\":{\"name\":\"int\"},\"block\":\"generated_quantities\"},{\"name\":\"log_lik\",\"type\":{\"name\":\"vector\",\"length\":" + std::to_string(log_lik_1dim__) + "},\"block\":\"generated_quantities\"}]"); + return std::string("[{\"name\":\"Beta_raw\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"parameters\"},{\"name\":\"L_Theta\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"parameters\"},{\"name\":\"sigma_theta\",\"type\":{\"name\":\"vector\",\"length\":" + std::to_string(K) + "},\"block\":\"parameters\"},{\"name\":\"Beta\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"transformed_parameters\"},{\"name\":\"Sigma\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"transformed_parameters\"},{\"name\":\"Rho\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"transformed_parameters\"},{\"name\":\"min_beep\",\"type\":{\"name\":\"int\"},\"block\":\"generated_quantities\"},{\"name\":\"log_lik\",\"type\":{\"name\":\"vector\",\"length\":" + std::to_string(log_lik_1dim__) + "},\"block\":\"generated_quantities\"},{\"name\":\"Sigma_chol\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"generated_quantities\"}]"); } inline std::string get_unconstrained_sizedtypes() const { - return std::string("[{\"name\":\"Beta_raw\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"parameters\"},{\"name\":\"L_Theta\",\"type\":{\"name\":\"vector\",\"length\":" + std::to_string(((K * (K - 1)) /2)) + "},\"block\":\"parameters\"},{\"name\":\"sigma_theta\",\"type\":{\"name\":\"vector\",\"length\":" + std::to_string(K) + "},\"block\":\"parameters\"},{\"name\":\"Beta\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"transformed_parameters\"},{\"name\":\"Sigma\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"transformed_parameters\"},{\"name\":\"Rho\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"transformed_parameters\"},{\"name\":\"min_beep\",\"type\":{\"name\":\"int\"},\"block\":\"generated_quantities\"},{\"name\":\"log_lik\",\"type\":{\"name\":\"vector\",\"length\":" + std::to_string(log_lik_1dim__) + "},\"block\":\"generated_quantities\"}]"); + return std::string("[{\"name\":\"Beta_raw\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"parameters\"},{\"name\":\"L_Theta\",\"type\":{\"name\":\"vector\",\"length\":" + std::to_string(((K * (K - 1)) /2)) + "},\"block\":\"parameters\"},{\"name\":\"sigma_theta\",\"type\":{\"name\":\"vector\",\"length\":" + std::to_string(K) + "},\"block\":\"parameters\"},{\"name\":\"Beta\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"transformed_parameters\"},{\"name\":\"Sigma\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"transformed_parameters\"},{\"name\":\"Rho\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"transformed_parameters\"},{\"name\":\"min_beep\",\"type\":{\"name\":\"int\"},\"block\":\"generated_quantities\"},{\"name\":\"log_lik\",\"type\":{\"name\":\"vector\",\"length\":" + std::to_string(log_lik_1dim__) + "},\"block\":\"generated_quantities\"},{\"name\":\"Sigma_chol\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(K) + ",\"cols\":" + std::to_string(K) + "},\"block\":\"generated_quantities\"}]"); } // Begin method overload boilerplate template inline void @@ -1007,8 +1158,8 @@ class model_VAR_LKJ_beep final : public model_base_crtp { const size_t num_params__ = (((K * K) + (K * K)) + K); const size_t num_transformed = emit_transformed_parameters * ((((K * K) + (K * K)) + (K * K))); - const size_t num_gen_quantities = emit_generated_quantities * ((1 + - log_lik_1dim__)); + const size_t num_gen_quantities = emit_generated_quantities * (((1 + + log_lik_1dim__) + (K * K))); const size_t num_to_write = num_params__ + num_transformed + num_gen_quantities; std::vector params_i; @@ -1026,8 +1177,8 @@ class model_VAR_LKJ_beep final : public model_base_crtp { const size_t num_params__ = (((K * K) + (K * K)) + K); const size_t num_transformed = emit_transformed_parameters * ((((K * K) + (K * K)) + (K * K))); - const size_t num_gen_quantities = emit_generated_quantities * ((1 + - log_lik_1dim__)); + const size_t num_gen_quantities = emit_generated_quantities * (((1 + + log_lik_1dim__) + (K * K))); const size_t num_to_write = num_params__ + num_transformed + num_gen_quantities; vars = std::vector(num_to_write, diff --git a/src/stanExports_VAR_LKJ_beep.o b/src/stanExports_VAR_LKJ_beep.o index c055872..781f32b 100644 Binary files a/src/stanExports_VAR_LKJ_beep.o and b/src/stanExports_VAR_LKJ_beep.o differ diff --git a/src/stanExports_VAR_wishart.o b/src/stanExports_VAR_wishart.o index 37a5514..2b2ba80 100644 Binary files a/src/stanExports_VAR_wishart.o and b/src/stanExports_VAR_wishart.o differ diff --git a/src/stanExports_VAR_wishart_beep.o b/src/stanExports_VAR_wishart_beep.o index 13f9976..853822d 100644 Binary files a/src/stanExports_VAR_wishart_beep.o and b/src/stanExports_VAR_wishart_beep.o differ diff --git a/src/tsnet.dll b/src/tsnet.dll index 9a48d67..aebd44b 100644 Binary files a/src/tsnet.dll and b/src/tsnet.dll differ diff --git a/tests/testthat/test-centrality.R b/tests/testthat/test-centrality.R index 2833d22..4feae54 100644 --- a/tests/testthat/test-centrality.R +++ b/tests/testthat/test-centrality.R @@ -48,8 +48,8 @@ test_that("plot_centrality returns expected output", { cent <- get_centrality(fitobj) result <- plot_centrality(cent) - # Check that the result is a ggplot object - expect_type(result, "list") + # check that the result is a ggplot object + expect_true(is_ggplot(result)) }) diff --git a/tsnet.Rproj b/tsnet.Rproj index 9031c9c..810fc87 100644 --- a/tsnet.Rproj +++ b/tsnet.Rproj @@ -1,4 +1,5 @@ Version: 1.0 +ProjectId: 386e5bc4-5cf4-4302-83f9-edb3ecf2018d RestoreWorkspace: Default SaveWorkspace: Default