library(AICcmodavg) #Used constructing AIC and BIC tables library(rsq) #Import Data Set df<-data.frame(stringsAsFactors=FALSE, Year = c(2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2015L, 2016L, 2014L, 2015L, 2016L, 2014L, 2015L, 2016L, 2014L, 2015L, 2016L, 2014L, 2015L, 2014L, 2015L, 2016L, 2014L, 2015L, 2016L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2014L, 2015L, 2016L, 2014L, 2015L, 2016L, 2014L, 2015L, 2016L, 2014L, 2015L, 2016L, 2015L, 2016L, 2014L, 2015L, 2016L, 2014L, 2015L, 2016L, 2014L, 2015L, 2016L), Count = c(24L, 59L, 79L, 58L, 56L, 21L, 24L, 42L, 33L, 23L, 36L, 5L, 13L, 40L, 61L, 32L, 23L, 11L, 20L, 52L, 63L, 40L, 39L, 18L, 8L, 31L, 42L, 26L, 35L, 20L, 16L, 34L, 66L, 38L, 50L, 20L, 15L, 30L, 33L, 15L, 19L, 12L, 6L, 38L, 36L, 31L, 37L, 19L, 2L, 5L, 8L, 12L, 3L, 12L, 3L, 4L, 3L, 16L, 2L, 21L, 19L, 7L, 16L, 13L, 9L, 1L, 2L, 7L, 3L, 5L, 23L, 6L, 5L, 26L, 16L, 2L, 40L, 81L, 110L, 94L, 107L, 41L, 22L, 79L, 74L, 89L, 47L, 23L, 10L, 53L, 77L, 76L, 30L, 45L, 45L, 83L, 132L, 130L, 82L, 63L, 14L, 59L, 79L, 95L, 79L, 21L, 15L, 87L, 112L, 87L, 90L, 54L, 23L, 57L, 49L, 38L, 51L, 32L, 14L, 90L, 99L, 95L, 58L, 87L, 5L, 15L, 61L, 52L, 18L, 61L, 4L, 10L, 18L, 7L, 22L, 8L, 22L, 40L, 24L, 13L, 12L, 16L, 10L, 2L, 9L, 2L, 12L, 19L, 21L, 7L, 29L, 24L, 11L), TT = c("G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T"), GRID = c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 433L, 433L, 680L, 680L, 680L, 681L, 681L, 681L, 840L, 840L, 840L, 841L, 841L, 842L, 842L, 842L, 843L, 843L, 843L, 850L, 850L, 850L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 433L, 433L, 433L, 680L, 680L, 680L, 681L, 681L, 681L, 840L, 840L, 840L, 841L, 841L, 842L, 842L, 842L, 843L, 843L, 843L, 850L, 850L, 850L) ) #Visualize histograms of the data a<-df[df$TT=="G",] b<-df[df$TT=="T",] hist(a$Count, col='blue', breaks = 10, xlab="# Captured in Ground Traps") hist(b$Count, col='red', add=F, breaks = 10, xlab="# Captured in Tree Traps") #GLM Model fits glm.set<-list() glm.set[[1]]<-glm(Count~ as.factor(TT), data=df, family = poisson()) glm.set[[2]]<-glm(Count~ 1, data=df, family = poisson()) #Model Summaries and Coefficients summary(glm.set[[1]]) coef(glm.set[[1]]) confint(glm.set[[1]]) #AIC table AIC<-aictab(glm.set, sort = TRUE, second.ord = FALSE) #Rsq rsq(glm.set[[1]], type= c('n')) #Nagelkerke, 1991 #Predicted Site- and Year-Specific Capture Frequencies newdat <- data.frame(TT= c("G","T")) ginv <- glm.set[[1]]$family$linkinv ## inverse link function prs <- predict(glm.set[[1]], newdata = newdat, type = "link", se.fit=TRUE) newdat$pred <- ginv(prs[[1]]) newdat$lo <- ginv(prs[[1]] - 1.96 * prs[[2]]) newdat$up <- ginv(prs[[1]] + 1.96 * prs[[2]])