In snsf-data/ERforResearch: Expected Rank for Research Evaluation

knitr::opts_chunk$set(echo = TRUE)
options(dplyr.summarise.inform = FALSE)

library(ERforResearch)
library(tidyverse)
library(here)
library(ggalluvial)

The complex methodology of the ER is rather straightforward to use, even without knowledge in building Bayesian hierarchical models. First, for illustration, we will load a mock data set in R. The function get_mock_data() does exactly that.

# Load the mock data
data <- get_mock_data()

The dataset comprises r data %>% pull(panel) %>% unique() %>% length() different panels, with r data %>% filter(panel == "p1") %>% pull(assessor) %>% unique() %>% length() voter and r data %>% filter(panel == "p1") %>% pull(proposal) %>% unique() %>% length() proposals each. Every one of these voters gives a score to each of the proposals. For the first illustrations, we will only use one panel (p1).

data_panel1 <- data %>% 
  filter(panel == "p1")

The aim of this package is to be able to rank the proposals in a certain call or panel. The most straightforward way to do this, is to use the average of the scores a proposal received from all the voters:

# First ranking based on simple means: 
data_panel1 %>% 
  group_by(proposal) %>% 
  summarise(avg = mean(num_grade)) %>% 
  arrange(-avg)

Let us assume, we had funding to accept five proposals in panel 1. In this case, we would accept and fund the proposals #1-5, with average scores ranging from 5 to 4.4. The methodology in the ERforResearch package allows us to incorporate the information on the voters behavior (are they more or less strict in giving grades to proposals?). Additionally to that, the Bayesian hierarchical models (BHM) can incorporate random variation in the grading of the proposals, for example, due to conflict of interests or other one proposal is not grades by every voter, \dots (Note that here no voter had a COI). This happens very regularly in the grant peer review process. For this reason we will use a JAGS model via the rjags package to build such a flexible model. The ERforResearch package includes a function to call a default JAGS BHM model:

?ERforResearch::get_default_jags_model

Then, we have the possibility to simply use that default model and sample from it using the get_mcmc_samples() function.

# !Just for illustration purposes we set the rhat threshold a bit higher.
# Make sure to leave it at the default value of 1.01 once the package-
# functionalities are understood.

# Estimate theta's 
mcmc_samples_object <- 
  get_mcmc_samples(
    data = data_panel1, id_proposal = "proposal",
    id_assessor = "assessor",
    grade_variable = "num_grade",
    path_to_jags_model = NULL, 
          # NULL means we use the default model
    seed = 6,
    rhat_threshold = 1.05,
    dont_bind = TRUE)

# It took the model .... to converge:
mcmc_samples_object$conv_status
c(mcmc_samples_object$n_adapt, 
  mcmc_samples_object$n_burnin,
  mcmc_samples_object$n_iter)

# With the latter object we can do some convergence diagnostics, 
bayesplot::mcmc_trace(mcmc_samples_object$samples$mcmc,
                      pars = paste0("proposal_intercept[", 1:15, "]"))
# Instead of the trace plots, we can directly look at the Summary in the object,
# which includes the rhat values (psrf), Monte-Carlo error (MCerr) and the
# effective sample size (SSeff) for each parameter, like so:
# mcmc_samples_object$summary

# If the model in path_to_jags_model is different from the default model with
# different names for the specific parameters, this has to be declared in the 
# arguments (theta_name, tau_name, sigma_name, tau_voter_name, rank_theta_name).

# We can further produce rank histograms
pl1 <- 
  bayesplot::mcmc_rank_hist(mcmc_samples_object$samples$mcmc,
                            pars = paste0("proposal_intercept[", 1:7, "]"))
pl2 <- 
  bayesplot::mcmc_rank_hist(mcmc_samples_object$samples$mcmc,
                            pars = paste0("proposal_intercept[", 8:15, "]"))

pl1

pl2
# For the further analyis we will simply use all three chains
# and therefore bind them here.
mcmc_samples <- coda::mcmc(do.call(rbind, mcmc_samples_object$samples$mcmc))

We do not have to separately sample from the model, as the functions that follow can easily call the get_mcmc_samples() themselves. If the same seed is given to the seed argument in the respective functions the results are reproducible.

Then, we use the plot_rankogram() function to plot the Rankograms of each proposal, which are the probabilities of each rank plotted against all possible ranks.

plot_rankogram(data = data_panel1,
               id_proposal = "proposal",
               id_assessor = "assessor",
               grade_variable = "num_grade",
               mcmc_samples = mcmc_samples_object
                  # if NULL the `get_mcmc_samples function is called from within
)

Next, we can use the same function with the parameter cumulative_rank_prob set to TRUE in order to plot the cumulative ranking probabilities against the possible ranks.

