# R-script about graphics for multivariate data 

#############
# data sets #
#############

# we will look at the following four data sets:

# bumpus data:
# all lengths are measured in millimeters
# alar length is length from wing tip to wing tip
# for humerus and sternum, see http://fsc.fernbank.edu/Birding/skeleton.htm
# the first 21 birds survived the storm, the others died
bumpus <- read.table("http://www.ndsu.nodak.edu/ndsu/doetkott/introsas/rawdata/bumpus.html", skip=20, nrows=49, col.names=c("id","total","alar","head","humerus","sternum"))
bumpus <- bumpus[,-1]

# take a first look at the data:
bumpus
dim(bumpus)
summary(bumpus)
boxplot(bumpus, main="Boxplot of Bumpus' data")
pairs(bumpus, panel=panel.smooth, 
      main="Scatterplot matrix of Bumpus' data")

# car data:
?mtcars
# engine displacement = total volume of air/fuel mixture an 
#        engine can draw in during one complete engine cycle
# low rear axel ratio means high towing power

# take a first look at the data:
mtcars
dim(mtcars)
summary(mtcars)
boxplot(mtcars[,1:7], main="Boxplot of car data")
pairs(mtcars[,1:7], panel=panel.smooth,
      main="Scatterplot matrix of car data")

# US judge data:
?USJudgeRatings

# take a first look at the data:
USJudgeRatings
dim(USJudgeRatings)
summary(USJudgeRatings)
boxplot(USJudgeRatings, main="Boxplot of US judge ratings")
pairs(USJudgeRatings, panel=panel.smooth, 
      main="Scatterplot matrix of US judge ratings")

# earth quake data:
?quakes

# take a first look at the data:
dim(quakes)
summary(quakes)
boxplot(quakes)
pairs(quakes, panel=panel.smooth, 
      main="Scatterplot matrix of earth quakes")


##############
# star plots #
##############

# using automatic scaling
stars(bumpus, labels=c(1:49), nrow=6, key.loc=c(20,-1), scale=T,
      main="Star plot of Bumpus' data") 

# re-doing the automatic scaling by hand:
bumpus.scaled.1 <- apply(bumpus, 2, function(x) (x-min(x))/(max(x)-min(x)))
stars(bumpus.scaled.1, labels=c(1:49), nrow=6, key.loc=c(20,-1), scale=F,
      main="Star plot of Bumpus' data") 

# a different way of scaling:
bumpus.scaled.2 <- apply(bumpus, 2, function(x) x/max(x))
stars(bumpus.scaled.2, labels=c(1:49), nrow=6, key.loc=c(20,-1), 
      scale=F, main="Star plot of Bumpus' data") 

# star plot with colors
stars(bumpus, labels=c(1:49), nrow=6, key.loc=c(20,-1), 
      main="Star plot of Bumpus' data (blue=alive, red=dead)", 
      col.stars=c(rep("blue",21),rep("red",28)))

# segment plot
palette(rainbow(12, s = 0.6, v = 0.75)) 
windows()   # open new plotting window (on Windows machines)
stars(bumpus, labels=c(1:49), nrow=6, key.loc=c(20,-1), 
      main="Segment plot of Bumpus' data", draw.segment=TRUE) 

# car data:
# using the option 'full=FALSE':
stars(mtcars[, 1:7], key.loc = c(14, 2), scale=T,
      main = "Motor Trend Cars", full = FALSE)

# segment plot:
stars(mtcars[, 1:7], key.loc = c(14,2), scale=T,
      main = "Motor Trend Cars", draw.segment=TRUE)
    
# US judge ratings (using our own scaling):
USJudge <- apply(USJudgeRatings, 2, function(x) x/max(x))
loc <- stars(USJudge, labels = NULL, scale = FALSE,
             radius = FALSE, frame.plot = TRUE,
             key.loc = c(13, 1.5), main = "Judge data", 
             len = 1.2)
# loc contains the centers of the segment plots

# write name of judge in the middle of each segment plot:
text(loc, abbreviate(row.names(USJudgeRatings)), 
     col = "blue", cex = 0.8, xpd = TRUE)


######################
# conditioning plots #
######################

# car data:
coplot(mpg ~ disp | as.factor(cyl), data = mtcars,
       panel = panel.smooth, rows=1)

# reproduce the plot ourselves:
windows()
par(mfrow=c(1,3))
attach(mtcars)
plot(disp[cyl==4], mpg[cyl==4], xlim=range(disp), 
     ylim=range(mpg), xlab="disp", ylab="mpg")
plot(disp[cyl==6], mpg[cyl==6], xlim=range(disp), 
     ylim=range(mpg), xlab="disp", ylab="mpg")
plot(disp[cyl==8], mpg[cyl==8], xlim=range(disp), 
     ylim=range(mpg), xlab="disp", ylab="mpg")

# earth quake data:

# condition on depth of the earth quake:
coplot(lat ~ long | depth, data = quakes, rows=1)

# determine the conditioning intervals ourselves:
given.depth <- co.intervals(quakes$depth, number=4, overlap=.1)
given.depth
coplot(lat ~ long | depth, data = quakes, 
       given.values=given.depth, rows=1)

# condition on magnitude of the earth quake:
coplot(lat ~ long | mag, data = quakes, rows=1)

# condition on depth and magnitude:
coplot(lat ~ long | depth*mag, data = quakes, rows=1)