code/data_input_v2_nz_land.R

In davan690/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://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")

abund.dat5 <-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))