1

I'm trying to find an equation to fit my data. I recognize the shape to be y=-exp(x) but nls(y~-a*exp(x*b)) with various parameter start values fail. y is negative so the "easy" fit of log(y)~log(a)+bx doesn't work well. I tried log(y+2)~a+bx to make everything positive but that didn't result in the right shape.

Can someone provide some help for fitting this data? Thanks!

my data

y=swediff

x=avgdate

dat2=structure(list(swediff = c(0.0379635202678687, 0.0845477936160927, 
0.146010217481196, 0.0416237104326292, 0.0659140490644253, 0.134535534695029, 
0.0095147654468483, 0.238456044233877, 0.276025694437364, 0.29435448415394, 
0.00301157777812485, 0.19171002685605, 0.277759059448242, 0.00400780564144798, 
0.342605838471721, 0.236804884903432, 0.151048712082562, 0.188620966368049, 
-0.0615972418208484, -0.00184933102124457, -0.0163171325413688, 
0.00370250929658511, 0.30014673206306, 0.135354035472228, 0.00699671782210069, 
0.0174510674253347, -0.0145499677497698, 0.0113155610814752, 
-0.0683884523999768, 0.20157093417998, 0.186320361855075, -0.115609443650563, 
0.069177592825418, -0.0161221161393796, 0.150181081582068, 0.0632121126749741, 
0.0769960292118834, 0.061783685314432, 0.0442014176783082, -0.00990798027657931, 
-0.00186219548019918, 0.0274216740478325, 0.118878480695049, 
0.0592089915185285, -0.00823096478874009, 0.120750948230554, 
0.278594307094423, -0.0111994006625954, 0.0379360193757585, 6.29460162030332e-05, 
0.0602068958909111, 0.173755367986025, 0.135902420389977, 0.124863098282806, 
0.0706487190649132, 0.284186487140106, 0.395475978216811, 0.207437141623493, 
0.125297332063683, 0.173198938943319, 0.286364692763195, 0.187181622294097, 
0.0155249739392275, 0.0284811383447792, 0.3302152196822, 0.279809724823687, 
0.137578862533244, 0.116541887214791, 0.0092271882470954, -0.244303611220416, 
-0.0553353805122039, 0.000651263069993122, 0.095154425644194, 
0.208713883186232, 0.188079755592783, 0.205982363732731, 0.0411534316934769, 
0.0245409797776046, 0.309614117523271, 0.355492037841145, 0.288275778246439, 
0.321807979017761, 0.287947161484914, 0.156217353669998, 0.172890522219893, 
-0.0179090707980069, 0.0179643473340213, 0.25572307647916, 0.129243837667331, 
-0.122599542543602, -0.250851649402414, -0.590393411432316, -0.358595360363101, 
-0.0550622410807962, 0.117220026030236, 0.164838184961244, -0.069042978776428, 
0.167357636244551, 0.196797358553372, 0.0360344219433584, 0.104666349576576, 
0.180262208126638, 0.310581785237572, -0.0391739945001232, 0.183767269585731, 
-0.152544955460898, -0.158059621319118, 0.0295497815037277, 0.279152618815129, 
0.449891855544769, 0.0647285940606174, 0.198100342057544, 0.180654262068054, 
-0.0276951482994757, -0.265654834430295, -0.391084042220014, 
-0.232183493444027, -0.350719544866729, -0.430609172084572, -0.360201519737938, 
-0.150370469022476, 0.193858381273031, 0.0324799885766501, 0.00919239862196636, 
0.0389874409767672, 0.0735492920963668, -0.0246395678120185, 
-0.456275023145791, -0.454966261045632, -0.449668204190322, 0.109736227977039, 
0.258795304502191, 0.097510631994483, 0.311711079270686, -0.214511654649866, 
-0.369907957164848, -1.10256371343854, -0.540460735663232, -0.524390818371818, 
-0.807243819575856, -0.821087579040873, -0.464518686689061, -0.00423706426563258, 
-0.215993248122225, -0.250111278782723, -0.388687939032414, -0.360549487118099, 
-0.650590724024303, -0.519445819136207, -0.202899217122569, 0.0586223761941996, 
0.216084208573312, 0.0706260521616783, 0.0756847641524316, 0.0770749593558848, 
0.117128900658644, 0.176013651395269, -0.0908276818873188, -0.264146913474664, 
-0.397227204225366, -0.520243226629073, -0.315371107607145, -0.317452140016727, 
-0.288852876616768, -0.341051691256037, -0.221530347691469, -0.459411134488964, 
-0.637940256805166, -0.27643564304224, -0.0101249183839163, 0.286080813914053, 
0.185985595857955, 0.0607903622684369, 0.0739324892640134, 0.0613910749095987, 
0.081924014991581, 0.191652099453867, 0.0671457700438239, -0.239590260546166, 
-0.153731321264683, -0.925001640233871, -0.768208799308804, -0.953874558691501, 
-0.879388655106864, -0.731569006141193, -0.80662535340525, -0.546541080560938, 
-0.0290726402048612, 0.0899201018225928, -0.427683149429864, 
-0.642914036196173, -0.00975702612958795, 0.108679758607406, 
0.115330348818079, 0.0636578050454333, 0.12389689549272, 0.0784572603799192, 
0.0984567282801977, -0.192624142725435, -0.512093030098511, -0.478389465413485, 
-0.621953473157772, -0.14258209361421, 0.215855068127397, 0.10933805943701, 
0.0371862858446161, 0.147743925023497, 0.110106354151748, 0.167891959372529, 
0.0869501603757304, 0.079154830708945, -0.0538071559599136, 0.329630218853886, 
-0.069082249729998, -0.168466429438044, -0.0413130450215624, 
-0.220680450761757, -0.778270577227636, -0.579688137311179, -0.377826042303253, 
-0.172521105687602, 0.0123755037495739, 0.146581091763792, 0.142823263739513, 
0.118485470250309, 0.123284157640724, 0.00262964701743046, -0.324425563039869, 
0.0535619813150857, 0.172683623868712, 0.127449390525177, 0.268301759367799, 
-0.00433449818031201, -0.371211977194441, -0.457630978639089, 
-0.40157582709021, -0.303037159227197, -0.664551261534142, -0.619530450371125, 
-0.310697158309397, -0.146761733349917, 0.0236095942741553, 0.00995413363968556, 
-0.220037507344027, -0.298201547051362, -0.228524953982665, 0.139997538150043, 
-0.0164089841500689, 0.0480316934068684, -0.0202200927067654, 
0.210784261132325, 0.235666847665425, 0.187862324743091, -0.0954145478677223, 
-0.198055375090692, -0.299881653134342, -0.693717473430792, -0.826874682094633, 
-0.537211297019843, 0.0727237791693081, 0.252505366529933, -0.317726197889352, 
-0.279117149104267, -0.0368974550397417, 0.292430819817608, 0.14796346005382, 
0.037310174880905, 0.111369531415814, 0.222644841684763, 0.232955939088023, 
0.0985704990997788, -0.045553014664072, -0.207731060229671, -0.229935960514404, 
-0.686597725938266, -0.247064364067584, -0.371521331373238, -0.430762168759812, 
-0.196945305918869, 0.215658033923907, 0.390228305203605, 0.300921043791892, 
0.043076468096275, -0.00079320334809857, 0.252832538725812, -0.0267246881907509, 
-0.0904732083287169, -0.254702182910177, -0.289199140111472, 
-0.903663808758846, -0.730480757062446, -0.399450175857138, -0.204946388704593, 
0.167241850751596, 0.108062173538109, 0.158279079258542, 0.110014597255262, 
-0.0307183404739954, -0.105426447364233, 0.133090339859816, -0.178600042668314, 
-0.271699090774754, -0.342264904547807, -1.09265736425287, -0.699172072891869, 
-0.0979254975665448, 0.268596423677819, 0.155665651298694, 0.0992990830653245, 
0.277629262139193, 0.176822125972084, -0.103178451026634, -0.461191269578747, 
-0.216966326644849, -0.1084886076781, -0.154140082268591, -0.652188752264078, 
-0.657623380576153, -0.0784736154290889, 0.0742726185810828, 
0.0709763425947463, 0.294940653726301, 0.293652811592831, 0.0640793559532362, 
-0.0367757730565124, -0.415018906491611, -0.538095246521145, 
-0.18261089894835, -0.566401860686253, -0.194554300183061, -0.0327973145856545, 
0.106726178908186, 0.0748878637317865, 0.0234611871797437, 0.218270627297708, 
-0.15143333931606, -0.028720599313675, -0.604914220546465, -0.484120047066584, 
-0.629905362306152, -0.363942052067613, -0.1558791539258, 0.0559939704995262, 
-0.267565269947471, -0.507879399175005, -0.163173363661096, -0.0301926306852119, 
0.0141783247953152, -0.078676708722438), avgdate = c(107.642857142857, 
108.142857142857, 111.928571428571, 101.857142857143, 80.1428571428571, 
78.1428571428571, 120.214285714286, 114.785714285714, 110.428571428571, 
112.5, 111.571428571429, 103.214285714286, 68.1428571428571, 
85, 79.6428571428571, 78.3571428571429, 100, 99.6428571428571, 
78.2142857142857, 116.285714285714, 82.5714285714286, 84.4285714285714, 
87, 65.7142857142857, 54.8571428571429, 77.2142857142857, 113.214285714286, 
99.7142857142857, 98.4285714285714, 116.428571428571, 128.928571428571, 
102.785714285714, 96.5, 121.071428571429, 115.214285714286, 116.5, 
114.785714285714, 53, 45.2142857142857, 39.5, 72.5714285714286, 
84.4285714285714, 112.571428571429, 74.3571428571429, 96.8571428571429, 
118.571428571429, 159.642857142857, 164.214285714286, 20.2857142857143, 
49.1428571428571, 59.7142857142857, 43.5, 60.8571428571429, 117, 
64.6428571428571, 31.2857142857143, 16.5714285714286, 16.5714285714286, 
32.3571428571429, 33.5, 42.0714285714286, 47.4285714285714, 68.7142857142857, 
113.071428571429, 135, 154.214285714286, 159.714285714286, 167.5, 
167.642857142857, 143.642857142857, 97.4285714285714, 40.0714285714286, 
22.2142857142857, 18.5714285714286, 39.6428571428571, 66.2857142857143, 
80.3571428571429, 56.6428571428571, 28.7142857142857, 30.8571428571429, 
19.8571428571429, 18.7857142857143, 85, 138.142857142857, 121.357142857143, 
104.785714285714, 64.4285714285714, 147.785714285714, 161.285714285714, 
171.214285714286, 173.357142857143, 168.214285714286, 170.5, 
172.5, 33.5, 47.2857142857143, 134.214285714286, 101.571428571429, 
66.7857142857143, 62.8571428571429, 75.8571428571429, 57.3571428571429, 
135.785714285714, 131.357142857143, 140.214285714286, 167.214285714286, 
171.857142857143, 157, 120.928571428571, 92.2857142857143, 90.2142857142857, 
127.5, 137.285714285714, 167.785714285714, 170.857142857143, 
169.714285714286, 180.714285714286, 174.857142857143, 191.071428571429, 
162.857142857143, 156.142857142857, 143.642857142857, 112.071428571429, 
155.571428571429, 155.928571428571, 156.714285714286, 168.142857142857, 
173.5, 180.428571428571, 175.642857142857, 152.642857142857, 
144.214285714286, 119, 141.428571428571, 164.714285714286, 175.142857142857, 
166.428571428571, 129.571428571429, 143.071428571429, 162.571428571429, 
NA, 187.357142857143, 167.285714285714, 159.642857142857, 152.785714285714, 
174.142857142857, 180.357142857143, 180.357142857143, 173.642857142857, 
177.642857142857, 170.357142857143, 158.928571428571, 141.357142857143, 
117.571428571429, 136.571428571429, 118.714285714286, 165.785714285714, 
169.5, 171.357142857143, 196.785714285714, 198.428571428571, 
195.428571428571, 193.642857142857, 185, 191.928571428571, NA, 
195.642857142857, 169, 168.642857142857, 162.928571428571, 147, 
136.928571428571, 93.6428571428571, 105.285714285714, 130.428571428571, 
139.714285714286, 145.214285714286, 160.928571428571, 164.714285714286, 
171.785714285714, 174.5, 185, 187.071428571429, 191.5, 190.214285714286, 
NA, 173.357142857143, 171.285714285714, 165.214285714286, 186.642857142857, 
178.714285714286, 169.928571428571, 132.428571428571, 133.357142857143, 
137.5, 141.357142857143, 152.642857142857, 151.428571428571, 
157.142857142857, 172.142857142857, 166.714285714286, 187.785714285714, 
179.857142857143, 164.428571428571, 141.357142857143, 140.357142857143, 
174.357142857143, 142.214285714286, 126.285714285714, 150.571428571429, 
154.571428571429, 128.857142857143, 141.785714285714, 163.5, 
169, 165.714285714286, 175.071428571429, 170.285714285714, 180.071428571429, 
180.5, 178.357142857143, 167.285714285714, 152.357142857143, 
136.714285714286, 154.5, 162.571428571429, 168.714285714286, 
164.428571428571, 164, 129.214285714286, 134.142857142857, 152.428571428571, 
169, 169.5, 174.857142857143, 174.857142857143, 177.285714285714, 
188.071428571429, 191.214285714286, 184.5, 176.214285714286, 
170.142857142857, 162.428571428571, 174, 173.428571428571, 167.785714285714, 
165.785714285714, 137.714285714286, 147, 145.642857142857, 141.571428571429, 
150.5, 158.428571428571, 172.428571428571, 174.285714285714, 
171.285714285714, 185.285714285714, 190.857142857143, 177.214285714286, 
166.571428571429, 167.571428571429, 169.428571428571, 159.357142857143, 
156.357142857143, 142.571428571429, 136.928571428571, 158.714285714286, 
159.071428571429, 146.928571428571, 149.214285714286, 152.285714285714, 
159.928571428571, 162.785714285714, 174, 179.285714285714, 181.571428571429, 
176.428571428571, 176.071428571429, 169.285714285714, 150.428571428571, 
134.785714285714, 144, 152.714285714286, 163.785714285714, 159.428571428571, 
172.714285714286, 164.642857142857, 167.857142857143, 165.642857142857, 
179.571428571429, 179.142857142857, 144.428571428571, 176.357142857143, 
152.285714285714, 132, 129.142857142857, 139.071428571429, 161.142857142857, 
171.785714285714, 161.642857142857, 171.357142857143, 172.571428571429, 
172.428571428571, 180.285714285714, 180.857142857143, 175.142857142857, 
164.5, 135.357142857143, 132, 144, 147.714285714286, 166.928571428571, 
178.571428571429, 176.071428571429, 173.785714285714, 169.857142857143, 
183.785714285714, 193.857142857143, 169.714285714286, 160.928571428571, 
144, 148.785714285714, 150, 160.428571428571, 157.785714285714, 
174.357142857143, 175.571428571429, 171.785714285714, 180.142857142857, 
180.857142857143, 172.928571428571, 168.214285714286, 158.5, 
162.214285714286, 166.928571428571, 169.071428571429, 168.5, 
178.428571428571, 188.285714285714, 178.142857142857, 175.785714285714, 
170.214285714286, 154.571428571429, 170.142857142857, 165.214285714286, 
168.642857142857, 167.214285714286, 154.071428571429, 168.357142857143
)), row.names = c(NA, -349L), class = c("tbl_df", "tbl", "data.frame"
), .Names = c("swediff", "avgdate"))
MichaelChirico
  • 33,841
  • 14
  • 113
  • 198
