Exploring the UK Gender Pay Gap with R...

The gender pay gap in the UK may not of primary interest to some biochemists but the Biochemical Society is interested in gender equality and the majority of biochemistry undergraduates are female... Here is the Biochemical Society's policy statement and here is something relevant from their blog.

A report about the gender pay gap was tweeted about today by @UKParliament (it seems it was published last year). There is an Excel file that goes with the report.

Today, I've been using R to explore some of the data and I have written a script below to make these graphs - the first two graphs from the report.

The way the data is presented makes me uncomfortable with men being paid more represented as a positive percentage and women being paid more being expressed as a negative percentage. I feel sure there is a better way....

Still, the data is interesting....

Here is the script:
START
library(RCurl)
library(readxl)
library(ggplot2)
library(reshape2)
library(ggthemes)

# this is the link to the data
link <- "http://researchbriefings.files.parliament.uk/documents/SN07068/data-tables.xlsx"

# the download.file() function downloads and saves the file with the name given
download.file(url=link,destfile="file.xlsx", mode="wb")

# then we can open the file and extract the data using the read_excel() function. 
data <- read_excel("file.xlsx", skip=3, col_names=TRUE)

str(data)
# shows that Year is characters 
data[,1] <- as.numeric(data[,1])   # change to number
data <- data[1:22,]    # get rid of seven rows of NAs.
names <- colnames(data)
names[1] <- "Year"
names[2] <- "All_employees"
colnames(data) <- names  # make column names easier to use
data[,2:4] <- data[,2:4]*100   # Excel stores percents as decimals

# reshape the data from wide to long format
data.melt <- melt(data, id.vars = "Year")
colnames(data.melt) <- c("Year", "empType", "gendGap")


# draw the graph
p1 <- ggplot(data.melt, aes(x=Year, 
                           y= gendGap, 
                           colour = empType)) + 
  geom_point() +   # draw the points
  geom_line(size=1) +  # draw the lines
  labs(color = "Employment Type") + # customizes the legend title
  ylab("Gender Gap (%)") + # y-label
  ggtitle("Gender Pay Gap, UK, 1997-2015") +   # graph title
  ylim(-10,30) + 
  xlim(1995, 2015) +
  geom_hline(yintercept = 0) +  # nice line at zero
  theme_bw()

p1 <- p1 + theme(legend.text=element_text(size = 12), # increase size of text
               legend.title=element_text(size = 12)) # and title

p1 <- p1 + theme(axis.title.y = element_text(size = 14 )) + 
  theme(axis.text = element_text(size = 12))

p1 # show the graph

# maybe you prefer a different theme.
p1 + theme_hc()

# maybe without a legend but with labels on the lines:
p1 <- p1 + theme(legend.position="none") + 
  geom_text(data = data.melt[which(data.melt$Year == "2013"),],
               aes(label = empType),
               vjust = -2)
p1

# draw the second graph with the age data... 
data2 <- read_excel("file.xlsx", sheet=2, skip=3, col_names=TRUE)
View(data2)
str(data2)
data2 <- data2[1:8,]

# multiply numbers by 100 to give percentages
data2[,2:4] <- data2[,2:4]*100

names <- colnames(data2)
names[1] <- "Ages"
names[2] <- "All_employees"
colnames(data2) <- names
data2.type <- data2[3:8,]

data2.melt <- melt(data2.type, id.vars = "Ages")
colnames(data2.melt) <- c("Ages", "empType", "gendGap")

g <- ggplot(data = data2.melt[7:18,], aes (x = Ages, y = gendGap, fill = empType))  
g <- g + geom_bar(stat="identity", position="dodge", width = 0.75) +
  ylim(-11,20)  +
  ylab("Gender Gap (%)") + # y-label
  xlab("Age") + # x-label
  ggtitle("Gender Pay Gap by Age, April 2015") +
  labs(fill = "Employment Type") + # customizes the legend title
  theme_hc() +
  theme(legend.position=c(0,1), # move to the top left
        legend.justification=c(0,1.5)) # move it in a bit
g  # show the graph...

Overview of gender breakdown across types of staff

I have spent a large part of the last month grappling with Athena SWAN data for the School of Medicine at Cardiff University. To apply for this equality award, the Athena SWAN Self Assessment Team needed to look at the gender breakdown of students and staff in the School.

Here is a bar chart I generated with the data to give an overview of the School:

