#+ warning=FALSE, message=FALSE
library(mosaic)
library(dplyr)
library(ggplot2)
library(scales)

#' ## HOWL SURVEY DATA 
howl<-read.csv(file="2017howlsurvey_data.csv")
howl$Hsdistance<-as.numeric(howl$Hsdistance)

#' Comparing observer error across experience levels (novice vs. experienced)
obsErr<-favstats(bearDiff~observer,data=howl)
obsErr$skill<-c("e","n","n","n","e","n","n","e") #e=experience,n=novice
nov<-obsErr[obsErr$skill=="n",] #filtering out experienced observer
mean(nov$mean) #mean error for novice observers
sd(nov$mean) #standard deviation for novice observer error
exp<-obsErr[obsErr$skill=="e",]
mean(exp$mean) #mean error for experienced observers
sd(exp$mean) #standard deviation for experienced observer error

#' Independent Samples T-Test comparing experienced vs. novice observers
t.test(mean~skill, data=obsErr)

#'Because there was no difference between experience levels, we combined observers for subsequent analysis
mean(obsErr$mean) #mean observer error
sd(~bearDiff,data=howl) #standard deviation of observer error

#'Percentile Bootstrap for 95% Confidence Interval for Observer Error
bootDist<-do(1000)*{
  bootSamp<-resample(howl$bearDiff)
  mean(bootSamp)
}
quantile(bootDist$result, probs=c(0.025,0.975)) #95% Confidence Interval

#'Quadratic regression to assess the relationship between observer error and distance from homesites
distModel<-lm(bearDiff~Hsdistance+I(Hsdistance^2),data=howl) 
summary(distModel)
par(mfrow=c(2,2))
plot(distModel) #'Checking diagnostics. Meets assumption of regression pretty well.
par(mfrow=c(1,1))

#'Figure showing relationship between observer error and distance from homesites
Fig2<-ggplot(howl, aes(x=Hsdistance, y=bearDiff))+geom_point(size=8)+
  geom_smooth(method="lm", formula=y~x+I(x^2),se=T,size=2)+
  theme_classic()+scale_x_continuous(name="Homesite Distance from Observer (m)")+
  scale_y_continuous(name="Observer Error (degrees)")+
  theme(axis.text = element_text(size=28),axis.title = element_text(size=28))+
  annotate("text", x=1460, y=30, label=paste("R^2 == ", 0.11), size=12, family="Times New Roman", parse=TRUE)+ #manually added text to graph because I couldn't get R2 to display
  annotate("text", x=1450, y=28.5, label="p<0.05", size=12, family="Times New Roman", parse=TRUE) #manually adding the p=value
Fig2

#'For saving figures as png
png(filename="Fig2_HowlSurvey.png", #make the png file 
    type="cairo",
    units="in", 
    width=14, 
    height=14, 
    pointsize=12, 
    res=300) #dpi

Fig2

dev.off()