mentorship.Rmd

---
title: "Mentorship"
author: "Daniel E. Weeks"
date: "`r format(Sys.time(), '%B %d, %Y, %R')`"
header-includes:
- \usepackage{pdflscape}
- \newcommand{\blandscape}{\begin{landscape}}
- \newcommand{\elandscape}{\end{landscape}}
output:
  pdf_document:
    toc: true
    number_sections: true
    toc_depth: 3
  html_document:
    df_print: paged
    toc: true
    number_sections: true
    toc_depth: '3'
    code_folding: show
#   md_extensions: +raw_attribute
  github_document:
    toc: true
    number_sections: true
    toc_depth: 3
---

```{r,echo=FALSE,message=FALSE,warning=FALSE}
require(knitr)
# Set so that long lines in R will be wrapped:
opts_chunk$set(tidy.opts=list(width.cutoff=80),tidy=TRUE)
```

```{r,echo=FALSE}
# In the event of a crash, save the working space on exit:
save_all <- function() {
  save.image("recover.RData")
}
options(error = save_all)
```


# Load Libraries

```{r load_libraries,message=FALSE}
library(tidyverse)
# library(tidylog)
library(ggExtra)
library(arsenal)
```

# Input directory and files 

```{r}
# Print the working directory
getwd()
```

# The AlShebli et al (2020) mentorship paper

Here I examine the data provided by the authors of this paper:

AlShebli B, Makovi K, Rahwan T. The association between early career informal mentorship in academic collaborations and junior author performance. Nat Commun. 2020 Nov 17;11(1):5855. doi: 10.1038/s41467-020-19723-8. PMID: 33203848.
https://pubmed.ncbi.nlm.nih.gov/33203848/

# The data

## The input data

The input data are from the 'bedoor/Mentorship' GitHub repository at

https://github.com/bedoor/Mentorship

Their repository is described as:

'This repository includes all data used in "The Association between Early Career Informal Mentorship in Academic Collaborations and Junior Author Performance".'

```{r}
LinesRead <- 100000
inFile <- "Mentorship/Repository_Data/Data_7yearcutoff.csv"
a <- read_csv(file = inFile, n_max = LinesRead)
names(a)
```

```{r}
total_records <- as.integer(system2("wc",
                                    args = c("-l",
                                             inFile,
                                             " | awk '{print $1}'"),
                                    stdout = TRUE))
total_records
```

## Disclaimer: For speed, only read in the first `r LinesRead` lines

To speed up the analyses, we read only the first `r LinesRead` lines of the input file, instead of reading all `r total_records` lines.

## Problem: No gender information for mentors

There does not appear to be any information in the provided file about the mentor's gender, so the provided data do not appear to be 'all the data'.

There is no data dictionary, so we are left to infer what is in each column of the data file from the column names and content.


# ProtegeFirstPubYear

# Distribution of ProtegeFirstPubYear

"Year of the protege’s first publication: The year in which the protege published their first mentored paper."

```{r, dpi=600, dev='png'}
summary(a$ProtegeFirstPubYear)
ggplot(data = a, mapping = aes(x = ProtegeFirstPubYear)) + geom_histogram() +
  ggtitle("Distribution of ProtegeFirstPubYear")

table(a$Pr0tegeGender)

ggplot(data = a, mapping = aes(x = ProtegeFirstPubYear)) + geom_histogram() +
  ggtitle("Distribution of ProtegeFirstPubYear") + facet_grid(~ Pr0tegeGender)

a$MaleProtegeFirstPubYear <- a$ProtegeFirstPubYear
a$MaleProtegeFirstPubYear[a$Pr0tegeGender=="female"] <- NA

a$FemaleProtegeFirstPubYear <- a$ProtegeFirstPubYear
a$FemaleProtegeFirstPubYear[a$Pr0tegeGender=="male"] <- NA

summary(a$MaleProtegeFirstPubYear)
summary(a$FemaleProtegeFirstPubYear)

p <- ggplot(data=a, aes(x=x) ) +
  geom_histogram(aes(x = MaleProtegeFirstPubYear, y = ..density..), fill="#69b3a2" ) +
  annotate("text", x=1950, y=0.03, label="MaleProtegeFirstPubYear", color="#69b3a2") +
  geom_histogram(aes(x = FemaleProtegeFirstPubYear, y = -..density..), fill= "#404080") +
  annotate("text", x=1950, y=-0.03, label="FemaleProtegeFirstPubYear", color="#404080") +
  xlab("ProtegeFirstPubYear") +
  ggtitle("Distributions of ProtegeFirstPubYear by gender")
p

p <- ggplot(data=a, aes(x=x) ) +
  geom_density(aes(x = MaleProtegeFirstPubYear, y = ..density..), fill="#69b3a2" ) +
  annotate("text", x=1950, y=0.03, label="MaleProtegeFirstPubYear", color="#69b3a2") +
  geom_density(aes(x = FemaleProtegeFirstPubYear, y = -..density..), fill= "#404080") +
  annotate("text", x=1950, y=-0.03, label="FemaleProtegeFirstPubYear", color="#404080") +
  xlab("ProtegeFirstPubYear") +
  ggtitle("Distributions of ProtegeFirstPubYear by gender")
p
```

