########## R script: fossil.Rs ##########

# For illustrating penalized spline fitting
# of ordinary scatterplot data. Note that the
#  R package "nlme" is required as well as the 
# file default.knots.r containing the function
# default.knots().

# Last changed: 02 APR 2004 by M.P. Wand

create.pdf <- F  # flag for creation of pdf file

library(nlme)
source("default.knots.r")

# Read in and extract data

fossil <- read.table("fossil.dat",header=T)

x <- fossil[,1]
y <- 10000*fossil[,2]

knots <- default.knots(x)

X <- model.matrix(~x)
Z <- outer(x,knots,"-")
Z <- Z*(Z>0)
group <- rep(1,length(x))

data.fr <- groupedData(y~x|group,data=data.frame(x,y))
        	          
fit <- lme(y~-1+X,random=pdIdent(~-1+Z),data=data.fr)

beta.hat <- fit$coef$fixed
u.hat <- unlist(fit$coef$random)

sig.eps.hat <- fit$sigma
sig.u.hat <- sig.eps.hat*exp(unlist(fit$modelStruct))

n.g <- 101
x.g <- seq(min(x),max(x),length=n.g)
X.g <- cbind(rep(1,n.g),x.g)

Z.g <- outer(x.g,knots,"-")
Z.g <- Z.g*(Z.g>0)

f.g <- X.g%*%beta.hat + Z.g%*%u.hat

if (create.pdf)
   pdf("fossil-lin.pdf",width=8)

plot(x,y,bty="l",xlab="age",ylab="10000(strontium ratio)",
     type="n",main="linear spline fit")

points(x,y)
lines(x.g,f.g,col="blue",lwd=2)

cat("Hit return to continue\n")
ans <- readline()

# Now do radial cubic fit

svd.Omega <- svd(abs(outer(knots,knots,"-"))^3)
matrix.sqrt.Omega <- t(svd.Omega$v %*% (t(svd.Omega$u)
                                  *sqrt(svd.Omega$d)))
Z <- t(solve(matrix.sqrt.Omega,t(abs(outer(x,knots,"-")^3))))

group <- rep(1,length(x))

data.fr <- groupedData(y~x|group,data=data.frame(x,y))
        	          
fit <- lme(y~-1+X,random=pdIdent(~-1+Z),data=data.fr)

beta.hat <- fit$coef$fixed
u.hat <- unlist(fit$coef$random)

sig.eps.hat <- fit$sigma
sig.u.hat <- sig.eps.hat*exp(unlist(fit$modelStruct))

n.g <- 101
x.g <- seq(min(x),max(x),length=n.g)
X.g <- cbind(rep(1,n.g),x.g)

Z.g <- t(solve(matrix.sqrt.Omega,t(abs(outer(x.g,knots,"-")^3))))

f.g <- X.g%*%beta.hat + Z.g%*%u.hat

if (create.pdf)
   pdf("fossil-cub.pdf",width=8)

plot(x,y,bty="l",xlab="age",ylab="10000(strontium ratio)",
     type="n",main="cubic spline fit")

points(x,y)
lines(x.g,f.g,col="blue",lwd=2)

########## End of fossil ##########