########## R script: sitka-spruce.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: 05 APR 2004 by M.P. Wand

use.pdf <- F

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

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

idnum <- sitka$idnum
order <- sitka$order
days <- sitka$days
log.size <- sitka$log.size
ozone <- sitka$ozone

ozone.char <- as.character(ozone)
ozone.char[ozone.char=="0"] <- "control"
ozone.char[ozone.char=="1"] <- "ozone treated"

# Plot raw sitka data

if (use.pdf)
   pdf("sitka-raw.pdf",width=8)

pobj <-
xyplot(log.size ~ days |ozone.char, data =sitka,groups = idnum,
        strip = function(...) strip.default(..., style = 1),
        subscripts = T, panel = function(x, y, subscripts,groups) 
        {
           panel.grid()
           panel.xyplot(x, y,pch=16,col="black")
           panel.superpose(x, y, subscripts, groups, type = "l",lwd=2)
        }
        ,xlab = "number of days", ylab ="log(size)",
         bg="white")
         
print(pobj)

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

knots <- default.knots(days)

# Construct the response vector and X matrix

y <- log.size
X <- cbind(rep(1,length(y)),days,ozone)

# Create the Z matrix structure.

Z.spline <- outer(days,knots,"-")
Z.spline <- Z.spline*(Z.spline>0)         

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

data.fr <- groupedData(y~-1+X|group,data=data.frame(y,X))

Z.block <- list(list(group=pdIdent(~Z.spline-1)),list(idnum=pdIdent(~1)))

Z.block <- unlist(Z.block,recursive=FALSE)

fit <- lme(y~-1+X,data=data.fr,random=Z.block)

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))

print(summary(fit)$tTable)
print(intervals(fit))

########## End of sitka-spruce ##########