beech-publication-wr: Small mammal dynamics in NZ beech forests

# data
CR.model.out <- readRDS("C://Code/data/final_cauchy_2_5.rds")

## ----import-data-previous-data-------------------------------------------
# Overall Data
# source("./Rcode/manuscript-source-code.R", echo = FALSE)

# data --------------------------------------------------------------------
# export data
plot.dat.all1 <- read_csv("C://Users/s435389/Dropbox/data/plot-all-data1.csv")%>%
          mutate(Valley = factor(valley, labels = c("Eglinton", "Hollyford")))
# out.final <- read_csv("./data/final-outputs.csv")
# abund.dat5 <- read_csv("./data/abundance5.csv")

plot.dat.all1


## ----reduce-previous-data------------------------------------------------

plot.dat.all1 %>%
  select(N, lcl.N, ucl.N, var, grid, Control, Valley, trip, grid.n)


## ----combining-n-rats----------------------------------------------------
# rats have only 80 estimates not 144
meanR <- read_csv("C://Users/s435389/Dropbox/data/mna_allrat.csv") %>%
          mutate(Valley = factor(valley, labels = c("Eglinton", "Hollyford")),
                 trip = trip.no)

meanR <- meanR %>%
  select (Valley, grid ,trip ,n) %>%
  # select (grid ,trip.no ,n) %>%
  mutate(rat.mna = n,
         lag.rat.mna = lag(n))

joined.rats <- full_join(plot.dat.all1, meanR, by = c("Valley", "trip", "grid"))

joined.rats


## ----extracting-rates----------------------------------------------------
# extracting rates of increase
# from scratch
# for now this is hand selected...
rate.dat <- data.frame(CR.model.out$summary[, c(1, 2, 3, 7)])

#renaming variables
names(rate.dat) <- c("mean.r", "se.r", "lcl.r", "ucl.r")
rate.dat$var <- rownames(rate.dat)

#reduce dataset to 136 or 144 rows
## right now it is raw data
rate.dat1 <- rate.dat %>%
  filter(substr(var, 1, 1) == "r")  %>%
  mutate(
    grid = case_when(
      str_detect(var, ",1]") ~ "egl M1",
      str_detect(var, ",2]") ~ "egl M2",
      str_detect(var, ",3]") ~ "egl MR1",
      str_detect(var, ",4]") ~ "egl MR2",
      str_detect(var, ",5]") ~ "hol M1",
      str_detect(var, ",6]") ~ "hol M2",
      str_detect(var, ",7]") ~ "hol MR1",
      str_detect(var, ",8]") ~ "hol MR2",
      TRUE ~ "other"
    ),
    trip = case_when(
      str_detect(string = var, pattern = "r\\[1,") ~ 1,
      str_detect(string = var, pattern = "r\\[2,") ~ 2,
      str_detect(string = var, pattern = "r\\[3,") ~ 3,
      str_detect(string = var, pattern = "r\\[4,") ~ 4,
      str_detect(string = var, pattern = "r\\[5,") ~ 5,
      str_detect(string = var, pattern = "r\\[6,") ~ 6,
      str_detect(string = var, pattern = "r\\[7,") ~ 7,
      str_detect(string = var, pattern = "r\\[8,") ~ 8,
      str_detect(string = var, pattern = "r\\[9,") ~ 9,
      str_detect(string = var, pattern = "r\\[10,") ~ 10,
      str_detect(string = var, pattern = "r\\[11,") ~ 11,
      str_detect(string = var, pattern = "r\\[12,") ~ 12,
      str_detect(string = var, pattern = "r\\[13,") ~ 13,
      str_detect(string = var, pattern = "r\\[14,") ~ 14,
      str_detect(string = var, pattern = "r\\[15,") ~ 15,
      str_detect(string = var, pattern = "r\\[16,") ~ 16,
      str_detect(string = var, pattern = "r\\[17,") ~ 17,
      str_detect(string = var, pattern = "r\\[18,") ~ 18,
      str_detect(string = var, pattern = "r\\[19,") ~ 19,
      str_detect(string = var, pattern = "r\\[20,") ~ 20
    ),
    grid.n = as.numeric(factor(grid)),
    grid = as.character(grid),
    valley = as.character(ifelse(grepl("egl", grid), "egl", "hol"))
  )

rate.dat2 <- rate.dat1 %>%
  select("mean.r",
         "se.r" ,
         "lcl.r" ,
         "ucl.r",
         "trip",
         "grid.n",
         "grid",
         "valley",
         "grid.n")

out.full.136 <- left_join(rate.dat2, joined.rats)

out.full.136$mean.r


## ----parameter-data------------------------------------------------------
# variables ---------------------------------------------------------------
# adding parameter estimates as final data
# add in the appropriate parameter estimates
# parameters
out.para <- data.frame(CR.model.out$summary[, c(1, 2, 3, 7)])
# glimpse(out.para)

names(out.para) <- c("mean.b", "se.b", "lcl.b", "ucl.b")

out.para$var <- rownames(out.para)
var_names <- rownames(out.para)
# glimpse(out.para)

out.para1 <- out.para %>%
  filter(substr(var, 1, 1) == "b") %>%
  mutate(
    Control = case_when(
      str_detect(var, "\\[1,") ~ "control",
      str_detect(var,  "\\[2,") ~ "no control",
      str_detect(var,  "\\[3,") ~ "control",
      TRUE ~ "other"
    ),
    Valley = case_when(
      str_detect(var, "\\[1,") ~ "egl",
      str_detect(var,  "\\[2,") ~ "hol",
      str_detect(var,  "\\[3,") ~ "hol",
      TRUE ~ "other"
    ),
    month = case_when(
      str_detect(var, ",1]") ~ "Feb",
      str_detect(var, ",2]") ~ "May",
      str_detect(var, ",3]") ~ "Aug",
      str_detect(var, ",4]") ~ "Nov",
      TRUE ~ "other"
    )
  )

## issue is with leveling in model and in months in data!!
# great example here!!

    # month = case_when(
    #   str_detect(var, ",1]") ~ "February",
    #   str_detect(var, ",2]") ~ "May",
    #   str_detect(var, ",3]") ~ "August",
    #   str_detect(var, ",4]") ~ "November",
    #   TRUE ~ "other"
    # )

# final parameter dataset -----------------------------------------------------------
# parameter outputs
out.dat1 <- out.full.136 %>%
  group_by(Valley, Control, month) %>%
  summarise(dat.mice.mean = mean(N),
    dat.seed.mean = mean(log.seed),
    dat.rat.mean = mean(lag.rat.mna, na.rm = TRUE),
    dat.mice.min = min(N),
    dat.seed.min = min(log.seed),
    dat.rat.min = min(lag.rat.mna, na.rm = TRUE),
    dat.mice.max = max(N),
    dat.seed.max = max(log.seed),
    dat.rat.max = max(lag.rat.mna)
  ) %>%
  ungroup() %>%
  select(Valley,
         Control,
         month,
         dat.seed.mean,
         dat.rat.mean,
         dat.mice.min,
         dat.seed.min,
         dat.rat.min,
         dat.mice.max,
         dat.seed.max,
         dat.rat.max)%>%
  droplevels()

# out.dat11 <- out.dat1 %>%
#           mutate(Valley = factor(valley, labels = c("Eglinton", "Hollyford")),
#                  trip = trip.no)
# names(plot.dat.all1)


## ------------------------------------------------------------------------
out.para2 <- out.para1[1:48,]

names(out.para2)

out.para3 <- out.para2 %>%
  mutate(Valley = factor(Valley, labels = c("Eglinton", "Hollyford")),
      Control = factor(Control, labels = c("Yes", "No")),
para = case_when(str_detect(var, "b0") ~ "b0",
      str_detect(var, "b.seed") ~ "b.seed",
      str_detect(var, "b.dens") ~ "b.dens",
      str_detect(var, "b.rat") ~ "b.rat"))

# # merge parameter estimates into out.r
out.final <-left_join(out.para3,out.dat1, by = c("Valley",  "Control", "month"))

#remove error estimates for lines
rm.errors <- c("se.b", "lcl.b", "ucl.b", "var")

out.final1 <- out.para3 %>%
                select(-rm.errors)

para.plot.dat <- out.final1 %>%
  mutate_if(is.factor, as.character) %>%
  # filter(month == "Feb") %>%
    droplevels() %>%

                  mutate(Estimate = mean.b)

### With para errors
para.plot.erors <- out.para3 %>%
  mutate_if(is.factor, as.character) %>%
  # filter(month == "Feb") %>%
    droplevels() %>%

                  mutate(Estimate = mean.b)

para.plot.dat$para <- factor(para.plot.dat$para, labels = c("Density", "Rats", "Seed", "Intercept"))

para.plot.dat1 <- para.plot.dat %>%
  mutate_if(is.factor, as.character) %>%
  # filter(month == "Feb") %>%
    droplevels()

glimpse(para.plot.dat1)


## ----include=FALSE-------------------------------------------------------
# creates the reduced parameter dataset for seeds

# 24 long (x2) 12 replicates
# not true
para.plot.dat1 <- para.plot.dat
# currently 48 long
# %>%
                    # filter(month == "May")

#now this in 12 deep
out.para1 <- para.plot.dat1 %>%
  select(mean.b, Control, Valley, para, month) %>%
  spread(key = para, value = mean.b)

# for each group I need to ....
#feb filtered (data)
out.dat1 <- out.full.136 %>%
  mutate(lag.sjt = log.seed,
         lag.N = lag(N)) %>%
  select(-var) %>%
    droplevels()

out.dat2 <-  out.full.136 %>%
      mutate(lag.sjt = log.seed,
         lag.N = lag(N)) %>%
  group_by(Control, Valley, month) %>%
    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
seed.plots <- left_join(out.dat2, out.para1, by = c("Control", "Valley", "month"))

