Starting values for 4 parameter NLS - Chapman Richards function
*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")
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.
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
|
show 2 more comments
*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")
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.
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
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 asexp(-k * t * {1/(1-m)})
which can be re-parameterized toexp(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
|
show 2 more comments
*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")
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.
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
*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")
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.
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
r parameters curve-fitting modeling nls
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 asexp(-k * t * {1/(1-m)})
which can be re-parameterized toexp(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
|
show 2 more comments
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 asexp(-k * t * {1/(1-m)})
which can be re-parameterized toexp(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
|
show 2 more comments
1 Answer
1
active
oldest
votes
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.
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
add a comment |
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
});
}
});
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
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
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.
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
add a comment |
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.
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
add a comment |
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.
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.
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
add a comment |
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
add a comment |
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.
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
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
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
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
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 asexp(-k * t * {1/(1-m)})
which can be re-parameterized toexp(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