Sunday, July 8, 2012

Fitting a dynamic model, and determining the number of parameters that can be fitted.

Let's suppose that we have the same dynamic model we presented before - that is, the Lorentz system of differential equations. Remember?

library(deSolve) # require this library
solveLorenz <- function(pars, times=seq(0,100,by=0.1)) {
derivs <- function(t, state, pars) { # returns rate of change
with(as.list(c(state, pars)), {
dX <- a*X + Y*Z
dY <- b * (Y-Z)
dZ <- -X*Y + c*Y - Z
return(list(c(dX, dY, dZ)))
}
)
}
state <- c(X = 1, Y = 1, Z = 1)
## ode solves the model by integration...
return(ode(y = state, times = times, func = derivs, parms = pars))
}
view raw gistfile1.r hosted with ❤ by GitHub
In order to perform a fitting we need to define an objective function of sort: this will then be minimised.
library(FME)
Objective <- function(x, parset = names(x)) {
guess_pars[parset] <- x
tout <- seq(0, 50, by = 0.5)
## output times
bullet <- solveLorenz(guess_pars, tout)
## Model cost
return(modCost(obs = target, model = bullet))
}
view raw gistfile1.r hosted with ❤ by GitHub
Now, the two most important pieces are put together. What we still need is some data to fit, and a starting point for the fitting. We get those in the next few lines, followed by the actual command to fit the guess to our actual target (note: target is explicitly called within the Objective function itself, in the modCost() line).
target_pars <- c(a = -9/3, b = -5, c = 30); target_pars
target <- solveLorenz(target_pars)
guess_pars<-c(a = -6/3, b = -8, c = 20); guess_pars
guess <- solveLorenz(guess_pars)
# print(system.time(Fit <- modFit(p = guess_pars, f = Objective))) # this works
# print(system.time(Fit <- modFit(p = guess_pars, f = Objective, method="SANN", control=list(maxit=100), lower=c(-5,-10,0)))) # this works too
# print(system.time(Fit <- modFit(p = guess_pars, f = Objective, method="SANN", control=list(maxit=10000)))) # this works too
print(system.time(Fit <- modFit(p = guess_pars, f = Objective, method="Nelder-Mead", control=list(maxit=1000)))) # this works too
# print(system.time(Fit <- modFit(p = guess_pars, f = Objective, method="Nelder-Mead", control=list(maxit=10000)))) # this works too
Fit$par
bullet_pars <- Fit$par
bullet <- solveLorenz(bullet_pars)
Now we're in business! After running the code above, the system suggests a value of coefficients fairly close to the target's owns. -3.6, -5.2, 27.7, compared to -3, -5, 30. let's have a look at the solutions: in blue, our target - in red, our initial guess, in green, our final fit (bullet).
pXY<-ggplot(as.data.frame(bullet)) +geom_path(aes(X, Y, alpha=Z), col='green') + opts(legend.position = "none")
pXY<-pXY +geom_path(data=as.data.frame(target), aes(X, Y, alpha=Z), col='blue')
pXY<-pXY +geom_path(data=as.data.frame(guess), aes(X, Y, alpha=Z), col='red')
print(pXY)
Ok, the fit isn't that good if you look at it. But it's better than nothing... If you want to try and get closer just increase the number of iterations allowed to the minimiser. But it isn't really necessary - we're going to move onto something more interesting. How do we know that the parameters we're fitting are sensible? Yes, they certainly look so in this case, but in real life we have no clue of what the real numbers may be, or we wouldn't be fitting to them in first place. Perhaps the trajectory we're fitting against is just too short, and doesn't have enough information to distinguish between an infinite number of radically different solutions, which fit well locally but diverge wildly later on... Pretty much assured to be the case with Lorentz's system, I'd say... So, we run a sensitivity analysis on the system - the authors of the FME package have made this very easy for us:
sf<-sensFun(func = Objective, parms = bullet_pars, varscale = 1)
Coll <- collin(sF)
the result of this is then fed to a collinearity analysys function, which - at least that's how I understand it, in 'non-statitiscian' speaks - checks that simultaneous changes in (a subset of) your parameters do not affect the model in a too similar (collinear) way. So, here's your result. If you print the 'Coll' variable, you'll see a table-like frame where a collinearity value is attached to each combination of parameters... Like this:
> Coll
a b c N collinearity
1 1 1 0 2 1.1
2 1 0 1 2 2.8
3 0 1 1 2 1.0
4 1 1 1 3 2.9
>
view raw 4_Coll_out hosted with ❤ by GitHub
Now, according to the authors of the package, values below 20 say that it's OK to trust a fitting of that particular combination. value above 20, are a no-no... This table is easy to look at but I thought that transforming it into a graph may come in handy, particularly for cases several parameters and combinations have to be explored. First off we 'melt' our Collinearity table using the reshape package:
require(reshape)
# select only the columns relative to parameter on/off, and the collinearity value
Coll_m<-melt(t(Coll[,c(1:length(guess_pars),ncol(Coll))]));
# split them in two different frames for easier handling
Coll_f<-Coll_m[Coll_m$X1=='collinearity',];
Coll_m<-Coll_m[Coll_m$X1!='collinearity',];
# add a categorical label
Coll_f$label[Coll_f$value<20]<- sprintf("%3.2g", Coll_f$value[Coll_f$value<20]); Coll_f$label[Coll_f$value>=20] <- ">= 20"; Coll_f$label[Coll_f$value>=100] <- "> 100";
Coll_f$flag<-'ID OK'; Coll_f$flag[Coll_f$value>=20] <- "ID KO";
Coll_f$flag <- factor(Coll_f$flag, levels=c('ID OK', 'ID KO'))
view raw 5_MeltColl.r hosted with ❤ by GitHub
As you may see, we also create some categorical (i.e. factor() variables which make or easy coloring of the table. Then, we can finally plot:
grid.newpage()
# create a tile matrix with 0 and 1 depending on which parameter is being fitted.
cl1<-ggplot(Coll_m) + geom_tile(aes(X1, X2, fill=value), , col='white') + scale_fill_gradient(low='white', high='black', limits=c(0, 1))
# fill the same matrix with tiles of two different colors
cl1<-cl1 + geom_text(aes(X1, X2, label=value, col=value)) + scale_color_gradient(low='black', high='white', limits=c(0, 1))
cl1<-cl1 + opts(title="Parameters' N-ways Identifiability", axis.title.x = theme_blank(), axis.title.y = theme_blank(), axis.text.y = theme_text(colour = NA), axis.ticks = theme_segment(colour = NA), axis.text.x = theme_text(size = 12), legend.position='none')
# now, on the right side, fill colors depending on the collinearity values
cl2<-ggplot(Coll_f) + geom_tile(aes(X1, X2, fill=flag), col='white') + scale_fill_manual(values=c('green','red'))
# finally print clipped collinearity values
cl2<-cl2 + geom_text(data=Coll_f, aes(X1, X2, label=label), col='black') + theme_bw() + guides(fill=guide_legend(title=NULL)) +
# and all options to make the graph look right
opts(title="", axis.title.x = theme_blank(), axis.title.y = theme_blank(), axis.text.y = theme_text(colour = NA), axis.ticks = theme_segment(colour = NA), axis.text.x = theme_text(size = 12), legend.position='none')
pushViewport(viewport(layout = grid.layout(6, 7)));
print(cl1, vp=vplayout(1:6,1:5))
print(cl2, vp=vplayout(1:6,6:7))
view raw 6_Coll_Plot.r hosted with ❤ by GitHub
And here's the final result:
This graph suggests, like the table from which is derived, that all parameters can be fit together (with these specific data), and a reliable fit can be found... If only my real-life example were so easy... Homeworks: Try adding some random noise to the XYZ values of the target, and see if the collinearity analysis changes. I'll try to update this post with a tentative walkthrough in the next few days. For now, enjoy!

No comments:

Post a Comment