基于terra包为时序栅格单元格拟合自定义回归模型的问题

问题回顾 #

需要基于terra包为时序栅格的每个单元格拟合回归模型，尝试用regress函数时出现斜率全为0的问题，示例代码如下：

# 创建示例时序栅格
d<-rast(lapply(1:10, FUN = function(x){
  rast(matrix(rnorm(100, mean=x), ncol=10))
}))
names(d)<-1:10

# 创建年份栅格（错误写法）
d1<-rast(lapply(names(d)[1:10], function(x) {rast(matrix(x, nrow=10, ncol=10))
}))
names(d1)<-1:10

# 执行回归后斜率全为0
plot(regress(d, d1))

错误原因 #

你创建的d1栅格中，年份是字符类型（names(d)返回字符串），而线性回归要求自变量为数值型。字符型自变量无法参与回归计算，导致斜率全部为0。

解决方案 #

1. 修正年份栅格的数值类型 #

将年份转为数值后再生成栅格：

d1 <- rast(lapply(as.numeric(names(d)), function(x) {
  rast(matrix(x, nrow=10, ncol=10))
}))
names(d1) <- names(d)

2. 用terra::regress执行标准线性回归 #

重新运行regress即可得到正确的斜率结果：

lm_result <- regress(d, d1, fun=lm)
plot(lm_result[["x"]]) # 查看斜率栅格

3. 适配mblm包的中位数回归 #

要拟合mblm的自定义模型，需编写一个返回模型参数的函数，供regress调用：

library(mblm)

# 自定义函数：输入因变量y（栅格时间序列）和自变量x（年份），返回模型系数
mblm_fit <- function(y, x) {
  # 过滤可能的NA值（可选，根据数据情况调整）
  valid_idx <- !is.na(y) & !is.na(x)
  if(sum(valid_idx) < 2) return(c(NA, NA))
  
  model <- mblm(y[valid_idx] ~ x[valid_idx])
  coef(model) # 返回截距和斜率两个参数
}

# 执行自定义回归
mblm_result <- regress(d, d1, fun=mblm_fit)
names(mblm_result) <- c("Intercept", "Slope")
plot(mblm_result[["Slope"]]) # 查看中位数回归的斜率栅格

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