---
title: "MH survey analysis"
output: 
  html_document:
    toc: true
    toc_float: true
    toc_depth: 3
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(MASS)
library(dplyr)
library(ggplot2)
library("nnet")
library("R.devices")
library("MHTdiscrete")
```

## Abbreviations used in the markdown/comments of this script
- LR: Logistic regression
- FLR: Factorial Logistic regreesion
- OLR: Ordinal logistic regression
- KW: Kruskal-Wallis test
- OR: Odds ratio
- CI: Confidence interval (two-sided 95% confidence)
- glm: Generalised linear model
- CS: career stage
- PNTA: Prefer not to answer


## A quick introduction to the statistical tests and choices
This section intends to give a high level explanation to the code that follows. The information provided in this section does NOT supercede the methods section of the report - please always regard to the information provided in the report as final.

### Independent variables
We checked mainly for the effect of 2 independent variables on the responses:
- Gender - men vs women (variable name: All_Gender_clean)
- Career stages - PhDs vs Postdocs vs Early career PIs vs mid-career PIs vs late-career PIs (variable: Supporter_CareerStage_clean)
We had respondents of other genders and career stages, but for this report, we chose to only analyse these groups because they had a reasonable and comparable number of samples. (If you'd like to look into other groups specifically, you're more than welcomed!)

For some questions we also checked the effect of whether someone is a supporter on their response (variable name: Provided_Support_clean). Supporter status is determined based on their response to the statemnet "I have provided support to someone who was doing research and who was struggling with their mental health". We compared respondents who answered "No" and one of the "yes" options to this statement.
For some questions we also wanted to see if early career PI stands up from other career stages. For this, we run either (1) Kruskal-Wallis test between early PI and others pooled together; or (2) pairwise comparison between early PIs and all other groups. In the second case, Sidak correction was applied posthoc. 
### Dependent variables, i.e. the response
A wide variety of questions were asked - we focussed our statistical analyses on two types of responses (or responses that could be logically recoded to one of the following two types):
- Categorical, and in most cases, binary, i.e. Yes/No
- Ordinal, i.e. Likert scale type responses

For categorical responses, we modelled response using (binomial) logistic regression. We ran 3 seperate models - one with genders as predictors, the second with career stages (CS) as a predictors, and the final with both as predictors. For gender, men is always the baseline indicator variable, and for CS, PhD is always the baseline indicator variable. The first two models provide the reported p-values, odds ratios (OR, obtained with exp(coef()) and confidence intervals of the OR (two-sided 95% confidenece, obtained with exp(confint())). This final factorial logistic regression allows us to check for interactions between gender and CS.

For ordinal responses, we modelled responses using ordinal logistic regression (OLR). We ran 2 seperate models - one with genders as predictors, the second with career stages (CS) as a predictors. For gender, men is always the baseline indicator variable, and for CS, PhD is always the baseline indicator variable. These give the ORs and CIs reported. As the ordinal responses we've obtained are usually not normally distributed (code not reported here, but normality was tested using Shapiro.test()), we chose to use the Kruskal-Wallis test to test for the effect of gender/CS on the ordinal responses. To check for any interactions between the two independent variables we also split the data into groups by one independent variable and used the chi-square test to check for the effect of the other independent variable (i.e. within all women, is there an effect of CS, etc). 

## Initialising
Load and inspect data structure

```{r load_data}
data=read.csv("../mh-data/cleandata2604 REDACTED.csv")
# colnames(data)
```

##Filtering the data
Looking at only 2 independent variables:

1. Gender (label: "All_Gender_clean"): Man -1, Woman -3

2. Career stage (label "Supporter_CareerStage_clean"): 
  * A PhD student-3; 
  * A postdoctoral researcher (postdoc) - 4; 
  * A group leader (less than 5 years’ experience)- 6
  * A group leader (5 to 10 years’ experience) -7
  * A group leader (10 or more years’ experience) -8
  
```{r filter_data}
genderCS=data[data$All_Gender_clean %in% c(1,3) & data$Supporter_CareerStage_clean %in% c(3,4,6,7,8), ]
genderCS$All_Gender_clean=as.factor(genderCS$All_Gender_clean)
genderCS$Supporter_CareerStage_clean=as.factor(genderCS$Supporter_CareerStage_clean)
dim(genderCS)

levels(genderCS$All_Gender_clean)=c("Men", "Women")
levels(genderCS$Supporter_CareerStage_clean)=c("PhD students", "Postdocs", "Group leaders (<5yr)", "Group leaders (5-10yr)", "Group leaders (>10yr)")
```

## the eLife colour palette

```{r}
eGreen="#346A2D"
eLime="#7DB441"
eBlue="#06589C"
eSky="#2997D4"
ePurple="#881350"
eFuschia="#D81F62"
eGrey="#666B6E"

ePalette=c(eGreen, eLime, eGrey, eSky, eBlue)
```

## Q28xQ04 Gender x Career stage (CS) (Section 2, Panel 1)
chi-square
```{r gender_CS}
gender=genderCS[,c("Supporter_CareerStage_clean","All_Gender_clean")]
dim(gender)[1] #this gives N
table(gender)
chisq.test(table(gender))
```

### Q28xQ04 Gender x Career stage (CS) (Graph 2.1.2)

```{r}

gender$Supporter_CareerStage_clean=factor(gender$Supporter_CareerStage_clean) #be careful when rerunning this line

# This als all other eps images will be output to your code directory inside the images sub-folder
#eps("images/CareerStage_gender.eps", width=1000, height=578)
graphdata = gender %>%
  group_by(Supporter_CareerStage_clean, All_Gender_clean) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))
ggplot(graphdata, aes(x=Supporter_CareerStage_clean, y=perc, fill=All_Gender_clean)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=20))+
  geom_text(aes(label=round(perc, digit=1)), size=8, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", size=3, title="During my supporting role, I was:") +
    guides(fill=guide_legend(reverse=TRUE)) +
  guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=c(eFuschia, ePurple))
ggsave("2-1-2.png")
```

## Q01 Providing support
Dependent variable: Provided_Support

* 1- No
* 2 - Yes, 1 individual
* 3 - Yes, 1-5
* 4- Yes, >5

### Treat Dependent variable as ordinal

```{r}
Provided_Support=genderCS[,c("Provided_Support", "Supporter_CareerStage_clean", "All_Gender_clean")]
Provided_Support=Provided_Support[Provided_Support$Provided_Support %in% seq(1,4), ] #remove NAs and some choices
dim(Provided_Support)[1] #N
Provided_Support$Provided_Support=as.factor(Provided_Support$Provided_Support)
table(Provided_Support)
```

#### Providing support x gender
```{r}
exp(coef(polr(Provided_Support~All_Gender_clean, data=Provided_Support, Hess=TRUE, method="logistic")))
exp(confint(polr(Provided_Support~All_Gender_clean, data=Provided_Support, Hess=TRUE, method="logistic")))
kruskal.test(Provided_Support ~ All_Gender_clean, data=Provided_Support)
```

#### Providing support x Career stage
```{r}
exp(coef(polr(Provided_Support~Supporter_CareerStage_clean, data=Provided_Support, Hess=TRUE, method="logistic")))
exp(confint(polr(Provided_Support~Supporter_CareerStage_clean, data=Provided_Support, Hess=TRUE, method="logistic")))
kruskal.test(Provided_Support~Supporter_CareerStage_clean, data=Provided_Support)
```

#### chisq for intersections [NOT REPORTED -- EMMY]

```{r}
#split by career stages, gender effect
chisq.test(table(Provided_Support[Provided_Support$Supporter_CareerStage_clean==levels(Provided_Support$Supporter_CareerStage_clean)[1],c("Provided_Support","All_Gender_clean")]))
chisq.test(table(Provided_Support[Provided_Support$Supporter_CareerStage_clean==levels(Provided_Support$Supporter_CareerStage_clean)[2],c("Provided_Support","All_Gender_clean")]))
chisq.test(table(Provided_Support[Provided_Support$Supporter_CareerStage_clean==levels(Provided_Support$Supporter_CareerStage_clean)[3],c("Provided_Support","All_Gender_clean")]))
chisq.test(table(Provided_Support[Provided_Support$Supporter_CareerStage_clean==levels(Provided_Support$Supporter_CareerStage_clean)[4],c("Provided_Support","All_Gender_clean")]))
chisq.test(table(Provided_Support[Provided_Support$Supporter_CareerStage_clean==levels(Provided_Support$Supporter_CareerStage_clean)[5],c("Provided_Support","All_Gender_clean")]))

#split by gender, career effect (same as KW above)
chisq.test(table(Provided_Support[Provided_Support$All_Gender_clean==levels(Provided_Support$All_Gender_clean)[1],c("Provided_Support","Supporter_CareerStage_clean")]))
chisq.test(table(Provided_Support[Provided_Support$All_Gender_clean==levels(Provided_Support$All_Gender_clean)[2],c("Provided_Support","Supporter_CareerStage_clean")]))
```

## Q02 Not_ProvidedSupport (Section 2, Panel 2)
Dependent variable:Not_ProvidedSupport

* 1: No
* 2: Yes, once
* 3: Yes, a few times
* 4: Yes, often
* 5: I'm not sure -removed
* 6: PNTA -removed

```{r Not_ProvidedSupport}
NotProvidedSupport=genderCS[,c("Not_ProvidedSupport", "All_Gender_clean", "Supporter_CareerStage_clean")]
summary(NotProvidedSupport) #there are some NAs
#clean
NotProvidedSupport=NotProvidedSupport[(NotProvidedSupport$Not_ProvidedSupport %in% seq(1:4)),]
dim(NotProvidedSupport)[1] #this gives N
#factorise
NotProvidedSupport[NotProvidedSupport$Not_ProvidedSupport==1, "Not_ProvidedSupport"]="No"
NotProvidedSupport[NotProvidedSupport$Not_ProvidedSupport==2, "Not_ProvidedSupport"]="Yes, once"
NotProvidedSupport[NotProvidedSupport$Not_ProvidedSupport==3, "Not_ProvidedSupport"]="Yes, a few times"
NotProvidedSupport[NotProvidedSupport$Not_ProvidedSupport==4, "Not_ProvidedSupport"]="Yes, often"
NotProvidedSupport$Not_ProvidedSupport=factor(NotProvidedSupport$Not_ProvidedSupport, levels=c("Yes, often", "Yes, a few times", "Yes, once",  "No"))
```

### Not provided support x gender graph
```{r}
graphdata = NotProvidedSupport %>%
  group_by(All_Gender_clean, Not_ProvidedSupport) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

#eps("images/NotProvidedSupport_gender.eps", width=1000, height=578) #[NOT REPORTED]
ggplot(graphdata, aes(x=All_Gender_clean, y=perc, fill=Not_ProvidedSupport)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=20))+
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="I have been in a situation when I thought someone in academia needed help, but I couldn’t or didn’t provide
support") +
    scale_x_discrete(limits = rev(levels(graphdata$Supporter_CareerStage_clean))) +
  guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
dev.off()
```

### Not provided support x CS graph (Graph 2.2.2)
```{r}
graphdata = NotProvidedSupport %>%
  group_by(Supporter_CareerStage_clean,Not_ProvidedSupport) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))
# This will save the image to your local code folder
#eps("images/NotProvidedSupport_CS.eps", width=1000, height=578)
ggplot(graphdata, aes(x=Supporter_CareerStage_clean, y=perc, fill=Not_ProvidedSupport)) +
  geom_bar(stat="identity") +
  theme_minimal() +
      theme(plot.title = element_text(hjust = 0.5), text=element_text(size=10), legend.position="bottom", legend.text = element_text(size=7))+
  geom_text(aes(label=round(perc, digit=1)), size=2.5, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="I have been in a situation when I thought someone in academia needed help, \nbut I couldn’t or didn’t provide
support") +
  guides(fill=guide_legend(reverse=TRUE)) +
  #coord_flip() +
  scale_fill_manual(name="", values=ePalette)
ggsave(dpi=700, "2-2-2-v3.png")
dev.off()

```

### Stats (Section 2, panel 2)
recode to Yes - 1 / No-0
Model: FLR
```{r}
# recoding response to binary variable
NotProvidedSupport$Not_ProvidedSupport_binary=NotProvidedSupport$Not_ProvidedSupport
levels(NotProvidedSupport$Not_ProvidedSupport_binary)=c(1, 1, 1, 0)
NotProvidedSupport$Not_ProvidedSupport_binary=as.numeric(NotProvidedSupport$Not_ProvidedSupport_binary)
NotProvidedSupport[NotProvidedSupport$Not_ProvidedSupport_binary=="2", "Not_ProvidedSupport_binary"]=0

# LR intersection between gender and CS
NotProvidedSupport_glm=glm(Not_ProvidedSupport_binary~Supporter_CareerStage_clean*All_Gender_clean, data=NotProvidedSupport, family="binomial")
summary(NotProvidedSupport_glm)

# LR gender
summary(glm(Not_ProvidedSupport_binary~All_Gender_clean, data=NotProvidedSupport, family="binomial")) #glm summary with p-values
exp(coef(glm(Not_ProvidedSupport_binary~All_Gender_clean, data=NotProvidedSupport, family="binomial"))) #OR (men as baseline)
exp(confint(glm(Not_ProvidedSupport_binary~All_Gender_clean, data=NotProvidedSupport, family="binomial"))) #CI

# LR Career stage
summary(glm(Not_ProvidedSupport_binary~Supporter_CareerStage_clean, data=NotProvidedSupport, family="binomial")) #glm summary with p-values
exp(coef(glm(Not_ProvidedSupport_binary~Supporter_CareerStage_clean, data=NotProvidedSupport, family="binomial"))) #OR (PhD as baseline)
exp(confint(glm(Not_ProvidedSupport_binary~Supporter_CareerStage_clean, data=NotProvidedSupport, family="binomial"))) #CI
```

## Q05 How often do they support someone who is officially their responsibility (Section 3, Panel 1)

Dependent variable: Supporter_Officially

* 1: No - not my official job
* 2: Yes - my official job

(remove 3- not sure, and 4 - PNTA)

Model: First individual factors (Gender, CS), then Factorial logistical regression Supporter_Officially ~ Supporter_CareerStage * All_Gender_clean) to check for interactions

```{r Supporter_Officially}
#subsetting data and recording response to binary
Officially=genderCS[,c("Supporter_Officially", "Supporter_CareerStage_clean", "All_Gender_clean")]
Officially=Officially[Officially$Supporter_Officially %in% seq(1,2), ] #remove NAs and some choices
dim(Officially)[1]
Officially[Officially$Supporter_Officially==1, "Supporter_Officially"]=0 #recoding No->0 and Yes->1 for binom modelling
Officially[Officially$Supporter_Officially==2, "Supporter_Officially"]=1

# LR gender
summary(glm(Supporter_Officially~All_Gender_clean, data=Officially, family="binomial")) #glm and p-value
exp(coef(glm(Supporter_Officially~All_Gender_clean, data=Officially, family="binomial"))) #OR (men as baseline)
exp(confint(glm(Supporter_Officially~All_Gender_clean, data=Officially, family="binomial"))) #CI

# LR CS
summary(glm(Supporter_Officially~Supporter_CareerStage_clean, data=Officially, family="binomial")) #glm and p-value
exp(coef(glm(Supporter_Officially~Supporter_CareerStage_clean, data=Officially, family="binomial"))) #OR (PhD as baseline)
exp(confint(glm(Supporter_Officially~Supporter_CareerStage_clean, data=Officially, family="binomial"))) #CI

# FLR gender-CS interactions
summary(glm(Supporter_Officially~Supporter_CareerStage_clean*All_Gender_clean, data=Officially, family="binomial"))
```

### How often do they support someone who is officially their responsibility x gender graph (Graph 3.1.2)
```{r}
Officially[Officially$Supporter_Officially==0, "Supporter_Officially"]="No"
Officially[Officially$Supporter_Officially==1, "Supporter_Officially"]="Yes"
Officially$Supporter_Officially=factor(Officially$Supporter_Officially)
#ImpactPersonal$Impact_Personal=factor(ImpactPersonal$Impact_Personal, levels=c("No","Yes"))
graphdata =Officially %>%
  group_by(All_Gender_clean, Supporter_Officially) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/Officially_gender.eps", width=1000, height=578)
ggplot(graphdata, aes(x=All_Gender_clean, y=perc, fill=Supporter_Officially)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16))+
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="It was part of my official job description to be in charge of the person I helped") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=c(eBlue, eGreen))
ggsave(dpi=1000, "3-1-2.png", limitsize = FALSE)
dev.off()
```
### How often do they support someone who is officially their responsibility x career graph (Graph 3.1.1)
```{r}
graphdata = Officially %>%
  group_by(Supporter_CareerStage_clean, Supporter_Officially) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/Officially_CS.eps", width=1000, height=578)
