#' Extract DCPO Results
#'
#' \code{extract_dcpo_results} is a convenience function that extracts 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 row_number
#' @importFrom tidyr pivot_longer
#' @importFrom purrr map_df
#' @importFrom tibble as_tibble rownames_to_column
#' @importFrom stringr str_replace
#'
#' @export
extract_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))
ktcodes <- dat %>%
dplyr::group_by(country) %>%
dplyr::summarize(first_yr = min(year),
last_yr = max(year))
res <- map_df(pars, function(par) {
if (par == "theta") {
rstan::extract(dcpo_output, pars = "theta") %>%
dplyr::first() %>%
purrr::array_branch(1) %>%
purrr::map(function(x) {
tibble::as_tibble(x, .name_repair = "universal") %>%
dplyr::mutate(tt = dplyr::row_number()) %>%
dplyr::left_join(tcodes, by = "tt") %>%
tidyr::pivot_longer(cols = starts_with("..."), names_to = "kk") %>%
dplyr::mutate(year = if_else(tt == 1,
as.integer(year),
as.integer(min(year, na.rm = TRUE) + tt - 1)),
kk = str_replace(kk, "...", "") %>% as.numeric()) %>%
dplyr::left_join(kcodes, by = "kk") %>%
dplyr::left_join(ktcodes, by = "country") %>%
dplyr::filter(year >= first_yr & year <= last_yr) %>%
dplyr::arrange(kk, tt) %>%
dplyr::select(-first_yr, -last_yr)
})
} else if (par == "sigma") {
rstan::extract(dcpo_output, pars = "sigma") %>%
dplyr::first() %>%
purrr::array_branch(1) %>%
purrr::map(function(x) {
tibble::as_tibble(x) %>%
dplyr::mutate(tt = row_number()) %>%
dplyr::left_join(tcodes, by = "tt") %>%
pivot_longer(cols = starts_with("V"), names_to = "kk") %>%
dplyr::mutate(year = if_else(tt == 1,
as.integer(year),
as.integer(min(year, na.rm = TRUE) + tt - 1)),
kk = str_replace(kk, "V", "") %>% as.numeric()) %>%
dplyr::left_join(kcodes, by = "kk") %>%
dplyr::left_join(ktcodes, by = "country") %>%
dplyr::filter(year >= first_yr & year <= last_yr) %>%
dplyr::arrange(kk, tt) %>%
dplyr::select(-first_yr, -last_yr)
})
} 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.