##### Installing Packages #####
install.packages("foreign")
install.packages("nFactors")
install.packages("psych")
install.packages("Hmisc")
install.packages("sjstats")
install.packages("lme4")
install.packages("lm.beta")
install.packages("effectsize")

##### Activating Packages #####
library(foreign)
library(nFactors)
library(psych)
library(Hmisc)
library(sjstats)
library(lme4)
library(lm.beta)
library(heplots)
library(effectsize)
library(lavaan)

##### Loading Data #####

#Setting Working directory
setwd("")#Please set working directory

#Loading Data into Work space
data<-read.spss("Reimagining Recruitment - complete_Coded.sav",
                use.value.labels = FALSE,
                to.data.frame = TRUE)

#Cutting data down to the relevant bits
Data<-data[,c(1,83,92,100,101,109,110,111,119,127:148,156:172,180:213,
              243:259,267,275:276,284:286,294:297,328:330,338,346:347)]

#Naming all Variables
colnames(Data)<-c("SN","University","Discipline","Role","Role_Text","Full_Time",
                  "Fixed_Term_Contract","Employment_Duration_Months",
                  "Nr_of_Contracts","Identity_1","Identity_2","Identity_3",
                  "Identity_4","Identity_5","Identity_6","Identity_7",
                  "Identity_8","Identity_9","Identity_10","Career_Attitude_1",
                  "Career_Attitude_2","Career_Attitude_3","Career_Attitude_4",
                  "Career_Attitude_5","Career_Attitude_6","Career_Attitude_7",
                  "Career_Attitude_8","Career_Attitude_9","Career_Attitude_10",
                  "Staying_in_Academia","Staying_in_Academia_text",
                  "WorkingStyle_1","WorkingStyle_2","WorkingStyle_3",
                  "WorkingStyle_4","WorkingStyle_5","WorkingStyle_6",
                  "WorkingStyle_7","WorkingStyle_8","WorkingStyle_9",
                  "ReceivedOpportunity_1","ReceivedOpportunity_2",
                  "ReceivedOpportunity_3","ReceivedOpportunity_4",
                  "ReceivedOpportunity_5","ReceivedOpportunity_6",
                  "DiversityLevel","DiversityLevel_text","DiversityClimate_1",
                  "DiversityClimate_2","DiversityClimate_3",
                  "DiversityClimate_4","DiversityClimate_5",
                  "DiversityClimate_6","DiversityClimate_7", "Affective_1",
                  "Affective_2","Affective_3","Affective_4","Affective_5",
                  "Affective_6","Affective_7","Affective_8",
                  "WorkplaceCulture_1","WorkplaceCulture_2",
                  "WorkplaceCulture_3","WorkplaceCulture_4",
                  "WorkplaceCulture_5","WorkplaceCulture_6",
                  "WorkplaceCulture_7","WorkplaceCulture_8",
                  "Harrassment_Self_1","Harrassment_Self_2",
                  "Harrassment_Self_3","Harrassment_Self_4",
                  "Harrassment_Self_5","Harrassment_Other_1",
                  "Harrassment_Other_2","Harrassment_Other_3",
                  "Harrassment_Other_4","Harrassment_Other_5",
                  "Harrassment_text","Offline_1","Offline_2","Offline_3",
                  "Offline_4","Offline_5","Offline_6","Offline_7","Offline_8",
                  "Online_1","Online_2","Online_3","Online_4","Online_5",
                  "Online_6","Online_7","Online_8","Offline_Online_text","Age",
                  "Sex_Gender","Sex_Gender_Text","Sex_Impact",
                  "Sexual_Orientation","Sexual_Orientation_Text",
                  "Sexual_Orientation_Impact","Parental_Education",
                  "Parental_Eduaction_Impact","Marital_Status","UK_Citizen",
                  "Ethnicity","Ethnicity_Text_Mixed","Ethnicity_Text_Other",
                  "Ethnicity_Impact","Disability")

