R/tune_temps.R

In duckmayr/bggum: Bayesian Estimation of Generalized Graded Unfolding Model Parameters

#' tune_temperatures
#'
#' Find Optimal Temperatures for the GGUM MCMCMC Sampler
#'
#' Atchadé, Roberts, and Rosenthal (2011) determine the optimal swap-acceptance
#' rate for Metropolis-coupled MCMC and provide an algorithm for building
#' optimal temperature schedules. We implement this algorithm in the context of
#' the GGUM to provide a temperature schedule that should result in
#' approximately 0.234 swap acceptance rate between adjacent chains.
#'
#' Please see the vignette (via \code{vignette("bggum")}) for a full in-depth
#' practical guide to Bayesian estimation of GGUM parameters.
#'
#' @param data An integer matrix giving the response by each
#'   respondent to each item
#' @param n_temps How many temperatures to make?
#' @param temp_tune_iterations How many iterations should the temperature
#'   tuning algorithm run for each temperature? (default is 5000)
#' @param n_draws How many draws should be used to determine each temperature?
#'   (specifying n_draws < temp_tune_iterations will result in an error;
#'   default is 2500).
#' @param K (Optional) A numeric vector with an element for each item giving
#'   the number of options for the item; if not provided, it is generated by
#'   taking the number of unique options observed in the data
#' @param proposal_sds (Optional) A list of length four where is element is a
#'   numeric vector giving standard deviations for the proposals;
#'   the first element should be a numeric vector with a standard deviation
#'   for the proposal for each respondent's theta parameter (the latent trait),
#'   the second a vector with a standard deviation for each item's alpha
#'   (discrimination) parameter, the third a vector with a standard deviation
#'   for each item's delta (location) parameter, and the fourth a vector with
#'   a standard deviation for each item's tau (option threshold) parameters.
#'   If not given, the standard deviations are all set to 1.0 before any
#'   tuning begins.
#' @param sd_tune_iterations A numeric vector of length one; if proposal
#'   standard deviations are not given, this provides the number of
#'   iterations to use to tune the proposals before the temperature finding
#'   algorithm begins (default is 5000)
#' @param theta_prior_params A numeric vector of length two;
#'   the mean and standard deviation of theta parameters' prior distribution
#'   (where the theta parameters have a normal prior; the default is 0 and 1)
#' @param alpha_prior_params A numeric vector of length four;
#'   the two shape parameters and a and b values for alpha parameters' prior
#'   distribution (where the alpha parameters have a four parameter beta prior;
#'   the default is 1.5, 1.5, 0.25, and 4)
#' @param delta_prior_params A numeric vector of length four;
#'   the two shape parameters and a and b values for delta parameters' prior
#'   distribution (where the delta parameters have a four parameter beta prior;
#'   the default is 2, 2, -5, and 5)
#' @param tau_prior_params A numeric vector of length four;
#'   the two shape parameters and a and b values for tau parameters' prior
#'   distribution (where the tau parameters have a four parameter beta prior;
#'   the default is 2, 2, -6, and 6)
#'
#' @return A numeric vector of temperatures
#'
#' @seealso \code{\link{ggumMCMC}}, \code{\link{ggumMC3}}
#'
#' @references Atchadé, Yves F., Gareth O. Roberts, and Jeffrey S. Rosenthal.
#'   2011.  \dQuote{Towards Optimal Scaling of Metropolis-Coupled Markov Chain
#'   Monte Carlo.} \emph{Statistics and Computing} 21(4): 555--68.
#'
#' @examples
#' ## NOTE: This is a toy example just to demonstrate the function, which uses
#' ## a small dataset and an unreasonably low number of sampling interations.
#' ## For a longer practical guide on Bayesian estimation of GGUM parameters,
#' ## please see the vignette ( via vignette("bggum") ).
#' ## We'll simulate data to use for this example:
#' set.seed(123)
#' sim_data <- ggum_simulation(100, 10, 2)
#' ## Now we can tune the temperature schedule:
#' ## (for the purposes of example, we use 100 iterations,
#' ## though in practice you would use much more)
#' temps <- tune_temperatures(data = sim_data$response_matrix, n_temps = 5,
#'                            temp_tune_iterations = 100, n_draws = 50,
#'                            sd_tune_iterations = 100)
#'
#' @export
tune_temperatures <- function(data, n_temps, temp_tune_iterations = 5000,
                              n_draws = 2500, K = NULL, proposal_sds = NULL,
                              sd_tune_iterations = 5000,
                              theta_prior_params = c(0.0, 1.0),
                              alpha_prior_params = c(1.5, 1.5, 0.25, 4.0),
                              delta_prior_params = c(2.0, 2.0, -5.0, 5.0),
                              tau_prior_params = c(2.0, 2.0, -6.0, 6.0)) {
    n <- nrow(data)
    m <- ncol(data)
    K <- integer(m)
    for ( j in 1:m ) {
        K[j] = length(unique(na.omit(data[ , j])))
    }
    if ( is.null(proposal_sds) ) {
        theta_init <- init_thetas(n, theta_prior_params[1],
                                  theta_prior_params[2])
        alpha_init <- init_alphas(m, alpha_prior_params[1],
                                  alpha_prior_params[2],
                                  alpha_prior_params[3],
                                  alpha_prior_params[4])
        delta_init <- init_deltas(m, delta_prior_params[1],
                                  delta_prior_params[2],
                                  delta_prior_params[3],
                                  delta_prior_params[4])
        tau_init <- init_taus(m, tau_prior_params[1], tau_prior_params[2],
                              tau_prior_params[3], tau_prior_params[4], K)
        proposal_sds <- tune_proposals(data, sd_tune_iterations, K, theta_init,
                                       alpha_init, delta_init, tau_init,
                                       theta_prior_params, alpha_prior_params,
                                       delta_prior_params, tau_prior_params)
    }
    return(.tune_temperatures(data, n_temps, temp_tune_iterations, n_draws,
                              n, m, K, proposal_sds,
                              theta_prior_params[1], theta_prior_params[2],
                              alpha_prior_params[1], alpha_prior_params[2],
                              alpha_prior_params[3], alpha_prior_params[4],
                              delta_prior_params[1], delta_prior_params[2],
                              delta_prior_params[3], delta_prior_params[4],
                              tau_prior_params[1], tau_prior_params[2],
                              tau_prior_params[3], tau_prior_params[4]))
}