# Visualization

fgl0 <- fgl[ ,-10] # omit type.
fgl.df <- data.frame(type = rep(fgl$type, 9), 
   y = as.vector(as.matrix(fgl0)), 
   meas = factor(rep(1:9, rep(214,9)), labels=names(fgl0)))
trellis.device()
stripplot(type ~ y | meas, data=fgl.df, scales=list(x="free"),
   strip=function(...) strip.default(style=1, ...), xlab="", 
	cex=0.5, pch=18, col=6)
dev.off()

fgl.col <- c("SkyBlue", "SlateBlue", "Orange", "Orchid", "Green", "HotPink")[fgl$type]
xgobi(fgl0, colors=fgl.col)
fgl1 <- fgl0[1:185,]
xgobi(fgl1, colors=fgl.col[1:185])

fgl.sam <- sammon(dist(as.matrix(fgl0[-40, ])))
fgl.iso <- isoMDS(dist(as.matrix(fgl0[-40, ])))


eqscplot(fgl.sam$points, type="n", xlab="", ylab="")
for(i in seq(along=levels(fgl$type))) {
	set <- fgl$type[-40] == levels(fgl$type)[i]
	points(fgl.sam$points[set,], pch=18, cex=0.6, col=2+i)
}
key(text=list(levels(fgl$type), col=3:8))

eqscplot(fgl.iso$points, type="n", xlab="", ylab="")
for(i in seq(along=levels(fgl$type))) {
	set <- fgl$type[-40] == levels(fgl$type)[i]
	points(fgl.iso$points[set,], pch=18, cex=0.6, col=2+i)
}
key(text=list(levels(fgl$type), col=3:8))

# crabs

crabs.grp <- c("B", "b", "O", "o")[rep(1:4, rep(50,4))]

trellis.device()
def <- trellis.par.get("superpose.symbol")
def$pch <- rep(16, 7)
def$cex <- rep(0.5, 7)
trellis.par.set("superpose.symbol", def) 
splom(~ crabs[ ,4:8], groups=crabs.grp, panel=panel.superpose,
    key = list(text = list(c("Blue male", "Blue female", 
                            "Orange Male", "Orange female")),
        points = Rows(trellis.par.get("superpose.symbol"), 1:4), 
        columns = 4)
   )

lcrabs <- log(crabs[,4:8])
lcrabs.pc <- predict(princomp(lcrabs))
trellis.par.set("superpose.symbol", def) 
splom(~ lcrabs.pc[, 1:3], groups = crabs.grp,
    panel = panel.superpose,
    key = list(text = list(c("Blue male", "Blue female", 
                            "Orange Male", "Orange female")),
        points = Rows(trellis.par.get("superpose.symbol"), 1:4),
        columns = 4)
   )

xgobi(lcrabs, colors=c("SkyBlue", "SlateBlue", "Orange", "Red")[rep(1:4, rep(50,4))])

cr1.iso <- isoMDS(dist(as.matrix(lcrabs)))
spcrabs <- scale(lcrabs.pc, , T)
spcr.iso <- isoMDS(dist(as.matrix(spcrabs)))
cr.scale <- 0.5 * log(crabs$CL * crabs$CW)
slcrabs <- lcrabs - cr.scale
cr2.iso <- isoMDS(dist(as.matrix(slcrabs)))

tmp <- apply(cr1.iso$points, 2, function(x)diff(range(x)))
xyplot(cr1.iso$points[,2] ~ cr1.iso$points[,1], groups = crabs.grp,
    panel = panel.superpose, aspect=tmp[2]/tmp[1], xlab="", ylab="",
       key = list(text = list(c("Blue male", "Blue female", 
                    "Orange Male", "Orange female")),
         points = Rows(trellis.par.get("superpose.symbol"), 1:4),
         columns = 4))

tmp <- apply(spcr.iso$points, 2, function(x)diff(range(x)))
xyplot(spcr.iso$points[,2] ~ spcr.iso$points[,1], groups = crabs.grp,
       panel = panel.superpose, aspect=tmp[2]/tmp[1], xlab="", ylab="",
       key = list(text = list(c("Blue male", "Blue female", 
                    "Orange Male", "Orange female")),
         points = Rows(trellis.par.get("superpose.symbol"), 1:4),
         columns = 2))

