如何解决R:将“group_by”与“nls()”一起使用
我有一个数据集,我想拟合按 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
,您可能不需要担心上限和下限,尽管它们仍然可以设置。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。