3
$\begingroup$

I am trying to model the relationship between forest age and individual tree mortality rate. The probability of mortality declines rapidly as forests go from being very young, and then creeps back up as they age (Juveniles more likely to die than intermediate aged trees. After this initial decline in mortality probability, mortality creeps back up as a function of age).

Here are some example data generated in R, and a plot of the relationship. I build the example data using two independently parameterized exponential functions, and then summing them together.

#age range, combine two curves to get J
age <- seq(1:200)
age.m1 <- 0.1    * exp(-age * .1) #mortality rate declines initially from juvenile stage
age.m2 <- 0.001 * exp(age* 0.015) #mortality rate creeps up as trees get older
age.m3 <- age.m1 + age.m2 #combine the two mortality rates
#generate some scatter in the data to plot
scatter <- age.m3 + rnorm(length(age.m3),0, 0.002)

#plot data and true relationship
plot(scatter, pch=16, ylab = 'tree mortality rate', xlab='forest age')
lines(smooth.spline(age.m3), lwd=3, col='purple')

The outcome looks like this: enter image description here

I'd like to fit a function that could capture this non-linear relationship, so that I could compare the parameters of these fits to different sets of data. Specifically, I'm asking:

  1. Which equation could fit an asymmetric curve like this using any statistical software?
  2. Can you demonstrate it works using the nls package, or similar, in R?
$\endgroup$
  • 1
    $\begingroup$ There isn't enough information here to provide definitive answers. What's needed are some indications of the shape of your data, of the error structure, and any relevant theoretical suggestions. To illustrate what I am requesting, take a look at the considerations involved in the solutions at stats.stackexchange.com/a/64039/919 and stats.stackexchange.com/a/148166/919, both of which solve problems that have exactly the characteristics you describe here (parameterized nonlinear U-shaped curves being fit using R). $\endgroup$ – whuber Jun 2 '16 at 19:48
  • 1
    $\begingroup$ I don't see any cross-post: are you perhaps suggesting one of my links does answer your question? If so, that's fine and we can point this thread at it. I see you did provide a narrative description that leads to a U-shaped curve, but a more useful description would be quantitative (although I realize that can be problematic to come by). We should bear in mind that the model needn't be perfect, but only good enough. However, if we have the opportunity to adopt a model that is likely to match the reality and isn't overly complicated, then we should try to do so. $\endgroup$ – whuber Jun 2 '16 at 20:00
  • 1
    $\begingroup$ It is--but I had understood it as being only a qualitative illustration rather than as a statement of a specific model. Important information like that shouldn't be buried in code! (Not everybody is conversant with every software platform, but it's reasonable to expect readers to understand English and standard mathematical notation.) Furthermore, once I saw that your simulated errors were inappropriate for mortality rates (they should have a Poisson-like distribution, not a Normal one) I stopped paying any more attention to those details, figuring they weren't meant to be studied closely. $\endgroup$ – whuber Jun 2 '16 at 20:06
  • 1
    $\begingroup$ What's wrong with the function you used to generate the data? In what way would it not already be a good answer to your question? $\endgroup$ – Glen_b Jun 3 '16 at 13:56
  • 1
    $\begingroup$ @Glen_b I generated the data by summing two exponential functions. If I attempt to fit two exponential functions in nls (or anything else), with identical format, it will throw a singularity error. $\endgroup$ – colin Jun 3 '16 at 14:32
3
$\begingroup$

You just need to choose some sensible options in nls. Here I used constrained nls, forcing the coefficients to have opposite signs. I ran this after your code above:

mfit=nls(scatter ~ a1*exp(-b1*age)+a2*exp(b2*age),
    start=list(a1=.01,a2=.02,b1=.04,b2=.04),lower=list(0,0,0,0),
    algorithm="port",trace=TRUE)
 lines(age,fitted(mfit),col="yellow2",lwd=2,lty=3)

enter image description here

The fitted curve is the light dashed line right up the middle of the purple curve

$\endgroup$
  • $\begingroup$ ah, I didn't realize you could force the coefficient signs that way, given that nls may explore a coefficient space that flips the sign. Great! $\endgroup$ – colin Jun 6 '16 at 0:52
  • 1
    $\begingroup$ Even without using formal constraints, with a bit of cleverness and reparameterization you may be able to achieve a similar effect (i.e. impose positivity), so similar models might be fitted in packages that don't have formal constraints as options. $\endgroup$ – Glen_b Jun 6 '16 at 1:17
4
$\begingroup$

FIRST APPROACH:

Hurst et al. use the following form for a similar problem (as far as I can see):

mortality = a + b * age * exp(c * age)

where a, b, and c are parameters. The call to nls would be along the lines

dat <- data.frame(age=age, mortality=scatter)
fit <- nls(
    mortality ~ a + b * age * exp(c*age), 
    data=dat, start=list(a=0.05, b=-0.001, c=-0.001)
)
plot(mortality~age, data=dat)
lines(predict(fit))

where dat is the data.frame with your data, including the columns mortality and age. The fit does not look very nice, though.

enter image description here

SECOND APPROACH:

To exploit the known functional relation, you can use constrained optimisation to get a reasonable fit:

ll <- function(par, dat) -sum(dnorm(dat[,"mortality"], 
        mean=par[1]*exp(par[2]*dat[,"age"]) +   
             par[3]*exp(par[4]*dat[,"age"]), 
        sd=par[5], log=TRUE)
)

fit_mle <- constrOptim(
    theta=c(1,-0.05,1, 0.05, 1), 
    f=ll, grad=NULL, 
    ui=diag(c(1,-1,1,1,1)), ci=rep(0,5), dat=dat
)

p <- fit_mle[["par"]]
y <- p[1]*exp(p[2]*age)+p[3]*exp(p[4]*age)
plot(mortality~age, data=dat)
lines(y)

The important constraints are that the second parameter needs to be negative and the fourth parameter positive.

enter image description here

$\endgroup$
  • $\begingroup$ this is progress! But yeah, the fit is not great, considering these are idealized data where we know the true functional form. $\endgroup$ – colin Jun 3 '16 at 14:33
  • $\begingroup$ Edited the answer with an second approach. $\endgroup$ – Karsten W. Jun 3 '16 at 15:26
  • $\begingroup$ you have nailed it! $\endgroup$ – colin Jun 3 '16 at 15:27
  • $\begingroup$ What's the difference between the optimization procedure in nls & constrOptim? (You can use constraints in both, I think.) $\endgroup$ – Scortchi - Reinstate Monica Jun 3 '16 at 18:45
  • $\begingroup$ Interesting! I did not know about constrains in nls $\endgroup$ – Karsten W. Jun 4 '16 at 9:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.