## Problem: ProtegeFirstPubYear time range is very large

in the first `r LinesRead` lines of this data set, the `ProtegeFirstPubYear` ranges from `r min(a$ProtegeFirstPubYear)` to `r max(a$ProtegeFirstPubYear)`.  It seems that it would be very difficult to properly compare the impact factors of papers published in `r `min(a$ProtegeFirstPubYear)`, well before the advent of team science and the current deluge of scientific publications, to those published as recently `r max(a$ProtegeFirstPubYear)`.


# NumYearsPostMentorship

"The number of years post mentorship: Since our dataset is up to Dec 31st 2019,
we are only able to calculate c5—the number of citations accumulated five years post publication—for papers published before 2015. Thus, given a protege whose first paper was published in year x, the number of years post mentorship is `2015 - (x + 6)`, bearing in mind that the mentorship period is 7 years and we do not include proteges with a gap of 5 years of more in their career history."

## Problem: Error in the equation for NumYearsPostMentorship

Note that the stated equation of `2015 - (x + 6)` for computing the number of years post mentorship is not what was used, as instead `2015 - (x + 2)`, where `x` is the year the protege's paper was first published, as we see here in the first `r LinesRead` lines of the input file:

```{r}
all.equal(a$NumYearsPostMentorship, 2015 -(a$ProtegeFirstPubYear + 2))
```

## Distribution of NumYearsPostMentorship

```{r, dpi=600, dev='png'}
summary(a$NumYearsPostMentorship)
ggplot(data = a, mapping = aes(x = NumYearsPostMentorship)) + geom_histogram() +
  ggtitle("Distribution of NumYearsPostMentorship")
```

```{r, dpi=600, dev='png'}
a$ProtegeGender <- as.factor(a$Pr0tegeGender)
p <- ggplot(data = a, mapping = aes(x = ProtegeFirstPubYear, y = NumYearsPostMentorship)) + geom_point() + ggtitle("NumYearsPostMentorship vs. ProtegeFirstPubYear")
ggMarginal(p, type="histogram") 

```

## Problem: some values of `NumYearsPostMentorship` are unrealistically large

In the first `r LinesRead` lines of the input file of this data set, there are values of `NumYearsPostMentorship` as large as `r max(a$NumYearsPostMentorship)`.  These values are unrealistically large. 

# AvgMentorsAcAges

## Distribution of AvgMentorsAcAges

"Average academic age of mentors: This is computed for any given protege by first
computing the academic age of each mentor in the year of their first publication with the protege, and then averaging these numbers over all the mentors."

"Given a scientist whose first paper was published in year x, the academic age of this scientist in year y is y-x."

```{r, dpi=600, dev='png'}
summary(a$AvgMentorsAcAges)
ggplot(data = a, mapping = aes(x = AvgMentorsAcAges)) + geom_histogram() +
  ggtitle("Distribution of AvgMentorsAcAges")

