In brouwern/avesdemazan: What the Package Does (One Line, Title Case)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(avesdemazan)

Diet_table <- ecuador_id %>%
  group_by(Band.Number, Diet) %>%
  summarize(Diet_group = length(Diet)) %>%
  group_by(Diet) %>%
  summarize(N = length(Diet))

kable(Diet_table, caption = "Number of individuals by diet")

ecuador$diet_insect <- ifelse(ecuador$Diet == "I", "Insectivore",
                              ifelse(ecuador$Diet == "N", "Nectarivore",
                                     "Omnivore"))

fit_nb2_corr_diet <- glmmTMB(N ~ 1 +                                                 # Null
                            diet_insect +
                            #diet_insect * time_cts +
                            Location +                                               # Location
                            time_cts +                                               # Continuous time variable
                            (1|Specie.Code:Location) +                               # Species-level intercept
                            (time_cts + 0|Specie.Code:Location) +
                            ar1(as.ordered(time_cts) + 0|Specie.Code) +              # Autocorrelation across time
                            offset(log(tot_net_hours)),
                          family = nbinom2,
                          data = ecuador)

Diet_tab_fit <- summary(fit_nb2_corr_diet)$coefficients$cond
rownames(Diet_tab_fit) <- c("Intercept", "Diet (Nectarivore)", "Diet (Omnivore)", "Location (MAIN)", "Location (MASE)", "Time")

#kable(Diet_tab_fit, caption = "Diet model fit on original scale")

# Transformation
vcov <- matrix(unlist(vcov(fit_nb2_corr_diet)), nrow = 6)

# LLAV
Diet_tab_fit2 <- data.frame(Estimate_Llav = c(Diet_tab_fit[1,1],
                                              Diet_tab_fit[1,1] + Diet_tab_fit[2,1],
                                              Diet_tab_fit[1,1] + Diet_tab_fit[3,1]))

Diet_tab_fit2$SE_Llav <- c(round(sqrt(vcov[1,1]), 3),
                           round(sqrt(vcov[1,1] + vcov[2,2] + 2*vcov[1,2]), 3),
                           round(sqrt(vcov[1,1] + vcov[3,3] + 2*vcov[1,3]), 3))

Diet_tab_fit2$LB_Llav <- round(exp(Diet_tab_fit2$Estimate_Llav - 1.96*Diet_tab_fit2$SE_Llav), 5)*1000
Diet_tab_fit2$Est_Llav <- round(exp(Diet_tab_fit2$Estimate_Llav), 5)*1000
Diet_tab_fit2$UB_Llav <- round(exp(Diet_tab_fit2$Estimate_Llav + 1.96*Diet_tab_fit2$SE_Llav), 5)*1000

Diet_tab_fit2$CI_Llav <- paste0("(", Diet_tab_fit2$LB_Llav, ", ", Diet_tab_fit2$Est_Llav, ", ", Diet_tab_fit2$UB_Llav, ")")

# MAIN
Diet_tab_fit2$Estimate_Main <- c(Diet_tab_fit[1,1] + Diet_tab_fit[4,1],
                                 Diet_tab_fit[1,1] + Diet_tab_fit[2,1] + Diet_tab_fit[4,1],
                                 Diet_tab_fit[1,1] + Diet_tab_fit[3,1] + Diet_tab_fit[4,1])

Diet_tab_fit2$SE_Main <- c(round(sqrt(vcov[1,1] + vcov[4,4] + 2*vcov[1,4]), 3),
                           round(sqrt(vcov[1,1] + vcov[2,2] + 2*vcov[1,2] + vcov[4,4] + 2*vcov[1,4] + 2*vcov[2,4]), 3),
                           round(sqrt(vcov[1,1] + vcov[3,3] + 2*vcov[1,3] + vcov[4,4] + 2*vcov[1,4] + 2*vcov[3,4]), 3))

Diet_tab_fit2$LB_Main <- round(exp(Diet_tab_fit2$Estimate_Main - 1.96*Diet_tab_fit2$SE_Main), 5)*1000
Diet_tab_fit2$Est_Main <- round(exp(Diet_tab_fit2$Estimate_Main), 5)*1000
Diet_tab_fit2$UB_Main <- round(exp(Diet_tab_fit2$Estimate_Main + 1.96*Diet_tab_fit2$SE_Main), 5)*1000

Diet_tab_fit2$CI_Main <- paste0("(", Diet_tab_fit2$LB_Main, ", ", Diet_tab_fit2$Est_Main, ", ", Diet_tab_fit2$UB_Main, ")")