ggplot(graphdata, aes(x=Supporter_CareerStage_clean, y=perc, fill=Supporter_Officially,)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16))+
  geom_text(aes(label=round(perc, digit=1)), size=3.5, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="It was part of my official job description \nto be in charge of the person I helped") +
    guides(fill=guide_legend(reverse=TRUE)) +
  scale_x_discrete(limits = rev(levels(graphdata$Supporter_CareerStage_clean))) +
  coord_flip() +
  scale_fill_manual(name="", values=c(eBlue, eGreen))
ggsave(dpi=1000, "3-1-1.png", limitsize = FALSE)
dev.off()
```

## Q06 MH at the time of support (Section 2, Panel 3)

Response (label: "Supporter_MH"): No-1 , Yes-2 (Not sure-3, and PNTA-4 are removed for binomial modelling)

Model: Factorial logistical regression Supporter_MH ~ Supporter_CareerStage * All_Gender_clean)

```{r Supporter_MH}
# Subsetting data and recoding response to binary
MH=genderCS[,c("Supporter_MH", "Supporter_CareerStage_clean", "All_Gender_clean")]
MH=MH[MH$Supporter_MH %in% seq(1,2), ] #remove NAs and some choices
dim(MH)[1] #this gives N
MH[MH$Supporter_MH==1, "Supporter_MH"]=0 #recoding No->0 and Yes->1 for binom modelling
MH[MH$Supporter_MH==2, "Supporter_MH"]=1
MH_glm=glm(Supporter_MH~Supporter_CareerStage_clean*All_Gender_clean, data=MH, family="binomial")
summary(MH_glm)

# LR gender
summary(glm(Supporter_MH~All_Gender_clean, data=MH, family="binomial")) #glm and p-values
exp(coef(glm(Supporter_MH~All_Gender_clean, data=MH, family="binomial"))) #OR (men as baseline)
exp(confint(glm(Supporter_MH~All_Gender_clean, data=MH, family="binomial"))) #CI

# LR CS
summary(glm(Supporter_MH~Supporter_CareerStage_clean, data=MH, family="binomial")) #glm and p-values
exp(coef(glm(Supporter_MH~Supporter_CareerStage_clean, data=MH, family="binomial"))) #OR (PhDs as baseline)
exp(confint(glm(Supporter_MH~Supporter_CareerStage_clean, data=MH, family="binomial"))) #CI
```

### MH at the time of support x gender graph (graph 2.3.3)

```{r}
# This will save the image to your local code folder
#eps("images/SupporterMH_gender.eps", width=1000, height=578)
MH[MH$Supporter_MH==0, "Supporter_MH"]="No"
MH[MH$Supporter_MH==1, "Supporter_MH"]="Yes"
MH$Supporter_MH=factor(MH$Supporter_MH)

graphdata = MH %>%
  group_by(All_Gender_clean, Supporter_MH) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))
ggplot(graphdata, aes(x=All_Gender_clean, y=perc, fill=Supporter_MH)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16))+
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="As I was providing support, I was struggling with my own mental health:") +
  guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=c(eBlue, eGreen))
ggsave(dpi=1000, "2-3-3-v2.png")
dev.off()
```

### MH at the time of support x CS graph (Graph 2.3.5)

```{r}
# This will save the image to your local code folder
#eps("images/SupporterMH_CS.eps", width=1000, height=578)
graphdata = MH %>%
  group_by(Supporter_CareerStage_clean, Supporter_MH) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))
ggplot(graphdata, aes(x=Supporter_CareerStage_clean, y=perc, fill=Supporter_MH)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16))+
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="As I was providing support, I was struggling with my own mental health") +
    guides(fill=guide_legend(reverse=TRUE)) +
  scale_x_discrete(limits = rev(levels(graphdata$Supporter_CareerStage_clean))) +
  coord_flip() +
  scale_fill_manual(name="", values=c(eBlue, eGreen))
ggsave(dpi=1000, "2-3-5-v1.png")
dev.off()
```

## Q07 Helping more than one person at the time (Section 2, Panel 3)

Dependent variable: Supporter_NumberReceivers

* 1: No (just one person)
* 2: Yes (more than one person)

Model: First individual factors (Gender, CS), then Factorial logistical regression Supporter_NumberReceivers ~ Supporter_CareerStage * All_Gender_clean) to check for interactions

```{r Supporter_NumberReceivers}
# Subsetting data and recoding response to binary
NumberReceivers=genderCS[,c("Supporter_NumberReceivers", "Supporter_CareerStage_clean", "All_Gender_clean")]
NumberReceivers=NumberReceivers[NumberReceivers$Supporter_NumberReceivers %in% seq(1,2), ] #remove NAs and some choices
dim(NumberReceivers)[1] #N
NumberReceivers[NumberReceivers$Supporter_NumberReceivers==1, "Supporter_NumberReceivers"]=0 #recoding No->0 and Yes->1 for binom modelling
NumberReceivers[NumberReceivers$Supporter_NumberReceivers==2, "Supporter_NumberReceivers"]=1

#LR gender
summary(glm(Supporter_NumberReceivers~All_Gender_clean, data=NumberReceivers, family="binomial")) #glm and p-value
exp(coef(glm(Supporter_NumberReceivers~All_Gender_clean, data=NumberReceivers, family="binomial"))) #OR (men as baseline)
exp(confint(glm(Supporter_NumberReceivers~All_Gender_clean, data=NumberReceivers, family="binomial"))) #CI

#LR CS
summary(glm(Supporter_NumberReceivers~Supporter_CareerStage_clean, data=NumberReceivers, family="binomial")) #glm and p-value
exp(coef(glm(Supporter_NumberReceivers~Supporter_CareerStage_clean, data=NumberReceivers, family="binomial"))) #OR (PhD as baseline)
exp(confint(glm(Supporter_NumberReceivers~Supporter_CareerStage_clean, data=NumberReceivers, family="binomial"))) #CI

#FLR gender-CS interactions
summary(glm(Supporter_NumberReceivers~Supporter_CareerStage_clean*All_Gender_clean, data=NumberReceivers, family="binomial"))
```
### graph for Helping more than one person at the time x gender (Graph 2.3.4)
```{r}
NumberReceivers[NumberReceivers$Supporter_NumberReceivers==0, "Supporter_NumberReceivers"]="No"
NumberReceivers[NumberReceivers$Supporter_NumberReceivers==1, "Supporter_NumberReceivers"]="Yes"
NumberReceivers$Supporter_NumberReceivers=factor(NumberReceivers$Supporter_NumberReceivers)
#ImpactPersonal$Impact_Personal=factor(ImpactPersonal$Impact_Personal, levels=c("No","Yes"))

# This will save the image to your local code folder
#eps("images/NumberReceivers_gender.eps", width=1000, height=578)
graphdata = NumberReceivers %>%
  group_by(All_Gender_clean, Supporter_NumberReceivers) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))
ggplot(graphdata, aes(x=All_Gender_clean, y=perc, fill=Supporter_NumberReceivers)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16))+
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="I was helping more than one person \nat the same time (in my professional and personal life)") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=c(eBlue, eGreen))
ggsave(dpi=1000, "2-3-4.png", limitsize=FALSE)
dev.off()
```

### graph for Helping more than one person at the time x CS (Graph 2.3.6)
```{r}
graphdata = NumberReceivers %>%
  group_by(Supporter_CareerStage_clean, Supporter_NumberReceivers) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/NumberReceivers_CS.eps", width=1000, height=578)
ggplot(graphdata, aes(x=Supporter_CareerStage_clean, y=perc, fill=Supporter_NumberReceivers)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16))+
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="I was helping more than one person at the same time \n(in my professional and personal life)") +
    guides(fill=guide_legend(reverse=TRUE)) +
  scale_x_discrete(limits = rev(levels(graphdata$Supporter_CareerStage_clean))) +
  
  coord_flip() +
  scale_fill_manual(name="", values=c(eBlue, eGreen))
ggsave(dpi=1000, "2-3-6.png", limitsize = FALSE)
dev.off()
```

## Q13 Missing resources for person supported (Section 3, Panel 4)

Filter down to people who picked at least 1 of the 6 options:

* Receiver_Missing_Personal
* Receiver_Missing_Institutions
* Receiver_Missing_Peers
* Receiver_Missing_Managers
* Receiver_Missing_Pro
* Receiver_Missing_OutPro

```{r}
#subsetting data and reocoding variables
Receiver_Missing=genderCS[(genderCS$Receiver_Missing_Personal==1 | genderCS$Receiver_Missing_Institutions==1 | genderCS$Receiver_Missing_Peers==1 | genderCS$ Receiver_Missing_Managers==1 | genderCS$Receiver_Missing_Pro==1 | genderCS$Receiver_Missing_OutPro==1), ]
Receiver_Missing=genderCS[,c("Receiver_Missing_Personal", "Receiver_Missing_Institutions", "Receiver_Missing_Peers", "Receiver_Missing_Managers", "Receiver_Missing_Pro", "Receiver_Missing_OutPro", "Supporter_CareerStage_clean","All_Gender_clean")]
dim(Receiver_Missing)[1]
#fill all na with 0
Receiver_Missing[is.na(Receiver_Missing)]=0
```

### Receiver_Missing_Personal (Yes/No) xCS xGender xCSxGender
Dependent variable - 

* 0 - didn't select this option 
* 1 - did select this option - i.e. Missing personal support

Model: Factorial logistical regression

```{r}
summary(glm(Receiver_Missing_Personal~All_Gender_clean, data=Receiver_Missing, family="binomial")) #LR gender
summary(glm(Receiver_Missing_Personal~Supporter_CareerStage_clean, data=Receiver_Missing, family="binomial")) #LR CS
summary(glm(Receiver_Missing_Personal~Supporter_CareerStage_clean*All_Gender_clean, data=Receiver_Missing, family="binomial")) #LR gender-CS interactions
```

### Receiver_Missing_Managers xCS xGender xCSxGender
Missing: Support from supervisors and managers
as above

```{r}
#LR gender
summary(glm(Receiver_Missing_Managers~All_Gender_clean, data=Receiver_Missing, family="binomial")) #glm and p-value
exp(coef(glm(Receiver_Missing_Managers~All_Gender_clean, data=Receiver_Missing, family="binomial"))) #OR (men as baseline)
exp(confint(glm(Receiver_Missing_Managers~All_Gender_clean, data=Receiver_Missing, family="binomial"))) #CI

#LR CS
summary(glm(Receiver_Missing_Managers~Supporter_CareerStage_clean, data=Receiver_Missing, family="binomial")) #glm and p-value
exp(coef(glm(Receiver_Missing_Managers~Supporter_CareerStage_clean, data=Receiver_Missing, family="binomial"))) #OR (PhD as baseline)
exp(confint(glm(Receiver_Missing_Managers~Supporter_CareerStage_clean, data=Receiver_Missing, family="binomial"))) #CI

#FLR gender x CS interactions
summary(glm(Receiver_Missing_Managers~Supporter_CareerStage_clean*All_Gender_clean, data=Receiver_Missing, family="binomial"))

```

### Receiver_Missing_Peers xCS xGender xCSxGender

Missing: Support from peers and colleagues
as above

```{r}
summary(glm(Receiver_Missing_Peers~All_Gender_clean, data=Receiver_Missing, family="binomial")) #LR Gender
summary(glm(Receiver_Missing_Peers~Supporter_CareerStage_clean, data=Receiver_Missing, family="binomial")) #LR CS
summary(glm(Receiver_Missing_Peers~Supporter_CareerStage_clean*All_Gender_clean, data=Receiver_Missing, family="binomial")) #LR gender-CS interactions
```

### Receiver_Missing_Institutions xCS xGender xCSxGender

Missing Institutional measures (e.g. extensions, intermission)
as above

```{r}
#LR gender
summary(glm(Receiver_Missing_Institutions~All_Gender_clean, data=Receiver_Missing, family="binomial")) #glm and p-value
exp(coef(glm(Receiver_Missing_Institutions~All_Gender_clean, data=Receiver_Missing, family="binomial"))) #OR (men as baseline)
exp(confint(glm(Receiver_Missing_Institutions~All_Gender_clean, data=Receiver_Missing, family="binomial"))) #CI

#LR CS
summary(glm(Receiver_Missing_Institutions~Supporter_CareerStage_clean, data=Receiver_Missing, family="binomial")) #glm and p-value
exp(coef(glm(Receiver_Missing_Institutions~Supporter_CareerStage_clean, data=Receiver_Missing, family="binomial"))) #OR (PhD as baseline)
exp(confint(glm(Receiver_Missing_Institutions~Supporter_CareerStage_clean, data=Receiver_Missing, family="binomial"))) #CI

#FLR gender-CS interactions
summary(glm(Receiver_Missing_Institutions~Supporter_CareerStage_clean*All_Gender_clean, data=Receiver_Missing, family="binomial"))
```

### Receiver_Missing_Pro xCS xGender xCSxGender

Missing Support from mental health professionals within the institution
as above

```{r}
summary(glm(Receiver_Missing_Pro~All_Gender_clean, data=Receiver_Missing, family="binomial")) #LR Gender
summary(glm(Receiver_Missing_Pro~Supporter_CareerStage_clean, data=Receiver_Missing, family="binomial")) #LR CS
summary(glm(Receiver_Missing_Pro~Supporter_CareerStage_clean*All_Gender_clean, data=Receiver_Missing,  family="binomial")) #FLR gender-CS interactions
```

### Receiver_Missing_OutPro xCS xGender xCSxGender

Missing Support from health professionals outside of the institution
as above

```{r}
summary(glm(Receiver_Missing_OutPro~All_Gender_clean, data=Receiver_Missing, family="binomial")) #LR gender
summary(glm(Receiver_Missing_OutPro~Supporter_CareerStage_clean, data=Receiver_Missing, family="binomial")) #LR CS
summary(glm(Receiver_Missing_OutPro~Supporter_CareerStage_clean*All_Gender_clean, data=Receiver_Missing, family="binomial")) #FLR gender-CS interactions
```
## Q14 Length of support (Section 3, panel 2)
Dependent variable: Lenght_Support (note the typo)
Recode

* 0- less than six months
  * 1: Less than six months
  * 4: I'm still supporting, <6 mo
* 1- more than six months
  * 2: 6mo-1yr
  * 3: >1yr
  * 5: I'm still supporting, >6 mo
(ignore 6- not sure and 7- PNTA)

Model: Factorial Logistical regression: Length_Support ~ Gender * CS
```{r}
#subsetting data and recoding response to binary
Length_Support=genderCS[,c("Lenght_Support", "Supporter_CareerStage_clean", "All_Gender_clean")]
colnames(Length_Support)[1]="Length_Support"
Length_Support=Length_Support[Length_Support$Length_Support %in% seq(1,5), ] #remove NAs and some choices
dim(Length_Support)[1] #this gives N
#recode
Length_Support[Length_Support$Length_Support==1, "Length_Support"]=0 
Length_Support[Length_Support$Length_Support==4, "Length_Support"]=0 
Length_Support[Length_Support$Length_Support==2, "Length_Support"]=1
Length_Support[Length_Support$Length_Support==3, "Length_Support"]=1
Length_Support[Length_Support$Length_Support==5, "Length_Support"]=1 