#Create Variable Look-up Screen
VariableLookup<-matrix(nrow = 84,ncol = 2)
VariableLookup[1,1]<-"Identity_1"
VariableLookup[2,1]<-"Identity_2"
VariableLookup[3,1]<-"Identity_3"
VariableLookup[4,1]<-"Identity_4"
VariableLookup[5,1]<-"Identity_5"
VariableLookup[6,1]<-"Identity_6"
VariableLookup[7,1]<-"Identity_7"
VariableLookup[8,1]<-"Identity_8"
VariableLookup[9,1]<-"Identity_9"
VariableLookup[10,1]<-"Identity_10"
VariableLookup[11,1]<-"Career_Attitude_1"
VariableLookup[12,1]<-"Career_Attitude_2"
VariableLookup[13,1]<-"Career_Attitude_3"
VariableLookup[14,1]<-"Career_Attitude_4"
VariableLookup[15,1]<-"Career_Attitude_5"
VariableLookup[16,1]<-"Career_Attitude_6"
VariableLookup[17,1]<-"Career_Attitude_7"
VariableLookup[18,1]<-"Career_Attitude_8"
VariableLookup[19,1]<-"Career_Attitude_9"
VariableLookup[20,1]<-"Career_Attitude_10"
VariableLookup[21,1]<-"WorkingStyle_1"
VariableLookup[22,1]<-"WorkingStyle_2"
VariableLookup[23,1]<-"WorkingStyle_3"
VariableLookup[24,1]<-"WorkingStyle_4"
VariableLookup[25,1]<-"WorkingStyle_5"
VariableLookup[26,1]<-"WorkingStyle_6"
VariableLookup[27,1]<-"WorkingStyle_7"
VariableLookup[28,1]<-"WorkingStyle_8"
VariableLookup[29,1]<-"WorkingStyle_9"
VariableLookup[30,1]<-"ReceivedOpportunity_1"
VariableLookup[31,1]<-"ReceivedOpportunity_2"
VariableLookup[32,1]<-"ReceivedOpportunity_3"
VariableLookup[33,1]<-"ReceivedOpportunity_4"
VariableLookup[34,1]<-"ReceivedOpportunity_5"
VariableLookup[35,1]<-"ReceivedOpportunity_6"
VariableLookup[36,1]<-"DiversityClimate_1"
VariableLookup[37,1]<-"DiversityClimate_2"
VariableLookup[38,1]<-"DiversityClimate_3"
VariableLookup[39,1]<-"DiversityClimate_4"
VariableLookup[40,1]<-"DiversityClimate_5"
VariableLookup[41,1]<-"DiversityClimate_6"
VariableLookup[42,1]<-"DiversityClimate_7"
VariableLookup[43,1]<-"Affective_1"
VariableLookup[44,1]<-"Affective_2"
VariableLookup[45,1]<-"Affective_3"
VariableLookup[46,1]<-"Affective_4"
VariableLookup[47,1]<-"Affective_5"
VariableLookup[48,1]<-"Affective_6"
VariableLookup[49,1]<-"Affective_7"
VariableLookup[50,1]<-"Affective_8"
VariableLookup[51,1]<-"WorkplaceCulture_1"
VariableLookup[52,1]<-"WorkplaceCulture_2"
VariableLookup[53,1]<-"WorkplaceCulture_3"
VariableLookup[54,1]<-"WorkplaceCulture_4"
VariableLookup[55,1]<-"WorkplaceCulture_5"
VariableLookup[56,1]<-"WorkplaceCulture_6"
VariableLookup[57,1]<-"WorkplaceCulture_7"
VariableLookup[58,1]<-"WorkplaceCulture_8"
VariableLookup[59,1]<-"Harrassment_Self_1"
VariableLookup[60,1]<-"Harrassment_Self_2"
VariableLookup[61,1]<-"Harrassment_Self_3"
VariableLookup[62,1]<-"Harrassment_Self_4"
VariableLookup[63,1]<-"Harrassment_Self_5"
VariableLookup[64,1]<-"Harrassment_Other_1"
VariableLookup[65,1]<-"Harrassment_Other_2"
VariableLookup[66,1]<-"Harrassment_Other_3"
VariableLookup[67,1]<-"Harrassment_Other_4"
VariableLookup[68,1]<-"Harrassment_Other_5"
VariableLookup[68,1]<-"Harrassment_Other_5"
VariableLookup[69,1]<-"Offline_1"
VariableLookup[70,1]<-"Offline_2"
VariableLookup[71,1]<-"Offline_3"
VariableLookup[72,1]<-"Offline_4"
VariableLookup[73,1]<-"Offline_5"
VariableLookup[74,1]<-"Offline_6"
VariableLookup[75,1]<-"Offline_7"
VariableLookup[76,1]<-"Offline_8"
VariableLookup[77,1]<-"Online_1"
VariableLookup[78,1]<-"Online_2"
VariableLookup[79,1]<-"Online_3"
VariableLookup[80,1]<-"Online_4"
VariableLookup[81,1]<-"Online_5"
VariableLookup[82,1]<-"Online_6"
VariableLookup[83,1]<-"Online_7"
VariableLookup[84,1]<-"Online_8"
VariableLookup[1,2]<-"Being an academic is an important part of my self-concept"
VariableLookup[2,2]<-"I am typical of most academics"
VariableLookup[3,2]<-"I have a lot in common with the avaerage academic"
VariableLookup[4,2]<-"I feel a strong sense of belonging as an academic"
VariableLookup[5,2]<-"I sometimes feel on the periphery as an academic"
VariableLookup[6,2]<-"Other academics recognise me as a valued member"
VariableLookup[7,2]<-"Sometimes, I think other academics doubt my credentials"
VariableLookup[8,2]<-"I am considered to have academic expertise"
VariableLookup[9,2]<-"I can always be myself around other academics"
VariableLookup[10,2]<-"I sometimes hide aspects of my identity to avoid negative attention from academics"
VariableLookup[11,2]<-"I like this Career too much to give it up"
VariableLookup[12,2]<-"If I could do it all over again, I would choose to work in a different profession"
VariableLookup[13,2]<-"Most of my academic career has been a natural progression from one step to the next"
VariableLookup[14,2]<-"I stayed in academia because it was the easiest option/I couldn't think of another option"
VariableLookup[15,2]<-"Sometimes I doubt whether I can succeed in an academic career"
VariableLookup[16,2]<-"Someone like me can succeed in an academic career"
VariableLookup[17,2]<-"I am satisfied with my chances of getting ahead in an academic career "
VariableLookup[18,2]<-"Academic Recruitment Processes select the best person for the job"
VariableLookup[19,2]<-"From my observations, there is bias in academic recruitment processes"
VariableLookup[20,2]<-"From my observations, there is bias in academic promotion processes"
VariableLookup[21,2]<-"I prefer to work in collaboration with others"
VariableLookup[22,2]<-"I prefer to work in independently and alone"
VariableLookup[23,2]<-"I feel comfortable working in a group"
VariableLookup[24,2]<-"I can challenge the ideas of people I work with"
VariableLookup[25,2]<-"I find it easy to articulate my ideas to strangers"
VariableLookup[26,2]<-"I find it easy to articulate my ideas to people more senior than me"
VariableLookup[27,2]<-"I rarely get credit for my work or ideas"
VariableLookup[28,2]<-"I thrive in a competitive environment"
VariableLookup[29,2]<-"I thrive in a diverse and collaborative environment"
VariableLookup[30,2]<-"I have opportunities to collaborate on publications"
VariableLookup[31,2]<-"I have opportunities to be involved in research grant applications"
VariableLookup[32,2]<-"I have the opportunity to attend external activities (e.g., conference attendence, workshops, meetings)"
VariableLookup[33,2]<-"I have opportunities to build professional and social networks"
VariableLookup[34,2]<-"I am encouraged to take on career development opportunities"
VariableLookup[35,2]<-"I have a clear understanding of how to advance in my academic career"
VariableLookup[36,2]<-"There is a fair division of rewards and opportunitites"
VariableLookup[37,2]<-"There is a fair division of desirable and undesirable tasks"
VariableLookup[38,2]<-"People with different backgrounds are treated fairly"
VariableLookup[39,2]<-"There are role models from diverse backgrounds"
VariableLookup[40,2]<-"Overall, I think people are treated equally on their merits"
VariableLookup[41,2]<-"More needs to be done to encourage diversity and inclusiveness"
VariableLookup[42,2]<-"My line manager/supervisor makes sure everyone has the same support and opportunities"
VariableLookup[43,2]<-"Happy vs Sad"
VariableLookup[44,2]<-"Angry vs Calm"
VariableLookup[45,2]<-"Disappointed vs Pleased"
VariableLookup[46,2]<-"Safe vs Intimidated"
VariableLookup[47,2]<-"Insecure vs Confident"
VariableLookup[48,2]<-"Hopeful vs Hopeless"
VariableLookup[49,2]<-"Accepted vs Rejected"
VariableLookup[50,2]<-"Included vs Excluded"
VariableLookup[51,2]<-"There is a spirit of collegiallity and support"
VariableLookup[52,2]<-"Collaborative working is valued and supported"
VariableLookup[53,2]<-"People tend to be very competitive"
VariableLookup[54,2]<-"A good work-life balance is encouraged"
VariableLookup[55,2]<-"You need to be part of a clique or aligned with a powerful individual or group to have influence"
VariableLookup[56,2]<-"I feel that I am listened to when I give my opinions (e.g., in a workplace meeting or forum)"
VariableLookup[57,2]<-"Social activities (e.g., parties or team building) are welcoming to all"
VariableLookup[58,2]<-"Offensive language and behaviours are not accepted (e.g., inappropriate images or comments, unwanted or degrading behaviour or interactions)"
VariableLookup[59,2]<-"Offensive Remarks, Jokes, Harrassment, or Bullying based on gender or sexuality"
VariableLookup[60,2]<-"Offensive Remarks, Jokes, Harrassment, or Bullying based on Ethnicity"
VariableLookup[61,2]<-"Offensive Remarks, Jokes, Harrassment, or Bullying based on Socio-Economic Background"
VariableLookup[62,2]<-"Offensive Remarks, Jokes, Harrassment, or Bullying based on Disability"
VariableLookup[63,2]<-"Offensive Remarks, Jokes, Harrassment, or Bullying based on Religion"
VariableLookup[64,2]<-"Offensive Remarks, Jokes, Harrassment, or Bullying based on someone else's gender or sexuality"
VariableLookup[65,2]<-"Offensive Remarks, Jokes, Harrassment, or Bullying based on someone else's Ethnicity"
VariableLookup[66,2]<-"Offensive Remarks, Jokes, Harrassment, or Bullying based on someone else's Socio-Economic Background"
VariableLookup[67,2]<-"Offensive Remarks, Jokes, Harrassment, or Bullying based on someone else's Disability"
VariableLookup[68,2]<-"Offensive Remarks, Jokes, Harrassment, or Bullying based on someone else's Religion"
VariableLookup[69,2]<-"It has been easy to advance my research interests"
VariableLookup[70,2]<-"I have often felt excluded or on the periphery"
VariableLookup[71,2]<-"I have found it easy to express my optinion and be listened to"
VariableLookup[72,2]<-"Other participants generally take me seriously and make me feel valued"
VariableLookup[73,2]<-"I have found it difficult o meet people with whom I might collaborate"
VariableLookup[74,2]<-"I have enjoyed the opportunity to catch up with people I know"
VariableLookup[75,2]<-"I have enjoyed the opportunity to meet new people and develeop my networks"
VariableLookup[76,2]<-"The experience has made me feel more positive about my future in academia"
VariableLookup[77,2]<-"It has been easy to advance my research interests"
VariableLookup[78,2]<-"I have often felt excluded or on the periphery"
VariableLookup[79,2]<-"I have found it easy to express my optinion and be listened to"
VariableLookup[80,2]<-"Other participants generally take me seriously and make me feel valued"
VariableLookup[81,2]<-"I have found it difficult o meet people with whom I might collaborate"
VariableLookup[82,2]<-"I have enjoyed the opportunity to catch up with people I know"
VariableLookup[83,2]<-"I have enjoyed the opportunity to meet new people and develeop my networks"
VariableLookup[84,2]<-"The experience has made me feel more positive about my future in academia"

