R：将“group_by”与“nls()”一起使用

Question

我有一个数据集，我想拟合按 4 个不同因素（主题、种族、目标和干扰因素）分组的 Gompertz 模型。 Gompertz 模型在应用于整个数据集时有效（即，不应用“group_by”）。 group_by 函数在我使用（更简单的）线性回归时起作用。但是，当我尝试将 group_by 与 Gompertz 模型一起使用时，我收到以下错误：

Error in chol2inv(object$m$Rmat()) : 
  element (3, 3) is zero, so the inverse cannot be computed
In addition: Warning messages:
1: In nls(yt ~ ymin + ymax * (exp(-exp((alpha * 2.718282/ymax) * (lambda -  :
  Convergence failure: false convergence (8)
2: In nls(yt ~ ymin + ymax * (exp(-exp((alpha * 2.718282/ymax) * (lambda -  :
  Convergence failure: singular convergence (7)

代码如下：

grouped_data = all_merged %>%
  group_by(subject,race,target,distractor)

gomp_fits = do(grouped_data, tidy(nls(yt ~ ymin+ymax*(exp(-exp((alpha* 2.718282/ymax)*(lambda-time)+1))), data = ., start = list(lambda = 0.480, alpha = 5.8, ymin = 0, ymax = 1.6),
                                 control = list(warnOnly = TRUE),
                                 algorithm = "port",
                                 lower = c(0,-Inf, -Inf, 0),
                                 upper= c(2, Inf, Inf, 2)))) 

谢谢！

Answer 1

TLDR

考虑 nlsLM，一个自启动的 Gompertz 模型或使用一种方法来计算起始值，在 group_modify 工作流中使用它。

可能是这样的（尽管可能不需要上限和下限

fit_gomp <- function(data, ...) {
    nlsLM(formula = y ~ SSgompertz(x, Asym, b2, b3),
          data = data,
          lower = c(0,-Inf, -Inf, 0),
          upper = c(2, Inf, Inf, 2),
          ...) %>% tidy()
}

data %>%
  group_by(subject, race, target, distractor) %>%
  group_modify(~ fit_qomp(data = .x), .keep = TRUE)

获取起始值

虽然我没有使用过 Gompertz 模型，但请考虑一下您是否可以找到一种方法以数学方式获取起始值。

例如，假设我想拟合一个二次平台模型（但是它只有 3 个起始参数）。首先，我有一个定义方程的函数，稍后将进入nls。

# y = b0 + b1x + b2x^2
# b0 = intercept
# b1 = slope
# b2 = quadratic term
# jp = join point = critical concentration

quadp <- function(x, b0, b1, jp) {
    b2 <- -0.5 * b1 / jp
    if_else(
        condition = x < jp,
        true  = b0 + (b1 * x) + (b2 * x * x),
        false = b0 + (b1 * jp) + (b2 * jp * jp)
    )
}

第二部分是制作拟合二次多项式的拟合函数，使用这些系数作为 nls 部分的起始值，并拟合 nls 模型。

fit_quadp <- function(data, ...) {
    # get starting values from simple quadratic
    start <- lm(y ~ poly(x, 2, raw = TRUE), data = data)
    start_values <- list(b0 = start$coef[[1]], # intercept
                         b1 = start$coef[[2]], # slope
                         jp = median(data$x)) # join-point
    
    # nls model that uses those starting values
    nlsLM(formula = y ~ quadp(x, b0, b1, jp),
        data = data,
        start = start_values,
        ...
    ) %>% tidy()
}

... 是在需要时为 nls.control 添加参数。

分析分组数据

至于分析分组数据，我使用group_modify()，因为它返回一个数据框，而group_map() 返回一个列表。所以我的基本工作流程如下：

dataset %>%
    group_by(grouping_variable_1, grouping_variable_2, ...) %>%
    group_modify(~ fit_quadp(data = .x), .keep = TRUE)

然后出现一个包含所有整洁统计信息的表格，因为函数中使用了tidy()。您可以考虑将try() 包裹在函数的nls() 部分周围，这样如果它在前两组中成功，但在第三组中，它仍然会继续，您仍然应该得到一些结果。

nlsLM()

另外，如果你想使用minpack.lm 中的nlsLM，那里的算法比nls() 中的算法更成功。有些人担心错误收敛，但我还没有在我的应用程序中看到它。同样使用nlsLM，您可能不需要担心上限和下限，尽管它们仍然可以设置。