R/summarize_dcpo_results.R

In DCPO: Dynamic Comparative Public Opinion

#' Extract DCPO Results
#'
#' \code{summarize_dcpo_results} is a convenience function that produces summary statistics of the main parameters of a DCPO stanfit object along with the relevant identifying information (country, year, question, and cutpoint).
#'
#' @param dcpo_input the data frame of survey items and marginals generated by \code{DCPOtools::dcpo_setup} previously passed to \code{DCPO::dcpo} to generate the stanfit object passed as \code{dcpo_output}
#' @param dcpo_output a stanfit object output by \code{DCPO::dcpo}
#' @param pars a character vector of parameter names to be summarized from the \code{DCPO} model: theta (mean public opinion), sigma (polarization in public opinion), alpha (question dispersion), beta (question-cutpoint difficulty), and/or delta (country-specific question bias)
#' @param probs a numeric vector of quantiles of interest; the default is c(.1, .9)
#'
#' @examples
#' \donttest{
#' out1 <- dcpo(demsup_data,
#'              chime = FALSE,
#'              chains = 1,
#'              iter = 300) # 1 chain/300 iterations for example purposes only; use defaults
#'
#' theta_results <- summarize_dcpo_results(dcpo_input = demsup_data,
#'                                         dcpo_output = out1,
#'                                         pars = "theta")
#' }
#'
#' @return a tibble
#'
#' @importFrom rstan summary
#' @importFrom dplyr mutate group_by summarize select first arrange
#' @importFrom purrr map_df
#' @importFrom tibble as_tibble rownames_to_column
#'
#' @export

summarize_dcpo_results <- function(dcpo_input,
                                 dcpo_output,
                                 pars = c("theta", "sigma", "alpha", "beta", "delta"),
                                 probs = c(.1, .9)) {

  question <- country <- year <- parameter <- kk <- tt <- qq <- rr <- NULL

  dat <- dcpo_input$data

  qcodes <- dat %>%
    dplyr::group_by(question) %>%
    dplyr::summarize(qq = first(qq) %>%
                       as.numeric())

  kcodes <- dat %>%
    dplyr::group_by(country) %>%
    dplyr::summarize(kk = first(kk) %>%
                       as.numeric())

  tcodes <- dat %>%
    dplyr::group_by(year) %>%
    dplyr::summarize(tt = first(tt))

  res <- map_df(pars, function(par) {
    if (par == "theta") {
      rstan::summary(dcpo_output, pars = "theta", probs = probs) %>%
        dplyr::first() %>%
        as.data.frame() %>%
        tibble::rownames_to_column("parameter") %>%
        tibble::as_tibble() %>%
        dplyr::mutate(tt = as.numeric(gsub("theta\\[(\\d+),\\d+\\]",
                                           "\\1",
                                           parameter)),
                      kk = as.numeric(gsub("theta\\[\\d+,(\\d+)\\]",
                                           "\\1",
                                           parameter))) %>%
        dplyr::left_join(kcodes, by = "kk") %>%
        dplyr::left_join(tcodes, by = "tt") %>%
        dplyr::arrange(kk, tt)
    } else if (par == "sigma") {
      rstan::summary(dcpo_output, pars = "sigma", probs = probs) %>%
        dplyr::first() %>%
        as.data.frame() %>%
        tibble::rownames_to_column("parameter") %>%
        tibble::as_tibble() %>%
        dplyr::mutate(tt = as.numeric(gsub("sigma\\[(\\d+),\\d+\\]",
                                           "\\1",
                                           parameter)),
                      kk = as.numeric(gsub("sigma\\[\\d+,(\\d+)\\]",
                                           "\\1",
                                           parameter))) %>%
        dplyr::left_join(kcodes, by = "kk") %>%
        dplyr::left_join(tcodes, by = "tt") %>%
        dplyr::arrange(kk, tt)
    } else if (par == "alpha") {
      rstan::summary(dcpo_output, pars = "alpha", probs = probs) %>%
        dplyr::first() %>%
        as.data.frame() %>%
        tibble::rownames_to_column("parameter") %>%
        tibble::as_tibble() %>%
        dplyr::mutate(qq = as.numeric(gsub("alpha\\[(\\d+)]",
                                           "\\1",
                                           parameter))) %>%
        dplyr::left_join(qcodes, by = "qq") %>%
        dplyr::arrange(qq)
    } else if (par == "beta") {
      rstan::summary(dcpo_output, pars = "beta", probs = probs) %>%
        dplyr::first() %>%
        as.data.frame() %>%
        tibble::rownames_to_column("parameter") %>%
        tibble::as_tibble() %>%
        dplyr::mutate(rr = as.numeric(gsub("beta\\[(\\d+),\\d+\\]",
                                           "\\1",
                                           parameter)),
                      qq = as.numeric(gsub("beta\\[(\\d+),\\d+\\]",
                                           "\\1",
                                           parameter)))%>%
        dplyr::left_join(qcodes, by = "qq") %>%
        dplyr::arrange(qq, rr)
    } else if (par == "delta") {
      rstan::summary(dcpo_output, pars = "delta", probs = probs) %>%
        dplyr::first() %>%
        as.data.frame() %>%
        tibble::rownames_to_column("parameter") %>%
        tibble::as_tibble() %>%
        dplyr::mutate(qq = as.numeric(gsub("delta\\[(\\d+),\\d+\\]",
                                           "\\1",
                                           parameter)),
                      kk = as.numeric(gsub("delta\\[\\d+,(\\d+)\\]",
                                           "\\1",
                                           parameter)))%>%
        dplyr::left_join(qcodes, by = "qq") %>%
        dplyr::left_join(kcodes, by = "kk") %>%
        dplyr::arrange(qq)
    }
  })
  return(res)
}