# predict estimates for seed lines
# old parameter names below (dirct from jags)

# seed.plot1 <- seed.plots %>%
#   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.plot1 <- seed.plots %>%
  mutate(pred.mean = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*M.rat),
         pred.min = Intercept + (Seed * min.seed) + (Density*M.dens) + (Rats*M.rat),
         pred.max = Intercept + (Seed * MAX.seed) + (Density*M.dens) + (Rats*M.rat))

# # r2
# # plot(pred~mean.r, data = seed.feb.plot1)
#
#
# # data in tidyverse for ggplot2
out.seed2 <- seed.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, month, pred.fact) %>%
  mutate(lag.sjt = ifelse(lag.sjt>0, lag.sjt, 0.1),
         treat = paste(Valley,Control))

# # dens-feb ----------------------------------------------------------------
# #spread
# dens.feb.plot1 <- feb.plot %>%
#   mutate(pred.min = b0 + (b.dens * dat.mice.min),
#          pred.max = (b0 + b.dens * dat.mice.max))
#
# #tidyverse
# out.dens.feb2 <- dens.feb.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))
#
# feb.plot <- out.dat2 %>%
#   filter(month == "Feb")
# # rat-feb -----------------------------------------------------------------
# #spread and predict
# rat.feb.plot1 <- feb.plot %>%
#   mutate(pred.min = b0 + (b.dens * dat.rat.min),
#          pred.max = (b0 + b.dens * dat.rat.max))
#
# #tidyverse
# out.rat.feb2 <- rat.feb.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))



## ----include=FALSE-------------------------------------------------------
# data ----------------------
### Seed lines (feb)

pC.plot.3 <- out.full.136

out.dat1 <- out.full.136 %>%
  mutate(lag.sjt = log.seed,
         lag.N = lag(N)) %>%
  select(-var) %>%
    droplevels()

# pc.plot.feb <-

ggplot(out.dat1, aes(y = mean.r, x = lag.sjt)) +

  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r), lwd = 0.75, alpha = 0.2, position=position_dodge(width=30), width = 0) +

  geom_point(aes(colour = Valley ,shape = Valley, fill = Control),
             stroke = 1.5, size = 2, alpha = 0.7) +

  # ggplot(out.seed.feb2, aes(y = mean.r, x = lag.sjt)) +
    geom_line(data = out.seed2, aes(y = mean.r, x = lag.sjt, group = treat),size = 0.75, alpha = 0.7) +
  geom_point(data = out.seed2, aes(y = mean.r, x = lag.sjt, shape = Valley, fill = Control, group = treat), size = 5) +

facet_grid(Control~month, scales = "free") +
  scale_shape_manual(name = "Valley",
                         labels = c("E", "H"),
                         values = c(25,21)) +

  scale_colour_manual(name = "Stoat control",
                      labels = c("Eglinton", "Hollyford", "Hollyford"),
                      values = c("darkgoldenrod","black", "black")) +
  scale_fill_manual(name = "Stoat control",
                    labels = c("Yes", "No", "Yes"),
                    values = c("darkgoldenrod","black", "darkgoldenrod")) +
  #
  # # scale_colour_manual(name = "Stoat control",
  # #                     labels = c("E-", "H+", "H-"),
  # #                     values = c("black","black", "black")) +
  # # scale_shape_manual(name = "Valley",
  # #                    labels = c("E", "H"),
  # #                    values = c(25,21)) +
  # # scale_fill_manual(name = "Stoat control",
  # #                   labels = c("E-", "H+", "H-"),
  # #                   values = c("white","black", "white")) +


xlab(expression(paste("Intake"," ", "Rate"," ","(",italic(S[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]),")"))) )+

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
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", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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"))


## ------------------------------------------------------------------------
# data for lines ----------------------------------------------------------
# Feb - seed
#feb filtered (parameters)
# flip dataset
out.para.feb <- para.plot.dat %>%
  filter(month == "Feb") %>%
  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 ....

#feb filtered (data)
out.dat1.feb <- out.dat1 %>%
  filter(month == "Feb") %>%
  mutate(lag.sjt = log.seed,
         lag.N = lag(N)) %>%
  # select(-var) %>%
    droplevels()

out.dat1.feb2 <-  out.full.136 %>%
  filter(month == "Feb") %>%
      mutate(lag.sjt = log.seed,
         lag.N = lag(N)) %>%
  group_by(Control, Valley,month) %>%
    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 = 0,
    min.rat = 0,
    min.r = min(mean.r, na.rm = TRUE)) %>%
  ungroup()

# #merge two datasets
feb.plot <- left_join(out.dat1.feb2, out.para.feb)

feb.plot1 <- feb.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.rat),
         pred.max = Intercept + (Seed * MAX.seed) + (Density*M.dens) + (Rats*M.rat))

# # data in tidyverse for ggplot2
out.seed.feb2 <- feb.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))


## ------------------------------------------------------------------------
# data ----------------------
### Seed lines (feb)

pD.plot.3 <- out.full.136 %>%
  filter(month == "Feb") %>%
  mutate(lag.sjt = log.seed,
         lag.N = lag(N))

pd.plot.feb <- ggplot(pD.plot.3, aes(y = mean.r, x = lag.sjt)) +

  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r), lwd = 0.75, alpha = 0.2, position=position_dodge(width=30), width = 0) +

  geom_point(aes(colour = Valley,shape = Valley, fill = Control),
             stroke = 1.5, size = 2, alpha = 0.7) +

  # ggplot(out.seed.feb2, aes(y = mean.r, x = lag.sjt)) +
    geom_line(data = out.seed.feb2, aes(y = mean.r, x = lag.sjt, group = treat),size = 0.75, alpha = 0.7) +
  geom_point(data = out.seed.feb2, aes(y = mean.r, x = lag.sjt, shape = Valley, fill = Control, group = treat), size = 5) +

   scale_color_manual(name = "Stoat Control",
                     values = c("white", "black", "white")) +

  scale_shape_manual(name = "Ecosystem",
                     values = c(24, 21)) +
  # scale_size_manual(name = "Rat Control", values = c(2.5, 3, 2.5)) +
  scale_fill_manual(
    name = "Stoat Control",
    values = c("cornflowerblue", "darkorange", "cornflowerblue")
  ) +
     # Remove fill legend and replace the fill legend using the newly created size
  guides(
    col = "none",
    size = guide_legend(override.aes = list(
  shape = c(15,0),alpha = 1
    )),
    shape = guide_legend(override.aes = list(
      shape = c(24, 21), size = 4
    )),
    fill = guide_legend(override.aes = list(
      col = c("cornflowerblue", "darkorange"),shape = c("square"),
      size = 4
    ))
  ) +

  #
  # # scale_colour_manual(name = "Stoat control",
  # #                     labels = c("E-", "H+", "H-"),
  # #                     values = c("black","black", "black")) +
  # # scale_shape_manual(name = "Valley",
  # #                    labels = c("E", "H"),
  # #                    values = c(25,21)) +
  # # scale_fill_manual(name = "Stoat control",
  # #                   labels = c("E-", "H+", "H-"),
  # #                   values = c("white","black", "white")) +


xlab(expression(paste("Intake"," ", "Rate"," ","(",italic(S[jt]),")" ))) +
# 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]),")"))) )+

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
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", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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"))
pd.plot.feb


## ------------------------------------------------------------------------
# data for lines ----------------------------------------------------------
# May - seed
#May filtered (parameters)
# flip dataset
out.para.May <- para.plot.dat %>%
  filter(month == "May") %>%
  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 ....

#May filtered (data)
out.dat1.May <- out.dat1 %>%
  filter(month == "May") %>%
  mutate(lag.sjt = log.seed,
         lag.N = lag(N)) %>%
  # select(-var) %>%
    droplevels()

out.dat1.May2 <-  out.full.136 %>%
  filter(month == "May") %>%
      mutate(lag.sjt = log.seed,
         lag.N = lag(N)) %>%
  group_by(Control, Valley, month) %>%
    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 = 0,
    min.rat = 0,
    min.r = min(mean.r, na.rm = TRUE)) %>%
  ungroup()

# #merge two datasets
May.plot <- left_join(out.dat1.May2, out.para.May)

May.plot1 <- May.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.rat),
         pred.max = Intercept + (Seed * MAX.seed) + (Density*M.dens) + (Rats*M.rat))

# # data in tidyverse for ggplot2
out.seed.May2 <- May.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))


## ------------------------------------------------------------------------
# data ----------------------
### Seed lines (May)

pD.plot.3 <- out.full.136 %>%
  filter(month == "May") %>%
  mutate(lag.sjt = log.seed,
         lag.N = lag(N))

pd.plot.May <- ggplot(pD.plot.3, aes(y = mean.r, x = lag.sjt)) +

  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r), lwd = 0.75, alpha = 0.2, position=position_dodge(width=30), width = 0) +

  geom_point(aes(colour = Valley,shape = Valley, fill = Control),
             stroke = 1.5, size = 2, alpha = 0.7) +

  # ggplot(out.seed.May2, aes(y = mean.r, x = lag.sjt)) +
    geom_line(data = out.seed.May2, aes(y = mean.r, x = lag.sjt, group = treat),size = 0.75, alpha = 0.7) +
  geom_point(data = out.seed.May2, aes(y = mean.r, x = lag.sjt, shape = Valley, fill = Control, group = treat), size = 5) +

   scale_color_manual(name = "Stoat Control",
                     values = c("white", "black", "white")) +

  scale_shape_manual(name = "Ecosystem",
                     values = c(24, 21)) +
  # scale_size_manual(name = "Rat Control", values = c(2.5, 3, 2.5)) +
  scale_fill_manual(
    name = "Stoat Control",
    values = c("cornflowerblue", "darkorange", "cornflowerblue")
  ) +
     # Remove fill legend and replace the fill legend using the newly created size
  guides(
    col = "none",
    size = guide_legend(override.aes = list(
  shape = c(15,0),alpha = 1
    )),
    shape = guide_legend(override.aes = list(
      shape = c(24, 21), size = 4
    )),
    fill = guide_legend(override.aes = list(
      col = c("cornflowerblue", "darkorange"),shape = c("square"),
      size = 4
    ))
  ) +

  #
  # # scale_colour_manual(name = "Stoat control",
  # #                     labels = c("E-", "H+", "H-"),
  # #                     values = c("black","black", "black")) +
  # # scale_shape_manual(name = "Valley",
  # #                    labels = c("E", "H"),
  # #                    values = c(25,21)) +
  # # scale_fill_manual(name = "Stoat control",
  # #                   labels = c("E-", "H+", "H-"),
  # #                   values = c("white","black", "white")) +


