在函数中使用nls()拟合指数衰减曲线时遇str2lang错误求助

尝试在函数中使用nls()对生物数据拟合指数衰减曲线（需多次执行该操作），但遇到如下错误：

Error in str2lang(x) : <text>:2:0: unexpected end of input
1: ~ 
   ^

测试数据集及代码如下：

dat_small <-  structure(list(actual_time = c(5, 9, 30, 59, 119, 171, 216), 
                             activity = c(7158, 7386, 5496, 3884, 1502, 819, 409)), row.names = c(NA, 
                                                                                                  -7L), class = c("tbl_df", "tbl", "data.frame"))
curve_fit <-  function(x, y, df) {
  #browser()
  # Mono-exponential model for activity using actual sampling time
  mod <-  lm(log(df[[y]]) ~ df[[x]], data = df) # get starting values from log-linear model
  C0 <-  as.numeric(exp(mod$coef[1])) # exponential of log-linear intercept
  lambda1 <-  as.numeric(abs(mod$coef[2])) # absolute value of slope of log-linear model
  nls_mod_mono <-  nls(df[[y]] ~ C0*exp(-lambda1*df[[x]]), start = c(C0 = C0, lambda1 = lambda1), data = df)
  summary(nls_mod_mono) # estimate model
  xNew <- seq(0, 240, length.out = 100) # new grid of times
  yNew <- predict(nls_mod_mono, list(x = xNew)) # predict activity
  dfNew <-  data.frame(x = xNew, y = yNew)
  p_mono <-  ggplot(dfNew, aes(x = x, y = y)) +
    geom_line() +
    geom_point(data = df, aes(x = x, y = y), size = 3) +
    xlab("Actual Sampling Time (mins)") + ylab("Activity (Bq)") +
    theme_bw(base_size = 20)
  mono_half_life <-  0.693/as.numeric(coef(nls_mod_mono)[2]) # half-life
  mono_auc <-  as.numeric(coef(nls_mod_mono)[1])/as.numeric(coef(nls_mod_mono)[2]) # AUC to infinity (can be estimated in mono case by simply dividing the intercept by lambda)
  results <-  c(p_mono, mono_half_life, mono_auc)
}
curve_fit("actual_time", "activity", dat_small)

错误核心原因是nls()无法解析直接用df[[x]]/df[[y]]写的公式，需要动态构建公式并修正几处代码逻辑：

关键修正点 #

1. 动态构建公式：用reformulate()和as.formula()根据传入的列名参数生成模型公式，让lm()和nls()能正确识别数据变量
2. 统一预测变量名：预测时使用与原数据列名一致的变量名，确保predict()匹配模型参数
3. 规范结果返回：用列表打包ggplot对象和数值结果，避免混合类型导致的混乱
4. ggplot动态变量引用：用.data[[x]]语法在ggplot中动态调用列名，符合tidyeval规范

修正后的完整代码 #

library(ggplot2)

dat_small <-  structure(list(actual_time = c(5, 9, 30, 59, 119, 171, 216), 
                             activity = c(7158, 7386, 5496, 3884, 1502, 819, 409)), row.names = c(NA, 
                                                                                                  -7L), class = c("tbl_df", "tbl", "data.frame"))

curve_fit <-  function(x, y, df) {
  # 动态构建log线性模型公式
  lm_formula <- reformulate(x, response = paste0("log(", y, ")"))
  mod <- lm(lm_formula, data = df)
  
  C0 <- as.numeric(exp(mod$coef[1]))
  lambda1 <- as.numeric(abs(mod$coef[2]))
  
  # 动态构建nls模型公式
  nls_formula <- as.formula(paste0(y, " ~ C0*exp(-lambda1*", x, ")"))
  nls_mod_mono <- nls(nls_formula, start = list(C0 = C0, lambda1 = lambda1), data = df)
  
  # 生成预测数据（变量名与原数据一致）
  xNew <- seq(0, 240, length.out = 100)
  pred_data <- data.frame(!!x := xNew)
  yNew <- predict(nls_mod_mono, newdata = pred_data)
  
  dfNew <- data.frame(x = xNew, y = yNew)
  p_mono <- ggplot(dfNew, aes(x = x, y = y)) +
    geom_line() +
    geom_point(data = df, aes(x = .data[[x]], y = .data[[y]]), size = 3) +
    xlab("Actual Sampling Time (mins)") + ylab("Activity (Bq)") +
    theme_bw(base_size = 20)
  
  mono_half_life <- 0.693/as.numeric(coef(nls_mod_mono)[2])
  mono_auc <- as.numeric(coef(nls_mod_mono)[1])/as.numeric(coef(nls_mod_mono)[2])
  
  # 用列表返回不同类型结果
  list(plot = p_mono, half_life = mono_half_life, auc = mono_auc)
}

# 运行函数并查看结果
result <- curve_fit("actual_time", "activity", dat_small)
print(result$plot)
cat("半衰期：", round(result$half_life, 2), "分钟\n")
cat("AUC：", round(result$auc, 2), "\n")

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