#Coding Factor (Ordinal or Nominal) Variables
#University
Data$University<-as.numeric(Data$University)
Data$University<-factor(Data$University,
                        levels = c(1:41),
                        labels = c("Bangor Univeristy","University of Bath",
                                   "University of Birmingham",
                                   "University of Bristol",
                                   "University of Cambridge",
                                   "Universidad Nacional de Colombia",
                                   "University of Edinburgh",
                                   "University of Exeter", 
                                   "Heriot Watt University",
                                   "University of Hertfordshire",
                                   "Imperial College London",
                                   "Kings College London",
                                   "Lancaster University","University of Leeds",
                                   "University of Leicester", 
                                   "University of Liverpool", 
                                   "Liverpool John Moores University",
                                   "London Metropolitan University",
                                   "University of Manchester",
                                   "Newcastle University",
                                   "University of Nottingham",
                                   "The Open University",
                                   "University of Oxford",
                                   "University of Plymouth",
                                   "University of Portsmouth",
                                   "Queen Mary University of London",
                                   "Queen's University Belfast",
                                   "University of Sheffield",
                                   "University of Southampton",
                                   "University of St. Andrews",
                                   "University of Stirling",
                                   "University of Strathclyde, Glasgow",
                                   "University of Surrey",
                                   "University of Sussex",
                                   "University of Swansea",
                                   "Trinity College Dublin",
                                   "Ulster University",
                                   "University College London",
                                   "University of Warwick",
                                   "University of the West of England",
                                   "University of York"))

#Discipline
Data$Discipline<-as.numeric(Data$Discipline)
Data$Discipline<-factor(Data$Discipline,
                        levels = c(1:8),
                        labels = c("Biological Sciences",
                                   "Computer Sciences",
                                   "Engineering",
                                   "Mathematical Sciences",
                                   "Physics",
                                   "Other Discipline",
                                   "Chemistry",
                                   "Earth Sciences"))

#Role
Data$Role<-factor(Data$Role,levels = c(1:11),
                  labels = c("PhD-Student",
                             "Post-Doc",
                             "Teaching Fellow",
                             "Research Fellow",
                             "Senior Teaching Fellow",
                             "Senior Research Fellow",
                             "Lecturer",
                             "Associate Professor/Reader",
                             "Professor",
                             "Other",
                             "Senior Lecturer"))

#Full Time
Data$Full_Time[which(Data$Full_Time==3)]<-3
Data$Full_Time<-factor(Data$Full_Time,levels = c(1:2),
                       labels = c("Full Time",
                                  "Part Time"))

#Contract
Data$Fixed_Term_Contract[which(Data$Fixed_Term_Contract>2)]<-NA
Data$Fixed_Term_Contract[which(Data$Fixed_Term_Contract==2)]<-0
Data$Fixed_Term_Contract<-factor(Data$Fixed_Term_Contract,
                                 levels = c(0,1),
                                 labels = c("Open Ended Contract",
                                            "Fixed Term Contract"))
#Employment Duration
Data$Employment_Duration_Months<-as.numeric(Data$Employment_Duration_Months)

#Number of Contracts
Data$Nr_of_Contracts<-as.numeric(Data$Nr_of_Contracts)

#Age
Data$Age<-as.numeric(Data$Age)

#Sex Gender
Data$Sex_Gender[which(Data$Sex_Gender>2)]<-3
Data$Sex_Gender[which(Data$Sex_Gender==2)]<-4
Data$Sex_Gender[which(Data$Sex_Gender==1)]<-2
Data$Sex_Gender[which(Data$Sex_Gender==4)]<-1
Data$Sex_Gender<-factor(Data$Sex_Gender, levels = c(1:3),
                        labels = c("Cis-Man","Cis-Woman","Non-Cis"))


#Sexual Orientation
Data$Sexual_Orientation[which(Data$Sexual_Orientation==7)]<-NA
Data$Sexual_Orientation[which(Data$Sexual_Orientation>2)]<-2
Data$Sexual_Orientation<-factor(Data$Sexual_Orientation,levels = c(1:2),
                                labels = c("Hetherosexual", 
                                           "Non-Hetherosexual"))

#Parental Orientation
Data$Parental_Education[which(Data$Parental_Education==6)]<-NA
Data$Parental_Education[which(Data$Parental_Education==1)]<-2
Data$Parental_Education[which(Data$Parental_Education>2)]<-1
Data$Parental_Education<-factor(Data$Parental_Education,
                                levels = c(1:2),
                                labels = c("Higher Education",
                                           "No Higher Education"))

#Marital Status
Data$Marital_Status[which(Data$Marital_Status==3)]<-NA
Data$Marital_Status<-factor(Data$Marital_Status,
                            levels = c(1:2),
                            labels = c("Not Married",
                                       "Married"))

#Citizenship
Data$UK_Citizen[which(Data$UK_Citizen==3)]<-NA
Data$UK_Citizen[which(Data$UK_Citizen==2)]<-0
Data$UK_Citizen<-factor(Data$UK_Citizen,
                        levels = c(0:1),
                        labels = c("Not a UK Citizen",
                                   "UK Citizen"))

#Ethnicity
Data$Ethnicity[which(Data$Ethnicity==10)]<-NA
Data$Ethnicity[which(Data$Ethnicity==2)]<-1
Data$Ethnicity[which(Data$Ethnicity>2)]<-2
Data$Ethnicity<-factor(Data$Ethnicity,
                       levels = c(1:2),
                       labels = c("White",
                                  "Non-White"))

#Disability
Data$Disability[which(Data$Disability==3)]<-NA
Data$Disability[which(Data$Disability==2)]<-0
Data$Disability<-factor(Data$Disability,
                        levels = c(0:1),
                        labels = c("No Disability",
                                   "Disability"))

#Career Stage
Data$Career_Stage[which(Data$Role=="PhD-Student")]<-1
Data$Career_Stage[which(Data$Role=="Post-Doc")]<-2
Data$Career_Stage[which(Data$Role=="Teaching Fellow")]<-3
Data$Career_Stage[which(Data$Role=="Research Fellow")]<-3
Data$Career_Stage[which(Data$Role=="Senior Teaching Fellow")]<-3
Data$Career_Stage[which(Data$Role=="Senior Research Fellow")]<-3
Data$Career_Stage[which(Data$Role=="Lecturer")]<-4
Data$Career_Stage[which(Data$Role=="Associate Professor/Reader")]<-5
Data$Career_Stage[which(Data$Role=="Professor")]<-5
Data$Career_Stage[which(Data$Role=="Other")]<-NA
Data$Career_Stage[which(Data$Role=="Senior Lecturer")]<-5

Data$Career_Stage<-factor(Data$Career_Stage,levels = c(1:5),
                          labels = c("PhD-Student",
                                     "Post-Doc",
                                     "Fellow",
                                     "Lecturer",
                                     "Senior Academic"))

Data$Staying_in_STEM_in_Academia[which(Data$Staying_in_Academia==1)]<-0
Data$Staying_in_STEM_in_Academia[which(Data$Staying_in_Academia==3)]<-0
Data$Staying_in_STEM_in_Academia[which(Data$Staying_in_Academia==2)]<-1
Data$Staying_in_STEM_in_Academia<-factor(Data$Staying_in_STEM_in_Academia,
                                         levels = c(0:1),
                                         labels = c("Not Staying",
                                                    "Staying"))