xlab(expression(paste("Intake"," ", "Rate"," ","(",italic(S[jt]),")" ))) +
# 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]),")"))) )+

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
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", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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"))
pd.plot.May


## ------------------------------------------------------------------------
# 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.dat1 %>%
  filter(month == "Aug") %>%
  mutate(lag.sjt = log.seed,
         lag.N = N) %>%
  # select(-var) %>%
    droplevels()

out.dat1.Aug2 <-  out.full.136 %>%
  filter(month == "Aug") %>%
      mutate(lag.sjt = log.seed,
         lag.N = N) %>%
  group_by(Control, Valley, month) %>%
    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 = 0,
    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)

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.rat),
         pred.max = Intercept + (Seed * MAX.seed) + (Density*M.dens) + (Rats*M.rat))

# # data in tidyverse for ggplot2
out.dens.Aug2 <- 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))


## ------------------------------------------------------------------------
# data ----------------------
### Seed lines (Aug)

pD.plot.3 <- out.full.136 %>%
  filter(month == "Aug") %>%
  mutate(lag.sjt = log.seed,
         lag.N = N)

pd.plot.Aug <- ggplot(pD.plot.3, aes(y = mean.r, x = lag.sjt)) +

  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r), lwd = 0.75, alpha = 0.2, position=position_dodge(width=30), width = 0) +

  geom_point(aes(colour = Control,shape = Valley, fill = Control),
             stroke = 1.5, size = 2, alpha = 0.7) +

  # ggplot(out.seed.Aug2, aes(y = mean.r, x = lag.sjt)) +
    geom_line(data = out.dens.Aug2, aes(y = mean.r, x = lag.sjt, group = treat),size = 0.75, alpha = 0.7) +
  geom_point(data = out.dens.Aug2, aes(y = mean.r, x = lag.sjt, shape = Valley, fill = Control, group = treat), size = 5) +

   scale_color_manual(name = "Stoat Control",
                     values = c("white", "black", "white")) +

  scale_shape_manual(name = "Ecosystem",
                     values = c(24, 21)) +
  # scale_size_manual(name = "Rat Control", values = c(2.5, 3, 2.5)) +
  scale_fill_manual(
    name = "Stoat Control",
    values = c("cornflowerblue", "darkorange", "cornflowerblue")
  ) +
     # Remove fill legend and replace the fill legend using the newly created size
  guides(
    col = "none",
    size = guide_legend(override.aes = list(
  shape = c(15,0),alpha = 1
    )),
    shape = guide_legend(override.aes = list(
      shape = c(24, 21), size = 4
    )),
    fill = guide_legend(override.aes = list(
      col = c("cornflowerblue", "darkorange"),shape = c("square"),
      size = 4
    ))
  ) +

  #
  # # scale_colour_manual(name = "Stoat control",
  # #                     labels = c("E-", "H+", "H-"),
  # #                     values = c("black","black", "black")) +
  # # scale_shape_manual(name = "Valley",
  # #                    labels = c("E", "H"),
  # #                    values = c(25,21)) +
  # # scale_fill_manual(name = "Stoat control",
  # #                   labels = c("E-", "H+", "H-"),
  # #                   values = c("white","black", "white")) +


xlab(expression(paste("Intake"," ", "Rate"," ","(",italic(S[jt]),")" ))) +
# 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]),")"))) )+

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
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", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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"))
pd.plot.Aug


## ------------------------------------------------------------------------
# data for lines ----------------------------------------------------------
# Nov - seed
#Nov filtered (parameters)
# flip dataset
out.para.Nov <- para.plot.dat %>%
  filter(month == "Nov") %>%
  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 ....

#Nov filtered (data)
out.dat1.Nov <- out.dat1 %>%
  filter(month == "Nov") %>%
  mutate(lag.sjt = log.seed,
         lag.N = N) %>%
  # select(-var) %>%
    droplevels()

out.dat1.Nov2 <-  out.full.136 %>%
  filter(month == "Nov") %>%
      mutate(lag.sjt = log.seed,
         lag.N = N) %>%
  group_by(Control, Valley, month) %>%
    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 = 0,
    min.rat = 0,
    min.r = min(mean.r, na.rm = TRUE)) %>%
  ungroup()

# #merge two datasets
Nov.plot <- left_join(out.dat1.Nov2, out.para.Nov)

Nov.plot1 <- Nov.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.rat),
         pred.max = Intercept + (Seed * MAX.seed) + (Density*M.dens) + (Rats*M.rat))

# # data in tidyverse for ggplot2
out.dens.Nov2 <- Nov.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))


## ------------------------------------------------------------------------
# data ----------------------
### Seed lines (Nov)

pD.plot.3 <- out.full.136 %>%
  filter(month == "Nov") %>%
  mutate(lag.sjt = log.seed,
         lag.N = N)

pd.plot.Nov <- ggplot(pD.plot.3, aes(y = mean.r, x = lag.sjt)) +

  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r), lwd = 0.75, alpha = 0.2, position=position_dodge(width=30), width = 0) +

  geom_point(aes(colour = Valley,shape = Valley, fill = Control),
             stroke = 1.5, size = 2, alpha = 0.7) +

  # ggplot(out.seed.Nov2, aes(y = mean.r, x = lag.sjt)) +
    geom_line(data = out.dens.Nov2, aes(y = mean.r, x = lag.sjt, group = treat),size = 0.75, alpha = 0.7) +
  geom_point(data = out.dens.Nov2, aes(y = mean.r, x = lag.sjt, shape = Valley, fill = Control, group = treat), size = 5) +

   scale_color_manual(name = "Stoat Control",
                     values = c("white", "black", "white")) +

  scale_shape_manual(name = "Ecosystem",
                     values = c(24, 21)) +
  # scale_size_manual(name = "Rat Control", values = c(2.5, 3, 2.5)) +
  scale_fill_manual(
    name = "Stoat Control",
    values = c("cornflowerblue", "darkorange", "cornflowerblue")
  ) +
     # Remove fill legend and replace the fill legend using the newly created size
  guides(
    col = "none",
    size = guide_legend(override.aes = list(
  shape = c(15,0),alpha = 1
    )),
    shape = guide_legend(override.aes = list(
      shape = c(24, 21), size = 4
    )),
    fill = guide_legend(override.aes = list(
      col = c("cornflowerblue", "darkorange"),shape = c("square"),
      size = 4
    ))
  ) +

  #
  # # scale_colour_manual(name = "Stoat control",
  # #                     labels = c("E-", "H+", "H-"),
  # #                     values = c("black","black", "black")) +
  # # scale_shape_manual(name = "Valley",
  # #                    labels = c("E", "H"),
  # #                    values = c(25,21)) +
  # # scale_fill_manual(name = "Stoat control",
  # #                   labels = c("E-", "H+", "H-"),
  # #                   values = c("white","black", "white")) +


xlab(expression(paste("Intake"," ", "Rate"," ","(",italic(S[jt]),")" ))) +
# 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]),")"))) )+

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
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", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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"))
pd.plot.Nov


## ----include=FALSE-------------------------------------------------------
# data for lines ----------------------------------------------------------
# Feb - seed
#feb filtered (parameters)
# flip dataset
para.plot.dat1 <- para.plot.dat

out.para1 <- para.plot.dat %>%
  select(mean.b, Control, Valley, para, month) %>%
  # filter(month == "Feb") %>%
    droplevels() %>%
  spread(key = para, value = mean.b)

# for each group I need to ....

#feb filtered (data)
out.dat1 <- out.full.136 %>%
  filter(month == "Feb") %>%
  select(-var) %>%
    droplevels()

