install.packages("R2admb")
install.packages("glmmADMB", 
                 repos=c("http://glmmadmb.r-forge.r-project.org/repos",
                         getOption("repos")),
                 type="source")


library(lme4)
library(glmmADMB)
library(readr)
library(readxl)




county2 <- read_csv("E:/Rabies/Data/Rabies Project Init - Edited.csv")

# Zone 2 subset

county2Z2 <- subset(county2, RiskZone==2)
County7z2<- as.factor(county2Z2$County)

# Zone 3 subset

county2Z3 <- subset(county2, RiskZone==3)
County7z3<- as.factor(county2Z3$County)

# Zones combined

county2Z23 <- subset(county2, RiskZone>=2)
County7z23<- as.factor(county2Z23$County)

#File for plot

county3 <- read_excel("Rabies/RabiesMelt.xlsx")


#Regression Model for zones combined

mod3a <- glmmadmb(NoRRVpos ~ Year + (1|County7z23) + offset (logPD),
                 data = county2Z23,
                 family= "nbinom",
                 link="log",
                 corStruct = "diag",
                 zeroInflation=FALSE)

summary(mod3a)


#Regression Model for zone 2


mod3b <- glmmadmb(NoRRVpos ~ Year + (1|County7z2) + offset (logPD),
                  data = county2Z2,
                  family= "nbinom",
                  link="log",
                  corStruct = "diag",
                  zeroInflation=FALSE )

summary(mod3b)


##Regression Model for zone 3



mod3c <- glmmadmb(NoRRVpos ~ Year + (1|County7z3) + offset (logPD),
                  data = county2Z3,
                  family= "nbinom",
                  link="log",
                  corStruct = "diag",
                  zeroInflation=FALSE)

summary(mod3c)



#Z score calculation for all hosts

zall= (-0.09 - -0.04)/sqrt(0.018^2 + 0.012^2)
zall
pvalueAll = pnorm(-abs(zall))
pvalueAll

#Breakout for raccoons only

#Regression Model for zones combined

mod4a <- glmmadmb(Racc ~ Year + (1|County7z23) + offset (logPD),
                 data = county2Z23,
                 family= "nbinom",
                 link="log",
                 corStruct = "diag",
                 zeroInflation=TRUE)

summary(mod4a)

#Regression Model for zone 2

mod4b <- glmmadmb(Racc ~ Year + (1|County7z2) + offset (logPD),
                  data = county2Z2,
                  family= "nbinom",
                  link="log",
                  corStruct = "diag",
                  zeroInflation=FALSE)

summary(mod4b)


##Regression Model for zone 3

mod4c <- glmmadmb(Racc ~ Year + (1|County7z3) + offset (logPD),
                  data = county2Z3,
                  family= "nbinom",
                  link="log",
                  corStruct = "diag",
                  zeroInflation=FALSE)

summary(mod4c)

#Z score calculation for raccoons


zrac = (-0.102 - -0.061)/sqrt(0.02^2 + 0.014^2)
zrac
pvaluerac = pnorm(-abs(zrac))
pvaluerac

#Breakout for nondomestic, non-raccoons

#Regression Model for zones combined


mod5a <- glmmadmb(NDNR ~ Year + (1|County7z23) + offset (logPD),
                 data = county2Z23,
                 family= "nbinom",
                 link="log",
                 corStruct = "diag",
                 zeroInflation=TRUE)
summary(mod5a)


#Regression Model for zone 2

mod5b <- glmmadmb(NDNR ~ Year + (1|County7z2) + offset (logPD),
                  data = county2Z2,
                  family= "nbinom",
                  link="log",
                  corStruct = "diag",
                  zeroInflation=FALSE)
summary(mod5b)


##Regression Model for zone 3

mod5c <- glmmadmb(NDNR ~ Year + (1|County7z3) + offset (logPD),
                  data = county2Z3,
                  family= "nbinom",
                  link="log",
                  corStruct = "diag",
                  zeroInflation=FALSE)
summary(mod5c)


#Z score calculation for NDNR

zndnr= (-0.032 - -0.015)/sqrt(0.025^2 + 0.018^2)
zndnr
pvaluendnr = pnorm(-abs(zndnr))
pvaluendnr

#Breakout for domestic animals

#Regression Model for zones combined

mod6a <- glmmadmb(Dom ~ Year + (1|County7z23) + offset (logPD),
                 data = county2Z23,
                 family= "nbinom",
                 link ="log",
                 corStruct = "diag",
                 zeroInflation=FALSE)

summary(mod6a)

#Regression Model for zone 2


mod6b <- glmmadmb(Dom ~ Year + (1|County7z2) + offset (logPD),
                  data = county2Z2,
                  family= "nbinom",
                  link ="log",
                  corStruct = "diag",
                  zeroInflation=FALSE)

summary(mod6b)

##Regression Model for zone 3


mod6c <- glmmadmb(Dom ~ Year + (1|County7z3) + offset (logPD),
                  data = county2Z3,
                  family= "nbinom",
                  link ="log",
                  corStruct = "diag",
                  zeroInflation=FALSE)

summary(mod6c)


#Z score calculation for domestics

zdom = (-0.018 -0.037)/sqrt(0.04^2 + 0.035^2)
zdom
pvaluedom = pnorm(-abs(zdom))
pvaluedom





library(ggplot2)


#Lineplot
ggplot(county3, aes(x=Year, y= Value, group=Type))+
  geom_line(aes(linetype=Type))+ 
  scale_linetype_manual(values=c(4,3,2,1),
                        labels=c("Domestic","NDNR","Raccoon","Total"))+
  labs(y="# RRV positive", linetype = "Animal type")+
  theme_bw()