How to determine the distribution of a parameter fit by nonlinear regression

Question

The example above shows enzyme kinetics -- enzyme velocity as a function of substrate concentration. The well-established Michaelis-Menten equation is:

$Y=V_{max} \cdot \dfrac{X}{K_m + X}$

• $X$ are the concentrations of substrate (set by the experimenter)
• $Y$ are the enzyme activities (measured by the experimenter)
• $V_{max}$ is the maximum enzyme velocity at high substrate concentrations. It is fit by nonlinear regression. It has the same units as Y and must be positive.
• $K_m$ is the Michaelis-Menten constant, which is the substrate concentration that leads to half-maximal velocity. Since it is a concentration, it must be positive. It is fit by nonlinear regression, and has the same units as X.

The left panel shows one simulated data set. $V_{max}$ was set to 84 and $K_m$ was set to 4. Each $Y$ value was computed from the equation above plus a random error (Gaussian, SD=12). I made the SD high to make the variation in $K_m$ pronounced. The curve was fit by nonlinear regression using the equation above to determine the $V_{max}$ and $K_m$. Since the residuals are assumed to be Gaussian (and for this example were simulated that way), the nonlinear regression minimizes the sum of the squared residuals.

The middle panel shows the values of $K_m$ fit by nonlinear regression for 100 such simulations. The asymmetry is clear.

The right panel shows the frequency distribution of $K_m$ determined from 10,000 simulations. The distribution was fit to both a normal distribution (red; fits poorly) and a lognormal distribution (blue; fits well). I think this demonstrates pretty clearly that the distribution of $K_m$ is lognormal, or at least it follows a distribution very similar to the lognormal distribution.

My questions are:

• For this example, can algebra and/or calculus prove that the distribution of $K_m$ values is lognormal (or prove it has some other distribution)?
• More generally, what method can be used to derive the distribution of any parameter fit by nonlinear regression?

Answer 1

This answer does not (yet) answer the question but should at least help to clarify what the question really is:

---

"fit by nonlinear regression" sounds like you are using the following model:

$\mathcal{Y}\sim \mathcal{N}(\mu=\frac{X}{X+K_m}, \sigma^2)$

(this assumes that there is no error in measuring the substrate concentration X; If this nevertheless a good model is another question)

The corresponding likelihood function given a sample $Y^N$ is:

$p_{\mathcal{Y^N}}(Y^N|K_m, \sigma, X^N) = \prod_{i=1}^Np_{\mathcal{N}}(Y^N|\mu=\frac{X^N_i}{X^N_i+K_m}, \sigma^2)$,

where $p_\mathcal{N}$ is the density of the normal.

and sounds like you are using maximum-likelihood to estimate $K_m$ (and $\sigma^2$).

(if this is a good approach is yet another question)

$ML_{\hat{K_m}}(X^N,Y^N) = \operatorname*{argmax}\limits_{K_m} \operatorname*{max}\limits_{\sigma} p_{\mathcal{Y^N}}(Y^N|K_m, \sigma, X^N)$

You then seem to sample $\mathcal{Y^N}$ for some fixed $X^N$, $K_m$ and $\sigma$

(Where $X^N$ might be your data while $K_m$ and $\sigma$ might be the estimate you obtained for your data with above ML approach)