Here's is the script I used to generate this data and some other graphs that explore using ggplot to make bar charts.

START
# overview of the School in terms of gender for our Athena SWAN application
library(readxl)
library(ggplot2)
library(reshape2)

# the data is here
link <- "https://raw.githubusercontent.com/brennanpincardiff/AthenaSWANBenchmarkData/master/asStaffOverview_comparedtype_20160505.xlsx"

# the download.file() function downloads and saves the file with the name given
download.file(url=link, destfile="file.xlsx", mode="wb")

# read in the data...
data <- read_excel("file.xlsx")

# calculate Female and Male as a percentage
data$tot <- data$f + data$m
data$Female <- (data$f / data$tot) *100 
data$Male <- (data$m / data$tot) *100

# first a simple plot of the Cardiff School of Medicine Data
# plot the number of females of each type of member of the School
# pull out the Cardiff data
data.c <- data[data$uni == "Cardiff",2:5]

# make the object with the Cardiff subset  
g <- ggplot(data.c,
       aes(x = type,
           y = f)) +
     geom_bar(stat="identity")+
     theme(axis.text.x = element_text(size = 12))

# show the graph
g

# The order of the bars is not what I want.
# take the order from the data file
types <- data.c$type

# apply this order to the graph
# http://stackoverflow.com/questions/5208679/order-bars-in-ggplot2-bar-graph
g + scale_x_discrete(limits = types)

# I want to put female and male on the same graph so...
# pull out and melt the Cardiff data so that have f and m for each type in long format
data.melt <- melt(data[data$uni == "Cardiff",2:5], 
                  id.var = c("type", "year"))

# put the data in an object
g <- ggplot(data.melt, 
       aes(x = type, 
           y = value,
           fill = variable)) + 
      scale_x_discrete(limits = types)

# show with geom_bar
# this is a stacked bar chart showing both female and male numbers
g + geom_bar(stat="identity")

# show female and male in separate bars using position arguement
g + geom_bar(stat="identity", position="dodge")

# have as stacked bar chart but normalised to one. 
g +  geom_bar(stat="identity", position="fill")

# I really want a percentage of each stacked 
# to allow different parts of the School to be compared
# remove the raw numbers
data.s <- data[, -c(3:6)]
data.melt.s <- melt(data.s[data.s$uni == "Cardiff",2:4], 
                  id.var = c("type"))  # id = or id.var are same
colnames(data.melt.s) <- c("type", "Gender", "Percent")

# make a new graph
p <- ggplot(data.melt.s, 
           aes(x = type, 
               y = Percent,
               fill = Gender)) +
      geom_bar(stat="identity") + 
      scale_x_discrete(limits = types) +
      ylab("Percentage of staff") +
      xlab("") +
      theme_bw() +
      theme(axis.text.x = element_text(size = 14))
p # show the graph

# adjust the colours using scale_fill_brewer()
p + scale_fill_brewer(palette = 16, direction=-1) +
  # and add a title
    ggtitle("Gender breakdown of different members of the School of Medicine") +
    theme(axis.title.y = element_text(size = 14 ))

# use data from Imperial & Leeds as comparison
# use the facet_wrap function to make all looks the same. 
# remove Prof & Support as we only have for Cardiff
data.s <- data.s[-5,]
types <- types[-5]

# melt the data into the format for the bar chart. 
data.melt.s <- melt(data.s, 
                    id.var = c("uni", "type"))
colnames(data.melt.s) <- c("uni", "type", "Gender", "Percent")

# this is a final plot with three Medical Schools (it's a bit big). 
p <- ggplot(data.melt.s, 
            aes(x = type, 
                y = Percent,
                group = uni,
                fill = Gender)) +
  geom_bar(stat="identity") +
  facet_wrap(~uni) +    # this function gives three plots
  scale_x_discrete(limits = types) +
  ylab("Percentage of staff") +
  xlab("") +
  theme_bw() +
  scale_fill_brewer(palette = 16, direction=-1) +
  ggtitle("Gender breakdown of three Schools of Medicine") +
  theme(strip.text.x = element_text(size = 14, colour = "black")) +
  theme(axis.text.x = element_text(size = 11))

p # show the plot

END of SCRIPT

I looked up these pages for help:

• http://stackoverflow.com/questions/3290330/facet-label-font-size
• http://www.sthda.com/english/wiki/ggplot2-axis-ticks-a-guide-to-customize-tick-marks-and-labels

