setwd ('C:\\Users\\kabur\\Dropbox\\R workshop\\2021-22 series\\session 4')

#1.	Upload the data "hm_data" 
hm_data <- read.csv ("hm_data.csv", header = T)

#2.Generate a new column in which you calculate grooming rates 
#(by dividing total number of grooming interactions by observation time).
hm_data$gr_rates<- hm_data$Tot_Gr/hm_data$Obs_time

#3.	Check for collinearity 
#between the continuous independent variables. Did you find any collinearity?
cor.test (hm_data$RI, hm_data$HM_agg_rates)
cor.test (hm_data$RI, hm_data$HM_sub_rates)
cor.test (hm_data$RI, hm_data$Prov_rates)
cor.test (hm_data$HM_agg_rates, hm_data$HM_sub_rates)
cor.test (hm_data$HM_agg_rates, hm_data$Prov_rates)
cor.test (hm_data$HM_sub_rates, hm_data$Prov_rates)

#the only variables that have a correlation coefficient above 0.7 are human-macaque
#aggression and submission. So we are going to run multiple regression without human-
#macaque submission

#4.Run linear regression to assess whether human-monkey interactions affect grooming rates, 
#making sure you include variables that are not correlated. 

reg_lm <- lm (gr_rates ~ RI + HM_agg_rates + Prov_rates, data = hm_data)
summary (reg_lm)

#5.Plot the relationship between grooming and rank, with the correct axis labels

ggplot(hm_data, aes(x=RI, y=gr_rates)) +
  geom_point(shape=1) +   
  geom_smooth(method=lm)+
  xlab ("Dominance rank") + ylab ("Grooming rates (sec/hr.)")+ 
  ggtitle ("Relationship between dominance rank and grooming rates")+
  theme(plot.title = element_text(color = "black", size = 20, hjust = .5), 
        axis.text= element_text(color = "black", size = 20, hjust =1, vjust = .5, face = "plain"), 
        axis.title = element_text(color = "black", size = 20, vjust = 2, face = "plain"),
        panel.background = element_rect(fill = "transparent",colour = NA), 
        axis.line = element_line(color="black", size = .2))

#6.	Now run the analysis with Generalized Linear Model 
#with Poisson distribution, rather than linear regression. 

glm_poi <-glm (Tot_Gr~RI + HM_agg_rates + Prov_rates, family=poisson, data=hm_data)
summary (glm_poi)
library(AER)
dispersiontest(glm_poi)


#7.	Run the analysis with Generalized Linear Model with Negative Binomial distribution.
glm_nb <-glm.nb (Tot_Gr~RI + HM_agg_rates + Prov_rates, data=hm_data)
summary (glm_nb)

#8.Plot the relationship between total number of grooming bouts and Provisioning rates
glm_nb <-glm.nb (Tot_Gr~Prov_rates, data=hm_data)
plot(Tot_Gr~Prov_rates, data=hm_data)
curve(predict(glm_nb,data.frame(Prov_rates=x),type="response"),add=TRUE, col="grey50",lty=2)