############################################################################
#Professional Boundaries NUS Data
###########################################################################
library(tidyverse)
library(psych)
library(MASS)
library(lmtest)
library(aod)
library(ggpubr)
library(finalfit)
library(naniar)
###########################################################################
#Reading in data and check it
###########################################################################
setwd("C:\\Users\\Alex\\Documents\\My Research\\Professional Boundaries data")

#Read in the data
Survey<-read.csv("Staff-studentSexualMisconductAndAttitudesTowardsProfessionalBoundariesUKHE.csv",header = TRUE,sep="," ,strip.white=TRUE)
#Checking data
str(Survey)
summary(Survey)

###########################################################################
#Creating Factors of Demographic Variables and Tidying Data
###########################################################################
Survey<-Survey%>%mutate(DEM1Gender=factor(DEM1Gender,levels = c(1,2,3,4),labels=c("Male","Female","Non-Binary","Prefer not to Say")),
                        DEM2GenderMatchBirth=factor(DEM2GenderMatchBirth,levels = c(1,2,3),labels = c("Yes","No","Prefer not to say")),
                        DEM3AgeBands=factor(DEM3AgeBands, levels = c(1,2,3,4,5,6,7,8,9), labels = c("Under 18", "18-24", "25-29","30-39","40-49","50-59","60-65","Over 65","Prefer not to say")),
                        DEM4SexualOrientation.=factor(DEM4SexualOrientation., levels=c(1,2,3,4,5,6),labels = c ("Heterosexual/Straight","Gay/Lesbian","Bisexual/Bi","Queer","Prefer not to say","In an other way")),
                        DEM5LearningDisabilityImpairnmentHealthCondition = factor(DEM5LearningDisabilityImpairnmentHealthCondition,levels = c(1,2,3,4),labels = c("No","Yes","Prefer not to say","Don't Know")),
                        DEM6EthnicityReduced=factor(DEM6EthnicityReduced, levels = c(1,2,3,4,5,6), labels = c("White","Black","Asian","Mixed","Other","Prefer not to say")),
                        DEM7Religion=factor(DEM7Religion, levels = c(1,2,3,4,5,6,7,8,9,10,11,12),labels = c("Christian","Buddhist","Hindu","Jewish","Sikh","Muslim","Spirtitual","Agnostic","Atheist","Any other religion","None","Prefer not to say")),
                        DEM8International.status = factor(DEM8International.status,levels = c(1,2,3), labels = c("Home Student","International Student","Prefer not to say")),
                        DEM9CurrentInstitution = factor(DEM9CurrentInstitution, levels = c(1,2,3),labels = c("University","Music Conservatoire","Specialist Arts College")),
                        DEM10ParentCarerGiverAttendedHEI = factor(DEM10ParentCarerGiverAttendedHEI, levels = c(1,2,3),labels = c("Yes","No","Not sure")),
                        DEM11YearStudy = factor(DEM11YearStudy, levels = c(1,2,3,4,5,6,7,8),labels = c("1st UG","2nd UG","3rd UG","4th UG","5th UG","Postgraduate Taught","Postgraduate Research","Other")),
                        DEM11ReducedYearStudy=factor(DEM11ReducedYearStudy, levels = c(1,2,3,4),labels = c("Undergraduate","Postgraduate Taught","Postgraduate Research","Other")),
                        DEM12CurrentPastStudentHEI=factor(DEM12CurrentPastStudentHEI,levels = c(1,2),labels = c("Current Student","Graduated Student")),
                        DEM13PostgraduateStudents= factor(DEM13PostgraduateStudents, levels = c(1,2),labels = c("Postgraduate Taught", "Postgradute Research")))

#Rescale Professional Boundaries. So 0 refers to Very Uncomfortable and 4 Very Comfortable. I.e. lots of 0 equals zero tolerance to professional boundaries blurring.
table(Survey$PBGettingdrunkstudents)
Survey<-Survey %>%mutate_at(vars(7:17),funs(dplyr::recode(.,`5`=0,`4`=1,`3`=2,`2`=3,`1`=4,.default = NaN)))
table(Survey$PBGettingdrunkstudents)

###########################################################################
# Sample Statistics
###########################################################################
#Gender
table(Survey$DEM1Gender)
prop.table(table(Survey$DEM1Gender))

