---
title: "Reproducibility Analysis for 'NIH Peer Review: Criterion Scores Completely Account for Racial Disparities in Overall Impact Scores'"
author: "Elena A. Erosheva, Sheridan Grant, Mei-Ching Chen, Mark D. Lindner, Richard K. Nakamura, Carole J. Lee"
date: "April 3, 2020"
output: html_document
editor_options: 
chunk_output_type: console
---

This R Markdown (RMD) document can be used to reproduce the main results in "NIH Peer Review: Criterion Scores Completely Account for Racial Disparities in Overall Impact Scores", by Elena A. Erosheva, Sheridan Grant, Mei-Ching Chen, Mark Lindner, Richard Nakamura, and Carole J. Lee. 

**Before compilation, download the data ("NIH-public-data_Erosheva-et-al.csv"), move it to a folder of your choice and set your working directory to that folder.**

# Environment Setup

```{r setup, echo=FALSE}
knitr::opts_chunk$set(cache = FALSE)
```

Required R packages:
```{r message=FALSE, warning=FALSE}
require(readxl)
require(readr)
require(dplyr)
require(lme4)
require(knitr)
require(kableExtra)
require(Hmisc)
```

```{r, echo=F}
options(digits = 3)
options(tibble.print_max = 15, tibble.print_min = 15)
```

# Data Handling

In this section, we read in the data, recode variables for clarity in the analyses, and define variable subsets and data subsets for use in the subsequent analyses.

Read data into tibble "long":
```{r, message=F}
long <- read_csv("NIH-public-data_Erosheva-et-al.csv")
n <- dim(long)[1]
```

Make sure factor columns are factor-type objects:
```{r}
long_factors <- c('GROUP_ID', 'PI_RACE', 'PI_ID', 'REVIEWER_ID',
                  'APPLICATION_ID', 'IRG', 'ADMIN_ORG', 'SRG',
                  'APPLICATION_TYPE', 'PI_GENDER', 'PI_ETHNICITY',
                  'CAREER_STAGE', 'DEG_CATEGORY', 'INSTITUTION_BIN')
long[long_factors] <- lapply(long[long_factors], as.factor)
```

We reorder the factors to set up baseline categories for convenience. For example, we make "White" the default level of the PI_RACE variable so that the PI_RACE coefficient represents the conditional expected score difference for a PI being black _relative to_ being white.
```{r}
# White is considered the default level of the race variable
long$GROUP_ID <- factor(long$GROUP_ID,
                        levels = c('Matched White',
                                   'Matched Black',
                                   'Random White',
                                   'All Black'))

# "Experienced" is the default career stage
long$CAREER_STAGE <- factor(long$CAREER_STAGE,
                            levels = c('Experienced',
                                       'ESI',
                                       'Non-ESI NI'))

# "Male" is considered the default gender in this study
long$PI_GENDER <- factor(long$PI_GENDER,
                             levels = c('Male', 'Female'))

# "Non-Hispanic" is considered the default ethnicity in this study
long$PI_ETHNICITY <- factor(long$PI_ETHNICITY,
                            levels = c('Non-Hispanic', 'Hispanic/Latino'))

# White is considered the default level of the race variable
long$PI_RACE <- factor(long$PI_RACE,
                       levels = c('White', 'Black'))

# PhD is considered the default degree
long$DEG_CATEGORY <- factor(long$DEG_CATEGORY,
                            levels = c('PHD', 'MD' ,'MD/PHD', 'Others'))
```


Since all SEPs (Special Emphasis Panels) are coded as a single SRG (a single level of the "SRG" variable), we recode this factor to be a concatentation of the IRG and SRG, so that SEPs in different IRGs are identified as different by the model.
```{r, warning=F}
long <- long %>% mutate(SRG = paste(SRG, IRG, sep = '_'))
long$SRG <- as.factor(long$SRG)
```

## Variables for Models

