R语言GAM平滑NDVI数据集得到直线而非曲线的问题求助

问题描述 #

我正在对NDVI数据集进行平滑处理以提取每日NDVI值，此前处理另一数据集时获得了良好的曲线拟合效果，但处理新数据集时，GAM模型拟合结果为直线，不符合预期。

原始代码 #

library(janitor)

ndvifile2 <- read.csv("https://www.dropbox.com/s/1q13uym7h2medtq/adj_Final_NDVI_Mar10.txt?dl=1" , sep = "\t", header = T) %>% clean_names()

yf2 <- ndvifile2$ndvi_p50
xf2 <- ndvifile2$julian

# 可视化数据集
plot(xf2, yf2, type = "l")

# 将字符型转换为数值型
y.new2 <- as.numeric(yf2)

# 计算每个日期距初始日的天数
days.since2 <- as.numeric(xf2- min(xf2))

# 可视化转换后的变量
plot(days.since2,y.new2,typ="l")

# 高斯模型平滑
library(mgcv)
m2 <- gam(y.new2 ~ s(days.since2,k=40))
plot(days.since2,y.new2,typ="l")

# 在指定时间区间内预测
E.y2 <- predict(m2,newdata=data.frame(days.since2=0:10957))
points(0:10957,E.y2,col="gold",typ="l")

问题原因分析 #

1. 缺失值干扰：ndvi_p50列可能包含非数值字符，转换为numeric类型时生成大量NA，GAM模型默认忽略缺失值，若剩余有效数据点过少或分布无明显趋势，模型会拟合出直线。
2. 模型参数不匹配：指定的k=40可能远超有效数据的实际波动复杂度，导致模型无法捕捉到曲线趋势。
3. 预测范围超限：0:10957的预测范围远大于原始数据的时间跨度，超出部分模型只能用恒定值填充，呈现直线形态。

解决方案与修正代码 #

步骤1：处理缺失值并检查有效数据 #

先过滤缺失值，确认有效数据的分布情况：

# 过滤缺失值
valid_data <- data.frame(days = days.since2[!is.na(y.new2)], ndvi = y.new2[!is.na(y.new2)])
# 查看有效数据量
cat("有效数据点数量：", nrow(valid_data), "\n")
# 绘制有效数据时序图
plot(valid_data$days, valid_data$ndvi, type = "l", main = "有效NDVI时间序列")

步骤2：调整GAM模型参数 #

根据有效数据的波动情况调整平滑项的k值（建议设置为小于有效数据点数量的合理值），使用限制性样条（bs="cr"）提升拟合稳定性：

library(mgcv)
# 拟合优化后的GAM模型
m2 <- gam(ndvi ~ s(days, k = 30, bs = "cr"), data = valid_data)
# 查看模型摘要，确认平滑项显著性
summary(m2)

步骤3：在数据范围内预测 #

仅在原始有效数据的时间区间内预测，避免超出范围的无效延伸：

# 生成与原始数据匹配的预测时间序列
new_days <- seq(min(valid_data$days), max(valid_data$days), by = 1)
# 预测NDVI值
E.y2 <- predict(m2, newdata = data.frame(days = new_days))
# 绘制拟合结果
plot(valid_data$days, valid_data$ndvi, type = "l", main = "NDVI平滑拟合结果", xlab = "距初始日天数", ylab = "NDVI")
lines(new_days, E.y2, col = "gold", lwd = 2)

完整修正代码 #

library(janitor)
library(mgcv)

# 读取并清洗数据
ndvifile2 <- read.csv("https://www.dropbox.com/s/1q13uym7h2medtq/adj_Final_NDVI_Mar10.txt?dl=1" , sep = "\t", header = T) %>% clean_names()

# 变量转换与预处理
yf2 <- ndvifile2$ndvi_p50
xf2 <- ndvifile2$julian
y.new2 <- as.numeric(yf2)
days.since2 <- as.numeric(xf2 - min(xf2))

# 过滤缺失值
valid_data <- data.frame(days = days.since2[!is.na(y.new2)], ndvi = y.new2[!is.na(y.new2)])

# 检查有效数据
cat("有效数据点数量：", nrow(valid_data), "\n")
plot(valid_data$days, valid_data$ndvi, type = "l", main = "有效NDVI时间序列")

# 拟合GAM模型
m2 <- gam(ndvi ~ s(days, k = 30, bs = "cr"), data = valid_data)
summary(m2)

# 预测并可视化
new_days <- seq(min(valid_data$days), max(valid_data$days), by = 1)
E.y2 <- predict(m2, newdata = data.frame(days = new_days))

plot(valid_data$days, valid_data$ndvi, type = "l", main = "NDVI平滑拟合结果", xlab = "距初始日天数", ylab = "NDVI")
lines(new_days, E.y2, col = "gold", lwd = 2)

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