#Variance Check
Data$SumVar<-apply(Data[,c(10:30,32:47,49:81,83:98)], 1, var,na.rm=TRUE)

#Manual exclusion check
Data$Flag<-0
Data$Flag[which(Data$SN==4)]<-1 #Ignore Survey
Data$Flag[which(Data$SN==187)]<-1 #Manager - not Academic Staff
Data$Flag[which(Data$SN==498)]<-1 #Master's Student
Data$Flag[which(Data$SN==538)]<-1 #Lab Manager - Not Academic Staff
Data$Flag[which(Data$SN==629)]<-1 #Professional Services Staff
Data$Flag[which(Data$SN==637)]<-1 #Technical Services Manager
Data$Sexual_Orientation[which(Data$SN==9)]<-NA #Prefers not to Answer
Data$Sexual_Orientation[which(Data$SN==136)]<-NA #Prefers not to Answer
Data$Sexual_Orientation[which(Data$SN==167)]<-NA #Prefers not to Answer

DataAnalysis<-Data[which(Data$Career_Stage!="Fellow"& 
                           is.na(Data$SumVar)==FALSE &
                           Data$Flag==0),]

DataAnalysis<-Data[which(is.na(Data$SumVar)==FALSE &
                           Data$Flag==0),]

##### Factor Analyses####

#Identity Variables
#Reversing Inverse Items
DataAnalysis$Identity_5<-8-DataAnalysis$Identity_5
DataAnalysis$Identity_7<-8-DataAnalysis$Identity_7
DataAnalysis$Identity_10<-8-DataAnalysis$Identity_10

RV<-DataAnalysis[,c(10:19)]
RV<-RV[complete.cases(RV),]

ev <- eigen(cor(RV))
ap <- parallel(subject=nrow(RV),
               var=ncol(RV),
               rep=100,cent=.05)
nS <- nScree(x=ev$values, 
             aparallel=ap$eigen$qevpea)
plotnScree(nS)

numFact=nS$Components[[1]]
cormatrix<-cor(RV, method = "pearson")

fit <- fa(r=cormatrix, nfactors=numFact, rotate="varimax")
print(fit, digits=2, cutoff=.3, sort=TRUE)

#Thoughts about Academia
#Reversing Inverse Items
DataAnalysis$Career_Attitude_2<-8-DataAnalysis$Career_Attitude_2
DataAnalysis$Career_Attitude_5<-8-DataAnalysis$Career_Attitude_5
DataAnalysis$Career_Attitude_9<-8-DataAnalysis$Career_Attitude_9
DataAnalysis$Career_Attitude_10<-8-DataAnalysis$Career_Attitude_10

RV<-DataAnalysis[,c(20:29)]
RV<-RV[complete.cases(RV),]

ev <- eigen(cor(RV))
ap <- parallel(subject=nrow(RV),
               var=ncol(RV),
               rep=100,cent=.05)
nS <- nScree(x=ev$values, 
             aparallel=ap$eigen$qevpea)
plotnScree(nS)

numFact=nS$Components[[1]]
cormatrix<-cor(RV, method = "pearson")

fit <- fa(r=cormatrix, nfactors=numFact, rotate="varimax")
print(fit, digits=2, cutoff=.3, sort=TRUE)

#Working Style
#Reversing Inverse Items
DataAnalysis$WorkingStyle_2<-8-DataAnalysis$WorkingStyle_2
DataAnalysis$WorkingStyle_7<-8-DataAnalysis$WorkingStyle_7
DataAnalysis$WorkingStyle_8<-8-DataAnalysis$WorkingStyle_8

RV<-DataAnalysis[,c(32:40)]
RV<-RV[complete.cases(RV),]

ev <- eigen(cor(RV))
ap <- parallel(subject=nrow(RV),
               var=ncol(RV),
               rep=100,cent=.05)
nS <- nScree(x=ev$values, 
             aparallel=ap$eigen$qevpea)
plotnScree(nS)

numFact=nS$Components[[1]]
cormatrix<-cor(RV, method = "pearson")

fit <- fa(r=cormatrix, nfactors=numFact, rotate="varimax")
print(fit, digits=2, cutoff=.3, sort=TRUE)

#Received Opportunity
RV<-DataAnalysis[,c(41:46)]
RV<-RV[complete.cases(RV),]

ev <- eigen(cor(RV))
ap <- parallel(subject=nrow(RV),
               var=ncol(RV),
               rep=100,cent=.05)
nS <- nScree(x=ev$values, 
             aparallel=ap$eigen$qevpea)
plotnScree(nS)

numFact=nS$Components[[1]]
cormatrix<-cor(RV, method = "pearson")

fit <- fa(r=cormatrix, nfactors=numFact, rotate="varimax")
print(fit, digits=2, cutoff=.3, sort=TRUE)

#Diversity Climate
DataAnalysis$DiversityClimate_6<-8-DataAnalysis$DiversityClimate_6

RV<-DataAnalysis[,c(49:55)]
RV<-RV[complete.cases(RV),]

ev <- eigen(cor(RV))
ap <- parallel(subject=nrow(RV),
               var=ncol(RV),
               rep=100,cent=.05)
nS <- nScree(x=ev$values, 
             aparallel=ap$eigen$qevpea)
plotnScree(nS)

numFact=nS$Components[[1]]
cormatrix<-cor(RV, method = "pearson")

fit <- fa(r=cormatrix, nfactors=numFact, rotate="varimax")
print(fit, digits=2, cutoff=.3, sort=TRUE)

#Semantic Differential
DataAnalysis$Affective_1<-8-DataAnalysis$Affective_1
DataAnalysis$Affective_4<-8-DataAnalysis$Affective_4
DataAnalysis$Affective_6<-8-DataAnalysis$Affective_6
DataAnalysis$Affective_7<-8-DataAnalysis$Affective_7
DataAnalysis$Affective_8<-8-DataAnalysis$Affective_8

RV<-DataAnalysis[,c(56:63)]
RV<-RV[complete.cases(RV),]

ev <- eigen(cor(RV))
ap <- parallel(subject=nrow(RV),
               var=ncol(RV),
               rep=100,cent=.05)
nS <- nScree(x=ev$values, 
             aparallel=ap$eigen$qevpea)
plotnScree(nS)

numFact=nS$Components[[1]]
cormatrix<-cor(RV, method = "pearson")

fit <- fa(r=cormatrix, nfactors=numFact, rotate="varimax")
print(fit, digits=2, cutoff=.3, sort=TRUE)

#Workplace Culture
DataAnalysis$WorkplaceCulture_3<-8-DataAnalysis$WorkplaceCulture_3
DataAnalysis$WorkplaceCulture_5<-8-DataAnalysis$WorkplaceCulture_5


RV<-DataAnalysis[,c(64:71)]
RV<-RV[complete.cases(RV),]

ev <- eigen(cor(RV))
ap <- parallel(subject=nrow(RV),
               var=ncol(RV),
               rep=100,cent=.05)
nS <- nScree(x=ev$values, 
             aparallel=ap$eigen$qevpea)
plotnScree(nS)

numFact=nS$Components[[1]]
cormatrix<-cor(RV, method = "pearson")

fit <- fa(r=cormatrix, nfactors=numFact, rotate="varimax")
print(fit, digits=2, cutoff=.3, sort=TRUE)

#Harassment Experience
RV<-DataAnalysis[,c(72:76)]
RV<-RV[complete.cases(RV),]

ev <- eigen(cor(RV))
ap <- parallel(subject=nrow(RV),
               var=ncol(RV),
               rep=100,cent=.05)
nS <- nScree(x=ev$values, 
             aparallel=ap$eigen$qevpea)