ggplot(data = a, mapping = aes(x = ProtegeFirstPubYear, y = AvgMentorsAcAges, col=ProtegeGender)) + geom_point() +
  ggtitle("ProtegeFirstPubYear vs. AvgMentorsAcAges")

ggplot(data = a, mapping = aes(x = ProtegeFirstPubYear, y = AvgMentorsAcAges)) + geom_point() + facet_grid(~ ProtegeGender) +
  ggtitle("ProtegeFirstPubYear vs. AvgMentorsAcAges by ProtegeGender")
```

## Problem: some `AvgMentorsAcAges` are unrealistically large

There are `AvgMentorsAcAges` as large as `r max(a$AvgMentorsAcAges)` in the first `r LinesRead` lines of this data set.  This is unrealistically large. 

Among proteges who published their first paper after 2000, there are `AvgMentorsAcAges` as large as `r max(a[a$ProtegeFirstPubYear > 2000, ]$AvgMentorsAcAges)` in this data set.  Should we be evaluating the mentorship effects of mentors who were born than 200 years ago?


# Avg_c5

## Distribution of Avg_c5

```{r, dpi=600, dev='png'}
summary(a$Avg_c5)
ggplot(data = a, mapping = aes(x = Avg_c5)) + geom_histogram() +
  ggtitle("Distribution of Avg_c5")

ggplot(data = a, mapping = aes(x = ProtegeFirstPubYear, y = Avg_c5, col=ProtegeGender)) + geom_point() +
  ggtitle("ProtegeFirstPubYear vs. Avg_c5")

ggplot(data = a, mapping = aes(x = ProtegeFirstPubYear, y = Avg_c5)) + geom_point() + facet_grid( ~ ProtegeGender) +
  ggtitle("ProtegeFirstPubYear vs. Avg_c5")
```


## Joint distribution of Avg_c5 and AvgMentorsAcAges

```{r, dpi=600, dev='png'}
ggplot(data = a, mapping = aes(y = Avg_c5, x = AvgMentorsAcAges, col=ProtegeGender)) + geom_point() + 
  ggtitle("Avg_c5 vs. AvgMentorsAcAges")

ggplot(data = a, mapping = aes(y = Avg_c5, x = AvgMentorsAcAges, col=ProtegeFirstPubYear)) + 
  geom_point() + 
  scale_color_gradient(low="red", high="blue") +
  facet_grid( ~ ProtegeGender) +
  ggtitle("Avg_c5 vs. AvgMentorsAcAges")
```

# Avg_c10

## Distribution of Avg_c10

```{r, dpi=600, dev='png'}
summary(a$Avg_c10)
ggplot(data = a, mapping = aes(x = Avg_c10)) + geom_histogram() +
  ggtitle("Distribution of Avg_c10")

ggplot(data = a, mapping = aes(x = ProtegeFirstPubYear, y = Avg_c10, col=ProtegeGender)) + geom_point() +
  ggtitle("ProtegeFirstPubYear vs. Avg_c10")

ggplot(data = a, mapping = aes(x = ProtegeFirstPubYear, y = Avg_c10)) + geom_point() + facet_grid( ~ ProtegeGender) +
  ggtitle("ProtegeFirstPubYear vs. Avg_c10")
```



## Problem with Avg_c10: Cannot compute a value Avg_c10 for recent publications

As the dataset is up to Dec 31st 2019, how can one compute a valid Avg_c10 for someone whose first publication year was 2012 or 2013?

```{r}
a %>% filter(ProtegeFirstPubYear == max(ProtegeFirstPubYear)-1) %>% select(ProtegeFirstPubYear, Avg_c5, Avg_c10)  %>% head(20)
```

# Avg_c5 vs Avg_c10

## Tables of Avg_c5 and Avg_c10 

```{r}

my_controls <- tableby.control(
  test = T,
  total = T,
  numeric.stats = c("meansd", "medianq1q3", "range", "Nmiss2"),
  cat.stats = c("countpct", "Nmiss2"),
  stats.labels = list(
    meansd = "Mean (SD)",
    medianq1q3 = "Median (Q1, Q3)",
    range = "Min - Max",
    Nmiss2 = "Missing"
  )
)

