6

I have been working on a curve fitting script which fits 3 exponentially modified Gaussians (EMGs) to a convolved curve. My base function is similar to a Gaussian distribution, but includes a third parameter (the first two being mu and sigma) which determines the weight of the exponential component of the function.

So overall, each EMG peak takes 3 parameters, plus an amplitude coefficient (in order to match experimental data with values > 1.0)

To deconvolve 3 EMG peaks the total # of parameters to be minimized is 3x4 = 12

In some cases the fit works nicely, but in many cases it fails to converge, and returns an error like this

Convergence failure: false convergence (8)

after only 50 or so iterations (well below the limit).

By using the trace option, I can see that the result is very close to matching the data. And by plotting the curve from my initial estimate values, it can also be seen that the starting parameters are within reasonable proximity to the data:

Data = black (noise added), initial = orange, final iteration before error = red

Here is a sample of my code, where I call nls():

t <- 0.05
fit <- nls(y ~ emgmix(a,b,c,d,a1,b1,c1,d1,a2,b2,c2,d2), 
  start = list(
    a=pk1coef[2],
    b=pk1coef[2],
    c=t,
    d=y[pk1c[2]]*40,

    a1=pk2coef[1],
    b1=pk2coef[2],
    c1=t,
    d1=y[pk2c[2]]*40,

    a2=pk3coef[1],
    b2=pk3coef[2],
    c2=t,
    d2=y[pk3c[2]]*40),

    lower=rep(0.001,12),

    control = list(maxiter = 1000),
    trace = TRUE,
    algorithm = "port",
)

so why am I getting an error when the algorithm seems to be succeeding?

0:     562831.45:  341.700  10.6000 0.0500000  27623.1  419.300  10.8000 0.0500000  2132.38  497.000  14.1000 0.0500000  1026.47
1:     405050.97:  341.603  10.5350 0.0508866  27738.3  419.883  10.7618 0.0471600  2065.57  498.294  14.0557 0.0465954  1057.21
2:     115191.71:  341.507  10.5354 0.0556600  27858.3  421.299  10.1276 0.0418391  1986.87  503.484  13.9263 0.0356617  1262.92
3:     38076.077:  342.417  11.2347 0.0632863  27377.3  420.770  14.8188 0.0546385  2213.08  505.655  18.1187 0.0495791  1407.27
4:     36401.368:  343.360  11.7864 0.0723805  26889.9  426.228  23.2991 0.115937  2330.60  507.362  26.3221 0.0784007  1706.85
5:     16394.715:  343.437  11.7838 0.0741048  26883.4  423.172  19.5157 0.154983  2482.43  519.106  27.3302 0.165639  1558.34
6:     12437.878:  343.449  11.7884 0.0743107  26868.4  426.309  21.3703 0.207416  2569.34  517.635  24.8692 0.263019  1512.44
7:     12248.298:  343.429  11.7789 0.0740482  26885.0  426.114  20.9106 0.213771  2551.15  516.084  24.6528 0.200320  1527.81
8:     12235.845:  343.430  11.7791 0.0740580  26884.1  426.175  20.9728 0.214606  2555.89  515.690  24.4315 0.192340  1523.57
9:     12230.776:  343.430  11.7794 0.0740656  26883.7  426.227  20.9872 0.217407  2556.37  515.362  24.3697 0.180294  1523.84
10:     12217.446:  343.432  11.7803 0.0740936  26881.7  426.645  21.0955 0.238821  2558.55  514.148  24.1081 0.139162  1524.57
11:     12185.853:  343.435  11.7813 0.0741224  26879.7  427.203  21.2201 0.274725  2561.21  513.228  23.8124 0.126246  1525.05
12:     12174.404:  343.436  11.7819 0.0741410  26878.4  427.589  21.2985 0.310384  2564.07  512.065  23.4146 0.106315  1524.86
13:     12161.212:  343.437  11.7826 0.0741606  26877.1  427.933  21.3557 0.351018  2565.29  512.085  23.3748 0.109496  1524.09
14:     12155.955:  343.437  11.7826 0.0741621  26876.9  428.243  21.3974 0.394982  2565.96  511.729  23.2536 0.104486  1524.67
15:     12152.489:  343.438  11.7827 0.0741652  26876.7  428.497  21.4262 0.441353  2566.25  511.693  23.2270 0.104343  1524.34
16:     12150.125:  343.438  11.7829 0.0741713  26876.3  428.714  21.4500 0.491154  2566.61  511.637  23.2104 0.103651  1524.53
17:     12148.632:  343.438  11.7829 0.0741714  26876.3  429.008  21.4756 0.569129  2566.55  511.659  23.2185 0.103983  1524.51
18:     12147.015:  343.438  11.7827 0.0741674  26876.5  429.225  21.4869 0.653321  2566.19  511.648  23.2186 0.103855  1524.68
19:     12145.989:  343.438  11.7828 0.0741677  26876.4  429.391  21.4974 0.738613  2566.22  511.659  23.2218 0.103998  1524.65
20:     12145.369:  343.438  11.7829 0.0741710  26876.2  429.533  21.5074 0.830413  2566.43  511.651  23.2199 0.103902  1524.69
21:     12145.021:  343.438  11.7829 0.0741707  26876.2  429.685  21.5152 0.947698  2566.43  511.656  23.2211 0.103965  1524.66
22:     12144.714:  343.438  11.7828 0.0741698  26876.3  429.815  21.5202  1.08360  2566.35  511.653  23.2208 0.103927  1524.70
23:     12144.463:  343.438  11.7828 0.0741698  26876.3  429.913  21.5239  1.22124  2566.36  511.656  23.2217 0.103966  1524.69
24:     12144.317:  343.438  11.7829 0.0741705  26876.2  429.992  21.5273  1.35908  2566.42  511.651  23.2198 0.103907  1524.69
25:     12144.214:  343.438  11.7829 0.0741712  26876.2  430.059  21.5299  1.50140  2566.47  511.654  23.2205 0.103943  1524.67
26:     12144.204:  343.438  11.7829 0.0741712  26876.2  430.059  21.5300  1.51704  2566.50  511.650  23.2189 0.103890  1524.67
27:     12144.204:  343.438  11.7829 0.0741713  26876.2  430.059  21.5303  1.51705  2566.51  511.650  23.2188 0.103891  1524.67
28:     12144.204:  343.438  11.7829 0.0741714  26876.2  430.059  21.5305  1.51706  2566.53  511.651  23.2185 0.103891  1524.65
29:     12144.204:  343.438  11.7829 0.0741714  26876.2  430.059  21.5305  1.51706  2566.53  511.651  23.2185 0.103891  1524.65
30:     12144.204:  343.438  11.7829 0.0741714  26876.2  430.059  21.5305  1.51706  2566.53  511.651  23.2185 0.103891  1524.65
31:     12144.204:  343.438  11.7829 0.0741714  26876.2  430.059  21.5305  1.51706  2566.53  511.651  23.2185 0.103891  1524.65
Ryan
  • 203
  • 2
  • 9
  • From `?nls()` : "Setting warnOnly = TRUE in the control argument (see nls.control) returns a non-converged object (since R version 2.5.0) which might be useful for further convergence analysis, but not for inference." – Bryan Goggin Jun 28 '16 at 01:01
  • Plot the SSE in dependence of the parameters; e.g., contour plots let you examine sensitivity of SSE against variation in two parameters. I expect you'll see that there is no clear minimum. – Roland Jun 28 '16 at 07:13
  • Hmm well unfortunately I have 12 parameters so there is little chance of being able to visualize their relationship to the SSE. But I see your point, because the SSE reaches a minimum, but the algorithm continues to modify the parameters with no resulting change in SSE, and declares a false convergence. Could it be that my step size is too small? please see my updated post for the traces output – Ryan Jun 28 '16 at 18:09
  • Also, since this is a composite function with 12 parameters, each parameter's weight is less than the algorithm may be designed for.. – Ryan Jun 28 '16 at 18:12
  • The SSE does not reach a minimum. As you know, the SSE hyper-plane in the parameter space can have all kinds of valleys, plateaus and other features, where there is no sufficient gradient for the algorithm to reach a (numerically) clear minimum. The more parameters you have, the more likely this is to happen. I'd usually consider 12 parameters to be too many for a non-linear model. But you could try different optimizers. The optimx package is handy for this. – Roland Jun 29 '16 at 07:19
  • Ok thanks, I will try this out. I've been using warnOnly recently with some success, but as you said, there is no guarantee that the result is a true minimum, and i have found that my manually identified starting parameters sometimes produce different results for the same data, presumably because the algorithm is getting caught up in a local valley – Ryan Jun 29 '16 at 08:05
  • Did you have a try with the nlsLM function from the minpack.lm package? That's usually a lot more robust than nls... For parameters that you want to constrain to positive values it's easiest for those to fit the log of those parameters. Hope this helps? – Tom Wenseleers Sep 02 '18 at 21:42
  • Well it has been years since I dealt with this problem, but the warnOnly flag got me far enough for for my needs at the time. Another post recommended using the optim method but I wasn't able to implement it. If I were to redo it I would probably reimplement everything in Python, using a general purpose regression library like scikitlearn, or even PyTorch. – Ryan Sep 04 '18 at 18:37