plotnScree(nS)

numFact=nS$Components[[1]]
cormatrix<-cor(RV, method = "pearson")

fit <- fa(r=cormatrix, nfactors=numFact, rotate="varimax")
print(fit, digits=2, cutoff=.3, sort=TRUE)


##### Creating Variables and Reliability ####
#Identifying as an Academic
DataAnalysis$IdentifyingAsAnAcademic<-rowMeans(DataAnalysis[,c(10:13)])
psych::alpha(DataAnalysis[,c(10:13)])

#Feeling Authentic
DataAnalysis$FeelingAuthentic<-rowMeans(DataAnalysis[,c(14,16,18:19)])
psych::alpha(DataAnalysis[,c(14,16,18:19)])

#External Validation
DataAnalysis$ExternalValidation<-rowMeans(DataAnalysis[c(15,17)])
psych::alpha(DataAnalysis[,c(15,17)])

#Identifying as an Academic Alternative
DataAnalysis$IdentifyingAsAnAcademicAlt<-rowMeans(DataAnalysis[,c(10:19)])
psych::alpha(DataAnalysis3[,c(10:19)])

#External Validation Alternative
DataAnalysis$ExternalValidationAlt<-rowMeans(DataAnalysis[c(13,15,17)])
psych::alpha(DataAnalysis[,c(13,15,17)])

#Confidence to Succeed
DataAnalysis$ConfidenceToSucceed<-rowMeans(DataAnalysis[c(20:21,24:26)])
psych::alpha(DataAnalysis[c(20:21,24:26)])

#Procedural Justice
DataAnalysis$ProceduralJustice<-rowMeans(DataAnalysis[c(27:29)])
psych::alpha(DataAnalysis[c(27:29)])

#Collaborative Style
DataAnalysis$CollaborativeStyle<-rowMeans(DataAnalysis[c(32:34,40)])
psych::alpha(DataAnalysis[c(32:34,40)])

#Speaking-Up Culture
DataAnalysis$SpeakingUpCulture<-rowMeans(DataAnalysis[c(35:37)])
psych::alpha(DataAnalysis[c(35:37)])

#Received Opportunities
DataAnalysis$ReceivedOpportunity<-rowMeans(DataAnalysis[c(41:46)])
psych::alpha(DataAnalysis[c(41:46)])

#Diversity Climate
DataAnalysis$DiversityClimate<-rowMeans(DataAnalysis[c(49:55)])
psych::alpha(DataAnalysis[,c(49:55)])

#Affective Workplace Appraisal
DataAnalysis$AffectiveWorkplaceAppraisal<-rowMeans(DataAnalysis[c(56:63)])
psych::alpha(DataAnalysis[,c(56:63)])

#Collaborative Culture
DataAnalysis$CollaborativeCulture<-rowMeans(DataAnalysis[c(64:71)])
psych::alpha(DataAnalysis[,c(64:71)])

#Harassment Experience
DataAnalysis$HarassmentExperience<-rowMeans(DataAnalysis[c(72:81)])
psych::alpha(DataAnalysis[,c(72:81)])

#Offline Good
DataAnalysis$Pro_Offline<-rowMeans(DataAnalysis[c(83:90)])
psych::alpha(DataAnalysis[,c(83:90)])

#Online Good
DataAnalysis$Pro_Online<-rowMeans(DataAnalysis[,c(91:98)])
psych::alpha(DataAnalysis[,c(91:98)])

#Diversity and Inclusivity
DataAnalysis$DiversityAndInclusivity<-rowMeans(DataAnalysis[c(49:51,53,55,64:71)])
psych::alpha(DataAnalysis[,c(49:55,64:71)])

#Diversity and Inclusivity Alt
DataAnalysis$PositiveWorkplaceClimate<-rowMeans(DataAnalysis[c(49:55,64:71,56:63,27:29)])
psych::alpha(DataAnalysis[,c(49:55,64:71,56:63,27:29)])


##### Reducing Collinearity by FA #####
RV<-DataAnalysis[,c(49:55,64:71,27:29)]
RV<-RV[complete.cases(RV),]

ev <- eigen(cor(RV))
ap <- parallel(subject=nrow(RV),
               var=ncol(RV),
               rep=100,cent=.05)
nS <- nScree(x=ev$values, 
             aparallel=ap$eigen$qevpea)
plotnScree(nS)

numFact=3
cormatrix<-cor(RV, method = "pearson")

fit <- fa(r=cormatrix, nfactors=3, rotate="varimax")
print(fit, digits=2, cutoff=.3, sort=TRUE)

##### Confirmatory FA #####

Model<-"WorkplaceCulture =~ WorkplaceCulture_1+WorkplaceCulture_2+WorkplaceCulture_3+WorkplaceCulture_4+WorkplaceCulture_5+WorkplaceCulture_6+WorkplaceCulture_7+WorkplaceCulture_8
DiversityClimate =~ DiversityClimate_1+DiversityClimate_2+DiversityClimate_3+DiversityClimate_4+DiversityClimate_5+DiversityClimate_6+DiversityClimate_7
CareerAttitude =~ Career_Attitude_8+Career_Attitude_9+Career_Attitude_10
Climate =~WorkplaceCulture+DiversityClimate+CareerAttitude"


Check<-cfa(Model, data = RV)
summary(Check, fit.measures = TRUE)

##### Analysis Prep ####

#Excluding all Fellows

#Exclude people who want to quit researching all together
DataAnalysis2<-DataAnalysis[which(DataAnalysis$Staying_in_Academia!=3),]

#No under 18s
DataAnalysis2[which(as.numeric(DataAnalysis2$Age)<18),]<-NA

#####Descriptive Statistics ####

# Gender
summary(as.factor(DataAnalysis2$Sex_Gender))

#Age
mean(DataAnalysis2$Age,na.rm=TRUE)
sd(DataAnalysis2$Age,na.rm=TRUE)

#Career Stage
summary(as.factor(DataAnalysis2$Career_Stage))

#Full Time
summary(as.factor(DataAnalysis2$Full_Time))

#Contract
summary(as.factor(DataAnalysis2$Fixed_Term_Contract))

#Ethnicity
round((summary(as.factor(DataAnalysis2$Ethnicity))/nrow(DataAnalysis2))*100,0)

#Citizenship
round((summary(as.factor(DataAnalysis2$UK_Citizen))/nrow(DataAnalysis2))*100,0)

#Sexual Orientation
round((summary(as.factor(DataAnalysis2$Sexual_Orientation))/nrow(DataAnalysis2))*100,0)

#Marital Status
round((summary(as.factor(DataAnalysis2$Marital_Status))/nrow(DataAnalysis2))*100,0)

#Disability 
round((summary(as.factor(DataAnalysis2$Disability))/nrow(DataAnalysis2))*100,0)

#Parental Education
round((summary(as.factor(DataAnalysis2$Parental_Education))/nrow(DataAnalysis2))*100,0)

#Staying in STEM Academia
summary(as.factor(DataAnalysis2$Staying_in_STEM_in_Academia))


#Received Opportunities
psych::alpha(DataAnalysis2[c(41:46)])
round(mean(DataAnalysis2$ReceivedOpportunity,na.rm = TRUE),2)
round(sd(DataAnalysis2$ReceivedOpportunity,na.rm = TRUE),2)

#Harassment Experiences
psych::alpha(DataAnalysis2[,c(72:81)])
round(mean(DataAnalysis2$HarassmentExperience,na.rm = TRUE),2)
round(sd(DataAnalysis2$HarassmentExperience,na.rm = TRUE),2)

# Collaborative Climate 
psych::alpha(DataAnalysis2[,c(64:71)])
round(mean(DataAnalysis2$CollaborativeCulture,na.rm = TRUE),2)
round(sd(DataAnalysis2$CollaborativeCulture,na.rm = TRUE),2)

