In Sz-Tim/gbPopMod: Mechanistic species distribution model for glossy buckthorn (*Frangula alnus*)

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

This is a wrapper for organizing and running sensitivity analyses to assess the impact of each parameter on the resulting buckthorn distribution.

if(!require(pacman)) {install.packages("pacman"); require(pacman)}
p_load("tidyverse", "magrittr", "stringr", "doSNOW", "here", "devtools",
       "fastmatch", "gganimate"); theme_set(theme_bw())
library(gbPopMod)
data(lc.rct)

##---
## set parameters
##---

  p <- "pr.s"
  # p.seq <- seq(0.01, 0.8, length.out=5)
  p.seq <- expand_cnpy(Op=c(0.1, 0.99), Cl=c(0.1, 0.99), length_out=2)
  # p.seq <- expand_LCs(OpI=c(0.3, 0.6), Oth=c(0.01, 0.2), Dec=c(0.2, 0.4),
                      # WP=c(0.2, 0.6), Evg=c(0.2, 0.4), Mxd=c(0.2, 0.4),
                      # length_out=3)
  n.sim <- 3
  g.p <- set_g_p(tmax=100, lc.r=30, lc.c=30, n.ldd=1, n.cores=3)
  control.p <- set_control_p()
  set.seed(225)


##--
## initialize landscape
##--

  # land cover
  lc.df <- lc.rct %>% 
    filter(y >= (max(lc.rct$y) - g.p$lc.r) & x <= g.p$lc.c) %>%
    mutate(id=row_number(), 
           id.in=min_rank(na_if(inbd*id, 0)))
  ggplot(lc.df, aes(x=x, y=-y, fill=inbd)) + geom_tile() +
    scale_fill_manual("Inbounds", values=c("gray", "#543005"))
  ngrid <- nrow(lc.df)
  ncell <- sum(lc.df$inbd)

  # populations
  # indexed 1:ngrid including rows for out-of-bounds cells
  N.init <- pop_init(ngrid, g.p, lc.df)


  # set directories
  sim.wd <- paste0("out/", ncell, "_t", g.p$tmax, "/")
  par.wd <- paste0(sim.wd, p, "/")
  parSet.wd <- paste0(par.wd, p.seq, "/")
  if(!dir.exists(here("out"))) dir.create(here("out"))
  if(!dir.exists(here(sim.wd))) dir.create(here(sim.wd))
  if(!dir.exists(here(par.wd))) dir.create(here(par.wd))


  # dispersal probabilities
  # indexed 1:ncell with only inbound cells: sdd.pr[,,,i] = ith inbound cell
  # row i = dplyr::filter(lc.df, inbds)[i,] = lc.df[id.in==i,]
  # cell indexes in sdd.pr[,,2,] correspond with lc.df$id
  # this extra complication avoids calculating many, many unnecessary cells
  sdd.pr <- sdd_set_probs(ncell, lc.df, g.p)
  saveRDS(sdd.pr, paste0(par.wd, "sdd_pr.rds"))

sdd.pr <- readRDS(paste0(par.wd, "sdd_pr.rds"))
for(j in 1:length(p.seq)) {
  # setup for particular parameter value
  set.seed(225)
  g.p[[p]] <- p.seq[[j]]

  # run n.sim simulations
  p.c <- makeCluster(g.p$n.cores); registerDoSNOW(p.c)
  out.j <- foreach(i=1:n.sim) %dopar% {
    gbPopMod::run_sim(ngrid, ncell, g.p, lc.df, sdd.pr, 
                      N.init, control.p, verbose=F)
  }
  stopCluster(p.c)

  # rearrange output
  ad.j <- map(out.j, ~.$N[,,max(g.p$age.f)]) %>% 
    unlist %>% array(., dim=c(ngrid, g.p$tmax+1, n.sim))
  sb.j <- map(out.j, ~.$N.sb) %>% 
    unlist %>% array(., dim=c(ngrid, g.p$tmax+1, n.sim))

  # store output
  if(!dir.exists(here(parSet.wd[j]))) dir.create(here(parSet.wd[j]))
  save_pars(parSet.wd[j], g.p, control.p)
  save_abundances(parSet.wd[j], ad.j, sb.j)
  rm(ad.j); rm(sb.j)

  # progress
  cat("Finished parameter set", j, "of", length(p.seq), "\n\n")
}

