Fitting a curve to data: inconsistent results

I have a dataset of nutrient budgets for a set of lakes. Values are compiled from the literature and I've made the data available online.

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from lmfit import Model

data_url = r"https://gist.githubusercontent.com/JamesSample/029f59471818929ef9fb87bde95169bc/raw/c5b66b554b9ff9d3c448a80c0b4104f410ab5185/lake_n_budgets.csv"
df = pd.read_csv(data_url)
df.head()

Lin is the flux of nutrients into the lake (tonnes/year); Lout is the flux from the outflow (tonnes/year); and tau is the water "residence time" (in years) i.e. the average time water spends circulating in the lake.

According to simple theory, these variables can be related as

where sigma and n are empirical constants that I'd like to find.

Attempt 1

I had a go at fitting this using lmfit

def load_out(Lin, tau, sigma=1, n=1):
    return Lin / (1 + (sigma * (tau**n)))


model = Model(load_out, independent_vars=["Lin", "tau"])
fit = model.fit(df["Lout"], Lin=df["Lin"], tau=df["tau"])
fit

which gives the following results

So far, so good.

However, in the literature it's common to think in terms of the "transmission factor" of the lakes, which is defined as

Combining this with the first equation gives

Attempt 2

As a "sense check" while learning to use lmfit, I had a go at fitting this version too

def transmission(tau, sigma=1, n=1):
    return 1 / (1 + (sigma * (tau**n)))


df["trans"] = df["Lout"] / df["Lin"]
model = Model(transmission, independent_vars=["tau"])
fit = model.fit(df["trans"], tau=df["tau"])
fit

The results are quite different: the best estimates for sigma and n are roughly twice what they were before, and the differences are larger than the uncertainty estimates on the parameters. Since this is essentially the same equation - with the same sigma and n - I'm a bit surprised by this.

Attempt 3

As a further check, I tried rearranging the equation, taking logs and then fitting a linear regression using statsmodels

df["y"] = (df["Lin"] / df["Lout"]) - 1
mod = smf.ols(formula='np.log10(y) ~ np.log10(tau)', data=df)
res = mod.fit()
print(res.summary())
print()
print("Estimate for sigma:", 10**res.params[0])
print("Estimate for n:", res.params[1])

This gives

Estimate for sigma: 0.757
Estimate for n: 0.361

which is broadly compatible with the results from lmfit in Attempt 2.

My question(s)

I'd like to understand why attempt 1 gives such different results, and the implications for my analysis.

• Am I missing something really obvious, or is it reasonable to be surprised by these differences?

• Should the uncertainty estimates on the parameters not be larger than the differences between simple rearrangements of the fitted equation?

If I create a synthetic dataset where the relationship matches exactly, all three approaches give identical results, so I think my code is working correctly. I therefore think the differences are due to characteristics/uncertainties in my dataset, but I'm still puzzled over why the difference is so large.

• Is there a bug in my lmfit code for attempt 1?
• Is attempt 1 fundamentally wrong in some way?

For attempt 1, I suppose I'm fitting a surface to a set of points in three dimensions, whereas in the other versions I'm fitting a curve in two dimensions. The latter is presumably easier/more robust with a small dataset like this (~170 data points), so I think attempts 2 and 3 probably give the "best" answer.

• Does the fact that attempt 1 is so different imply this model is just not appropriate for my dataset?

Thank you! :-)