#Diversity Climate
psych::alpha(DataAnalysis2[,c(49:55)])
round(mean(DataAnalysis2$DiversityClimate,na.rm = TRUE),2)
round(sd(DataAnalysis2$DiversityClimate,na.rm = TRUE),2)

#Speaking Up Climate
psych::alpha(DataAnalysis2[c(35:37)])
round(mean(DataAnalysis2$SpeakingUpCulture,na.rm = TRUE),2)
round(sd(DataAnalysis2$SpeakingUpCulture,na.rm = TRUE),2)
             
#Procedural Justice
psych::alpha(DataAnalysis2[c(27:29)])
round(mean(DataAnalysis2$ProceduralJustice,na.rm = TRUE),2)
round(sd(DataAnalysis2$ProceduralJustice,na.rm = TRUE),2)

#Academic Identity
psych::alpha(DataAnalysis2[,c(10:19)])
round(mean(DataAnalysis2$IdentifyingAsAnAcademicAlt,na.rm = TRUE),2)
round(sd(DataAnalysis2$IdentifyingAsAnAcademicAlt,na.rm = TRUE),2)

#Confidence to Succeed
psych::alpha(DataAnalysis2[c(20:21,24:26)])
round(mean(DataAnalysis2$ConfidenceToSucceed,na.rm = TRUE),2)
round(sd(DataAnalysis2$ConfidenceToSucceed,na.rm = TRUE),2)

#Affective Workplace Evaluation
psych::alpha(DataAnalysis2[,c(56:63)])
round(mean(DataAnalysis2$AffectiveWorkplaceAppraisal,na.rm = TRUE),2)
round(sd(DataAnalysis2$AffectiveWorkplaceAppraisal,na.rm = TRUE),2)

#Collaborative Working style
psych::alpha(DataAnalysis2[c(32:34,40)])
round(mean(DataAnalysis2$CollaborativeStyle,na.rm = TRUE),2)
round(sd(DataAnalysis2$CollaborativeStyle,na.rm = TRUE),2)

#Diversity and Inclusivity
psych::alpha(DataAnalysis2[c(49:51,53,55,64:71)])
round(mean(DataAnalysis2$DiversityAndInclusivity,na.rm = TRUE),2)
round(sd(DataAnalysis2$DiversityAndInclusivity,na.rm = TRUE),2)

#Positive Workplace Climate
psych::alpha(DataAnalysis2[c(49:55,64:71,56:63,27:29)])
round(mean(DataAnalysis2$PositiveWorkplaceClimate,na.rm = TRUE),2)
round(sd(DataAnalysis2$PositiveWorkplaceClimate,na.rm = TRUE),2)

##### Analysis ####

#Standardizing Everything
DataAnalysis2[,c(120:137)]<-scale(DataAnalysis2[,c(120:137)],scale = TRUE,center = TRUE)

#Log Age - including centering of the log
DataAnalysis2$Age[DataAnalysis2$Age<18]<-NA
DataAnalysis2$LogAge<-log(DataAnalysis2$Age)
DataAnalysis2$LogAge<-scale(DataAnalysis2$LogAge,scale = TRUE,center = TRUE)

#Correlations
DataAnalysis2$Staying_in_STEM_in_Academia<-as.numeric(DataAnalysis2$Staying_in_STEM_in_Academia)-1
cor(DataAnalysis2[,c(129,133,132,130,128,126,123,125,131,127,137,117)],use = "pairwise.complete.obs")

#Intraclass Correlation (ICC) via constructing a RI only model.

set.seed(1)
outICC<-matrix(nrow=5000,ncol = 1)

mod0<-glmer(Staying_in_STEM_in_Academia~1+(1|University),
            data = sub,
            family = "binomial",
            na.action = na.omit)

ICCM<-performance::icc(mod0)


DataAnalysis3<-DataAnalysis2
DataAnalysis3[which(as.numeric(DataAnalysis2$Sex_Gender)>2),]<-NA
DataAnalysis3$Sex_Gender<-factor(as.numeric(DataAnalysis3$Sex_Gender),
                                 levels = c(1,2),
                                 labels = c("Cis-Man",
                                            "Cis-Woman"))

##### VIF Check ######
ViFCheck1<-glm(Staying_in_STEM_in_Academia ~
            Career_Stage +
            Fixed_Term_Contract +
            Sex_Gender +
            Sexual_Orientation +
            Parental_Education +
            Ethnicity +
            Disability +
            LogAge +
            HarassmentExperience +
            ReceivedOpportunity +
            ProceduralJustice +
            CollaborativeCulture+
            DiversityClimate+
            AffectiveWorkplaceAppraisal+
            ConfidenceToSucceed+
            IdentifyingAsAnAcademicAlt+
            CollaborativeStyle,
          data=DataAnalysis2,
          family = "binomial")

car::vif(ViFCheck1)

ViFCheck2<-glm(Staying_in_STEM_in_Academia ~
                 Career_Stage +
                 Sex_Gender +
                 Sexual_Orientation +
                 Parental_Education +
                 Ethnicity +
                 Disability +
                 HarassmentExperience +
                 ReceivedOpportunity +
                 PositiveWorkplaceClimate+
                 ProceduralJustice+
                 ConfidenceToSucceed+
                 IdentifyingAsAnAcademicAlt+
                 CollaborativeStyle,
               data=DataAnalysis2,
               family = "binomial")

car::vif(ViFCheck2)

ViFCheck3<-glm(Staying_in_STEM_in_Academia ~
                 Career_Stage +
                 Sex_Gender +
                 Sexual_Orientation +
                 Parental_Education +
                 Ethnicity +
                 Disability +
                 HarassmentExperience +
                 ReceivedOpportunity +
                 PositiveWorkplaceClimate+
                 ConfidenceToSucceed+
                 IdentifyingAsAnAcademicAlt+
                 CollaborativeStyle,
               data=DataAnalysis2,
               family = "binomial")

car::vif(ViFCheck3)

#####model 1 - just demographic variables####

mod1<-glm(Staying_in_STEM_in_Academia~
            Career_Stage+
            Sex_Gender+
            Sexual_Orientation+
            Parental_Education+
            Ethnicity+
            Disability,
          data=DataAnalysis2,
          family = "binomial")
summary(mod1)

#without career stage, gender becomes significant
mod1alt<-glm(Staying_in_STEM_in_Academia~
            Sex_Gender+
            Sexual_Orientation+
            Parental_Education+
            Ethnicity+
            Disability,
          data=DataAnalysis2,
          family = "binomial")
summary(mod1alt)

#Exponents for Odds-Ratios
exp(cbind(OR=coef(mod1),confint(mod1)))

#Exponents for Odds-Ratios
exp(cbind(OR=coef(mod1alt),confint(mod1alt)))

######Model 2 - Experiences #####

mod2<-glm(Staying_in_STEM_in_Academia~
            Career_Stage+
            Sex_Gender+
            Sexual_Orientation+
            Parental_Education+
            Ethnicity+
            Disability+
            ReceivedOpportunity+
            HarassmentExperience,
          data=DataAnalysis2,
          family = "binomial")
summary(mod2)

#Exponents for Odds-Ratios
exp(cbind(OR=coef(mod2),confint(mod2)))

######Model 3 - Workplace Perceptions #####

mod3<-glm(Staying_in_STEM_in_Academia~
            Career_Stage+
            Sex_Gender+
            Sexual_Orientation+
            Parental_Education+
            Ethnicity+
            Disability+
            HarassmentExperience+
            ReceivedOpportunity+
            PositiveWorkplaceClimate,
          data=DataAnalysis2,
          family = "binomial")
summary(mod3)
exp(cbind(OR=coef(mod3),confint(mod3)))

######Model 4 - Self Perceptions #####

