在R循环中调用ns()拟合lme模型时出现变量长度不匹配错误

问题重现 #

运行以下循环代码，尝试用ns(x, df=j)拟合带随机效应的模型时触发错误：

library(splines)
library(nlme)

# 生成数据
x <- 1:100
xid <- sort(rep.int(1:10, 10))
y <- x + rnorm(10)

# 循环遍历样条节点数
for (j in 1:5) {
  # 拟合j个节点的样条模型
  fit <- lme(y ~ ns(x, j), random = ~1|xid)
  # 打印模型摘要
  print(summary(fit))
}

得到错误：

Error in model.frame.default(formula = ~y + x + j, data = <environment>,  : variable lengths differ (found for 'j')

注意：使用无分组的lm模型时运行正常，直接将j替换为1-5的具体整数时也能正常运行。

问题根源 #

lme的公式解析机制与lm存在差异：

• lm会优先识别公式中ns(x,j)里的j为当前环境中的迭代变量（长度为1的标量），不会将其当作数据变量处理。
• lme在解析公式时，会把j当作需要从数据环境中读取的变量，而循环中的j是长度为1的标量，与y、x等长度为100的向量长度不匹配，因此触发模型框构建的长度错误。

解决方案 #

方案1：动态构建公式 #

用as.formula()将j的当前值嵌入到公式字符串中，让lme直接解析出固定的节点数：

library(splines)
library(nlme)

x <- 1:100
xid <- sort(rep.int(1:10, 10))
y <- x + rnorm(10)

for (j in 1:5) {
  # 动态生成公式
  formula_str <- paste0("y ~ ns(x, ", j, ")")
  fit <- lme(as.formula(formula_str), random = ~1|xid)
  print(summary(fit))
}

方案2：预计算样条矩阵 #

先在循环中计算好样条基矩阵，再将其作为变量代入模型，避免公式解析时的歧义：

library(splines)
library(nlme)

x <- 1:100
xid <- sort(rep.int(1:10, 10))
y <- x + rnorm(10)

for (j in 1:5) {
  # 预计算样条基矩阵
  spline_mat <- ns(x, j)
  # 将矩阵合并到数据框
  model_data <- data.frame(y, spline_mat, xid)
  # 拟合模型
  fit <- lme(y ~ ., data = model_data, random = ~1|xid)
  print(summary(fit))
}

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