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 )