search

kaizadp / field_lab_q10

0 stars 0 forks source link

SIDb Q10, R10 calculations #1

Closed kaizadp closed 3 years ago

kaizadp commented 3 years ago

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

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

https://github.com/kaizadp/field_lab_q10/blob/309b282a73a4317687a07304244b19a9e81b5b30/code/1-initial_processing.R#L208-L225

bpbond commented 3 years ago

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 3 years ago

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 3 years ago

@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 3 years ago

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 3 years ago

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 3 years ago

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 3 years ago

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:

https://github.com/kaizadp/field_lab_q10/blob/309b282a73a4317687a07304244b19a9e81b5b30/code/1-initial_processing.R#L220-L225

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

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 3 years ago

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