p.c <- makeCluster(g.p$n.cores); registerDoSNOW(p.c)
foreach(j=1:length(p.seq), .combine=rbind) %dopar% {
  require(tidyverse); require(stringr); require(gbPopMod)
  # setup
  p.j <- paste0(p, ": ", p.seq[j])
  ad.j <- readRDS(paste0(parSet.wd[j], "abund_ad.rds"))
  sb.j <- readRDS(paste0(parSet.wd[j], "abund_sb.rds"))
  g.p <- readRDS(paste0(parSet.wd[j], "pars_glbl.rds"))

  # munge data
  ## cell summaries
  ad.mn <- apply(ad.j, 1:2, mean)
  cell.j <- cbind(lc.df, ad.mn) %>% as.tibble %>%
    gather(year, N.adult, (1:ncol(ad.mn)) + ncol(lc.df)) %>%
    rename(ad_Ab=N.adult) %>%
    mutate(ad_sd=c(apply(ad.j, 1:2, sd)),
           ad_pP=c(apply(ad.j > 0, 1:2, mean)),
           sb_Ab=c(apply(log(sb.j+1), 1:2, mean)),
           sb_sd=c(apply(log(sb.j+1), 1:2, sd)),
           sb_pP=c(apply(sb.j > 0, 1:2, mean)),
           ad_L5=c(ad.mn > 0))
  cell.j$ad_L5 <- cell.j$ad_L5 + c(ad.mn > 5)
  cell.j$year <- str_pad(cell.j$year, 3, "left", "0")
  cell.j$p <- p
  if(length(p.seq[[j]])==1) {
    cell.j$p.j <- p.seq[j]
  } else {
    cell.j$p.j <- as.character(p.seq[j])
    cell.j$p.j.OpI <- p.seq[[j]][1]
    cell.j$p.j.Oth <- p.seq[[j]][2]
    cell.j$p.j.Dec <- p.seq[[j]][3]
    cell.j$p.j.WP <- p.seq[[j]][4]
    cell.j$p.j.Evg <- p.seq[[j]][5]
    cell.j$p.j.Mxd <- p.seq[[j]][6]
  }
  ## landscape summaries
  occ.s.ad <- apply(ad.j[lc.df$inbd,,]>0, 2:3, mean)*100
  occ.s.sb <- apply(sb.j[lc.df$inbd,,]>0, 2:3, mean)*100
  less5.s <- apply(ad.j[lc.df$inbd,,]>0 & ad.j[lc.df$inbd,,]<=5, 2:3, mean)*100
  K.ag <- round(as.matrix(lc.df[,4:9]) %*% g.p$K)
  K.s <- apply(ad.j[lc.df$inbd,,]==K.ag[lc.df$inbd,], 2:3, mean)/occ.s.ad*10000
  grid.j <- tibble(year=unique(cell.j$year),
                   pOcc_ad_mn=apply(occ.s.ad, 1, mean),
                   pOcc_ad_sd=apply(occ.s.ad, 1, sd),
                   pOcc_sb_mn=apply(occ.s.sb, 1, mean),
                   pOcc_sb_sd=apply(occ.s.sb, 1, sd),
                   pL5_mn=apply(less5.s, 1, mean),
                   pL5_sd=apply(less5.s, 1, sd),
                   pK_Occ_mn=apply(K.s, 1, mean),
                   pK_Occ_sd=apply(K.s, 1, sd),
                   p=p)
  if(length(p.seq[[j]])==1) {
    grid.j$p.j <- p.seq[j]
  } else {
    grid.j$p.j <- as.character(p.seq[j])
    grid.j$p.j.OpI <- p.seq[[j]][1]
    grid.j$p.j.Oth <- p.seq[[j]][2]
    grid.j$p.j.Dec <- p.seq[[j]][3]
    grid.j$p.j.WP <- p.seq[[j]][4]
    grid.j$p.j.Evg <- p.seq[[j]][5]
    grid.j$p.j.Mxd <- p.seq[[j]][6]
  }


  # save data summaries
  write_csv(cell.j, paste0(parSet.wd[j], "cell_j.csv"))
  write_csv(grid.j, paste0(parSet.wd[j], "grid_j.csv"))
  rm(ad.j); rm(sb.j)

  # save plots
  make_plots_final_t(parSet.wd[j], g.p, filter(cell.j, year==max(cell.j$year)),
                     p.j, 8, 6)
  make_plots_gifs(parSet.wd[j], g.p, cell.j, p.j)
  if(j==1) make_plots_lc(sim.wd, lc.df)
  paste("Finished parameter set", j, "of", length(p.seq))
}
stopCluster(p.c)

# grid summary plots
grid.sum <- map_df(parSet.wd, ~(read_csv(paste0(., "grid_j.csv")))) %>%
  mutate(year=as.numeric(year)) 
if(length(p.seq[[1]])==1) {
  grid.sum %<>% mutate(p.j=as.factor(p.j)) 
} else {
  grid.sum %<>% mutate_at(vars(p.j:p.j.Mxd), as.factor)
}
make_plots_gridSummary(par.wd, grid.sum, byLC=(length(p.seq[[1]])>1))

Sz-Tim/gbPopMod documentation built on Dec. 7, 2020, 1:07 p.m.