and then apply above ML estimator ( let's call it ), thus sampling a random variable $\mathcal{\hat{K_m}}$ whose distribution you are asking about (and which you are plotting). There are legit reasons to desire an explicit form of this distribution; for example, to construct confidence intervals for your estimation of $K_m$.

However since this distribution is not (symmetric and uni-modal) it's yet another question which is the best way to construct a confidence interval given this distribution

Note, however, that this distribution is NOT the posterior distribution of nor a likelihood function for $K_m$ and thus probably not what you desired when you said "the distribution of a parameter".

the likelihood function is trivial to obtain (look at logLik for your model in R) while the posterior requires you to choose a prior (the empirical distribution of $K_m$ values in databases might be a good choice)

Anyway, let's see how far we get. Let's start by expressing it as compound distribution using the the distribution of $Y^N$ that we know:

$p_{\mathcal{\hat{K_m}}} (\hat{K_M})=\int_{ \{Y^N|\hat{K_M}=ML_{\hat{K_m}}(X^N,Y^N)\}} p_{\mathcal{Y^N}}(Y^N) \mathrm{d} Y^N$

This contains $ML_{\hat{K_m}}(X^N,Y^N)$ for which we might be able to find and algebraic expression for: $ML_{\hat{K_m}}(X^N,Y^N) = \operatorname*{argmax}\limits_{K_m} \operatorname*{max}\limits_{\sigma} \prod_{i=1}^Np_{\mathcal{N}}(Y^N_i|\mu=\frac{X^N_i}{X^N_i+K_m}, \sigma^2)$

$ = \operatorname*{argmax}\limits_{K_m} \operatorname*{max}\limits_{\sigma} \sum_{i=1}^N\log(p_{\mathcal{N}}(Y^N_i|\mu=\frac{X^N_i}{X^N_i+K_m}, \sigma^2))$

$ = \operatorname*{argmax}\limits_{K_m} \operatorname*{max}\limits_{\sigma} \sum_{i=1}^N\log(\frac{1}{\sqrt{2\pi\sigma^2}}) - \frac{\left(Y^N_i-\frac{X^N_i}{X^N_i+K_m}\right)^2}{2\sigma^2}$

$ = \operatorname*{argmin}\limits_{K_m} \sum_{i=1}^N \left(Y^N_i-\frac{X^N_i}{X^N_i+K_m}\right)^2$

$ 0 = \left.\frac{\mathrm{d}}{\mathrm{d} K_m} \sum_{i=1}^N \left(Y^N_i-\frac{X^N_i}{X^N_i+K_m}\right)^2\right|_\hat{K_m}$ $ = \sum_{i=1}^N \left.\frac{\mathrm{d}}{\mathrm{d} K_m} \left(Y^N_i-\frac{X^N_i}{X^N_i+K_m}\right)^2\right|_\hat{K_m}$ $ = \sum_{i=1}^N \frac{X^N_i(\hat{K_m}Y^N_i+X^N_i(Y^N_i-1))}{(\hat{K_m}+X^N_i)^3}$

From where I don't know how to continue.

---

I'm still in the progress of refining this answer please find below a current draft to decide if it's worth your bounty:

In this answer I assume $V_{max}$ is known to be (without loss of generality) 1. As confirmed in the comments you are using the following model:

$\mathcal{Y}\sim \mathcal{N}(\mu=\frac{X}{X+K_m}, \sigma^2)$

The corresponding likelihood function is

$L(K_m, \sigma) = p_{\mathcal{Y^N}}(Y^N|K_m, \sigma, X^N) = \prod_{i=1}^Np_{\mathcal{N}}(Y^N|\mu=\frac{X^N_i}{X^N_i+K_m}, \sigma^2)$,

where $p_\mathcal{N}$ is the density of the normal distribution.

Now, you would like to know the distribution of a random variable $\mathcal{\hat{K_m}}$ that is the maximum likelihood estimate,

$ML_{\hat{K_m}}(X^N,Y^N) = \operatorname*{argmax}\limits_{K_m} \operatorname*{max}\limits_{\sigma} p_{\mathcal{Y^N}}(Y^N|K_m, \sigma, X^N)$ $ = \operatorname*{argmax}\limits_{K_m} \operatorname*{max}\limits_{\sigma} \prod_{i=1}^Np_{\mathcal{N}}(Y^N_i|\mu=\frac{X^N_i}{X^N_i+K_m}, \sigma^2)$

$ = \operatorname*{argmax}\limits_{K_m} \operatorname*{max}\limits_{\sigma} \sum_{i=1}^N\log(p_{\mathcal{N}}(Y^N_i|\mu=\frac{X^N_i}{X^N_i+K_m}, \sigma^2))$

$ = \operatorname*{argmax}\limits_{K_m} \operatorname*{max}\limits_{\sigma} \sum_{i=1}^N\log(\frac{1}{\sqrt{2\pi\sigma^2}}) - \frac{\left(Y^N_i-\frac{X^N_i}{X^N_i+K_m}\right)^2}{2\sigma^2}$

$ = \operatorname*{argmin}\limits_{K_m} \sum_{i=1}^N \left(Y^N_i-\frac{X^N_i}{X^N_i+K_m}\right)^2$,

obtained for draws of draws of size $N$ from $\mathcal{Y}$, $\mathcal{Y^N}$, for any $N$, $X^N$, $\sigma$.

You then sampled $K_m$ for some fixed $K$, $X^N$, $K_m$ and $\sigma$ by first sampling $\mathcal{Y^N}$ accordingly and then applying above ML estimator. Based on this, you think that $\mathcal{K_m}$ follows a log normal distribution.

It is known that, for any differentiable function $f: \mathbb{R}^N \to \mathbb{R}$ and $\mathcal{Y} = f(\mathcal{X})$,
$p_\mathcal{Y}(y) = \int_x \delta(f(x)-y) p_\mathcal{X}(x)\mathrm{d}x$ , where $\delta$ is the Dirac delta.

And that for any monotonic function $g: \mathbb{R} \to \mathbb{R}$ and $\mathcal{Y} = f(\mathcal{X})$,
$p_\mathcal{Y}(y) = p_\mathcal{X}(g^{-1}(y)) \left|\frac{\mathrm{d}}{\mathrm{d}y} g^{-1}(y) \right|$

We can use this to try to derive a closed form for the density of the distribution of $\mathcal{\hat{K_m}}$:

$p_{\mathcal{\hat{K_m}}}(\hat{K_m})=\int \delta (\hat{K_m}-ML_{\hat{K_m}}(X^N,Y^N)) p_{\mathcal{Y^N}}(Y^N) \mathrm{d} Y^N$

$\overset{\tiny{\text{if i'm lucky}}}{=}\int \delta(\frac{\mathrm{d}}{\mathrm{d} \hat{K_m}} \sum_{i=1}^N \left(Y^N_i-\frac{X^N_i}{X^N_i+\hat{K_m}}\right)^2) p_{\mathcal{Y^N}}(Y^N) \mathrm{d} Y^N$

$=\int \delta(\sum_{i=1}^N \frac{X^N_i(\hat{K_m}Y^N_i+X^N_i(Y^N_i-1))}{(\hat{K_m}+X^N_i)^3}) p_{\mathcal{Y^N}}(Y^N) \mathrm{d} Y^N$

But I don't how to find a simpler form for that.

For $N=1$ this is a bit simpler:

$p_{\mathcal{\hat{K_m}}}(\hat{K_m})=p_\mathcal{Y}(g^{-1}(\hat{K_m})) \left|\frac{\mathrm{d}}{\mathrm{d}\hat{K_m}} g^{-1}(\hat{K_m}) \right| = p_\mathcal{Y}(\frac{X}{X+\hat{K_m}}) \left|\frac{\mathrm{d}}{\mathrm{d}\hat{K_m}} \frac{X}{X+\hat{K_m}} \right|= p_\mathcal{Y}(\frac{X}{X+\hat{K_m}}) \left|- \frac{X}{(X+\hat{K_m})^2} \right|= p_{\mathcal{N}}(\frac{X}{X+\hat{K_m}}|\mu=\frac{X^N_i}{X^N_i+K_m}, \sigma^2) \frac{X}{(X+\hat{K_m})^2} $

Where I used: $ML_{\hat{K_m}}(X^N,Y^N) = \operatorname*{argmin}\limits_{K_m}\left(y-\frac{x}{x+K_m}\right)^2 \Leftrightarrow 0 =\frac{x(\hat{K_m}y+x(y-1))}{(\hat{K_m}+x)^3} \land (\text{further conditions})$ which solves $\hat{K_m}=x(\frac{1}{y}-1)$.

For $N=2$ the explicit form of $ML_{K_m}$ has quite a few more terms

In any case, this shows that $p_{\mathcal{\hat{K_m}}}(\hat{K_m})$ is not log normal (but might converge to it (before converging to normal)).

Answer 2

My questions are:

• For this example, can algebra and/or calculus prove that the distribution of Km values is lognormal (or prove it has some other distribution)?
• More generally, what method can be used to derive the distribution of any parameter fit by nonlinear regression?

• The Km values can not be exactly lognormal. This is because in your problem formulation negative values can occur as the maximum likelihood estimate (yes the negative values do not make sense, but neither do normal distributed errors, which can cause the negative Km values). Of course, the lognormal might still be a reasonable approximation.

  A more rigorous 'proof' that the distribution can not be exactly lognormal is given below in the special case with measurements in two points. For that case it is possible/easy to compute the estimates explicitly and express the sample distribution of the estimates.

• Below I describe a method that derives an approximate distribution by not performing a normal approximation to the $K_m$ parameter directly, but instead on two other parameters from which a different approximate sample distribution for $K_m$ is derived.

  The second part in the following, improving it even more, is very experimental. It shows a very reasonable fit, but I do not have a proof for it. I have to look into that further. But I thought it was interesting to share.

1 Different parameterization

I can re-express the Michaelis-Menten equation as a generalized linear model (using the normal family with inverse as link function):

$$y \sim N\left( \frac{1}{\beta_0+\beta_1 z},\sigma^2 \right)$$

Where

• $z = 1/x$ the inverse of your variable $x$ for the substrate concentrate
• $\beta_0 = 1/V_{max}$ the inverse of your enzyme velocity parameter
• $\beta_1 = K_m/V_{max}$ the ratio of your half-maximal and velocity parameters

The parameters $\beta_i$ will be approximately multivariate normal distributed. Then the distribution of $K_m = \beta_1/\beta_0$ is the ratio of two correlated normal variables.

When we compute this then we get a slightly more reasonable fit

set.seed(1)

### parameters
a = 10
b = 5
n <- 10^5

### two arrays of sample distribution of parameters
am <- rep(0,n)
bm <- rep(0,n)

### perform n times a computation to view te sample distribution
for (i in 1:n) {
  x <-seq(0,40,5)
  y <- a*x/(x+b)+rnorm(length(x),0,1)
  mod <- nls(y ~ ae * x/(x+be), start = list(ae=a,be=b))
  am[i] <- coef(mod)[1]
  bm[i] <- coef(mod)[2]
}

### histogram
hist(bm, breaks = seq(-2,30,0.3), freq = 0 , xlim = c(0,20), ylim = c(0,0.20),
     main = "histogram compared with \n two normal approximations",
     xlab = "Km", cex.main = 1)

### fit with normal approximation
s <- seq(0,22,0.01)
lines(s,dnorm(s,mean(bm),var(bm)^0.5))

### fit with ratio of normal approximation
w <- fw(s,mean(bm/am),mean(1/am),var(bm/am)^0.5,var(1/am)^0.5,cor(1/am,bm/am))
lines(s,w,col=2)

legend(20,0.20,
       c("normal approximation",
         "normal ratio approximation"),
       xjust = 1, cex = 0.7, col = c(1,2), lty = 1 )

Here we used the following function to compute the ratio of two correlated normal distributions (see also here). It is based on: Hinkley D.V., 1969, On the Ratio of Two Correlated Normal Random Variables, Biometrica vol. 56 no. 3.

## X1/X2      
fw <- function(w,mu1,mu2,sig1,sig2,rho) {
  #several parameters
  aw <- sqrt(w^2/sig1^2 - 2*rho*w/(sig1*sig2) + 1/sig2^2)
  bw <- w*mu1/sig1^2 - rho*(mu1+mu2*w)/(sig1*sig2)+ mu2/sig2^2
  c <- mu1^2/sig1^2 - 2 * rho * mu1 * mu2 / (sig1*sig2) + mu2^2/sig2^2
  dw <- exp((bw^2 - c*aw^2)/(2*(1-rho^2)*aw^2))
  
  # output from Hinkley's density formula
  out <- (bw*dw / ( sqrt(2*pi) * sig1 * sig2 * aw^3)) * (pnorm(bw/aw/sqrt(1-rho^2),0,1) - pnorm(-bw/aw/sqrt(1-rho^2),0,1)) + 
    sqrt(1-rho^2)/(pi*sig1*sig2*aw^2) * exp(-c/(2*(1-rho^2)))
  
  out
}
fw <- Vectorize(fw)

In the above computation, we estimated the covariance matrix for the sample distribution of the parameters $\beta_0$ and $\beta_1$ by simulating many samples. In practice, when you only have a single sample, you could be using an estimate of the variance based on the observed information matrix (for instance when you use in R the glm function, then you can obtain estimates for the covariance, based on the observed information matrix by using the vcov function) .

2 Improving normal approximation for parameter $\beta_1$

The above result, using $K_m = \beta_1/\beta_0$ is still not great because the normal approximation for the parameter $\beta_1$ is not perfect. However, with some trial and error, I found that a scaled noncentral t-distribution is a very good fit (I have some intuitive idea about it but I can not yet explain so well why, let alone proof it).

h <- hist(bm/am, breaks = seq(-2,3,0.02), freq = 0 , xlim = c(-0.2,1.3), ylim = c(0,3),
     main = "histogram compared with normal and t-distribution",
     xlab = expression(beta[1]), cex.main = 1)

### fitting a normal distribution
s <- seq(0,22,0.001)
lines(s,dnorm(s,mean(bm/am),var(bm/am)^0.5))

### fitting a t-distribution to the histogram
xw <- h$mids
yw <- h$density
wfit <- nls(yw ~ dt(xw*a, df, ncp)*a, start = list(a=2,df=1, ncp = 0.5),
            control = nls.control(tol = 10^-5, maxiter = 10^5),
            algorithm = 'port',
            lower = c(0.1,0.1,0.1))
wfit
lines(xw,predict(wfit),col = 2)

legend(1.3,3,
       c("normal approximation",
         "t-distribution approximation"),
       xjust = 1, cex = 0.7, col = c(1,2), lty = 1 )

Special case with measurements in two points

If you measure in only two points $x=s$ and $x = t$, then you could reparameterize the curve in terms of the values in those two points $y(s)$ and $y(t)$. The parameter $K_m$ will be

$$K_m = \frac{y(t)-y(s)}{y(s)/s-y(t)/t}$$

Since estimates of $y(t)$ and $y(s)$ will be independent and normally distributed the sample distribution of the estimate of $K_m$ will be the ratio of two correlated normal distributions.

The computation below illustrates this with a perfect match.

The fit with a lognormal distribution is actually not so bad either (and I needed to use some extreme parameters to make the difference clearly visible). There might be a connection between a product/ratio distribution and the lognormal distribution. It is similar to this question/answer where you have a variable that is a product of several terms. This is the same as the exponent of the sum of the log of those terms. That sum might be approximately normal distributed if either you have a lot of terms or when you have a few terms that are already approximately normal distributed.

$$K_m = e^{\log(K_m/V_{max}) - \log(1/V_{max})}$$

set.seed(1)

### parameters
a = 50
b = 5
n <- 10^5
t = 2
s = 4

### two arrays of sample distribution of parameters
am <- rep(0,n)
bm <- rep(0,n)

### perform n times a computation to view the sample distribution
x <- c(t,s)
for (i in 1:n) {
  y <- a*x/(x+b)+rnorm(length(x),0,1)
  mod <- lm(1/y ~ 1+I(1/x))
  am[i] <- 1/coef(mod)[1]
  bm[i] <- coef(mod)[2]/coef(mod)[1]
}

### histogram
h <- hist(bm, breaks = c(-10^5,seq(-100,100,0.2),10^5), freq = 0 , xlim = c(0,15), ylim = c(0,0.30),
          main = "special case of measurement in two points",
          xlab = "Km", cex.main = 1)

### plotting fit with lognormal distribution
xw <- h$mids
yw <- h$density
wfit <- nls(yw ~ dlnorm(xw, mu, sd), start = list(mu = log(5), sd = 0.5),
            control = nls.control(tol = 10^-5, maxiter = 10^5),
            algorithm = 'port',
            lower = c(0.1,0.1))
wfit
lines(xw,predict(wfit),col = 1)


### plotting ratio distribution
### means, sigma and distribution
y1 = a*s/(b+s)
y2 = a*t/(b+t)
cc = -(1/s + 1/t)/sqrt(1+1)/sqrt(1/t^2+1/s^2)

lines(ts,fw(ts, mu1 = y2-y1   , 
            mu2 = y1/s-y2/t,
            sig1 = sqrt(1+1), 
            sig2 = sqrt(1/t^2+1/s^2),
            rho  = cc  ),
      col  = 2)  


legend(15,0.3,
       c("ratio distribution", "fit with lognormal"),
       xjust = 1, cex = 0.7, col = c(2,1), lty = 1 )