2_audience_genderPosition.Rmd

---
title: "Gender bias in audience of seminars and career position"
author: "Melina Leite, Júlia Barreto & Isabella Romitelli"
date: "`r format(Sys.time(), '%d de %B de %Y')`"
output:
  rmdformats::readthedown:
    highlight: kate
    self_contained: true
    thumbnails: false
    lightbox: true
    gallery: false
  pdf_document:
    highlight: tango
    toc: yes
---
```{r setup, echo=FALSE, warning=FALSE, message=FALSE}
library(knitr)
library(tidyverse); library(cowplot); library(patchwork)
theme_set(theme_cowplot())
library(janitor); library(lubridate)
library(DHARMa); library(ggeffects)
library(bbmle); library(MASS);
library(performance)
library(FactoMineR); library(factoextra)
library(glmmTMB)


opts_chunk$set(fig.align = 'center', warning = FALSE, message = FALSE, error = FALSE, echo=T, cache=F)
options(formatR.arrow = TRUE, width = 90, help_type = "html")
```

# Data

Data description and summary statitiscs in script `0_data_summary`.

```{r}
load("data_clean/data_audience.Rdata")
data <- data_a
```

# Modeling

The audience (number of attendants in the seminar) is modeled by the gender (male, female), the academic level (student, postdoc, professor) of the speaker, and affirmative actions (before, after).

Negative binomial mixed-effects model with year as random intercept.

```{r}
mg0 <- glmmTMB(audience_n~ 1 + (1|year), data=data, family= nbinom2)
mg1 <- glmmTMB(audience_n~ gender + (1|year), data=data, family= nbinom2)
mg2 <- glmmTMB(audience_n~ position_cat + (1|year), data=data, family= nbinom2)
mg3 <- glmmTMB(audience_n~ affirm_action + (1|year), data=data, family= nbinom2)

mg4 <- glmmTMB(audience_n~ gender + position_cat + (1|year), 
               data=data, family= nbinom2)
mg5 <- glmmTMB(audience_n~ gender + affirm_action + (1|year), 
               data=data, family= nbinom2)
mg6 <- glmmTMB(audience_n~ affirm_action + position_cat+ (1|year), 
               data=data, family= nbinom2)

mg7 <- glmmTMB(audience_n~ gender * position_cat+ (1|year), 
               data=data, family= nbinom2)
mg8 <- glmmTMB(audience_n~ gender * affirm_action+ (1|year), 
               data=data, family= nbinom2)
mg9 <- glmmTMB(audience_n~ affirm_action * position_cat+ (1|year), 
               data=data, family= nbinom2)

mg10 <- glmmTMB(audience_n~ gender + position_cat + affirm_action + (1|year), 
               data=data, family= nbinom2)
mg11 <- glmmTMB(audience_n~ gender * position_cat + affirm_action + (1|year), 
               data=data, family= nbinom2)
mg12 <- glmmTMB(audience_n~ gender + position_cat * affirm_action + (1|year), 
               data=data, family= nbinom2)

mg13 <- glmmTMB(audience_n~ gender * position_cat * affirm_action + (1|year), 
               data=data, family= nbinom2)

AICtab(mg0,mg1,mg2, mg3, mg4,mg5,mg6,mg7,mg8,mg9,mg10,mg11,mg12,mg13, base=T, 
       weights=T) %>% 
  kable(digits=2)
```


## Residual diagnostic

Using the `DHARMa` package.

The two most plausible models presented a satisfactory residual diagnostic.

```{r}
plot(simulateResiduals(mg11))
```

```{r}
plot(simulateResiduals(mg10))
```

## Models result

The two equally plausible models for the audience included gender, academic level and affirmative actions as predictors, with the difference that the best fitted model includes an interaction of gender and academic level.

```{r echo=TRUE}
summary(mg11)
performance::r2(mg11)
```