Dominik
  • 782
  • 7
  • 27
  • The RHS of `y~-a*exp(x*b)` cannot exceed 0... – MichaelChirico Sep 12 '16 at 01:44
  • 1
    You do not yet know how to describe an algorithm. You need statistics or mathematics help. – IRTFM Sep 12 '16 at 02:09
  • yes good point. do you have a suggestion for an alternative eqn? – Dominik Sep 12 '16 at 02:19
  • perhaps someone should migrate this to math or crossvalidated, if they think its more appropriate there? I don't know if I can – Dominik Sep 12 '16 at 02:25
  • Looking at `log(y + 2)` vs `x`, the plot wasn't linear so I could see why the equation above didn't provide a good fit. Perhaps consider fitting a piecewise linear model to your data (maybe this is what you mean by hockey stick decay?). Copy paste this code to see what a piecewise linear model would look like when you create a knot at 150: `ggplot(data = dat2, aes(x = avgdate, y = swediff, group = avgdate < 150)) + geom_point() + geom_smooth(method = 'lm', se = F)` – jav Sep 12 '16 at 03:01
  • @jav I actually took a look at strucchange and it recommended two breaks. A single break is optimal around 140. I might go with that but I'm evaluating what form is most meaningful. – Dominik Sep 12 '16 at 03:30