Defining variable subsets for easier specification of models:
```{r}
criteria <- c('SIGNIFICANCE_INIT', 'INVESTIGATOR_INIT', 'INNOVATION_INIT',
              'APPROACH_INIT', 'ENVIRONMENT_INIT')

ID_clusters <- c('PI_ID', 'REVIEWER_ID')

org_clusters <- c('ADMIN_ORG', 'IRG', 'SRG')

matching <- c('APPLICATION_TYPE', 'AMENDED', 'PI_GENDER', 'PI_ETHNICITY',
              'CAREER_STAGE', 'DEG_CATEGORY', 'INSTITUTION_BIN', 'IRG')

# Application- and applicant-specific variables for the public data are all matching variables except for IRG, which is a structural variable 
app_app_vars <- matching[1:7]
```

Here, we define the matched and random subsets of the data by filtering the rows of "long" based on "GROUP_ID". Then, we make the appropriate variables in "d_matched" and "d_random" factor types.
```{r}
# Filter based on GROUP_ID
d_matched <- long[long$GROUP_ID %in% c('Matched White', 'Matched Black'),]
d_random <- long[long$GROUP_ID != 'Matched White',]

# Ensure factor variables have factor data type, remove unnecessary GROUP_ID factor levels after filtering
d_matched$GROUP_ID <- factor(d_matched$GROUP_ID)
d_random$GROUP_ID <- factor(d_random$GROUP_ID)
for (i in 1:(dim(long)[2])) {
  if (is.factor(long[,i][[1]])) {
    d_matched[,i][[1]] <- factor(d_matched[,i][[1]])
    d_random[,i][[1]] <- factor(d_random[,i][[1]])
  }
}
```

# Racial Disparity Models


We now generate the results of the racial disparity analysis as reported in Table S9 of the Reproducibility section of the supplement. These results are comparable to the full-data results in Table 5 of the main paper, discussed in the section titled "Racial Disparity in Preliminary Overall Impact Scores."

Model specifications:
```{r}
# Model 1: only structural variables and PI Race
structural <- paste('IMPACT_INIT ~ PI_RACE',
                    ' + (1|PI_ID)',
                    ' + (1|REVIEWER_ID)',
                    ' + IRG',
                    ' + ADMIN_ORG',
                    ' + (1|SRG)',
                    sep = ' ')

# Model 2: structural variables, application-/applicant-specific covariates, and PI race
structural_app_app <- paste(structural, ' + ',
                            paste(app_app_vars, collapse = ' + '))

# Model 3: structural variables, criterion scores, and PI race
structural_crit <- paste(structural, ' + ',
                         paste(criteria, collapse = ' + '))

# Model 4: structural variables, application-/applicant-specific covariates, criterion scores, and PI race
structural_app_app_crit <- paste(structural, ' + ',
                                 paste(app_app_vars, collapse = ' + '),
                                 ' + ',
                                 paste(criteria, collapse = ' + '))
```

"Bobyqua" optimization algorithm used throughout:
```{r}
optimizer <- lmerControl(optimizer='bobyqa')
```

We utilize a helper function for computing two-sided p-values based on a normal approximation (justified by the approximate model df of 8000+), which is the following:
```{r}
twoSideP <- function(t) {
  return(2 * pnorm(-abs(t)))
}
```

Then we estimate the matched subset racial disparity models, below: 
```{r}
structural_matched <- lmer(as.formula(structural),
                           data = d_matched, 
                           REML = T,
                           control = optimizer)

structural_app_app_matched <- lmer(as.formula(structural_app_app),
                                   data = d_matched,
                                   REML = T,
                                   control = optimizer)

structural_crit_matched <- lmer(as.formula(structural_crit),
                                data = d_matched,
                                REML = T,
                                control = optimizer)

structural_app_app_crit_matched <- lmer(as.formula(structural_app_app_crit),
                                        data = d_matched,
                                        REML = T,
                                        control = optimizer)
```

