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

Question

我有一个数据集，我想要拟合一个Gompertz模型，该模型由4个不同的因素(主题、种族、目标和分心者)组成。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()。您可以考虑在函数的nls()部分包含一个try()，这样，如果它在前两个组上成功，但在第三个组上成功，它仍然会继续，您应该仍然会得到一些结果。

nlsLM()

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