主页 / coding / 问题 / 79366909
Accepted
LucaS
Asked: 2025-01-18 17:54:01 +0800 CST

尝试在函数内使用 nls() 时，str2lang(x) 出现错误

  • 772

我正在尝试使用函数将指数衰减曲线拟合到一些生物数据nls()（我想多次执行此操作）。

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)

2 个回答

  1. 在这种情况下，尝试formula使用。nlsactivity ~ C0 * exp(-lambda1 * actual_time)

    您可以使用as.formula(paste(y, "~ C0*exp(-lambda1*", x, ")", sep = ""))、 或 ，sprintf如下所示。

    要使用中指定的名称来命名参数list中的变量，请使用。newdata=predictxsetNames

    请注意，使用df[[x]]是可能的，因为 ist 仍在使用lm。但此时该data=参数是多余的。

    可能as.numeric会造成混淆，因为实际上已经是了"numeric"。请unname改用。

    curve_fit <- function(x, y, df, plot=FALSE) {
 ## 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 <- unname(exp(mod$coef[1]))  ## exponential of log-linear intercept
 lambda1 <- unname(abs(mod$coef[2]))  ## absolute value of slope of log-linear model
 fo <- as.formula(sprintf('%s ~ C0*exp(-lambda1*%s)', y, x))
 nls_mod_mono <- nls(fo, start=c(C0=C0, lambda1=lambda1), data=df)
 xNew <- seq(0, 240, length.out=100)  ## new grid of times
 yNew <- predict(nls_mod_mono, newdata= list(xNew) |> setNames(x))  ## predict activity
 fit_df <- data.frame(x=xNew, y=yNew)
 mono_half_life <- 0.693/coef(nls_mod_mono)[2]  ## half-life
 mono_auc <- coef(nls_mod_mono)[1]/coef(nls_mod_mono)[2]  ## AUC to infinity (can be estimated in mono case by simply dividing the intercept by lambda)
 list(fit_df=fit_df, 
      mono_half_life=unname(mono_half_life), 
      mono_auc=unname(mono_auc))
}

> res <- curve_fit("actual_time", "activity", dat_small)
> res$mono_half_life
[1] 53.25588
> res$mono_auc
[1] 615006.6

    我建议在函数外部绘制曲线以避免特征蔓延。

    par(mar=c(5, 4, 2, 2))
plot(y ~ x, res$fit_df, type='l', xlab='Actual Sampling Time (mins)',
     ylab="Activity (Bq)", las=1)

    在此处输入图片描述

    • 4
  2. Best Answer

    对 SO 的问题应该包含最少的代码，重点放在手头的问题上。这里的问题是如何避免调用中的错误nls，因此之后的所有内容都不相关，我们将其省略。（但是，我们确实指出，aes在 ggplot2 中接受未加引号的变量，并且由于x和y是字符，因此请像这样写aes(.data[[x]], .data[[y]])。）

    1）一个简单的解决方法可能只是通过定义具有固定 x 和 y 列名的数据框然后使用它来避免该问题。

    curve_fit1 <- function(x, y, df) {
  data <- data.frame(x = df[[x]], y = df[[y]])
  co <- coef(lm(log(y) ~ x, data))
  st <- c(C0 = exp(co[[1]]), lambda1 = -co[[2]])
  nls(y ~ C0 * exp(-lambda1 * x), data = data, start = st) 
}

curve_fit1("actual_time", "activity", dat_small)

    给予

    Nonlinear regression model
  model: y ~ C0 * exp(-lambda1 * x)
   data: data
       C0   lambda1 
8.003e+03 1.301e-02 
 residual sum-of-squares: 270020

Number of iterations to convergence: 4 
Achieved convergence tolerance: 4.228e-06

    2)更简单的方法是使用NLS.expoDecaystatforbiology 包中的自启动函数，在这种情况下不需要起始值。输出如 (1) 所示。

    library(statforbiology)

curve_fit2 <- function(x, y, df) {
  data <- data.frame(x = df[[x]], y = df[[y]])
  nls(y ~ NLS.expoDecay(x, C0, lambda1), data)
}

curve_fit2("actual_time", "activity", dat_small)

    3)我们可以交替使用substitute将变量代入公式。这可以适用于 (1) 或 (2)。我们将展示它如何应用于 (2)。

    library(statforbiology)

curve_fit3 <- function(x, y, df) {
  fo <- substitute(y ~ NLS.expoDecay(x, C0, lambda1), 
    list(x = as.name(x), y = as.name(y)))
  nls(fo, df)
}

curve_fit3("actual_time", "activity", dat_small)

    给予

    Nonlinear regression model
  model: activity ~ NLS.expoDecay(actual_time, C0, lambda1)
   data: df
       C0   lambda1 
8.003e+03 1.301e-02 
 residual sum-of-squares: 270020

Number of iterations to convergence: 4 
Achieved convergence tolerance: 4.228e-06

    4)另一个选择是将 drc 包与 statforbiology 中的 DRC.exponDecay 一起使用。在这种情况下，接受的公式的简单形式允许我们使用reformulate来生成它，此外，drc从drmhas生成的类对象coef和plot其他方法——尝试methods(class = "drc")。

    library(drc)
library(statforbiology)

curve_fit4 <- function(x, y, df) {
  drm(reformulate(x, y), data = df, fct = DRC.expoDecay())
}

res <- curve_fit4("actual_time", "activity", dat_small)
plot(res)

    截屏

    • 4

相关问题