#Gender as assigned at birth
table(Survey$DEM2GenderMatchBirth)
prop.table(table(Survey$DEM2GenderMatchBirth))

#Age
table(Survey$DEM3AgeBands)
prop.table(table(Survey$DEM3AgeBands))

# Year of Study
table(Survey$DEM11ReducedYearStudy)
prop.table(table(Survey$DEM11ReducedYearStudy))

#Sexual Orientation
table(Survey$DEM4SexualOrientation.)
prop.table(table(Survey$DEM4SexualOrientation.))

#Ethnicity/Race
table(Survey$DEM6Ethnicity)
prop.table(table(Survey$DEM6Ethnicity))
table(Survey$DEM6EthnicityReduced)
prop.table(table(Survey$DEM6EthnicityReduced))

#Learning Disability
table(Survey$DEM5LearningDisabilityImpairnmentHealthCondition)
prop.table(table(Survey$DEM5LearningDisabilityImpairnmentHealthCondition))

#Home and International Students
table(Survey$DEM8International.status)
prop.table(table(Survey$DEM8International.status))

###########################################################################
# Reducing data to only variables required for analysis and creating factors.
###########################################################################
PBDF<-Survey%>%dplyr::select("PBGettingdrunkstudents":"PBRomanticRelationshipsStudents","DEM1Gender",
                             "DEM4SexualOrientation.","DEM6EthnicityReduced","DEM8International.status","DEM11ReducedYearStudy")%>%droplevels()

##############################################################################
#Professional Boundaries Factor Analysis
##############################################################################
# How many factors to select
fa.parallel(PBDF[,1:11],fm="pa",fa="fa",
            main="Parallel Analysis Scree Plots",
            n.iter=20) #Suggests 2 clear factors and two marginal other factors.
vss(PBDF[,1:11],title="Very Simple Structure of Professional Boundaries", plot = TRUE)
Scree<-scree(PBDF[,1:11],factors=TRUE,main="Scree plot",hline=NULL,add=FALSE)

#Scree, Velicer MAP and VSS all apear to indicate two factors.
KMO(PBDF[,1:11])#KMO test of sampling adequacy
cor<-corr.test(PBDF[,1:11])
bartlet<-cortest.bartlett(cor$r,n=1492)#Bartlett test of spehericity

#Conducting factor analysis
#Unrotated
fa(PBDF[,1:11],nfactors=2, fm="pa")
#Rotated using promax rotation due to correlation and dropping two items Iniviting to dinner and supervision at house.
PBFactor<-fa(PBDF[,c(1,3,4,5,6,7,8,10,11)],nfactors=2,rotate="promax", scores="tenBerge", fm="pa")
PBFactor
PBFactor$STATISTIC 
PBFactor$dof 
PBFactor$PVAL

#Calculating Cronbach Alpha for whole and subscales.
alpha(PBDF[,1:11])#OVerall Professional Boundaries
alpha(PBDF[,c(5,6,7,10,11)]) #ComfortSexualInteraction .94 alpha
alpha(PBDF[,c(1,3,4,8)]) #ComfortPersonalInteraction .78 alpha

#Sum Scores
PBDF<-PBDF%>% rowwise() %>% mutate(SUMComfortSexualInteraction= sum(c_across(c("PBAskingDate", "PBTellingAttractedYou","PBCommentingBody","PBSexualRelationsStudents", "PBRomanticRelationshipsStudents")), na.rm = FALSE),
                                   SUMComfortPersonalInteraction= sum(c_across(c("PBGettingdrunkstudents","PBAddingFacebook","PBMessageFacebookWhatsapp","PBMeetingsOutsideTimetable")), na.rm = FALSE))
ggplot(PBDF, aes(x=SUMComfortSexualInteraction)) + geom_histogram(binwidth=1, fill="#69b3a2", color="#e9ecef", alpha=0.9)
ggplot(PBDF, aes(x=SUMComfortPersonalInteraction)) + geom_histogram(binwidth=1, fill="#69b3a2", color="#e9ecef", alpha=0.9)

