[R] questions regarding spline functions

Dylan Beaudette dylan.beaudette at gmail.com
Tue Aug 1 20:06:16 CEST 2006

Posting reply to the list as well:

On Tuesday 01 August 2006 07:02, Myers, Brent wrote:
> Dylan,
> We are working on the same problem. I'm shocked that I hadn't seen the
> paper you posted and I am most grateful for that. I try to pay attention
> to McBratney's work too!

Hi Brent, glad to hear from a fellow soil scientist.
I originally found the McBratney paper in the book: Upscaling and Downscaling 
Methods for Environmental Research.

> The originator of the area-preserving spline is DeBoor, in A Practical
> Guide to Splines, Chapter 8. The DeBoor algorithm is available in the
> Matlab Spline Toolkit, which he wrote. Too pricey for me however. I have
> been unsuccessful implementing it thus far in S-Plus, the reference is
> not very clear in the construction of the coefficient matrix. I can send
> you my notes on that. I've paused on that for now, more below. Notice
> that these are quadratic splines (due to the number of constraints) and
> that they are not so 'flexible' in representing soils data.

Interesting. I will lookup the original author, and see if I can make any 
sense of it. This is not critical, merely a bit of initial exploration on 
horizon depth function characterizations for large data sets. 

> There are several other functions available from a CERN library
> available also in their Physics Analyst Workstation software (GPL). I
> haven't tried them. Realize also that some of the various density
> estimation algorithms in S/R are spline fits of this nature.

I will have to read up on those...

> I also considered the possibility of boundary data to define slope at
> some of the knots, but only about 1/3 of my data have good horizon
> descriptions. Linear stairstep with sloped lines across boundaries are
> trivial, but are a marked improvement over mid-horizon data.

I too have seen similar results. However, I cannot quite get the bs() function 
to use the "knots" that I have in mind: note that this may be an operator 
error... Does it make sense to use the horizon boundaries as knots?

> I am sure you are thinking about this...one big problem with area
> preserving splines is that they will not represent the true distribution
> of the soil property within the horizon if there are peaks. You
> undoubtedly realize that since the vector of pedogenesis is anisotropic,
> and that splines cannot represent this anisotropy. This is a real
> problem for me in that it is critical to identify the location of minima
> and maxima within the soil profile, more important than the property
> value in my case. Their location is grossly distorted by a spline fit to
> a mid-horizon data point. The mid-point data is better. I am taking a
> separate approach to solve that problem, and then fit the area
> preserving spline with knots defined at these locations. We need a
> general algorithm to handle this data!!!! It would be very useful for
> the soil science community!

Indeed. If you wouldn't mind, I would be interested in hearing how this 
progresses. While we are not in a position to develop such a technique, 
outside of something rather simple, we sure would be able to put it to use.

A quick summary in code and with an attached image of the various methods I 
have tried thus far, with one simple, and one complicated clay depth profile:

#example of two pedons on one plot:
z.1 <- c(0,2,18,24,68,160,170,192,200)
z.2 <- c(0,3,14,26,70,108,145,170,226,240)

x.1 <- mid(z.1)
x.2 <- mid(z.2)

#clay pct
y.1 <- c(0,1,2,2,4,7,6,1)
y.2 <- c(0,5,3,3,3,2,27,3,5)

fm.1 <- lm(y.1 ~ bs(x.1, df=5) )
fm.2 <- lm(y.2 ~ bs(x.2, df=5) )

isfm.1 <- interpSpline(y.1 ~ x.1)
isfm.2 <- interpSpline(y.2 ~ x.2)

new_x.1 <-  seq(range(x.1)[1], range(x.1)[2], len = 200)
new_x.2 <-  seq(range(x.2)[1], range(x.2)[2], len = 200)

new_y.1 <- predict(fm.1, data.frame(x.1=new_x.1) )
new_y.2 <- predict(fm.2, data.frame(x.2=new_x.2) )

#interpSpline() method
#note different predict() syntax
ispline_y.1 <- predict(isfm.1, new_x.1)
ispline_y.2 <- predict(isfm.2, new_x.2)


#simple pedon example
plot(y.1 ~ x.1, xlab="Depth", ylab="Percent Clay", type="b", pch=16, lwd=2, 
main="Simple Pedon")
#spline examples
lines(spline(y.1 ~ x.1), col="blue", lty=2)
lines(new_x.1, new_y.1, col='red', lty=2)
lines(ispline_y.1, col='green', lty=2)

