1_speakers_genderPosition.Rmd

---
title: "Gender bias in speakers and career position"
author: "Melina Leite & Júlia Barreto"
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(DHARMa); library(ggeffects)
library(bbmle);
library(performance)
library(glmmTMB)
library(here)

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 statistics in script `0_data_summary`.

```{r}
load(here("data_clean", "data_proportion.Rdata"))
load(here("data_clean", "data_pop_ppge.Rdata"))
data <- data_p
```


Creating dummy column to indicate if the speaker is a female (1) or not (0)
```{r}
data$fem <- 1
data$fem[data$gender == "M"] <- 0 
```


# Modeling

Proportions of female speakers by academic level, before and after affirmative actions.

Mixed effects model, including year as random intercept.

Binomial distribution.

```{r, echo=T}
mod0 <- glmmTMB(fem ~ 1 + (1|year), 
              family=binomial, data = data)
mod1 <- glmmTMB(fem ~ affirm_action  + (1|year), 
              family=binomial, data = data)
mod2 <- glmmTMB(fem ~ position_cat + (1|year), 
              family=binomial, data = data)
mod3 <- glmmTMB(fem ~ position_cat + affirm_action + (1|year), 
              family=binomial, data = data)
mod4 <- glmmTMB(fem ~ position_cat*affirm_action + (1|year), 
              family=binomial, data = data)

kable(AICtab(mod0,mod1,mod2,mod3,mod4, base=T, weights=T), digits=2)
```

## Residual diagnostic of the selected models

Using the `DHARMa` package.

The two most plausible models presented a satisfactory residual diagnostic.

```{r}
plot(simulateResiduals(mod4))
plot(simulateResiduals(mod2))
```

## Models results

Predicting the proportion of female speakers fixing the population gender ration at 1 (1:1).

```{r}
summary(mod4)
performance::r2(mod4)
my4 <- ggpredict(mod4, terms=c("position_cat","affirm_action"))
```


```{r fig prop, echo=F}
suma <- data %>% count(position_cat, affirm_action,fem) %>%
  mutate(position_cat = fct_relevel(position_cat, "student", "postdoc", "professor"))
prs <- as.data.frame(my4)

ggplot(suma, aes(x = position_cat, y = fem,col = affirm_action))+
  geom_point(aes(size=n), position = position_dodge(0.6), alpha = 0.2, show_guides = F) +
  scale_size(range = c(1,10), breaks = c(3,10,20,60)) +
  geom_pointrange(data = prs, aes(x = x, y = predicted, col = group,
                                ymax = conf.high, ymin = conf.low), 
             position = position_dodge(0.6)) +
  geom_hline(yintercept = 0.5, linetype = "dashed", col = 'gray') +
  scale_color_manual(name="Affirmative \n actions", 
                     values = c("goldenrod", "green4")) +
  ylab("Proportion of female speakers")+
  xlab("Academic level") +
  theme(text = element_text(size=20),
        axis.text = element_text(size=18))
ggsave(here("figures", "FIG_2_prop_female_speakers_model.jpeg"), width=9, height = 6)  
```

```{r}
summary(mod2)
performance::r2(mod2)
my2 <- ggpredict(mod2, terms=c( "position_cat"))
```


```{r, echo=F}
prs2 <- as.data.frame(my2)
ggplot(prs2, aes(x=x, y=predicted, ymin=conf.low, ymax=conf.high)) +
  geom_pointrange()+
   geom_hline(yintercept = 0.5, linetype="dashed", col='gray') +
  ylab("Proportion of female speakers")+
  xlab("Academic level")+ ylim(0,1)+
  theme(text = element_text(size=20),
        axis.text = element_text(size=18))
ggsave(here("figures","FIG_S3_prop_female_speakers_model2.jpeg"), width=9, height = 6)  
```

# Extra analysis: speakers of the PPGE community

Subset analysis with only speakers of the PPGE community in order to compare the proportion of female speakers with the proportion of female in the population source.

Data form the population of speakers is the proporion of female in each academic level and each year in the PPGE community.

Including the information of the PPGE community in the dataset
```{r}
data <- data %>% left_join(pop2[,c(1:2,9:10)], 
                           by=c("year", "position_cat"="category")) %>%
  mutate(position_cat = fct_relevel(position_cat, "student","postdoc",
                                    "professor"))
```

Selecting only speakers of the PPGE community.
```{r}
datappge <- data %>% filter(department == "Ecologia", institute == "IB",
                            university == "USP")
dim(datappge)
```

```{r}
table(datappge$position_cat)
table(datappge$affirm_action)
```



## Modeling

```{r}
mnullppge<- glmmTMB(fem ~ 1 + (1|year), 
                  family = binomial, data = datappge)
mod0ppge <- glmmTMB(fem ~ 1 + propFcat + (1|year), 
                  family = binomial, data = datappge)
mod1ppge <- glmmTMB(fem ~ affirm_action + propFcat + (1|year), 
                  family = binomial, data = datappge)
mod2ppge <- glmmTMB(fem ~ position_cat + propFcat + (1|year), 
                  family = binomial, data = datappge)
mod3ppge <- glmmTMB(fem ~ position_cat + affirm_action + propFcat + (1|year), 
                  family = binomial, data = datappge)
kable(AICtab(mod0ppge, mod1ppge, mod2ppge, mod3ppge, mnullppge,
             base = T, weights = T), digits = 2)
```



Residual diagnostic of 2 the best-fitting models.

```{r}
plot(simulateResiduals(mod0ppge))
plot(simulateResiduals(mod2ppge))
```


## Model's Predictions

## Model with only the population-level proportion of females

Solid line indicates the model's prediction.
Dashed line indicates the proportions 50% of females for each proportion - relationship 1:1. 
Dotted horizontal line indicates the 50% of speakers.

```{r, fig.height=6, fig.width=6}
plot(ggpredict(mod0ppge, terms=c("propFcat[all]"))) +
    theme_cowplot() + ggtitle("") +
  geom_abline(slope=1, intercept = 0, linetype="dashed") +
  geom_hline(yintercept = 0.5, linetype="dotted")+
  xlab("Proportion of female academics in the PPGE community") +
  ylab("Proportion of female speakers")
ggsave(here("figures", "FIG_S2_propFemalePPGE.jpeg"), width=9, height = 6)
```

## Model with academic level and the population-level gender ratio

```{r, fig.height=6, fig.width=6}
plot(ggpredict(mod2ppge, terms=c("propFcat[all]", "position_cat"))) +
  theme_cowplot() + ggtitle("") +
  geom_hline(yintercept = 0.5, linetype="dotted") +
   geom_abline(slope=1, intercept = 0, linetype="dashed", alpha=0.4) +
  xlab("Proportion of female academics in the PPGE community") +
  ylab("Proportion of female speakers")
#ggsave(here("figures", "FIG_S2_propFemalePPGE.jpeg"), width=9, height = 6)
```