|
Fitting an Age model One of the basic tasks in analyzing a
sedimentological record is the construction of an age model, or an age vs.
depth relationship that can be used to assign ages to samples between
horizons that have been dated by, for example, radiocarbon. The
following illustration uses Sarah Millspaugh's data from her Geology
paper (Millspaugh, et al., 2000, Geology 28:211-214).
Plot the data
plot(calAge ~ Depth)
First, consider a straight-line fit
# Cygnet
Lake age models
attach(cygage)
# linear (straight line)
cygage.lm1 <- lm(calAge ~ Depth)
# plot the fitted line
plot(calAge ~ Depth)
abline(cygage.lm1)
segments(Depth, fitted(cygage.lm1), Depth, calAge)
# examine the model
summary(cygage.lm1)
oldpar <- par(mfrow = c(2, 2))
plot(cygage.lm1)
par(oldpar)
Second and third-degree polynomial age models. Note the
incorporation of observations weights.
# quadratic (2nd-degree
polynomial (parabola))
cygage.lm2 <- lm(calAge ~ Depth + I(Depth^2), weight=Weight)
# plot the fitted curve
plot(calAge ~ Depth)
lines(spline(fitted(cygage.lm2) ~ Depth), col="green")
segments(Depth, fitted(cygage.lm2), Depth, calAge, col="green")
# examine the model
summary(cygage.lm2)
oldpar <- par(mfrow = c(2, 2))
plot(cygage.lm2)
par(oldpar)
# cubic (3rd-degree polynomial)
cygage.lm3 <- lm(calAge ~ Depth + I(Depth^2) + I(Depth^3),
weight=Weight)
# plot the fitted curve
plot(calAge ~ Depth)
lines(spline(fitted(cygage.lm3) ~ Depth), col="red")
segments(Depth, fitted(cygage.lm3), Depth, calAge, col="red")
# examine the model
summary(cygage.lm3)
oldpar <- par(mfrow = c(2, 2))
plot(cygage.lm3)
par(oldpar)
Interpolate ages for new samples, illustrating the use of the
predict() function
# get ages for new depths
cygage.pred3 <- predict(cygage.lm3, newdata=cygdepths)
plot(calAge ~ Depth)
lines(cygage.pred3, col="blue")
Now try a couple of other curve-fitting methods
# lowess
cygage.lo1 <- loess(calAge ~ Depth, span=0.8)
# plot the fitted curve
plot(calAge ~ Depth)
hat <- predict(cygage.lo1)
lines(Depth[order(Depth)], hat[order(Depth)], col="red")
segments(Depth, fitted(cygage.lo1), Depth, calAge, col="red")
# examine the model
summary(cygage.lo1)
# smoothing spline
cygage.spline1 <- smooth.spline(calAge ~ Depth )
# plot the fitted curve
plot(calAge ~ Depth)
hat <- predict(cygage.spline1)
lines(spline(hat$x, hat$y))
segments(Depth, hat$y, Depth, calAge, col="orange")
# examine the model
cygage.spline1
|