R/figures/pc-plot-code.R

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

# data for lines ----------------------------------------------------------
# Aug - seed
#Aug filtered (parameters)
# flip dataset
out.para.Aug <- para.plot.dat %>%
  filter(month == "Aug") %>%
  select(mean.b, Control, Valley, para) %>%
  droplevels() %>%
  spread(key = para, value = mean.b)

#this gives us 12 times lines to estimate from


# for each group I need to ....

#Aug filtered (data)
out.dat1.Aug <- out.full.136 %>%
  mutate(Control = factor(control, labels = c("Yes", "No")),
         Valley = factor(valley, labels = c("Eglinton", "Hollyford")),
         Date = as.Date(true.date),
         lag.sjt = lag(log.seed),
         lag.N = lag(N))%>%
  filter(month == "Aug") %>%
  droplevels()

out.dat1.Aug2 <-  out.dat1.Aug %>%
  group_by(Control, Valley) %>%
  summarise(
    M.seed = mean(lag.sjt, na.rm = TRUE),
    M.dens = mean(lag.N, na.rm = TRUE),
    M.rat = mean(lag.rat.mna, na.rm = TRUE),

    MAX.rat = max(lag.rat.mna, na.rm = TRUE),
    MAX.seed = max(lag.sjt, na.rm = TRUE),
    MAX.dens = max(lag.N, na.rm = TRUE),
    MAX.r = max(mean.r, na.rm = TRUE),

    min.seed = 0,
    min.dens = min(lag.N, na.rm = TRUE),
    min.rat = 0,
    min.r = min(mean.r, na.rm = TRUE)) %>%
  ungroup()

# #merge two datasets
Aug.plot <- left_join(out.dat1.Aug2 ,out.para.Aug)

# predict estimates for seed lines
# seed.Aug.plot1 <- Aug.plot %>%
#   mutate(pred.mean = b0 + (b.seed * M.seed) + (b.dens*M.dens) + (b.rat*M.rat),
#          pred.min = b0 + (b.seed * min.seed) + (b.dens*min.dens) + (b.rat*min.rat),
#          pred.max = b0 + (b.seed * MAX.seed) + (b.dens*MAX.dens) + (b.rat*MAX.rat))
seed.Aug.plot1 <- Aug.plot %>%
  mutate(pred.mean = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*M.rat),
         pred.min = Intercept + (Seed * min.seed) + (Density*M.dens) + (Rats*M.dens),
         pred.max = Intercept + (Seed * MAX.seed) + (Density*M.dens) + (Rats*M.dens))

# # r2
# # plot(pred~mean.r, data = seed.Aug.plot1)
#
#
# # data in tidyverse for ggplot2
out.seed.Aug2 <- seed.Aug.plot1 %>%
  gather(key = pred.fact, value = mean.r, pred.min:pred.max) %>%
  mutate(lag.sjt = ifelse(pred.fact == "pred.min", min.seed, MAX.seed)) %>%
  select(Valley, Control, mean.r, lag.sjt) %>%
  mutate(treat = paste(Valley, Control),
         lag.sjt = ifelse(lag.sjt>0, lag.sjt, 0))
#
# # dens-Aug ----------------------------------------------------------------
# #spread
# dens.Aug.plot1 <- Aug.plot %>%
#   mutate(pred.min = b0 + (b.dens * dat.mice.min),
#          pred.max = (b0 + b.dens * dat.mice.max))
#
# #tidyverse
# out.dens.Aug2 <- dens.Aug.plot1 %>%
#   gather(key = pred.fact, value = mean.r, pred.min:pred.max) %>%
#   mutate(lag.sjt = ifelse(pred.fact == "pred.min", dat.mice.min, dat.mice.max)) %>%
#   select(valley, control, mean.r, lag.sjt) %>%
#   mutate(treat = paste(valley, control))
#
#
# # rat-Aug -----------------------------------------------------------------
# #spread and predict
# rat.Aug.plot1 <- Aug.plot %>%
#   mutate(pred.min = b0 + (b.dens * dat.rat.min),
#          pred.max = (b0 + b.dens * dat.rat.max))
#
# #tidyverse
# out.rat.Aug2 <- rat.Aug.plot1 %>%
#   gather(key = pred.fact, value = mean.r, pred.min:pred.max) %>%
#   mutate(lag.sjt = ifelse(pred.fact == "pred.min", dat.rat.min, dat.rat.max)) %>%
#   select(valley, control, mean.r, lag.sjt) %>%
#   mutate(treat = paste(valley, control))
# source("./R/figures/pd-plot-code.R")
# pc.plot

pc.plot <- ggplot(out.dat1.Aug, aes(y = mean.r, x = log.seed)) +
  #error bars
  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r),
                lwd = 0.75,
                alpha = 0.2,
                position=position_dodge(width=30),
                width = 0) +
  #background points
  geom_point(aes(colour = Control,
                 shape = Valley,
                 fill = Conditions),
             stroke = 1.001, size = 2,
             alpha = 0.2) +
  geom_line(data = out.seed.Aug2,
            aes(y = mean.r, x = lag.sjt, group = treat),
            size = 0.75,
            alpha = 0.7) +
geom_point(data = out.seed.Aug2, aes(y = mean.r, x = lag.sjt, shape = Valley, fill = Control, group = treat), size = 5)+
  #Aesthics
  scale_color_manual(name = "Stoat removal",
                     values = c("black", "black", "black", "black")) +

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

  scale_fill_manual(name = "Stoat control",
                    values = c("white", "black","white", "black")) +

  guides(fill = "none",
         col = guide_legend(override.aes = list(
           shape = c(15,0), alpha = 1
         )),
         shape = guide_legend(override.aes = list(
           shape = c(24, 21), alpha = 1, size = 2
         ))) +
  # xlab(expression(paste("Mouse"," ", "Density"," ","(",italic(N[jt]),")" ))) +
  # xlab(expression(atop(paste("Minimum"," ", "number"," "),paste("of", " ", "rats"," ","(",italic(R[jt]),")"))) )+
  # xlab(expression(atop(paste("Minimum"," ", "number"," "),paste("of", " ", "rats"," ","(",italic(R[jt]),")"))) )+
  xlab(expression(paste("Avaliable seed"," ", "(", italic(Seed[jt]), ")"))) +
  ylab(expression(atop(paste("Rate"," ", "of"," ",
                             "increase"),paste(" ", "of"," ",
                                               "mice"," ","(",italic(r[jt]),")"))) ) +
  geom_hline(yintercept = 0, lty = 2, size = 0.5) +
  #theme
  theme_tufte() +
  theme_bw() +
  theme(strip.background = element_blank(),
        strip.text.y = element_blank(),

        plot.title = element_text(hjust = 0, size=24, family = "Times", color="black"),
        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.ticks.x = element_line(size = 1),
        axis.ticks.y = element_line(size = 1),
        axis.line.x = element_line(size = 1),
        axis.line.y = element_line(size = 1),

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

pc.plot
# Save plot
# export plot for example vignette
# png("C://Code/beech-publication-wr/figs/seeded.png")
# pc.plot
# dev.off()