out.dat2 <-  out.full.136 %>%
  group_by(Control, Valley, month) %>%
    summarise(
    M.seed = mean(log.seed, na.rm = TRUE),
    M.dens = mean(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(log.seed, na.rm = TRUE),
    MAX.dens = max(N, na.rm = TRUE),
    MAX.r = max(mean.r, na.rm = TRUE),

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

# #merge two datasets
dens.plots <- bind_cols(out.dat2, out.para1)

# predict estimates for seed lines
# seed.plot1 <- seed.plots %>%
#   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))
dens.plot1 <- dens.plots %>%
  mutate(pred.mean = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*M.rat),
         pred.min = Intercept + (Seed * M.seed) + (Density*min.dens) + (Rats*M.rat),
         pred.max = Intercept + (Seed * M.seed) + (Density*MAX.dens) + (Rats*M.rat))

# # r2
# # plot(pred~mean.r, data = seed.feb.plot1)
#
#
# # data in tidyverse for ggplot2
out.dens2 <- dens.plot1 %>%
  gather(key = pred.fact, value = mean.r, pred.min:pred.max) %>%
  mutate(lag.N = ifelse(pred.fact == "pred.min", min.dens, MAX.dens)) %>%
  select(Valley, Control, mean.r, lag.N, month, pred.fact) %>%
  mutate(lag.N = ifelse(lag.N>0, lag.N, 0.1),
         treat = paste(Valley,Control))

# dens-feb ----------------------------------------------------------------
#spread
#NOT WORKING
# dens.feb.plot1 <- out.dat1 %>%
#   mutate(pred.min = b0 + (b.dens * dat.mice.min),
#          pred.max = (b0 + b.dens * dat.mice.max))
#
# #tidyverse
# out.dens.feb2 <- dens.feb.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-feb -----------------------------------------------------------------
# #spread and predict
# rat.feb.plot1 <- feb.plot %>%
#   mutate(pred.min = b0 + (b.dens * dat.rat.min),
#          pred.max = (b0 + b.dens * dat.rat.max))
#
# #tidyverse
# out.rat.feb2 <- rat.feb.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))



## ----include=FALSE-------------------------------------------------------
# data ----------------------

pC.plot.3 <- out.full.136

out.dat1 <- out.full.136 %>%
  mutate(lag.sjt = log.seed,
         lag.N = lag(N)) %>%
  select(-var) %>%
    droplevels()

ggplot(out.dat1, aes(y = mean.r, x = lag.N)) +

  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r), lwd = 0.75, alpha = 0.2, position=position_dodge(width=30), width = 0) +

  geom_point(aes(colour = Valley,shape = Valley, fill = Control),
             stroke = 1.5, size = 2, alpha = 0.7) +

  # ggplot(out.seed.feb2, aes(y = mean.r, x = lag.sjt)) +
    geom_line(data = out.dens2, aes(y = mean.r, x = lag.N, group = treat),size = 0.75, alpha = 0.7) +
  geom_point(data = out.dens2, aes(y = mean.r, x = lag.N, shape = Valley, fill = Control, group = treat), size = 5) +

facet_grid(Control~month, scales = "free") +
  scale_shape_manual(name = "Valley",
                         labels = c("E", "H"),
                         values = c(25,21)) +

  scale_colour_manual(name = "Stoat control",
                      labels = c("Eglinton", "Hollyford", "Hollyford"),
                      values = c("darkgoldenrod","black", "black")) +
  scale_fill_manual(name = "Stoat control",
                    labels = c("Yes", "No", "Yes"),
                    values = c("darkgoldenrod","black", "darkgoldenrod")) +
  #
  # # scale_colour_manual(name = "Stoat control",
  # #                     labels = c("E-", "H+", "H-"),
  # #                     values = c("black","black", "black")) +
  # # scale_shape_manual(name = "Valley",
  # #                    labels = c("E", "H"),
  # #                    values = c(25,21)) +
  # # scale_fill_manual(name = "Stoat control",
  # #                   labels = c("E-", "H+", "H-"),
  # #                   values = c("white","black", "white")) +


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]),")"))) )+

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
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", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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"))



## ------------------------------------------------------------------------
# data for lines ----------------------------------------------------------
# Feb - seed
#feb filtered (parameters)
# flip dataset
out.para.feb <- para.plot.dat %>%
  filter(month == "Feb") %>%
  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 ....

#feb filtered (data)
out.dat1.feb <- out.dat1 %>%
  filter(month == "Feb") %>%
  mutate(lag.sjt = log.seed,
         lag.N = lag(N)) %>%
  # select(-var) %>%
    droplevels()

out.dat1.feb2 <-  out.full.136 %>%
  filter(month == "Feb") %>%
      mutate(lag.sjt = log.seed,
         lag.N = lag(N)) %>%
  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 = 0,
    min.rat = 0,
    min.r = min(mean.r, na.rm = TRUE)) %>%
  ungroup()

# #merge two datasets
feb.plot <- left_join(out.dat1.feb2, out.para.feb)

feb.plot1 <- feb.plot %>%
  mutate(pred.mean = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*M.rat),
         pred.min = Intercept + (Seed * M.seed) + (Density*min.dens) + (Rats*M.dens),
         pred.max = Intercept + (Seed * M.seed) + (Density*MAX.dens) + (Rats*M.dens))

# # data in tidyverse for ggplot2
out.dens.feb2 <- feb.plot1 %>%
  gather(key = pred.fact, value = mean.r, pred.min:pred.max) %>%
  mutate(lag.N = ifelse(pred.fact == "pred.min", min.dens, MAX.dens)) %>%
  select(Valley, Control, mean.r, lag.N) %>%
  mutate(treat = paste(Valley, Control),
         lag.N = ifelse(lag.N>0, lag.N, 0.1))


## ------------------------------------------------------------------------
# data ----------------------
### Seed lines (feb)

pD.plot.3 <- out.full.136 %>%
  filter(month == "Feb") %>%
  mutate(lag.sjt = log.seed,
         lag.N = lag(N))

pd.plot.feb <- ggplot(pD.plot.3, aes(y = mean.r, x = lag.N)) +

  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r), lwd = 0.75, alpha = 0.2, position=position_dodge(width=30), width = 0) +

  geom_point(aes(colour = Valley,shape = Valley, fill = Control),
             stroke = 1.5, size = 2, alpha = 0.7) +

  # ggplot(out.seed.feb2, aes(y = mean.r, x = lag.sjt)) +
    geom_line(data = out.dens.feb2, aes(y = mean.r, x = lag.N, group = treat),size = 0.75, alpha = 0.7) +
  geom_point(data = out.dens.feb2, aes(y = mean.r, x = lag.N, shape = Valley, fill = Control, group = treat), size = 5) +

   scale_color_manual(name = "Stoat Control",
                     values = c("white", "black", "white")) +

  scale_shape_manual(name = "Ecosystem",
                     values = c(24, 21)) +
  # scale_size_manual(name = "Rat Control", values = c(2.5, 3, 2.5)) +
  scale_fill_manual(
    name = "Stoat Control",
    values = c("cornflowerblue", "darkorange", "cornflowerblue")
  ) +
     # Remove fill legend and replace the fill legend using the newly created size
  guides(
    col = "none",
    size = guide_legend(override.aes = list(
  shape = c(15,0),alpha = 1
    )),
    shape = guide_legend(override.aes = list(
      shape = c(24, 21), size = 4
    )),
    fill = guide_legend(override.aes = list(
      col = c("cornflowerblue", "darkorange"),shape = c("square"),
      size = 4
    ))
  ) +

  #
  # # scale_colour_manual(name = "Stoat control",
  # #                     labels = c("E-", "H+", "H-"),
  # #                     values = c("black","black", "black")) +
  # # scale_shape_manual(name = "Valley",
  # #                    labels = c("E", "H"),
  # #                    values = c(25,21)) +
  # # scale_fill_manual(name = "Stoat control",
  # #                   labels = c("E-", "H+", "H-"),
  # #                   values = c("white","black", "white")) +


# xlab(expression(paste("Intake"," ", "Rate"," ","(",italic(S[jt]),")" ))) +
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]),")"))) )+

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
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", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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"))
pd.plot.feb


## ------------------------------------------------------------------------
# data for lines ----------------------------------------------------------
# May - seed
#May filtered (parameters)
# flip dataset
out.para.May <- para.plot.dat %>%
  filter(month == "May") %>%
  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 ....

#May filtered (data)
out.dat1.May <- out.full.136 %>%
  filter(month == "May") %>%
  mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N)) %>%
  select(-var) %>%
    droplevels()

out.dat1.May2 <-  out.full.136 %>%
  filter(month == "May") %>%
      mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N)) %>%
  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 = 0,
    min.rat = 0,
    min.r = min(mean.r, na.rm = TRUE)) %>%
  ungroup()

# #merge two datasets
May.plot <- left_join(out.dat1.May2 ,out.para.May)

# predict estimates for seed lines
# seed.May.plot1 <- May.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))
May.plot1 <- May.plot %>%
  mutate(pred.mean = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*M.rat),
         pred.min = Intercept + (Seed * M.seed) + (Density*min.dens) + (Rats*M.dens),
         pred.max = Intercept + (Seed * M.seed) + (Density*MAX.dens) + (Rats*M.dens))

# # r2
# # plot(pred~mean.r, data = seed.May.plot1)
#
#
# # data in tidyverse for ggplot2
out.dens.May2 <- May.plot1 %>%
  gather(key = pred.fact, value = mean.r, pred.min:pred.max) %>%
  mutate(lag.N = ifelse(pred.fact == "pred.min", min.dens, MAX.dens)) %>%
  select(Valley, Control, mean.r, lag.N) %>%
  mutate(treat = paste(Valley, Control),
         lag.N = ifelse(lag.N>0, lag.N, 0.1))
#
# # dens-May ----------------------------------------------------------------
# #spread
# dens.May.plot1 <- May.plot %>%
#   mutate(pred.min = b0 + (b.dens * dat.mice.min),
#          pred.max = (b0 + b.dens * dat.mice.max))
#
# #tidyverse
# out.dens.May2 <- dens.May.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-May -----------------------------------------------------------------
# #spread and predict
# rat.May.plot1 <- May.plot %>%
#   mutate(pred.min = b0 + (b.dens * dat.rat.min),
#          pred.max = (b0 + b.dens * dat.rat.max))
#
# #tidyverse
# out.rat.May2 <- rat.May.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))



## ------------------------------------------------------------------------
# data ----------------------
### Seed lines (May)

pD.plot.3 <- out.full.136 %>%
  filter(month == "May") %>%
  mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N))

