Re: [R] Non-Linear Regression Help
I have a fair bit of experience with both nls and rating curves. This is not a nls() problem, this is a model problem. The power law rating curve favored by hydrologists would not apply to your data as it's based on the idea that a log-log plot of discharge vs. stage, or state+a in your case is a straight line, statistical assumptions notwithstanding. A log-log plot of your data, plot(discharge~stage,data=yourdata,pch=19,log='xy') clear is not a straight line. The very large discharge at stage=6.53 vs. stage=6.32 says one of two things: 1) there is an error in the data (perhaps the 2592.05 should be 592.05) or 2) the river channel geometry has changed dramatically, as in overtopping its banks or perhaps there's a smaller central channel set into a larger flood channel, similar to the LA river of the movies. The way forward is 1) recheck your data or 2) recheck your data and use a two-piece model with one piece for stage <= 6.32 and a second piece for stage > 6.32. For this second approach to work, you'll need more data than you have given us here. BTW, nls() should work fine if the model/data combination meet the requirements of 1) the model 'fits' the data, 2) the residuals are NIID(0,sigma^2), the parameters C, a, and n are identifiable from the data (should be if the last point is excluded). As always, you'll need good starting values for the parameters (get them from a log-log plot). You may find, based on the residuals, that linear regression (lm, glm) are most appropriate so that the errors meet the criteria of constant variance. If none of this makes sense, buy and study the book Nonlinear regression analysis: Its applications, D. M. Bates and D. G. Watts, Wiley, New York, 1988. ISBN 0471-816434. The nls() application is the easy part. Good luck David Stevens On 5/4/2017 4:58 PM, Zachary Shadomy wrote: > I am having some errors come up in my first section of code. I have no > issue in plotting the points. Is there an easier method for creating a > non-linear regression using C*(x+a)^n. The .txt file is named > stage_discharge with the two variables being stage and discharge. > The data is a relatively small file listed below: > > stage discharge > 6.53 2592.05 > 6.32 559.5782 > 5.96 484.2151 > 4.99 494.7527 > 3.66 456.0778 > 0.51 291.13 > > > > > >> power.nls<-nls(stage_dischargee$discharge~C*(stage_discharge$stage+a)^n, > data=stage_discharge, start=list(C=4, a=0, n=1)) >> C<-coef(power.nls)["C"] >> a<-coef(power.nls)["a"] >> n<-coef(power.nls)["n"] >> plot(stage_discharge$stage, stage_discharge$discharge, pch=17, cex=1.25, > ylab='Discharge (cfs )', xlab='Stage (ft)', font.lab=2, main='Boone Creek\n > St age-Discharge Curve') >> curve(C*(x+a)^n, add=TRUE, col="red") > [[alternative HTML version deleted]] > > __ > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- David K Stevens, P.E., Ph.D. Professor and Head, Environmental Engineering Civil and Environmental Engineering Utah Water Research Laboratory 8200 Old Main Hill Logan, UT 84322-8200 435 797 3229 - voice 435 797 1363 - fax david.stev...@usu.edu [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Non-Linear Regression Help
If you insist on using nls() for anything that you don't understand
extremely well, you will end up with frustration. nls() uses the same
method K F Gauss used (with good understanding of the details) over
200 years ago. The Gauss-Newton approach inside works very well and
efficiently for problems where the assumptions are met, and terribly
most other times. But nls() does have some nice "extras", and rather
than rewrite all the code, we have a wrapnls() function for the codes
in the much more modern nlsr package. It tries (and mostly succeeds) in
getting analytic derivatives in cases like this. Note that nls(), when
you output the diagnostic, tells you that it is having trouble with
the numeric derivative.
I did the following:
1) made a csv file from the data in the posting (Shadomy17.csv)
2) edited the nls() call and added trace and try()
3) ran nlxb() from nlsr. Note that it uses a lot of iterations -- the
problem is quite close to singular. The singular values have NOTHING
to do with the individual parameters. Their display position is just
convenient. Together they show that the ratio of largest / smallest sv
(a measure of the condition number) is very large -- an ill-conditioned
problem. Now we know this -- there's no guessing and hand-waving.
Best, JN
Here's the (rather verbose) output
> library(nlsr)
> shadomy <- read.csv("./Shadomy17.csv")
> power.nls<-try(nls(discharge~C*(stage+a)^n, data=shadomy, start=list(C=4,
a=0, n=1), trace=TRUE))
7585285 : 4 0 1
> print(power.nls)
[1] "Error in numericDeriv(form[[3L]], names(ind), env) : \n Missing value or an infinity produced when evaluating the
model\n"
attr(,"class")
[1] "try-error"
attr(,"condition")
> tmp <- readline("Try a better approach")
> p.nlxb <- nlxb(discharge~C*(stage+a)^n, data=shadomy, start=list(C=4, a=0,
n=1), trace=TRUE)
formula: discharge ~ C * (stage + a)^n
lower:[1] -Inf -Inf -Inf
upper:[1] Inf Inf Inf
$watch
[1] FALSE
$phi
[1] 1
$lamda
[1] 1e-04
$offset
[1] 100
$laminc
[1] 10
$lamdec
[1] 4
$femax
[1] 1
$jemax
[1] 5000
$rofftest
[1] TRUE
$smallsstest
[1] TRUE
vn:[1] "discharge" "C" "stage" "a" "n"
Finished masks check
datvar:[1] "discharge" "stage"
Data variable discharge :[1] 2592.05 559.58 484.22 494.75 456.08 291.13
Data variable stage :[1] 6.53 6.32 5.96 4.99 3.66 0.51
trjfn:
function (prm)
{
if (is.null(names(prm)))
names(prm) <- names(pvec)
localdata <- list2env(as.list(prm), parent = data)
eval(residexpr, envir = localdata)
}
no weights
lower:[1] -Inf -Inf -Inf
upper:[1] Inf Inf Inf
Start:lamda: 1e-04 SS= 7585285 at C = 4 a = 0 n = 1 1 / 0
lamda: 0.001 SS= Inf at C = -843.56 a = 228.63 n = 123.33 2 / 1
lamda: 0.01 SS= Inf at C = -515.08 a = 148.03 n = 84.714 3 / 1
lamda: 0.1 SS= 9.0129e+100 at C = -50.129 a = 37.016 n = 29.648 4 / 1
lamda: 1 SS= 8.5013e+47 at C = 58.103 a = 30.986 n = 13.954 5 / 1
lamda: 10 SS= 4.2141e+31 at C = 49.209 a = 45.025 n = 8.0734 6 / 1
lamda: 100 SS= 9.4088e+10 at C = 17.369 a = 15.101 n = 2.9571 7 / 1
< print(p.nlxb)
nlsr object: x
residual sumsquares = 2468891 on 6 observations
after 31Jacobian and 51 function evaluations
namecoeff SE tstat pval gradient
JSingval
C17.2692 1404 0.0123 0.991 -733.5
4112
a -0.509779 60.67 -0.008403 0.9938 35772
188.7
n2.44601 31.750.07704 0.9434 -126335
0.6456
> sink()
On 2017-05-04 06:58 PM, Zachary Shadomy wrote:
I am having some errors come up in my first section of code. I have no
issue in plotting the points. Is there an easier method for creating a
non-linear regression using C*(x+a)^n. The .txt file is named
stage_discharge with the two variables being stage and discharge.
The data is a relatively small file listed below:
stage discharge
6.53 2592.05
6.32 559.5782
5.96 484.2151
4.99 494.7527
3.66 456.0778
0.51 291.13
power.nls<-nls(stage_dischargee$discharge~C*(stage_discharge$stage+a)^n,
data=stage_discharge, start=list(C=4, a=0, n=1))
C<-coef(power.nls)["C"]
a<-coef(power.nls)["a"]
n<-coef(power.nls)["n"]
plot(stage_discharge$stage, stage_discharge$discharge, pch=17, cex=1.25,
ylab='Discharge (cfs )', xlab='Stage (ft)', font.lab=2, main='Boone Creek\n
St age-Discharge Curve')
curve(C*(x+a)^n, add=TRUE, col="red")
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the
Re: [R] Non-Linear Regression Help
Hi I am not an expert in nonlinear regression, but your data seems to be rather weird. Last five points has almost linear relationship, the first one is several times higher. If there is no error in your data, I doubt that you can model it by power equation. Cheers Petr > -Original Message- > From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of Zachary > Shadomy > Sent: Friday, May 5, 2017 12:58 AM > To: r-help@r-project.org > Subject: [R] Non-Linear Regression Help > > I am having some errors come up in my first section of code. I have no > issue in plotting the points. Is there an easier method for creating a > non-linear regression using C*(x+a)^n. The .txt file is named > stage_discharge with the two variables being stage and discharge. > The data is a relatively small file listed below: > > stage discharge > 6.53 2592.05 > 6.32 559.5782 > 5.96 484.2151 > 4.99 494.7527 > 3.66 456.0778 > 0.51 291.13 > > > > > > > power.nls<- > nls(stage_dischargee$discharge~C*(stage_discharge$stage+a)^n, > data=stage_discharge, start=list(C=4, a=0, n=1)) > > C<-coef(power.nls)["C"] > > a<-coef(power.nls)["a"] > > n<-coef(power.nls)["n"] > > plot(stage_discharge$stage, stage_discharge$discharge, pch=17, cex=1.25, > ylab='Discharge (cfs )', xlab='Stage (ft)', font.lab=2, main='Boone Creek\n > St age-Discharge Curve') > > curve(C*(x+a)^n, add=TRUE, col="red") > > [[alternative HTML version deleted]] > > __ > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting- > guide.html > and provide commented, minimal, self-contained, reproducible code. Tento e-mail a jakékoliv k němu připojené dokumenty jsou důvěrné a jsou určeny pouze jeho adresátům. Jestliže jste obdržel(a) tento e-mail omylem, informujte laskavě neprodleně jeho odesílatele. Obsah tohoto emailu i s přílohami a jeho kopie vymažte ze svého systému. Nejste-li zamýšleným adresátem tohoto emailu, nejste oprávněni tento email jakkoliv užívat, rozšiřovat, kopírovat či zveřejňovat. Odesílatel e-mailu neodpovídá za eventuální škodu způsobenou modifikacemi či zpožděním přenosu e-mailu. V případě, že je tento e-mail součástí obchodního jednání: - vyhrazuje si odesílatel právo ukončit kdykoliv jednání o uzavření smlouvy, a to z jakéhokoliv důvodu i bez uvedení důvodu. - a obsahuje-li nabídku, je adresát oprávněn nabídku bezodkladně přijmout; Odesílatel tohoto e-mailu (nabídky) vylučuje přijetí nabídky ze strany příjemce s dodatkem či odchylkou. - trvá odesílatel na tom, že příslušná smlouva je uzavřena teprve výslovným dosažením shody na všech jejích náležitostech. - odesílatel tohoto emailu informuje, že není oprávněn uzavírat za společnost žádné smlouvy s výjimkou případů, kdy k tomu byl písemně zmocněn nebo písemně pověřen a takové pověření nebo plná moc byly adresátovi tohoto emailu případně osobě, kterou adresát zastupuje, předloženy nebo jejich existence je adresátovi či osobě jím zastoupené známá. This e-mail and any documents attached to it may be confidential and are intended only for its intended recipients. If you received this e-mail by mistake, please immediately inform its sender. Delete the contents of this e-mail with all attachments and its copies from your system. If you are not the intended recipient of this e-mail, you are not authorized to use, disseminate, copy or disclose this e-mail in any manner. The sender of this e-mail shall not be liable for any possible damage caused by modifications of the e-mail or by delay with transfer of the email. In case that this e-mail forms part of business dealings: - the sender reserves the right to end negotiations about entering into a contract in any time, for any reason, and without stating any reasoning. - if the e-mail contains an offer, the recipient is entitled to immediately accept such offer; The sender of this e-mail (offer) excludes any acceptance of the offer on the part of the recipient containing any amendment or variation. - the sender insists on that the respective contract is concluded only upon an express mutual agreement on all its aspects. - the sender of this e-mail informs that he/she is not authorized to enter into any contracts on behalf of the company except for cases in which he/she is expressly authorized to do so in writing, and such authorization or power of attorney is submitted to the recipient or the person represented by the recipient, or the existence of such authorization is known to the recipient of the person represented by the recipient. __ R-help@r-project.org mailing lis
[R] Non-Linear Regression Help
I am having some errors come up in my first section of code. I have no issue in plotting the points. Is there an easier method for creating a non-linear regression using C*(x+a)^n. The .txt file is named stage_discharge with the two variables being stage and discharge. The data is a relatively small file listed below: stage discharge 6.53 2592.05 6.32 559.5782 5.96 484.2151 4.99 494.7527 3.66 456.0778 0.51 291.13 > power.nls<-nls(stage_dischargee$discharge~C*(stage_discharge$stage+a)^n, data=stage_discharge, start=list(C=4, a=0, n=1)) > C<-coef(power.nls)["C"] > a<-coef(power.nls)["a"] > n<-coef(power.nls)["n"] > plot(stage_discharge$stage, stage_discharge$discharge, pch=17, cex=1.25, ylab='Discharge (cfs )', xlab='Stage (ft)', font.lab=2, main='Boone Creek\n St age-Discharge Curve') > curve(C*(x+a)^n, add=TRUE, col="red") [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