#LR gender
summary(glm(Length_Support~All_Gender_clean, data=Length_Support, family="binomial")) #glm and p-value

#LR CS
summary(glm(Length_Support~Supporter_CareerStage_clean, data=Length_Support, family="binomial")) #glm and p-value

#FLR gender-CS interactions
summary(glm(Length_Support~Supporter_CareerStage_clean*All_Gender_clean, data=Length_Support, family="binomial")) 

#effect of gender on response
kruskal.test(Length_Support ~ All_Gender_clean, data=Length_Support)

#effect of CS on response
kruskal.test(Length_Support ~ Supporter_CareerStage_clean, data=Length_Support)

```


(This alternative is not reported)
Alternative recode

* 0- less than six months
  * 1: Less than six months
  * 4: I'm still supporting, <6 mo
* 1- 6 mo to 1 yr
  * 2: 6mo-1yr
* 2- >1 yr
  * 3: >1yr

ignore 5

```{r}
Length_Support=genderCS[,c("Lenght_Support", "Supporter_CareerStage_clean", "All_Gender_clean")]
colnames(Length_Support)[1]="Length_Support"
Length_Support=Length_Support[Length_Support$Length_Support %in% seq(1,4), ] #remove NAs and some choices
dim(Length_Support)[1] #this gives N
#recode
Length_Support[Length_Support$Length_Support==1, "Length_Support"]=0 
Length_Support[Length_Support$Length_Support==4, "Length_Support"]=0 
Length_Support[Length_Support$Length_Support==2, "Length_Support"]=1
Length_Support[Length_Support$Length_Support==3, "Length_Support"]=2

summary(glm(Length_Support~All_Gender_clean, data=Length_Support, family="gaussian"))
summary(glm(Length_Support~Supporter_CareerStage_clean, data=Length_Support, family="gaussian"))
summary(glm(Length_Support~Supporter_CareerStage_clean*All_Gender_clean, data=Length_Support, family="gaussian"))
```


## Q17c Impact_draining (Section 4, Panel 3 and Section 4, Panel 4)
Response: Impact_draining (1-5), PNTA(6) is removed

Model: Ordinal logistical regression

```{r imapact_draining}
#subset
ImpactDraining=genderCS[,c("Impact_Draining", "All_Gender_clean", "Supporter_CareerStage_clean")]
summary(ImpactDraining) #there are some NAs
#clean
ImpactDraining=ImpactDraining[(ImpactDraining$Impact_Draining %in% seq(1:5)),]
dim(ImpactDraining)[1] #this gives N
#factorise
ImpactDraining$Impact_Draining=as.factor(ImpactDraining$Impact_Draining)
```

### Impact_draining x gender stats (Section 4, panel 3)
```{r impact_draining_OR}
exp(coef(polr(Impact_Draining~All_Gender_clean, data=ImpactDraining, Hess=TRUE, method="logistic"))) #OR (men as baseline)
exp(confint(polr(Impact_Draining~All_Gender_clean, data=ImpactDraining, Hess=TRUE, method="logistic"))) #CI
kruskal.test(Impact_Draining ~ All_Gender_clean, data=ImpactDraining) #effect of gender on response
```

### Impact Draining x CS (Section 4, panel 4)

```{r}
exp(coef(polr(Impact_Draining~Supporter_CareerStage_clean, data=ImpactDraining, Hess=TRUE, method="logistic"))) # OR (PhD as baseline)
exp(confint(polr(Impact_Draining~Supporter_CareerStage_clean, data=ImpactDraining, Hess=TRUE, method="logistic"))) #CI
kruskal.test(Impact_Draining ~ Supporter_CareerStage_clean, data=ImpactDraining) #effect of CS on response

# multiway kruskal wallis
## CSpairs=combinations(n=5, r=2, v=levels(ImpactDraining$Supporter_CareerStage_clean), repeats.allowed = F)
# Compare between early PIs and other CSs only
## this requires the package "MHTdiscrete"
CSpairs=cbind(rep("Group leaders (<5yr)", 4), levels(ImpactDraining$Supporter_CareerStage_clean)[c(1,2,4,5)])
ImpactDrainingCSKW=data.frame(CS1=c(), CS2=c(), X2=c(), N=c(), df=c(), pVal=c(), adjpVal=c()) 
for (i in 1:dim(CSpairs)[1]) {
  comparedData=ImpactDraining[(ImpactDraining$Supporter_CareerStage_clean==CSpairs[i,1]|ImpactDraining$Supporter_CareerStage_clean==CSpairs[i,2]),]
  KW=kruskal.test(Impact_Draining ~ Supporter_CareerStage_clean, data=comparedData) #effect of CS on response
  ImpactDrainingCSKW=rbind(ImpactDrainingCSKW, data.frame(CS1=CSpairs[i,1], CS2=CSpairs[i,2], X2=as.numeric(KW$statistic), N=dim(comparedData)[1], df=as.numeric(KW$parameter), pVal=as.numeric(KW$p.value)))
}

ImpactDrainingCSKW$adjpVal=Sidak.p.adjust(ImpactDrainingCSKW$pVal)
ImpactDrainingCSKW
```

### Impact Draining x early PI vs others (PhD/postdoc/midPI/latePI)
```{r}
# code a new variable based on genderCS$Supporter_CareerStage_clean, 1= early PI, 0=PhD/postdoc/midPI/latePI
genderCS$earlyPI=rep("Group leaders (<5yr)", nrow(genderCS))
genderCS[(genderCS$Supporter_CareerStage_clean!="Group leaders (<5yr)"),"earlyPI"]="Other career stages"
genderCS$earlyPI=as.factor(genderCS$earlyPI)

#redo subsetting
ImpactDrainingPI=genderCS[,c("Impact_Draining", "All_Gender_clean", "earlyPI")]
summary(ImpactDrainingPI) #there are some NAs
#clean
ImpactDrainingPI=ImpactDrainingPI[(ImpactDrainingPI$Impact_Draining %in% seq(1:5)),]
dim(ImpactDrainingPI)[1] #this gives N
#factorise
ImpactDrainingPI$Impact_Draining=as.factor(ImpactDrainingPI$Impact_Draining)

#kruskal wallis
kruskal.test(Impact_Draining ~ earlyPI, data=ImpactDrainingPI)
```


### Impact_draining x gender graph (graph 4.3.1)
```{r}
levels(ImpactDraining$Impact_Draining)=c("Strongly disagree", "Somewhat disagree", "Neither agree or disagree", "Somewhat agree", "Strongly agree")
ImpactDraining$Impact_Draining=factor(ImpactDraining$Impact_Draining, levels=c("Strongly agree", "Somewhat agree", "Neither agree or disagree", "Somewhat disagree", "Strongly disagree"))
graphdata = ImpactDraining %>%
  group_by(All_Gender_clean, Impact_Draining) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))
#graphdata$perc=0
#graphdata[graphdata$All_Gender_clean=="Women",]=mutate(graphdata[graphdata$All_Gender_clean=="Women",], perc=n/sum(n))
#graphdata[graphdata$All_Gender_clean=="Men",]=mutate(graphdata[graphdata$All_Gender_clean=="Men",], perc=n/sum(n))

# This will save the image to your local code folder
#eps("images/ImpactDraining_gender.eps", width=1000, height=578)
ggplot(graphdata, aes(x=All_Gender_clean, y=perc, fill=Impact_Draining)) +
  geom_bar(stat="identity") +
  theme_minimal() +
      theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16), legend.position="bottom", legend.text = element_text(size=7))+
  geom_text(aes(label=round(perc, digit=1)), size=3, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="The experience was emotionally draining or stressful") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
ggsave(dpi=1000, "4-3-1.png", limitsize = FALSE)
#dev.off()
```



###  Impact_draining x career graph (graph 4.4.1)
```{r}
#levels(ImpactDraining$Impact_Draining)=c("Strongly disagree", "Somewhat disagree", "Neither agree or #disagree", "Somewhat agree", "Strongly agree")
#ImpactDraining$Impact_Draining=factor(ImpactDraining$Impact_Draining, levels=c("Strongly agree", "Somewhat agree", "Neither agree or disagree", "Somewhat disagree", "Strongly disagree"))
graphdata = ImpactDraining %>%
  group_by(Supporter_CareerStage_clean, Impact_Draining) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))
# This will save the image to your local code folder
#eps("images/ImpactDraining_CS.eps", width=1000, height=578)
ggplot(graphdata, aes(x=rev(Supporter_CareerStage_clean), y=rev(perc), fill=rev(Impact_Draining))) +
  geom_bar(stat="identity") +
  theme_minimal() +
      theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16), legend.position="bottom", legend.text = element_text(size=7))+
  geom_text(aes(label=round(rev(perc), digit=1)), size=3.5, position=position_stack(vjust=0.5), color="white") +
  scale_x_discrete(limits = rev(levels(ImpactDraining$Supporter_CareerStage_clean))) +
  labs(x="", y="Percentage", title="The experience was emotionally draining or stressful") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
ggsave(dpi=1000, "4-4-1.png", limitsize = FALSE)
#dev.off()
```

It looks like career stage 6 (early PIs) are the only group more likely to find it more draining than PhDs (career stage 3, also the intercept here).

We may be able to Chi-square here to show early career PIs are significantly different from other career stages? My concern is big differences in sample sizes

### Chi squares for interactions between CS and gender [NOT REPORTED]

```{r}
#split by career stages, gender effect
chisq.test(table(ImpactDraining[ImpactDraining$Supporter_CareerStage_clean==levels(ImpactDraining$Supporter_CareerStage_clean)[1],c("Impact_Draining","All_Gender_clean")]))
chisq.test(table(ImpactDraining[ImpactDraining$Supporter_CareerStage_clean==levels(ImpactDraining$Supporter_CareerStage_clean)[2],c("Impact_Draining","All_Gender_clean")]))
chisq.test(table(ImpactDraining[ImpactDraining$Supporter_CareerStage_clean==levels(ImpactDraining$Supporter_CareerStage_clean)[3],c("Impact_Draining","All_Gender_clean")]))
chisq.test(table(ImpactDraining[ImpactDraining$Supporter_CareerStage_clean==levels(ImpactDraining$Supporter_CareerStage_clean)[4],c("Impact_Draining","All_Gender_clean")]))
chisq.test(table(ImpactDraining[ImpactDraining$Supporter_CareerStage_clean==levels(ImpactDraining$Supporter_CareerStage_clean)[5],c("Impact_Draining","All_Gender_clean")]))
#split by gender, career effect (same as KW above)
chisq.test(table(ImpactDraining[ImpactDraining$All_Gender_clean==levels(ImpactDraining$All_Gender_clean)[1],c("Impact_Draining","Supporter_CareerStage_clean")]))
chisq.test(table(ImpactDraining[ImpactDraining$All_Gender_clean==levels(ImpactDraining$All_Gender_clean)[2],c("Impact_Draining","Supporter_CareerStage_clean")]))
```

## Q17d Impact Time (Seciont 4, Panel 3)

```{r}
#Subsetting
ImpactTime=genderCS[,c("Impact_Time", "All_Gender_clean", "Supporter_CareerStage_clean")]
#clean
ImpactTime=ImpactTime[(ImpactTime$Impact_Time %in% seq(1:5)),]
dim(ImpactTime)[1] #N
#factorise
ImpactTime$Impact_Time=as.factor(ImpactTime$Impact_Time)
table(ImpactTime$Impact_Time)
```

### Impact Time x gender (Section 4, Panel 3)
```{r}
exp(coef(polr(Impact_Time~All_Gender_clean, data=ImpactTime, Hess=TRUE, method="logistic"))) #OR (men as baseline)
exp(confint(polr(Impact_Time~All_Gender_clean, data=ImpactTime, Hess=TRUE, method="logistic"))) #CI
kruskal.test(Impact_Time ~ All_Gender_clean, data=ImpactTime) #effect of gender on response
```

### Impact Time x Career stage
```{r}
exp(coef(polr(Impact_Time~Supporter_CareerStage_clean, data=ImpactTime, Hess=TRUE, method="logistic"))) #OR (PhD as baseline)
kruskal.test(Impact_Time ~ Supporter_CareerStage_clean, data=ImpactTime) #effect of CS on response
```

### Impact Time x gender graph (graph 4.3.2)
```{r}
levels(ImpactTime$Impact_Time)=c("Strongly disagree", "Somewhat disagree", "Neither agree or disagree", "Somewhat agree", "Strongly agree")
ImpactTime$Impact_Time=factor(ImpactTime$Impact_Time, levels=c("Strongly agree", "Somewhat agree", "Neither agree or disagree", "Somewhat disagree", "Strongly disagree"))
graphdata = ImpactTime %>%
  group_by(All_Gender_clean, Impact_Time) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))
# This will save the image to your local code folder
#eps("images/ImpactTime_gender.eps", width=1000, height=578)
ggplot(graphdata, aes(x=All_Gender_clean, y=perc, fill=Impact_Time)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16), legend.position="bottom", legend.text = element_text(size=7))+
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="The experience took a lot of my time") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
ggsave(dpi=1000, "4-3-2.png", limitsize = FALSE)
dev.off()
```

### Impact Time x Career stage graph
```{r}
graphdata = ImpactTime %>%
  group_by(Supporter_CareerStage_clean, Impact_Time) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))
# This will save the image to your local code folder
#eps("images/ImpactTime_CS.eps", width=1000, height=578)
ggplot(graphdata, aes(x=Supporter_CareerStage_clean, y=perc, fill=Impact_Time)) +
  geom_bar(stat="identity") +
  theme_minimal() +
      theme(plot.title = element_text(hjust = 0.5), text=element_text(size=20), legend.position="bottom", legend.text = element_text(size=7))+
  scale_x_discrete(limits = rev(levels(ImpactTime$Supporter_CareerStage_clean))) +
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="The experience took a lot of my time") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
dev.off()
```


## Q17e Impact work (Section 4, Panel 3 and Section 4, Panel 4)
```{r}
#Subset
ImpactWork=genderCS[,c("Impact_Work", "All_Gender_clean", "Supporter_CareerStage_clean")]
#clean
ImpactWork=ImpactWork[(ImpactWork$Impact_Work %in% seq(1:5)),]
dim(ImpactWork)[1]
#factorise
ImpactWork$Impact_Work=as.factor(ImpactWork$Impact_Work)
table(ImpactWork$Impact_Work)
``` 

### Impact work x gender (Section 4, Panel 3)
```{r}
exp(coef(polr(Impact_Work~All_Gender_clean, data=ImpactWork, Hess=TRUE, method="logistic"))) #OR (men as baseline)
exp(confint(polr(Impact_Work~All_Gender_clean, data=ImpactWork, Hess=TRUE, method="logistic"))) #CI
kruskal.test(Impact_Work ~ All_Gender_clean, data=ImpactWork) #effect of gender on response
```

### Impact work x Career stage (Section 4, Panel 4)
```{r}
exp(coef(polr(Impact_Work~Supporter_CareerStage_clean, data=ImpactWork, Hess=TRUE, method="logistic"))) #OR (PhD as baseline)
exp(confint(polr(Impact_Work~Supporter_CareerStage_clean, data=ImpactWork, Hess=TRUE, method="logistic"))) #CI
kruskal.test(Impact_Work ~ Supporter_CareerStage_clean, data=ImpactWork) #effect of CS on response