# pc.plot.May <-

  ggplot(pD.plot.3, aes(y = mean.r, x = lag.N)) +
  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r), lwd = 0.75, alpha = 0.2, position=position_dodge(width=30), width = 0) +

  geom_point(aes(colour = Valley,shape = Valley, fill = Control),
             stroke = 1.5, size = 2, alpha = 0.7) +

  # ggplot(out.seed.feb2, aes(y = mean.r, x = lag.sjt)) +
    geom_line(data = out.dens.May2, aes(y = mean.r, x = lag.N, group = treat),size = 0.75, alpha = 0.7) +
  geom_point(data = out.dens.May2, aes(y = mean.r, x = lag.N, shape = Valley, fill = Control, group = treat), size = 5) +

   scale_color_manual(name = "Stoat Control",
                     values = c("white", "black", "white")) +

  scale_shape_manual(name = "Ecosystem",
                     values = c(24, 21)) +
  # scale_size_manual(name = "Rat Control", values = c(2.5, 3, 2.5)) +
  scale_fill_manual(
    name = "Stoat Control",
    values = c("cornflowerblue", "darkorange", "cornflowerblue")
  ) +
     # Remove fill legend and replace the fill legend using the newly created size
  guides(
    col = "none",
    size = guide_legend(override.aes = list(
  shape = c(15,0),alpha = 1
    )),
    shape = guide_legend(override.aes = list(
      shape = c(24, 21), size = 4
    )),
    fill = guide_legend(override.aes = list(
      col = c("cornflowerblue", "darkorange"),shape = c("square"),
      size = 4
    ))
  ) +

  #
  # # scale_colour_manual(name = "Stoat control",
  # #                     labels = c("E-", "H+", "H-"),
  # #                     values = c("black","black", "black")) +
  # # scale_shape_manual(name = "Valley",
  # #                    labels = c("E", "H"),
  # #                    values = c(25,21)) +
  # # scale_fill_manual(name = "Stoat control",
  # #                   labels = c("E-", "H+", "H-"),
  # #                   values = c("white","black", "white")) +


# xlab(expression(paste("Intake"," ", "Rate"," ","(",italic(S[jt]),")" ))) +
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]),")"))) )+

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
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", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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"))
# pd.plot.May


## ------------------------------------------------------------------------
# data for lines ----------------------------------------------------------
# Aug - density
#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 %>%
  filter(month == "Aug") %>%
  mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N)) %>%
  select(-var) %>%
    droplevels()

out.dat1.Aug2 <-  out.full.136 %>%
  filter(month == "Aug") %>%
      mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N)) %>%
  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 = 0,
    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))
Aug.plot1 <- Aug.plot %>%
  mutate(pred.mean = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*M.rat),
         pred.min = Intercept + (Seed * M.seed) + (Density*min.dens) + (Rats*M.dens),
         pred.max = Intercept + (Seed * M.seed) + (Density*MAX.dens) + (Rats*M.dens))

# # r2
# # plot(pred~mean.r, data = seed.Aug.plot1)
#
#
# # data in tidyverse for ggplot2
out.dens.Aug2 <- Aug.plot1 %>%
  gather(key = pred.fact, value = mean.r, pred.min:pred.max) %>%
  mutate(lag.N = ifelse(pred.fact == "pred.min", min.dens, MAX.dens)) %>%
  select(Valley, Control, mean.r, lag.N) %>%
  mutate(treat = paste(Valley, Control),
         lag.N = ifelse(lag.N>0, lag.N, 0.1))



## ------------------------------------------------------------------------
# data ----------------------
### Seed lines (Aug)

pD.plot.3 <- out.full.136 %>%
  filter(month == "Aug") %>%
  mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N))

# pc.plot.Aug <-

  ggplot(pD.plot.3, aes(y = mean.r, x = lag.N)) +

  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r), lwd = 0.75, alpha = 0.2, position=position_dodge(width=30), width = 0) +

  geom_point(aes(colour = Valley,shape = Valley, fill = Control),
             stroke = 1.5, size = 2, alpha = 0.7) +

  # ggplot(out.seed.Aug2, aes(y = mean.r, x = lag.sjt)) +
    geom_line(data = out.dens.Aug2, aes(y = mean.r, x = lag.N, group = treat),size = 0.75, alpha = 0.7) +
  geom_point(data = out.dens.Aug2, aes(y = mean.r, x = lag.N, shape = Valley, fill = Control, group = treat), size = 5) +

 scale_color_manual(name = "Stoat Control",
                     values = c("white", "black", "white")) +

  scale_shape_manual(name = "Ecosystem",
                     values = c(24, 21)) +
  # scale_size_manual(name = "Rat Control", values = c(2.5, 3, 2.5)) +
  scale_fill_manual(
    name = "Stoat Control",
    values = c("cornflowerblue", "darkorange", "cornflowerblue")
  ) +
     # Remove fill legend and replace the fill legend using the newly created size
  guides(
    col = "none",
    size = guide_legend(override.aes = list(
  shape = c(15,0),alpha = 1
    )),
    shape = guide_legend(override.aes = list(
      shape = c(24, 21), size = 4
    )),
    fill = guide_legend(override.aes = list(
      col = c("cornflowerblue", "darkorange"),shape = c("square"),
      size = 4
    ))
  ) +

  #
  # # scale_colour_manual(name = "Stoat control",
  # #                     labels = c("E-", "H+", "H-"),
  # #                     values = c("black","black", "black")) +
  # # scale_shape_manual(name = "Valley",
  # #                    labels = c("E", "H"),
  # #                    values = c(25,21)) +
  # # scale_fill_manual(name = "Stoat control",
  # #                   labels = c("E-", "H+", "H-"),
  # #                   values = c("white","black", "white")) +


# xlab(expression(paste("Intake"," ", "Rate"," ","(",italic(S[jt]),")" ))) +
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]),")"))) )+

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
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", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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"))
# pd.plot.May


## ------------------------------------------------------------------------
# data for lines ----------------------------------------------------------
# Nov - seed
#Nov filtered (parameters)
# flip dataset
out.para.Nov <- para.plot.dat %>%
  filter(month == "Nov") %>%
  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 ....

#Nov filtered (data)
out.dat1.Nov <- out.full.136 %>%
  filter(month == "Nov") %>%
  mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N)) %>%
  select(-var) %>%
    droplevels()

out.dat1.Nov2 <-  out.full.136 %>%
  filter(month == "Nov") %>%
      mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N)) %>%
  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 = 0,
    min.rat = 0,
    min.r = min(mean.r, na.rm = TRUE)) %>%
  ungroup()

# #merge two datasets
Nov.plot <- left_join(out.dat1.Nov2 ,out.para.Nov)

# predict estimates for seed lines
# seed.Nov.plot1 <- Nov.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))
Nov.plot1 <- Nov.plot %>%
  mutate(pred.mean = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*M.rat),
         pred.min = Intercept + (Seed * M.seed) + (Density*min.dens) + (Rats*M.dens),
         pred.max = Intercept + (Seed * M.seed) + (Density*MAX.dens) + (Rats*M.dens))

# # r2
# # plot(pred~mean.r, data = seed.Nov.plot1)
#
#
# # data in tidyverse for ggplot2
out.dens.Nov2 <- Nov.plot1 %>%
  gather(key = pred.fact, value = mean.r, pred.min:pred.max) %>%
  mutate(lag.N = ifelse(pred.fact == "pred.min", min.dens, MAX.dens)) %>%
  select(Valley, Control, mean.r, lag.N) %>%
  mutate(treat = paste(Valley, Control),
         lag.N = ifelse(lag.N>0, lag.N, 0.1))
#
# # dens-Nov ----------------------------------------------------------------
# #spread
# dens.Nov.plot1 <- Nov.plot %>%
#   mutate(pred.min = b0 + (b.dens * dat.mice.min),
#          pred.max = (b0 + b.dens * dat.mice.max))
#
# #tidyverse
# out.dens.Nov2 <- dens.Nov.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-Nov -----------------------------------------------------------------
# #spread and predict
# rat.Nov.plot1 <- Nov.plot %>%
#   mutate(pred.min = b0 + (b.dens * dat.rat.min),
#          pred.max = (b0 + b.dens * dat.rat.max))
#
# #tidyverse
# out.rat.Nov2 <- rat.Nov.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))



## ------------------------------------------------------------------------
# data ----------------------
### Seed lines (Nov)

pD.plot.3 <- out.full.136 %>%
  filter(month == "Nov") %>%
  mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N))

# pc.plot.Nov <-

  ggplot(pD.plot.3, aes(y = mean.r, x = lag.N)) +

  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r), lwd = 0.75, alpha = 0.2, position=position_dodge(width=30), width = 0) +

  geom_point(aes(colour = Valley,shape = Valley, fill = Control),
             stroke = 1.5, size = 2, alpha = 0.7) +

  # ggplot(out.seed.Nov2, aes(y = mean.r, x = lag.sjt)) +
    geom_line(data = out.dens.Nov2, aes(y = mean.r, x = lag.N, group = treat),size = 0.75, alpha = 0.7) +
  geom_point(data = out.dens.Nov2, aes(y = mean.r, x = lag.N, shape = Valley, fill = Control, group = treat), size = 5) +

# facet_wrap(~month, scales = "free")+
  scale_shape_manual(name = "Valley",
                         labels = c("E", "H"),
                         values = c(25,21)) +

  scale_colour_manual(name = "Stoat control",
                      labels = c("Eglinton", "Hollyford", "Hollyford"),
                      values = c("darkgoldenrod","black", "black")) +
  scale_fill_manual(name = "Stoat control",
                    labels = c("Yes", "No", "Yes"),
                    values = c("darkgoldenrod","black", "darkgoldenrod")) +
  #
  # # scale_colour_manual(name = "Stoat control",
  # #                     labels = c("E-", "H+", "H-"),
  # #                     values = c("black","black", "black")) +
  # # scale_shape_manual(name = "Valley",
  # #                    labels = c("E", "H"),
  # #                    values = c(25,21)) +
  # # scale_fill_manual(name = "Stoat control",
  # #                   labels = c("E-", "H+", "H-"),
  # #                   values = c("white","black", "white")) +


# xlab(expression(paste("Intake"," ", "Rate"," ","(",italic(S[jt]),")" ))) +
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]),")"))) )+

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
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", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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"))


## ----include = FALSE-----------------------------------------------------
# data for lines ----------------------------------------------------------
# Feb - rats
#feb filtered (parameters)
# flip dataset
para.plot.dat1 <- para.plot.dat

