关于terra::regress获取单元格斜率P值及扩展非参数方法的技术问询

1. 获取线性回归斜率的P值 #

terra::regress函数默认仅返回截距和斜率两个结果图层，不会直接输出斜率的P值。要获取每个单元格的P值，你可以通过自定义回归函数结合terra::app()来实现：

• 核心思路：针对每个单元格的时间序列值（即多层栅格对应位置的数值向量），调用lm()完成线性回归，从模型结果中提取截距、斜率以及斜率的P值。
• 示例代码：

library(terra)

# 构造示例多层SpatRaster（模拟10年气候数据，2x2单元格）
r <- rast(nrow=2, ncol=2, nlyr=10)
values(r) <- rnorm(ncell(r)*nlyr(r), mean=10, sd=2) + 0.5*rep(1:10, each=ncell(r))

# 自定义函数：输入单元格时间序列，输出截距、斜率、斜率P值
cell_lm <- function(x) {
  if (any(is.na(x))) return(c(NA, NA, NA))  # 处理缺失值
  model <- lm(x ~ seq_along(x))
  coefs <- coef(model)
  slope_p <- summary(model)$coefficients[2, 4]
  return(c(intercept=coefs[1], slope=coefs[2], slope_p=slope_p))
}

# 应用函数到每个单元格
result_lm <- app(r, cell_lm)
names(result_lm) <- c("intercept", "slope", "slope_p")

运行后result_lm会包含三个图层：截距、斜率、斜率的P值。

2. 扩展到非参数趋势指标（Theil-Sen、Mann-Kendall） #

完全可以扩展到非参数方法，同样通过terra::app()结合专业统计包实现：

（1）Theil-Sen斜率 #

可以使用zyp包的zyp.sen()函数，它能稳健估计趋势斜率：

library(zyp)

cell_sen <- function(x) {
  if (any(is.na(x))) return(NA)
  model <- zyp.sen(x ~ seq_along(x))
  return(coef(model)[2])  # 提取Theil-Sen斜率
}

result_sen <- app(r, cell_sen)
names(result_sen) <- "theil_sen_slope"

（2）Mann-Kendall趋势检验 #

使用trend包的mk.test()函数，提取趋势的显著性P值：

library(trend)

cell_mk <- function(x) {
  if (any(is.na(x))) return(NA)
  test <- mk.test(x)
  return(c(mk_stat=test$statistic, mk_p=test$p.value))
}

result_mk <- app(r, cell_mk)
names(result_mk) <- c("mk_statistic", "mk_p_value")

如果需要同时获取Theil-Sen斜率和Mann-Kendall检验结果，也可以合并到一个自定义函数中，一次性输出多个指标图层。

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