Starting values for 4 parameter NLS - Chapman Richards function












1















*Note - I have read several of the posts on how to find starting values for NLS - however, I have not found one with an equation of this form (i.e. 4 parameters, exponent raised to a power)



I am struggling tremendously to find suitable starting values for the Chapman Richards equation, which is commonly used in forestry to model tree growth.



y(t) = α * (1 - β * exp(-k * t)^{1/(1-m)})


I typically try to find initial values by plotting a line with set parameters, and then tweaking it to fit the data more closely (Image 1). After this I would use the parameters in the function:



initial.test <- chapmanRichards(seq(0:15),42,0.95,0.28, 0.67)
plot(age,topHeight,type="p",xlab="year since planting",ylab="Dom height (m)", xlim = c(0,20), ylim = c(0, 50))
lines(seq(0:15),initial.test,col="red")


enter image description here



nls(topHeight ~ chapmanRichards(age,a,b,k,m),start=list(a=42,b=0.95,k=0.28,m=0.67))


In this case, the program is able to fit the curve with the starting values provided. The problem, however, is when the data is a bit noisy, and after 2 hours of fiddling with the initial test values, I still can't find good enough starting values (Image 2 shows a few attempts on another dataset.



enter image description here



Can anyone advise on what a good way would be to find suitable starting values? I have thought of creating a matrix that basically runs a sequence for each of the parameters and looping the nls with those starting values, but not sure how the code would look. Any other advice would be greatly appreciated!



PS - would this be something more suited to Excel - solver?










share|improve this question

























  • Isn't the Chapman-Richards model y(t) = α * (1 - β * exp(-k * t))^{1/(1-m)}?

    – Roland
    Nov 20 '18 at 8:11






  • 1





    The way I have it is straight from the {growthmodels} package in R. There are supposedly other versions that one can use, such as the 3-parameter derivative to model growth rate.

    – Benjamin van Heerden
    Nov 20 '18 at 8:58








  • 1





    Well, I think you should look in some actual publications and not some random R package. exp(-k * t)^{1/(1-m)} is the same as exp(-k * t * {1/(1-m)}) which can be re-parameterized to exp(a * t), i.e., you have two parameters that cannot be estimated independently.

    – Roland
    Nov 20 '18 at 10:29








  • 1





    @Roland that is true, I have been working through Burkhart and Tome (2012) and the same 3 parameter form is given as in the blog that. It must be b0(1-e(b1*t))^b2

    – Benjamin van Heerden
    Nov 23 '18 at 7:02






  • 1





    @JamesPhillips Thanks! The blog is really useful. I notice it uses the function in the same format as presented in Burkhart and Tome (2012).

    – Benjamin van Heerden
    Nov 23 '18 at 7:03
















1















*Note - I have read several of the posts on how to find starting values for NLS - however, I have not found one with an equation of this form (i.e. 4 parameters, exponent raised to a power)



I am struggling tremendously to find suitable starting values for the Chapman Richards equation, which is commonly used in forestry to model tree growth.



y(t) = α * (1 - β * exp(-k * t)^{1/(1-m)})


I typically try to find initial values by plotting a line with set parameters, and then tweaking it to fit the data more closely (Image 1). After this I would use the parameters in the function:



initial.test <- chapmanRichards(seq(0:15),42,0.95,0.28, 0.67)
plot(age,topHeight,type="p",xlab="year since planting",ylab="Dom height (m)", xlim = c(0,20), ylim = c(0, 50))
lines(seq(0:15),initial.test,col="red")


enter image description here



nls(topHeight ~ chapmanRichards(age,a,b,k,m),start=list(a=42,b=0.95,k=0.28,m=0.67))


In this case, the program is able to fit the curve with the starting values provided. The problem, however, is when the data is a bit noisy, and after 2 hours of fiddling with the initial test values, I still can't find good enough starting values (Image 2 shows a few attempts on another dataset.



enter image description here



Can anyone advise on what a good way would be to find suitable starting values? I have thought of creating a matrix that basically runs a sequence for each of the parameters and looping the nls with those starting values, but not sure how the code would look. Any other advice would be greatly appreciated!



PS - would this be something more suited to Excel - solver?










