Fitting an asymptope curve using NLS in R

Question

I'm trying to fit a model using nls to some data

data <- structure(list(x = c(0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 
                     0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 
                     0.19, 0.2, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 
                     0.3, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 
                     0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 
                     0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6, 0.61, 0.62, 
                     0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.7, 0.71, 0.72, 0.73, 
                     0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.8, 0.81, 0.82, 0.83, 0.84, 
                     0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 
                     0.96, 0.97, 0.98, 0.99, 1), 
               y = c(0, 0.485920598224734, 0.579084967320261, 
                                                       0.629311870308531, 0.66413220309958, 0.688711933099295, 0.708245568969946, 
                                                       0.724788081171333, 0.738425093472615, 0.749807346519394, 0.760647315694837, 
                                                       0.769951765276708, 0.778616890715529, 0.786471444472986, 0.793875046379542, 
                                                       0.800645032394326, 0.806792819019893, 0.812866398378857, 0.818266404087108, 
                                                       0.823723492308131, 0.828706795673146, 0.833244855438536, 0.837942746239689, 
                                                       0.842309558466764, 0.846539372663185, 0.850717812598111, 0.854388218169364, 
                                                       0.85805291548934, 0.862094357393613, 0.865770471216143, 0.86910408996204, 
                                                       0.872403459200274, 0.875571538659132, 0.878585495333505, 0.881485286982333, 
                                                       0.884584867425864, 0.887478950823415, 0.890293118703085, 0.893267117618518, 
                                                       0.895978536975198, 0.89874703884465, 0.901190170391301, 0.903810257727545, 
                                                       0.906407512058681, 0.908776436338728, 0.911179610126438, 0.913645574678197, 
                                                       0.916151496988897, 0.918457630504895, 0.920558266974912, 0.922852983988355, 
                                                       0.924925079201986, 0.926900134143905, 0.928966521106259, 0.930821702771356, 
                                                       0.932751091703057, 0.934743271398807, 0.936695493335617, 0.938704797785199, 
                                                       0.940520021691355, 0.942449410623056, 0.944150469503667, 0.945914319148329, 
                                                       0.947643919285327, 0.949436310186374, 0.951114536061877, 0.952821303193767, 
                                                       0.954670776607586, 0.956343294231812, 0.957907355081771, 0.959505665439393, 
                                                       0.961024060279134, 0.962559579872706, 0.964140765476496, 0.965573536547079, 
                                                       0.966994891115107, 0.968530410708679, 0.970071638553529, 0.971550075634329, 
                                                       0.972954305448526, 0.974472700288267, 0.975876930102463, 0.977189827896224, 
                                                       0.978691097982133, 0.980049661786112, 0.981516682364358, 0.982886662670891, 
                                                       0.984171019208266, 0.985409709735423, 0.986694066272797, 0.988006964066558, 
                                                       0.989222821588606, 0.990358763592773, 0.99155749636099, 0.992716271370266, 
                                                       0.993983503153809, 0.995324942203956, 0.996375260438965, 0.997591117961013, 
                                                       0.998909724006051, 1)), row.names = c(NA, -101L), class = c("tbl_df", 
                                                                                                                   "tbl", "data.frame"))

The data looks like this

I'm trying to fit an asymptotic model, due to the nature of the data, we know this should converge at a point, and the value it converges at is useful information.

I tried fitting a model using ax^b+c and it doesn't seems to work.

> model =  nls(y~a*x^(b),
+              data = data,
+              start = list(a=-10, b=-0.2))
Error in numericDeriv(form[[3L]], names(ind), env, central = nDcentral) : 
  Missing value or an infinity produced when evaluating the model

From what I've gathered the suggests the models don't fit the data. Is this true? I'm generating the data from another source during different time periods, and other time periods produce data which I can model just fine and the data follows a similar shape for example,

