## ----setup, include=FALSE, echo=FALSE----------------------------------------- library(multiScaleR) knitr::opts_chunk$set(collapse = TRUE, comment = "#>") pkg_extdata <- function(...) { src_path <- normalizePath(file.path("..", "inst", "extdata", ...), winslash = "/", mustWork = FALSE) if (file.exists(src_path)) { return(src_path) } inst_path <- system.file("extdata", ..., package = "multiScaleR") if (nzchar(inst_path) && file.exists(inst_path)) { return(inst_path) } stop("Could not locate required extdata file: ", paste(..., collapse = "/")) } vignette_cache <- readRDS(pkg_extdata("vignette_cache.rds")) # Slow fitted objects used for printed summaries and plots are loaded from # package extdata. Regenerate them with data-raw/generate_vignette_cache.R. ## ----echo=FALSE--------------------------------------------------------------- plot_kernel(prob = 0.9, sigma = 100, kernel = "gaus", scale_dist = F, add_label = F) ## ----echo=FALSE--------------------------------------------------------------- plot_kernel(prob = 0.9, sigma = 50, kernel = "exp", scale_dist = F, add_label = F) ## ----echo=FALSE--------------------------------------------------------------- plot_kernel(prob = 0.9, sigma = 250, beta = 5, kernel = "expow", scale_dist = F, add_label = F) ## ----echo=FALSE--------------------------------------------------------------- plot_kernel(prob = 0.9, sigma = 300, kernel = "fixed", scale_dist = F, add_label = F) ## ----------------------------------------------------------------------------- library(multiScaleR) ## Read in data data("landscape_counts") dat <- landscape_counts data("surv_pts") pts <- vect(surv_pts) land_rast <- terra::rast(system.file("extdata", "landscape.tif", package = 'multiScaleR')) ## ----------------------------------------------------------------------------- summary(dat) pts land_rast plot(land_rast) ## Plot with points plot(land_rast$land1) plot(pts, add = T, pch = 19) ## ----------------------------------------------------------------------------- kernel_inputs <- kernel_prep(pts = pts, raster_stack = land_rast, max_D = 1700, kernel = 'gaussian', verbose = FALSE) kernel_inputs ## ----eval = FALSE------------------------------------------------------------- # kernel_inputs <- kernel_prep(pts = pts, # raster_stack = land_rast, # max_D = 1700, # kernel = 'gaussian', # store_cell_data = FALSE, # verbose = FALSE) ## ----------------------------------------------------------------------------- df <- data.frame(dat, kernel_inputs$kernel_dat) str(df) ## ----------------------------------------------------------------------------- mod0 <- glm(counts ~ site + land1 + land2 + land3, family = poisson(), data = df) summary(mod0) ## ----opt1, eval=FALSE--------------------------------------------------------- # opt1 <- multiScale_optim(fitted_mod = mod0, # kernel_inputs = kernel_inputs) ## ----load-cached-opt1, include=FALSE------------------------------------------ opt1 <- vignette_cache$guide$opt1 ## ----opt1-screened-start, eval=FALSE------------------------------------------ # set.seed(1) # opt1_screen <- multiScale_optim( # fitted_mod = mod0, # kernel_inputs = kernel_inputs, # start_strategy = "screen", # screen_n_sigma = 4, # screen_n_jitter = 2, # screen_maxit = 5 # ) ## ----------------------------------------------------------------------------- summary(opt1) ## ----eval=FALSE--------------------------------------------------------------- # summary(opt1, profile = TRUE) ## ----opt2--------------------------------------------------------------------- ## New model mod0_2 <- glm(counts ~ site + land1 + land2, family = poisson(), data = df) ## ----eval=FALSE--------------------------------------------------------------- # ## Optimize # opt2 <- multiScale_optim(fitted_mod = mod0_2, # kernel_inputs = kernel_inputs) ## ----load-cached-opt2, include=FALSE------------------------------------------ opt2 <- vignette_cache$guide$opt2 ## ----------------------------------------------------------------------------- summary(opt2) ## ----eval=FALSE--------------------------------------------------------------- # summary(opt2, profile = TRUE) ## ----------------------------------------------------------------------------- ## Kernel function plot(opt2) ## Kernel function; 99% contribution plot(opt2, prob = 0.99) ## ----------------------------------------------------------------------------- plot_marginal_effects(opt2) ## ----sigma-profile, fig.cap = "AICc profile across sigma for each spatial covariate. The red dashed line marks the optimized sigma."---- n_profile_cores <- min(length(opt2$scale_est$Mean), max(1L, parallel::detectCores(logical = FALSE))) prof2 <- vignette_cache$guide$prof2 plot(prof2) ## ----kernel_raster------------------------------------------------------------ rast_opt <- terra::unwrap(vignette_cache$guide$rast_opt) plot(rast_opt) ## ----other_kernel, eval=FALSE------------------------------------------------- # ## Negative Exponential # exp_inputs <- kernel_prep(pts = pts, # raster_stack = land_rast, # max_D = 1700, # kernel = 'exp', # verbose = FALSE) # # ## Exponential Power # expow_inputs <- kernel_prep(pts = pts, # raster_stack = land_rast, # max_D = 1700, # kernel = 'expow', # verbose = FALSE) # # ## Fixed width buffer # fixed_inputs <- kernel_prep(pts = pts, # raster_stack = land_rast, # max_D = 1700, # kernel = 'fixed', # verbose = FALSE) ## ----eval=FALSE--------------------------------------------------------------- # opt_exp <- multiScale_optim(fitted_mod = mod0_2, # kernel_inputs = exp_inputs) # # ## Starting values needed # opt_expow <- multiScale_optim(fitted_mod = mod0_2, # kernel_inputs = expow_inputs, # par = c(500/expow_inputs$unit_conv, # 500/exp_inputs$unit_conv, # 5,10)) # # opt_fixed <- multiScale_optim(fitted_mod = mod0_2, # kernel_inputs = fixed_inputs) ## ----eval=FALSE--------------------------------------------------------------- # mod_list <- list(opt2, # opt_exp, # opt_expow, # opt_fixed) # # ## AIC table # aic_tab(mod_list) # # ## BIC table # bic_tab(mod_list) ## ----eval=FALSE--------------------------------------------------------------- # ## Landscape only effect # mod0_3 <- glm(counts ~ land1 + land2, # family = poisson(), # data = df) # # ## Landscape 1 only effect # mod0_4 <- glm(counts ~ land1, # family = poisson(), # data = df) # # ## Landscape 2 only effect # mod0_5 <- glm(counts ~ land2, # family = poisson(), # data = df) # # ## Landscape 3 only effect # mod0_6 <- glm(counts ~ land3, # family = poisson(), # data = df) # # ## Site only effect # ## No multiScaleR optimization # mod0_7 <- glm(counts ~ site, # family = poisson(), # data = df) ## ----eval=FALSE--------------------------------------------------------------- # opt3 <- multiScale_optim(fitted_mod = mod0_3, # kernel_inputs = kernel_inputs) # opt4 <- multiScale_optim(fitted_mod = mod0_4, # kernel_inputs = kernel_inputs) # opt5 <- multiScale_optim(fitted_mod = mod0_5, # kernel_inputs = kernel_inputs) # opt6 <- multiScale_optim(fitted_mod = mod0_6, # kernel_inputs = kernel_inputs) ## ----eval=FALSE--------------------------------------------------------------- # mod_list2 <- list(opt1, opt2, opt3, opt4, opt5, opt6, mod0_7) # # aic_tab(mod_list2) # # bic_tab(mod_list2) ## ----------------------------------------------------------------------------- rs <- sim_rast(user_seed = 999, dim = 250) plot(rs) ## ----------------------------------------------------------------------------- rs <- terra::subset(rs, c(1,4)) s_dat <- sim_dat_unmarked(alpha = 1, beta = c(0.75,-0.75), kernel = 'gaussian', sigma = c(75, 150), n_points = 75, n_surv = 5, det = 0.5, type = 'count', raster_stack = rs, max_D = 550, user_seed = 555) plot(s_dat$df$y ~ s_dat$df$bin1) plot(s_dat$df$y ~ s_dat$df$cont2) ## ----------------------------------------------------------------------------- library(unmarked) kernel_inputs <- kernel_prep(pts = s_dat$pts, raster_stack = rs, max_D = 550, kernel = 'gaus', verbose = FALSE) umf <- unmarkedFramePCount(y = s_dat$y, siteCovs = kernel_inputs$kernel_dat) ## Base unmarked model mod0_umf.p <- unmarked::pcount(~1 ~bin1 + cont2, data = umf, K = 100) ## ----------------------------------------------------------------------------- ## The pcount optimization refits a K = 100 N-mixture model at every optimizer ## step (~4 min), so a small precomputed result -- built from this same ## simulated data -- is loaded here to keep the vignette build fast. Every other ## model in this vignette is fit live. Regenerate the snapshot with ## data-raw/generate_pcount_vignette_object.R, or run it live with: ## opt_umf.p <- multiScale_optim(mod0_umf.p, kernel_inputs) opt_umf.p <- readRDS(system.file("extdata", "opt_umf_p.rds", package = "multiScaleR")) ## ----------------------------------------------------------------------------- summary(opt_umf.p) plogis(opt_umf.p$opt_mod@estimates@estimates$det@estimates[[1]]) ## ----------------------------------------------------------------------------- plot(opt_umf.p) plot_marginal_effects(opt_umf.p) ## ----------------------------------------------------------------------------- ## Project model rast_scale.center <- kernel_scale.raster(raster_stack = rs, multiScaleR = opt_umf.p, scale_center = TRUE, clamp = TRUE, pct_mx = 0.00) plot(rast_scale.center) ab.mod_pred <- terra::predict(rast_scale.center, opt_umf.p$opt_mod, type = 'state') plot(ab.mod_pred) ## ----eval=FALSE--------------------------------------------------------------- # s_dat.occ <- sim_dat_unmarked(alpha = 0.25, # beta = c(-0.75,1), # kernel = 'gaussian', # sigma = c(225, 75), # n_points = 250, # n_surv = 5, # det = 0.5, # type = 'occ', # raster_stack = rs, # max_D = 800, # user_seed = 555) # # plot(s_dat.occ$df$y ~ s_dat.occ$df$bin1) # plot(s_dat.occ$df$y ~ s_dat.occ$df$cont2) ## ----eval=FALSE--------------------------------------------------------------- # kernel_inputs <- kernel_prep(pts = s_dat.occ$pts, # raster_stack = rs, # max_D = 800, # kernel = 'gaus', # verbose = FALSE) # # ## Occupancy frame # umf <- unmarkedFrameOccu(y = s_dat.occ$y, # siteCovs = kernel_inputs$kernel_dat) # # ## Base unmarked model # (mod0_umf.occ <- unmarked::occu(~1 ~bin1 + cont2, # data = umf)) ## ----eval=FALSE--------------------------------------------------------------- # opt_umf.occ <- multiScale_optim(fitted_mod = mod0_umf.occ, # kernel_inputs = kernel_inputs) ## ----eval=FALSE--------------------------------------------------------------- # summary(opt_umf.occ) # plogis(opt_umf.occ$opt_mod@estimates@estimates$det@estimates[[1]]) ## ----eval=FALSE--------------------------------------------------------------- # ## Project model # rast_scale.center <- kernel_scale.raster(raster_stack = rs, # multiScaleR = opt_umf.occ, # scale_center = TRUE, # clamp = T) # plot(rast_scale.center) # # occ.mod_pred <- terra::predict(rast_scale.center, # opt_umf.occ$opt_mod, # type = 'state') # # plot(occ.mod_pred) ## ----eval=FALSE--------------------------------------------------------------- # set.seed(12345) # rs <- sim_rast(user_seed = 999, # dim = 500, resolution = 30) # rs <- terra::subset(rs, c(1,3)) # # ## Zero-inflated parameters # zi_alpha <- -0.25 # zi_beta <- -1 # # n_points <- 400 # alpha <- 0.5 # beta <- 0.75 # kernel <- 'gaussian' # sigma <- c(400, # For main effect, bin1 # 200) # For zero-inflated, cont1 # StDev <- 2 # Negative binomial dispersion # # ## Generate random points # bnd <- buffer(vect(ext(rs[[1]])), -1000) # pts <- spatSample(x = rs[[1]], # size = n_points, # method = 'random', # ext = ext(bnd), # replace = F, # as.points = T) # # kernel_out <- kernel_prep(pts = pts, # raster_stack = rs, # max_D = 1500, # kernel = kernel, # sigma = sigma, # verbose = FALSE) # # ## ZINB simulation # zi_prob <- plogis(zi_alpha + zi_beta*kernel_out$kernel_dat$cont1) # # # Simulate zero-inflated counts # # 1 = excess zero, 0 = from Poisson # is_zero <- rbinom(length(zi_prob), size = 1, prob = zi_prob) # # obs <- exp(alpha + beta*kernel_out$kernel_dat$bin1) # # obs.zinb <- ifelse(is_zero, 0, # rnbinom(length(is_zero), # mu = obs, # size = StDev)) # # dat <- data.frame(cnt = obs.zinb, # kernel_out$kernel_dat) # with(dat, plot(cnt ~ bin1)) ## ----eval=FALSE--------------------------------------------------------------- # library(pscl) # zinb_mod <- pscl::zeroinfl(cnt ~ bin1 | cont1, # dist = 'negbin', # data = dat) # # kernel_inputs <- kernel_prep(pts = pts, # kernel = kernel, # max_D = 1500, # raster_stack = rs, # verbose = FALSE) ## ----eval=FALSE--------------------------------------------------------------- # zinb_opt <- multiScale_optim(zinb_mod, # kernel_inputs) ## ----eval=FALSE--------------------------------------------------------------- # summary(zinb_opt) # plot(zinb_opt) # plot_marginal_effects(zinb_opt) ## ----eval=FALSE--------------------------------------------------------------- # ## With fitted model # kernel_dist(opt2) # # ## Distance of 95% kernel # kernel_dist(opt2, prob = 0.95) # ## ----------------------------------------------------------------------------- kernel_dist(kernel = "gaussian", sigma = 100, prob = 0.9) kernel_dist(kernel = "exp", sigma = 100, prob = 0.9) kernel_dist(kernel = "expow", sigma = 100, beta = 5, prob = 0.9) ## ----------------------------------------------------------------------------- plot_kernel(kernel = 'exp', sigma = 50) plot_kernel(kernel = 'expow', beta = 5, sigma = 50) plot_kernel(kernel = 'gaussian', sigma = 50) plot_kernel(kernel = 'gaussian', sigma = 50) plot_kernel(kernel = 'gaussian', sigma = 50, scale_dist = F) plot_kernel(kernel = 'gaussian', sigma = 50, add_label = F) ## ----error=FALSE, message=FALSE----------------------------------------------- r_sim1 <- sim_rast(dim = 100, resolution = 30, user_seed = 555) r_sim2 <- sim_rast(dim = 100, resolution = 30, autocorr_range1 = 100, autocorr_range2 = 1, sill = 10, user_seed = 555) plot(r_sim1) plot(r_sim2) ## ----eval=FALSE--------------------------------------------------------------- # s_dat <- sim_dat(alpha = 0.25, # beta = c(0.75, -0.75), # kernel = 'gaussian', # sigma = c(350, 200), # type = 'count', # n_points = 100, # raster_stack = terra::subset(r_sim1, c(1,4)), # min_D = 250, # user_seed = 999) ## ----eval=FALSE--------------------------------------------------------------- # plot(s_dat$df$y ~ s_dat$df$bin1) # plot(s_dat$df$y ~ s_dat$df$cont2) ## ----eval=FALSE--------------------------------------------------------------- # sim_mod <- glm(y ~ bin1 + cont2, # family = 'poisson', # data = s_dat$df) # # kernel_inputs <- kernel_prep(pts = s_dat$pts, # raster_stack = r_sim1, # max_D = 1500, # kernel = 'gaussian', # verbose = FALSE) ## ----eval=FALSE--------------------------------------------------------------- # sim_opt <- multiScale_optim(fitted_mod = sim_mod, # kernel_inputs = kernel_inputs) ## ----eval=FALSE--------------------------------------------------------------- # summary(sim_opt) # # plot(sim_opt) # # plot(kernel_scale.raster(raster_stack = r_sim1, # multiScaleR = sim_opt)) # # plot_marginal_effects(sim_opt)