###########################################################################
# Checking Missingness ahead of Demographic Analysis
###########################################################################
CheckMissing<-PBDF%>%dplyr::select("SUMComfortSexualInteraction", "SUMComfortPersonalInteraction","DEM1Gender",
                                   "DEM4SexualOrientation.","DEM6EthnicityReduced","DEM8International.status","DEM11ReducedYearStudy")
CheckMissing%>% ff_glimpse()

CheckMissing %>%missing_plot()

CheckMissing %>% missing_pattern()

#Check missingness for each of the DV and their IV's
mcar_test(CheckMissing[,2:7])#Non significant indicating the data is missing completely at random MCAR
mcar_test(CheckMissing[,c(1,3,4,5,6,7)]) #Non significant indicating the data is missing completely at random MCAR

#Given the small number of missing values and MCAR list wise deletion is used.
PBDF<-na.omit(PBDF)#loses 1492 to 1448 (i.e. 44 (2.94%))
##############################################################################
#Professional Boundaries Demographic Analysis on Un/Comfort with Personal and Sexual interactions
##############################################################################
#Negative binomial for Comfort with Sexual Interaction 
summary(PBComfortSexualInteractionNB<-glm.nb(SUMComfortSexualInteraction ~ DEM1Gender+DEM4SexualOrientation.+DEM6EthnicityReduced+DEM8International.status+DEM11ReducedYearStudy, data=PBDF))
#Incident ratio's for the model
exp(coef(PBComfortSexualInteractionNB))
#Testing main effects of Gender, Home/International and Level of study.
m2 <- update(PBComfortSexualInteractionNB, . ~ . - DEM1Gender)
anova(PBComfortSexualInteractionNB, m2)#Testing gender

m3 <- update(PBComfortSexualInteractionNB, . ~ . - DEM8International.status)
anova(PBComfortSexualInteractionNB, m3)#Home/International

m4 <- update(PBComfortSexualInteractionNB, . ~ . - DEM11ReducedYearStudy)
anova(PBComfortSexualInteractionNB, m4)#level of study

m5 <- update(PBComfortSexualInteractionNB, . ~ . - DEM4SexualOrientation.)
anova(PBComfortSexualInteractionNB, m5)#Sexual orientation

m6 <- update(PBComfortSexualInteractionNB, . ~ . - DEM6EthnicityReduced)
anova(PBComfortSexualInteractionNB, m6)#Sexual orientation

library(lm.beta)
summary(PBSUMComfortPersonalInteraction<-lm(SUMComfortPersonalInteraction~ DEM1Gender+DEM4SexualOrientation.+DEM6EthnicityReduced+DEM8International.status+DEM11ReducedYearStudy, data=PBDF))
lm.beta(PBSUMComfortPersonalInteraction)

# diagnostic plots
layout(matrix(c(1,2,3,4),2,2)) # optional 4 graphs/page
plot(PBSUMComfortPersonalInteraction)

anova(PBSUMComfortPersonalInteraction)#anova table
confint(PBSUMComfortPersonalInteraction, level=0.95) # CIs for model parameters

#Compare poisson model to Negative binomial
summary.glm(SUMComfortSexualInteractionPoi <- glm(SUMComfortSexualInteraction ~ DEM1Gender+DEM4SexualOrientation.+DEM6EthnicityReduced+DEM8International.status+DEM11ReducedYearStudy,family = "poisson", data = PBDF))
lrtest(SUMComfortSexualInteractionPoi, PBComfortSexualInteractionNB)
pchisq(2 * (logLik(PBComfortSexualInteractionNB) - logLik(SUMComfortSexualInteractionPoi)), df = 1, lower.tail = FALSE)
logLik(PBComfortSexualInteractionNB)
logLik(SUMComfortSexualInteractionPoi)

#AIC 10727 Poisson
#AIC 6096 Negative Binomial
#Both the likelhood ratio test and AIC suggest the NB is a substantially better fit than the poisson model.
###########################################################################
# Likert Graph Professional Boundaries
###########################################################################

#library likert package now as it masks recode above which I have dealt with but just to be safe.
library(likert)