my_controls2 <- tableby.control(
  numeric.stats = c( "medianq1q3"),
  stats.labels = list(
    medianq1q3 = "Median (Q1, Q3)"
  )
)

t.test(Avg_c5 ~ ProtegeGender, data=a)

t.test(Avg_c10 ~ ProtegeGender, data=a)

t1 <- tableby(ProtegeGender ~ ., data = a %>% select(ProtegeGender, Avg_c5, Avg_c10))
```

```{r, results='asis'}
summary(t1, title = "Table of Avg_c5 and Avg_c10 by ProtegeGender")
```

```{r}
t1 <- tableby(ProtegeGender ~ ., data = a %>% select(ProtegeGender, Avg_c5, Avg_c10), control = my_controls2)
```

```{r, results='asis'}
summary(t1, title = "Table of Avg_c5 and Avg_c10 by ProtegeGender")
```

```{r}
t1 <- tableby(ProtegeGender ~ ., data = a %>% select(ProtegeGender, Avg_c5, Avg_c10), control = my_controls)
```

```{r, results='asis'}
summary(t1, title = "Table of Avg_c5 and Avg_c10 by ProtegeGender")
```

## Joint distribution of Avg_c5 and Avg_c10

```{r, dpi=600, dev='png'}
ggplot(data = a, mapping = aes(x = Avg_c5, y = Avg_c10, col=ProtegeGender)) + geom_point() + 
  geom_abline(intercept = 0, slope = 1, col="red") +
  ggtitle("Avg_c5 vs. Avg_c10")

ggplot(data = a, mapping = aes(x = Avg_c5, y = Avg_c10, col=ProtegeFirstPubYear)) + 
  geom_point() + 
  scale_color_gradient(low="red", high="blue") +
  geom_abline(intercept = 0, slope = 1, col="red") +
  facet_grid( ~ ProtegeGender) +
  ggtitle("Avg_c5 vs. Avg_c10")
```

# numMentors

## Distribution of numMentors

```{r, dpi=600, dev='png'}
summary(a$numMentors)
ggplot(data = a, mapping = aes(x = numMentors)) + geom_histogram() +
  ggtitle("Distribution of numMentors")

ggplot(data = a, mapping = aes(x = ProtegeFirstPubYear, y = numMentors, col=ProtegeGender)) + geom_point() +
  ggtitle("ProtegeFirstPubYear vs. numMentors")

ggplot(data = a, mapping = aes(x = ProtegeFirstPubYear, y = numMentors)) + geom_point() + facet_grid( ~ ProtegeGender) +
  ggtitle("ProtegeFirstPubYear vs. numMentors")
```

## Problem: some `numMentors` values are too large

How could someone possibly really have more than a handful of mentors?

In the first `r LinesRead` lines of this data set, there are individuals with `r max(a$numMentors)` mentors!

The Supplement states that:

"Whenever a junior scientist publishes a paper with a senior scientist, we consider the former to be a protege, and the latter to be a mentor, as long as they authored at least one paper with 20 or less co-authors and share the same discipline and US-based affiliation."

Oh, but we are really measuring co-authorship at the same US-based affiliation on joint papers with 20 or less co-authors. 

# Summary

Here I examine the data provided by the authors of this paper:

AlShebli B, Makovi K, Rahwan T. The association between early career informal mentorship in academic collaborations and junior author performance. Nat Commun. 2020 Nov 17;11(1):5855. doi: 10.1038/s41467-020-19723-8. PMID: 33203848.
https://pubmed.ncbi.nlm.nih.gov/33203848/

For speed, I examined only the first `r LinesRead` lines of the data file `r inFile` from the 'bedoor/Mentorship' GitHub repository at

https://github.com/bedoor/Mentorship

Their repository is described by the authors as:

'This repository includes all data used in "The Association between Early Career Informal Mentorship in Academic Collaborations and Junior Author Performance".'

## Problem: No gender information for mentors

There does not appear to be any information in the provided file about the mentor's gender, so the provided data do not appear to be 'all the data'.

## Problem: ProtegeFirstPubYear time range is very large

In the first `r LinesRead` lines of this data set, the `ProtegeFirstPubYear` ranges from `r min(a$ProtegeFirstPubYear)` to `r max(a$ProtegeFirstPubYear)`.  It seems that it would be very difficult to properly compare the impact factors of papers published in `r min(a$ProtegeFirstPubYear)`, well before the advent of team science and the current deluge of scientific publications, to those published as recently `r max(a$ProtegeFirstPubYear)`.

```{r, dpi=600, dev='png', echo=FALSE}
ggplot(data = a, mapping = aes(x = ProtegeFirstPubYear, y = Avg_c5)) + geom_point() + facet_grid( ~ ProtegeGender) +
  ggtitle("ProtegeFirstPubYear vs. Avg_c5")
