# Fitting a power law to the relationship between two variables

I have a small set of data that I need to fit a curve to (see image and data below). algaeConc    adultDens
2.934961    6665.145469
2.921015    6665.144453
2.910769    6665.142377
2.908331    6665.138925
2.907658    6665.134515
2.907278    6665.128736
2.907130    6665.121034
2.907160    6665.117992
2.907326    6665.124547
2.907594    6665.131266


The data came out of a simulation model for which I was stressing a parameter (filtration rate) then recording what happened to my variables (algae concentration and adult density). Now, I need to fit a curve to this, and I have to admit I never saw this kind of decline before. So, I would like some suggestions about how to fit a curve to this. It seems to me that a power law might be needed, but I have never worked with those before and am not sure how to go about it.

I've spent all day reading about power laws and trying to fit it on R but only managed to confuse myself, so any help will be highly appreciated.

Edit: The X axis is inverted to match the stress on the filtration parameter (increasing stress, caused a decrease in algae concentration).

• Have you tried drawing what power law relationships look like? – Glen_b Apr 19 '17 at 23:18
• I have. It seems to me that the formula below should be the correct one, but I'm not sure: adult density = a * algae concentration^(-b) Right now I'm trying to fit this formula in R using the poweRlaw package while I wait to see some cyber opinions, but no success. It think the last value in the table is bugging up the whole thing, because it rises instead of dropping further. – Renata Dis Apr 19 '17 at 23:20
• I've only just noticed that you have reversed your x-axis in the plot. Why would you do that without saying you had done so? Indeed why do it at all? Just showing the plot without being very clear you had flipped the axis is very likely to mislead people. If you want to lead answerers into giving bad advice, that's a good way to do it. Please unflip the plot or give a compelling reason in your text for flipping it around. – Glen_b Apr 19 '17 at 23:51
• It could also be quite useful to people who might answer if you'd paste the data you plotted (in ASCII form) into your post. After you paste it in, you'll need to select it with the code tool at the top of the edit window to indent it 4 spaces so it maintains its line-by-line format – Glen_b Apr 19 '17 at 23:56
• Thank you Glen. I have implemented the changes you asked for. Also, I'm sorry about the X axis, I plotted this way right away because of the stress test, and completly forgot to point it out. – Renata Dis Apr 20 '17 at 0:08

Here's the unflipped version of your plot: No choices of $a$ and $b$ in $y=ax^{-b}$ (where $y$ is adult density and $x$ is algae concentration) will fit that shape of relationship; they don't increase and flatten off up to a positive asymptote: (with negative $a$ you can get it to increase to an asymptote, but the fit will be all negative)

Indeed, the range of your $x$-variable is so narrow (relative to its mean say) that a typical power curve will look effectively like a straight line.

A quick power-law diagnostic plot:

If you don't have intuition for power laws, to see whether a power law is likely to be a reasonable description, plot log y against log x (it doesn't matter which base of logs). It should look close to linear on that scale.

In your case we can see immediately that won't happen (indeed the coefficients of variation are so small that the log-scale plot looks a lot like the first plot above aside from the axis scales).

There are a number of nonlinear functions used for this kind of "saturation" relationship we see in the plot of the data (dozens of them, really). It's best to proceed from theoretical notions of what is happening (which might lead to a differential equation, say, and then to a functional form) than to guess at functional forms; indeed that sort of theoretical understanding is where most of the common increasing-to-an-asymptote type models arise.

Here's a few examples of nonlinear models that increase to an asymptote (at least for some parameter values); most of them will fit your data poorly:

\begin{eqnarray} E(y) &=& ax/(b+x) \\ E(y) &=& a(1-\exp(-bx)) \\ E(y) &=& a (1 - \exp(-b(x - c))) \\ E(y) &=& a+(c-a)\exp(-bx) \\ E(y) &=& a - c(x-d)^{-b} \\ E(y) &=& a \exp(-b/x) \end{eqnarray}

There are many others besides. You should be looking at subject-area knowledge to develop a model, though -- trawling through such a list of models to find one that fits effectively ruins any statistical inference on the same data, so unless you want to go collect another set of data (or indeed possibly several more), don't do that. Strictly speaking you should have a model in place before you even collect data.