share|improve this question

























  • Isn't the Chapman-Richards model y(t) = α * (1 - β * exp(-k * t))^{1/(1-m)}?

    – Roland
    Nov 20 '18 at 8:11






  • 1





    The way I have it is straight from the {growthmodels} package in R. There are supposedly other versions that one can use, such as the 3-parameter derivative to model growth rate.

    – Benjamin van Heerden
    Nov 20 '18 at 8:58








  • 1





    Well, I think you should look in some actual publications and not some random R package. exp(-k * t)^{1/(1-m)} is the same as exp(-k * t * {1/(1-m)}) which can be re-parameterized to exp(a * t), i.e., you have two parameters that cannot be estimated independently.

    – Roland
    Nov 20 '18 at 10:29








  • 1





    @Roland that is true, I have been working through Burkhart and Tome (2012) and the same 3 parameter form is given as in the blog that. It must be b0(1-e(b1*t))^b2

    – Benjamin van Heerden
    Nov 23 '18 at 7:02






  • 1





    @JamesPhillips Thanks! The blog is really useful. I notice it uses the function in the same format as presented in Burkhart and Tome (2012).

    – Benjamin van Heerden
    Nov 23 '18 at 7:03














1












1








1


1






*Note - I have read several of the posts on how to find starting values for NLS - however, I have not found one with an equation of this form (i.e. 4 parameters, exponent raised to a power)



I am struggling tremendously to find suitable starting values for the Chapman Richards equation, which is commonly used in forestry to model tree growth.



y(t) = α * (1 - β * exp(-k * t)^{1/(1-m)})


I typically try to find initial values by plotting a line with set parameters, and then tweaking it to fit the data more closely (Image 1). After this I would use the parameters in the function:



initial.test <- chapmanRichards(seq(0:15),42,0.95,0.28, 0.67)
plot(age,topHeight,type="p",xlab="year since planting",ylab="Dom height (m)", xlim = c(0,20), ylim = c(0, 50))
lines(seq(0:15),initial.test,col="red")


enter image description here



nls(topHeight ~ chapmanRichards(age,a,b,k,m),start=list(a=42,b=0.95,k=0.28,m=0.67))


In this case, the program is able to fit the curve with the starting values provided. The problem, however, is when the data is a bit noisy, and after 2 hours of fiddling with the initial test values, I still can't find good enough starting values (Image 2 shows a few attempts on another dataset.



enter image description here



Can anyone advise on what a good way would be to find suitable starting values? I have thought of creating a matrix that basically runs a sequence for each of the parameters and looping the nls with those starting values, but not sure how the code would look. Any other advice would be greatly appreciated!



PS - would this be something more suited to Excel - solver?










share|improve this question
















*Note - I have read several of the posts on how to find starting values for NLS - however, I have not found one with an equation of this form (i.e. 4 parameters, exponent raised to a power)



I am struggling tremendously to find suitable starting values for the Chapman Richards equation, which is commonly used in forestry to model tree growth.



y(t) = α * (1 - β * exp(-k * t)^{1/(1-m)})


I typically try to find initial values by plotting a line with set parameters, and then tweaking it to fit the data more closely (Image 1). After this I would use the parameters in the function:



initial.test <- chapmanRichards(seq(0:15),42,0.95,0.28, 0.67)
plot(age,topHeight,type="p",xlab="year since planting",ylab="Dom height (m)", xlim = c(0,20), ylim = c(0, 50))
lines(seq(0:15),initial.test,col="red")


enter image description here



nls(topHeight ~ chapmanRichards(age,a,b,k,m),start=list(a=42,b=0.95,k=0.28,m=0.67))


In this case, the program is able to fit the curve with the starting values provided. The problem, however, is when the data is a bit noisy, and after 2 hours of fiddling with the initial test values, I still can't find good enough starting values (Image 2 shows a few attempts on another dataset.



enter image description here



Can anyone advise on what a good way would be to find suitable starting values? I have thought of creating a matrix that basically runs a sequence for each of the parameters and looping the nls with those starting values, but not sure how the code would look. Any other advice would be greatly appreciated!



PS - would this be something more suited to Excel - solver?







r parameters curve-fitting modeling nls






share|improve this question















share|improve this question













share|improve this question




share|improve this question








edited Nov 20 '18 at 11:50









toti08

1,74331523




1,74331523










asked Nov 20 '18 at 7:24









Benjamin van HeerdenBenjamin van Heerden

61




