R语言adehabitatHR包kerneloverlapHR(PHR)异常值问题排查

问题核心 #

使用adehabitatHR包的kerneloverlapHR()函数，通过PHR方法计算多组动物家域重叠时，出现以下异常：

• 三组及以上个体时，结果矩阵的对角线值（i=j）严重偏离1（出现极小值或高达84的极端值）
• 重叠计算结果与可视化的家域分布完全不符（如模拟数据中无可见重叠，但计算显示27%的重叠率）
• 持续出现警告：In vi * aj : longer object length is not a multiple of shorter object length
• 尝试设置same4all=TRUE无改善，且核函数构建、家域面积估计均正常

PHR方法的定义为：动物j的利用分布（UD）在动物i家域内的积分，公式为 $\text{PHR}_{i,j} = \iint_i \text{UD}_j(x,y) dx dy$，结果应在0-1之间，i=j时近似为1。

问题原因 #

1. 网格分辨率不统一：kernelUD()默认会为每个个体生成独立的网格（基于各自的数据分布范围），当个体数据量或空间分布差异较大时，不同个体的UD网格维度、范围不一致，导致kerneloverlapHR()计算积分时出现维度不匹配，引发警告和错误的积分结果。
2. 函数调用大小写错误：原代码中调用的kerneloverlaphr()（小写HR）并非官方标准函数名，正确函数名为kerneloverlapHR()（大写HR），虽然R可能兼容小写调用，但易引发潜在逻辑问题。
3. 家域范围定义差异：PHR方法默认基于100%的UD范围计算积分，而可视化使用的是95%家域多边形，若UD尾部延伸超出95%范围，会导致计算值与可视化感知的重叠情况不符，但模拟数据中的极端值核心原因仍是网格不统一。

解决方案 #

1. 统一核密度计算网格 #

强制所有个体使用相同的空间网格范围和分辨率，确保UD矩阵维度一致，消除维度不匹配的警告，保证积分计算的准确性。

2. 使用标准函数名调用 #

替换kerneloverlaphr()为官方标准的kerneloverlapHR()函数。

3. 可选：自定义家域范围计算PHR #

若需基于95%（而非100%）家域计算重叠，可手动提取家域多边形后，计算UD在该范围内的积分。

修正后的验证代码 #

# 加载依赖包
library(tidyverse)
library(hms)
library(adehabitatHR)
library(sp)
library(sf)

# 模拟数据生成
date <- seq(from = as.Date("2023/01/01"), to = as.Date("2023/01/31"), by = "day")
time <- as_hms(seq(from = as_hms("06:00:00"), to = as_hms("17:00:00"), length.out = 500))
id <- seq(1, 3)

example <- data.frame(
  expand_grid(date, time, id)
) |> 
  arrange(id, date, time)

# 模拟家域坐标分布
set.seed(3)
lat_range <- list("1" = rnorm(n = 500, 
                            mean = 0.43, 
                            sd = 0.0038),
                  "2" = rnorm(n = 500,
                            mean = 0.40,
                            sd = 0.0043),
                  "3" = rnorm(n = 300,
                            mean = 0.39,
                            sd = 0.0040))
lon_range <- list("1" = rnorm(n = 500, 
                              mean = 23.01, 
                              sd = 0.0025),
                  "2" = rnorm(n = 500,
                              mean = 22.97,
                              sd = 0.0029),
                  "3" = rnorm(n = 500,
                              mean = 22.99,
                              sd = 0.0029))

example$lat <- numeric(46500)
example$lon <- numeric(46500)
for(id in seq(1, 3)){
  example$lat <- replace(example$lat, example$id == id, sample(lat_range[[id]], size = 15500, replace = TRUE))
  example$lon <- replace(example$lon, example$id == id, sample(lon_range[[id]], size = 15500, replace = TRUE))
}

# 转换为空间数据并投影
example.sp <- data.frame(id = example$id, x = example$lat, y = example$lon)
coordinates(example.sp) <- ~x+y
proj4string(example.sp) <- CRS( "+init=epsg:4326" )
example.sp <- spTransform(example.sp, CRS("+proj=utm +zone=34 +datum=WGS84 +units=m +no_defs"))

# 生成统一网格：基于所有个体的坐标范围设置固定分辨率
x_range <- range(example.sp@coords[,1])
y_range <- range(example.sp@coords[,2])
unified_grid <- expand.grid(x = seq(x_range[1], x_range[2], length.out = 50),
                            y = seq(y_range[1], y_range[2], length.out = 50))
coordinates(unified_grid) <- ~x+y
proj4string(unified_grid) <- proj4string(example.sp)

# 使用统一网格构建核密度
example.kernelref <- kernelUD(example.sp, h = "href", grid = unified_grid)
example.kernel.poly <- getverticeshr(example.kernelref, percent = 95, unin = "m", unout = "km2") 

# 查看家域面积
example.kernel.poly$area

# 使用标准函数计算PHR重叠
overlap_example <- kerneloverlapHR(example.kernelref, method = "PHR")
print(overlap_example)

自定义95%家域PHR计算（可选） #

若需严格基于95%家域范围计算重叠，可手动实现积分：

# 初始化PHR矩阵
phr_95 <- matrix(nrow = 3, ncol = 3, dimnames = list(paste0("id",1:3), paste0("id",1:3)))
# 计算网格单元面积
cell_area <- (x_range[2]-x_range[1])/49 * (y_range[2]-y_range[1])/49

for(i in 1:3){
  # 获取个体i的95%家域多边形
  poly_i <- example.kernel.poly[i,]
  # 遍历所有个体j
  for(j in 1:3){
    # 提取个体j的UD数据
    ud_j <- example.kernelref[[j]]$ud
    # 筛选出家域i内的UD网格点
    pts_in_i <- over(ud_j, poly_i)
    ud_in_i <- ud_j[!is.na(pts_in_i),]
    # 计算积分：UD值总和乘以网格单元面积
    phr_95[i,j] <- sum(ud_in_i$ud) * cell_area
  }
}
print(phr_95)

效果验证 #

• 统一网格后，维度不匹配的警告消失
• 结果矩阵对角线值接近1（因数值精度略有偏差，但无极端值）
• 模拟数据中无重叠的个体，PHR值趋近于0，与可视化结果一致

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