user2359280
user2359280

Reputation: 1

Problems calculating nls() in R

I'm trying to fit point data with polynomial functions within the program R and want to see which order makes the best fit. I use the nonlinear regression model nls(function, data, start), the function being y~a*(x+c)^b.

What I want to do is to calculate the cross-sectional area of a valley fill deposit. For this reason I need to model how the underlying valley base might look like. I already have the linear scenario and now want to try some polynomial scenarios.

My dataset represents the valley profile and looks like this (with x_combsl being the distance and y_combsl the elevation):

    x_combsl  y_combsl    
1   0.0000000 120.9095
2   0.9904867 120.5066
3   1.9809735 120.0947
4   2.9714602 119.6811
5   3.9619470 119.2492
6   4.9524337 118.8483
7   5.9429204 118.4866
8   6.9334072 118.1120
9   7.9238939 117.7750
10  8.9143806 117.2833
11  9.9048674 116.7698
12 10.8953541 116.2841
13 11.8858409 115.8285
14 12.8763276 115.2949
15 13.8668143 114.6750
16 14.8573011 114.1301
17 15.8477878 113.6537
18 16.8382746 113.2016
19 17.8287613 112.8163
20 18.8192480 112.4945
21 19.8097348 112.0304
22 20.8002215 111.2370
23 21.7907082 110.7463
24 22.7811950 110.2954
25 23.7716817 109.6715
26 24.7621685 109.1829
27 44.5719032 109.0435
28 45.5623900 109.4721
29 46.5528767 110.0491
30 47.5433634 110.4394
31 48.5338502 110.6832
32 49.5243369 111.0763
33 50.5148237 111.8376
34 51.5053104 112.6162
35 52.4957971 113.2467
36 53.4862839 113.8065
37 54.4767706 114.4694
38 55.4672573 114.9547
39 56.4577441 115.4724
40 57.4482308 116.0013
41 58.4387176 116.4606
42 59.4292043 117.0797
43 60.4196910 117.7074
44 61.4101778 118.2127    
45 62.4006645 118.7544
46 63.3911512 119.3134
47 64.3816380 119.9159
48 65.3721247 120.5462
49 66.3626115 121.0418
50 67.3530982 121.5350
51 68.3435849 122.0184
52 69.3340717 122.5490
53 70.3245584 123.1162
54 71.3150452 123.6437

When I try to generate the model I get the following error message:

data<-data.frame(x_combsl, y_combsl)

fit_nls<-nls(y_combsl~a*(x_combsl + c)^b, data=data, start=list(a=1, b=2, c=35))
Error in numericDeriv(form[[3]], names(ind), env) :
   Missing value or an Infinity produced when evaluating the model

Any ideas what might cause the problem? The peak of the polynom is very likely to be around c=35, so the start parameters don't seem to cause problems. Might it be the very linear alignment of the points? I tried it with fewer data points [24:31] but I got the same error message.

I'm not very experienced in programming with R, so specific answers would be great.

Upvotes: 0

Views: 4326

Answers (1)

G. Grothendieck
G. Grothendieck

Reputation: 269346

A few points:

  • the curve is approximately symmetric about roughly 35 so the starting c should equal -35, not +35.
  • given the symmetry b should be even, 2, 4, 6, ... .
  • the model needs the flexibility to shift up and down so add a parameter for that
  • with the last point we have two linear parameters so use algorithm="plinear" to avoid having to come up with starting values for them. Also we are fixing b so we now only have c to give a starting value for and we already determined that we would use -35 for that.

Thus try the following and perhaps additional runs with b=2 and b=6 as well:

nls(y_combsl ~ cbind(1, (x_combsl + c)^4), data, alg = "plinear", start = list(c = -35))

Of courese we are really now down to a polynomial so we might change the model to a full polynomial and use lm:

lm(y_combsl ~ poly(x_combsl, 4), data)

Upvotes: 1

Related Questions