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

编程问答 2022-04-24

问题描述

我有一个数据集，我想拟合按 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,0),upper= c(2,Inf,2))))

谢谢！

解决方法

TLDR

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

也许是这样的（虽然上限和下限可能不是必需的

fit_gomp <- function(data,...) {
    nlsLM(formula = y ~ SSgompertz(x,Asym,b2,b3),data = data,lower = c(0,-Inf,0),upper = c(2,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,jp),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()

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

broomdplyrdplyrnlsrrtidyr