0

Still a novice with R.

Background

I have the data from Li et al. 2003 paper "Belowground biomass dynamics in the Carbon Budget Model of the Canadian Forest Sector".(https://cdnsciencepub.com/doi/10.1139/x02-165) I am trying to recreate equation 6 in R so that I can produce the same parameter estimates. (eq. 6), however, Pf has a maximum value of 0.426, as values greater will produce a 0.

To impose a maximum for the function we can impose limits on the parameters. The easiest is if I can reparameterize the function in terms of the maximum. In the case that RB>0 and the exponential term is decaying, we have for the function Pf(RB)=a+bexp(cRB) the maximum

m = max(Pf) = a+b

So this could be rewritten as

a = m − b

and fit the equation

Pf(RB) = m − b + b·exp(c·RB) = m(1 − b·exp(c·RB))

I would like to code this using nls in r but I don't know how I would go about doing this with the upper limit. I am trying to produce the model in R so I can get the same parameters generated in table 3 equation number 6. Any help would be greatly appreciated.

enter image description here

enter image description here

dww
  • 30,425
  • 5
  • 68
  • 111
James
  • 33
  • 5
  • Any chance you can post the data or a link to it? – Ben Bolker Feb 25 '22 at 23:29
  • https://docs.google.com/spreadsheets/d/1A6Jkm-zHb6pE7QeLvCmFF3_l70xJ9wiEFsQHiijUVPA/edit?usp=sharing Does this work? – James Feb 25 '22 at 23:37
  • https://docs.google.com/spreadsheets/d/1A6Jkm-zHb6pE7QeLvCmFF3_l70xJ9wiEFsQHiijUVPA/edit?usp=sharing this link might be better if anyone else wants to look at the data. – James Feb 25 '22 at 23:51
  • from your last equation, wouldn't `m − b + bexp(cRB) = m - b(1 − exp(cRB))` (instead of `m(1 − bexp(cRB))`) – user20650 Feb 26 '22 at 02:30
  • Please provide enough code so others can better understand or reproduce the problem. – Community Feb 26 '22 at 13:21

1 Answers1

1

We can use nls as follows.

fit = nls(fine_prop ~ m * (1 - b * exp(c*total_root)),
    data = fr,
    start = list(m=0.1, b = 2, c=-0.05))

coefficients(fit)
#          m           b           c 
# 0.06851410 -5.06265496 -0.05369425 

plot(fr)
x = seq(0, max(fr$total_root), length.out=100)
lines(x, predict(fit, newdata = data.frame(total_root = x)))

enter image description here

Equivalently, we can also put the equation in the same format as in the published article, rather than the form in OP. The we have:

fit = nls(fine_prop ~ a + b * exp(c*total_root),
    data = fr,
    start = list(a=1, b = 1, c=-0.1))

coefficients(fit)
#          a          b          c 
#  0.0685124  0.3468601 -0.0536925

The data:

fr = structure(list(total_root = c(28, 47.5, 104.5, 62, 59.8304, 50.9335, 
40.1412, 43.3121, 131.0057, 24.74, 137.71, 25.004, 88.1, 88.1, 
57.6, 57.6, 54.16, 82.82, 88.6, 134.68, 20.92, 25.004, 141.31, 
143.8, 204.04, 104.88, 172.8, 122.62, 157.7, 18.184, 13.538, 
10.4, 27, 36.7, 44.2, 53.9, 78.6, 47.5, 6.6, 10.1, 14.6, 16.8, 
18.3, 8, 9.5, 6.2, 14.1, 15.8, 21.6, 29.1, 33.2, 41, 45, 46, 
59, 37.6, 9.3, 15, 22.7, 43.2, 37.1, 51, 28.6, 60, 75.2, 12.2, 
43.8, 43.3, 6.8, 6, 7.1, 8.9, 94.3, 100.5, 44.9, 9.98, 9.17, 
7.38, 20.9, 9.6, 11.48, 22.68, 27.63, 27.63, 15.19, 27.78, 26.63, 
12.44, 16.7, 0.5, 3.13), fine_prop = c(0.033796, 0.037486, 0.033818, 
0.033774, 0.10147, 0.032788, 0.068508, 0.045253, 0.005412, 0.373484, 
0.092876, 0.196049, 0.030647, 0.051078, 0.144097, 0.182292, 0.182422, 
0.051195, 0.059594, 0.113306, 0.291587, 0.141337, 0.115562, 0.0758, 
0.053911, 0.075324, 0.075231, 0.08563, 0.061509, 0.125825, 0.130743, 
0.203846, 0.37037, 0.174387, 0.144796, 0.079777, 0.075064, 0.042316, 
0.090909, 0.059406, 0.047945, 0.039881, 0.040984, 0.0875, 0.063158, 
0.091935, 0.061702, 0.063924, 0.057407, 0.052234, 0.052711, 0.047317, 
0.042667, 0.043913, 0.033898, 0.035106, 0.26129, 0.28, 0.129956, 
0.306019, 0.190566, 0.117647, 0.158042, 0.106667, 0.142553, 0.142623, 
0.025571, 0.010624, 0.633824, 0.616667, 0.56338, 0.278652, 0.059597, 
0.013433, 0.088864, 0.305611, 0.314068, 0.298103, 0.398086, 0.491667, 
0.408537, 0.260582, 0.299312, 0.103511, 0.357472, 0.236501, 0.241833, 
0.276527, 0.353293, 0.4, 0.389776)), row.names = c(NA, -91L), class = "data.frame")
dww
  • 30,425
  • 5
  • 68
  • 111
  • Amazing!, Thank you very much! This is very helpful. Just a few questions to clarify. Where do the "start = list(m=0.1, b = 2, c=-0.05))" number come from? Just curious because I noticed the coeffients that are in table 3 equation 6 don't match the output of the model in the answer. – James Feb 26 '22 at 01:17
  • Looking at the table in the photo above is it possible they didn't use a non linear model? Maybe I am wrong, but table 3 equation number 6 has an R squared value which I take as meaning the model was linear? If it is indeed linear I am not sure what the original authors did. – James Feb 26 '22 at 01:28
  • The coefficients are close to the eq. 6. If we put it into same form we have 0.07 + 0.35 exp(-0.054 RB). Not sure why they got slightly different. Could be a few reasons (e.g. perhaps they used a weighted fit, or excluded some outliers, etc.) If the exact method is not reported in the paper you may need to contact the authors if the small difference really matters. The start values were just guesses that are close enough to converge on the solution. – dww Feb 26 '22 at 04:36
  • Wouldn't necessarily read too much into the R2. You are correct that it is not strictly an appropriate measure for non-linear models. Nonetheless a lot of people report R2 for these in any case. – dww Feb 26 '22 at 04:38
  • Thank you again @dww that clears up a lot. Sorry for my ignorance but I see where 0.07 and -0.054 come from but where does 0.35 in the equation come from if the coefficient b=-5.06265496? – James Feb 26 '22 at 05:48
  • algebra. multiply out the brackets in your equation, and the 2nd term is `m * -b`. i.e. 0.0685 * 5.06 – dww Feb 26 '22 at 12:21