1 Answers1

5

It's little wonder it doesn't converge. Your model doesn't allow the fit to go above 0.

You need an asymptote parameter in your nls model rather than forcing it to be 0.

If you try something like 

nls(y~c-a*exp(x*b),start=list(a=.1,b=.1,c=0.2))

you get a reasonable fit:

Nonlinear fit to data with asymptote parameter added

A statistician will probably ask you why you need a functional form; typical nonparametric regression / smoothing approaches* should describe the relationship as well or slightly better; you could use cross-validation to choose the df.

* I'd expect that a natural cubic spline with say 3 df plus the constant term would do pretty well (it has an extra parameter so that would be little surprise if it does better than the nls exponential fit)

You can do a broken line fit easily with the package segmented:

 library(segmented)
 linmdl <- lm(y~x)
 segmentedmdl <- segmented(linmdl, seg.Z = ~x, psi=150)
 summary(segmentedmdl)

        ***Regression Model with Segmented Relationship(s)***

Call: 
segmented.lm(obj = linmdl, seg.Z = ~x, psi = 150)

Estimated Break-Point(s):
    Est.  St.Err 
154.173   1.962 

Meaningful coefficients of the linear terms:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.1518234  0.0427016   3.555  0.00043 ***
x           -0.0003317  0.0003784  -0.877  0.38134    
U1.x        -0.0191520  0.0015109 -12.676       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1909 on 342 degrees of freedom
Multiple R-Squared: 0.5835,  Adjusted R-squared: 0.5799 