As a biochemist, where would you like to work from an Athena SWAN point of view...

As part of my role in the School of Medicine at Cardiff University, I am helping to prepare an application for an Athena SWAN Award. As part of this, I have prepared a list of Athena SWAN Awards across the UK over the last few years and looking at other applications to learn more about the process. The data is on github. I have been using this data to find Medical Schools that have Athena SWAN awards to use as comparison data for our application and to use as a style guide.

Athena SWAN is a quality mark that focuses on equality. As such, departments and Institutions with Athena SWAN awards should be good places to work. Furthermore the better the award, the better the institution as a work place. I decided to analyse the Department names from the point of view of a biochemist.

Here is a graph, I produced (similar to one I produced before and published here):

To do the analysis, I made a user-defined function using function(). The function is called countAwardType(). To use the function, you put in an argument - the name of a department or institution in the brackets: countAwardType("Biochemistry").

There are a two Gold Awards for "Bio" Departments. These are:

• Queen’s University Belfast, School of Biological Sciences
• University of York, Department of Biology

There are three for Chemistry:

• University of York, Department of Chemistry
• Imperial College London, Department of Chemistry
• School of Chemistry, University of Edinburgh

Oxford University has the most awards - see below for Russell group graph.

SCRIPT START
# As a biochemist, where would you like to work for from an Athena SWAN point of view?

library(RCurl)
library(readxl)

# this is the link to the data
link <- "https://raw.githubusercontent.com/brennanpincardiff/AthenaSWANBenchmarkData/master/ASAwards/AthenaSWANAwardList2012_201520160316.xlsx"

# the download.file() function downloads and saves the file with the name given
download.file(url=link, destfile="file.xlsx", mode="wb")

# read in the file with the read_excel() function
awards<- read_excel("file.xlsx", col_names = TRUE)

award.type <- c("Gold", "Silver", "Bronze")

# I have created a function to count the different types of Award - Gold, Silver & Bronze
# and then return a data.frame with specific names
countAwardType <- function(x){
  
  # to test this function  make x <- c("Medicine") and go through it line by line. 
  
  # list of the different types of Athena SWAN awards
  award.type <- c("Gold", "Silver", "Bronze")
  
  # find the pattern supplied in the function as x e.g. "Medicine"
  # this creates a data.frame. 
  x.Awards <- awards[grep(x, awards$org),]
  
  # use split() function which separates the data.frame into a list by award level
  # there is a data.frame within each list.
  # this can be called by the dollar sign - $. 
  x.split <- split(x.Awards, x.Awards$Level)
  
  # count the number of rows in each list 
  gold <- nrow(x.split$Gold)
  silver <- nrow(x.split$Silver)
  bronze <- nrow(x.split$Bronze)
  
  # if there are no entries then a NULL value will be returned 
  # zero is better therefore recode the NULL values
  if(is.null(gold) == TRUE) {gold <- 0}
  if(is.null(silver) == TRUE) {silver <- 0}
  if(is.null(bronze) == TRUE) {bronze <- 0}
  
  # combine these variables into a vector
  aW <- c(gold, silver, bronze)
  
  # combine the vectors into a data.frame
  awardsDepts.x <- data.frame(x, aW, award.type)
  
  # rename the columns
  colnames(awardsDepts.x) <- c("Dept.Type", "Award.Count", "Award.Type")
  
  # return the answer. 
  return(awardsDepts.x)
}

# with this function, I can answer the following questions....

# How many Departments of "Biochemistry" have Athena SWAN awards?
# I can answer this using the grep() function 
# this searches for the word "Biochemistry" in the org column from the awards data.frame
# it creates a data.frame called biochemAwards
biochemAwards <- awards[grep("Biochemistry", awards$org),]
nrow(biochemAwards)  # answer 5

# or
# I can use my function countAwardType()
countAwardType("Biochemistry")

# returns this output:
# Dept.Type Award.Count Award.Type
# 1 Biochemistry           0       Gold
# 2 Biochemistry           1     Silver
# 3 Biochemistry           4     Bronze


# As a biochemist you might be willing to work in any "Bio" department
# these include "Biochemistry" of course. 
# or a "Chem" Department 
# or a "Medic" Department like I do...
# or perhaps a Pharmacology department "Pharm"

# lets make a list. 

biochemistry.options <- c("Biochemistry", "Bio", "Chem", "Medic", "Pharm", "Immun")

