R语言计算有机材料碳半衰期SSasymp函数奇异梯度报错求助

首先我不确定该问题更适配StackOverflow还是CrossValidated板块，最终判断其更偏向编码实现而非底层统计原理，因此发布在此处，若归类有误请告知。

我正在计算多种有机材料中的碳半衰期：将材料为期数周培养，定期测定释放的CO₂量，再将累计释放的CO₂体积（单位mL）转换为碳质量（单位mg），以此计算每个采样间隔内各样品剩余碳含量。数据规律为易矿化碳消耗阶段出现初始碳含量下降（不同材料下降幅度存在差异），后续碳损失速率明显收窄。

我尝试使用渐近线设为0的SSasymp函数计算各样品半衰期，相关代码和示例数据如下：

dat<-structure(list(Item = c("litter", "woodlt10", "litter", "woodlt10", "chargt10", "woodlt10", "litter", "chargt10", "chargt10"), `0` = c(161.4599767, 178.78608, 154.3154933, 179.5406033, 177.9216, 185.262, 150.8786667, 195.4312667, 227.50085), `1` = c(161.0021445, 178.3139851, 153.6009328, 179.2539234, 177.8203349, 185.262, 150.7417449, 195.358527, 227.3496655), `2.5` = c(158.8259128, 177.5134301, 152.5134086, 178.6545425, 177.7889754, 184.3638163, 149.216371, 195.358527, 227.3496655), `4.5` = c(156.5077532, 175.921231, 151.4148628, 177.7793692, 177.4767007, 183.2183622, 147.201998, 195.0909267, 227.0262222), `6.5` = c(154.7131141, 174.9474735, 150.4432374, 177.1403608, 177.2406706, 182.4578207, 146.234637, 194.8740861, 226.7688705), `9.5` = c(153.2392748, 174.0175268, 149.3042064, 176.5575212, 176.8846807, 181.7943539, 145.5862023, 194.6301544, 226.4292793), `13` = c(152.0007445, 173.2103072, 148.4350239, 176.0309575, 176.5002673, 181.1742383, 145.0268347, 194.4425645, 226.3546808), `16.5` = c(150.9846197, 172.6132263, 147.6816509, 175.5924338, 176.3494115, 180.7311843, 144.555467, 194.3811803, 226.2060901), `21` = c(150.2721712, 172.192254, 147.2036125, 175.3900685, 176.341071, 180.4498045, 144.2670636, 194.355281, 226.1714114), `25.5` = c(149.6342556, 171.7482415, 146.6502626, 175.1314172, 176.2993861, 180.0476477, 143.9400763, 194.2879702, 226.1714114), `30.5` = c(149.0119875, 171.2716008, 146.1358666, 174.8327655, 176.1848876, 179.7659473, 143.5427987, 194.2297192, 226.0823399), `36.5` = c(148.5402568, 170.8086499, 145.6660173, 174.5592093, 176.0286056, 179.5362906, 143.2190717, 194.1430492, 225.949889), `43` = c(148.0427195, 170.2820678, 145.1835833, 174.1679759, 175.8830912, 179.218831, 142.8504933, 194.0126381, 225.76894), `49.5` = c(147.7386827, 170.0050513, 144.8519388, 173.8786241, 175.7664341, 178.9888063, 142.5957979, 193.9764544, 225.6975125), `56.5` = c(147.4501476, 169.7254062, 144.5900736, 173.6467626, 175.6701446, 178.805284, 142.3922732, 193.9764544, 225.6401166), `64.5` = c(147.0743873, 169.3696494, 144.2525808, 173.2666537, 175.5422531, 178.5399775, 142.1173998, 193.8920486, 225.5513622), `73` = c(146.7558811, 169.0445058, 143.940297, 172.871404, 175.4054422, 178.2874291, 141.7951601, 193.7639492, 225.3946395), `81` = c(146.6028383, 168.9047583, 143.8443744, 172.6848769, 175.4054422, 178.1929838, 141.664276, 193.7639492, 225.3946395), `88.5` = c(146.3920556, 168.7163201, 143.7024872, 172.488525, 175.3520018, 178.0604067, 141.4846825, 193.7430944, 225.3551643), `99.5` = c(146.1854778, 168.5426061, 143.5068639, 172.3002049, 175.2961331, 177.9321711, 141.290387, 193.7237412, 225.2926565)), row.names = c(1L, 2L, 4L, 6L, 7L, 11L, 12L, 16L, 38L), class = "data.frame")