61













  • Isn't the Chapman-Richards model y(t) = α * (1 - β * exp(-k * t))^{1/(1-m)}?

    – Roland
    Nov 20 '18 at 8:11






  • 1





    The way I have it is straight from the {growthmodels} package in R. There are supposedly other versions that one can use, such as the 3-parameter derivative to model growth rate.

    – Benjamin van Heerden
    Nov 20 '18 at 8:58








  • 1





    Well, I think you should look in some actual publications and not some random R package. exp(-k * t)^{1/(1-m)} is the same as exp(-k * t * {1/(1-m)}) which can be re-parameterized to exp(a * t), i.e., you have two parameters that cannot be estimated independently.

    – Roland
    Nov 20 '18 at 10:29








  • 1





    @Roland that is true, I have been working through Burkhart and Tome (2012) and the same 3 parameter form is given as in the blog that. It must be b0(1-e(b1*t))^b2

    – Benjamin van Heerden
    Nov 23 '18 at 7:02






  • 1





    @JamesPhillips Thanks! The blog is really useful. I notice it uses the function in the same format as presented in Burkhart and Tome (2012).

    – Benjamin van Heerden
    Nov 23 '18 at 7:03



















  • Isn't the Chapman-Richards model y(t) = α * (1 - β * exp(-k * t))^{1/(1-m)}?

    – Roland
    Nov 20 '18 at 8:11






  • 1





    The way I have it is straight from the {growthmodels} package in R. There are supposedly other versions that one can use, such as the 3-parameter derivative to model growth rate.

    – Benjamin van Heerden
    Nov 20 '18 at 8:58








  • 1





    Well, I think you should look in some actual publications and not some random R package. exp(-k * t)^{1/(1-m)} is the same as exp(-k * t * {1/(1-m)}) which can be re-parameterized to exp(a * t), i.e., you have two parameters that cannot be estimated independently.

    – Roland
    Nov 20 '18 at 10:29








  • 1





    @Roland that is true, I have been working through Burkhart and Tome (2012) and the same 3 parameter form is given as in the blog that. It must be b0(1-e(b1*t))^b2

    – Benjamin van Heerden
    Nov 23 '18 at 7:02






  • 1





    @JamesPhillips Thanks! The blog is really useful. I notice it uses the function in the same format as presented in Burkhart and Tome (2012).

    – Benjamin van Heerden
    Nov 23 '18 at 7:03

















Isn't the Chapman-Richards model y(t) = α * (1 - β * exp(-k * t))^{1/(1-m)}?

– Roland
Nov 20 '18 at 8:11





Isn't the Chapman-Richards model y(t) = α * (1 - β * exp(-k * t))^{1/(1-m)}?

– Roland
Nov 20 '18 at 8:11




1




1





The way I have it is straight from the {growthmodels} package in R. There are supposedly other versions that one can use, such as the 3-parameter derivative to model growth rate.

– Benjamin van Heerden
Nov 20 '18 at 8:58







The way I have it is straight from the {growthmodels} package in R. There are supposedly other versions that one can use, such as the 3-parameter derivative to model growth rate.

– Benjamin van Heerden
Nov 20 '18 at 8:58






1




1





Well, I think you should look in some actual publications and not some random R package. exp(-k * t)^{1/(1-m)} is the same as exp(-k * t * {1/(1-m)}) which can be re-parameterized to exp(a * t), i.e., you have two parameters that cannot be estimated independently.

– Roland
Nov 20 '18 at 10:29







Well, I think you should look in some actual publications and not some random R package. exp(-k * t)^{1/(1-m)} is the same as exp(-k * t * {1/(1-m)}) which can be re-parameterized to exp(a * t), i.e., you have two parameters that cannot be estimated independently.

– Roland
Nov 20 '18 at 10:29






1




1





@Roland that is true, I have been working through Burkhart and Tome (2012) and the same 3 parameter form is given as in the blog that. It must be b0(1-e(b1*t))^b2

– Benjamin van Heerden
Nov 23 '18 at 7:02





@Roland that is true, I have been working through Burkhart and Tome (2012) and the same 3 parameter form is given as in the blog that. It must be b0(1-e(b1*t))^b2

– Benjamin van Heerden
Nov 23 '18 at 7:02




1




1





@JamesPhillips Thanks! The blog is really useful. I notice it uses the function in the same format as presented in Burkhart and Tome (2012).

– Benjamin van Heerden
Nov 23 '18 at 7:03





@JamesPhillips Thanks! The blog is really useful. I notice it uses the function in the same format as presented in Burkhart and Tome (2012).

– Benjamin van Heerden
Nov 23 '18 at 7:03












1 Answer
1






active

oldest

votes


















2














