First off, let me say that I'm not certain whether this belongs more in StackOverfloww or CrossValidated. In the end, I figured it was more a question of coding/execution than underlying statistical concepts, so I elected to put it here. If it belongs elsewhere, let me know.
I am attempting to calculate the half-life of the carbon in several types of organic materials. I incubated these materials for several weeks and measured the CO2 emitted over time. I then converted the cumulative mL CO2 emitted to mg C, which allowed me to calculate the C remaining in each of the samples at each sampling interval. As you'll note, there is an initial drop (whose degree differs between the materials) as the easily mineralized C is consumed, after which C loss is limited.
I then attempted to calculate the half-life for each sample using the SSAsymp function with 0 as the asymptote. I've included the code and some sample data below:
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)
However, I keep getting an error:
ERROR : singular gradient matrix at initial parameter estimates
which I presume is because I'm setting the asymptote to zero? Or because the response value doesn't change enough over time? Can anyone identify a solution, either in terms of my existing code or another way of calculating half-life?
Many thanks for sifting through my crude code and taking the time to help me out!