We also estimate the random subset racial disparity models, given below: 
```{r}
structural_random <- lmer(as.formula(structural),
                           data = d_random, 
                           REML = T,
                           control = optimizer)

structural_app_app_random <- lmer(as.formula(structural_app_app),
                                   data = d_random,
                                   REML = T,
                                   control = optimizer)

structural_crit_random <- lmer(as.formula(structural_crit),
                               data = d_random,
                               REML = T,
                               control = optimizer)

structural_app_app_crit_random <- lmer(as.formula(structural_app_app_crit),
                                       data = d_random,
                                       REML = T,
                                       control = optimizer)
```

```{r, include=F}

tabDisparity_matched <- matrix(nrow = 8, ncol = 4)

# Note that the PI_RACE coefficient is the second coefficient (after the intercept), hence the [2] index. The columns of the model summary coefficients object are [,1]: coefficient estimate; [,2]: standard error; [,3]: t-value.

# Coefficient estimates
tabDisparity_matched[1,] <- c(fixef(structural_matched)[2],
                      fixef(structural_app_app_matched)[2],
                      fixef(structural_crit_matched)[2],
                      fixef(structural_app_app_crit_matched)[2])

# Standard error estimates
tabDisparity_matched[2,] <- c(summary(structural_matched)$coefficients[2,2],
                      summary(structural_app_app_matched)$coefficients[2,2],
                      summary(structural_crit_matched)$coefficients[2,2],
                      summary(structural_app_app_crit_matched)$coefficients[2,2])

# P-values
tabDisparity_matched[3,] <- c(twoSideP(summary(structural_matched)$coefficients[2,3]),
                      twoSideP(summary(structural_app_app_matched)$coefficients[2,3]),
                      twoSideP(summary(structural_crit_matched)$coefficients[2,3]),
                      twoSideP(summary(structural_app_app_crit_matched)$coefficients[2,3]))

# Random effect standard deviations
tabDisparity_matched[5:8,] <- cbind(data.frame(VarCorr(structural_matched))$sdcor,
                            data.frame(VarCorr(structural_app_app_matched))$sdcor,
                            data.frame(VarCorr(structural_crit_matched))$sdcor,
                            data.frame(VarCorr(structural_app_app_crit_matched))$sdcor)
```

Effect sizes are computed by dividing coefficient estimates by residual standard deviations, per one of _multiple valid options for mixed-effects models_ as detailed by Hedges (2007). We encourage readers to focus on coefficient estimates themselves (first line of table) rather than effect sizes, since many different effect sizes can be computed for mixed-effects models, each with their own subtle interpretations.
```{r, include=F}
# Effect size is coefficient estimate divided by residual standard deviation
tabDisparity_matched[4,] <- tabDisparity_matched[1,] / tabDisparity_matched[8,]
```

```{r, include=F}
tabDisparity_random <- matrix(nrow = 8, ncol = 4)

# Coefficient estimates
tabDisparity_random[1,] <- c(fixef(structural_random)[2],
                     fixef(structural_app_app_random)[2],
                     fixef(structural_crit_random)[2],
                     fixef(structural_app_app_crit_random)[2])

# Standard error estimates
tabDisparity_random[2,] <- c(summary(structural_random)$coefficients[2,2],
                      summary(structural_app_app_random)$coefficients[2,2],
                      summary(structural_crit_random)$coefficients[2,2],
                      summary(structural_app_app_crit_random)$coefficients[2,2])

# P-values
tabDisparity_random[3,] <- c(twoSideP(summary(structural_random)$coefficients[2,3]),
                      twoSideP(summary(structural_app_app_random)$coefficients[2,3]),
                      twoSideP(summary(structural_crit_random)$coefficients[2,3]),
                      twoSideP(summary(structural_app_app_crit_random)$coefficients[2,3]))

# Random effect standard deviations
tabDisparity_random[5:8,] <- cbind(data.frame(VarCorr(structural_random))$sdcor,
                            data.frame(VarCorr(structural_app_app_random))$sdcor,
                            data.frame(VarCorr(structural_crit_random))$sdcor,
                            data.frame(VarCorr(structural_app_app_crit_random))$sdcor)

# Effect sizes
tabDisparity_random[4,] <- tabDisparity_random[1,] / tabDisparity_random[8,]
```