out.para1 <- para.plot.dat %>%
  select(mean.b, Control, Valley, para, month) %>%
    droplevels() %>%
  spread(key = para, value = mean.b)

# for each group I need to ....

#feb filtered (data)
out.dat1 <- out.full.136 %>%
  mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N),
         lag.R = lag.rat.mna) %>%
  select(-var) %>%
    droplevels()

out.dat2 <-  out.full.136 %>%
      mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N),
         lag.R = lag.rat.mna) %>%
  group_by(Control, Valley, month) %>%
    summarise(
    M.seed = mean(lag.sjt, na.rm = TRUE),
    M.dens = mean(lag.N, na.rm = TRUE),
    M.rat = mean(lag.R, na.rm = TRUE),

    MAX.rat = max(lag.R, 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
dens.plots <- left_join(out.dat2, out.para1, by = c("Control", "Valley", "month"))

# predict estimates for seed lines
# seed.plot1 <- seed.plots %>%
#   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))
dens.plot1 <- dens.plots %>%
  mutate(pred.mean = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*M.rat),
         pred.min = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*min.rat),
         pred.max = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*MAX.rat))

# # r2
# # plot(pred~mean.r, data = seed.feb.plot1)
#
#
# # data in tidyverse for ggplot2
out.dens2 <- dens.plot1 %>%
  gather(key = pred.fact, value = mean.r, pred.min:pred.max) %>%
  mutate(lag.R = ifelse(pred.fact == "pred.min", min.rat, MAX.rat)) %>%
  select(Valley, Control, mean.r, lag.R, month, pred.fact) %>%
  mutate(lag.R = ifelse(lag.R>0, lag.R, 0.1),
         treat = paste(Valley,Control))


#
# # dens-feb ----------------------------------------------------------------
# #spread
# dens.feb.plot1 <- feb.plot %>%
#   mutate(pred.min = b0 + (b.dens * dat.mice.min),
#          pred.max = (b0 + b.dens * dat.mice.max))
#
# #tidyverse
# out.dens.feb2 <- dens.feb.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-feb -----------------------------------------------------------------
# #spread and predict
# rat.feb.plot1 <- feb.plot %>%
#   mutate(pred.min = b0 + (b.dens * dat.rat.min),
#          pred.max = (b0 + b.dens * dat.rat.max))
#
# #tidyverse
# out.rat.feb2 <- rat.feb.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))



## ----include = FALSE-----------------------------------------------------
# data ----------------------

pC.plot.3 <- out.full.136

out.dat1 <- out.full.136 %>%
  mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N),
         lag.R = lag.rat.mna) %>%
  select(-var) %>%
    droplevels()

ggplot(out.dat1, aes(y = mean.r, x = lag.R)) +

  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r), lwd = 0.75, alpha = 0.2, position=position_dodge(width=30), width = 0) +

  geom_point(aes(colour = Valley,shape = Valley, fill = Control),
             stroke = 1.5, size = 2, alpha = 0.7) +

  # ggplot(out.seed.feb2, aes(y = mean.r, x = lag.sjt)) +
    geom_line(data = out.dens2, aes(y = mean.r, x = lag.R, group = treat),size = 0.75, alpha = 0.7) +
  geom_point(data = out.dens2, aes(y = mean.r, x = lag.R, shape = Valley, fill = Control, group = treat), size = 5) +

facet_grid(Control~month, scales = "free") +
  scale_shape_manual(name = "Valley",
                         labels = c("E", "H"),
                         values = c(25,21)) +

  scale_colour_manual(name = "Stoat control",
                      labels = c("Eglinton", "Hollyford", "Hollyford"),
                      values = c("darkgoldenrod","black", "black")) +
  scale_fill_manual(name = "Stoat control",
                    labels = c("Yes", "No", "Yes"),
                    values = c("darkgoldenrod","black", "darkgoldenrod")) +
  #
  # # scale_colour_manual(name = "Stoat control",
  # #                     labels = c("E-", "H+", "H-"),
  # #                     values = c("black","black", "black")) +
  # # scale_shape_manual(name = "Valley",
  # #                    labels = c("E", "H"),
  # #                    values = c(25,21)) +
  # # scale_fill_manual(name = "Stoat control",
  # #                   labels = c("E-", "H+", "H-"),
  # #                   values = c("white","black", "white")) +


# 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]),")"))) )+

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
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", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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"))



## ------------------------------------------------------------------------
# data for lines ----------------------------------------------------------
# Feb - rats
#feb filtered (parameters)
# flip dataset
para.plot.dat1 <- para.plot.dat %>%
  filter(month == "Feb")

out.para1 <- para.plot.dat %>%
  filter(month == "Feb") %>%
  select(mean.b, Control, Valley, para, month) %>%
    droplevels() %>%
  spread(key = para, value = mean.b)

# for each group I need to ....

#feb filtered (data)
out.dat1 <- out.full.136 %>%
  filter(month == "Feb") %>%
  mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N),
         lag.R = lag.rat.mna) %>%
  select(-var) %>%
    droplevels()

out.dat2 <-  out.full.136 %>%
  filter(month == "Feb") %>%
      mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N),
         lag.R = lag.rat.mna) %>%
  group_by(Control, Valley, month) %>%
    summarise(
    M.seed = mean(lag.sjt, na.rm = TRUE),
    M.dens = mean(lag.N, na.rm = TRUE),
    M.rat = mean(lag.R, na.rm = TRUE),

    MAX.rat = max(lag.R, 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
dens.plots <- left_join(out.dat2, out.para1, by = c("Control", "Valley", "month"))

# predict estimates for seed lines
# seed.plot1 <- seed.plots %>%
#   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))
dens.plot1 <- dens.plots %>%
  mutate(pred.mean = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*M.rat),
         pred.min = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*min.rat),
         pred.max = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*MAX.rat))

# # r2
# # plot(pred~mean.r, data = seed.feb.plot1)
#
#
# # data in tidyverse for ggplot2
out.dens2 <- dens.plot1 %>%
  gather(key = pred.fact, value = mean.r, pred.min:pred.max) %>%
  mutate(lag.R = ifelse(pred.fact == "pred.min", min.rat, MAX.rat)) %>%
  select(Valley, Control, mean.r, lag.R, month, pred.fact) %>%
  mutate(lag.R = ifelse(lag.R>0, lag.R, 0.1),
         treat = paste(Valley,Control))


## ------------------------------------------------------------------------
# data ----------------------

pC.plot.3 <- out.full.136 %>%
  filter(month == "Feb")

out.dat1 <- out.full.136 %>%
  filter(month == "Feb") %>%
  mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N),
         lag.R = lag.rat.mna) %>%
  select(-var) %>%
    droplevels()

ggplot(out.dat1, aes(y = mean.r, x = lag.R)) +

  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r), lwd = 0.75, alpha = 0.2, position=position_dodge(width=30), width = 0) +

  geom_point(aes(colour = Valley,shape = Valley, fill = Control),
             stroke = 1.5, size = 2, alpha = 0.7) +

  # ggplot(out.seed.feb2, aes(y = mean.r, x = lag.sjt)) +
    geom_line(data = out.dens2, aes(y = mean.r, x = lag.R, group = treat),size = 0.75, alpha = 0.7) +
  geom_point(data = out.dens2, aes(y = mean.r, x = lag.R, shape = Valley, fill = Control, group = treat), size = 5) +

 scale_color_manual(name = "Stoat Control",
                     values = c("white", "black", "white")) +

  scale_shape_manual(name = "Ecosystem",
                     values = c(24, 21)) +
  # scale_size_manual(name = "Rat Control", values = c(2.5, 3, 2.5)) +
  scale_fill_manual(
    name = "Stoat Control",
    values = c("cornflowerblue", "darkorange", "cornflowerblue")
  ) +
     # Remove fill legend and replace the fill legend using the newly created size
  guides(
    col = "none",
    size = guide_legend(override.aes = list(
  shape = c(15,0),alpha = 1
    )),
    shape = guide_legend(override.aes = list(
      shape = c(24, 21), size = 4
    )),
    fill = guide_legend(override.aes = list(
      col = c("cornflowerblue", "darkorange"),shape = c("square"),
      size = 4
    ))
  ) +

  #
  # # scale_colour_manual(name = "Stoat control",
  # #                     labels = c("E-", "H+", "H-"),
  # #                     values = c("black","black", "black")) +
  # # scale_shape_manual(name = "Valley",
  # #                    labels = c("E", "H"),
  # #                    values = c(25,21)) +
  # # scale_fill_manual(name = "Stoat control",
  # #                   labels = c("E-", "H+", "H-"),
  # #                   values = c("white","black", "white")) +


# xlab(expression(paste("Intake"," ", "Rate"," ","(",italic(S[jt]),")" ))) +
# xlab(expression(paste("Mouse"," ", "Density"," ","(",italic(N[jt]),")" ))) +
xlab(expression(atop(paste("Minimum"," ", "number"," "),paste("of", " ", "rats"," ","(",italic(R[jt]),")"))) )+

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
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", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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"))
# pd.plot.May



## ------------------------------------------------------------------------
# data for lines ----------------------------------------------------------
# May - rats
#May filtered (parameters)
# flip dataset
para.plot.dat1 <- para.plot.dat %>%
  filter(month == "May")

out.para1 <- para.plot.dat %>%
  filter(month == "May") %>%
  select(mean.b, Control, Valley, para, month) %>%
    droplevels() %>%
  spread(key = para, value = mean.b)

# for each group I need to ....

#May filtered (data)
out.dat1 <- out.full.136 %>%
  filter(month == "May") %>%
  mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N),
         lag.R = lag.rat.mna) %>%
  select(-var) %>%
    droplevels()