Convergence attained in 2 iterations with relative change 1.425652e-16

As you see it placed the break at about 154, producing this fit:

fit of segmented regression

Glen_b
  • 7,883
  • 2
  • 37
  • 48
  • Thanks. I tried a c- but rather blindly and didn't give good enough starting values. Now that I understand what the c means that makes sense. – Dominik Sep 12 '16 at 03:33
  • I did try a gam fit and it fit pretty well but had a tail going up on the left. Might still be worth using though. Did you have a particular smoothing method recommendation? – Dominik Sep 12 '16 at 03:34
  • 2
    @Domink I've added a hockey-stick fit (where the break is optimized automatically and should be considered a parameter of the model) as well. – Glen_b Sep 12 '16 at 03:35
  • 1
    It partly depends on your needs but `library(splines);lm(y~ns(x,df=3))` should do fine. Or you might look at smoothing splines/penalized splines.. local linear regression or loess(/lowess) might also do quite well... it's hard to give advice without really understanding what you need to achieve. – Glen_b Sep 12 '16 at 03:38
  • oops looks like we were writing at the same time. Thanks for the additional suggestions! Segmented looks more useful than strucchange too. – Dominik Sep 12 '16 at 03:41
  • I think I need a functional form with parmeters that I can compare and interpret across sites. the exp and segmented forms do that well – Dominik Sep 12 '16 at 03:45