```{r, echo=F}
# Formatting tabDisparity with row/column names
tabDisparity_rownames <- c('Race Coefficient', 'Standard Error', 'P Value',
                   'Effect Size',
                   'Reviewer Intercept Std. Dev.', 
                   'PI Intercept Std. Dev.',
                   'SRG Intercept Std. Dev.',
                   'Residual Std. Dev.')

disparity_models <- c('Structural', 
                      'Structural + Applicant/Application-Specific',
                      'Structural + Criteria',
                      'Structural + Applicant/Application-Specific + Criteria')

rownames(tabDisparity_matched) <- tabDisparity_rownames
colnames(tabDisparity_matched) <- disparity_models

rownames(tabDisparity_random) <- tabDisparity_rownames
colnames(tabDisparity_random) <- disparity_models
```

**Table S9, racial disparity models: matched subset analysis, public-use data.**
```{r, echo=F}
tabDisparity_matched %>%
  kable('html', row.names = T, digits = 3) %>%
  kable_styling(full_width = F) %>%
  column_spec(1, border_right = T)
```
Race coefficient estimates, their effect sizes, and variance components estimates from four hierarchical linear models for preliminary overall impact scores fit on n = 7471 reviews of 2566 applications (matched subset). Model 1 controls for structural covariates; Model 2 controls for structural and matching covariates; Model 3 controls for structural covariates and criterion scores; Model 4 controls for structural, matching covariates, and criterion scores. Coefficient estimates for control variables are not shown.

**Table S9, racial disparity models: random subset analysis, public-use data.**
```{r, echo=F}
tabDisparity_random %>%
  kable('html', row.names = T, digits = 3) %>%
  kable_styling(full_width = F) %>%
  column_spec(1, border_right = T)
```
Race coefficient estimates, their effect sizes, and variance components estimates from four hierarchical linear models for preliminary overall impact scores fit on n = 8595 reviews of 3045 applications (random subset). Model 1 controls for structural covariates; Model 2 controls for structural and matching covariates; Model 3 controls for structural covariates and criterion scores; Model 4 controls for structural, matching covariates, and criterion scores. Coefficient estimates for control variables are not shown.

# Commensuration Bias Models

Next, we generate the results of the commensuration analysis as reported in Table S10 and Figures S3 and S4 of the Reproducibility section of the supplement. These results are comparable to the full-data results in Table 6 of the main paper, discussed in the section titled "Commensuration Model for Preliminary Overall Impact Scores" and elaborated on in the supplement's "Commensuration Practices" section. "NA" in tables refers to random effect standard deviation estimates not having standard errors or p-values.

```{r}
# Commensuration model includes structural variables, application-/applicant-specific covariates, criterion scores, PI race, and PI race-criterion score interactions
commensuration <- paste(structural_app_app_crit,
                        ' + (',
                        paste(criteria, collapse = ' + '), ') * PI_RACE')

# Model fitting
commensuration_matched <- lmer(as.formula(commensuration),
                               data = d_matched,
                               REML = T,
                               control = optimizer)

commensuration_random <- lmer(as.formula(commensuration),
                              data = d_random,
                              REML = T,
                              control = optimizer)
```

```{r, include=F}
tabCommens_matched <- matrix(ncol = 3, nrow = 15)
p_matched <- length(fixef(commensuration_matched))

coefs_matched <- summary(commensuration_matched)$coefficients

# Coefficient estimates and standard errors
tabCommens_matched[1:5,1:2] <- coefs_matched[(p_matched-9):(p_matched-5),1:2]
tabCommens_matched[6,1:2] <- coefs_matched[2,1:2]
tabCommens_matched[7:11,1:2] <- coefs_matched[(p_matched-4):p_matched,1:2]

# P-values
tabCommens_matched[1:5,3] <- sapply(coefs_matched[(p_matched-9):(p_matched-5),3],
                              twoSideP)
tabCommens_matched[6,3] <- sapply(coefs_matched[2,3],
                            twoSideP)
tabCommens_matched[7:11,3] <- sapply(coefs_matched[(p_matched-4):p_matched,3],
                               twoSideP)

# Random effect standard deviations
tabCommens_matched[12:15,1] <- data.frame(VarCorr(commensuration_matched))$sdcor
```