items_likert <-PBDF %>%mutate_at(vars(c(1,3,4,5,6,7,8,10,11)),funs(factor(.,levels = c(0,1,2,3,4),labels=c("Very Uncomfortable","Somewhat Uncomfortable","Neutral","Somewhat Comfortable","Very Comfortable"))))
items_likert <-data.frame(items_likert[,c(1,3,4,5,6,7,8,10,11)])
names(items_likert) <- c(
  B1aGettingdrunkstudents="Getting Drunk with students",
  B1cAddingFacebook="Adding you social media",
  B1dMessageFacebookWhatsapp="Sending private messages on social media",
  B1eAskingDate="Asking you out on a date",
  B1fTellingAttractedYou="Telling you they are attracted to you",
  B1gCommentingBody="Commenting on your body",
  B1hMeetingsOutsideTimetable="Arranging meetings outside of academic term time",
  B1jSexualRelationsStudents="Having sexual relations with students",
  B1kRomanticRelationshipsStudents="Having romantic relationships with students")
ItemsLikert<-likert(items_likert)
summary(ItemsLikert)
plot(ItemsLikert,center =3)

rm(items_likert,ItemsLikert)


###########################################################################
#Graphing Professional Boundaries and Demographic Factors
###########################################################################
#Graph Sexualised Interaction Significant Predictors.
#####################################################
#Break down of gender on average response to Sexualised Interactions Subscale
AverageGenderSI<- PBDF%>%pivot_longer(cols=c("PBAskingDate", "PBTellingAttractedYou","PBCommentingBody","PBSexualRelationsStudents", "PBRomanticRelationshipsStudents"), names_to = "RapeAcceptance", values_to = "value")%>%
  group_by(DEM1Gender) %>%summarise(mean.value = mean(value,na.rm = TRUE),median=median(value,na.rm = TRUE), sd.value = sd(value,na.rm = TRUE), count = n(),se.mean = sd.value/sqrt(count))

GenderSI<-ggplot(AverageGenderSI ,aes(x=DEM1Gender, y=mean.value, fill=DEM1Gender))+geom_bar(stat='identity',position ='dodge')+xlab('Gender')+
  ylab('Average Comfort with Sexualised Interactions')+theme_classic()+geom_errorbar(aes(ymin=mean.value-se.mean,ymax=mean.value+se.mean), position= position_dodge(width=0.9),width=.5)+scale_y_continuous(limits = c(0, 4),breaks=seq(0,10,1))+theme(legend.position = "none")

#Break down of gender on average response to Sexualised Interactions Subscale
AverageHomeSI<- PBDF%>%pivot_longer(cols=c("PBAskingDate", "PBTellingAttractedYou","PBCommentingBody","PBSexualRelationsStudents", "PBRomanticRelationshipsStudents"), names_to = "RapeAcceptance", values_to = "value")%>%
  group_by(DEM8International.status) %>%summarise(mean.value = mean(value,na.rm = TRUE),median=median(value,na.rm = TRUE), sd.value = sd(value,na.rm = TRUE), count = n(),se.mean = sd.value/sqrt(count))

HomeSI<-ggplot(AverageHomeSI ,aes(x=DEM8International.status, y=mean.value, fill=DEM8International.status))+geom_bar(stat='identity',position ='dodge')+xlab('Home/International Students')+
  ylab('Average Comfort with Sexualised Interactions')+theme_classic()+geom_errorbar(aes(ymin=mean.value-se.mean,ymax=mean.value+se.mean), position= position_dodge(width=0.9),width=.5)+scale_y_continuous(limits = c(0, 4),breaks=seq(0,10,1))+theme(legend.position = "none")
library(ggpubr)
ggarrange(GenderSI, HomeSI,
          labels = c("A", "B"),
          ncol = 2, nrow = 1)

#####################################################
#Graph Personalised Interaction Significant Predictors.
#####################################################
#Break down of gender on average response to Personalised Interactions Subscale
AverageGenderPI<- PBDF%>%pivot_longer(cols=c("PBGettingdrunkstudents","PBAddingFacebook","PBMessageFacebookWhatsapp","PBMeetingsOutsideTimetable"), names_to = "RapeAcceptance", values_to = "value")%>%
  group_by(DEM1Gender) %>%summarise(mean.value = mean(value,na.rm = TRUE),median=median(value,na.rm = TRUE), sd.value = sd(value,na.rm = TRUE), count = n(),se.mean = sd.value/sqrt(count))

