[R] questions regarding spline functions

[Dylan Beaudette]

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&_user=4421&_handle=V-WA-A-W-AU-MsSAYWW-UUA-U-AACZEDZWBC-AACVCCDUBC-BZAYUEWB-AU-U&_fmt=summary&_coverDate=08%2F31%2F1999&_rdoc=3&_orig=browse&_srch=%23toc%235807%231999%23999089998%23108393!&_cdi=5807&view=c&_acct=C000059598&_version=1&_urlVersion=0&_userid=4421&md5=488f1e114d8d64265ff65506e9587e71

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
530.754.7341
-------------- next part --------------
A non-text attachment was scrubbed...
Name: spline1.png
Type: image/png
Size: 3953 bytes
Desc: not available
Url : https://stat.ethz.ch/pipermail/r-help/attachments/20060731/5a5b9df6/attachment.png