## ----setup, include=FALSE, echo=FALSE----------------------------------------- library(multiScaleR) # Define the base URL for the extdata directory base_url <- "https://raw.githubusercontent.com/wpeterman/multiScaleR/master/inst/extdata/" f <- c("opt_exp.RData", "opt_expow.RData", "opt_fixed.RData", "opt_umf_counts.RData", "opt_umf_occ.RData", "opt1.RData", "opt2.RData", "opt3.RData", "opt4.RData", "opt5.RData", "opt6.RData", "sim_opt.RData", "zinb.RData") # Download each file for (file in f) { download.file( url = paste0(base_url, file), destfile = file.path(tempdir(), file), mode = "wb" ) } rdat <- list.files(tempdir(), pattern = "*.RData", full.names = T) invisible(lapply(rdat, load, .GlobalEnv)) ## ----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 ## ----------------------------------------------------------------------------- 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) ## ----------------------------------------------------------------------------- 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) ## ----------------------------------------------------------------------------- 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."---- prof2 <- profile_sigma(opt2, n_pts = 15, verbose = FALSE) plot(prof2) ## ----kernel_raster------------------------------------------------------------ rast_opt <- kernel_scale.raster(raster_stack = land_rast, multiScaleR = opt2) plot(rast_opt) ## ----other_kernel, eval=FALSE------------------------------------------------- # ## Negative Exponential # ## Note: these kernel_prep objects are pre-loaded in the setup chunk. # ## eval=FALSE keeps the vignette from re-running these computationally # ## intensive calls during knitting. # 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) ## ----------------------------------------------------------------------------- mod_list <- list(opt2, opt_exp, opt_expow, opt_fixed) ## AIC table aic_tab(mod_list) ## BIC table bic_tab(mod_list) ## ----------------------------------------------------------------------------- ## 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) ## ----------------------------------------------------------------------------- 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) ## ----eval=FALSE--------------------------------------------------------------- # opt_umf.p <- multiScale_optim(fitted_mod = mod0_umf.p, # kernel_inputs = kernel_inputs, # n_cores = 8) ## ----------------------------------------------------------------------------- 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) ## ----------------------------------------------------------------------------- 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) ## ----------------------------------------------------------------------------- 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, # n_cores = 8) ## ----------------------------------------------------------------------------- summary(opt_umf.occ) plogis(opt_umf.occ$opt_mod@estimates@estimates$det@estimates[[1]]) ## ----------------------------------------------------------------------------- ## 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) ## ----------------------------------------------------------------------------- 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)) ## ----------------------------------------------------------------------------- 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, # n_cores = 4) ## ----------------------------------------------------------------------------- summary(zinb_opt) plot(zinb_opt) plot_marginal_effects(zinb_opt) ## ----------------------------------------------------------------------------- ## 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) ## ----------------------------------------------------------------------------- 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) ## ----------------------------------------------------------------------------- plot(s_dat$df$y ~ s_dat$df$bin1) plot(s_dat$df$y ~ s_dat$df$cont2) ## ----------------------------------------------------------------------------- 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) ## ----------------------------------------------------------------------------- summary(sim_opt) plot(sim_opt) plot(kernel_scale.raster(raster_stack = r_sim1, multiScaleR = sim_opt)) plot_marginal_effects(sim_opt)