```{r}
myg11 <- ggpredict(mg11, terms=c("position_cat","gender", "affirm_action"))
prs <- as.data.frame(myg11) %>% rename(affirm_action = facet)
colnames(prs)[1] <- "position_cat"
```


```{r, echo=F}
ggplot(data, aes(x=position_cat, y=audience_n)) +
  geom_point(aes(col=gender), position = position_dodge(0.6), alpha=0.3,
             size=3,show.legend = F) +
  facet_grid(~affirm_action, labeller = 
               as_labeller(c(before = "Before Affirmative action",
                             after =  "After Affirmative action"))) +
  scale_color_manual(values = c("#6D57CF","#FCA532")) +
  scale_fill_manual(name="Gender", values = c("#6D57CF","#FCA532")) +
  geom_pointrange(data=prs, aes(x=position_cat, y=predicted,fill=group,
                                ymax=conf.high, ymin=conf.low), alpha=1,
             position=position_dodge(0.6), size=1.2, shape=21, col="black") +
    xlab("Academic level") + ylab("Audience (N)")
#ggsave("figures/FIG_3a_alternative.jpeg", width=8, height = 4)  
```


```{r echo=TRUE}
summary(mg10)
performance::r2(mg10)
```

```{r}
myg10 <- ggpredict(mg10, terms=c("position_cat","gender", "affirm_action"))
pr10 <- as.data.frame(myg10) %>% rename(affirm_action = facet)
colnames(pr10)[1] <- "position_cat"
```

```{r Fig_S3, echo=F}
ggplot(data, aes(x=position_cat, y=audience_n)) +
  geom_point(aes(col=gender), position = position_dodge(0.6), alpha=0.3,
             size=3,show.legend = F) +
  facet_grid(~affirm_action, labeller = 
               as_labeller(c(before = "Before Affirmative action",
                             after =  "After Affirmative action"))) +
  scale_color_manual(values = c("#6D57CF","#FCA532")) +
  scale_fill_manual(name="Gender", values = c("#6D57CF","#FCA532")) +
  geom_pointrange(data=pr10, aes(x=position_cat, y=predicted,fill=group,
                                ymax=conf.high, ymin=conf.low), alpha=1,
             position=position_dodge(0.6), size=1.2, shape=21, col="black") +
    xlab("Academic level") + ylab("Audience (N)") 

ggsave("figures/FIG_S3_audience_speakers_model2.jpeg", width=8, height = 4)  
```


# Only professors - productivity metrics

Investigating if differences in productivity between male and female professors and researches are related to the audience.

We measured productivity publication metrics from Google Scholar for professors and researchers. 

Creating productivity index using PCA 1st axis from metrics.

## PCA productivity metrics

```{r}
dp <- data %>% filter(!is.na(data$total_citation_n),
                      !is.na(data$nature_index_count))
colnames(dp)[22:29] <- c("tot_cit", "h", "i10", "most_cit", "cit_cum",
                        "career_Y", "nature_count", "nature_share")
table(dp$gender, dp$affirm_action)
```

Productivity publication metrics

```{r}
pca1 <- PCA(dp[, c(22:29)], graph=F)
```


```{r, eval=F}
p1 <- fviz_pca_biplot(pca1, col.ind = dp$gender, addEllipses=TRUE,
                      col.ind.sub="none",  geom="point",
                      repel = TRUE) +
  geom_vline(xintercept = 0, linetype="dashed") + 
  geom_hline(yintercept = 0, linetype="dashed")+
  scale_color_manual(name="Gender",values = c("#6D57CF","#FCA532"))+
  scale_shape(name="Gender")+
  scale_fill_manual(name="Gender",values = c("#6D57CF","#FCA532"))+
  ggtitle("Productivity metrics") +
  xlab("PC1 (52%)") + ylab("PC2 (21%)") +
  theme_cowplot()

p1
#("figures/pca_biplot.jpeg", width=6, height = 6)  
```