1 Answers1

4

I was having a similar problem so what I did was "force" new iterations by using the warnOnly=T, which would result in actual estimates, then use those estimates as new starting values in a second nls(). This is what my code ended up looking like:

a_start2 = 40
b_start2 = 200
p_start2 = 16.5 * mean(no.stage[Position == i & Stage == j & Year == k]) + 74.167

subset1 = which(Position == i & Stage == j & Year == k)

m2 = nls(
  Percent ~ ((a) / sqrt(2*b*pi)) * exp(-(((DAFB - p)^2) / (2*b))),
  start = list(a = a_start2, b = b_start2, p = p_start2),
  control = list(maxiter = 50000, minFactor = 1/2000, warnOnly = TRUE),
  algorithm = "port",
  lower = list(a = 0.1, b = 100, p = -100),
  upper = list(a = 200, b = 800, p = 400),
  subset = subset1
)

print(summary(m2))


a_start3 = coef(summary(m2))["a", "Estimate"]
b_start3 = coef(summary(m2))["b", "Estimate"]
p_start3 = coef(summary(m2))["p", "Estimate"]

m3 = nls(
  Percent ~ ((a) / sqrt(2*b*pi)) * exp(-(((DAFB - p)^2) / (2*b))),
  start = list(a = a_start3, b = b_start3, p = p_start3),
  control = list(maxiter = 50000, minFactor = 1/2000, warnOnly = TRUE),
  algorithm = "port",
  lower = list(a = 0.1, b = 100, p = -100),
  upper = list(a = 200, b = 800, p = 400),
  subset = subset1
)
MichaelChirico
  • 33,841
  • 14
  • 113
  • 198
  • This works sometimes but not always. It seems that changing the starting values can suppress the error message even if the final values are the same but sometimes the algorithm can converge with different starting values but not when taking the final values as starting values. I have no idea why we have that... – Jean Paul Aug 06 '20 at 14:25