8
$\begingroup$

Suppose I'm working with the following model

$y_i = \alpha(1-\exp(-\beta t_i))+\gamma(1-\exp(-\delta t_i)) + \varepsilon_i$.

The $\varepsilon_i$ are i.i.d. gaussian with zero mean and I'm trying to find best fit values of $\alpha,\beta,\gamma,\delta$.

For concreteness, say this is a model for the total amount of some bacterial species with two subspecies who grow in time according to the first and second terms on the RHS, but we only measure the total population. Note: this is not the actual setting, but suffices for the question.

The model is not identifiable in the usual sense since I can always just swap $\alpha$ and $\gamma$ for example and get exactly the same density/likelihood.

As you might expect, if I run an MCMC on this I end up with horribly wide posteriors, and any nonlinear least squares approach is incredibly sensitive to initial guesses -- we have a big plateau in the likelihood function.

Better experimental design is not an option at this stage -- obviously measuring the subspecies separately would be the best option.

Is there anything I can do with this problem, or is better experimental design the only option?

$\endgroup$
9
  • 1
    $\begingroup$ Why do you have two terms with the same structure? As you said, they are misidentified. $\endgroup$ – Anna SdTC Feb 5 '17 at 3:34
  • 1
    $\begingroup$ They represent two 'sub-quantities' of the measured quantity. I did not formulate this model, I was given it to try and fit to some data and the people who gave it to me said they don't want to change it. I wanted to make sure there was nothing I could do before telling them it was not possible. $\endgroup$ – Will Feb 5 '17 at 3:45
  • 1
    $\begingroup$ I don't think you can based on this model. Can you do some sort of two-step estimation in which in the first step there is only one of these exponential terms and in the second step you compute the "sub-quantitities"? $\endgroup$ – Anna SdTC Feb 5 '17 at 3:48
  • 1
    $\begingroup$ That might work. I definitely need to fix some parameters if I'm going to get anywhere at all and that might be a good way to do it. What if I put an exponential prior (with 'big' rate parameter) on either $\alpha$ or $\gamma$ (they are non-negative)? Thinking along the lines of trying to regularise somehow... $\endgroup$ – Will Feb 5 '17 at 3:58
  • 1
    $\begingroup$ For general information, as nonlinear least-squares methods there is "Levenberg-Marquardt" (LM) and the "trust-region reflective" (TRR), have you tried one of them? I know that the TRR allows you to set boundaries on parameters whereas LM don't. This said, for your problem, I think you have to have something more that allows you to distinguish both terms. Either you can identify measurements with only one population and do a two-step estimation as Anna noted, or change your model to account for the specificities of each subspecies. Otherwise I don't see the usefulness of this model. $\endgroup$ – Learn_and_Share Feb 5 '17 at 11:43
5
$\begingroup$

There is no identifiability problem, except in the trivial sense that any one particular model can have two descriptions. The real problem appears to be the difficulty in fitting the model--but that's due to how the models are parameterized rather than lack of identifiability.

This problem has an equally trivial solution: declare, without any loss of generality, that $\beta \ge \delta$. If you want to be really fussy, also insist that if $\beta=\delta$, then $\alpha \ge \gamma$.

Unfortunately, that requires any procedure to fit the model to respect these constraints. Introducing a constraint here is not so bad though, because the application is such that obviously all parameters are non-negative anyway: the parameter space already has sharp boundaries. Including one more constraint won't force any changes in how we go about fitting the model.


One well-known method to convert a constrained optimization into an unconstrained one is to reparameterize the problem so that in the new parameter space the boundaries are pushed out to infinity. There are many ways to accomplish that here. A consideration of what the parameters mean will guide us. In particular, $\nu = \alpha + \gamma$ is the maximum attained by the function $$t\to g(t; \alpha,\beta,\gamma,\delta) = \alpha\left(1 - e^{-\beta t}\right) + \gamma\left(1 - e^{-\delta t}\right)$$ for $t \ge 0$. Given $\nu$, then necessarily $0 \le \alpha\le \nu$ and $\gamma = \nu - \alpha$. When non-negative values sum to a fixed whole, it often works to parameterize their proportions of the whole in terms of angles: let one proportion be the squared cosine and the other be the squared sine. Furthermore, a simple way to ensure $\nu$, $\beta$, and $\delta$ are positive is to make them exponentials--that is, use their logarithms as parameters. Finally, to enforce $\delta \le \beta$, set $\delta$ to be the squared cosine of some angle times $\beta$. Thus we might reparameterize the problem by fitting the function

$$t \to f(t;n,a,b,d) = e^n\left(1 - \cos(a)^2\exp\left(-e^{b} t\right) - \sin(a)^2\exp\left(-e^{b\cos(d)^2} t\right)\right).$$

From estimates of these parameters (which, by the way, are not "identifiable" due to the ambiguity in the angles $a$ and $d$) you can recover the original ones as

$$\eqalign{ \alpha &= e^n\cos(a)^2 \\ \beta &=e^b \\ \gamma &=e^n\sin(a)^2 \\ \delta &= e^b\cos(d)^2. }$$

Properties of the exponential and trig functions assure all the constraints hold: $\alpha \gt 0$, $\beta \ge \delta \gt 0$, and $\gamma \gt 0$. (Since double precision floats can become astronomically small, there is no practical distinction between $\gt$ and $\ge$ in these constraints.)

In this well-defined sense the model is identifiable even though the parameters used to fit it are not identifiable.


