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

county2 <- read_csv("~/Rabies/Rabies Project Init.csv")
county3 <- read_excel("Rabies/RabiesMelt.xlsx")
summary (county2)

County7<- as.factor(county2$County)

#Linear-mixed model
mod1 <- glmmadmb(NoRRVpos ~ Year + (Year|County7) ,
                  data = county2,
                  family= "gaussian",
                  link="identity")

#Poisson model
mod2 <- glmmadmb(NoRRVpos ~ Year + (Year|County7) + offset(logPD),
                 data = county2,
                 family= "poisson",
                 link="log",
                 corStruct = "diag")

#Zero inflated Poisson model
mod3 <- glmmadmb(NoRRVpos ~ Year + (Year|County7) + offset (logPD),
                 data = county2,
                 family= "poisson",
                 link="log",
                 corStruct = "diag",
                 zeroInflation=TRUE)





#Breakout for raccoons only
mod4 <- glmmadmb(Racc ~ Year + (Year|county7) + offset (logPD),
                 data = county2,
                 family= "poisson",
                 link="log",
                 corStruct = "diag",
                 zeroInflation=TRUE)

#Breakout for nondomestic, non-raccoons
mod5 <- glmmadmb(NDNR ~ Year + (Year|county7) + offset (logPD),
                 data = county2,
                 family= "poisson",
                 link="log",
                 corStruct = "diag",
                 zeroInflation=TRUE)

#Breakout for domestic animals
mod6 <- glmmadmb(Dom ~ Year + (Year|county7) + offset (logPD),
                 data = county2,
                 family= "poisson",
                 link ="log",
                 corStruct = "diag",
                 zeroInflation=TRUE)



summary(mod1)
summary(mod2)
summary(mod3)
summary(mod4)
summary(mod5)
summary(mod6)

#Histogram
ggplot(county2, aes(x=NoRRVpos)) +
  geom_histogram(binwidth = 1) +
  expand_limits(x=50)+
  labs(x="Count of RRV positive animals in a given county and year", y="Frequency of count",
       title="Figure 1: Histogram of total RRV positives by county (n=55) in West Virginia, 2000-2015.")+
  theme_bw()

#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", title="Figure 2: Numbers of RRV positive animals by animal type, West Virginia, 2000-2015.", linetype = "Animal type")+
  theme_bw()