tmp <- apply(cr2.iso$points, 2, function(x)diff(range(x)))
xyplot(cr2.iso$points[,2] ~ cr2.iso$points[,1], groups = crabs.grp,
    panel = panel.superpose, aspect=tmp[2]/tmp[1], xlab="", ylab="",
       key = list(text = list(c("Blue male", "Blue female", 
                    "Orange Male", "Orange female")),
         points = Rows(trellis.par.get("superpose.symbol"), 1:4),
         columns = 2))
dev.off()

# Density estimation
library(logsplin) # on Windows
x <- seq(8, 35, 0.2)
gal <- galaxies/1000
plot(range(x), c(0, 0.35), type="n",
     xlab="velocity of galaxy", ylab="density")
rug(gal)
lines(x, dlogspline(x, logspline.fit(gal)),col=8)
lines(density(gal, n=200, window="gaussian", width=width.SJ(gal)), col=3)
abline(h=0, col=2)

library(ksmooth)
plot(range(x), c(0, 0.35), type="n",
     xlab="velocity of galaxy", ylab="density")
lines(locpoly(gal, degree=1, bandwidth=0.625), col=8)
lines(locpoly(gal, degree=2, bandwidth=0.625), col=3)
rug(gal)
abline(h=0, col=2)

library(locfit, first=T)
plot(range(x), c(0, 0.35), type="n",
     xlab="velocity of galaxy", ylab="density")
rug(gal)
abline(h=0, col=2)
lines(locfit(~ gal, flim=c(8, 35)), m=500, col=8)

fit <- locfit(~ gal, alpha=0.45, flim=c(8, 35),
        kern="gauss", ev="grid", mg=200)
lines(fit, m=500, col=3)
galaxies.bin <- data.frame(velocity=seq(8, 35, 0.5),
    count=hist(gal, breaks=seq(7750, 35250, 500)/1000,
               plot=F)$counts)
fit2 <- locfit(count ~ velocity, data=galaxies.bin, 
               weights=rep(82*0.5, nrow(galaxies.bin)), 
               alpha=c(0, 0, 2), family="poisson")
lines(fit2, m=500, col=4)


# Smoothing

attach(wtloss)
plot(Days, Weight, type="p",ylab="Weight (kg)")
Wt.lbs <- pretty(range(Weight*2.205))
lines(smooth.spline(Days, Weight), col=4)
lines(loess.smooth(Days, Weight, degree=2, family="symmetric"), col=6)
detach()

attach(mcycle)

plot(times, accel, main="Polynomial regression")
lines(times, fitted(lm(accel ~ poly(times, 6))), col=8)
lines(times, fitted(lm(accel ~ poly(times, 12))), col=3)
legend(40, -100, c("degree=6", "degree=12"), lty=1, col=c(8,3), bty="n")

plot(times, accel, main="Natural splines")
lines(times, fitted(lm(accel ~ ns(times, df=5))), col=8)
lines(times, fitted(lm(accel ~ ns(times, df=10))), col=3)
lines(times, fitted(lm(accel ~ ns(times, df=20))), col=4)
legend(40, -100, c("df=5", "df=10", "df=20"), lty=1, col=c(8,3,4), bty="n")

acc <- unique(mcycle$times)
acc.ns <- ns(acc, df=10)
matplot(acc, acc.ns, type="l", lty=1)
title(main="ns basis")
rug(attr(acc.ns, "knots"))


plot(times, accel, main="Smoothing splines")
lines(smooth.spline(times, accel))
plot(times, accel, main="Lowess")
lines(lowess(times, accel), col=8)
lines(lowess(times, accel, 0.2), col=3)
legend(40, -100, c("default span", "f = 0.2"), lty=1, col=c(8,3), bty="n")

plot(times, accel, main ="ksmooth")
lines(ksmooth(times, accel,"normal", bandwidth=5), col=8)
lines(ksmooth(times, accel,"normal", bandwidth=2), col=3)
legend(40, -100, c("bandwidth=5", "bandwidth=2"), lty=1, col=c(8,3), bty="n")

