如何在R中用PLS模型平均系数对栅格图像进行预测？

问题背景 #

我有含100个观测值的数据集（50个X变量+1个Y变量），已完成100次置换划分（60-40训练/测试集）并训练PLS模型，得到100组含截距的系数。需要将这些系数的平均值应用到栅格图像，生成预测的Y值栅格。

完整实现代码与步骤 #

library(pls)
library(raster)

# 1. 加载数据并初始化系数存储
my_data <- read.csv("file1.csv", header = TRUE)
# 初始化矩阵：100次迭代 × 51个系数（50个X变量+截距）
output.coefs <- matrix(NA, nrow = 100, ncol = 51)

# 2. 100次置换训练PLS并保存系数
for (i in 1:100) {
  # 固定样本量为整数（避免小数问题）
  index <- sample(1:nrow(my_data), size = floor(0.6*nrow(my_data)))
  cal_data <- my_data[index,]
  
  # 训练PLS模型
  pls_model <- plsr(Y~., data = cal_data, ncomp = 10, validation = "LOO")
  
  # 提取含截距的系数并存储
  output.coefs[i,] <- as.vector(coef(pls_model, ncomp = 10, intercept = TRUE))
}

# 3. 计算平均系数并保存
mean_coefs <- colMeans(output.coefs, na.rm = TRUE)
# 拆分截距和X变量系数
intercept_val <- mean_coefs[1]
x_var_coefs <- mean_coefs[-1]
# 保存平均系数到CSV
write.csv(mean_coefs, "mean_pls_coefs.csv", row.names = FALSE)

# 4. 加载栅格并验证匹配性
# 注意：需加载50波段的栅格（对应50个X变量），用stack而非raster
my_raster_stack <- stack("example_raster1")

# 检查波段数与X变量数量是否一致
if(nlayers(my_raster_stack) != length(x_var_coefs)){
  stop("栅格波段数与X变量数量不匹配，请检查数据！")
}

# 5. 应用系数生成预测Y栅格
# 计算公式：Y = 截距 + Σ(栅格波段i × 对应X变量系数i)
predicted_y_raster <- intercept_val + sum(my_raster_stack * x_var_coefs)

# 6. 保存预测结果
writeRaster(predicted_y_raster, "predicted_Y_raster.tif", overwrite = TRUE)

关键说明 #

• 系数存储初始化：原代码未初始化output.coefs，需提前创建矩阵避免赋值报错
• 样本量处理：用floor()确保训练集样本数为整数，避免sample()函数报错
• 栅格加载：必须用stack()加载多波段栅格（50个波段对应50个X变量），单波段栅格不符合需求
• 系数匹配：务必保证栅格波段顺序与数据集中X变量的顺序完全一致，否则结果会出错
• 栅格运算：利用raster包的向量化运算，高效完成多波段与系数的线性组合

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