dat.a<-data.frame(t(dat))
dat.a$Days<-rownames(dat.a)
colnames(dat.a)<-dat.a[1,]
colnames(dat.a)<-paste("X",seq(1:ncol(dat.a)),sep="")
dat.a<-dat.a[-1,]
library(dplyr)
dat.a<-mutate_all(dat.a, function(x) as.numeric(as.character(x)))
storage <- list()
for(i in names(dat.a)){
tryCatch({
    storage[[i]] <- log(2)/exp(coefficients(nls(dat.a[,i] ~ SSasymp(X10, 0.0001, R0, lrc), data=dat.a))[3])
}, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})} 
library(dplyr)
halflives<-melt(storage)
samplelist<-data.frame(matrix(NA, nrow = 10, ncol = 1))
samplelist$L1<-colnames(dat.a)
halflives<-merge(samplelist,halflives,by="L1",all=TRUE)
library(readr)
halflives$ord<-parse_number(halflives$L1)
halflives <- halflives[order(halflives$ord),]
colnames(halflives)<-c("L1","junk","halflife","ord")
halflives<-halflives[-10,]
halflives$material<-dat$Item

aggregate(x = halflives$halflife, by = list(halflives$material), FUN = mean)

运行时持续报错：

ERROR : singular gradient matrix at initial parameter estimates

我猜测报错原因是渐近线设为0，或是观测周期内响应值变化幅度不足？请问是否有现有代码的修复方案，或是其他计算碳半衰期的可行方法？

修复方案 #

报错原因 #

报错核心来自两个问题：

1. 你固定了SSasymp的渐近线为0.0001，但你的数据中尤其是chargt10类样品，整个观测周期碳损失不到1%，远没到趋近于0的阶段，模型无法识别合理的衰减速率参数，导致梯度矩阵奇异
2. 基础nls和SSasymp的默认初始值估计逻辑是针对完整衰减曲线设计的，对衰减幅度极小的数据初始值估计失效

方案1：调整拟合逻辑，使用更稳健的拟合算法 #

改用minpack.lm包的nlsLM拟合，对初始值敏感度更低、收敛性更好，同时将渐近线改为待估参数，避免人为固定导致的拟合失效：

# 安装加载必要包
install.packages(c("minpack.lm", "reshape2"))
library(minpack.lm)
library(reshape2)
library(dplyr)
library(readr)

# 原有数据处理逻辑保留
dat<-structure(list(Item = c("litter", "woodlt10", "litter", "woodlt10", "chargt10", "woodlt10", "litter", "chargt10", "chargt10"), `0` = c(161.4599767, 178.78608, 154.3154933, 179.5406033, 177.9216, 185.262, 150.8786667, 195.4312667, 227.50085), `1` = c(161.0021445, 178.3139851, 153.6009328, 179.2539234, 177.8203349, 185.262, 150.7417449, 195.358527, 227.3496655), `2.5` = c(158.8259128, 177.5134301, 152.5134086, 178.6545425, 177.7889754, 184.3638163, 149.216371, 195.358527, 227.3496655), `4.5` = c(156.5077532, 175.921231, 151.4148628, 177.7793692, 177.4767007, 183.2183622, 147.201998, 195.0909267, 227.0262222), `6.5` = c(154.7131141, 174.9474735, 150.4432374, 177.1403608, 177.2406706, 182.4578207, 146.234637, 194.8740861, 226.7688705), `9.5` = c(153.2392748, 174.0175268, 149.304206