Search code examples
rnon-linear-regressionnls

Coding non linear regression, Exponential decay

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

Solution

  • 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")