out.dat2 <-  out.full.136 %>%
  filter(month == "May") %>%
      mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N),
         lag.R = lag.rat.mna) %>%
  group_by(Control, Valley, month) %>%
    summarise(
    M.seed = mean(lag.sjt, na.rm = TRUE),
    M.dens = mean(lag.N, na.rm = TRUE),
    M.rat = mean(lag.R, na.rm = TRUE),

    MAX.rat = max(lag.R, 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
dens.plots <- left_join(out.dat2, out.para1, by = c("Control", "Valley", "month"))

# predict estimates for seed lines
# seed.plot1 <- seed.plots %>%
#   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))
dens.plot1 <- dens.plots %>%
  mutate(pred.mean = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*M.rat),
         pred.min = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*min.rat),
         pred.max = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*MAX.rat))

# # r2
# # plot(pred~mean.r, data = seed.May.plot1)
#
#
# # data in tidyverse for ggplot2
out.dens2 <- dens.plot1 %>%
  gather(key = pred.fact, value = mean.r, pred.min:pred.max) %>%
  mutate(lag.R = ifelse(pred.fact == "pred.min", min.rat, MAX.rat)) %>%
  select(Valley, Control, mean.r, lag.R, month, pred.fact) %>%
  mutate(lag.R = ifelse(lag.R>0, lag.R, 0.1),
         treat = paste(Valley,Control))


## ------------------------------------------------------------------------
# data ----------------------

pC.plot.3 <- out.full.136 %>%
  filter(month == "May")

out.dat1 <- out.full.136 %>%
  filter(month == "May") %>%
  mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N),
         lag.R = lag.rat.mna) %>%
  select(-var) %>%
    droplevels()

ggplot(out.dat1, aes(y = mean.r, x = lag.R)) +

  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r), lwd = 0.75, alpha = 0.2, position=position_dodge(width=30), width = 0) +

  geom_point(aes(colour = Valley,shape = Valley, fill = Control),
             stroke = 1.5, size = 2, alpha = 0.7) +

  # ggplot(out.seed.May2, aes(y = mean.r, x = lag.sjt)) +
    geom_line(data = out.dens2, aes(y = mean.r, x = lag.R, group = treat),size = 0.75, alpha = 0.7) +
  geom_point(data = out.dens2, aes(y = mean.r, x = lag.R, shape = Valley, fill = Control, group = treat), size = 5) +

 scale_color_manual(name = "Stoat Control",
                     values = c("white", "black", "white")) +

  scale_shape_manual(name = "Ecosystem",
                     values = c(24, 21)) +
  # scale_size_manual(name = "Rat Control", values = c(2.5, 3, 2.5)) +
  scale_fill_manual(
    name = "Stoat Control",
    values = c("cornflowerblue", "darkorange", "cornflowerblue")
  ) +
     # Remove fill legend and replace the fill legend using the newly created size
  guides(
    col = "none",
    size = guide_legend(override.aes = list(
  shape = c(15,0),alpha = 1
    )),
    shape = guide_legend(override.aes = list(
      shape = c(24, 21), size = 4
    )),
    fill = guide_legend(override.aes = list(
      col = c("cornflowerblue", "darkorange"),shape = c("square"),
      size = 4
    ))
  ) +

  #
  # # scale_colour_manual(name = "Stoat control",
  # #                     labels = c("E-", "H+", "H-"),
  # #                     values = c("black","black", "black")) +
  # # scale_shape_manual(name = "Valley",
  # #                    labels = c("E", "H"),
  # #                    values = c(25,21)) +
  # # scale_fill_manual(name = "Stoat control",
  # #                   labels = c("E-", "H+", "H-"),
  # #                   values = c("white","black", "white")) +


# xlab(expression(paste("Intake"," ", "Rate"," ","(",italic(S[jt]),")" ))) +
# xlab(expression(paste("Mouse"," ", "Density"," ","(",italic(N[jt]),")" ))) +
xlab(expression(atop(paste("Minimum"," ", "number"," "),paste("of", " ", "rats"," ","(",italic(R[jt]),")"))) )+

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
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", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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"))
# pd.plot.May



## ------------------------------------------------------------------------
# data for lines ----------------------------------------------------------
# Aug - rats
#Aug filtered (parameters)
# flip dataset
para.plot.dat1 <- para.plot.dat %>%
  filter(month == "Aug")

out.para1 <- para.plot.dat %>%
  filter(month == "Aug") %>%
  select(mean.b, Control, Valley, para, month) %>%
    droplevels() %>%
  spread(key = para, value = mean.b)

# for each group I need to ....

#Aug filtered (data)
out.dat1 <- out.full.136 %>%
  filter(month == "Aug") %>%
  mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N),
         lag.R = lag.rat.mna) %>%
  select(-var) %>%
    droplevels()

out.dat2 <-  out.full.136 %>%
  filter(month == "Aug") %>%
      mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N),
         lag.R = lag.rat.mna) %>%
  group_by(Control, Valley, month) %>%
    summarise(
    M.seed = mean(lag.sjt, na.rm = TRUE),
    M.dens = mean(lag.N, na.rm = TRUE),
    M.rat = mean(lag.R, na.rm = TRUE),

    MAX.rat = max(lag.R, 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
dens.plots <- left_join(out.dat2, out.para1, by = c("Control", "Valley", "month"))

# predict estimates for seed lines
# seed.plot1 <- seed.plots %>%
#   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))
dens.plot1 <- dens.plots %>%
  mutate(pred.mean = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*M.rat),
         pred.min = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*min.rat),
         pred.max = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*MAX.rat))

# # r2
# # plot(pred~mean.r, data = seed.Aug.plot1)
#
#
# # data in tidyverse for ggplot2
out.dens2 <- dens.plot1 %>%
  gather(key = pred.fact, value = mean.r, pred.min:pred.max) %>%
  mutate(lag.R = ifelse(pred.fact == "pred.min", min.rat, MAX.rat)) %>%
  select(Valley, Control, mean.r, lag.R, month, pred.fact) %>%
  mutate(lag.R = ifelse(lag.R>0, lag.R, 0.1),
         treat = paste(Valley,Control))


## ------------------------------------------------------------------------
# data ----------------------

pC.plot.3 <- out.full.136 %>%
  filter(month == "Aug")

out.dat1 <- out.full.136 %>%
  filter(month == "Aug") %>%
  mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N),
         lag.R = lag.rat.mna) %>%
  select(-var) %>%
    droplevels()

ggplot(out.dat1, aes(y = mean.r, x = lag.R)) +

  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r), lwd = 0.75, alpha = 0.2, position=position_dodge(width=30), width = 0) +

  geom_point(aes(colour = Valley,shape = Valley, fill = Control),
             stroke = 1.5, size = 2, alpha = 0.7) +

  # ggplot(out.seed.Aug2, aes(y = mean.r, x = lag.sjt)) +
    geom_line(data = out.dens2, aes(y = mean.r, x = lag.R, group = treat),size = 0.75, alpha = 0.7) +
  geom_point(data = out.dens2, aes(y = mean.r, x = lag.R, shape = Valley, fill = Control, group = treat), size = 5) +

 scale_color_manual(name = "Stoat Control",
                     values = c("white", "black", "white")) +

  scale_shape_manual(name = "Ecosystem",
                     values = c(24, 21)) +
  # scale_size_manual(name = "Rat Control", values = c(2.5, 3, 2.5)) +
  scale_fill_manual(
    name = "Stoat Control",
    values = c("cornflowerblue", "darkorange", "cornflowerblue")
  ) +
     # Remove fill legend and replace the fill legend using the newly created size
  guides(
    col = "none",
    size = guide_legend(override.aes = list(
  shape = c(15,0),alpha = 1
    )),
    shape = guide_legend(override.aes = list(
      shape = c(24, 21), size = 4
    )),
    fill = guide_legend(override.aes = list(
      col = c("cornflowerblue", "darkorange"),shape = c("square"),
      size = 4
    ))
  ) +

  #
  # # scale_colour_manual(name = "Stoat control",
  # #                     labels = c("E-", "H+", "H-"),
  # #                     values = c("black","black", "black")) +
  # # scale_shape_manual(name = "Valley",
  # #                    labels = c("E", "H"),
  # #                    values = c(25,21)) +
  # # scale_fill_manual(name = "Stoat control",
  # #                   labels = c("E-", "H+", "H-"),
  # #                   values = c("white","black", "white")) +


# xlab(expression(paste("Intake"," ", "Rate"," ","(",italic(S[jt]),")" ))) +
# 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]),")"))) )+

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
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", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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"))
# pd.plot.Aug



## ------------------------------------------------------------------------
# data for lines ----------------------------------------------------------
# Nov - rats
#Nov filtered (parameters)
# flip dataset
para.plot.dat1 <- para.plot.dat %>%
  filter(month == "Nov")

out.para1 <- para.plot.dat %>%
  filter(month == "Nov") %>%
  select(mean.b, Control, Valley, para, month) %>%
    droplevels() %>%
  spread(key = para, value = mean.b)

# for each group I need to ....

#Nov filtered (data)
out.dat1 <- out.full.136 %>%
  filter(month == "Nov") %>%
  mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N),
         lag.R = lag.rat.mna) %>%
  select(-var) %>%
    droplevels()

out.dat2 <-  out.full.136 %>%
  filter(month == "Nov") %>%
      mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N),
         lag.R = lag.rat.mna) %>%
  group_by(Control, Valley, month) %>%
    summarise(
    M.seed = mean(lag.sjt, na.rm = TRUE),
    M.dens = mean(lag.N, na.rm = TRUE),
    M.rat = mean(lag.R, na.rm = TRUE),

    MAX.rat = max(lag.R, 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
dens.plots <- left_join(out.dat2, out.para1, by = c("Control", "Valley", "month"))

# predict estimates for seed lines
# seed.plot1 <- seed.plots %>%
#   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))
dens.plot1 <- dens.plots %>%
  mutate(pred.mean = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*M.rat),
         pred.min = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*min.rat),
         pred.max = Intercept + (Seed * M.seed) + (Density*M.dens) + (Rats*MAX.rat))

