R/figures/pb-plot-code.R

In davan690/beech-publication-wr: Small mammal dynamics in NZ beech forests

# Prediction A plot
# Oct 2019
# Anthony Davidson
require(extrafont)
# library(jtools)
## ------------------------------------------------------------------------
high.abund.dat <- plot.dat.all1 %>%
  # mutate(N = log(N) %>%  # for when using out.N%>%
  filter(trip.no == 6 | trip.no == 7 | trip.no == 16 | trip.no == 17 )  %>%
  mutate(n.seed.l = N/log.seed,
         n.seed.c = N/cum.seed) %>%
  droplevels()  %>%
  mutate(n.seed.l = ifelse(n.seed.l > 10000 , 0, n.seed.l),
         n.seed.c = ifelse(n.seed.c > 10000 , 0, n.seed.c))
# filter(high.abund.dat, trip.no == 15)$true.date

# summary for plots
# create summarised datasets of need
high.cond.sum <- high.abund.dat %>%
  group_by(Control,Valley,Conditions) %>%
  summarise(N.count = n(),
            N = mean(N, rm.na = TRUE),
            tvalue = qt(p = 0.025, df = N.count - 1, lower.tail = FALSE),
            high.se = mean(se.N),
            lcl.high.tv = N - (tvalue*high.se),
            ucl.high.tv = N + (tvalue*high.se),
            lcl.high = N - (1.96*high.se),
            ucl.high = N + (1.96*high.se))

high.con <- high.abund.dat %>%
  group_by(Control)  %>%
  summarise(N.count = n(),
            N = mean(n.seed.l, rm.na = TRUE),
            tvalue = qt(p = 0.025, df = N.count - 1, lower.tail = FALSE),
            high.se = mean(se.N),
            lcl.high.tv = N - (tvalue*high.se),
            ucl.high.tv = N + (tvalue*high.se),
            lcl.high = N - (1.96*high.se),
            ucl.high = N + (1.96*high.se))

high.con1 <- high.con %>%
  mutate(lcl.high = ifelse(lcl.high < 0, 0, lcl.high),
         Conditions = NA,
         valley = NA)

high.mod3 <- glm(N ~ Control + Valley + Conditions + Conditions:Control, family = "gaussian", data = high.abund.dat)
summary.high.mod3 <- summary(high.mod3)
high.mod3.fun <- anova(high.mod3, test = "F")

# parameters for extraction
output.dens <- plot_summs(high.mod3, scale = TRUE, plot.distributions = TRUE, inner_ci_level = .9)

high.plot.data <- high.abund.dat %>%
  select(n.seed.l, valley, Control, Conditions) %>%
  mutate(Control = factor(Control, labels = c("Yes", "No")),
         Conditions = factor(Conditions, labels = c("No", "Yes")),
         Valley = factor(valley, labels = c("Eglinton", "Hollyford")),
         N = n.seed.l)

plot
high.plot.out <- ggplot(data = high.plot.data, aes(y = N, x = Conditions))+
  # data
  geom_jitter(aes(y = N, x = Conditions, col = Control, shape = Valley, fill = Control),
              stroke = 1.05,
              alpha = 0.7, size = 4, width = 0.3) +

  geom_point(data = high.con1, aes(y = N, x = Control), shape = "square",col= "grey60",  size = 5) +

  geom_errorbar(data = high.con1, aes(x = Control, ymin = lcl.high, ymax = ucl.high,
                                      width = 0),
                lwd = 1,
                 col= "grey60") +

  #labels
  geom_vline(xintercept = 1.5, lty = 2) +
  xlab("Rat removal") +
  ylab(expression(atop(paste("Relative "," ","mouse "," ", "abundance"," "),
                       paste(log(italic(N[jt]))/log(italic(S[jt])))))) +
  #Aesthics
  scale_color_manual(name = "Stoat removal",
                     values = c("black", "black")) +

  scale_shape_manual(name = "Ecosystem",
                     values = c(24, 21)) +

  scale_fill_manual(name = "Rat removal",
                    values = c("white", "black")) +

  guides(fill = "none",
         col = guide_legend(override.aes = list(
           shape = c(0,15), alpha = 1
         )),
         shape = guide_legend(override.aes = list(
           shape = c(24, 21), size = 4
         ))) +

  #theme
  theme_classic()+
  theme_bw() +
  theme(strip.background = element_blank(),
        strip.text.y = element_blank(),

        plot.title = element_text(hjust = 0, size=24, family = "Times", color="black", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "right",
        legend.key = element_blank(),
        legend.background = element_rect(fill="white", size=1),
        legend.key.size=unit(1,"cm"),
        legend.text = element_text(colour = "black", size =10, family = "Times"),
        legend.title = element_text(colour = "black", size =10, family = "Times"),

        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.spacing = unit(2, "lines"),
        panel.border = element_blank(),

        axis.title.y = element_text(colour = "black",size =20, family = "Times", angle = 90),
        axis.title.x = element_text(colour = "black", size =20, family = "Times"),
        axis.text.y=element_text(colour = "black",size = 20, family = "Times"),
        axis.text.x = element_text(colour = "black", size =20, family = "Times"),

        axis.line.y = element_line(size = 1),
        axis.line.x = element_blank(),
        axis.ticks.x = element_blank(),

        strip.text = element_text(face="bold",colour = "black",size =14, family = "Times"))

# high.plot.out