## Sunday, July 7, 2013

### Estimating ODE's parameters

In a previous post I used R to solve Ordinary Differential Equations (ODE). Since I'm still more familiar with Matlab to work with ODE, I decided to go one step further and learn how to use R to estimate parameters in ODE.

In this short tutorial, I will use the ODE presented here which quantifies the salt concentration in a well-stirred tank:

\begin{align} &\frac{dx}{dt}=a\frac{150-x(t)}{200}\\[10pt] \end{align}

With $x(0) = 20$

The analytic solution is:

\begin{align} &x(t)=150-130e^{-at/200} \end{align}

Step 1

Generate data using the analytic solution and add some random noise. Note that the ODE uses one parameter that has been fixed at $$a = 0.75$$.

t = seq(1, 800, by = 10)
a = 0.75  ## Fixed parameter used to simulate data.

## Simulate data and add noise.
conc.observed = 150 - 130 * exp(-a * t/200) + rnorm(length(t), sd = 5)

## Plot
plot(t, conc.observed, pch = 21, bg = "gray", ylim = c(20, 180), xlab = "Time",
ylab = "Salt (kg)", axes = F)
box(bty = "l")
axis(1, tck = -0.02)
axis(2, tck = -0.02, las = 2)

legend("bottomright", legend = c("Observed data"), pch = 21, pt.bg = "gray",
bty = "n")


Step 2

Write a function that will be used to solve the ODE. While we there, solve it with fixed parameter $$a = 0.75$$.

library(deSolve)

salttank = function(t, state, parameters) {
with(as.list(c(state, parameters)), {

# rate of change
dX = a * ((150 - X)/200)

# return the rate of change
list(c(dX))
})
}

## Solve the ODE. Not necessary at this point.

## Define the initial value for the state variable
state = c(X = 20)

## Time range.
times = seq(1, 800, by = 10)

## Parameters
parameters = c(a = 0.75)

conc.modeled = ode(y = state, times = times, func = salttank, parms = parameters,
method = "rk4")


Plot the results.

plot(t, conc.observed, pch = 21, bg = "gray", ylim = c(20, 180), xlab = "Time",
ylab = "Salt (kg)", axes = F)
box(bty = "l")
axis(1, tck = -0.02)
axis(2, tck = -0.02, las = 2)
lines(conc.modeled[, "time"], conc.modeled[, "X"], col = "red")

legend("bottomright", legend = c("Observed data", "True solution (a = 0.75)"),
lty = c(NA, 1), pch = c(21, NA), col = c("black", "red"), pt.bg = "gray",
bty = "n")


Step 3

This is where the magic happens. To fit parameters, I will use nls.lm() from the minpack.lm package. In its simplest form, the function uses par which is a list of guessed parameters and fn a function used to minimize the sum square of the residuals using the Levenberg-Marquardt algorithm.

Now we have to write a function (here I named it ssq) that will use ODE parameters as input (in this case only $$a$$) and produces a residuals vector as output.

ssq = function(params) {

## Range and initial condition as before.
times = seq(1, 800, by = 10)
state = c(X = 20)

## Resolve the ODE.
out = ode(y = state, times = times, func = salttank, parms = params, method = "rk4")

## modeled - observed
ssq = out[, "X"] - conc.observed
}


Step 4

Finally, we can estimate the parameter $$a$$.


library(minpack.lm)

params.guessed = c(a = 1)
params.fitted = nls.lm(par = params.guessed, fn = ssq)

summary(params.fitted)

##
## Parameters:
##   Estimate Std. Error t value Pr(>|t|)
## a   0.7438     0.0113      66   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.93 on 79 degrees of freedom
## Number of iterations to termination: 5
## Reason for termination: Relative error in the sum of squares is at most ftol'.


We see from the summary of params.fitted that the estimated parameter is $$a = 0.7438$$ which is obviously very close to $$a=0.75$$ used to produce the observed data.

Step 5

Finally, simulate the data using the estimated parameter and plot the results.

## Simulate using fitted parameters
params = coef(params.fitted)

conc.modeled3 = ode(y = state, times = times, func = salttank, parms = params,
method = "rk4")

plot(t, conc.observed, pch = 21, bg = "gray", ylim = c(20, 180), xlab = "Time",
ylab = "Salt (kg)", axes = F)
box(bty = "l")
axis(1, tck = -0.02)
axis(2, tck = -0.02, las = 2)
lines(conc.modeled[, "time"], conc.modeled[, "X"], col = "red")
lines(conc.modeled3[, "time"], conc.modeled3[, "X"], col = "blue")

legend("bottomright", legend = c("Observed data", "True solution (a = 0.75)",
paste("Simulated with fitted parameters (a = ", sprintf("%2.2f", params),
")", sep = "")), lty = c(NA, 1, 1), pch = c(21, NA, NA), col = c("black",
"red", "blue"), pt.bg = "gray", bty = "n")
`

1. Super,
Je pensais que tu aurais fait cette modÃ©lisation pour la croissance d'une population de levure de biÃ¨re! Peut Ãªtre pour une prochaine.

Sinon quand est-ce que tu t'attaque Ã  PARAFAC?

Bravo pour le travail de ton blog!
Thomas

2. Merci Thomas!

L'idÃ©e de la biÃ¨re serait trÃ¨s intÃ©ressante. D'ailleurs, cela ferait un excellent sujet pour un cours Ã  la maÃ®trise avec les Ã©tudiants. Suivre la cinÃ©tique des levures dans une expÃ©rience de fermentation.

Bon courage,
Phil

3. Many thanks for making the truthful effort to explain this. I feel very strong about it and would like to read more. If you
Qassim University
.,

1. Thank you for your comment. I have some posts coming soon with more detailed examples.

4. Thank you for the example. I am also more familiar with matlab, maple and mathematica. Do you know how to fit also the initial condition ? i.e state = c(X = 20), intead of 20 a new parameter X0 ? Thanks.