# Determination of maximum log-likelihood of nonlinear model for calculation of Aikaike IC

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?)

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}$$

• 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. – S. Michael May 8 '18 at 16:04
• 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? – S. Michael May 22 '18 at 17:35