```

## Problem: Error in the equation for NumYearsPostMentorship

Note that the stated equation of `2015 - (x + 6)` for computing the number of years post mentorship is not what was used, as instead `2015 - (x + 2)`, where `x` is the year the protege's paper was first published.

## Problem: some values of `NumYearsPostMentorship` are unrealistically large

In the first `r LinesRead` lines of the input file of this data set, there are values of `NumYearsPostMentorship` as large as `r max(a$NumYearsPostMentorship)`.  These values are unrealistically large. 

```{r, dpi=600, dev='png', echo=FALSE}
a$ProtegeGender <- as.factor(a$Pr0tegeGender)
p <- ggplot(data = a, mapping = aes(x = ProtegeFirstPubYear, y = NumYearsPostMentorship)) + geom_point() + ggtitle("NumYearsPostMentorship vs. ProtegeFirstPubYear")
ggMarginal(p, type="histogram") 

```

## Problem: some `AvgMentorsAcAges` are unrealistically large

There are `AvgMentorsAcAges` as large as `r max(a$AvgMentorsAcAges)` in the first `r LinesRead` lines of this data set.  This is unrealistically large. 

Among proteges who published their first paper after 2000, there are `AvgMentorsAcAges` as large as `r max(a[a$ProtegeFirstPubYear > 2000, ]$AvgMentorsAcAges)` in this data set.  Should we be evaluating the mentorship effects of mentors who were born than 200 years ago?

```{r, dpi=600, dev='png', echo=FALSE}
ggplot(data = a, mapping = aes(x = ProtegeFirstPubYear, y = AvgMentorsAcAges)) + geom_point() + facet_grid(~ ProtegeGender) +
  ggtitle("ProtegeFirstPubYear vs. AvgMentorsAcAges by ProtegeGender")
```

## Problem with Avg_c10: Cannot compute a value `Avg_c10` for recent publications

As the dataset is up to Dec 31st 2019, how can one compute a valid `Avg_c10` for someone whose first publication year was 2012 or 2013?

```{r, dpi=600, dev='png', echo=FALSE}
ggplot(data = a, mapping = aes(x = ProtegeFirstPubYear, y = Avg_c10)) + geom_point() + facet_grid( ~ ProtegeGender) +
  ggtitle("ProtegeFirstPubYear vs. Avg_c10")
```

## Problem: some `numMentors` values are too large

How could someone possibly really have more than a handful of mentors?

In the first `r LinesRead` lines of this data set, there are individuals with `r max(a$numMentors)` mentors!

The Supplement states that:

"Whenever a junior scientist publishes a paper with a senior scientist, we consider the former to be a protege, and the latter to be a mentor, as long as they authored at least one paper with 20 or less co-authors and share the same discipline and US-based affiliation."

Oh, but we are really measuring **co-authorship** at the same US-based affiliation on joint papers with 20 or less co-authors. 

```{r, dpi=600, dev='png', echo=FALSE}
ggplot(data = a, mapping = aes(x = ProtegeFirstPubYear, y = numMentors)) + geom_point() + facet_grid( ~ ProtegeGender) +
  ggtitle("ProtegeFirstPubYear vs. numMentors")
```


# Generating GitHub Markdown

```{r}
library(rmarkdown)
library(here)
finalize <- function() {
  rmd <- here("mentorship.Rmd")
  rmarkdown::render(rmd, "github_document")
}
```

# Session Information

```{r}
sessionInfo()
```