[R] Non-Linear Regression Help
J C Nash
profjcnash at gmail.com
Fri May 5 15:59:20 CEST 2017
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")
<simpleError in numericDeriv(form[[3L]], names(ind), env): Missing value or an infinity produced when evaluating the model>
> 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] 10000
$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)
}
<bytecode: 0x3263280>
<environment: 0x3021208>
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
<<lamda: 40 SS= 7139465 at C = 5.6661 a = 1.9085 n = 1.2421 8 / 1
<<lamda: 16 SS= 6321710 at C = 7.4018 a = 3.5556 n = 1.3955 9 / 2
<<lamda: 6.4 SS= 5015080 at C = 9.3512 a = 5.1077 n = 1.5166 10 / 3
<<lamda: 2.56 SS= 3724863 at C = 11.195 a = 6.3333 n = 1.6023 11 / 4
<<lamda: 1.024 SS= 3144435 at C = 12.47 a = 6.9964 n = 1.6516 12 / 5
<<lamda: 0.4096 SS= 3044065 at C = 13.076 a = 7.0845 n = 1.6763 13 / 6
<<lamda: 0.16384 SS= 3016569 at C = 13.402 a = 6.7057 n = 1.6978 14 / 7
<<lamda: 0.065536 SS= 2965611 at C = 13.887 a = 5.7652 n = 1.739 15 / 8
<<lamda: 0.026214 SS= 2874080 at C = 14.836 a = 4.0457 n = 1.8266 16 / 9
<<lamda: 0.010486 SS= 2769844 at C = 16.1 a = 1.9237 n = 1.9871 17 / 10
<<lamda: 0.0041943 SS= 2639613 at C = 15.821 a = 0.1568 n = 2.2672 18 / 11
lamda: 0.041943 SS= 1.7977e+308 at C = 12.901 a = -1.5552 n = 2.7895 19 / 12
lamda: 0.41943 SS= 1.7977e+308 at C = 16.976 a = -0.52795 n = 2.4402 20 / 12
<<lamda: 0.16777 SS= 2550653 at C = 16.551 a = 0.15159 n = 2.3095 21 / 12
<<lamda: 0.067109 SS= 2524756 at C = 16.778 a = -0.066675 n = 2.3521 22 / 13
lamda: 0.67109 SS= 1.7977e+308 at C = 17.157 a = -0.52924 n = 2.4441 23 / 14
<<lamda: 0.26844 SS= 2517716 at C = 16.855 a = -0.12164 n = 2.3641 24 / 14
<<lamda: 0.10737 SS= 2501124 at C = 16.986 a = -0.2586 n = 2.3913 25 / 15
lamda: 1.0737 SS= 1.7977e+308 at C = 17.264 a = -0.55996 n = 2.454 26 / 16
<<lamda: 0.4295 SS= 2496748 at C = 17.03 a = -0.29226 n = 2.3988 27 / 16
<<lamda: 0.1718 SS= 2486194 at C = 17.117 a = -0.37629 n = 2.4163 28 / 17
lamda: 1.718 SS= 1.7977e+308 at C = 17.307 a = -0.56916 n = 2.4578 29 / 18
<<lamda: 0.68719 SS= 2483488 at C = 17.143 a = -0.39692 n = 2.421 30 / 18
<<lamda: 0.27488 SS= 2476879 at C = 17.199 a = -0.44853 n = 2.4322 31 / 19
lamda: 2.7488 SS= 1.7977e+308 at C = 17.323 a = -0.57068 n = 2.459 32 / 20
<<lamda: 1.0995 SS= 2475207 at C = 17.214 a = -0.46124 n = 2.4351 33 / 20
<<lamda: 0.4398 SS= 2471092 at C = 17.249 a = -0.49305 n = 2.4422 34 / 21
lamda: 4.398 SS= 1.7977e+308 at C = 17.329 a = -0.56992 n = 2.4593 35 / 22
<<lamda: 1.7592 SS= 2470058 at C = 17.259 a = -0.5009 n = 2.444 36 / 22
lamda: 17.592 SS= 1.7977e+308 at C = 17.28 a = -0.52057 n = 2.4484 37 / 23
<<lamda: 7.0369 SS= 2469799 at C = 17.261 a = -0.50286 n = 2.4444 38 / 23
<<lamda: 2.8147 SS= 2469154 at C = 17.267 a = -0.50778 n = 2.4455 39 / 24
lamda: 28.147 SS= 1.7977e+308 at C = 17.28 a = -0.52008 n = 2.4483 40 / 25
<<lamda: 11.259 SS= 2468993 at C = 17.268 a = -0.50901 n = 2.4458 41 / 25
lamda: 112.59 SS= 1.7977e+308 at C = 17.272 a = -0.51208 n = 2.4465 42 / 26
<<lamda: 45.036 SS= 2468952 at C = 17.269 a = -0.50931 n = 2.4459 43 / 26
lamda: 450.36 SS= 1.7977e+308 at C = 17.27 a = -0.51008 n = 2.4461 44 / 27
<<lamda: 180.14 SS= 2468942 at C = 17.269 a = -0.50939 n = 2.4459 45 / 27
<<lamda: 72.058 SS= 2468917 at C = 17.269 a = -0.50958 n = 2.446 46 / 28
lamda: 720.58 SS= 1.7977e+308 at C = 17.27 a = -0.51006 n = 2.4461 47 / 29
<<lamda: 288.23 SS= 2468911 at C = 17.269 a = -0.50963 n = 2.446 48 / 29
<<lamda: 115.29 SS= 2468895 at C = 17.269 a = -0.50975 n = 2.446 49 / 30
lamda: 1152.9 SS= 1.7977e+308 at C = 17.27 a = -0.51005 n = 2.4461 50 / 31
<<lamda: 461.17 SS= 2468891 at C = 17.269 a = -0.50978 n = 2.446 51 / 31
> print(p.nlxb)
nlsr object: x
residual sumsquares = 2468891 on 6 observations
after 31 Jacobian and 51 function evaluations
name coeff SE tstat pval gradient JSingval
C 17.2692 1404 0.0123 0.991 -733.5 4112
a -0.509779 60.67 -0.008403 0.9938 35772 188.7
n 2.44601 31.75 0.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 at 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.
>
More information about the R-help
mailing list