########## R script: ozone-calif.Rs ##########

# For illustrating penalized spline fitting
# of an additive model. 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 an extract data

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

pairs(ozone.calif)

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

y <- ozone.calif$ozone

x1 <- ozone.calif$inversion.base.height
x2 <- ozone.calif$daggot.pressure.gradient
x3 <- ozone.calif$inversion.base.temp

# Set up matrices for additive model

knots.1 <- default.knots(x1,15)
knots.2 <- default.knots(x2,15)
knots.3 <- default.knots(x3,15)

K.1 <- length(knots.1)
K.2 <- length(knots.2)
K.3 <- length(knots.3)


re.block.inds <- list(1:K.1,(K.1+1):(K.1+K.2),
                     (K.1+K.2+1):(K.1+K.2+K.3))         

X <- cbind(rep(1,length(y)),x1,x2,x3)

Z.1 <- outer(x1,knots.1,"-")
Z.1 <- Z.1*(Z.1>0)

Z.2 <- outer(x2,knots.2,"-")
Z.2 <- Z.2*(Z.2>0)

Z.3 <- outer(x3,knots.3,"-")
Z.3 <- Z.3*(Z.3>0)

Z <- cbind(Z.1,Z.2,Z.3)


# Fit model using lme() and extract coefficient estimates

Z.block  <-  list()

for (i in 1:length(re.block.inds))
   Z.block[[i]] <- as.formula(paste("~Z[,c(",paste(
             re.block.inds[[i]],collapse=","),")]-1"))

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

data.fr <- groupedData(y~-1+X|group,data=data.frame(y,X,Z))
                  
fit <- lme(y~-1+X,random=list(group=pdBlocked(Z.block,
           pdClass=rep("pdIdent",length(Z.block)))),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))

# Draw plots

if (create.pdf)
   pdf("calif-ozone.pdf",width=8)

par(mfrow=c(1,3))

n.g <- 101
ylim.val <- c(0,30)

# First function

x1.g <- seq(min(x1),max(x1),length=n.g)
x2.g <- rep(mean(x2),length=n.g)
x3.g <- rep(mean(x3),length=n.g)

X1.g <- x1.g 
Z1.g <- outer(x1.g,knots.1,"-")
Z1.g <- Z1.g*(Z1.g>0)

X2.g <- x2.g
Z2.g <- outer(x2.g,knots.2,"-")
Z2.g <- Z2.g*(Z2.g>0)

X3.g <- x3.g
Z3.g <- outer(x3.g,knots.3,"-")
Z3.g <- Z3.g*(Z3.g>0)

X.g <- cbind(rep(1,n.g),X1.g,X2.g,X3.g)
Z.g <- cbind(Z1.g,Z2.g,Z3.g)

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

plot(x1.g,f1.g,col="red",lwd=4,type="l",bty="l",
     xlab="inversion base height",ylab="ozone",ylim=ylim.val)
rug(x1)

# Second function

x1.g <- rep(mean(x1),length=n.g)
x2.g <- seq(min(x2),max(x2),length=n.g)
x3.g <- rep(mean(x3),length=n.g)

X1.g <- x1.g
Z1.g <- outer(x1.g,knots.1,"-")
Z1.g <- Z1.g*(Z1.g>0)

X2.g <- x2.g
Z2.g <- outer(x2.g,knots.2,"-")
Z2.g <- Z2.g*(Z2.g>0)

X3.g <- x3.g
Z3.g <- outer(x3.g,knots.3,"-")
Z3.g <- Z3.g*(Z3.g>0)

X.g <- cbind(rep(1,n.g),X1.g,X2.g,X3.g)
Z.g <- cbind(Z1.g,Z2.g,Z3.g)

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

plot(x2.g,f2.g,col="blue",lwd=4,type="l",bty="l",
     xlab="Daggot pressure gradient",ylab="ozone",ylim=ylim.val)
rug(x2)

# Third function

x1.g <- rep(mean(x1),length=n.g)
x2.g <- rep(mean(x2),length=n.g)
x3.g <- seq(min(x3),max(x3),length=n.g)

X1.g <- x1.g
Z1.g <- outer(x1.g,knots.1,"-")
Z1.g <- Z1.g*(Z1.g>0)

X2.g <- x2.g
Z2.g <- outer(x2.g,knots.2,"-")
Z2.g <- Z2.g*(Z2.g>0)

X3.g <- x3.g
Z3.g <- outer(x3.g,knots.3,"-")
Z3.g <- Z3.g*(Z3.g>0)

X.g <- cbind(rep(1,n.g),X1.g,X2.g,X3.g)
Z.g <- cbind(Z1.g,Z2.g,Z3.g)

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

plot(x3.g,f3.g,col="orange",lwd=4,type="l",bty="l",
     xlab="inversion base temperature",ylab="ozone",ylim=ylim.val)
 
rug(x3)

########## End of ozone-calif ##########