如何在R包amt中模拟稳态利用分布（替代旧函数方案）

问题背景 #

尝试复现Signer等人2019年论文的代码时，发现当前amt 0.2.1.0版本已移除habitat_kernel()、movement_kernel()和simulate_ud()函数。除了从旧版本包提取这些函数外，可通过以下两种方法实现结合移动核与生境核的稳态利用分布模拟，最终生成预测图。

替代方案1：手动实现核心模拟逻辑 #

核心思路 #

稳态利用分布是生境偏好权重与动物移动扩散过程的结合，手动实现分为三步：计算生境栅格的偏好权重、基于步长/转角分布模拟移动路径、统计栅格访问频率得到分布。

代码实现 #

library(amt)
library(terra)
library(circular)

# 加载数据
data("deer")
data("sh_forest")

# 预处理数据（与原代码一致）
data_ssf <- deer |>
  steps_by_burst() |>
  random_steps(n_control = 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_))

# 拟合SSF模型
m1 <- data_ssf |> 
  fit_ssf(case_ ~ forest*cos_ta + forest*log_sl + strata(step_id_))

# 提取关键参数
# 1. 步长Weibull分布参数
sl_params <- sl_distr(m1)$params
shape <- sl_params$shape
scale <- sl_params$scale

# 2. 转角von Mises分布浓度参数（基于模型系数推导）
ta_coef <- coef(m1)[grepl("cos_ta", names(coef(m1)))]
kappa <- exp(sum(ta_coef))

# 3. 生境偏好系数
forest_coef <- coef(m1)["forestforest"]

# 计算生境栅格权重并标准化
sh_forest$habitat_weight <- exp(forest_coef * (sh_forest$forest == 1))
sh_forest$habitat_weight <- sh_forest$habitat_weight / sum(sh_forest$habitat_weight, na.rm = TRUE)

# 转换为terra栅格方便空间操作
habitat_r <- rast(sh_forest)[["habitat_weight"]]
count_r <- habitat_r
values(count_r) <- 0

# 模拟移动路径生成稳态UD
set.seed(123)
n_steps <- 1e6 # 模拟步数，可根据精度需求调整
current_point <- as.numeric(data_ssf[1, c("x1_", "y1_")])

for (i in 1:n_steps) {
  # 模拟步长与转角
  sl <- rweibull(1, shape = shape, scale = scale)
  ta <- rvonmises(1, mu = 0, kappa = kappa)
  
  # 计算下一个点坐标
  next_x <- current_point[1] + sl * cos(ta)
  next_y <- current_point[2] + sl * sin(ta)
  
  # 检查点是否在栅格范围内，并用Metropolis-Hastings采样接受/拒绝
  if (!is.na(extract(habitat_r, c(next_x, next_y)))) {
    current_weight <- extract(habitat_r, current_point)
    next_weight <- extract(habitat_r, c(next_x, next_y))
    accept_prob <- min(1, next_weight / current_weight)
    
    if (runif(1) < accept_prob) {
      current_point <- c(next_x, next_y)
    }
    
    # 记录当前点访问计数
    cell_idx <- cellFromXY(count_r, current_point)
    count_r[cell_idx] <- count_r[cell_idx] + 1
  }
}

# 标准化计数得到稳态利用分布
ssud <- count_r / sum(values(count_r), na.rm = TRUE)

# 绘制预测图
plot(ssud, main = "稳态利用分布（结合移动核与生境核）")
plot(rast(sh_forest)[["forest"]], add = TRUE, alpha = 0.3)

替代方案2：利用amt预测函数+空间卷积 #

核心思路 #

先用amt的predict()生成生境偏好栅格，再结合移动扩散核进行卷积操作，模拟移动对生境利用的扩散效应，效率更高适合大尺度区域。

代码实现 #

library(amt)
library(terra)
library(stats)

# 加载并预处理数据（同前）
data("deer")
data("sh_forest")
data_ssf <- deer |>
  steps_by_burst() |>
  random_steps(n_control = 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_))
m1 <- data_ssf |> fit_ssf(case_ ~ forest*cos_ta + forest*log_sl + strata(step_id_))

# 1. 生成生境偏好预测栅格
pred_grid <- sh_forest |> 
  mutate(forest = factor(forest, levels = 1:0, labels = c("forest", "non-forest")),
         cos_ta = 0, log_sl = 0) # 固定转角与步长的均值

# 预测生境偏好相对概率
pred_ssf <- predict(m1, newdata = pred_grid, type = "response")
pred_r <- rast(pred_ssf, type = "xyz")
names(pred_r) <- "habitat_preference"

# 2. 构建移动扩散核
sl_params <- sl_distr(m1)$params
# 生成步长与转角样本，转换为坐标偏移
sl_vals <- rweibull(1000, shape = sl_params$shape, scale = sl_params$scale)
ta_vals <- rvonmises(1000, mu = 0, kappa = exp(coef(m1)[grepl("cos_ta", names(coef(m1)))]))
dx <- sl_vals * cos(ta_vals)
dy <- sl_vals * sin(ta_vals)

# 创建移动核栅格
kernel_size <- ceiling(max(abs(dx), abs(dy)))
kernel_grid <- expand.grid(x = -kernel_size:kernel_size, y = -kernel_size:kernel_size)
kernel_grid$weight <- dnorm(sqrt(kernel_grid$x^2 + kernel_grid$y^2), mean = 0, sd = sd(sl_vals))
kernel_r <- rast(kernel_grid, type = "xyz")

# 3. 卷积生境偏好与移动核，得到考虑移动的利用分布
ssud <- terra::convolve(pred_r, kernel_r, normalize = TRUE)

# 绘制预测图
plot(ssud, main = "稳态利用分布（卷积生境偏好与移动核）")
plot(rast(sh_forest)[["forest"]], add = TRUE, alpha = 0.3)

关键注意事项 #

• 模拟步数或核样本量越大，稳态分布精度越高，但计算耗时越长；
• 转角分布参数可根据模型交互项细化，比如不同生境类型下设置不同的浓度参数；
• 手动模拟更贴近论文原始逻辑，卷积方法更适合大尺度研究场景。

内容的提问来源于stack exchange，提问作者kflorko