# loop through the list - make a data.frame for plotting
awardsDepts <- NULL
for(i in 1:length(biochemistry.options)){
  awardsDepts <- rbind(awardsDepts, countAwardType(biochemistry.options[i])) 
}

# re-order the data.frame
awardsDepts$Award.Type <- factor(awardsDepts$Award.Type, levels = c("Bronze", "Silver", "Gold"))

# let's draw a graph
# make the first version of the plot
p <- ggplot(awardsDepts, 
            aes(x=Dept.Type, 
                y=Award.Count, 
                fill=Award.Type)) + 
            geom_bar(stat="identity") +  # makes the barplot
     scale_fill_manual(values=c("#956C3E", 
                                "#757576", 
                                "#A28D30")) + 
     # colours as per AS website - see below
     xlab("") +  # no need for the "Depts" x-axis
     theme_few() # nice clean theme  

# I want to add a y-axis label
# and increase the size of the text
p <- p + ylab("Number of Awards") +
  theme(axis.title.y = element_text(size = 14 )) + 
  theme(axis.text = element_text(size = 14))

# make the legend a bit bigger and move it
p <- p + theme(legend.position=c(0.9,0.85), # top right
               legend.text=element_text(size = 12), # inc title size
               legend.title=element_text(size = 12)) # labels

# have a look at the bar chart
p + ggtitle("Departments Types with Athena SWAN Awards")

# Don't forget some of the quality Research Institutes: 
# Babraham and Institue of Cancer Research both have Silver Awards. 


# by way of contrast let's look at some other department types:
departments <- c("Bio", "Chem", "Medic", "Engin", "Comp", "Math", "Psych", "Physic")
awardsDepts <- NULL
for(i in 1:length(departments)){
  awardsDepts <- rbind(awardsDepts, countAwardType(departments[i])) 
}
awardsDepts$Award.Type <- factor(awardsDepts$Award.Type, levels = c("Bronze", "Silver", "Gold"))

# push in the new data
d <- p %+% awardsDepts
d + ggtitle("Departments Types with Athena SWAN Awards")

# we can also compare Institutions
russell <- c("University of Birmingham", "University of Bristol", "University of Cambridge", 
             "Cardiff University", "Durham University", "University of Edinburgh",
             "University of Exeter", "University of Glasgow", "Imperial College London",
             "King’s College London", "University of Leeds", "University of Liverpool",
             "London School of Economics and Political Science", "University of Manchester",
             "Newcastle University", "University of Nottingham", "University of Oxford",
             "Queen Mary, University of London", "Queen’s University Belfast", 
             "University of Sheffield", "University of Southampton", "University College London",
             "University of Warwick", "University of York")

awardsDepts <- NULL
for(i in 1:length(russell)){
  awardsDepts <- rbind(awardsDepts, countAwardType(russell[i])) 
}
awardsDepts$Award.Type <- factor(awardsDepts$Award.Type, levels = c("Bronze", "Silver", "Gold"))

# push in the new data
inst <- p %+% awardsDepts
inst + ylim(0,50) +  # extend the y axis
    theme(legend.position=c(0.1,0.85)) +   # move the legend
    theme(axis.text.x=element_text(angle=90, hjust=1, vjust=0.5 )) +                    # vertical names
    ggtitle("Russell group Universities with Athena SWAN Awards")

# finally Who's got Gold?
awards[grep("Gold", awards$Level),]
Type New_Renew Level Year Month                                                       org
# 70  Dept   Renewal  Gold 2015 Apr University of York, Department of Chemistry
# 360 Dept       N/A  Gold 2013 Nov Queen’s University Belfast, School of Psychology
# 361 Dept       N/A  Gold 2013 Nov University of Cambridge, Department of Physics
# 362 Dept       N/A  Gold 2013 Nov University of York, Department of Biology
# 452 Dept       N/A  Gold 2013 Apr Imperial College London, Department of Chemistry
# 539 Dept       New  Gold 2012   Nov Queen’s University Belfast, School of Biological Sciences
# 540 Dept       New  Gold 2012   Nov School of Chemistry, University of Edinburgh


# colours from AthenaSWAN website
# https://athena-swan.medschl.cam.ac.uk/wp-content/uploads/2014/03/Athena-SWAN_style-guide_October-2012.pdf
# Bronze: R149 G108 B62
# Silver: R117 G117 B118
# Gold: R162 G141 B48
# Convert RGB to hex http://www.rgbtohex.net/