# MASE
Diet_tab_fit2$Estimate_Mase <- c(Diet_tab_fit[1,1] + Diet_tab_fit[5,1],
                                 Diet_tab_fit[1,1] + Diet_tab_fit[2,1] + Diet_tab_fit[5,1],
                                 Diet_tab_fit[1,1] + Diet_tab_fit[3,1] + Diet_tab_fit[5,1])

Diet_tab_fit2$SE_Mase <- c(round(sqrt(vcov[1,1] + vcov[5,5] + 2*vcov[1,5]), 3),
                           round(sqrt(vcov[1,1] + vcov[2,2] + 2*vcov[1,2] + vcov[5,5] + 2*vcov[1,5] + 2*vcov[2,5]), 3),
                           round(sqrt(vcov[1,1] + vcov[3,3] + 2*vcov[1,3] + vcov[5,5] + 2*vcov[1,5] + 2*vcov[3,5]), 3))

Diet_tab_fit2$LB_Mase <- round(exp(Diet_tab_fit2$Estimate_Mase - 1.96*Diet_tab_fit2$SE_Mase), 5)*1000
Diet_tab_fit2$Est_Mase <- round(exp(Diet_tab_fit2$Estimate_Mase), 5)*1000
Diet_tab_fit2$UB_Mase <- round(exp(Diet_tab_fit2$Estimate_Mase + 1.96*Diet_tab_fit2$SE_Mase), 5)*1000

Diet_tab_fit2$CI_Mase <- paste0("(", Diet_tab_fit2$LB_Mase, ", ", Diet_tab_fit2$Est_Mase, ", ", Diet_tab_fit2$UB_Mase, ")")

Diet_tab_fit2$Diet <- c("Insectivore", "Nectarivore", "Omnivore")

plot_diet <- data.frame(Estimate = c(Diet_tab_fit2$Est_Llav, Diet_tab_fit2$Est_Main, Diet_tab_fit2$Est_Mase),
                        LB = c(Diet_tab_fit2$LB_Llav, Diet_tab_fit2$LB_Main, Diet_tab_fit2$LB_Mase),
                        UB = c(Diet_tab_fit2$UB_Llav, Diet_tab_fit2$UB_Main, Diet_tab_fit2$UB_Mase), 
                        Diet = rep(Diet_tab_fit2$Diet, 3),
                        Habitat = c(rep("Secondary Forest (LLAV)", 3), 
                                    rep("Introduced Forest (MAIN)", 3),
                                    rep("Primary Forest (MASE)", 3)),
                        Y1 = c(3.15, 3, 2.85,
                              2.15, 2, 1.85,
                              1.15, 1, 0.85),
                        Y2 = c(0.85, 1.85, 2.85,
                               1.15, 2.15, 3.15,
                               1, 2, 3))

plot_diet2 <- plot_diet[,c(5,4,1,2,3)]
colnames(plot_diet2)[4:5] <- c("Lower 95% Estimate", "Upper 95% Estimate")
kable(plot_diet2, caption = "Diet model baseline capture rates per 1000 net hours")

The diet model indicates that there is a significant difference in baseline capture rates between the insectivores and nectarivores. The baseline capture rates for insectivores and omnivores are not statistically different. There is insufficient statistical evidence to suggest that capture rates are changing differently for insectivores, nectarivores, and omnivores (i.e. all diet groups share the estimated 2.83\% decrease in capture rate from year to year; this decrease in capture rate starts from the estimated baseline capture rate). Baseline capture rates per 1000 net hours are summarized in the table above and in the figures on the next page. The time effect is very similar to our original model. We predict a 2.83\% decrease in capture rate from year to year (95\% CI: 0.9-4.7\% decrease).

ggplot(aes(x = Estimate, y = Y1, color = Diet), data = plot_diet) +
  geom_point(size = 1.5) +
  geom_segment(aes(x = LB, xend = UB, y = Y1, yend = Y1), size = 1) +
  theme_bw() +
  scale_y_continuous(breaks = 1:3, labels = c("Primary Forest (MASE)", "Introduced Forest (MAIN)", "Secondary Forest (LLAV)")) +
  ylab("Habitat\n") +
  xlab("Predicted captures per 1000 net hours") +
  ggtitle("Predicted capture rates by diet")

ggplot(aes(x = Estimate, y = Y2, color = Habitat), data = plot_diet) +
  geom_point(size = 1.5) +
  geom_segment(aes(x = LB, xend = UB, y = Y2, yend = Y2), size = 1) +
  theme_bw() +
  scale_y_continuous(breaks = 1:3, labels = c("Insectivore", "Nectarivore", "Omnivore")) +
  ylab("Diet\n") +
  xlab("Predicted captures per 1000 net hours") +
  ggtitle("Predicted capture rates by diet")