########## R script: sitkaNew ##########

# For illustrating improved use of lme() for
# fitting additive mixed models in R. Note that
# SemiPar 1.0 does NOT use this improvement,
# but it will be incorporated into future
# releases of SemiPar. 

# Last changed: 05 DEC 2005

library(SemiPar)
data(sitka) ; attach(sitka)

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

num.knots.val <- 7
knots <- default.knots(days,num.knots=num.knots.val)

# Create the spline component of the Z matrix.

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

# The component of the Z matrix corresponding to the random
# intercept does not need to be specified and can be handled
# through the identification numbers stored in "id.num", as
# follows.

n <- length(y)
group.vec <- factor(rep(1,n))
Z.block <- list(list(group.vec=pdIdent(~Z.spline-1)),
                list(id.num=pdIdent(~1)))

Z.block <- unlist(Z.block,recursive=FALSE)
data.fr <- groupedData(y~-1+X|group.vec,
                       data=data.frame(y,X,Z.spline,id.num))

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

# Extract the estimated parameters.

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

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

# Print summary table of linear predictors.

print(round(summary(fit)$tTable[1:2,],4))

# Plot the estimated effect of "days".

ng <- 101
days.g <- seq(min(days),max(days),length=101)
X.g <- cbind(rep(1,ng),rep(mean(ozone),ng),days.g)
Z.spline.g <- t(solve(matrix.sqrt.Omega,t(abs(outer(days.g,knots,"-")^3))))

f.hat.g <- as.vector(X.g%*%beta.hat + 
                     Z.spline.g%*%u.hat[1:num.knots.val])

get(getOption("device"))()
plot(days.g,f.hat.g,bty="l",xlab="days",ylab="log(size)",
     type="l",lwd=3,col="maroon")

########## End of sitkaNew ##########