# # r2
# # plot(pred~mean.r, data = seed.Nov.plot1)
#
#
# # data in tidyverse for ggplot2
out.dens2 <- dens.plot1 %>%
  gather(key = pred.fact, value = mean.r, pred.min:pred.max) %>%
  mutate(lag.R = ifelse(pred.fact == "pred.min", min.rat, MAX.rat)) %>%
  select(Valley, Control, mean.r, lag.R, month, pred.fact) %>%
  mutate(lag.R = ifelse(lag.R>0, lag.R, 0.1),
         treat = paste(Valley,Control))


## ------------------------------------------------------------------------
# data ----------------------

pC.plot.3 <- out.full.136 %>%
  filter(month == "Nov")

out.dat1 <- out.full.136 %>%
  filter(month == "Nov") %>%
  mutate(lag.sjt = lag(log.seed),
         lag.N = lag(N),
         lag.R = lag.rat.mna) %>%
  select(-var) %>%
    droplevels()

ggplot(out.dat1, aes(y = mean.r, x = lag.R)) +

  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r), lwd = 0.75, alpha = 0.2, position=position_dodge(width=30), width = 0) +

  geom_point(aes(colour = Valley,shape = Valley, fill = Control),
             stroke = 1.5, size = 2, alpha = 0.7) +

  # ggplot(out.seed.Nov2, aes(y = mean.r, x = lag.sjt)) +
    geom_line(data = out.dens2, aes(y = mean.r, x = lag.R, group = treat),size = 0.75, alpha = 0.7) +
  geom_point(data = out.dens2, aes(y = mean.r, x = lag.R, shape = Valley, fill = Control, group = treat), size = 5) +

 scale_color_manual(name = "Stoat Control",
                     values = c("white", "black", "white")) +

  scale_shape_manual(name = "Ecosystem",
                     values = c(24, 21)) +
  # scale_size_manual(name = "Rat Control", values = c(2.5, 3, 2.5)) +
  scale_fill_manual(
    name = "Stoat Control",
    values = c("cornflowerblue", "darkorange", "cornflowerblue")
  ) +
     # Remove fill legend and replace the fill legend using the newly created size
  guides(
    col = "none",
    size = guide_legend(override.aes = list(
  shape = c(15,0),alpha = 1
    )),
    shape = guide_legend(override.aes = list(
      shape = c(24, 21), size = 4
    )),
    fill = guide_legend(override.aes = list(
      col = c("cornflowerblue", "darkorange"),shape = c("square"),
      size = 4
    ))
  ) +

  #
  # # scale_colour_manual(name = "Stoat control",
  # #                     labels = c("E-", "H+", "H-"),
  # #                     values = c("black","black", "black")) +
  # # scale_shape_manual(name = "Valley",
  # #                    labels = c("E", "H"),
  # #                    values = c(25,21)) +
  # # scale_fill_manual(name = "Stoat control",
  # #                   labels = c("E-", "H+", "H-"),
  # #                   values = c("white","black", "white")) +


# xlab(expression(paste("Intake"," ", "Rate"," ","(",italic(S[jt]),")" ))) +
# xlab(expression(paste("Mouse"," ", "Density"," ","(",italic(N[jt]),")" ))) +
xlab(expression(atop(paste("Minimum"," ", "number"," "),paste("of", " ", "rats"," ","(",italic(R[jt]),")"))) ) +

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
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", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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"))
# pd.plot.Nov



## ------------------------------------------------------------------------
para.rats <- para.plot.erors %>%
              filter(para == "b.rat") %>%
                mutate(month = factor(month, levels = c("Feb", "May", "Aug", "Nov"), labels = c("February", "May", "August", "November")))



para.plot.rat <- ggplot(para.rats, aes(y = mean.b, x = month, col = Control, shape = Valley, fill = Control)) +
  geom_errorbar(aes(ymin = lcl.b, ymax = ucl.b), lwd = 0.75, alpha = 1, width = 0.1, position =position_dodge(width = 1)) +

  geom_abline(intercept = 0, slope = 0) +
  geom_point(size =4, alpha = 0.7, position = position_dodge(width = 1)) +
  # facet_wrap(month~para, scales = "free") +
  facet_wrap(~month, scales = "free") +
  # scale_y_continuous(limits = c(-1,1)) +
   scale_color_manual(name = "Stoat Control",
                     values = c("cornflowerblue", "darkorange", "cornflowerblue")) +

  scale_shape_manual(name = "Ecosystem",
                     values = c(24, 21)) +
  scale_size_manual(name = "Rat Control", values = c(2.5, 3, 2.5)) +
  scale_fill_manual(name = "Stoat Control",
    values = c("cornflowerblue", "darkorange", "cornflowerblue")) +
     # Remove fill legend and replace the fill legend using the newly created size
  guides(
    col = "none",
    size = guide_legend(override.aes = list(
  shape = c(15,0),alpha = 1
    )),
    shape = guide_legend(override.aes = list(
      shape = c(24, 21), size = 4
    )),
    fill = guide_legend(override.aes = list(
      col = c("cornflowerblue", "darkorange"),shape = c("square"),
      size = 4
    ))
  ) +


# xlab(expression(paste("Intake"," ", "Rate"," ","(",italic(S[jt]),")" ))) +
xlab("") +
# xlab(expression(atop(paste("Minimum"," ", "number"," "),paste("of", " ", "rats"," ","(",italic(R[jt]),")"))) )+
# xlab(expression(atop(paste("Minimum"," ", "number"," "),paste("of", " ", "rats"," ","(",italic(R[jt]),")"))) )+

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
theme_tufte() +
  theme_bw() +
  theme(strip.background = element_blank(),
        strip.text.y = element_blank(),

        plot.title = element_blank(),
        plot.subtitle= element_blank(),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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 =14, family = "Times", angle = 90),
        axis.title.x = element_text(colour = "black", size =14, family = "Times"),
        axis.text.y=element_text(colour = "black",size = 14, family = "Times"),
        axis.text.x = element_text(colour = "black", size =14, family = "Times"),

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

        strip.text = element_blank())
para.plot.rat



## ----pc-seed-plot--------------------------------------------------------
pd.plot.Aug

# Save plot
# export plot for example vignette
png("./figs/fig-5-1.png")
pd.plot.Aug
dev.off()


## ------------------------------------------------------------------------
knitr::include_graphics(c("./figs/fig-5-1.png"), dpi = 400)


## ----pd-dens-plot--------------------------------------------------------
pd.plot.feb

# Save plot
# export plot for example vignette
png("./figs/fig-6-1.png")
pd.plot.feb
dev.off()


## ------------------------------------------------------------------------
knitr::include_graphics(c("./figs/fig-6-1.png"))


## ----rat-plot------------------------------------------------------------
rat.plot.all <- ggplot(out.dat1, aes(y = mean.r, x = lag.R)) +

  geom_errorbar(aes(ymin = lcl.r, ymax = ucl.r), lwd = 0.75, alpha = 0.2, position=position_dodge(width=30), width = 0) +

  geom_point(aes(colour = Control,shape = Valley, fill = Control),
             stroke = 1.5, size = 2, alpha = 0.7) +
    geom_line(data = out.dens2, aes(y = mean.r, x = lag.R, group = treat),size = 0.75, alpha = 0.7) +
  geom_point(data = out.dens2, aes(y = mean.r, x = lag.R, shape = Valley, fill = Control, group = treat), size = 5) +
 scale_color_manual(name = "Stoat Control",
                     values = c("white", "black", "white")) +

  scale_shape_manual(name = "Ecosystem",
                     values = c(24, 21)) +
  scale_size_manual(name = "Rat Control", values = c(2.5, 3, 2.5)) +
  scale_fill_manual(name = "Stoat Control",
    values = c("cornflowerblue", "darkorange", "cornflowerblue")) +
     # Remove fill legend and replace the fill legend using the newly created size
  guides(
    col = "none",
    size = guide_legend(override.aes = list(
  shape = c(15,0),alpha = 1
    )),
    shape = guide_legend(override.aes = list(
      shape = c(24, 21), size = 4
    )),
    fill = guide_legend(override.aes = list(
      col = c("cornflowerblue", "darkorange"),shape = c("square"),
      size = 4
    ))
  ) +

# 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]),")"))) )+

ylab(expression(atop(paste("Rate"," ", "of"," ",
                           "increase"),paste(" ", "of"," ",
                                             "mice"," ","(",italic(r[jt]),")"))) ) +
# scale_y_continuous(expand = c(0,0.01),breaks = seq(-4,4,1)) +
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", margin = margin(t = 10, b = 10)),
        plot.subtitle=element_text(size=16, face="italic", color="black"),

        legend.position = "none",
        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 =16, family = "Times"),
        legend.title = element_text(colour = "black", size =16, 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"))

rat.plot.all
davan690/beech-publication-wr documentation built on March 29, 2020, 11:09 a.m.
rdrr.io home R language documentation Run R code online
CRAN packages Bioconductor packages R-Forge packages GitHub packages
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
davan690/beech-publication-wr
Small mammal dynamics in NZ beech forests

R/old-scripts/Davidson_2019_Model_wrangling.R
In davan690/beech-publication-wr: Small mammal dynamics in NZ beech forests

R Package Documentation

Browse R Packages

We want your feedback!

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

R/old-scripts/Davidson_2019_Model_wrangling.R In davan690/beech-publication-wr: Small mammal dynamics in NZ beech forests

R Package Documentation

Browse R Packages

We want your feedback!

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

R/old-scripts/Davidson_2019_Model_wrangling.R
In davan690/beech-publication-wr: Small mammal dynamics in NZ beech forests