data <- structure(list(x = c(0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 
                     0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 
                     0.19, 0.2, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 
                     0.3, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 
                     0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 
                     0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6, 0.61, 0.62, 
                     0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.7, 0.71, 0.72, 0.73, 
                     0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.8, 0.81, 0.82, 0.83, 0.84, 
                     0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 
                     0.96, 0.97, 0.98, 0.99, 1), y = c(0.387696586431811, 0.50631394232361, 
                                                       0.573625136575837, 0.620130450617488, 0.653935039565606, 0.680528424328709, 
                                                       0.702201768036288, 0.720147005264378, 0.73551633943648, 0.748852762970566, 
                                                       0.760487368804423, 0.771175048836208, 0.780975399794722, 0.78935205112075, 
                                                       0.796907591961063, 0.804330695626262, 0.811118100850909, 0.81746184153892, 
                                                       0.823328808396517, 0.828639539118631, 0.833526470880376, 0.838598814687283, 
                                                       0.843386418567692, 0.848240241035659, 0.852398768334271, 0.856808926265603, 
                                                       0.860755554084031, 0.864616097738635, 0.868152170314207, 0.871648511737245, 
                                                       0.875105122007748, 0.878171042611661, 0.881203853921796, 0.884395589842069, 
                                                       0.887593947621097, 0.890540674767407, 0.893235771280999, 0.896036817534682, 
                                                       0.898765023342052, 0.901360791974307, 0.90402277919412, 0.906492732510015, 
                                                       0.909022282554713, 0.911326689401715, 0.913558255802404, 0.915922259378207, 
                                                       0.917922060722445, 0.92012713968811, 0.922193159619905, 0.924504188325663, 
                                                       0.926484124093633, 0.92848392543787, 0.930477104923352, 0.932159057047313, 
                                                       0.934238320696619, 0.936125550442009, 0.938045889481177, 0.939840413203986, 
                                                       0.941654802503063, 0.943429460649604, 0.945190875078635, 0.947071482965268, 
                                                       0.948746813230474, 0.950534715094527, 0.952263020229778, 0.95381915703738, 
                                                       0.955467999867563, 0.957037380392676, 0.958514054895209, 0.960143032149124, 
                                                       0.961705790815482, 0.963202330894282, 0.964698870973082, 0.966195411051882, 
                                                       0.967658841836904, 0.969162003774459, 0.970506241101877, 0.972049134191968, 
                                                       0.973373505943118, 0.974658146541734, 0.976075224315465, 0.977406217925372, 
                                                       0.978863026851637, 0.980081448862696, 0.981425686190114, 0.982551402178591, 
                                                       0.983955236234811, 0.985067708505778, 0.986425189550707, 0.987656855279277, 
                                                       0.98886865543158, 0.990080455583882, 0.991530642651392, 0.992802039532497, 
                                                       0.994033705261067, 0.995066715226964, 0.996291759096778, 0.997523424825349, 
                                                       0.998702115683872, 1)), row.names = c(NA, -100L), class = c("tbl_df", 
                                                                                                                   "tbl", "data.frame"))

model2 =  nls(y~a*x^(b)+c,
+              data = data,
+              start = list(a=-10, b=-0.2, c=1.5))

This works fine.

I need a model that can fit both of these. Am I using the wrong formula or is there a better alternative to NLS? Thanks

Answer 1

The model seems to work fine for the first data set if you use better starting values. For this model, you can get decent starting values by doing a log-log regression:

m1 <- lm(log(y) ~ log(x), data1, subset = x > 0)
model =  nls(y~a*x^(b),
             data = data1,
             start = list(a=exp(coef(m1)[1]), b=coef(m1)[2]))

"Better starting values" is usually the answer to this kind of problem. It can also help to supply nls with the function gradients, which you can do as follows (in this case we have to drop the x=0 data point again because the derivative is NaN there ...)

power_law <- deriv(~a*x^b, c("a", "b"),
                   function.arg = c("a", "b", "x"))

model =  nls(y~power_law(a, b, x),
             data = subset(data1, x>0),
             start = list(a=exp(coef(m1)[1]), b=coef(m1)[2]))

As I'm sure I have commented elsewhere, you can also fit this model via

model <- glm(y~log(x), family = gaussian(link="log"), data = data1)

(although the intercept term will be log(a) rather than a)