## multiway kruskal wallis
# Compare between early PIs and other CSs onl
## this requires the package "MTHdiscrete", and having run the similar analysis in the Impact draining section
ImpactWorkCSKW=data.frame(CS1=c(), CS2=c(), X2=c(), N=c(), df=c(), pVal=c()) 
for (i in 1:dim(CSpairs)[1]) {
  comparedData=ImpactWork[(ImpactWork$Supporter_CareerStage_clean==CSpairs[i,1]|ImpactWork$Supporter_CareerStage_clean==CSpairs[i,2]),]
  KW=kruskal.test(Impact_Work ~ Supporter_CareerStage_clean, data=comparedData) #effect of CS on response
  ImpactWorkCSKW=rbind(ImpactWorkCSKW, data.frame(CS1=CSpairs[i,1], CS2=CSpairs[i,2], X2=as.numeric(KW$statistic), N=dim(comparedData)[1], df=as.numeric(KW$parameter), pVal=as.numeric(KW$p.value)))
}
ImpactWorkCSKW$adjpVal=Sidak.p.adjust(ImpactWorkCSKW$pVal)
ImpactWorkCSKW
```

### Impact work x PI vs others
```{r}
#redo subsetting
ImpactWorkPI=genderCS[,c("Impact_Work", "All_Gender_clean", "earlyPI")]
summary(ImpactWorkPI) #there are some NAs
#clean
ImpactWorkPI=ImpactWorkPI[(ImpactWorkPI$Impact_Work %in% seq(1:5)),]
dim(ImpactWorkPI)[1] #this gives N
#factorise
ImpactWorkPI$Impact_Work=as.factor(ImpactWorkPI$Impact_Work)

#kruskal wallis
kruskal.test(Impact_Work ~ earlyPI, data=ImpactWorkPI)
```


### Impact work x gender graph (Graph 4.3.4)
```{r}
levels(ImpactWork$Impact_Work)=c("Strongly disagree", "Somewhat disagree", "Neither agree or disagree", "Somewhat agree", "Strongly agree")
ImpactWork$Impact_Work=factor(ImpactWork$Impact_Work, levels=c("Strongly agree", "Somewhat agree", "Neither agree or disagree", "Somewhat disagree", "Strongly disagree"))
graphdata = ImpactWork %>%
  group_by(All_Gender_clean, Impact_Work) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))
# This will save the image to your local code folder
#eps("images/ImpactWork_gender.eps", width=1000, height=578)
ggplot(graphdata, aes(x=All_Gender_clean, y=perc, fill=Impact_Work)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16), legend.position="bottom", legend.text = element_text(size=7))+
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="The experience negatively impacted my academic work") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
ggsave(dpi=1000, "4-3-4 GOOD.png", limitsize = FALSE)
dev.off()
```

### Impact work x Career stage graph (Graph 4.4.2)
```{r}
graphdata = ImpactWork %>%
  group_by(Supporter_CareerStage_clean, Impact_Work) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/ImpactWork_CS.eps", width=1000, height=578)
ggplot(graphdata, aes(x=Supporter_CareerStage_clean, y=perc, fill=rev(Impact_Work))) +
  geom_bar(stat="identity") +
  theme_minimal() +
      theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16), legend.position="bottom", legend.text = element_text(size=7))+
  scale_x_discrete(limits = rev(levels(ImpactWork$Supporter_CareerStage_clean))) +
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="The experience negatively impacted my academic work") +
    guides(fill=guide_legend(reverse=FALSE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
ggsave(dpi=1000, "4-4-2.png", limitsize = FALSE)
dev.off()
```
### chisq for intersections [NOT REPORTED]

```{r}
#split by career stages, gender effect
chisq.test(table(ImpactWork[ImpactWork$Supporter_CareerStage_clean==levels(ImpactWork$Supporter_CareerStage_clean)[1],c("Impact_Work","All_Gender_clean")]))
chisq.test(table(ImpactWork[ImpactWork$Supporter_CareerStage_clean==levels(ImpactWork$Supporter_CareerStage_clean)[2],c("Impact_Work","All_Gender_clean")]))
chisq.test(table(ImpactWork[ImpactWork$Supporter_CareerStage_clean==levels(ImpactWork$Supporter_CareerStage_clean)[3],c("Impact_Work","All_Gender_clean")]))
chisq.test(table(ImpactWork[ImpactWork$Supporter_CareerStage_clean==levels(ImpactWork$Supporter_CareerStage_clean)[4],c("Impact_Work","All_Gender_clean")]))
chisq.test(table(ImpactWork[ImpactWork$Supporter_CareerStage_clean==levels(ImpactWork$Supporter_CareerStage_clean)[5],c("Impact_Work","All_Gender_clean")]))
#split by gender, career effect (same as KW above)
chisq.test(table(ImpactWork[ImpactWork$All_Gender_clean==levels(ImpactWork$All_Gender_clean)[1],c("Impact_Work","Supporter_CareerStage_clean")]))
chisq.test(table(ImpactWork[ImpactWork$All_Gender_clean==levels(ImpactWork$All_Gender_clean)[2],c("Impact_Work","Supporter_CareerStage_clean")]))
```

## Q17f Impact personal (Section 4, Panel 3 and Section 4, Panel 4)
```{r}
ImpactPersonal=genderCS[,c("Impact_Personal", "All_Gender_clean", "Supporter_CareerStage_clean")]
#clean
ImpactPersonal=ImpactPersonal[(ImpactPersonal$Impact_Personal %in% seq(1:5)),]
dim(ImpactPersonal)[1]
#factorise
ImpactPersonal$Impact_Personal=as.factor(ImpactPersonal$Impact_Personal)
table(ImpactPersonal$Impact_Personal)
``` 

###  Impact personal x gender stats (Section 4, Panel 3)
```{r}
exp(coef(polr(Impact_Personal~All_Gender_clean, data=ImpactPersonal, Hess=TRUE, method="logistic"))) #OR (men as baseline)
exp(confint(polr(Impact_Personal~All_Gender_clean, data=ImpactPersonal, Hess=TRUE, method="logistic"))) #CI
kruskal.test(Impact_Personal ~ All_Gender_clean, data=ImpactPersonal) #effect of gender on response
```

### Impact personal x Career stage (Section 4, Panel 4)
```{r}
exp(coef(polr(Impact_Personal~Supporter_CareerStage_clean, data=ImpactPersonal, Hess=TRUE, method="logistic"))) #OR (PhD as baseline)
exp(confint(polr(Impact_Personal~Supporter_CareerStage_clean, data=ImpactPersonal, Hess=TRUE, method="logistic"))) #CI
kruskal.test(Impact_Personal ~ Supporter_CareerStage_clean, data=ImpactPersonal) #effect of CS on response
```

###  Impact personal x gender graph (Graph 4.3.3)
````{r}
levels(ImpactPersonal$Impact_Personal)=c("Strongly disagree", "Somewhat disagree", "Neither agree or disagree", "Somewhat agree", "Strongly agree")
ImpactPersonal$Impact_Personal=factor(ImpactPersonal$Impact_Personal, levels=c("Strongly agree", "Somewhat agree", "Neither agree or disagree", "Somewhat disagree", "Strongly disagree"))
# This will save the image to your local code folder
#eps("images/ImpactPersonal_gender.eps", width=1000, height=578)
graphdata = ImpactPersonal %>%
  group_by(All_Gender_clean, Impact_Personal) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))
ggplot(graphdata, aes(x=All_Gender_clean, y=perc, fill=Impact_Personal)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16), legend.position="bottom", legend.text = element_text(size=7))+
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="The experience had an impact on my personal life") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
ggsave(dpi=1000, "4-3-3.png", limitsize = FALSE)
dev.off()
```

### Impact personal x Career stage (Graph 4.4.3)
```{r}
#split by gender
kruskal.test(Impact_Personal ~ Supporter_CareerStage_clean, data=ImpactPersonal[ImpactPersonal$All_Gender_clean==levels(ImpactPersonal$All_Gender_clean)[1],])
kruskal.test(Impact_Personal ~ Supporter_CareerStage_clean, data=ImpactPersonal[ImpactPersonal$All_Gender_clean==levels(ImpactPersonal$All_Gender_clean)[2],])

#split by career stages, gender effect
chisq.test(table(ImpactPersonal[ImpactPersonal$Supporter_CareerStage_clean==levels(ImpactPersonal$Supporter_CareerStage_clean)[1],c("Impact_Personal","All_Gender_clean")]))
chisq.test(table(ImpactPersonal[ImpactPersonal$Supporter_CareerStage_clean==levels(ImpactPersonal$Supporter_CareerStage_clean)[2],c("Impact_Personal","All_Gender_clean")]))
chisq.test(table(ImpactPersonal[ImpactPersonal$Supporter_CareerStage_clean==levels(ImpactPersonal$Supporter_CareerStage_clean)[3],c("Impact_Personal","All_Gender_clean")]))
chisq.test(table(ImpactPersonal[ImpactPersonal$Supporter_CareerStage_clean==levels(ImpactPersonal$Supporter_CareerStage_clean)[4],c("Impact_Personal","All_Gender_clean")]))
chisq.test(table(ImpactPersonal[ImpactPersonal$Supporter_CareerStage_clean==levels(ImpactPersonal$Supporter_CareerStage_clean)[5],c("Impact_Personal","All_Gender_clean")]))
#split by gender, career effect (same as KW above)
chisq.test(table(ImpactPersonal[ImpactPersonal$All_Gender_clean==levels(ImpactPersonal$All_Gender_clean)[1],c("Impact_Personal","Supporter_CareerStage_clean")]))
chisq.test(table(ImpactPersonal[ImpactPersonal$All_Gender_clean==levels(ImpactPersonal$All_Gender_clean)[2],c("Impact_Personal","Supporter_CareerStage_clean")]))

graphdata = ImpactPersonal %>%
  group_by(Supporter_CareerStage_clean, Impact_Personal) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/ImpactPersonal_CS.eps", width=1000, height=578)
ggplot(graphdata, aes(x=Supporter_CareerStage_clean, y=perc, fill=Impact_Personal)) +
  geom_bar(stat="identity") +
  theme_minimal() +
      theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16), legend.position="bottom", legend.text = element_text(size=7))+
  scale_x_discrete(limits = rev(levels(ImpactPersonal$Supporter_CareerStage_clean))) +
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="The experience had an impact on my personal life") +
    guides(fill=guide_legend(reverse=TRUE)) +
 coord_flip() +
  scale_fill_manual(name="", values=ePalette)
ggsave(dpi=1000, "4-4-3-TRY.png", limitsize = FALSE)
#dev.off()
```


## Q12 Professional support for the receiver (Section 3, Panel 2)

Dependent variable: Receiver_ProfessionalHelp

Original question: During the time I was supporting them, the person received at least some professional help
* 1: No
* 2: Yes, it started before I started supporting them
* 3: Yes, it started after I started supporting them
* 4: Yes, but I'm unsure when it started
(remove 5- not sure, and 6 - PNTA)

We would like to check if respondents supported someone who doesn't have professional support at the time the support was offered.

hence recode to:
* 0: No(1), and (3) Yes it started after
* 1: Yes, it started before (2)
removing 2- unsure.

Model: First individual factors (Gender, CS), then Factorial logistical regression Supporter_NumberReceivers ~ Supporter_CareerStage * All_Gender_clean) to check for interactions

```{r Receiver_ProfessionalHelp}

# subsetting data and recoding response to binary
Receiver_ProfessionalHelp=genderCS[,c("Receiver_ProfessionalHelp", "Supporter_CareerStage_clean", "All_Gender_clean")]
#recode andclean
Receiver_ProfessionalHelp[(Receiver_ProfessionalHelp$Receiver_ProfessionalHelp %in% c(1,3)),"Receiver_ProfessionalHelp"]=0
Receiver_ProfessionalHelp[(Receiver_ProfessionalHelp$Receiver_ProfessionalHelp %in% c(2)),"Receiver_ProfessionalHelp"]=1
Receiver_ProfessionalHelp=Receiver_ProfessionalHelp[Receiver_ProfessionalHelp$Receiver_ProfessionalHelp %in% seq(0,1), ] #remove NAs and some choices
dim(Receiver_ProfessionalHelp)[1]

#LR gender
summary(glm(Receiver_ProfessionalHelp~All_Gender_clean, data=Receiver_ProfessionalHelp, family="binomial")) #glm and p-value

#LR CS
summary(glm(Receiver_ProfessionalHelp~Supporter_CareerStage_clean, data=Receiver_ProfessionalHelp, family="binomial")) #glm and p-value
exp(coef(glm(Receiver_ProfessionalHelp~Supporter_CareerStage_clean, data=Receiver_ProfessionalHelp, family="binomial"))) #OR (PhD as baseline)
exp(confint(glm(Receiver_ProfessionalHelp~Supporter_CareerStage_clean, data=Receiver_ProfessionalHelp, family="binomial"))) #CI

#FLR gender-CS interactions
summary(glm(Receiver_ProfessionalHelp~Supporter_CareerStage_clean*All_Gender_clean, data=Receiver_ProfessionalHelp, family="binomial"))
```

### Professional support for the receiver x (PhD+ Postdoc) vs PIs
Code new variable PI:
- 0: PhD+ postdoc
- 1: PIs of all stages
```{r}
genderCS$PI="Group leaders"
genderCS[(genderCS$Supporter_CareerStage_clean==levels(genderCS$Supporter_CareerStage_clean)[1] | genderCS$Supporter_CareerStage_clean==levels(genderCS$Supporter_CareerStage_clean)[2]),'PI']="PhDs and Postdocs"
genderCS$PI=as.factor(genderCS$PI)
genderCS$PI=relevel(genderCS$PI, "PhDs and Postdocs") #relevel so that "PhDs and postdocs"" is baseline
summary(genderCS$PI)
# redo subsetting
# subsetting data and recoding response to binary
Receiver_ProfessionalHelp=genderCS[,c("Receiver_ProfessionalHelp", "Supporter_CareerStage_clean", "PI", "All_Gender_clean")]
#recode andclean
Receiver_ProfessionalHelp[(Receiver_ProfessionalHelp$Receiver_ProfessionalHelp %in% c(1,3)),"Receiver_ProfessionalHelp"]=0
Receiver_ProfessionalHelp[(Receiver_ProfessionalHelp$Receiver_ProfessionalHelp %in% c(2)),"Receiver_ProfessionalHelp"]=1
Receiver_ProfessionalHelp=Receiver_ProfessionalHelp[Receiver_ProfessionalHelp$Receiver_ProfessionalHelp %in% seq(0,1), ] #remove NAs and some choices
dim(Receiver_ProfessionalHelp)[1]
#LR
summary(glm(Receiver_ProfessionalHelp~PI, data=Receiver_ProfessionalHelp, family="binomial")) #glm and p-value
exp(coef(glm(Receiver_ProfessionalHelp~PI, data=Receiver_ProfessionalHelp, family="binomial"))) #OR (PhD as baseline)
exp(confint(glm(Receiver_ProfessionalHelp~PI, data=Receiver_ProfessionalHelp, family="binomial"))) #CI
### Professional support for the receiver x gender graph [NOT REPORTED]

```{r}
Receiver_ProfessionalHelp[Receiver_ProfessionalHelp$Receiver_ProfessionalHelp==0, "Receiver_ProfessionalHelp"]="No"
Receiver_ProfessionalHelp[Receiver_ProfessionalHelp$Receiver_ProfessionalHelp==1, "Receiver_ProfessionalHelp"]="Yes"
Receiver_ProfessionalHelp$Receiver_ProfessionalHelp=factor(Receiver_ProfessionalHelp$Receiver_ProfessionalHelp)
#ImpactPersonal$Impact_Personal=factor(ImpactPersonal$Impact_Personal, levels=c("No","Yes"))
graphdata =Receiver_ProfessionalHelp %>%
  group_by(All_Gender_clean, Receiver_ProfessionalHelp) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/ReceiverProfessionalHelp_gender.eps", width=1000, height=578)
ggplot(graphdata, aes(x=All_Gender_clean, y=perc, fill=Receiver_ProfessionalHelp)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=20))+
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="During the time I was supporting them, the person received at least some professional help") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=c(eBlue, eGreen))
dev.off()
```

### Professional support for the receiver career graph 
```{r}
graphdata = Receiver_ProfessionalHelp %>%
  group_by(Supporter_CareerStage_clean, Receiver_ProfessionalHelp) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))