plot(times, accel, main ="supsmu")
lines(supsmu(times, accel), col=8)
lines(supsmu(times, accel, bass=3), col=3)
legend(40, -100, c("default", "bass=3"), lty=1, col=c(8,3), bty="n")

library(ksmooth)
plot(times, accel, main="locpoly")
lines(locpoly(times, accel, bandwidth=dpill(times,accel)), col=8)
lines(locpoly(times, accel, bandwidth=dpill(times,accel), degree=2), col=3)
legend(40, -100, c("linear", "quadratic"), lty=1, col=c(8,3), bty="n")
detach()

library(locfit, first=T)
fit <- locfit(accel ~ times, alpha = 0.3, data=mcycle)
plot(fit, se.fit=T, get.data=T, main="locfit")


#=== Beyond linear regression =================

# Additive models

rock.lm <- lm(log(perm) ~ area + peri + shape, rock)
summary(rock.lm)
library(mda)
attach(rock)
rock.bruto <- bruto(cbind(area,peri,shape), log(perm))
rock.bruto$type
detach()

set.seed(123)
cpus0 <- cpus[, 2:8]  # excludes names, authors' predictions
for(i in 1:3) cpus0[,i] <- log10(cpus0[,i])
samp <- sample(1:209, 100)
Xin <- as.matrix(cpus0[samp,1:6])
cpus.bruto <- bruto(Xin, log10(cpus0[samp,7]))
# examine the fitted functions
par(mfrow=c(3,2))
Xp <- matrix(sapply(cpus0[samp, 1:6], mean), 100, 6, byrow=T)
for(i in 1:6) {
  xr <- sapply(cpus0, range)
  Xp1 <- Xp; Xp1[,i] <- seq(xr[1,i], xr[2,i], len=100)
  Xf <- predict(cpus.bruto, Xp1)
  plot(Xp1[ ,i], Xf, xlab=names(cpus0)[i], ylab="", type="l")
}

cpus.mars <- mars(Xin, log10(cpus0[samp,7]))
showcuts(cpus.mars)
# examine the fitted functions
Xp <- matrix(sapply(cpus0[samp, 1:6], mean), 100, 6, byrow=T)
for(i in 1:6) {
  xr <- sapply(cpus0, range)
  Xp1 <- Xp; Xp1[,i] <- seq(xr[1,i], xr[2,i], len=100)
  Xf <- predict(cpus.mars, Xp1)
  plot(Xp1[ ,i], Xf, xlab=names(cpus0)[i], ylab="", type="l")
}
par(mfrow=c(1,1))
cpus.mars2 <- mars(Xin, log10(cpus0[samp,7]), degree=2)
showcuts(cpus.mars2)
cpus.mars6 <- mars(Xin, log10(cpus0[samp,7]), degree=6)
showcuts(cpus.mars6)

# Response transformation

attach(cpus0)
cpus.ace <- ace(cpus0[, 1:6], perf)
par(mfrow=c(1,2))
plot(perf, cpus.ace$ty)
plot(log10(perf), cpus.ace$ty)
par(mfrow=c(2,3))
for(i in 1:6)
  plot(cpus0[, i], cpus.ace$tx[, i], xlab=names(cpus0[i]), ylab="")

cpus.avas <- avas(cpus0[, 1:6], perf)
par(mfrow=c(1,2))
plot(perf, cpus.avas$ty)
plot(log10(perf), cpus.avas$ty)
par(mfrow=c(2,3))
for(i in 1:6)
  plot(cpus0[, i], cpus.avas$tx[, i], xlab=names(cpus0[i]), ylab="")
par(mfrow=c(1,1))
detach()

# Projection pursuit regression

library(ppr)
attach(rock)
rock1 <- data.frame(area=area/10000, peri=peri/10000,
                    shape=shape, perm=perm)
detach()
rock.ppr <- ppr(log(perm) ~ area + peri + shape, data=rock1,
                  nterms=2, max.terms=5)