For the analysis specific for professor talks (N=`r length(dp$id)`), the PCA results show that all the productivity metrics for professors were highly correlated (Figure 2B) with the first axis (52% of variance explained) while the institution indexes composed the second PCA axis (21% of variation explained).

Extracting PCA 2 first axes
```{r}
dp$pc1 <- pca1$ind$coord[,1]
dp$pc2 <- pca1$ind$coord[,2]
```

## Modeling
 
```{r}
m0 <- glmmTMB(audience_n ~ 1 + affirm_action + (1|year), data=dp, family=nbinom2)
m1 <- glmmTMB(audience_n ~ gender + affirm_action +(1|year), data=dp,
              family=nbinom2)
m2 <- glmmTMB(audience_n ~ pc1 + + affirm_action + (1|year), data=dp,
              family=nbinom2)
m3 <- glmmTMB(audience_n ~ gender + pc1 + affirm_action + (1|year), 
              data=dp, family=nbinom2)
m4 <- glmmTMB(audience_n ~ gender*pc1 + affirm_action + (1|year),
              data=dp, family=nbinom2)

AICtab(m0,m1,m2,m3,m4,
       base=T, weights=T) %>% kable(digits=2)
```

## Residual diagnostic

Using the `DHARMa` package.

The two most plausible models presented a satisfactory residual diagnostic.

```{r}
plot(simulateResiduals(m3))
plot(simulateResiduals(m4))
```

## Model results

```{r echo=TRUE}
summary(m3)
performance::r2(m3)
```

We used the first PCA axis as predictors together with gender to explain the professor's audience, and found that, as expected, audience increases with productivity index (first PCA axis) but female professors still presented on average audience *1.4* times smaller than male professors.


```{r, echo=F}
my3 <- ggpredict(m3, terms=c("pc1","gender", "affirm_action[before]")) %>% as.data.frame()
ggplot(my3, aes(x=x, y=predicted, col=group)) +
  geom_ribbon(aes(ymin=conf.low,ymax=conf.high, fill=group), alpha=0.3,
             colour = NA) +
    geom_line(size=0.5)+
  scale_color_manual(name="Gender",values = c("#6D57CF","#FCA532"))+
  scale_fill_manual(name="Gender",values = c("#6D57CF","#FCA532"))+
   theme_cowplot() + ggtitle("") +
  ylab("Audience (N)") + xlab("Productivity index (PC1 axis)")+
  geom_point(data=dp, aes(x=pc1, y=audience_n, col=gender), alpha=0.6)
```

```{r}
my3 <- ggpredict(m3, terms=c("gender"))
plot(my3)
```


```{r}
my4 <- ggpredict(m4, terms=c("pc1","gender")) %>% as.data.frame()
ggplot(my4, aes(x=x, y=predicted, col=group)) +
  geom_ribbon(aes(ymin=conf.low,ymax=conf.high, fill=group), alpha=0.3,
             colour = NA) +
    geom_line()+
  scale_color_manual(name="Gender",values = c("#6D57CF","#FCA532"))+
  scale_fill_manual(name="Gender",values = c("#6D57CF","#FCA532"))+
   theme_cowplot() + ggtitle("") +
  ylab("Audience (N)") + xlab("Productivity index (PC1 axis)")+
  geom_point(data=dp, aes(x=pc1, y=audience_n, col=gender), alpha=0.6)
#ggsave("figures/audience_professor.jpeg", width=9, height = 6)  
```

# Figure 3 audience