As @Roland pointed out in the comments the parameters in the equation shown in the question are not identifiable so assuming the equation is as he showed:



y = a * (1 - b * exp(-k * t))^{1/(1-m)}


take the log of both sides:



log(y) ~ log(a) + (1/(1-m)) * log(1 - b * exp(-k*t))


and let log(a) = A, 1/(1-m) = M and b = exp(k*B) giving:



log(y) ~ A + M * log(1 - exp(k*(B-t))


Since B is an offset and k is a scaling we can estimate them as B = mean(t) and k = 1/sd(t). Using algorithm = "plinear" we can avoid starting values for the linear parameters (A and M) provided we specify the right hand side as a matrix such that A times the first column plus M times the second column would give the predicted value. Thus we have:



st <- list(B = mean(t), k = 1/sd(t))
fm0 <- nls(log(y) ~ cbind(1, log(1 - exp(k*(B - t)))), start = st,
algorithm = "plinear")


and then back transform the coefficients so obtained to get the starting values for running the final nls.



Also note that nls2 in the nls2 package can evaluate the model on a grid or at a random set of points to get starting values.






share|improve this answer


























  • This is great! Thanks. I will play around with this a bit and see how it fits our data (lots of species and sites to work through!) Will post here once I have some useful feedback.

    – Benjamin van Heerden
    Nov 23 '18 at 7:10











Your Answer






StackExchange.ifUsing("editor", function () {
StackExchange.using("externalEditor", function () {
StackExchange.using("snippets", function () {
StackExchange.snippets.init();
});
});
}, "code-snippets");

StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "1"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);

StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});

function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: true,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: 10,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});


}
});














draft saved

draft discarded


















StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f53388112%2fstarting-values-for-4-parameter-nls-chapman-richards-function%23new-answer', 'question_page');
}
);

Post as a guest















Required, but never shown

























1 Answer
1






active

oldest

votes








1 Answer
1






active

oldest

votes









active

oldest

votes






active

oldest

votes









2














As @Roland pointed out in the comments the parameters in the equation shown in the question are not identifiable so assuming the equation is as he showed:



y = a * (1 - b * exp(-k * t))^{1/(1-m)}


take the log of both sides:



log(y) ~ log(a) + (1/(1-m)) * log(1 - b * exp(-k*t))


and let log(a) = A, 1/(1-m) = M and b = exp(k*B) giving:



