Submitted by dylan on Tue, 2010-02-23 21:42.

Stumbled upon an excellent example of how to perform numerical integration in R. Below is an example of piece-wise linear and spline fits to FTIR data, and the resulting computed area under the curve. With a high density of points, it seems like the linear approximation is most efficient and sufficiently accurate. With very large sequences, it may be necessary to adjust the value passed to the `subdivisions` argument of `integrate()`. Strangely, larger values seem to solve problems encountered with large datasets...

**FTIR Spectra Integration**

**Implementation**

# numerical integration in R

# example based on: http://tolstoy.newcastle.edu.au/R/help/04/10/6138.html

# sample data: FTIR spectra

x <- read.csv(url('http://casoilresource.lawr.ucdavis.edu/drupal/files/fresh_li_material.CSV'), header=FALSE)[100:400,]

names(x) <- c('wavenumber','intensity')

# fit a piece-wise linear function

fx.linear <- approxfun(x$wavenumber, x$intensity)

# integrate this function, over the original limits of x

Fx.linear <- integrate(fx.linear, min(x$wavenumber), max(x$wavenumber))

# fit a smooth spline, and return a function describing it

fx.spline <- splinefun(x$wavenumber, x$intensity)

# integrate this function, over the original limits of x

Fx.spline <- integrate(fx.spline, min(x$wavenumber), max(x$wavenumber))

# visual check, linear and spline fits shifted up for clarity

par(mar=c(0,0,0,0))

plot(x, type = "p", las=1, axes=FALSE, cex=0.5, ylim=c(0,0.12))

lines(x$wavenumber, fx.linear(x$wavenumber) + 0.01, col=2)

lines(x$wavenumber, fx.spline(x$wavenumber) + 0.02, col=3)

grid(nx=10, col=grey(0.5))

legend(x=615, y=0.11, legend=c('original','linear','spline'), col=1:3, pch=c(1,NA,NA), lty=c(NA, 1, 1), bg='white')

# results are pretty close

data.frame(method=c('linear', 'spline'), area=c(Fx.linear$value, Fx.spline$value), error=c(Fx.linear$abs.error,Fx.spline$abs.error))

method area error

1 linear 27.71536 0.0005727738

2 spline 27.71527 0.0030796529

`splinefun()` can also compute derivatives

par(mar=c(0,0,0,0), mfcol=c(2,1))

plot(x, type = "l", lwd=2, axes=FALSE)

grid(nx=10, col=grey(0.5))

plot(x$wavenumber, fx.spline(x$wavenumber, deriv=1), type='l', axes=FALSE)

lines(x$wavenumber, fx.spline(x$wavenumber, deriv=2), col='red')

grid(nx=10, col=grey(0.5))

abline(h=0, lty=2)

legend('topright', legend=c('1st derivative','2nd derivative'), lty=1, col=1:2, bg='white')

**Numerical Estimation of Derivatives**

## integration

Dylan, thanks for showing the use of approxfun and finding derivatives. But when I tried it with my own data, it didn't work quite right when trying to plot the original, spline, and piecewise fits on top of one another. The original data had 6 humps in series, but the smoothers got stretched out, and only captured four of the humps before running outside of the plot. Somehow the horizontal scale got stretched for the smoothers. thanks, Paul

## RE:

Hi Paul,

Without seeing your data and your specific incantation of R code, I can't offer much help. Posting an example on the R-help mailing list with reproducible code is usually the best way to get answers fast. I would also check for missing values in your x-variable, and working with a subset of your data first. Good luck!