GenderPI<-ggplot(AverageGenderPI ,aes(x=DEM1Gender, y=mean.value, fill=DEM1Gender))+geom_bar(stat='identity',position ='dodge')+xlab('Gender')+
  ylab('Average Comfort with Personalised Interactions')+theme_classic()+geom_errorbar(aes(ymin=mean.value-se.mean,ymax=mean.value+se.mean), position= position_dodge(width=0.9),width=.5)+scale_y_continuous(limits = c(0, 4),breaks=seq(0,10,1))+theme(legend.position = "none")
GenderPI

#Break down of Personalised Orientation on average response to Personalised Interactions Subscale
AverageSexualOPI<- PBDF%>%pivot_longer(cols=c("PBGettingdrunkstudents","PBAddingFacebook","PBMessageFacebookWhatsapp","PBMeetingsOutsideTimetable"), names_to = "RapeAcceptance", values_to = "value")%>%
  group_by(DEM4SexualOrientation.) %>%summarise(mean.value = mean(value,na.rm = TRUE),median=median(value,na.rm = TRUE), sd.value = sd(value,na.rm = TRUE), count = n(),se.mean = sd.value/sqrt(count))

SexualOrientationPI<-ggplot(AverageSexualOPI ,aes(x=DEM4SexualOrientation., y=mean.value, fill=DEM4SexualOrientation.))+geom_bar(stat='identity',position ='dodge')+xlab('Sexual Orientation')+
  ylab('Average Comfort with Personalised Interactions')+theme_classic()+geom_errorbar(aes(ymin=mean.value-se.mean,ymax=mean.value+se.mean), position= position_dodge(width=0.9),width=.5)+scale_y_continuous(limits = c(0, 4),breaks=seq(0,10,1))+theme(legend.position = "none")
SexualOrientationPI

#Break down of Ethnicity on average response to Personalised Interactions Subscale
AverageEthnicityPI<- PBDF%>%pivot_longer(cols=c("PBGettingdrunkstudents","PBAddingFacebook","PBMessageFacebookWhatsapp","PBMeetingsOutsideTimetable"), names_to = "RapeAcceptance", values_to = "value")%>%
  group_by(DEM6EthnicityReduced) %>%summarise(mean.value = mean(value,na.rm = TRUE),median=median(value,na.rm = TRUE), sd.value = sd(value,na.rm = TRUE), count = n(),se.mean = sd.value/sqrt(count))

EthnicityPI<-ggplot(AverageEthnicityPI ,aes(x=DEM6EthnicityReduced, y=mean.value, fill=DEM6EthnicityReduced))+geom_bar(stat='identity',position ='dodge')+xlab('Ethnicity')+
  ylab('Average Comfort with Personalised Interactions')+theme_classic()+geom_errorbar(aes(ymin=mean.value-se.mean,ymax=mean.value+se.mean), position= position_dodge(width=0.9),width=.5)+scale_y_continuous(limits = c(0, 4),breaks=seq(0,10,1))+theme(legend.position = "none")
EthnicityPI

#Break down of home/interational on average response to Sexualised Interactions Subscale
AverageHOMEPI<- PBDF%>%pivot_longer(cols=c("PBGettingdrunkstudents","PBAddingFacebook","PBMessageFacebookWhatsapp","PBMeetingsOutsideTimetable"), names_to = "RapeAcceptance", values_to = "value")%>%
  group_by(DEM8International.status) %>%summarise(mean.value = mean(value,na.rm = TRUE),median=median(value,na.rm = TRUE), sd.value = sd(value,na.rm = TRUE), count = n(),se.mean = sd.value/sqrt(count))

HomeSI<-ggplot(AverageHOMEPI ,aes(x=DEM8International.status, y=mean.value, fill=DEM8International.status))+geom_bar(stat='identity',position ='dodge')+xlab('Home/International Students')+
  ylab('Average Comfort with Personalised Interactions')+theme_classic()+geom_errorbar(aes(ymin=mean.value-se.mean,ymax=mean.value+se.mean), position= position_dodge(width=0.9),width=.5)+scale_y_continuous(limits = c(0, 4),breaks=seq(0,10,1))+theme(legend.position = "none")
HomeSI

ggarrange(GenderPI,SexualOrientationPI,EthnicityPI,HomeSI,
          labels = c("A", "B","C","D"),
          ncol = 2, nrow = 2)