Indeed there further consequences to using data-generated models / data-generated hypotheses that I haven't touched on there; there are a number of posts on site that discuss the effect of this kind of approach to choosing models.

In the situation where there's simply no subject-area knowledge to work form, and where the experimenter can collect any amount of data easily, such curve fitting strategies (trawling through a list) may have some value in helping to come up with plausible models -- as long as we take great care to make sure that the choice isn't just an artifact of our one set of data and we take similarly great care not to let the model-selection process affect our ultimate estimation and inference.

For example if we take enough data that we can do out of sample testing (i.e. we can see if the model performs well on data we didn't use to pick the model) and other people repeat the experiment (so we can see whether the same curve tends to be picked up across multiple independently-run experiments - that it's robust to the small changes in the setup that naturally occur when one goes to a different lab with their own pieces of equipment and different people running the experiment).

Once that's all done and the model actually can be seen to "work" in some sense, we can then sample an entirely new data set on which to perform estimation and inference (again, where possible, using out-of-sample assessment to avoid over-reliance on asymptotic approximations and on assuming the model is actually the data-generating process when computing thinks like intervals and standard errors). If we can't do something like that, such list trawling can be quite dangerous, leading to models that tell us something about the particulars of our data and our modelling process but little of value about the thing we want to model.

• Thanks for the explanation Glen. Indeed I underestimated the complexity of fitting a curve to this set of data. It seemed to be a straightforward case when I saw the scatter plot, but now I see that's not the case.Since I dont fully understand these formulations, I think I need to discuss the results further with my supervisors instead of hunting for a fit that I won't be able to discuss due to my inexperience. Once again, thanks. – Renata Dis Apr 20 '17 at 17:57
• A good idea, I think. A literature search may help identify one or two potential models (and perhaps likely flaws), so you may be able to go into such a discussion with more context than you currently have. – Glen_b Apr 20 '17 at 22:44

I fit the ASCII text data you posted to over 600 known, named equations with four or less parameters using a genetic algorithm for initial parameter estimation, and the best candidates were mostly hyperbolic type equations. Here is an example that I consider to be a good candidate equation:

y = x / (a + bx) + Offset

a =  6.8885679248813901E+05
b = -2.3699198001685570E+05
Offset =  6.6651456685672019E+03

Degrees of freedom (error): 7
Degrees of freedom (regression): 2
Chi-squared: 3.09776595618e-05
R-squared: 0.964515690092
Model F-statistic: 95.1351435065
Model F-statistic p-value: 8.41641753468e-06
Model log-likelihood: 49.2347614958
AIC: -9.24695229916
BIC: -9.15617677126
Root Mean Squared Error (RMSE): 0.00176004714601

a = 6.8885679248813901E+05
std err: 2.86714E+10
t-stat: 4.06822E+00
p-stat: 4.75990E-03
95% confidence intervals: [2.88464E+05, 1.08925E+06]

b = -2.3699198001685570E+05
std err: 3.39266E+09
t-stat: -4.06877E+00
p-stat: 4.75659E-03
95% confidence intervals: [-3.74723E+05, -9.92607E+04]

Offset = 6.6651456685672019E+03
std err: 2.13899E-06
t-stat: 4.55728E+06
p-stat: 0.00000E+00
95% confidence intervals: [6.66514E+03, 6.66515E+03]

Coefficient Covariance Matrix
[  6.47885451e+15  -2.22866161e+15  -4.28731703e+07]
[ -2.22866161e+15   7.66637461e+14   1.47483491e+07]
[ -4.28731703e+07   1.47483491e+07   4.83345536e-01]

• What the output's account of degrees of freedom misses is the "researcher degrees of freedom" used up in the search for a model across 600 options. The p-values here don't carry their usual properties or meaning (the same goes for standard errors and confidence intervals and a number of other quantities) -- as a simulation study would readily show. Indeed by ignoring the effect of the search on their properties, almost all of the information in the output is biased or misleading in one way or another (at least if interpreted in the usual way). ... the downvote wasn't mine, by the way – Glen_b Apr 21 '17 at 1:44