legend(8,6, legend=c('soil data', 'splines()', 'interpSpline()', 'bs()'), 
col=c('black', 'blue', 'red', 'green'), lty=c(1,1,1,1), lwd=c(2,1,1,1), 

#complex pedon:
plot(y.2 ~ x.2, xlab="Depth", ylab="Percent Clay", type="b", pch=16, lwd=2, 
main="Complex Pedon")
#spline examples
lines(spline(y.2 ~ x.2), col="blue", lty=2)
lines(new_x.2, new_y.2, col='red', lty=2)
lines(ispline_y.2, col='green', lty=2)

legend(8.4,24, legend=c('soil data', 'splines()', 'interpSpline()', 'bs()'), 
col=c('black', 'blue', 'red', 'green'), lty=c(1,1,1,1), lwd=c(2,1,1,1), 



> Brent
> D. Brenton Myers
> Graduate Fellow
> University of Missouri
> Soil Environmental and Atmospheric Sciences
> 269 Agriculture Engineering Building
> Columbia, MO 65211
> Work: (573) 882-1146
> Home: (573) 446-6856
> Cell: (573) 881-6855
> email: myersdb at missouri.edu
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Dylan Beaudette
> Sent: Monday, July 31, 2006 6:43 PM
> Subject: [R] questions regarding spline functions
> Greetings,
> A couple general questions regarding the use of splines to interpolate
> depth
> profile data.
> Here is an example of a set of depths, with associated attributes for a
> given
> soil profile, along  with a function for calculating midpoints from a
> set of
> soil horizon boundaries:
> #calculate midpoints:
> mid <- function(x) {
> for( i in 1:length(x)) {
>  if( i > 1) {
>    a[i] = (x[i] - x[i-1]) / 2 + x[i-1]
>   }
>  }
> #reurn the results
> a[which(!is.na(a))]
> }
> #horizon depth bounds
> z <- c(0,2,18,24,68,160,170,192,200)
> #horizon midpoints, associated with horizon attribute
> x <- mid(z)
> #clay pct
> y <- c(0,1,2,2,4,7,6,1)
> #plot them
> plot(y ~ x, xlab="Depth", ylab="Percent Clay", type="s")
> points(y ~ x, cex=0.5, pch=16)
> These point pairs usually represent a trend with depth, which I would
> like to
> model with splines - or some similar approach, as they have been found
> to
> work better than other methods such as a fitted polynomial.
> Using the B Spline function from the 'splines' package, it is possible
> to fit
> a model of some property with depth based on the bs() function:
> #natual, B-Splines
> library(splines)
> #fit a b-spline model:
> fm <- lm(y ~ bs(x, df=5) )
> I am able to predict a soil property with depth, at unsampled locations
> with
> this model with:
> new_x <-  seq(range(x)[1], range(x)[2], len = 200)
> #predict attribute at unsampled depths:
> new_y <- predict(fm, data.frame(x=new_x) )
> #plot the predicted attribute at the unsampled depths
> lines(new_x, new_y, col='red')
> This tends to work fairly well (see attached), but I am wondering if I
> can use
> the 'knots' parameter in the bs() function for incorporation of horizon
> boundary information into the creation of the spline. Moreover, it would
> be
> nice if the spline-based model 'fm' would predict a set of values with
> similar mean and range as the original data points: i.e
> #summary of clay values from original data:
> summary(y)
>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
>   0.000   1.000   2.000   2.875   4.500   7.00
> #see above
> summary(new_y)
>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
> -0.05786  2.09500  3.13200  3.62800  5.17100  7.08700
> This is based on an article I read :
> http://www.sciencedirect.com/science?_ob=ArticleURL&_udi=B6V67-3WWRDYY-3
> UEWB-AU-U&_fmt=summary&_coverDate=08%2F31%2F1999&_rdoc=3&_orig=browse&_s
> rch=%23toc%235807%231999%23999089998%23108393!&_cdi=5807&view=c&_acct=C0
> 00059598&_version=1&_urlVersion=0&_userid=4421&md5=488f1e114d8d64265ff65
> 506e9587e71
> where the author talks about a so-called 'equal-area quadratic smoothing
> spline' approach to describing a soil property depth function.
> Unfortunately
> the author did not provide sample code....
> Any thoughts / input would be greatly appreciated!
> Cheers,

Dylan Beaudette
Soils and Biogeochemistry Graduate Group
University of California at Davis
-------------- next part --------------
A non-text attachment was scrubbed...
Name: spline2.png
Type: image/png
Size: 33301 bytes
Desc: not available
Url : https://stat.ethz.ch/pipermail/r-help/attachments/20060801/4f119153/attachment.png 

More information about the R-help mailing list