mod4<-glm(Staying_in_STEM_in_Academia ~
            Career_Stage +
            Sex_Gender +
            Sexual_Orientation +
            Parental_Education +
            Ethnicity +
            Disability +
            HarassmentExperience +
            ReceivedOpportunity +
            PositiveWorkplaceClimate+
            ConfidenceToSucceed+
            IdentifyingAsAnAcademicAlt+
            CollaborativeStyle,
          data=DataAnalysis2,
          family = "binomial")
summary(mod4)
exp(cbind(OR=coef(mod4),confint(mod4)))

##### Mediation ####

set.seed(1)
out<-matrix(nrow = 5000,ncol = 8)
colnames(out)<-c("Indirect Effect Opportunity", 
                 "Indirect Effect Harassment", 
                 "CPrime",
                 "C","a11","d221","d32","b3")

for (i in 1:5000) 
  {
  sub<-DataAnalysis2[sample(1:nrow(DataAnalysis2), replace = TRUE), ]
  
  cpath<-glm(Staying_in_STEM_in_Academia ~
               Career_Stage +
               Sexual_Orientation +
               Parental_Education +
               Ethnicity +
               Disability +
               ProceduralJustice +
               IdentifyingAsAnAcademicAlt +
               CollaborativeStyle+
               Sex_Gender,
             data=sub,
             family = "binomial")
  
  dofc<-coef(cpath)[13]*(sqrt(3)/pi)
  rofc<-dofc/(sqrt((dofc^2)+4))
  
  a1path<-lm(ReceivedOpportunity ~
               Career_Stage +
               Sexual_Orientation +
               Parental_Education +
               Ethnicity +
               Disability +
               IdentifyingAsAnAcademicAlt +
               CollaborativeStyle+
               Sex_Gender,
             data=sub)
  
  rofa11<-sqrt(eta_squared(a1path)[08,2]) #the sqrt of eta squared is approximate to r.
  if (coef(a1path)[12]<0) {
    rofa11<-rofa11*-1
  }
  
  a2path<-lm(HarassmentExperience ~
               Career_Stage +
               Sexual_Orientation +
               Parental_Education +
               Ethnicity +
               Disability +
               IdentifyingAsAnAcademicAlt +
               CollaborativeStyle+
               Sex_Gender,
             data=sub)
  
  rofa12<-sqrt(eta_squared(a2path)[08,2])
  if (coef(a2path)[12]<0) {
    rofa12<-rofa12*-1
  }
  
  d1path<-lm(PositiveWorkplaceClimate ~
               Career_Stage +
               Sexual_Orientation +
               Parental_Education +
               Ethnicity +
               Disability +
               IdentifyingAsAnAcademicAlt +
               CollaborativeStyle+
               ReceivedOpportunity+
               HarassmentExperience+
               Sex_Gender,
             data=sub)
  
  rofd211<-sqrt(eta_squared(d1path)[8,2]) #the sqrt of eta squared is approximate to r.
  if (coef(d1path)[12]<0) {
    rofd211<-rofd211*-1
  }
  rofd212<-sqrt(eta_squared(d1path)[9,2]) #the sqrt of eta squared is approximate to r.
  if (coef(d1path)[13]<0) {
    rofd212<-rofd212*-1
  }
  rofa2<-sqrt(eta_squared(d1path)[10,2])
  if (coef(d1path)[14]<0) {
    rofa2<-rofa2*-1
  }
  
  d2path<-lm(ConfidenceToSucceed ~
               Career_Stage +
               Sexual_Orientation +
               Parental_Education +
               Ethnicity +
               Disability +
               IdentifyingAsAnAcademicAlt +
               CollaborativeStyle +
               ReceivedOpportunity +
               HarassmentExperience +
               PositiveWorkplaceClimate+
               Sex_Gender,
             data=sub)
  
  rofd32<-sqrt(eta_squared(d2path)[10,2])
  if (coef(d2path)[14]<0) {
    rofd211<-rofd211*-1
  }
  
  rofd311<-sqrt(eta_squared(d2path)[8,2])
  if (coef(d2path)[12]<0) {
    rofd311<-rofd311*-1
  }
  
  rofd312<-sqrt(eta_squared(d2path)[9,2])
  if (coef(d2path)[13]<0) {
    rofd312<-rofd312*-1
  }
  
  rofa3<-sqrt(eta_squared(d2path)[11,2])
  if (coef(d2path)[15]<0) {
    rofa3<-rofa3*-1
  }
  
  bpath<-glm(Staying_in_STEM_in_Academia ~
               Career_Stage +
               Sexual_Orientation +
               Parental_Education +
               Ethnicity +
               Disability +
               IdentifyingAsAnAcademicAlt+
               CollaborativeStyle +
               HarassmentExperience +
               ReceivedOpportunity +
               PositiveWorkplaceClimate+
               ConfidenceToSucceed+
               Sex_Gender,
             data=sub,
             family = "binomial")
  
  dofb3<-coef(bpath)[15]*(sqrt(3)/pi)
  rofb3<-dofb3/(sqrt((dofb3^2)+4))
  
  dofb2<-coef(bpath)[14]*(sqrt(3)/pi)
  rofb2<-dofb2/(sqrt((dofb2^2)+4))
  
  dofb11<-coef(bpath)[13]*(sqrt(3)/pi)
  rofb11<-dofb11/(sqrt((dofb11^2)+4))
  
  dofb12<-coef(bpath)[12]*(sqrt(3)/pi)
  rofb12<-dofb12/(sqrt((dofb12^2)+4))
  
  dofcprime<-coef(bpath)[16]*(sqrt(3)/pi)
  rofcprime<-dofcprime/(sqrt((dofcprime^2)+4))
  

  out[i,1]<-rofa11*rofd211*rofd32*rofb3
  out[i,2]<-rofa12*rofd212*rofd32*rofb3
  out[i,3]<-rofcprime
  out[i,4]<-rofc
  out[i,5]<-coef(a1path)[12]
  out[i,6]<-coef(d1path)[12]
  out[i,7]<-coef(d2path)[14]
  out[i,8]<-coef(bpath)[15]

print(paste(i))
}


##### Exploratory Gender Discrimination Analysis ####

set.seed(1)
out<-matrix(nrow = 5000,ncol = 3)
colnames(out)<-c("Indirect Effect","CPrime","C")