### Professional support for the receiver, PhDs& postdocs vs PIs graph (New graph 3.2.1?) 
```{r}
graphdata = Receiver_ProfessionalHelp %>%
  group_by(PI, Receiver_ProfessionalHelp) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))
# This will save the image to your local code folder
ggplot(graphdata, aes(x=PI, y=perc, fill=Receiver_ProfessionalHelp)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16), legend.text = element_text(size=7))+ 
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="During the time I was supporting them, the person received at least some professional help") +
    guides(fill=guide_legend(reverse=TRUE)) +
  scale_x_discrete(limits = rev(levels(graphdata$PI))) +
  coord_flip() +
  scale_fill_manual(name="", values=c(eBlue, eGreen))
ggsave(dpi=1000, "3-2-1.png", limitsize = FALSE)
dev.off()

```

## Q17a Comfortable_withMH (Seciont 6, Panel 1)
Model: ordinal logsitic regression
```{r}
#Subset
ComfortablewithMH=genderCS[,c("Comfortable_withMH", "All_Gender_clean", "Supporter_CareerStage_clean")]
#clean
ComfortablewithMH=ComfortablewithMH[(ComfortablewithMH$Comfortable_withMH %in% seq(1:5)),]
dim(ComfortablewithMH)[1]
#factorise
ComfortablewithMH$Comfortable_withMH=as.factor(ComfortablewithMH$Comfortable_withMH)
table(ComfortablewithMH$Comfortable_withMH)
```

### Comfortable_withMH x gender
```{r}
exp(coef(polr(Comfortable_withMH~All_Gender_clean, data=ComfortablewithMH, Hess=TRUE, method="logistic"))) #OR (men as baseline)
exp(confint(polr(Comfortable_withMH~All_Gender_clean, data=ComfortablewithMH, Hess=TRUE, method="logistic"))) #CI
kruskal.test(Comfortable_withMH ~ All_Gender_clean, data=ComfortablewithMH) #effect of gender on response
```

### Comfortable_withMH x CS
```{r}
exp(coef(polr(Comfortable_withMH~Supporter_CareerStage_clean, data=ComfortablewithMH, Hess=TRUE, method="logistic"))) #OR (PhD as baseline)
exp(confint(polr(Comfortable_withMH~Supporter_CareerStage_clean, data=ComfortablewithMH, Hess=TRUE, method="logistic"))) #CI
kruskal.test(Comfortable_withMH ~ Supporter_CareerStage_clean, data=ComfortablewithMH) #effect of CS on baseline
```

### Comfortable_withMH x gender graph (Graph 6.1.3)
```{r}
levels(ComfortablewithMH$Comfortable_withMH)=c("Strongly disagree", "Somewhat disagree", "Neither agree or disagree", "Somewhat agree", "Strongly agree")
ComfortablewithMH$Comfortable_withMH=factor(ComfortablewithMH$Comfortable_withMH, levels=c("Strongly agree", "Somewhat agree", "Neither agree or disagree", "Somewhat disagree", "Strongly disagree"))

graphdata = ComfortablewithMH %>%
  group_by(All_Gender_clean, Comfortable_withMH) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/ComfortablewithMH_gender.eps", width=1000, height=578)
ggplot(graphdata, aes(x=rev(All_Gender_clean), y=rev(perc), fill=rev(Comfortable_withMH))) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16), legend.position="bottom", legend.text = element_text(size=7))+
  geom_text(aes(label=round(rev(perc), digit=1)), size=3.5, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="I was comfortable discussing mental health problems") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
ggsave(dpi=1000, "6-1-3.png", limitsize = FALSE)
dev.off()
```

### Comfortable_withMH x CS graph
```{r}
graphdata = ComfortablewithMH %>%
  group_by(Supporter_CareerStage_clean, Comfortable_withMH) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/ComfortablewithMH_CS.eps", width=1000, height=578)
ggplot(graphdata, aes(x=Supporter_CareerStage_clean, y=perc, fill=Comfortable_withMH)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=20))+
  scale_x_discrete(limits = rev(levels(ComfortablewithMH$Supporter_CareerStage_clean))) +
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="I was comfortable discussing mental health problems") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
dev.off()
```

## Q17b Confident_withMH (Section 6, Panel 1)
Model: ordinal logsitic regression
```{r}
#Subset
ConfidentwithMH=genderCS[,c("Confident_withMH", "All_Gender_clean", "Supporter_CareerStage_clean")]
#clean
ConfidentwithMH=ConfidentwithMH[(ConfidentwithMH$Confident_withMH %in% seq(1:5)),]
dim(ConfidentwithMH)[1]
#factorise
ConfidentwithMH$Confident_withMH=as.factor(ConfidentwithMH$Confident_withMH)
table(ConfidentwithMH$Confident_withMH)
```

### Confident_withMH x gender
```{r}
exp(coef(polr(Confident_withMH~All_Gender_clean, data=ConfidentwithMH, Hess=TRUE, method="logistic"))) #OR (men as baseline)
exp(confint(polr(Confident_withMH~All_Gender_clean, data=ConfidentwithMH, Hess=TRUE, method="logistic"))) #CI
kruskal.test(Confident_withMH ~ All_Gender_clean, data=ConfidentwithMH) #effect of gender on response
```

### Confident_withMH x CS
```{r}
exp(coef(polr(Confident_withMH~Supporter_CareerStage_clean, data=ConfidentwithMH, Hess=TRUE, method="logistic"))) #OR (PhD as baseline)
exp(confint(polr(Confident_withMH~Supporter_CareerStage_clean, data=ConfidentwithMH, Hess=TRUE, method="logistic"))) #CI
kruskal.test(Confident_withMH ~ Supporter_CareerStage_clean, data=ConfidentwithMH) #effect of CS on response

## multiway kruskal wallis
## this requries the package 'gtools'
## CSpairs=combinations(n=5, r=2, v=levels(ConfidentwithMH$Supporter_CareerStage_clean), repeats.allowed = F)
# Compare between early PIs and other CSs only
## this requires the package "MTHdiscrete"
CSpairs=cbind(rep("Group leaders (<5yr)", 4), levels(ImpactDraining$Supporter_CareerStage_clean)[c(1,2,4,5)])
ConfidentwithMHCSKW=data.frame(CS1=c(), CS2=c(), X2=c(), N=c(), df=c(), pVal=c()) 
for (i in 1:dim(CSpairs)[1]) {
  comparedData=ConfidentwithMH[(ConfidentwithMH$Supporter_CareerStage_clean==CSpairs[i,1]|ConfidentwithMH$Supporter_CareerStage_clean==CSpairs[i,2]),]
  KW=kruskal.test(Confident_withMH ~ Supporter_CareerStage_clean, data=comparedData) #effect of CS on response
  ConfidentwithMHCSKW=rbind(ConfidentwithMHCSKW, data.frame(CS1=CSpairs[i,1], CS2=CSpairs[i,2], X2=as.numeric(KW$statistic), N=dim(comparedData)[1], df=as.numeric(KW$parameter), pVal=as.numeric(KW$p.value)))
}

ConfidentwithMHCSKW$adjpVal=Sidak.p.adjust(ConfidentwithMHCSKW$pVal)
ConfidentwithMHCSKW
```

### Confident_with MH x early PI vs others
Make sure you run the corresponding "impact draining" chunk (under heading "### Impact Draining x early PI vs others (PhD/postdoc/midPI/latePI)") before running this (need to have early PI variable in genderCS)
```{r}
#redo subsetting
ConfidentwithMHPI=genderCS[,c("Confident_withMH", "All_Gender_clean", "earlyPI")]
summary(ConfidentwithMHPI) #there are some NAs
#clean
ConfidentwithMHPI=ConfidentwithMHPI[(ConfidentwithMHPI$Confident_withMH %in% seq(1:5)),]
dim(ConfidentwithMHPI)[1] #this gives N
#factorise
ConfidentwithMHPI$Confident_withMH=as.factor(ConfidentwithMHPI$Confident_withMH)

#kruskal wallis
kruskal.test(Confident_withMH ~ earlyPI, data=ConfidentwithMHPI)
```



### Confident_withMH x gender graph
```{r}
levels(ConfidentwithMH$Confident_withMH)=c("Strongly disagree", "Somewhat disagree", "Neither agree or disagree", "Somewhat agree", "Strongly agree")
ConfidentwithMH$Confident_withMH=factor(ConfidentwithMH$Confident_withMH, levels=c("Strongly agree", "Somewhat agree", "Neither agree or disagree", "Somewhat disagree", "Strongly disagree"))
graphdata = ConfidentwithMH %>%
  group_by(All_Gender_clean, Confident_withMH) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/ConfidentwithMH_gender.eps", width=1000, height=578)
ggplot(graphdata, aes(x=All_Gender_clean, y=perc, fill=Confident_withMH)) +
  geom_bar(stat="identity") +
  theme_minimal() +
     theme(plot.title = element_text(hjust = 0.5), text=element_text(size=20))+
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="I was confident in my ability to do or say the right thing") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
dev.off()
```

### Confident_withMH x CS graph (Graph 6.1.4)
```{r}
graphdata = ConfidentwithMH %>%
  group_by(Supporter_CareerStage_clean, Confident_withMH) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/ConfidentwithMH_CS.eps", width=1000, height=578)
ggplot(graphdata, aes(x=Supporter_CareerStage_clean, y=perc, fill=Confident_withMH)) +
  geom_bar(stat="identity") +
  theme_minimal() +
      theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16), legend.position="bottom", legend.text = element_text(size=7))+
  scale_x_discrete(limits = rev(levels(ConfidentwithMH$Supporter_CareerStage_clean))) +
  geom_text(aes(label=round(perc, digit=1)), size=3, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="I was confident in my ability to do or say the right thing") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
ggsave(dpi=1000, "6-1-4 GOOD.png", limitsize = FALSE)
dev.off()
```

## Q18a Impact Rewarding (Seciont 4, Panel 1)
Model: ordinal logistic regression

```{r Impact_rewarding}
#subsetting
ImpactRewarding=genderCS[,c("Impact_Rewarding", "All_Gender_clean", "Supporter_CareerStage_clean")]
#clean
ImpactRewarding=ImpactRewarding[(ImpactRewarding$Impact_Rewarding %in% seq(1:5)),]
dim(ImpactRewarding)[1]
#factorise
ImpactRewarding$Impact_Rewarding=as.factor(ImpactRewarding$Impact_Rewarding)
table(ImpactRewarding$Impact_Rewarding)
```

### gender
```{r}
exp(coef(polr(Impact_Rewarding~All_Gender_clean, data=ImpactRewarding, Hess=TRUE, method="logistic"))) #OR (men as baseline)
kruskal.test(Impact_Rewarding ~ All_Gender_clean, data=ImpactRewarding) #effect of gender on response
```

### CS
```{r}
exp(coef(polr(Impact_Rewarding~Supporter_CareerStage_clean, data=ImpactRewarding, Hess=TRUE, method="logistic"))) #OR (men as baseline)
kruskal.test(Impact_Rewarding ~ Supporter_CareerStage_clean, data=ImpactRewarding) #effect of CS on response
```


## Q18b Impact Positive Difference
Model: ordinal logistic regression
```{r impact_positivedifference}
#Subsetting
ImpactPositiveDifference=genderCS[,c("Impact_PositiveDifference", "All_Gender_clean", "Supporter_CareerStage_clean")]
#clean
ImpactPositiveDifference=ImpactPositiveDifference[(ImpactPositiveDifference$Impact_PositiveDifference %in% seq(1:5)),]
dim(ImpactPositiveDifference)[1]
#factorise
ImpactPositiveDifference$Impact_PositiveDifference=as.factor(ImpactPositiveDifference$Impact_PositiveDifference)
table(ImpactPositiveDifference$Impact_PositiveDifference)
```

### Impact Positive Difference x gender
```{r}
exp(coef(polr(Impact_PositiveDifference~All_Gender_clean, data=ImpactPositiveDifference, Hess=TRUE, method="logistic"))) #OR (men as baseline)
kruskal.test(Impact_PositiveDifference ~ All_Gender_clean, data=ImpactPositiveDifference) #effect of gender on response
```

### Impact Positive Difference x CS
```{r}
exp(coef(polr(Impact_PositiveDifference~Supporter_CareerStage_clean, data=ImpactPositiveDifference, Hess=TRUE, method="logistic"))) #OR PhD as baseline
kruskal.test(Impact_PositiveDifference ~ Supporter_CareerStage_clean, data=ImpactPositiveDifference) #effect of CS on response
```

## Q20d Support Receiver (Section 4, Panel 1)
Model: ordinal logistic regression
```{r support_receiver}
#subset
SupportReceiver=genderCS[,c("Support_Receiver", "All_Gender_clean", "Supporter_CareerStage_clean")]
#clean
SupportReceiver=SupportReceiver[(SupportReceiver$Support_Receiver %in% seq(1:5)),]
dim(SupportReceiver)[1]
#factorise
SupportReceiver$Support_Receiver=as.factor(SupportReceiver$Support_Receiver)
table(SupportReceiver$Support_Receiver)
```

### gender
```{r}
exp(coef(polr(Support_Receiver~All_Gender_clean, data=SupportReceiver, Hess=TRUE, method="logistic"))) #OR (men as baseline)
kruskal.test(Support_Receiver ~ All_Gender_clean, data=SupportReceiver) # effect of gender on response
```

### CS
```{r}
exp(coef(polr(Support_Receiver~Supporter_CareerStage_clean, data=SupportReceiver, Hess=TRUE, method="logistic"))) #OR (PhD as baseline)
kruskal.test(Support_Receiver ~ Supporter_CareerStage_clean, data=SupportReceiver) #effect of CS on response

```

## Q18c Impact Beacon (Section 4, Panel 1)
Model: ordinal logistic regression
```{r}
#subset
ImpactBeacon=genderCS[,c("Impact_Beacon", "All_Gender_clean", "Supporter_CareerStage_clean")]
#clean
ImpactBeacon=ImpactBeacon[(ImpactBeacon$Impact_Beacon %in% seq(1:5)),]
dim(ImpactBeacon)[1]
#factorise
ImpactBeacon$Impact_Beacon=as.factor(ImpactBeacon$Impact_Beacon)
table(ImpactBeacon$Impact_Beacon)
```

