4
$\begingroup$

I set up a kinetic model consisting of a system of ODE's. These ODE's are solved with an ode-solver (ode45) for varying parameters $k_1$ & $k_2$.

The parameters are estimated using lsqnonlin (for nonlinear least-squares (nonlinear data-fitting) problems) which minimizes the "difference" between experimental and model data. The dataset consists of 180 observations from 6 experiments.

To evaluate model quality with the resulting $k_1$ & $k_2$, I would like to calculate the Akaike Information criterion (AIC) (to compare to other models with more or less parameters). For AIC I use the following formula:

$$AIC = 2p - 2\ln(L)$$

where $p$ is the number of parameters and $\ln(L)$ is the maximum log-likelihood of estimated model.

From Spiess & Neumeyer: An evaluation of R² as an inadequate measure for nonlinear models in pharmacological and biochemical research I got the following formula for $\ln(L)$:

lnL = 0.5 * (-nRes * (log(2*pi) + 1 - log(nRes) + log(sum((yExp(:)-yModel(:)).^2))));

nRes - number of observations,
yExp - experimantal data,
yModel - model data.

And here are my questions (especially Q3-Q5 in EDITs unanswered):

  1. Is this formula for max log-likelihood okay to use in this context? I feel like the "punishment" for more parameters is quite weak if I'm getting AIC values of e.g. -2343 and -2199 for a model with 3 and 2 parameters, respectively.

  2. Should I look at a different parameter than AIC to compare the models?

Edit 1: Trying to give a bit more context on the data.

6 Experiments were conducted. During each, 10 samples at different time points were taken with the first one being the starting conditions. In each of the 60 samples, the concentration of 3 species were measured, leading to the total of 180 observations. Different Kinetic models for the concentrations of these species are fitted to the data for different kinetic parameters.

Question 3: Is the 180 correct to use as nRes?

Edit 2: Regarding the discussion with measure_theory on why log(nRes) is in the formula of $\ln(L)$.

The original formula would be:

lnL = 0.5 * (-nRes * (log(2*pi) + 1 + log(σ_e²)));

with Maximum likelihood estimator σ_e²

σ_e²= RSS/nRes

and

log(RSS/nRes) = log(RSS) - log(nRes)

So I think the formula should be okay, if the normality assumption for the residuals is correct, which I hope it is...

Question 4: Does someone have any input for question 3 in Edit 1?

Question 5: I'm still wondering about the number of parameters that I have to use for $p$. Of course all parameters estimated in the model are counted, but I read I have to count σ². Is this correct? (number of estimated parameters +1?)

$\endgroup$
2
$\begingroup$

The log-likelihood equation you write is (almost) the log-likelihood equation for a normally distributed distribution with $\sigma = 1$.

The only thing that's throwing me off is the $-\log(nRes)$ term... Is that an error?

What you're basically assuming is that the difference between your model, $f(x)$, and the actual, $y$, is normally distributed conditional on the observed data, $x$.

i.e. $y_t-f(x_t) = \epsilon_t\sim \mathcal{N(0,1)}$

This is probably a safe assumption to make, but probably requires a bit more context to make sure.

You're right that the punishment appears weak, based on the the values of AIC, but consider the fact that you're only estimating 2-3 parameters for 180 observations of data. That's a lot of degrees of freedom to work with! The purpose of the penalty is to avoid over-fitting, and a ratio of 2-3 parameters to 180 obs. is most likely far from that.

I think you're on the right track. All you're assuming is $y_t-f(x_t) = \epsilon_t\sim \mathcal{N(0,1)}$ for different functional model forms, $f(\cdot)$, and you're using the information contained in the likelihood function to guide your model selection process.

Edit:

The inclusion of $\log(nRes)$ is a bit odd to have in the likelihood function, since if you carry the $nRes$ outside the parenthesis through you get $\log(nRes^{nRes})$, which is probably contributing to why you're getting such high AIC values.

Nethertheless, as long as all the models you're comparing have the same value of $nRes$, the results will cancel out comparatively.

e.g.

$$L_a = 0.5[-nRes(\log(2\pi) + 1 - \log(nRes) + \sum_t(y_t - f_a(x_t))^2)]= C\sum_t(y_t - f_a(x_t))^2$$

And taking the ratio of any other potential model $f_b$, the constant term's will cancel.

e.g. $$\frac{L_a}{L_b} = \frac{C\sum_t(y_t - f_a(x_t))^2}{C\sum_t(y_t - f_b(x_t))^2} = \frac{\sum_t(y_t - f_a(x_t))^2}{\sum_t(y_t - f_b(x_t))^2} $$

$\endgroup$
  • $\begingroup$ Thank you very much for your answer. I added some more context on my data in the edit of my initial question. The - log(nRes) term was no error according to the publication I got the formula from. $\endgroup$ – S. Michael May 8 '18 at 16:04
  • $\begingroup$ Thanks again for your input! So as long as nRes stays constant, I could pretty much just use the sum-term and forget the rest, if I only want to compare the models. Which would be almost like using . So maybe I should look for a different log-likelihood function. Regarding my edit in the original post and the added question 3, do you think the assumption of normally distributed residuals is correct? $\endgroup$ – S. Michael May 22 '18 at 17:35

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.