rock.ppr
summary(rock.ppr)
par(mfrow=c(3,2))
plot(rock.ppr, main="default")
plot(update(rock.ppr, bass=5), main="bass=5")
rock.ppr2 <- update(rock.ppr, sm.method="gcv", gcvpen=2)
plot(rock.ppr2, main="Smoothing splines")
summary(rock.ppr2)
summary(rock1) # to find the ranges of the variables
Xp <- expand.grid(area=seq(0.1,1.2,0.05),
                  peri=seq(0,0.5,0.02), shape=0.2)
par(mfrow=c(1,1))
trellis.device()
rock.grid <- cbind(Xp,fit=predict(rock.ppr2, Xp))
wireframe(fit ~ area+peri, rock.grid, screen=list(z=160,x=-60),
          aspect=c(1,0.5), drape=T)
dev.off()


# ==== Neural networks =====

x <- 1:250 *pi/250
f <- 4.26*(exp(-x) - 4 * exp(-2*x) + 3 * exp(-3*x))
set.seed(111)
y <- rnorm(250, f, 0.2)

plot(x, y, type="n")
points(x, y, cex=0.5, pch=16)
lines(x, lm(y ~ poly(x,6))$fitted, col=3)
lines(x, lm(y ~ poly(x,10))$fitted, col=4)
lines(x, lm(y ~ poly(x,16))$fitted, col=5)

plot(x, y, type="n")
points(x, y, cex=0.5, pch=16)
# try each of these fits several times
lines(x, nnet(y ~x, size=2, skip=T, linout=T, maxit=250)$fitted, col=3)
lines(x, nnet(y ~x, size=4, skip=T, linout=T, maxit=250)$fitted, col=4)
lines(x, nnet(y ~x, size=8, skip=T, linout=T, maxit=250)$fitted, col=5)

# now add weight decay
lines(x, nnet(y ~x, size=8, decay=0.01, skip=T, linout=T, maxit=250)$fitted, col=6)



attach(rock)
area1 <- area/10000; peri1 <- peri/10000
rock1 <- data.frame(perm, area=area1, peri=peri1, shape)
rock.nn <- nnet(log(perm) ~ area + peri + shape, data=rock1,
                size=3, decay=1e-3, linout=T, skip=T, maxit=1000, Hess=T)
summary(rock.nn)
sum((log(perm) - predict(rock.nn))^2)
detach()
eigen(rock.nn$Hess, T)$values


Xp <- expand.grid(area=seq(0.1,1.2,0.05),
                  peri=seq(0,0.5,0.02), shape=0.2)
trellis.device()
rock.grid <- cbind(Xp,fit=predict(rock.nn, Xp))
wireframe(fit ~ area + peri, rock.grid, screen=list(z=160,x=-60),
          aspect=c(1,0.5), drape=T)
dev.off()


pltnn("Logistic regression")
plt.bndry(size=0, col=2)

#par(mfrow=c(2,2))
pltnn("Size = 2")
set.seed(1); plt.bndry(size=2, col=2)
set.seed(3); plt.bndry(size=2, col=3); plt.bndry(size=2, col=4)

pltnn("Size = 2, lambda = 0.001")
set.seed(1); plt.bndry(size=2, decay=0.001, col=2)
set.seed(2); plt.bndry(size=0, decay=0.001, col=4)

pltnn("Size = 2, lambda = 0.01")
set.seed(1); plt.bndry(size=2, decay=0.01, col=2)
set.seed(2); plt.bndry(size=2, decay=0.01, col=4)

pltnn("Size = 5, 20  lambda = 0.01")
set.seed(2); plt.bndry(size=5, decay=0.01, col=1)
set.seed(2); plt.bndry(size=20, decay=0.01, col=2)

par(mfrow=c(1,1))
pltnn("Size=3, lambda=0.01")
Z <- matrix(0, nrow(cushT), ncol(tpi))
for(iter in 1:20) {
    set.seed(iter)
    cush.nn <- nnet(cush, tpi,  skip=T, softmax=T, size=3,   
        decay=0.01, maxit=1000, trace=F)
    Z <- Z + predict(cush.nn, cushT)
    cat("final value", format(round(cush.nn$value,3)), "\n")
    b1(predict(cush.nn, cushT), col=2, lwd=0.5)
    guiLocator(0)
}
pltnn("Averaged")
b1(Z, lwd=3)