Although one could use MCMC, if the purpose is just to fit the curve then it's more straightforward to use a numerical solver such as Newton-Raphson. The trick is to find a good starting value. The maximum of the $y_i$ would be a slight overestimate of $e^n$; so start, perhaps with $n=\log(\max(y_i)/2)$. You might begin with $a=\pi/4$, supposing each component makes a substantial contribution to the whole. Make some reasonable guesses about $e^b$ and $e^d$ based on expected decay rates. For instance, if the range of $t$ is reasonable, then take $b$ to be some fraction of the largest $t$ and perhaps arbitrarily pick $d=\pi/4$; maybe use a smaller starting value. (You will often get different values for the parameter estimates depending on these choices, but typically they will not appreciably affect the function $f$ itself.)

In many circumstances, this approach works strikingly well. Except when the variance of the errors is the same size as $\max{y_i}$ or larger (where it will be hard to discern any signal at all without a large amount of data), the fit works even with tiny amounts of data: all it needs is four.

Note that regardless of how the model is fit, there ordinarily will be huge uncertainty in the parameters: this family of curves is essentially a tiny perturbation of the two-parameter exponential family $t\to Ae^{-Bt}$. In many circumstances, then, two of the parameters (corresponding to the amplitude $A$ and longest decay rate $B$) can be identified with reasonable precision but the other two parameters, which reflect small variations from this exponential shape, will usually be highly uncertain.

Figure

The figure shows an example of a challenging fit. The underlying curve is shown in black. Ultimately it reaches a maximum of $4/3$, very slowly. Only $24$ data points are available, plotted as gray dots. The standard deviation of the random errors is $1/2$, a sizable proportion of that maximum. Many of the errors were positive, causing the fitted curve in red to be a little higher. The two exponential components of the fitted curve are shown as dashed and dotted gray lines. One shows a rapid rise to a threshold of $1/3$ by the time $t=1$; the other reflects the other exponential rising to its threshold of $1$. (You will have little hope of reproducing that sharp "shoulder" near $t=1$ until you have a $1000$ data points or more: try it out by varying n in the code below.)

Your success in any particular problem will depend on the magnitudes of the errors; the range of values of $t$ that are sampled; how those values are spaced; how many values are available; and choice of starting values. Nevertheless this appears to be a tractable problem in general, with solutions that can be obtained rapidly. Moreover, any maximum likelihood fitter will proceed similarly to minimize the sum of squares of residuals--and will, in addition, provide confidence regions for the parameters.

This is the R code I used to test this proposal. It will reproduce the figure and is easily modified--change the values of the variables at the beginning--to study data that look like any you might have.

#
# Describe the underlying model
#
set.seed(17)
alpha <- 1
beta <- 2
gamma <- 1/3
delta <- 1/10
sigma <- 1/2  # Error SD.
n <- 24
x.max <- 20   # Largest value of t.
#
# The original parameterization.
#
g <- function(x, alpha, beta, gamma, delta) {
  alpha * (1 - exp(-beta * x)) + gamma * (1 - exp(-delta * x))
}
#
# The re-parameterization.  `f.1` and `f.2` are the two exponential components.
#
f <- function(x, nu, t.a, log.b, t.d) {
  n <- exp(nu)
  a <- cos(t.a)^2
  alpha <- n*a
  gamma <- n*(1-a)
  beta <- exp(log.b)
  delta <- cos(t.d)^2 * beta
  n - alpha * exp(-beta * x) - gamma * exp(-delta * x)
}
f.1 <- function(x, nu, t.a, log.b, t.d) {
  n <- exp(nu)
  a <- cos(t.a)^2
  alpha <- n*a
  beta <- exp(log.b)
  alpha * (1 - exp(-beta * x))
}
f.2 <- function(x, nu, t.a, log.b, t.d) {
  n <- exp(nu)
  a <- cos(t.a)^2
  gamma <- n*(1-a)
  beta <- exp(log.b)
  delta <- cos(t.d)^2 * beta
  gamma * (1 - exp(-delta * x))
}
#
# The objective to minimize is the mean squared residual.
# This is equivalent to finding the MLE for Gaussian errors.
#
obj <- function(theta, x, y) {
  crossprod(y - f(x, theta[1], theta[2], theta[3], theta[4])) / length(x)
}
#
# Create data and plot them.
#
x <- seq(0, x.max, length.out=n)
y <- g(x, alpha, beta, gamma, delta) + rnorm(length(x), 0, sigma)
plot(x,y, pch=16, col="#00000040", xlab="t")
#
# Fit the curve.
#
theta <- c(nu=log(max(y)/2), t.a=pi/4, log.b=log(max(x)/10), t.d=pi/4)
fit <- nlm(obj, theta, x=x, y=y, gradtol=1e-14)
theta.hat <- fit$estimate
#
# Plot relevant curves.
#
curve(g(x, alpha, beta, gamma, delta), add=TRUE, lwd=2)
curve(f(x, theta.hat[1], theta.hat[2], theta.hat[3], theta.hat[4]),
      add=TRUE, col="Red", lwd=2)
curve(f.1(x, theta.hat[1], theta.hat[2], theta.hat[3], theta.hat[4]),
      add=TRUE, col="Gray", lty=2, lwd=2)
curve(f.2(x, theta.hat[1], theta.hat[2], theta.hat[3], theta.hat[4]),
      add=TRUE, col="Gray", lty=3, lwd=2)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.