### Impact Beacon x gender
```{r}
exp(coef(polr(Impact_Beacon~All_Gender_clean, data=ImpactBeacon, Hess=TRUE, method="logistic"))) #OR (men as baseline)
exp(confint(polr(Impact_Beacon~All_Gender_clean, data=ImpactBeacon, Hess=TRUE, method="logistic"))) #CI

kruskal.test(Impact_Beacon ~ All_Gender_clean, data=ImpactBeacon) #effect of gender on response
```

### Impact Beacon x CS
```{r}
exp(coef(polr(Impact_Beacon~Supporter_CareerStage_clean, data=ImpactBeacon, Hess=TRUE, method="logistic"))) #OR (PhD as baseline)
kruskal.test(Impact_Beacon ~ Supporter_CareerStage_clean, data=ImpactBeacon) #effect of CS on response
```

## Q20a Support Institutions (Section 5, Panel 2)
Model: ordinal logistic regression
```{r}
SupportInstitution=genderCS[,c("Support_Institution", "All_Gender_clean", "Supporter_CareerStage_clean")]
#clean
SupportInstitution=SupportInstitution[(SupportInstitution$Support_Institution %in% seq(1:5)),]
dim(SupportInstitution)[1]
#factorise
SupportInstitution$Support_Institution=as.factor(SupportInstitution$Support_Institution)
table(SupportInstitution$Support_Institution)
```

### Support Institutions x gender
```{r}
exp(coef(polr(Support_Institution~All_Gender_clean, data=SupportInstitution, Hess=TRUE, method="logistic"))) #OR (men as baseline)
exp(confint(polr(Support_Institution~All_Gender_clean, data=SupportInstitution, Hess=TRUE, method="logistic"))) #CI
kruskal.test(Support_Institution ~ All_Gender_clean, data=SupportInstitution) #effect of gender on response
```

### Support Institutions x CS
```{r}
exp(coef(polr(Support_Institution~Supporter_CareerStage_clean, data=SupportInstitution, Hess=TRUE, method="logistic"))) # OR (PhD as baseline)
exp(confint(polr(Support_Institution~Supporter_CareerStage_clean, data=SupportInstitution, Hess=TRUE, method="logistic"))) #CI 
kruskal.test(Support_Institution ~ Supporter_CareerStage_clean, data=SupportInstitution) #effect on CS on response
```

### Support Institution x gender graph (Graph 5.2.4)
```{r}
levels(SupportInstitution$Support_Institution)=c("Strongly disagree", "Somewhat disagree", "Neither agree or disagree", "Somewhat agree", "Strongly agree")
SupportInstitution$Support_Institution=factor(SupportInstitution$Support_Institution, levels=c("Strongly agree", "Somewhat agree", "Neither agree or disagree", "Somewhat disagree", "Strongly disagree"))
graphdata = SupportInstitution %>%
  group_by(All_Gender_clean, Support_Institution) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/SupportInstitution_gender.eps", width=1000, height=578)
ggplot(graphdata, aes(x=All_Gender_clean, y=perc, fill=Support_Institution)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16), legend.position="bottom", legend.text = element_text(size=7))+
  geom_text(aes(label=round(perc, digit=1)), size=3, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="I feel supported or valued by my institution \n(e.g. managers, department) for the help I was providing") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
ggsave(dpi=1000, "5-2-4.png", limitsize = FALSE)
dev.off()
```

### Support instituions x CS graph (Graph 5.2.5)
```{r}
graphdata = SupportInstitution %>%
  group_by(Supporter_CareerStage_clean, Support_Institution) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/SupportInstitution_CS.eps", width=1000, height=578)
ggplot(graphdata, aes(x=Supporter_CareerStage_clean, y=perc, fill=Support_Institution)) +
  geom_bar(stat="identity") +
  theme_minimal() +
      theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16), legend.position="bottom", legend.text = element_text(size=7))+
  scale_x_discrete(limits = rev(levels(SupportInstitution$Supporter_CareerStage_clean))) +
  geom_text(aes(label=round(perc, digit=1)), size=3, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="I feel supported or valued by my institution (e.g. managers, department) for the help I was providing") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
ggsave(dpi=1000, "5-2-5.png", limitsize = FALSE)
dev.off()
```

## Q20b Support colleagues (Section 5, Panel 2)
Model: ordinal logistic regression
```{r}
#subset
SupportColleagues=genderCS[,c("Support_Colleagues", "All_Gender_clean", "Supporter_CareerStage_clean")]
#clean
SupportColleagues=SupportColleagues[(SupportColleagues$Support_Colleagues %in% seq(1:5)),]
dim(SupportColleagues)[1]
#factorise
SupportColleagues$Support_Colleagues=as.factor(SupportColleagues$Support_Colleagues)
table(SupportColleagues$Support_Colleagues)
```

### Support colleagues x gender 
```{r}
exp(coef(polr(Support_Colleagues~All_Gender_clean, data=SupportColleagues, Hess=TRUE, method="logistic"))) #OR (men as baseline)
exp(confint(polr(Support_Colleagues~All_Gender_clean, data=SupportColleagues, Hess=TRUE, method="logistic"))) #CI
kruskal.test(Support_Colleagues ~ All_Gender_clean, data=SupportColleagues) #effect of gender on response
```

### Support colleagues x CS
```{r}
exp(coef(polr(Support_Colleagues~Supporter_CareerStage_clean, data=SupportColleagues, Hess=TRUE, method="logistic"))) #OR (PhD as baseline)
exp(confint(polr(Support_Colleagues~Supporter_CareerStage_clean, data=SupportColleagues, Hess=TRUE, method="logistic"))) #CI
kruskal.test(Support_Colleagues ~ Supporter_CareerStage_clean, data=SupportColleagues) #effect of CS on response
```

### Support colleagues x gender graph
```{r}
levels(SupportColleagues$Support_Colleagues)=c("Strongly disagree", "Somewhat disagree", "Neither agree or disagree", "Somewhat agree", "Strongly agree")
SupportColleagues$Support_Colleagues=factor(SupportColleagues$Support_Colleagues, levels=c("Strongly agree", "Somewhat agree", "Neither agree or disagree", "Somewhat disagree", "Strongly disagree"))
graphdata = SupportColleagues %>%
  group_by(All_Gender_clean, Support_Colleagues) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/SupportColleagues_gender.eps", width=1000, height=578)
ggplot(graphdata, aes(x=All_Gender_clean, y=perc, fill=Support_Colleagues)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=20))+
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="I felt supported or valued by my colleagues (e.g. other members of the lab) for the help I was providing") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
dev.off()
```

### Support colleagues x CS graph
```{r}
graphdata = SupportColleagues %>%
  group_by(Supporter_CareerStage_clean, Support_Colleagues) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/SupportColleagues_CS.eps", width=1000, height=578)
ggplot(graphdata, aes(x=Supporter_CareerStage_clean, y=perc, fill=Support_Colleagues)) +
  geom_bar(stat="identity") +
  theme_minimal() +
      theme(plot.title = element_text(hjust = 0.5), text=element_text(size=20))+
  scale_x_discrete(limits = rev(levels(SupportColleagues$Supporter_CareerStage_clean))) +
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="I felt supported or valued by my colleagues (e.g. other members of the lab) for the help I was providing") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
dev.off()
```

## Q20c Support Personal (Section 5, Panel 2)
Model: ordinal logistic regression
```{r}
#subset
SupportPersonal=genderCS[,c("Support_Personal", "All_Gender_clean", "Supporter_CareerStage_clean")]
#clean
SupportPersonal=SupportPersonal[(SupportPersonal$Support_Personal %in% seq(1:5)),]
dim(SupportPersonal)[1]
#factorise
SupportPersonal$Support_Personal=as.factor(SupportPersonal$Support_Personal)
table(SupportPersonal$Support_Personal)
```

### Support Personal x gender
```{r}
exp(coef(polr(Support_Personal~All_Gender_clean, data=SupportPersonal, Hess=TRUE, method="logistic"))) #OR (men as baseline)
kruskal.test(Support_Personal ~ All_Gender_clean, data=SupportPersonal) #effect of gender on response
```

### Support Personal x CS
```{r}
exp(coef(polr(Support_Personal~Supporter_CareerStage_clean, data=SupportPersonal, Hess=TRUE, method="logistic"))) #OR (PhD as baseline)
exp(confint(polr(Support_Personal~Supporter_CareerStage_clean, data=SupportPersonal, Hess=TRUE, method="logistic"))) #CI
kruskal.test(Support_Personal ~ Supporter_CareerStage_clean, data=SupportPersonal) #effect of CS on response
```

### Support Personal x gender graph
```{r}
levels(SupportPersonal$Support_Personal)=c("Strongly disagree", "Somewhat disagree", "Neither agree or disagree", "Somewhat agree", "Strongly agree")
SupportPersonal$Support_Personal=factor(SupportPersonal$Support_Personal, levels=c("Strongly agree", "Somewhat agree", "Neither agree or disagree", "Somewhat disagree", "Strongly disagree"))

# This will save the image to your local code folder
#eps("images/SupportPersonal_gender.eps", width=1000, height=578)
graphdata = SupportPersonal %>%
  group_by(All_Gender_clean, Support_Personal) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))
ggplot(graphdata, aes(x=All_Gender_clean, y=perc, fill=Support_Personal)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=20))+
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="I felt supported or valued by my friends, partners and family outside of academia for the help I was providing") +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
dev.off()
```

### Support Personal x CS graph
```{r}
graphdata = SupportPersonal %>%
  group_by(Supporter_CareerStage_clean, Support_Personal) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/SupportPersonal_CS.eps", width=1000, height=578)
ggplot(graphdata, aes(x=Supporter_CareerStage_clean, y=perc, fill=Support_Personal)) +
  geom_bar(stat="identity") +
  theme_minimal() +
      theme(plot.title = element_text(hjust = 0.5), text=element_text(size=20))+
  scale_x_discrete(limits = rev(levels(SupportPersonal$Supporter_CareerStage_clean))) +
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="I felt supported or valued by my friends, partners and family outside of academia for the help I was providing") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)
dev.off()
```

## Q21 Need and access to emotional support (Section 5, Panel 1)

Filter down to people who picked at least 1 of the 6 meaningful options (i.e. not PNTA or unsure):

* Support_Emotional_No
* Support_Emotional_NoFind
* Support_Emotional_Institution
* Support_Emotional_Colleagues
* Support_Emotional_Personal
* Support_Emotional_Other_Dummy

```{r}
Support_Emotional=genderCS[!is.na(genderCS$Support_Emotional_No==1 | genderCS$Support_Emotional_NoFind==1 | genderCS$Support_Emotional_Institution==1 | genderCS$Support_Emotional_Colleagues==1 | genderCS$Support_Emotional_Personal==1 | genderCS$Support_Emotional_Other_Dummy==1),]
Support_Emotional=Support_Emotional[,c("Support_Emotional_No", "Support_Emotional_NoFind", "Support_Emotional_Institution","Support_Emotional_Colleagues" , "Support_Emotional_Personal", "Supporter_CareerStage_clean","All_Gender_clean")]
dim(Support_Emotional)[1]
#fill all na with 0
Support_Emotional[is.na(Support_Emotional)]=0
```
### Support_Emotional_No (Yes/No) x Gender // xCS 
Dependent variable - 

* 0 - didn't select this option 
* 1 - did select this option - i.e. DON'T NEED EMOTIONAL SUPPORT

Model: Factorial logistical regression

```{r}
#LR gender
summary(glm(Support_Emotional_No~All_Gender_clean, data=Support_Emotional, family="binomial")) #glm and p-value
exp(coef(glm(Support_Emotional_No~All_Gender_clean, data=Support_Emotional, family="binomial"))) #OR (men as baseline)
exp(confint(glm(Support_Emotional_No~All_Gender_clean, data=Support_Emotional, family="binomial"))) #CI

#LR CS
summary(glm(Support_Emotional_No~Supporter_CareerStage_clean, data=Support_Emotional, family="binomial")) #glm and p-value
exp(coef(glm(Support_Emotional_No~Supporter_CareerStage_clean, data=Support_Emotional, family="binomial"))) #OR (PhDs as baseline)
exp(confint(glm(Support_Emotional_No~Supporter_CareerStage_clean, data=Support_Emotional, family="binomial"))) #CI

#FLR gender-CS interactions
summary(glm(Support_Emotional_No~Supporter_CareerStage_clean*All_Gender_clean, data=Support_Emotional, family="binomial"))
```

### Support_Emotional_NoFind x Gender // xCS
As above

```{r}
#LR gender
summary(glm(Support_Emotional_NoFind~All_Gender_clean, data=Support_Emotional, family="binomial")) #glm and p-value
exp(coef(glm(Support_Emotional_NoFind~All_Gender_clean, data=Support_Emotional, family="binomial"))) #OR (men as baseline)
exp(confint(glm(Support_Emotional_NoFind~All_Gender_clean, data=Support_Emotional, family="binomial")))#CI

#LR CS
summary(glm(Support_Emotional_NoFind~Supporter_CareerStage_clean, data=Support_Emotional, family="binomial")) #glm and p-value
exp(coef(glm(Support_Emotional_NoFind~Supporter_CareerStage_clean, data=Support_Emotional, family="binomial"))) #OR (PhD as baseline)
exp(confint(glm(Support_Emotional_NoFind~Supporter_CareerStage_clean, data=Support_Emotional, family="binomial"))) #CI

#FLR gender-CS interactions
summary(glm(Support_Emotional_NoFind~Supporter_CareerStage_clean*All_Gender_clean, data=Support_Emotional, family="binomial"))
```

### Support_Emotional_Institution x Gender // xCS
As above
```{r}
#LR gender
summary(glm(Support_Emotional_Institution~All_Gender_clean, data=Support_Emotional, family="binomial")) #glm and p-value

#LR CS
summary(glm(Support_Emotional_Institution~Supporter_CareerStage_clean, data=Support_Emotional, family="binomial")) #glm and p-value
exp(coef(glm(Support_Emotional_Institution~Supporter_CareerStage_clean, data=Support_Emotional, family="binomial"))) #OR (PhD as baseline)
exp(confint(glm(Support_Emotional_Institution~Supporter_CareerStage_clean, data=Support_Emotional, family="binomial"))) #CI

#FLR gender-CS interactions
summary(glm(Support_Emotional_Institution~Supporter_CareerStage_clean*All_Gender_clean, data=Support_Emotional, family="binomial")) #glm and p-value
```

### Support_Emotional_Colleagues x Gender // xCS
As above
```{r}
#LR gender
summary(glm(Support_Emotional_Colleagues~All_Gender_clean, data=Support_Emotional, family="binomial")) #glm and p-value
exp(coef(glm(Support_Emotional_Colleagues~All_Gender_clean, data=Support_Emotional, family="binomial"))) #OR (men as baseline)
exp(confint(glm(Support_Emotional_Colleagues~All_Gender_clean, data=Support_Emotional, family="binomial"))) #CI

#LR CS
summary(glm(Support_Emotional_Colleagues~Supporter_CareerStage_clean, data=Support_Emotional, family="binomial")) #glm and p-value
exp(coef(glm(Support_Emotional_Colleagues~Supporter_CareerStage_clean, data=Support_Emotional, family="binomial"))) #OR (PhD as baseline)
exp(confint(glm(Support_Emotional_Colleagues~Supporter_CareerStage_clean, data=Support_Emotional, family="binomial"))) #CI

#FLR gender-CS interactions
summary(glm(Support_Emotional_Colleagues~Supporter_CareerStage_clean*All_Gender_clean, data=Support_Emotional, family="binomial")) #glm and p-values
```

