R语言For循环报错：Error in qr.default(t(const))问题修复请求

问题根源分析 #

报错Error in qr.default(t(const)) : NA/NaN/Inf in foreign function call (arg 1)来自三个核心问题：

1. 无效节点组合：你生成了包含仅1个节点的组合，此时边界节点会被设为相同值，导致样条基矩阵线性相关，触发奇异矩阵错误。
2. 权重中的异常值：data1[, "Ec_x"] / data1[, "mu_x"]计算权重时，若mu_x（死亡率）为0会产生Inf，传入lm()后引发矩阵运算失败。
3. 冗余组合浪费资源：原逻辑生成大量无效组合，既拖慢计算又增加出错概率。

修复步骤与完整代码 #

1. 过滤无效节点组合 #

只保留至少2个节点的组合，确保边界节点有效且覆盖合理范围：

# 生成至少2个节点的组合，排除单节点情况
knot_combinations <- unlist(lapply(2:length(the.knots), combn, x = the.knots, simplify = FALSE), recursive = FALSE)
hyper_param <- expand.grid(knots = knot_combinations, Error.Rate = 0)

2. 处理权重异常值 #

替换权重中的Inf和NA为合理值，避免破坏线性回归计算：

weights_vec <- data1[, "Ec_x"] / data1[, "mu_x"]
# 替换Inf为最大有限权重，NA为权重均值
weights_vec[is.infinite(weights_vec)] <- max(weights_vec[is.finite(weights_vec)], na.rm = TRUE)
weights_vec[is.na(weights_vec)] <- mean(weights_vec[is.finite(weights_vec)], na.rm = TRUE)

3. 添加错误捕获机制 #

用tryCatch包裹拟合步骤，避免单个组合出错导致整个循环中断：

完整修复代码 #

library(demography)
library(splines)

# 从人类死亡率数据库读取日本死亡率数据
JPNmort <- hmd.mx("JPN","username","password")

age = 2:110

data1= cbind(1:109, 
             JPNmort[["rate"]][["total"]][age, "2016"], 
             JPNmort[["pop"]][["total"]][age, "2016"])
data2 = cbind(1:109, 
              JPNmort[["rate"]][["total"]][age, "2017"],
              JPNmort[["pop"]][["total"]][age, "2017"])

# 格式化数据行和列
column.names = c("Age", "mu_x", "Ec_x")
rownames(data1) = NULL
rownames(data2) = NULL

colnames(data1) = column.names
colnames(data2) = column.names

the.knots <- c(10,20,30,40,50,60,70,80,90,100)
# 过滤无效组合：仅保留至少2个节点的情况
knot_combinations <- unlist(lapply(2:length(the.knots), combn, x = the.knots, simplify = FALSE), recursive = FALSE)
hyper_param <- expand.grid(knots = knot_combinations, Error.Rate = 0)

# 初始化存储结果的列表
my.basis <- vector("list", nrow(hyper_param))
my.spline <- vector("list", nrow(hyper_param))
# 获取数据年龄范围，限制边界节点范围
age_range <- range(data1[, "Age"])

for (i in 1:nrow(hyper_param)) {
  my.knots <- as.numeric(unlist(hyper_param[i, 1]))
  
  # 生成样条基，确保边界节点在数据年龄范围内
  my.basis[[i]] <- ns(data1[, "Age"], 
                      knots = setdiff(my.knots, c(my.knots[1], my.knots[length(my.knots)])), 
                      Boundary.knots = pmax(pmin(c(my.knots[1], my.knots[length(my.knots)]), age_range[2]), age_range[1])
  )
  
  # 处理权重中的Inf/NA
  weights_vec <- data1[, "Ec_x"] / data1[, "mu_x"]
  weights_vec[is.infinite(weights_vec)] <- max(weights_vec[is.finite(weights_vec)], na.rm = TRUE)
  weights_vec[is.na(weights_vec)] <- mean(weights_vec[is.finite(weights_vec)], na.rm = TRUE)
  
  # 尝试拟合模型，捕获错误避免循环中断
  fit_result <- tryCatch({
    lm(data1[, "mu_x"] ~ my.basis[[i]], weights = weights_vec)
  }, error = function(e) {
    message(paste("Skipping combination", i, ":", e$message))
    return(NULL)
  })
  
  my.spline[[i]] <- fit_result
  
  # 仅在拟合成功时计算误差
  if (!is.null(fit_result)) {
    hyper_param$Error.Rate[i] <- mean((data2[, "mu_x"] - fitted(fit_result))^2, na.rm = TRUE)
  } else {
    hyper_param$Error.Rate[i] <- NA
  }
}

关键修复总结 #

• 排除单节点组合，避免边界节点重复引发的奇异矩阵问题。
• 清理权重中的异常值，确保线性回归输入有效。
• 加入错误捕获，防止单个组合的问题中断整个循环。
• 限制边界节点在数据年龄范围内，避免样条外推异常。

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