```{r,echo=F, fig.height=9, fig.width=9}
prs <- as.data.frame(myg11) %>% rename(affirm_action = facet,
                                        position_cat=x)
f1 <- ggplot(data, aes(x=position_cat, y=audience_n)) +
  geom_point(aes(col=gender), position = position_dodge(0.6), alpha=0.3,
             size=3,show.legend = F) +
  facet_grid(~affirm_action, labeller = 
               as_labeller(c(before = "Before Affirmative action",
                             after =  "After Affirmative action"))) +
  scale_color_manual(values = c("#6D57CF","#FCA532")) +
  #scale_x_discrete(labels = c("Before", "After"))+
  #scale_y_log10()+
  scale_fill_manual(name="Gender", values = c("#6D57CF","#FCA532")) +
  geom_pointrange(data=prs, aes(x=position_cat, y=predicted,fill=group,
                                ymax=conf.high, ymin=conf.low), alpha=1,
             position=position_dodge(0.6), size=1.2, shape=21, col="black") +
    xlab("Academic level") + ylab("Audience (N)") +  labs(tag="a)")
  
my3 <- ggpredict(m3, terms=c("pc1","gender"))
my3 <- as.data.frame(my3)

f2 <- ggplot(my3, aes(x=x, y=predicted, col=group)) +
  geom_ribbon(aes(ymin=conf.low,ymax=conf.high, fill=group), alpha=0.3,
             colour = NA) +
    geom_line(size=0.9)+
  scale_color_manual(name="Gender",values = c("#6D57CF","#FCA532"))+
  scale_fill_manual(name="Gender",values = c("#6D57CF","#FCA532"))+
   theme_cowplot() + ggtitle("") +
  ylab("Audience (N)") + xlab("Productivity index (PC1 axis)")+
  geom_point(data=dp, aes(x=pc1, y=audience_n, col=gender), alpha=0.6,
             size=2)+
  theme(legend.position="none",
         plot.title = element_text(size=12, vjust=2, hjust=0.5)) +
  labs(tag="c)", title="Professors' audience")

p1<- fviz_pca_biplot(pca1, col.ind = dp$gender, addEllipses=TRUE,
                      col.ind.sub="none",  geom="point",
                      repel = TRUE) +
  facet_grid(.~.)+
  geom_vline(xintercept = 0, linetype="dashed") + 
  geom_hline(yintercept = 0, linetype="dashed")+
  scale_color_manual(name="gender",values = c("#6D57CF","#FCA532"))+
  scale_shape(name="gender")+
  scale_fill_manual(name="gender",values = c("#6D57CF","#FCA532"))+
  labs(title="Professors' productivity", tag="b)") +
  xlab("PC1 (52%)") + ylab("PC2 (21%)") +
  theme_cowplot() +
  theme(legend.position="none",
        plot.title = element_text(size=12, vjust=2, hjust=0.5))

f1/(p1+f2) +plot_layout(guides="collect", heights = c(0.9,1))
ggsave("figures/FIG_3_audience.jpeg", width=10, height = 9)  
```


<!--T
```{r, echo=F, eval=F}
TEST
p1alt <- fviz_pca_biplot(pca1, col.ind = dp$gender, addEllipses=TRUE,
                      col.ind.sub="none",  geom="point",labelsize =3,
                      col.var = "#3B3B3B",
                      repel = TRUE) +
  facet_grid(.~.)+
  xlim(-5.8,7.5)+
  geom_vline(xintercept = 0, linetype="dashed", col="azure4") + 
  geom_hline(yintercept = 0, linetype="dashed", col="azure4")+
  scale_color_manual(name="gender",values = c("#6D57CF","#FCA532"))+
  scale_shape(name="gender")+
  scale_fill_manual(name="gender",values = c("#6D57CF","#FCA532"))+
  ggtitle("")+
  #labs(title="Professors' productivity") +
  xlab("PC1 (52%)") + ylab("PC2 (21%)") +
  theme_cowplot() +
  theme(legend.position="none",
        panel.background = element_rect(fill="white"),
        text = element_text(size=9),
        axis.text = element_text(size=9),
        plot.title = element_text(size=12, vjust=2, hjust=0.5))


f2 + ggtitle("") + labs(tag="")+ 
  theme(legend.position="right",
        plot.title = element_text(size=12, vjust=2, hjust=0.5)) +
 inset_element(p1alt,
               0.49, 0.5, 1.03, 1.14)
ggsave("figures/FIG_3_bc_alternative.jpeg", width=8, height=6)
```