```{r, include=F}
tabCommens_random <- matrix(ncol = 3, nrow = 15)
p_random <- length(fixef(commensuration_random))

coefs_random <- summary(commensuration_random)$coefficients

# Coefficient estimates and standard errors
tabCommens_random[1:5,1:2] <- coefs_random[(p_random-9):(p_random-5),1:2]
tabCommens_random[6,1:2] <- coefs_random[2,1:2]
tabCommens_random[7:11,1:2] <- coefs_random[(p_random-4):p_random,1:2]

# P-values
tabCommens_random[1:5,3] <- sapply(coefs_random[(p_random-9):(p_random-5),3],
                             twoSideP)
tabCommens_random[6,3] <- sapply(coefs_random[2,3],
                           twoSideP)
tabCommens_random[7:11,3] <- sapply(coefs_random[(p_random-4):p_random,3],
                              twoSideP)

# Random effect standard deviations
tabCommens_random[12:15,1] <- data.frame(VarCorr(commensuration_random))$sdcor
```

```{r, echo=F}
# Formatting tabCommens with row/column names
tabCommens_rownames <- c(criteria,
                   'PI Race = Black',
                   sapply(criteria, 
                          function(x) paste(x, '*', 'PI Race = Black', 
                                            sep = ' ')),
                   'Reviewer Intercept Std. Dev.', 
                   'PI Intercept Std. Dev.',
                   'SRG Intercept Std. Dev.',
                   'Residual Std. Dev.')
                   
tabCommens_colnames <- c('Estimate', 'Std. Error', 'P-value')

rownames(tabCommens_matched) <- tabCommens_rownames
colnames(tabCommens_matched) <- tabCommens_colnames

rownames(tabCommens_random) <- tabCommens_rownames
colnames(tabCommens_random) <- tabCommens_colnames
```

**Table S10, commensuration bias model: matched subset analysis, public-use data.**
```{r, echo=F}
tabCommens_matched %>%
  kable('html', row.names = T, digits = 3) %>%
  kable_styling(full_width = F) %>%
  column_spec(1, border_right = T)
```
Preliminary criterion, race, commensuration (race-criterion interaction) coefficients, and variance components estimates for preliminary overall impact scores on n = 7471 reviews of 2566 applications (matched subset). Control variables (coefficient estimates are not shown) are the matching variables.

**Table S10, commensuration bias model: random subset analysis, public-use data.**
```{r, echo=F}
tabCommens_random %>%
  kable('html', row.names = T, digits = 3) %>%
  kable_styling(full_width = F) %>%
  column_spec(1, border_right = T)
```
Preliminary criterion, race, commensuration (race-criterion interaction) coefficients, and variance components estimates for preliminary overall impact scores on n = 8595 reviews of 3045 applications (random subset). Control variables (coefficient estimates are not shown) are the matching variables.

## Commensuration Practices Analysis

The following code computes the expected differences, under a given model, in preliminary overall impact scores between black and white applicants who are assumed to have identical values on all application- and applicant-specific covariates except race. Figures S3 and S4, reproduced below, display histograms of these expected differences under the matched and random subset analyses for the public-use data. The expected differences are computed as E[Black-White], so that a positive expected difference is a worse expected score for the black applicant.
```{r}
# coefficient estimates
criteria_coef_matched <- summary(commensuration_matched)$coefficients[(p_matched-9):(p_matched-5),1]
interaction_coef_matched <- summary(commensuration_matched)$coefficients[(p_matched-4):p_matched,1]
race_coef_matched <- summary(commensuration_matched)$coefficients[2]

criteria_coef_random <- summary(commensuration_random)$coefficients[(p_random-9):(p_random-5),1]
interaction_coef_random <- summary(commensuration_random)$coefficients[(p_random-4):p_random,1]
race_coef_random <- summary(commensuration_random)$coefficients[2]

# Functions that compute the expected difference in preliminary overall impact score, black-white, for given criterion scores for a white and a black applicant
matched_score <- function(c_white, c_black) {
  dif <- c_black%*%(criteria_coef_matched + interaction_coef_matched) -
        c_white%*%criteria_coef_matched + race_coef_matched
}

random_score <- function(c_white, c_black) {
  dif <- c_black%*%(criteria_coef_random + interaction_coef_random) -
        c_white%*%criteria_coef_random + race_coef_random
}

# For each black application, compute the expected score difference between that application and an identical application from a white applicant
black_scores_matched <- apply(long[long$GROUP_ID %in% c('Matched Black', 'All Black'),criteria],
                              1, function(x) matched_score(x,x))
black_scores_random <- apply(long[long$GROUP_ID %in% c('Matched Black', 'All Black'),criteria],
                             1, function(x) random_score(x,x))
```

