r - 在 R 中计算碳随时间的半衰期

Question

首先，我要说的是，我不确定这是否更属于 StackOverfloww 还是 CrossValidated。最后，我认为这更多的是编码/执行的问题，而不是基本的统计概念，所以我选择把它放在这里。如果它属于其他地方，请告诉我。

我正在尝试计算几种有机 Material 中碳的半衰期。我将这些 Material 孵化了几周，并测量了一段时间内排放的二氧化碳。然后，我将累积排放的二氧化碳毫升数转换为毫克 C，这使我能够计算每个采样间隔每个 sample 中剩余的 C。您会注意到，随着容易矿化的碳被消耗，碳会出现初始下降(不同 Material 的下降程度不同)，之后碳的损失就会受到限制。

然后，我尝试使用 SSAsymp 函数(以 0 作为渐近线)计算每个样本的半衰期。我在下面包含了代码和一些示例数据:

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

我认为这是因为我将渐近线设置为零？或者因为响应值随着时间的推移变化不够大？任何人都可以根据我现有的代码或计算半衰期的其他方式确定解决方案吗？

非常感谢您仔细检查我的粗糙代码并花时间帮助我!

Answer 1

以天为单位的半衰期可以计算为 log2(y) 与天数的斜率倒数的负数。

dat2 <- t(dat[-1])
days <- as.numeric(rownames(dat2))
colnames(dat2) <- paste0("X", 1:ncol(dat2))
fm <- lm(log2(dat2) ~ days)
-1/coef(fm)[2, ] # half lives
##        X1        X2        X3        X4        X5        X6        X7        X8        X9 
##  816.2918 1306.5804 1063.3762 1854.8503 4882.6257 1914.4168 1275.1235 8489.2436 7799.1691 

为了存在恒定的半衰期，log2(y) 与天数的关系图应该是一条直线；然而，如下图所示，在问题中显示的数据范围内，情况似乎并非如此，因为早期和后期似乎有不同的半衰期。您可能希望将这些天分为两个或更多段，并分别计算每个段的半衰期。我们将天数范围分为两部分。我们将首先使用 10 天作为分界线，但我们将使用其结果来运行 nls 来拟合两条线——一条在截止日期之前，另一条在不假设任何特定截止日期之后。两条线的方程为 log2(y) = b1 + m1 * days 和 log2(y) = b2 + m2 * days。因为线条上方的区域是凸的，所以我们可以将这两者中的最大值作为在任何点使用的值。通过对 T0 求解方程 b1 + m1 * T0 = b2 + m2 * T0 即可获得截止点。

# starting values
fm1 <- lm(log2(dat2) ~ days, subset = days < 10)
co1 <- coef(fm1)
fm2 <- lm(log2(dat2) ~ days, subset = days >= 10)
co2 <- coef(fm2)

# calculate list of nls objects
fmList <- lapply(1:ncol(dat2), function(i) {
  st <- list(b1 = co1[1,i], m1 = co1[2,i], b2 = co2[1,i], m2 = co2[2,i])
  nls(log2(dat2)[, i] ~ pmax(b1 + m1 * days, b2 + m2 * days), start = st)
})

# plot
fits <- sapply(fmList, fitted)
matplot(days, log2(dat2), col = 1, cex = 0.5)
matplot(days, fits, type = "l", col = 1, lty = 1, add = TRUE)

# calculate statistics (half lives, etc.)
stats <- sapply(fmList, function(fm) {
  co <- as.list(coef(fm))
  with(co, c(co, half1 = -1/m1, half2 = -1/m2, T0 = (b1 - b2) / (m2 - m1)))
})
stats

给予:

      [,1]          [,2]          [,3]          [,4]          [,5]         
b1    7.334025      7.48148       7.267299      7.487643      7.475502     
m1    -0.008487284  -0.004336574  -0.005028027  -0.002671511  -0.0009155759
b2    7.242063      7.435454      7.210472      7.461353      7.46409      
m2    -0.0005902792 -0.0004453948 -0.0005338854 -0.0003596622 -0.0001180285
half1 117.8233      230.5968      198.8852      374.3201      1092.209     
half2 1694.113      2245.199      1873.061      2780.387      8472.529     
T0    11.6452       11.8284       12.64488      11.37197      14.30921     
      [,6]          [,7]          [,8]          [,9]         
b1    7.534028      7.236452      7.610853      7.829888     
m1    -0.003162845  -0.006031012  -0.0006520529 -0.0007230756
b2    7.500386      7.180011      7.60358       7.82288      
m2    -0.0002895582 -0.0004227974 -6.616548e-05 -8.131695e-05
half1 316.1711      165.8096      1533.618      1382.981     
half2 3453.538      2365.199      15113.62      12297.56     
T0    11.70858      10.06396      12.41357      10.92091  

检查

只是为了仔细检查假设我们有

x <- c(0, 3, 6)
y <- c(8, 4, 2)

其半衰期显然为 3。那么

-1/coef(lm(log2(y) ~ x))[[2]]
## [1] 3

已添加

关于评论，假设双指数表示您的意思是由 pracma::mexpfit 拟合的模型

library(pracma)

# fmList has one element per dat2 column with fit parameters, etc.
fmList <- lapply(1:ncol(dat2), function(i) {
  mexpfit(days, dat2[, i], -seq(2))
})

# plots
fits <- sapply(fmList, with,
  a0 + a[1] * exp(b[1] * days) + a[2] * exp(b[2] * days)
)
matplot(days, dat2, col = 1, cex = 0.5)
matplot(days, fits, type = "l", col = 1, lty = 1, add = TRUE)