在R中如何将By类nls模型列表转换为数据框并提取系数？

模型定义 #

要拟合的双池指数衰减模型为：

mass_loss ~ m1*exp(-k*cdi*t) + m2*exp(-k2*cdi*t)

对应的自定义R函数：

two_pool <- function(m1,k1,cdi,t,m2,k2){(m1*exp(-k1*cdi*t)+
                                           m2*exp(-k2*cdi*t))
}

数据准备 #

decomp_control <- 
  structure(list(site_number = c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 3L, 3L, 3L, 3L, 5L, 5L, 5L, 5L, 2L, 2L, 2L, 2L, 3L, 3L, 
3L, 3L), site_code = c("CCPp1a", "CCPp1b", "CCPp1c", "CCPp1d", 
"CCPp1e", "CCPp1f", "CCPp2a", "CCPp2b", "CCPp2c", "CCPp2d", "CCPp3a", 
"CCPp3b", "CCPp3c", "CCPp3d", "CCPp5a", "CCPp5b", "CCPp5c", "CCPp5d", 
"CMPp2a", "CMPp2b", "CMPp2c", "CMPp2d", "CMPp3a", "CMPp3b", "CMPp3c", 
"CMPp3d"), days_between = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L), mass_remaining = c(100, 100, 100, 100, 100, 100, 100, 100, 
100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 
100, 100, 100, 100, 100), cdi_mean = c(0.636176396, 0.636176396, 
0.636176396, 0.636176396, 0.636176396, 0.636176396, 0.699725486, 
0.699725486, 0.699725486, 0.699725486, 0.664880103, 0.664880103, 
0.664880103, 0.664880103, 0.679079501, 0.679079501, 0.679079501, 
0.679079501, 0.681828353, 0.681828353, 0.681828353, 0.681828353, 
0.705206121, 0.705206121, 0.705206121, 0.705206121), m2 = c(43.65475, 
43.65475, 43.65475, 43.65475, 43.65475, 43.65475, 43.98935, 43.98935, 
43.98935, 43.98935, 44.08055, 44.08055, 44.08055, 44.08055, 44.3788, 
44.3788, 44.3788, 44.3788, 45.22985, 45.22985, 45.22985, 45.22985, 
45.41225, 45.41225, 45.41225, 45.41225), m1 = c(56.34525, 56.34525, 
56.34525, 56.34525, 56.34525, 56.34525, 56.01065, 56.01065, 56.01065, 
56.01065, 55.91945, 55.91945, 55.91945, 55.91945, 55.6212, 55.6212, 
55.6212, 55.6212, 54.77015, 54.77015, 54.77015, 54.77015, 54.58775, 
54.58775, 54.58775, 54.58775)), row.names = c(NA, -26L), class = "data.frame")

当前问题修正与优化 #

1. 分组拟合代码错误修正 #

你当前的by循环存在核心错误：拟合时使用了全局数据decomp_Control而非分组后的子集x，导致每个分组都用全部数据拟合，而非对应站点的数据。修正后的代码：

models <- by(decomp_control, decomp_control$site_code, function(x) {
  fm <- nls(mass_remaining ~ two_pool(m1,k1,cdi_mean,days_between,
                                      m2,k2),
            start = list(k1 = 0.0121, k2 = 1.696e-12), data = x)
  fm
})

2. 高效系数提取方案 #

直接对nls模型对象执行unlist无法正确提取系数，推荐两种可靠方法：

方法1：使用broom包快速生成整洁数据框（推荐）

broom包的tidy()函数可以直接将模型系数转换为标准数据框，结合dplyr和purrr实现分组拟合+系数提取的一站式处理：

library(dplyr)
library(purrr)
library(broom)

# 按站点分组拟合模型并提取系数
model_coefs <- decomp_control %>%
  group_split(site_code) %>%
  map(~nls(mass_remaining ~ two_pool(m1,k1,cdi_mean,days_between,m2,k2),
           start = list(k1 = 0.0121, k2 = 1.696e-12), data = .x)) %>%
  map2_df(unique(decomp_control$site_code), ~tidy(.x) %>% mutate(site_code = .y))

# 查看结果
model_coefs

方法2：手动提取系数（无依赖方案）

如果不想额外安装包，可以手动提取每个模型的系数并组合：

# 提取每个模型的系数，转换为数据框
coef_list <- lapply(models, function(mod) {
  data.frame(t(coef(mod)))
})

# 合并为统一数据框，并添加站点标识
coef_df <- dplyr::bind_rows(coef_list, .id = "site_code")

重要提示 #

你的数据中days_between列全为0，会导致模型中的指数项exp(-k1*cdi*t)和exp(-k2*cdi*t)均等于1，模型退化为mass_remaining = m1 + m2，而你的数据中m1 + m2恰好等于100（与mass_remaining一致），因此拟合出的k1、k2值不具备实际生物学意义，建议检查数据是否包含完整的时间序列。

内容的提问来源于stack exchange，提问作者Daniel Fishburn