### Support_Emotional_Personal x Gender // xCS
As above
```{r}
#LR gender
summary(glm(Support_Emotional_Personal~All_Gender_clean, data=Support_Emotional, family="binomial")) #glm and p-value
exp(coef(glm(Support_Emotional_Personal~All_Gender_clean, data=Support_Emotional, family="binomial"))) #OR (men as baseline)
exp(confint(glm(Support_Emotional_Personal~All_Gender_clean, data=Support_Emotional, family="binomial"))) #CI

#LR CS
summary(glm(Support_Emotional_Personal~Supporter_CareerStage_clean, data=Support_Emotional, family="binomial")) #glm and p-value
exp(coef(glm(Support_Emotional_Personal~Supporter_CareerStage_clean, data=Support_Emotional, family="binomial"))) #OR (PhD as baseline)
exp(confint(glm(Support_Emotional_Personal~Supporter_CareerStage_clean, data=Support_Emotional, family="binomial"))) #CI

#FLR gender-CS interactions
summary(glm(Support_Emotional_Personal~Supporter_CareerStage_clean*All_Gender_clean, data=Support_Emotional, family="binomial")) #glm and p-values
```

## Q15 Supporter action advocacy (Section 3, Panel 3)

Answered to Q15: When supporting the person (select all that applies), I advocated for the person within my institution

Exclude everyone that only answered PNTA and/or unsure, or did not check any options at all

```{r}
#subset
Supporter_Actions_Advocacy=genderCS[!is.na(genderCS$Supporter_Actions_OpenUp==1|genderCS$Supporter_Actions_Resources==1|genderCS$Supporter_Actions_HelpResources==1|genderCS$Supporter_Actions_Emotional==1|genderCS$Supporter_Actions_Advice==1|genderCS$Supporter_Actions_Work==1|genderCS$Supporter_Actions_Personal==1|genderCS$Supporter_Actions_Crisis==1|genderCS$Supporter_Actions_Adjustments==1|genderCS$Supporter_Actions_Advocacy==1|genderCS$Supporter_Actions_Other_Dummy==1),c("Supporter_Actions_Advocacy", "Supporter_CareerStage_clean", "All_Gender_clean")]

dim(Supporter_Actions_Advocacy)[1] #N
#fill all na with 0
Supporter_Actions_Advocacy[is.na(Supporter_Actions_Advocacy)]=0

#factorise
Supporter_Actions_Advocacy$Supporter_Actions_Advocacy=as.factor(Supporter_Actions_Advocacy$Supporter_Actions_Advocacy)
table(Supporter_Actions_Advocacy$Supporter_Actions_Advocacy)
```

### Supporter_Action_Advocacy x gender // x CS

```{r}
#LR gender
summary(glm(Supporter_Actions_Advocacy~All_Gender_clean, data=Supporter_Actions_Advocacy, family="binomial")) #glm and p-value
exp(coef(glm(Supporter_Actions_Advocacy~All_Gender_clean, data=Supporter_Actions_Advocacy, family="binomial"))) #OR (men as baseline)
exp(confint(glm(Supporter_Actions_Advocacy~All_Gender_clean, data=Supporter_Actions_Advocacy, family="binomial"))) #CI

#LR CS
summary(glm(Supporter_Actions_Advocacy~Supporter_CareerStage_clean, data=Supporter_Actions_Advocacy, family="binomial")) #glm and p-value
exp(coef(glm(Supporter_Actions_Advocacy~Supporter_CareerStage_clean, data=Supporter_Actions_Advocacy, family="binomial"))) #OR (PhDs as baseline)
exp(confint(glm(Supporter_Actions_Advocacy~Supporter_CareerStage_clean, data=Supporter_Actions_Advocacy, family="binomial"))) #CI

#FLR gender-CS interactions
summary(glm(Supporter_Actions_Advocacy~Supporter_CareerStage_clean*All_Gender_clean, data=Supporter_Actions_Advocacy, family="binomial"))
```



## Q22 Needing practical information (Section 6, Panel 2)
Dependent variable: Support_Practical

* 1: No, did not need info
* 2: No, already had info
* 3: Yes, need info

(remove 4- not sure, 5 - PNTA and 6- other)

3 approaches:

* chi square (but cannot check for intersect)
* turn the 3 options into 3 x [yes/no]s
* multinomial logistic regression

FINAL APPROACH (immediately following here)
recode dependent variable:
- No - 1 and 2
- Yes - 3
Run as binomial logistic regression
### Needing practical information x Gender x CS xGenderxCS (Section 6, Panel 2)
```{r}
#Subset
Support_Practical=genderCS[,c("Support_Practical", "Supporter_CareerStage_clean", "All_Gender_clean")]
#Recode 
Support_Practical=Support_Practical[Support_Practical$Support_Practical %in% seq(1,3),] #remove options 4&5
Support_Practical[Support_Practical$Support_Practical %in% c(1,2), "Support_Practical"]=0
Support_Practical[Support_Practical$Support_Practical==3, "Support_Practical"]=1
# LR gender
summary(glm(Support_Practical~All_Gender_clean, data=Support_Practical, family="binomial")) #glm and p-value
exp(coef(glm(Support_Practical~All_Gender_clean, data=Support_Practical, family="binomial"))) #OR (men as baseline)
exp(confint(glm(Support_Practical~All_Gender_clean, data=Support_Practical, family="binomial"))) #CI
# LR CS
summary(glm(Support_Practical~Supporter_CareerStage_clean, data=Support_Practical, family="binomial")) #glm and p-value
exp(coef(glm(Support_Practical~Supporter_CareerStage_clean, data=Support_Practical, family="binomial")))
exp(confint(glm(Support_Practical~Supporter_CareerStage_clean, data=Support_Practical, family="binomial")))
#FLR gender-CS interactions
summary(glm(Support_Practical~Supporter_CareerStage_clean*All_Gender_clean, data=Support_Practical, family="binomial")) #glm and p-value
### Needing practical information x CS pairwise
```{r}
# Compare between early PIs and other CSs only
## this requires the package "MHTdiscrete"
CSpairs=cbind(rep("Group leaders (<5yr)", 4), levels(ImpactDraining$Supporter_CareerStage_clean)[c(1,2,4,5)])
Support_PracticalCSKW=data.frame(CS1=c(), CS2=c(), X2=c(), N=c(), df=c(), pVal=c()) 
for (i in 1:dim(CSpairs)[1]) {
  comparedData=Support_Practical[(Support_Practical$Supporter_CareerStage_clean==CSpairs[i,1]|Support_Practical$Supporter_CareerStage_clean==CSpairs[i,2]),]
  comparedData$Support_Practical=as.factor(comparedData$Support_Practical)
  comparedData$Supporter_CareerStage_clean=droplevels(comparedData$Supporter_CareerStage_clean) #drop levels to avoid incorrect chisq approx
 # KW=kruskal.test(Support_Practical ~ Supporter_CareerStage_clean, data=comparedData) #effect of CS on response
  chisq=chisq.test(table(comparedData[,c("Support_Practical", "Supporter_CareerStage_clean")]))
#  Support_PracticalCSKW=rbind(Support_PracticalCSKW, data.frame(CS1=CSpairs[i,1], CS2=CSpairs[i,2], X2=as.numeric(KW$statistic), df=as.numeric(KW$parameter), pVal=as.numeric(KW$p.value)))
    Support_PracticalCSKW=rbind(Support_PracticalCSKW, data.frame(CS1=CSpairs[i,1], CS2=CSpairs[i,2], X2=as.numeric(chisq$statistic), N=dim(comparedData)[1], df=as.numeric(chisq$parameter), pVal=as.numeric(chisq$p.value)))
}
Support_PracticalCSKW$adjpVal=Sidak.p.adjust(Support_PracticalCSKW$pVal)
Support_PracticalCSKW
Support_Practical=genderCS[,c("Support_Practical", "Supporter_CareerStage_clean", "All_Gender_clean")]
chisq.test(table(Support_Practical[,c("Support_Practical", "All_Gender_clean")]))
chisq.test(table(Support_Practical[,c("Support_Practical", "Supporter_CareerStage_clean")]))
```


### Approach 2: turn the 3 options into 3 x [yes/no]s [not reported]
```{r}
#Subset
Support_Practical=genderCS[,c("Support_Practical", "Supporter_CareerStage_clean", "All_Gender_clean")]
#Factorise
Support_Practical$Support_Practical=as.factor(Support_Practical$Support_Practical)
#recode and clean
Support_Practical=cbind(sapply(levels(Support_Practical$Support_Practical), function(x) as.integer(x==Support_Practical$Support_Practical)), Support_Practical[2:3])
Support_Practical=Support_Practical[,c(2,3,4,7,8)] #remove some options
#rename columns
colnames(Support_Practical)=c("NoNeed", "Had", "Need", "Supporter_CareerStage_clean", "All_Gender_clean")
```


#### Choosing option 1 - No, I did not need the information
```{r}
# LR gender
summary(glm(NoNeed~All_Gender_clean, data=Support_Practical, family="binomial")) #glm and p-value
exp(coef(glm(NoNeed~All_Gender_clean, data=Support_Practical, family="binomial"))) #OR (men as baseline)
exp(confint(glm(NoNeed~All_Gender_clean, data=Support_Practical, family="binomial"))) #CI

# LR CS
summary(glm(NoNeed~Supporter_CareerStage_clean, data=Support_Practical, family="binomial")) #glm and p-value
exp(coef(glm(NoNeed~Supporter_CareerStage_clean, data=Support_Practical, family="binomial")))
exp(confint(glm(NoNeed~Supporter_CareerStage_clean, data=Support_Practical, family="binomial")))

#FLR gender-CS interactions
summary(glm(NoNeed~Supporter_CareerStage_clean*All_Gender_clean, data=Support_Practical, family="binomial")) #glm and p-value
```

#### Choosing option 2 - No, I have already received information before
```{r}
# LR gender
summary(glm(Had~All_Gender_clean, data=Support_Practical, family="binomial")) #glm and p-value
exp(coef(glm(Had~All_Gender_clean, data=Support_Practical, family="binomial"))) #OR (men as baseline)
exp(confint(glm(Had~All_Gender_clean, data=Support_Practical, family="binomial"))) #CI

#LR CS
summary(glm(Had~Supporter_CareerStage_clean, data=Support_Practical, family="binomial")) #glm and p-value
exp(coef(glm(Had~Supporter_CareerStage_clean, data=Support_Practical, family="binomial"))) #OR (PhD as baseline)
exp(confint(glm(Had~Supporter_CareerStage_clean, data=Support_Practical, family="binomial"))) #CI

#FLR gender-CS interactions
summary(glm(Had~Supporter_CareerStage_clean*All_Gender_clean, data=Support_Practical, family="binomial")) #glm and p-value
```

#### Choosing option 3 - Yes I needed information
```{r}
#LR Gender
summary(glm(Need~All_Gender_clean, data=Support_Practical, family="binomial")) #glm and p-value
exp(coef(glm(Need~All_Gender_clean, data=Support_Practical, family="binomial"))) #OR (men as baseline)
exp(confint(glm(Need~All_Gender_clean, data=Support_Practical, family="binomial"))) #CI

#LR CS
summary(glm(Need~Supporter_CareerStage_clean, data=Support_Practical, family="binomial")) #glm and p-value
exp(coef(glm(Need~Supporter_CareerStage_clean, data=Support_Practical, family="binomial"))) #OR (PhD as baseline)
exp(confint(glm(Need~Supporter_CareerStage_clean, data=Support_Practical, family="binomial"))) #CI

#FLR gender-CS interactions
summary(glm(Need~Supporter_CareerStage_clean*All_Gender_clean, data=Support_Practical, family="binomial")) #glm and p-value

```

### Needing practical information x gender graph
```{r}
Support_Practical=genderCS[,c("Support_Practical", "Supporter_CareerStage_clean", "All_Gender_clean")]
#Subset
Support_Practical=genderCS[,c("Support_Practical", "Supporter_CareerStage_clean", "All_Gender_clean")]
#Recode 
Support_Practical=Support_Practical[Support_Practical$Support_Practical %in% seq(1,3),] #remove options 4&5
Support_Practical[Support_Practical$Support_Practical %in% c(1,2), "Support_Practical"]="No"
Support_Practical[Support_Practical$Support_Practical==3, "Support_Practical"]="Yes"

#ImpactPersonal$Impact_Personal=factor(ImpactPersonal$Impact_Personal, levels=c("No","Yes"))
graphdata =Support_Practical %>%
  group_by(All_Gender_clean, Support_Practical) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/SupportPractical_gender.eps", width=1000, height=578)
ggplot(graphdata, aes(x=All_Gender_clean, y=perc, fill=Support_Practical)) +
  geom_bar(stat="identity") +
  theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16), legend.position="bottom", legend.text = element_text(size=7))+
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="
       I needed practical information \n(e.g. advice, guidelines, training...) \nto help me in my supporting role.") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=c(eBlue, eGreen))
ggsave(dpi=1000, "6-2-2 (not reported).png", limitsize = FALSE)
dev.off()
```

### Needing practical information x CS graph (Graph 6.2.2)

```{r}
graphdata = Support_Practical %>%
  group_by(Supporter_CareerStage_clean, Support_Practical) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

# This will save the image to your local code folder
#eps("images/SupportPractical_CS.eps", width=1000, height=578)
ggplot(graphdata, aes(x=Supporter_CareerStage_clean, y=perc, fill=Support_Practical)) +
  geom_bar(stat="identity") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5), text=element_text(size=12), axis.title.x = element_text(size = 8), axis.title.y = element_text(size = 16), legend.text = element_text(size=7))+
  geom_text(aes(label=round(perc, digit=1)), size=3, position=position_stack(vjust=0.5), color="white") +
  scale_x_discrete(limits = rev(levels(graphdata$Supporter_CareerStage_clean))) +
  labs(x="", y="Percentage", title="I needed practical information \n(e.g. advice, guidelines, training...) to help me \nin my supporting role.") + 
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=c(eBlue, eGreen))
ggsave(dpi=1000, "6-2-2.png", limitsize = FALSE)
dev.off()
```

### Model: Multinomial logistic regression (not reported)
Makes the most sense to use "Yes, I needed information" (option 3) as the baseline. This tests whether gender/CS has an effect of a respondent choosing the other two responses ("No, I did not need information" and "No, I had already received information before") over the baseline. 

```{r Support_Practical}
#subset
Support_Practical=genderCS[,c("Support_Practical", "Supporter_CareerStage_clean", "All_Gender_clean")]
Support_Practical=Support_Practical[Support_Practical$Support_Practical %in% seq(1,3), ] #remove NAs and some choices
dim(Support_Practical)[1] #N
Support_Practical$Support_Practical=as.factor(Support_Practical$Support_Practical) #factorise
#reorder factor to make "Yes, I needed information" the baseline
Support_Practical$Support_Practical=relevel(Support_Practical$Support_Practical, "3")
levels(Support_Practical$Support_Practical)= c("Yes, I needed information", "No, I did not need information", "No, I had already received information before")

#multinomial model -> p-values (2-tailed z test)