plot_rankogram(data = data_panel1,
               id_proposal = "proposal",
               id_assessor = "assessor",
               grade_variable = "num_grade",
               cumulative_rank_prob = TRUE,
               mcmc_samples = mcmc_samples_object)

The SUCRA is the area under the cumulative ranking probability plot. The higher the SUCRA value, the higher the likelihood of a proposal being in the top rank or one of the top ranks.

SUCRA_results <- get_sucra(data = data_panel1,
                           id_proposal = "proposal",
                           id_assessor = "assessor",
                           grade_variable = "num_grade",
                           mcmc_samples = mcmc_samples_object)
SUCRA_results$sucra %>% round(2)

Finally we will compute the expected ranks (ER), which can directly be linked to the SUCRA. The lower the ER the better the evaluation of the proposal. The best possible ER is 1 for rank 1.

ER_results <- get_er_from_jags(data = data_panel1,
                               id_proposal = "proposal",
                               id_assessor = "assessor",
                               grade_variable = "num_grade",
                               mcmc_samples = mcmc_samples_object)

ER_results$rankings %>% 
  mutate_if(is.numeric, function(x) round(x, 2)) %>% 
  arrange(er)

The latter, can also be plotted using the plotting_er_results function. As can be seen in the figure, the five best proposals are still the same, but the ordering of all proposals changed. The ER actually does a comparative ranking, by comparing every proposal to every other proposal.

plotting_er_results(ER_results, title = "Panel 1", how_many_fundable = 5,
                    draw_funding_line = FALSE)

To better see, whether the funding line (FL) (here: between the fifth and the sixth proposal) goes through a cluster, we use the posterior distribution of the ER. The ERs of the competing proposals are plotted together with their 50% credible intervals (CrI). A provisional FL is defined at the ER of the last fundable (here the fifth) proposal. The proposals with a 50% CrI not crossing the FL are directly accecpted or rejected, depending whether their ER is lower or higher than the FL. The proposals with a 50% CrI crossing the FL will constitute the random selection group. In the case of panel 1 of the mock data four proposals will be randomly selected from a group of five.

plot_er_distributions(mcmc_samples_object, 
                      n_proposals = data_panel1 %>%
                        summarise(n = n_distinct(proposal)) %>% pull(), 
                      number_fundable = 5)

The flexibility of the method of the ER allows us to not only look at the panels individually but also rank all the proposals from the three panels in one large ranking. To do this, we can simply use the default model, but have to give a parameter for id_section.

writeLines(model_with_sections <- 
    "model{
      # Likelihood:
      for (i in 1:n) { # i is not the application but the review
          grade[i] ~ dnorm(mu[i], inv_sigma2)
                # inv_sigma2 is precision (1 / variance)
          mu[i] <- overall_mean + proposal_intercept[num_proposal[i]] +
          assessor_intercept[num_assessor[i]] + panel_intercept[num_panel[i]]
      }
      # Ranks:
      rank_theta[1:n_proposal] <- rank(-proposal_intercept[])
      # Priors:
      for (j in 1:n_proposal){
        proposal_intercept[j] ~ dnorm(0, inv_tau_proposal2)
      }
      for (l in 1:n_assessors){
        assessor_intercept[l] ~ dnorm(nu[l], inv_tau_assessor2)
      }
      for (l in 1:n_assessors){
        nu[l] ~ dnorm(0, 100) # 1/(0.1^2) = 100
      }
      for (l in 1:n_panels){
        panel_intercept[l] ~ dnorm(0, inv_tau_panel2)
      }
      sigma ~ dunif(0.000001, 2)
      inv_sigma2 <- pow(sigma, -2)
      inv_tau_proposal2 <- pow(tau_proposal, -2)
      tau_proposal ~ dunif(0.000001, 2)
      inv_tau_assessor2 <- pow(tau_assessor, -2)
      tau_assessor ~ dunif(0.000001, 2)
      inv_tau_panel2 <- pow(tau_panel, -2)
      tau_panel ~ dunif(0.000001, 2)
    }")
cat(model_with_sections,
    file = here("jags_model_with_sections.txt"))

Then we use the latter model to directly plot the ER:

ER_results_all <- get_er_from_jags(data = data,
                                   id_proposal = "proposal",
                                   id_assessor = "assessor",
                                   id_panel = "panel",
                                   grade_variable = "num_grade",
                                   tau_name_panel = "tau_panel",
                                   n_chains = 4, n_iter = 5000,
                                   n_burnin = 1000, n_adapt = 1000,
                                   rhat_threshold = 1.1)

plotting_er_results(ER_results_all, title = "All three panels",
                    draw_funding_line = FALSE,
                    how_many_fundable = 15)