for (i in 1:5000) 
{
  sub<-DataAnalysis2[sample(1:nrow(DataAnalysis2), replace = TRUE), ]
  
  cpath<-glm(Staying_in_STEM_in_Academia ~
               Career_Stage +
               Fixed_Term_Contract +
               Sexual_Orientation +
               Sex_Gender +
               Parental_Education +
               Ethnicity +
               Disability +
               LogAge +
               ReceivedOpportunity +
               SpeakingUpCulture +
               ProceduralJustice +
               IdentifyingAsAnAcademicAlt +
               CollaborativeStyle+
               HarassmentExperience,
             data=sub,
             family = "binomial")
  summary(cpath)
  
  if (is.na(coef(cpath)[19])) {
    
    dofc<-coef(cpath)[18]*(sqrt(3)/pi)
    rofc<-dofc/(sqrt((dofc^2)+4))
    
    apath<-lm(DiversityAndInclusivity ~
                Career_Stage +
                Fixed_Term_Contract +
                LogAge +
                Sexual_Orientation +
                Parental_Education +
                Ethnicity +
                Disability +
                ReceivedOpportunity+
                SpeakingUpCulture +
                ProceduralJustice +
                IdentifyingAsAnAcademicAlt +
                CollaborativeStyle+
                HarassmentExperience*Sex_Gender,
              data=sub)
    summary(apath)
    rofa1<-sqrt(eta_squared(apath)[13,2]) #the sqrt of eta squared is approximate to r.
    if (coef(apath)[17]<0) {
      rofa1<-rofa1*-1
    }
    
    rofa3<-sqrt(eta_squared(apath)[15,2])
    if (coef(apath)[19]<0) {
      rofa3<-rofa3*-1
    }
    dpath<-lm(ConfidenceToSucceed ~
                Career_Stage +
                Fixed_Term_Contract +
                LogAge +
                Sexual_Orientation +
                Parental_Education +
                Ethnicity +
                Disability +
                ReceivedOpportunity+
                SpeakingUpCulture +
                ProceduralJustice +
                IdentifyingAsAnAcademicAlt +
                CollaborativeStyle+
                HarassmentExperience+
                DiversityAndInclusivity*Sex_Gender,
              data=sub)
    summary(dpath)
    rofa2<-sqrt(eta_squared(dpath)[13,2]) #the sqrt of eta squared is approximate to r.
    if (coef(dpath)[17]<0) {
      rofa2<-rofa2*-1
    }
    rofd1<-sqrt(eta_squared(dpath)[14,2]) #the sqrt of eta squared is approximate to r.
    if (coef(dpath)[18]<0) {
      rofd1<-rofd1*-1
    }
    rofd2<-sqrt(eta_squared(dpath)[16,2])
    if (coef(dpath)[20]<0) {
      rofd2<-rofd2*-1
    }
    
    bpath<-glm(Staying_in_STEM_in_Academia ~
                 Career_Stage +
                 Fixed_Term_Contract +
                 Sexual_Orientation +
                 Parental_Education +
                 Ethnicity +
                 Disability +
                 LogAge +
                 ReceivedOpportunity +
                 SpeakingUpCulture +
                 ProceduralJustice +
                 IdentifyingAsAnAcademicAlt +
                 CollaborativeStyle +
                 HarassmentExperience +
                 DiversityAndInclusivity +
                 ConfidenceToSucceed*Sex_Gender,
               data=sub,
               family = "binomial")
    summary(bpath)
    
    dofcprime<-coef(bpath)[17]*(sqrt(3)/pi)
    rofcprime<-dofcprime/(sqrt((dofcprime^2)+4))
    
    dofb2<-coef(bpath)[19]*(sqrt(3)/pi)
    rofb2<-dofb2/(sqrt((dofb2^2)+4))
    
    dofb1<-coef(bpath)[18]*(sqrt(3)/pi)
    rofb1<-dofb1/(sqrt((dofb1^2)+4))
    
    dofb3<-coef(bpath)[21]*(sqrt(3)/pi)
    rofb3<-dofb3/(sqrt((dofb3^2)+4))
    
  } else {
  
    dofc<-coef(cpath)[19]*(sqrt(3)/pi)
    rofc<-dofc/(sqrt((dofc^2)+4))
  
    apath<-lm(DiversityAndInclusivity ~
                Career_Stage +
                Fixed_Term_Contract +
                LogAge +
                Sexual_Orientation +
                Parental_Education +
                Ethnicity +
                Disability +
                ReceivedOpportunity+
                SpeakingUpCulture +
                ProceduralJustice +
                IdentifyingAsAnAcademicAlt +
                CollaborativeStyle+
                HarassmentExperience*Sex_Gender,
              data=sub)
    summary(apath)
    rofa1<-sqrt(eta_squared(apath)[13,2]) #the sqrt of eta squared is approximate to r.
    if (coef(apath)[17]<0) {
      rofa1<-rofa1*-1
    }
    rofa3<-sqrt(eta_squared(apath)[15,2])
    if (coef(apath)[20]<0) {
      rofa3<-rofa3*-1
    }
    dpath<-lm(ConfidenceToSucceed ~
                Career_Stage +
                Fixed_Term_Contract +
                LogAge +
                Sexual_Orientation +
                Parental_Education +
                Ethnicity +
                Disability +
                ReceivedOpportunity+
                SpeakingUpCulture +
                ProceduralJustice +
                IdentifyingAsAnAcademicAlt +
                CollaborativeStyle+
                HarassmentExperience+
                DiversityAndInclusivity*Sex_Gender,
              data=sub)
    summary(dpath)
    rofa2<-sqrt(eta_squared(dpath)[13,2]) #the sqrt of eta squared is approximate to r.
    if (coef(dpath)[17]<0) {
      rofa2<-rofa2*-1
    }
    rofd1<-sqrt(eta_squared(dpath)[14,2]) #the sqrt of eta squared is approximate to r.
    if (coef(dpath)[18]<0) {
      rofd1<-rofd1*-1
    }
    rofd2<-sqrt(eta_squared(dpath)[16,2])
    if (coef(dpath)[21]<0) {
      rofd2<-rofd2*-1
    }
  
    bpath<-glm(Staying_in_STEM_in_Academia ~
                 Career_Stage +
                 Fixed_Term_Contract +
                 Sexual_Orientation +
                 Parental_Education +
                 Ethnicity +
                 Disability +
                 LogAge +
                 ReceivedOpportunity +
                 SpeakingUpCulture +
                 ProceduralJustice +
                 IdentifyingAsAnAcademicAlt +
                 CollaborativeStyle +
                 HarassmentExperience +
                 DiversityAndInclusivity +
                 ConfidenceToSucceed*Sex_Gender,
              data=sub,
              family = "binomial")
    summary(bpath)
  
    dofcprime<-coef(bpath)[17]*(sqrt(3)/pi)
    rofcprime<-dofcprime/(sqrt((dofcprime^2)+4))
  
    dofb2<-coef(bpath)[19]*(sqrt(3)/pi)
    rofb2<-dofb2/(sqrt((dofb2^2)+4))
  
    dofb1<-coef(bpath)[18]*(sqrt(3)/pi)
    rofb1<-dofb1/(sqrt((dofb1^2)+4))
  
    dofb3<-coef(bpath)[22]*(sqrt(3)/pi)
    rofb3<-dofb3/(sqrt((dofb3^2)+4))
  }
  out[i,1]<-rofa11*rofd211*rofd32*rofb3
  out[i,2]<-rofa12*rofd212*rofd32*rofb3
  out[i,2]<-rofcprime
  out[i,3]<-rofc
  
  print(paste(i))
}

##### Additional Simple Slopes Gender Analysis #####

al<-DataAnalysis2[which(DataAnalysis2$Career_Stage=="PhD-Student"),]

modGenderPhD<-glm(Staying_in_STEM_in_Academia
            ~Sex_Gender +
            Fixed_Term_Contract+
            Sexual_Orientation+
            Parental_Education+
            Ethnicity+
            Disability+
            LogAge,
          data=al,
          family = "binomial")
summary(modGenderPhD)

#Exponents for Odds-Ratios
exp(cbind(OR=coef(modGenderPhD),confint(modGenderPhD)))

al<-DataAnalysis2[which(DataAnalysis2$Career_Stage=="Fellow"),]

modGenderPhD2<-glm(Staying_in_STEM_in_Academia ~
                     Career_Stage +
                     Fixed_Term_Contract +
                     Sex_Gender +
                     Sexual_Orientation +
                     Parental_Education +
                     Ethnicity +
                     Disability +
                     LogAge +
                     HarassmentExperience +
                     ReceivedOpportunity +
                     IdentifyingAsAnAcademicAlt,
                   data=DataAnalysis2,
                   family = "binomial")
summary(modGenderPhD2)

#Exponents for Odds-Ratios
exp(cbind(OR=coef(modGenderPhD2),confint(modGenderPhD2)))

apath<-lm( DiversityAndInclusivity~
             Career_Stage+Sex_Gender +
             Fixed_Term_Contract +
             Sexual_Orientation +
             Parental_Education +
             Ethnicity +
             Disability +
             LogAge+
             HarassmentExperience+
             ReceivedOpportunity+
            IdentifyingAsAnAcademicAlt+
             CollaborativeStyle+
             ProceduralJustice+
             SpeakingUpCulture,
           data=DataAnalysis2)
summary(apath)