# Effect of gender
SP_gender=multinom(Support_Practical~All_Gender_clean, data=Support_Practical)
summary(SP_gender)
# Z-statistics and p-value using Wald tests
SP_gender_z=summary(SP_gender)$coefficients/summary(SP_gender)$standard.errors #Z statistics
(1-pnorm(abs(SP_gender_z), 0, 1))*2 #p-values
exp(coef(SP_gender)) #effect size (baseline IV is men, and baseline response is "yes")
exp(confint(SP_gender)) 

# Effect of career stage
SP_CS=multinom(Support_Practical~Supporter_CareerStage_clean, data=Support_Practical)
summary(SP_CS)
# Z-statistics and p-value using Wald tests
SP_CS_z=summary(SP_CS)$coefficients/summary(SP_CS)$standard.errors #Z statistics
(1-pnorm(abs(SP_CS_z), 0, 1))*2 #p-values
exp(coef(SP_CS)) #effect size
exp(confint(SP_CS))

```

### multinomial LR - Career stages pairwise analyses
using multinomial to compare between all career stages, pairwise

```{r}
#reverse CSpairs column order to not have early PIs on the baseline
CSpairs=CSpairs[,c(2,1)]
#new dataframe
SupportPracticalCSMLR=data.frame(CS1=c(), CS2=c(), response=c(), z=c(), edf=c(), pVal=c(), OR=c(), CIlower=c(), CIupper=c()) 
for (i in 1:dim(CSpairs)[1]) {
  comparedData=Support_Practical[(Support_Practical$Supporter_CareerStage_clean==CSpairs[i,1]|Support_Practical$Supporter_CareerStage_clean==CSpairs[i,2]),]
  comparedData$Supporter_CareerStage_clean=droplevels(comparedData$Supporter_CareerStage_clean) #drop levels to avoid levels with 0 values producing NaNs
  comparedData$Supporter_CareerStage_clean=relevel(comparedData$Supporter_CareerStage_clean, CSpairs[i,1]) #put the baseline first in factor level!
  model=multinom(Support_Practical~Supporter_CareerStage_clean, data=comparedData)
  # Z-statistics and p-value using Wald tests
  z=summary(model)$coefficients/summary(model)$standard.errors #Z statistics- from this line onwards, code produces 4 values, one for each non-baseline response, and intercept and CS2 over CS1. We ignore the intercept values (recording column 2)
  pVal=(1-pnorm(abs(z), 0, 1))*2 #p-values
  or=exp(coef(model)) #effect size
  CI=exp(confint(model)) #CI - [3,4] for first repsonse, [7,8] for second
  SupportPracticalCSMLR=rbind(SupportPracticalCSMLR, 
                              data.frame(CS1=CSpairs[i,1], CS2=CSpairs[i,2], response=rownames(z)[1], z=z[1,2], df=model$edf, pVal=pVal[1,2], OR=or[1,2], CIlower=CI[3], CIupper=CI[4]),
                              data.frame(CS1=CSpairs[i,1], CS2=CSpairs[i,2], response=rownames(z)[2], z=z[2,2], df=model$edf, pVal=pVal[2,2], OR=or[2,2], CIlower=CI[7], CIupper=CI[8]))
}

SupportPracticalCSMLR
```


## Q32c All_MH

### Q32c All_MH x gender x Supporter Status (no career stage here)

Looking at 2 independent variables:

1. Gender (label: "All_Gender_clean"): Man -1, Woman -3
  
2. Supporter/ Not supporter (label "Provided_Support"):

recode to new variable "Provided_Support_clean"

* 1- No - recode to "No"
* 2 - Yes, to one individual - recode to "Yes"
* 3 - Yes, to between one and five individuals - recode to "Yes"
* 4 - Yes, to more than five individuals - recode to "Yes"
* 5 - I am not sure - remove
* 6 - Prefer not to answer - remove

```{r}
#Filtering by gender and supporter status
genderSupport=data[data$All_Gender_clean %in% c(1,3) & data$Provided_Support %in% seq(1,4),]

#recode Supporter/not Supporter - Provided_Support_clean
genderSupport$Provided_Support_clean=genderSupport$Provided_Support
genderSupport[genderSupport$Provided_Support==1, "Provided_Support_clean"]="No"
genderSupport[genderSupport$Provided_Support %in% seq(2,4), "Provided_Support_clean"]="Yes"

genderSupport$All_Gender_clean=as.factor(genderSupport$All_Gender_clean)
genderSupport$Provided_Support_clean=as.factor(genderSupport$Provided_Support_clean)

levels(genderSupport$All_Gender_clean)=c("Men", "Women")
dim(genderSupport)
``` 

## Q01xQ28 Supporter vs gender (Section 2, Panel 1)
Chi-square & Logistical regression

```{r}
gender=genderSupport[,c("Provided_Support_clean","All_Gender_clean")]
dim(gender)[1] #this gives N
table(gender)
chisq.test(table(gender)) #not reported

summary(glm(Provided_Support_clean~All_Gender_clean, data=genderSupport, family="binomial")) 
exp(coef(glm(Provided_Support_clean~All_Gender_clean, data=genderSupport, family="binomial"))) #OR (men as baseline)
exp(confint(glm(Provided_Support_clean~All_Gender_clean, data=genderSupport, family="binomial"))) #CI

```



### filtering by response to All_MH
ordinal logistical regression: 1-5 strongly disagree - strongly agree

```{r}
AllMH=genderSupport[,c("All_MH", "All_Gender_clean", "Provided_Support_clean")]
#clean
AllMH=AllMH[(AllMH$All_MH %in% seq(1:5)),]
dim(AllMH)[1]
#factorise
AllMH$All_MH=as.factor(AllMH$All_MH)
table(AllMH$All_MH)
``` 

### All_MH x gender (Section 1, last paragraph)
```{r}
exp(coef(polr(All_MH~All_Gender_clean, data=AllMH, Hess=TRUE, method="logistic"))) #OR (men as baseline)
exp(confint(polr(All_MH~All_Gender_clean, data=AllMH, Hess=TRUE, method="logistic"))) #CI
kruskal.test(All_MH ~ All_Gender_clean, data=AllMH) #effect of gender on response
```

### All_MH x Supporter
```{r}
exp(coef(polr(All_MH~Provided_Support_clean, data=AllMH, Hess=TRUE, method="logistic"))) #OR (non-supporters are baseline)
kruskal.test(All_MH ~Provided_Support_clean, data=AllMH) #effect of supporter status on response
```

### All_MH x gender graph
```{r}
##reordering factors order for the graph
levels(AllMH$All_MH)=c("Strongly disagree", "Somewhat disagree", "Neither agree or disagree", "Somewhat agree", "Strongly agree")
AllMH$All_MH=factor(AllMH$All_MH, levels=c("Strongly agree", "Somewhat agree", "Neither agree or disagree", "Somewhat disagree", "Strongly disagree"))

#graph split by gender
graphdata = AllMH %>%
  group_by(All_Gender_clean, All_MH) %>%
  summarize(n=n()) %>%
  mutate(perc=n*100/sum(n))

#tiff("AllMH_gender.tiff", width=1000, height=578)
ggplot(graphdata, aes(x=All_Gender_clean, y=perc, fill=All_MH)) +
  geom_bar(stat="identity") +
  theme_minimal() +
     theme(plot.title = element_text(hjust = 0.5), text=element_text(size=20))+
  geom_text(aes(label=round(perc, digit=1)), size=4, position=position_stack(vjust=0.5), color="white") +
  labs(x="", y="Percentage", title="I have gone through times in my life when my mental health was poor") +
    guides(fill=guide_legend(reverse=TRUE)) +
  coord_flip() +
  scale_fill_manual(name="", values=ePalette)

#dev.off()


```

### All_MH x Supporter
```{r}
exp(coef(polr(All_MH~Provided_Support_clean, data=AllMH, Hess=TRUE, method="logistic"))) #OR (non-supporters are baseline)
kruskal.test(All_MH ~Provided_Support_clean, data=AllMH) #effect of supporter status on response
```

### Chisq for intersections [NOT REPORTED]
```{r}
# #split by supporter status, gender effect
# summary(AllMH$Provided_Support)
# chisq.test(table(AllMH[AllMH$Provided_Support_clean==levels(AllMH$Provided_Support_clean)[1],c("All_MH","All_Gender_clean")]))
# chisq.test(table(AllMH[AllMH$Provided_Support_clean==levels(AllMH$Provided_Support_clean)[2],c("All_MH","All_Gender_clean")]))
# #split by gender, supporter status (same as KW above)
# chisq.test(table(AllMH[AllMH$All_Gender_clean==levels(AllMH$All_Gender_clean)[1],c("All_MH","Provided_Support_clean")]))
# chisq.test(table(AllMH[AllMH$All_Gender_clean==levels(AllMH$All_Gender_clean)[2],c("All_MH","Provided_Support_clean")]))
```


## Q32d All_UnderstandingMH (Section 1, last paragraph)
ordinal logistical regression: 1-5 strongly disagree - strongly agree

```{r}
AllUnderstandingMH=genderSupport[,c("All_UnderstandingMH", "All_Gender_clean", "Provided_Support_clean")]
#clean
AllUnderstandingMH=AllUnderstandingMH[(AllUnderstandingMH$All_UnderstandingMH %in% seq(1:5)),]
dim(AllUnderstandingMH)[1]
#factorise
AllUnderstandingMH$All_UnderstandingMH=as.factor(AllUnderstandingMH$All_UnderstandingMH)
table(AllUnderstandingMH$All_UnderstandingMH)
``` 

### All_UnderstandingMH x gender

```{r}
exp(coef(polr(All_UnderstandingMH~All_Gender_clean, data=AllUnderstandingMH, Hess=TRUE, method="logistic"))) #OR (men as baseline)
exp(confint(polr(All_UnderstandingMH~All_Gender_clean, data=AllUnderstandingMH, Hess=TRUE, method="logistic"))) #CI
kruskal.test(All_UnderstandingMH ~ All_Gender_clean, data=AllUnderstandingMH) #effect of gender on response
```

### All_UnderstandingMH x Supporter
```{r}
exp(coef(polr(All_UnderstandingMH~Provided_Support_clean, data=AllUnderstandingMH, Hess=TRUE, method="logistic"))) #OR (non-supporters as baseline)
kruskal.test(All_UnderstandingMH ~Provided_Support_clean, data=AllUnderstandingMH) #effect of supporter status on response
```


<!-- ### Chisq for intersections [NOT REPORTED] -->
<!-- ```{r} -->
<!-- summary(AllUnderstandingMH$Provided_Support) -->
<!-- #split by supporter status, gender effect -->
<!-- chisq.test(table(AllUnderstandingMH[AllUnderstandingMH$Provided_Support_clean==levels(AllUnderstandingMH$Provided_Support_clean)[1],c("All_UnderstandingMH","All_Gender_clean")])) -->
<!-- chisq.test(table(AllUnderstandingMH[AllUnderstandingMH$Provided_Support_clean==levels(AllUnderstandingMH$Provided_Support_clean)[2],c("All_UnderstandingMH","All_Gender_clean")])) -->
<!-- #split by gender, supporter status (same as KW above) -->
<!-- chisq.test(table(AllUnderstandingMH[AllUnderstandingMH$All_Gender_clean==levels(AllUnderstandingMH$All_Gender_clean)[1],c("All_UnderstandingMH","Provided_Support_clean")])) -->
<!-- chisq.test(table(AllUnderstandingMH[AllUnderstandingMH$All_Gender_clean==levels(AllUnderstandingMH$All_Gender_clean)[2],c("All_UnderstandingMH","Provided_Support_clean")])) -->
<!-- ``` -->


### Q32c All_MH x career stage
Looking at the single indepedent variable - present career stage. This is because this section focuses on understanding the population that filled in the survey

1. CS (label: "All_CurrentRole_clean"): PhD - 3, postdoc -4, earlyPI - 6, midPI- 7, late PI -8
drop all other stages


```{r}

AllCS=data[data$All_CurrentRole_clean %in% c(3,4,6,7,8), ] #subset
AllCS$All_CurrentRole_clean=as.factor(AllCS$All_CurrentRole_clean) #factorise
dim(AllCS) 

levels(AllCS$All_CurrentRole_clean)=c("PhD students", "Postdocs", "Group leaders (<5yr)", "Group leaders (5-10yr)", "Group leaders (>10yr)")
```


Filter All_MH response as before
```{r}
#subset
AllMH_CS=AllCS[,c("All_MH", "All_CurrentRole_clean")]
#clean
AllMH_CS=AllMH_CS[(AllMH_CS$All_MH %in% seq(1:5)),]
dim(AllMH_CS)[1] #N
#factorise
AllMH_CS$All_MH=as.factor(AllMH_CS$All_MH)
table(AllMH_CS$All_MH)
``` 

### All_MH x CS (Section 1, last paragraph)
ordinal logistic regression, single Indepedent variable
```{r}
exp(coef(polr(All_MH~All_CurrentRole_clean, data=AllMH_CS, Hess=TRUE, method="logistic"))) #OR (men as baseline)
exp(confint(polr(All_MH~All_CurrentRole_clean, data=AllMH_CS, Hess=TRUE, method="logistic"))) #CI
kruskal.test(All_MH ~ All_CurrentRole_clean, data=AllMH_CS) #effect of CS on response
```

# Q1x Q32 Supporter vs CS
chi-square (not reported) and logistic regression
```{r}
#subset
CSSupport=AllCS[,c("Provided_Support","All_CurrentRole_clean")]
CSSupport$Provided_Support_clean=CSSupport$Provided_Support
CSSupport=CSSupport[!is.na(CSSupport$Provided_Support),] #remove NA
CSSupport[CSSupport$Provided_Support==1, 'Provided_Support_clean']=0
CSSupport[CSSupport$Provided_Support %in% seq(2,4), 'Provided_Support_clean']=1
CSSupport=CSSupport[CSSupport$Provided_Support_clean %in% c(0,1), c("Provided_Support_clean", "All_CurrentRole_clean")]

dim(CSSupport)[1] #this gives N

table(CSSupport)
chisq.test(table(CSSupport)) #not reported

summary(glm(Provided_Support_clean~All_CurrentRole_clean, data=CSSupport, family="binomial")) 
exp(coef(glm(Provided_Support_clean~All_CurrentRole_clean, data=CSSupport, family="binomial"))) #OR (PhD as baseline)
exp(confint(glm(Provided_Support_clean~All_CurrentRole_clean, data=CSSupport, family="binomial"))) #CI
```

## Q32d All_UnderstandingMH  x career stage (Section 1, last paragraph)
ordinal logistical regression: 1-5 strongly disagree - strongly agree

```{r}
AllUnderstandingMH_CS=AllCS[,c("All_UnderstandingMH", "All_CurrentRole_clean")]
#clean
AllUnderstandingMH_CS=AllUnderstandingMH_CS[(AllUnderstandingMH_CS$All_UnderstandingMH %in% seq(1:5)),]
dim(AllUnderstandingMH_CS)[1]
#factorise
AllUnderstandingMH_CS$All_UnderstandingMH=as.factor(AllUnderstandingMH_CS$All_UnderstandingMH)
table(AllUnderstandingMH_CS$All_UnderstandingMH)

exp(coef(polr(All_UnderstandingMH~All_CurrentRole_clean, data=AllUnderstandingMH_CS, Hess=TRUE, method="logistic"))) #OR (PhD as baseline)
exp(confint(polr(All_UnderstandingMH~All_CurrentRole_clean, data=AllUnderstandingMH_CS, Hess=TRUE, method="logistic"))) #CI
kruskal.test(All_UnderstandingMH ~ All_CurrentRole_clean, data=AllUnderstandingMH_CS) #effect of CS on response
```
