### R code from vignette source 'smoothing.Rnw'

###################################################
### code chunk number 1: zero
###################################################
options(show.signif.stars=FALSE,width=66)
set.seed(10131985)


###################################################
### code chunk number 2: kernel1
###################################################
library(alr3)
# kernel regression
par(pch=20, lwd=2)
oldpar <- par(mfrow=c(1,3))
with(oldfaith, {
  plot(Interval ~ Duration, main="bandwidth=5")
  lines(ksmooth(Duration, Interval, "normal", 5))
  plot(Interval ~ Duration, main="bandwidth=20")
  lines(ksmooth(Duration, Interval, "normal", 20))
  plot(Interval ~ Duration, main="bandwidth=40")
  lines(ksmooth(Duration, Interval, "normal", 40))
  })
par(oldpar)


###################################################
### code chunk number 3: kernel1a
###################################################
library(sm)
system.time(with(oldfaith, hcv(Duration, Interval)))
with(oldfaith, hcv(Duration, Interval, display="lines"))


###################################################
### code chunk number 4: kernel2
###################################################
system.time(with(oldfaith,
            print(h.select(Duration, Interval, method="cv"))))
with(oldfaith, sm.regression(Duration, Interval, method="cv"))


###################################################
### code chunk number 5: smoothing.Rnw:65-66 (eval = FALSE)
###################################################
## install.packages("rpanel", dependencies=TRUE)


###################################################
### code chunk number 6: smoothing.Rnw:71-73 (eval = FALSE)
###################################################
## library(rpanel)
## with(oldfaith, sm.regression(Duration, Interval, method="cv", panel=TRUE))


###################################################
### code chunk number 7: loess1
###################################################

(m1 <- loess(Interval~Duration, oldfaith))
with(oldfaith,{
  plot(Interval ~ Duration,
     main="loess, default is f=.75, deg=1")
  lines(spline(Duration, fitted(m1)),lty=1)
  lines(spline(Duration, fitted(update(m1, span=.15))),lty=2, lwd=2)
  lines(spline(Duration, fitted(update(m1, span=.95))),lty=3, lwd=2)
  lines(spline(Duration, fitted(update(m1, degree=1))),lty=4, lwd=2, col="red")
  })
legend("topleft", c("Default", "f=.15", "f=.95", "deg=1"), lty=1:4,
  col=c("black", "black", "black", "red"))


###################################################
### code chunk number 8: smoothing.Rnw:103-104
###################################################
predict(m1,newdata=data.frame(Duration=50*(1:8)))


###################################################
### code chunk number 9: smsplines1
###################################################
print(sm2 <- with(oldfaith, smooth.spline(Duration, Interval)))
with(oldfaith, {
  plot(Interval ~ Duration, main="Smoothing Splines")
  lines(sm2, lty=1, lwd=2)
  lines(update(sm2, spar=.1), lty=2, lwd=2)
  lines(update(sm2, spar=2.5), lty=3, lwd=2)
  })
legend("topleft", c("Default", "spar=.1", "spar=2.5"), lwd=2, lty=1:3, inset=0.02)


###################################################
### code chunk number 10: smoothing.Rnw:132-135
###################################################
mid <- with(oldfaith, (min(Duration) + max(Duration))/2)
StdD <- with(oldfaith, (Duration - mid)/max(Duration - mid))
x <- seq(-1, 1, length=201)


###################################################
### code chunk number 11: poly1
###################################################
par(mfrow=c(1, 2))
plot(x, x, type="l", lty=2, xlab="Std. X", ylab="y",
  main="Polynomial basis", ylim=c(-1, 1), xlim=c(-1, 1.4))
abline(h=0, lty=1)
lines(x,x^2, lty=3)
lines(x,x^3, lty=4)
lines(x,x^4, lty=5)
legend("bottomright", paste("Degree", 0:4), lty=1:5, cex=0.8, inset=0.02)
lm1 <- lm(Interval ~ poly(StdD, 4, raw=TRUE), oldfaith)
plot(Interval ~ StdD, oldfaith, main="Quartic fit",
  ylim=c(-20, 100), xlim=c(-1, 1.4))
lines(sort(StdD), predict(lm1, newdata=data.frame(StdD=sort(StdD))))
b <- coef(lm1)
#plot(x,b[2]*x,type="l",lty=2,xlab="Std. X",ylab="y",main="Polynomial basis",
#  ylim=c(-25,75))
lines(x,b[2]*x,lty=2)
abline(h=b[1],lty=1)
lines(x,b[3]*x^2, lty=3)
lines(x,b[4]*x^3, lty=4)
lines(x,b[5]*x^4, lty=5)
legend("bottomright", paste("Degree",0:4),lty=1:5, cex=.8)


###################################################
### code chunk number 12: rspline1
###################################################
library(splines)
matplot(x, bs(x, df=6), type="l", ylab="", col=1)


###################################################
### code chunk number 13: rspline2
###################################################
bspline.predictors <- bs(StdD, 6)
print(bspline.predictors[1:5,], digits=3)
sm1 <- lm(Interval ~ bs(Duration, 6),oldfaith)
plot(Interval ~ Duration, oldfaith, main="Regression splines")
or <- order(oldfaith$Duration)
lines(predict(sm1)[or]~Duration[or], oldfaith)


###################################################
### code chunk number 14: rspline3
###################################################
data(Orthodont, package="nlme") # get data from nlme package, don't load package
library(lattice)
print(xyplot(distance~age|Sex, data=Orthodont, group=Subject, type=c("g","l")))


###################################################
### code chunk number 15: smoothing.Rnw:190-197
###################################################
library(lme4)
m1 <- lmer(distance ~ Sex*bs(age, 3) + (1 + bs(age, 3)|Subject),
  data=Orthodont)
m2 <- lmer(distance ~ bs(age, 3) + (1 + bs(age, 3)|Subject),
  data=Orthodont)
anova(m1, m2)
AIC(m1, m2)


###################################################
### code chunk number 16: rspline4
###################################################
# fitted value plot
ages <- seq(8, 14, length=100)
fixef(m1)
newages <- model.matrix(~ 1 + bs(ages, 3))
fitmale <- newages %*% fixef(m1)[c(1, 3:5)]
fitfemale <- newages[,c(1, 1, 2, 3, 4, 2, 3, 4)] %*% fixef(m1)
plot(ages, fitmale, type="l", lwd=2, col="red", ylim=c(min(fitmale, fitfemale),
      max(fitmale, fitfemale)), xlab="Age", ylab="Fit distance",
      main="Fitted distance")
lines(ages, fitfemale, type="l", lwd=2, col="blue", lty=2)
legend("topleft", inset=.01, c("Male", "female"), lty=c(1, 2),
      lwd=2, col=c("red", "blue"))


###################################################
### code chunk number 17: rspline5
###################################################
# look at the derivative:
source("http://www.stat.umn.edu/~sandy/courses/8053/Data/mybs.R")
dmale <- mybs(ages, 3, deriv=1) %*% fixef(m1)[3:5]
dfemale <- mybs(ages, 3, deriv=1)[, c(1, 2, 3, 1, 2, 3)] %*% fixef(m1)[3:8]
plot(ages, dmale, type="l", lwd=2, col="red", ylim=c(min(dmale, dfemale),
      max(dmale, dfemale)), xlab="Age", ylab="Deriv. of Fit distance",
      main="Derivative of fitted distance")
lines(ages, dfemale, type="l", lwd=2, col="blue", lty=2)
legend("topleft", inset=.01, c("Male", "female"), lty=c(1,2),
      lwd=2, col=c("red", "blue"))