# 2-class synthetic data problem


attach(synth.tr)
synth.set()
synth.fit(size=0)

# try this several times
synth.fit(size=6, maxit=500)

synth.fit(size=6, decay=0.001, col=3, maxit=500)

synth.fit(size=6, decay=0.01, col=7, maxit=500)

synth.aver(size=3, decay=0.001, naver=10, maxit=500)



#======== Trees ==========================

cush.tr <- tree(tp ~ Tetrahydrocortisone + Pregnanetriol, Cf)
cush.tr
plot(cush[,1], cush[,2], type="n", main="Tree",
  xlab="Tetrahydrocortisone", ylab = "Pregnanetriol")
  for(il in 1:4) {
    set <- Cushings$Type == levels(Cushings$Type)[il]
    text(log(Cushings[set, 1]), log(Cushings[set, 2]),
         as.character(Cushings$Type[set]), col = 2 + il)
  }
par(cex=1.5); partition.tree(cush.tr, add=T); par(cex=1)


ir.tr <- tree(Species ~., ird)
summary(ir.tr)
ir.tr
plot(ir.tr)
text(ir.tr, all=T)
ir.tr1 <- snip.tree(ir.tr)
ir.tr1
summary(ir.tr1)

par(pty="s")
plot(ird[, 3],ird[, 4], type="n", xlab="petal length", ylab="petal width")
text(ird[, 3], ird[, 4], as.character(ird$Species), cex=0.7)
par(cex=2);partition.tree(ir.tr1, add=T);par(cex=1)
par(pty="m")

fgl.tr <- tree(type ~ ., fgl)
summary(fgl.tr)
plot(fgl.tr);  text(fgl.tr, all=T, cex=0.5)
fgl.tr1 <- snip.tree(fgl.tr)
tree.screens()
plot(fgl.tr1)
tile.tree(fgl.tr1, fgl$type)
close.screen(all = T)

par(mfrow=c(1,2), pty="s")
set.seed(123)
fgl.cv <- cv.tree(fgl.tr,, prune.tree)
for(i in 2:5)
  fgl.cv$dev <- fgl.cv$dev + cv.tree(fgl.tr,, prune.tree)$dev
fgl.cv$dev <- fgl.cv$dev/5
plot(fgl.cv, order="decreasing")

set.seed(123)
fgl.cv <- cv.tree(fgl.tr,, prune.misclass)
for(i in 2:5)
  fgl.cv$dev <- fgl.cv$dev + cv.tree(fgl.tr,, prune.misclass)$dev
fgl.cv$dev <- fgl.cv$dev/5
fgl.cv
plot(fgl.cv, order="decreasing")

par(mfrow=c(1,1), pty="m")

# RPart ==============

library(rpart)
ir.rp <- rpart(Species ~ ., data=ird, cp=0, minsplit=5, maxsurrogate=0)
ir.rp
plot(ir.rp);text(ir.rp)

plotcp(ir.rp)
printcp(ir.rp)
print(ir.rp, cp=0.015)
print(ir.rp, cp=0.094)

set.seed(123)
fgl.rp <- rpart(type ~ ., fgl, cp=0.001)
plotcp(fgl.rp)
printcp(fgl.rp)
print(fgl.rp, cp=0.02)
fgl.rp2 <- prune(fgl.rp, cp=0.02)
plot(fgl.rp2); text(fgl.rp2)



#======== Refining lda  ===================

tp
lda(cush, tp, CV=T)

guiModify("GraphSheet", Name = guiGetGSName(),
	NumImageColors = "2",
	NumImageShades = "16",
	ImageColorTable = "108,108,108|255,255,255|255,255,255|0,255,255|255,255,255|255,255,255|255,255,255|255,255,255|255,255,255|255,255,255|255,255,255|255,255,255|255,255,255|255,255,255|255,255,255|255,255,255")

cush.lda <- lda(cush, tp); predplot(cush.lda, "LDA")

cush.qda <- qda(cush, tp)
predplot(cush.qda, "QDA")
predplot(cush.qda, "QDA (debiased)", method="debiased")
predplot(cush.qda, "QDA (predictive)", method="predictive")