kaizadp/field_lab_q10

SIDb Q10, R10 calculations

Closed this issue · 9 comments

kaizadp commented

@bpbond I'm having trouble calculating Q10 and R10 from the SIDb data.
I've taken a crack at the function/equation, need help.

  • I took equations from Meyer et al. 2018 GBC. They may not be the ones we want, especially Q10.
  • How do I assign starting values to the parameters for curve-fitting? I always have trouble with that.
  • Each study (ID) has multiple respiration time points (time). Do we calculate the equation parameters per study, or per time-point in each study?
    • when I use group_by(ID), I get this error message:
      Error in nlsModel(formula, mf, start, wts) : singular gradient matrix at initial parameter estimates
    • when I use group_by(ID, time), my output is my starting a and b values, not calculated values.

P.S. respiration units are currently inconsistent across studies. I will work on that next, make it consistent with SRDB units.

field_lab_q10/code/1-initial_processing.R

Lines 208 to 225 in 309b282

calculate_sidb_q10_r10 = function(sidb_timeseries_clean){
# using equations from Meyer et al. 2018. https://doi.org/10.1002/2017GB005644
fit_q10_parameters = function(sidb_timeseries_clean){
curve.nls = nlsLM(response ~ a*exp(temperature * b),
start=list(a=1,
b=0),
data = sidb_timeseries_clean)
coef1 = coef(curve.nls) %>% as.data.frame()
coef_transpose = coef1 %>% transpose() %>% `colnames<-`(rownames(coef1))
coef_transpose
}
a =
sidb_timeseries_clean %>%
group_by(ID) %>%
do(fit_q10_parameters(.))
}

bpbond commented

That Meyer et al. paper is using a classic Q10 formulation. It's not the best, but it's simple and a fine starting point.

Start with a toy example (note I'm not familiar with the nlsLM() you're using above, so just using the basic nls()):

# Construct sample data
library(tibble)
R10 = 1.1
Q10 = 2.2
x <- tibble(temp = 0:20,
            resp = Q10 ^ ((temp - 10)/10))

# Add some noise
set.seed(1)
x$resp <- x$resp + rnorm(nrow(x), sd = 0.2)
plot(x)

Rplot

m <- nls(resp ~ a * exp(temp * b),
         start = list(a = 0.5, b = 0.1),
         data = x)

q10_est <- exp(10 * coef(m)["b"])
r10_est <- predict(m, newdata = tibble(temp = 10))

Here we recover 2.27 for Q10 (instead of the 'true' value 2.2) and 1.03 for R10 (instead of 1.1).

To your question above, re starting values. Well, a is the respiration rate at 0 C, so using some low value should work; you could also use quantile() to find the lowest 5% (or something) of the data and use that. For b, the 0.1 that I used above corresponds to a Q10 of 2.7, so that seemed like a reasonable guess.

bpbond commented

Each study (ID) has multiple respiration time points (time). Do we calculate the equation parameters per study, or per time-point in each study?

Here does "time-point" mean separate incubation, basically? If so yes I'd think the latter (per time-point in each study).

kaizadp commented

@bpbond

When I use nls(), I get an error. I tried different starting values, but no improvement.

Error in numericDeriv(form[[3L]], names(ind), env) : 
  Missing value or an infinity produced when evaluating the model

Apparently, a way around this is to use minpack.lm::nlsLM(), see here. But when I do that, the output file has my starting values for a and b throughout.

bpbond commented

Need to narrow down where this error is occurring, and why. For all data? Can you give me exact steps to reproduce it?

kaizadp commented

Line 163 onwards. I loaded the SIDb data, combined all the studies and cleaned, then ran the fit_q10_parameters() function. It works when I run it for the entire dataset, but not when I group, or run it on individual studies, eg for the "Andrews" studies.

     sidb_timeseries_clean %>%  
     filter(grepl("Andrews", ID) %>%  
     do(fit_q10_parameters(.))

edit: forgot to tag you, @bpbond !

bpbond commented

First, I'm already in this thread, so you don't need to re-tag me every time.

Second, trying to run your code and am hitting error after error. I'll open a PR with the first fixes I did but maybe we can touch base again via Teams, and I can watch you run this? Because it's not working for me. Thanks!

kaizadp commented

were the errors in the "misc code" section (line 304-356)? If so, my bad. that was miscellaneous code I was using to build the functions, and would not work otherwise. I have commented those lines.

If it's an error in other places, let's meet on Teams.

I do get an error for the final piece of calculate_sidb_q10_r10, where I am unable to fit the equation:

field_lab_q10/code/1-initial_processing.R

Lines 220 to 225 in 309b282

a =
sidb_timeseries_clean %>%
group_by(ID) %>%
do(fit_q10_parameters(.))
} 

Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates
bpbond commented

Sigh. Okay, now the file sources without error, as no code calls calculate_sidb_q10_r10().

How do I trigger the error above?

kaizadp commented

I was running individual pieces within the function all along. Just added a line to call the function.
516e914

  • 👍1