**Figure S3, matched subset analysis, public-use data.**
```{r,echo=F}
# Histogram of the expected score differences
hist(black_scores_matched, 
     xlab = 'Score Difference, Black - White', 
     # ylab = 'Number of Reviews of Applications',
     ylab = 'Scaled Frequency/Density',
     main = 'Distribution of Expected Impact Score Differences\nfor Black Applications, Matched Subset Commensuration Model',
     # at = seq(-1, 1, by = 0.1),
     yaxt = 'n',
     breaks = 10,
     xlim = c(-1,1))

# Axis ticks and numbering
minor.tick(nx = 5, ny = 1, tick.ratio = 0.5)
axis(2, at=seq(0,840,210), labels=seq(0,1,.25))

# Curves illustrating the variability of the random effects
x_grid <- seq(-1, 1, by = 0.01)
lines(x_grid, 
      600*dnorm(x_grid,
                sd = sqrt(summary(commensuration_matched)$varcor$REVIEWER_ID[1])),
      col = 'red')
lines(x_grid,
      1175*dnorm(x_grid,
                 sd = summary(commensuration_matched)$sigma),
      col = 'blue')

# Legend placement
legend(0.4,800, 
       legend = c('Reviewer Variability', 'Residual Variability', 'Commensuration Bias'),
       col = c('red', 'blue', 'black'),
       cex = 0.7, pch = 15)
```

Distribution of estimated expected preliminary overall impact score differences due to commensuration (histogram) an distributions of reviewer intercepts (red line) and model residuals (blue line), under the matched subset commensuration model.

**Figure S4, random subset analysis, public-use data.**
```{r,echo=F}
# Histogram of the expected score differences
hist(black_scores_random, 
     xlab = 'Score Difference, Black - White', 
     # ylab = 'Number of Reviews of Applications',
     ylab = 'Scaled Frequency/Density',
     main = 'Distribution of Expected Impact Score Differences\nfor Black Applications, Random Subset Commensuration Model',
     # at = seq(-1, 1, by = 0.1),
     yaxt = 'n',
     breaks = 10,
     xlim = c(-1,1))

# Axis ticks and numbering
minor.tick(nx = 5, ny = 1, tick.ratio = 0.5)
axis(2, at=seq(0,840,210), labels=seq(0,1,.25))

# Curves illustrating the variability of the random effects
x_grid <- seq(-1, 1, by = 0.01)
lines(x_grid, 
      580*dnorm(x_grid,
                sd = sqrt(summary(commensuration_random)$varcor$REVIEWER_ID[1])),
      col = 'red')
lines(x_grid,
      1200*dnorm(x_grid,
                 sd = summary(commensuration_random)$sigma),
      col = 'blue')

# Legend placement
legend(0.4,800, 
       legend = c('Reviewer Variability', 'Residual Variability', 'Commensuration Bias'),
       col = c('red', 'blue', 'black'),
       cex = 0.7, pch = 15)
```

Distribution of estimated expected preliminary overall impact score differences due to commensuration (histogram) an distributions of reviewer intercepts (red line) and model residuals (blue line), under the random subset commensuration model.