log(y) ~ A + M * log(1 - exp(k*(B-t))


Since B is an offset and k is a scaling we can estimate them as B = mean(t) and k = 1/sd(t). Using algorithm = "plinear" we can avoid starting values for the linear parameters (A and M) provided we specify the right hand side as a matrix such that A times the first column plus M times the second column would give the predicted value. Thus we have:



st <- list(B = mean(t), k = 1/sd(t))
fm0 <- nls(log(y) ~ cbind(1, log(1 - exp(k*(B - t)))), start = st,
algorithm = "plinear")


and then back transform the coefficients so obtained to get the starting values for running the final nls.



Also note that nls2 in the nls2 package can evaluate the model on a grid or at a random set of points to get starting values.






share|improve this answer


























  • This is great! Thanks. I will play around with this a bit and see how it fits our data (lots of species and sites to work through!) Will post here once I have some useful feedback.

    – Benjamin van Heerden
    Nov 23 '18 at 7:10
















2














As @Roland pointed out in the comments the parameters in the equation shown in the question are not identifiable so assuming the equation is as he showed:



y = a * (1 - b * exp(-k * t))^{1/(1-m)}


take the log of both sides:



log(y) ~ log(a) + (1/(1-m)) * log(1 - b * exp(-k*t))


and let log(a) = A, 1/(1-m) = M and b = exp(k*B) giving:



log(y) ~ A + M * log(1 - exp(k*(B-t))


Since B is an offset and k is a scaling we can estimate them as B = mean(t) and k = 1/sd(t). Using algorithm = "plinear" we can avoid starting values for the linear parameters (A and M) provided we specify the right hand side as a matrix such that A times the first column plus M times the second column would give the predicted value. Thus we have:



st <- list(B = mean(t), k = 1/sd(t))
fm0 <- nls(log(y) ~ cbind(1, log(1 - exp(k*(B - t)))), start = st,
algorithm = "plinear")


and then back transform the coefficients so obtained to get the starting values for running the final nls.



Also note that nls2 in the nls2 package can evaluate the model on a grid or at a random set of points to get starting values.






share|improve this answer


























  • This is great! Thanks. I will play around with this a bit and see how it fits our data (lots of species and sites to work through!) Will post here once I have some useful feedback.

    – Benjamin van Heerden
    Nov 23 '18 at 7:10














2












2








2







As @Roland pointed out in the comments the parameters in the equation shown in the question are not identifiable so assuming the equation is as he showed:



y = a * (1 - b * exp(-k * t))^{1/(1-m)}


take the log of both sides:



log(y) ~ log(a) + (1/(1-m)) * log(1 - b * exp(-k*t))


and let log(a) = A, 1/(1-m) = M and b = exp(k*B) giving:



log(y) ~ A + M * log(1 - exp(k*(B-t))


Since B is an offset and k is a scaling we can estimate them as B = mean(t) and k = 1/sd(t). Using algorithm = "plinear" we can avoid starting values for the linear parameters (A and M) provided we specify the right hand side as a matrix such that A times the first column plus M times the second column would give the predicted value. Thus we have:



st <- list(B = mean(t), k = 1/sd(t))
fm0 <- nls(log(y) ~ cbind(1, log(1 - exp(k*(B - t)))), start = st,
algorithm = "plinear")


and then back transform the coefficients so obtained to get the starting values for running the final nls.



Also note that nls2 in the nls2 package can evaluate the model on a grid or at a random set of points to get starting values.






share|improve this answer















As @Roland pointed out in the comments the parameters in the equation shown in the question are not identifiable so assuming the equation is as he showed:



y = a * (1 - b * exp(-k * t))^{1/(1-m)}


take the log of both sides:



log(y) ~ log(a) + (1/(1-m)) * log(1 - b * exp(-k*t))


and let log(a) = A, 1/(1-m) = M and b = exp(k*B) giving:



log(y) ~ A + M * log(1 - exp(k*(B-t))


Since B is an offset and k is a scaling we can estimate them as B = mean(t) and k = 1/sd(t). Using algorithm = "plinear" we can avoid starting values for the linear parameters (A and M) provided we specify the right hand side as a matrix such that A times the first column plus M times the second column would give the predicted value. Thus we have:



st <- list(B = mean(t), k = 1/sd(t))
fm0 <- nls(log(y) ~ cbind(1, log(1 - exp(k*(B - t)))), start = st,
algorithm = "plinear")


and then back transform the coefficients so obtained to get the starting values for running the final nls.



Also note that nls2 in the nls2 package can evaluate the model on a grid or at a random set of points to get starting values.







share|improve this answer














share|improve this answer



share|improve this answer








edited Nov 21 '18 at 22:51

























answered Nov 21 '18 at 22:38









G. GrothendieckG. Grothendieck

146k9129233




146k9129233













  • This is great! Thanks. I will play around with this a bit and see how it fits our data (lots of species and sites to work through!) Will post here once I have some useful feedback.

    – Benjamin van Heerden
    Nov 23 '18 at 7:10



















  • This is great! Thanks. I will play around with this a bit and see how it fits our data (lots of species and sites to work through!) Will post here once I have some useful feedback.

    – Benjamin van Heerden
    Nov 23 '18 at 7:10

















This is great! Thanks. I will play around with this a bit and see how it fits our data (lots of species and sites to work through!) Will post here once I have some useful feedback.

– Benjamin van Heerden
Nov 23 '18 at 7:10





This is great! Thanks. I will play around with this a bit and see how it fits our data (lots of species and sites to work through!) Will post here once I have some useful feedback.

– Benjamin van Heerden
Nov 23 '18 at 7:10


















draft saved

draft discarded




















































Thanks for contributing an answer to Stack Overflow!


  • Please be sure to answer the question. Provide details and share your research!

But avoid



  • Asking for help, clarification, or responding to other answers.

  • Making statements based on opinion; back them up with references or personal experience.


To learn more, see our tips on writing great answers.




draft saved


draft discarded














StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f53388112%2fstarting-values-for-4-parameter-nls-chapman-richards-function%23new-answer', 'question_page');
}
);

Post as a guest















Required, but never shown





















































Required, but never shown














Required, but never shown












Required, but never shown







Required, but never shown

































Required, but never shown














Required, but never shown












Required, but never shown







Required, but never shown







Popular posts from this blog

Can a sorcerer learn a 5th-level spell early by creating spell slots using the Font of Magic feature?

Does disintegrating a polymorphed enemy still kill it after the 2018 errata?

A Topological Invariant for $pi_3(U(n))$