tests/example2nd.R

In TIMP: Fitting Separable Nonlinear Models in Spectroscopy and Microscopy

## Example in modeling second order kinetics, by
## David Nicolaides.

## On simulated data.

##############################
## load TIMP
##############################

library("TIMP")

##############################
## SIMULATE DATA
##############################

## set up the Example problem, a la in-situ UV-Vis spectroscopy of a simple
## reaction.
## A + 2B -> C + D, 2C -> E

cstart <- c(A = 1.0, B = 0.8, C = 0.0, D = 0.0, E = 0.0)
times <- c(seq(0,2, length=21), seq(3,10, length=8))
k <- c(kA = 0.5, k2C = 1)

## stoichiometry matrices as per
## Puxty, G., Maeder, M., and Hungerbuhler, K. (2006) Tutorial on the fitting
## of kinetics models to mulivariate spectroscopic measurements with
## non-linear least-squares regression, Chemometrics and Intelligent
## Laboratory Systems 81, 149-164.

rsmatrix <- c(1,2,0,0,0,0,0,2,0,0)
smatrix <- c(-1,-2,1,1,0,0,0,-2,0,1)
concentrations <- calcD(k, times, cstart, rsmatrix, smatrix)

wavelengths <- seq(500, 700, by=2)
spectra <- matrix(nrow = length(wavelengths), ncol = length(cstart))
location <- c(550, 575, 625, 650, 675)
delta <- c(10, 10, 10, 10, 10)
spectra[, 1] <- exp( - log(2) * (2 * (wavelengths - location[1])/delta[1])^2)
spectra[, 2] <- exp( - log(2) * (2 * (wavelengths - location[2])/delta[2])^2)
spectra[, 3] <- exp( - log(2) * (2 * (wavelengths - location[3])/delta[3])^2)
spectra[, 4] <- exp( - log(2) * (2 * (wavelengths - location[4])/delta[4])^2)
spectra[, 5] <- exp( - log(2) * (2 * (wavelengths - location[5])/delta[5])^2)

sigma <- .001
Psi_q <- concentrations %*% t(spectra) + sigma *
  rnorm(dim(concentrations)[1] * dim(spectra)[1])

## store the simulated data in an object of class "dat"
kinetic_data <- dat(psi.df=Psi_q , x = times, nt = length(times),
 x2 = wavelengths, nl = length(wavelengths))

##############################
## DEFINE MODEL 
##############################

## starting values
kstart <- c(kA = 1, k2C = 0.5)

## model definition for 2nd order kinetics
kinetic_model <- initModel(mod_type = "kin", seqmod = FALSE, kinpar = kstart,
                           numericalintegration = TRUE,
                           initialvals = cstart, 
                           reactantstoichiometrymatrix = rsmatrix, 
                           stoichiometrymatrix = smatrix )

##############################
## FIT INITIAL MODEL 
## adding constraints to non-negativity of the
## spectra via the opt option nnls=TRUE
##############################

kinetic_fit <- fitModel(data=list(kinetic_data), modspec = list(kinetic_model),
                        opt = kinopt(nnls = TRUE, iter=80,
                          selectedtraces = seq(1,kinetic_data@nl,by=2)))

## look at estimated parameters

parEst(kinetic_fit)

## make a png of various results

## concentrations 

conRes <- getX(kinetic_fit)

matplot(times, conRes, type="b", col=1,pch=21, bg=1:5, xlab="time (sec)",
        ylab="concentrations", main="Concentrations (2nd order kinetics)")
                        
                        
## spectra 

specRes <- getCLP(kinetic_fit)

matplot(wavelengths, specRes, type="b", col=1,pch=21, bg=1:5,
        xlab="wavelength (nm)",
        ylab="amplitude", main="Spectra")


## see help(getResults) for how to get more results information from
## kinetic_fit