Code For Graph Gallery

# R code to generate graphs on www.milbo.users.sonic.net/gallery

#--- plotmo ---------------------------------------------------------------

library(earth)
data(ozone1)
a <- earth(O3 ~ ., data=ozone1, degree=2)
plotmo(a, caption="plotmo example 1: MARS model")

library(gam)
data(airquality)
airquality <- na.omit(airquality)  # plotmo doesn't know how to deal with NAs
plotmo(gam(Ozone^(1/3) ~ lo(Solar.R) + lo(Wind, Temp), data=airquality),
       caption="plotmo example 2: GAM model",
       pch=20, pt.col="blue",
       type2="image", ngrid2=100, image.col=terrain.colors(100))

plotmo(lm(log(O3) ~ humidity*temp, data=ozone1),
       caption="plotmo example 3", level=.95, clip=FALSE)

#--- plotd ----------------------------------------------------------------

library(earth)
data(etitanic)
fit <- earth(survived ~ ., data=etitanic, degree=2, glm=list(family=binomial))

plotd(fit, main="plotd example 1", type="response")

plotd(fit, main="plotd example 2", type="response", vline.col="gray",
      err.col=c("slategray1","slategray3","pink"))

plotd(fit, main="plotd example 3", type="response", hist=TRUE,
      legend.pos=c(.25,220))

plotd(fit, main="plotd example 4", type="class", hist=TRUE,
      labels=TRUE, xlab="", breaks=4)

#--- plotpc ---------------------------------------------------------------

library(plotpc)
data(iris)
x <- iris[,c(3,4)] # select Petal.Length and Petal.Width
plotpc(x, main="plotpc example 1\nPrincipal component histograms\n")

# example with some parameters and showing densities
plotpc(x,
       main="plotpc example 2\nPrincipal component densities\n",
       hist=FALSE,          # plot densities not histograms
       gp.points=gpar(col=2, cex=.6), # small red points
       adjust=.5,           # finer resolution in the density plots
       gp.axis=gpar(lty=3), # gpar of axes
       heightx=0,           # don't display x histogram
       heighty=0,           # don't display y histogram
       text1="",            # text above hist for 1st principal component
       text2="",            # text above hist for 2nd principal component
       axis.len2=4,         # length of 2nd principal axis (in std devs)
       offset1=2.5,         # offset of component 1 density plot
       offset2=5)           # offset of component 2 density plot

# example showing principal axes
x <- iris[iris$Species=="versicolor",c(3,4)]
vp <- plotpc(x,
       main="plotpc example 3\nPrincipal axes with confidence ellipse",
       sd.ellipse=2,                       # ellipse at two standard devs
       heightx=0, heighty=0, height1=0, height2=0, # no histograms
       gp.ellipse=gpar(col=1),             # ellipse in black
       axis.lenx=4, axis.leny=5,           # lengthen horiz and vertical axes
       axis.len1=4, gp.axis1=gpar(col=2),  # lengthen pc1 axis, draw in red
       axis.len2=8, gp.axis2=gpar(col=2))  # lengthen pc2 axis, draw in red

pushViewport(vp) # add text to the graph
un <- function(x) unit(x, "native")
grid.text("PC1", x=un(2.2), y=un(.6),   gp=gpar(cex=.8, font=2, col=2))
grid.text("PC2", x=un(3.9), y=un(2.35), gp=gpar(cex=.8, font=2, col=2))
grid.text("X1",  x=un(2.2), y=un(1.4),  gp=gpar(cex=.8, font=2))
grid.text("X2",  x=un(4.3), y=un(2.5),  gp=gpar(cex=.8, font=2))
popViewport()

# example comparing linear regression to principal axis
x <- iris[iris$Species=="setosa",c(3,4)]
library(alr3) # get banknote data
data(banknote)
x <- banknote[101:200,4:5] # select Forged, Top and Bottom
vp <- plotpc(x,
       main="plotpc example 4\nRegression lines and\nfirst principal component",
       heightx=0, heighty=0, height1=0, height2=0, # no histograms
       gp.points=gpar(col="steelblue"),      # color of points
       axis.len1=4,  gp.axis1=gpar(col=2, lwd=2),
       axis.len2=0,                          # no pc2 axis
       yonx=TRUE, xony=TRUE)                 # display regression lines

pushViewport(vp) # add text to the principal component line
grid.text("PC1", x=unit(6, "native"), y=unit(13.5, "native"),
          gp=gpar(col=2, cex=.8, font=2))
popViewport()

# Flury and Riedwyl "Multivar Stats A Practical Approach" Figure 7.1
library(plotpc)
library(alr3) # get banknote data
data(banknote)
x <- banknote[,4:5] # select Top and Bottom
plotpc(x, xrange=24,
       main="plotpc example 5\nVarious projections\n(Based on Flury and Riedwyl Figure 7.1)",
       height=-2,  # reverse directions of histograms to match Flury and Riedwyl
       gp.points=gpar(col="steelblue", cex=.5),
       breaks=12, heightx=0, heighty=0, height1=0, height2=0,
       axis.len1=0, axis.len2=0,
       angle3=30, angle4=102, angle5=174, angle6=246, angle7=318)

#--- plotld ---------------------------------------------------------------

library(plotpc)
data(iris)
x <- iris[, -5] # -5 to drop Species
plotld(x, main="principal component loadings")