# \donttest{
# RSF -------------------------------------------------------
# Fit an RSF, then calculate log-RSS to visualize results.
# Load packages
library(ggplot2)
# Load data
data("amt_fisher")
amt_fisher_covar <- get_amt_fisher_covars()
# Prepare data for RSF
rsf_data <- amt_fisher |>
filter(name == "Lupe") |>
make_track(x_, y_, t_) |>
random_points() |>
extract_covariates(amt_fisher_covar$elevation) |>
extract_covariates(amt_fisher_covar$popden) |>
extract_covariates(amt_fisher_covar$landuse) |>
mutate(lu = factor(landuse))
# Fit RSF
m1 <- rsf_data |>
fit_rsf(case_ ~ lu + elevation + popden)
# Calculate log-RSS
# data.frame of x1s
x1 <- data.frame(lu = factor(50, levels = levels(rsf_data$lu)),
elevation = seq(90, 120, length.out = 100),
popden = mean(rsf_data$popden))
# data.frame of x2 (note factor levels should be same as model data)
x2 <- data.frame(lu = factor(50, levels = levels(rsf_data$lu)),
elevation = mean(rsf_data$elevation),
popden = mean(rsf_data$popden))
# Calculate (use se method for confidence interval)
logRSS <- log_rss(object = m1, x1 = x1, x2 = x2, ci = "se")
# Plot
ggplot(logRSS$df, aes(x = elevation_x1, y = log_rss)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "gray80") +
geom_line() +
xlab(expression("Elevation " * (x[1]))) +
ylab("log-RSS") +
ggtitle(expression("log-RSS" * (x[1] * ", " * x[2]))) +
theme_bw()
# SSF -------------------------------------------------------
# Fit an SSF, then calculate log-RSS to visualize results.
# Load data
data(deer)
sh_forest <- get_sh_forest()
# Prepare data for SSF
ssf_data <- deer |>
steps_by_burst() |>
random_steps(n = 15) |>
extract_covariates(sh_forest) |>
mutate(forest = factor(forest, levels = 1:0,
labels = c("forest", "non-forest")),
cos_ta = cos(ta_),
log_sl = log(sl_))
# Fit an SSF (note model = TRUE necessary for predict() to work)
m2 <- ssf_data |>
fit_clogit(case_ ~ forest + strata(step_id_), model = TRUE)
# Calculate log-RSS
# data.frame of x1s
x1 <- data.frame(forest = factor(c("forest", "non-forest")))
# data.frame of x2
x2 <- data.frame(forest = factor("forest", levels = levels(ssf_data$forest)))
# Calculate
logRSS <- log_rss(object = m2, x1 = x1, x2 = x2, ci = "se")
# Plot
ggplot(logRSS$df, aes(x = forest_x1, y = log_rss)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0.25) +
geom_point(size = 3) +
xlab(expression("Forest Cover " * (x[1]))) +
ylab("log-RSS") +
ggtitle(expression("log-RSS" * (x[1] * ", " * x[2]))) +
theme_bw()
# }
Run the code above in your browser using DataLab