## Chapter 15: Graphs in R 

## Sec 15.1: Hardcopy graphics devices 
library(DAAG)
## Loading required package: lattice
pdf(file="fossilfuel.pdf", width=6.5, height=6.5, pointsize=18) 
  # For pdf() and postscript(), heights and widths are in inches  
plot(carbon ~ year, data=fossilfuel, pch=16) 
dev.off() 
## pdf 
##   2
par(pty="s")
plot(carbon ~ year, data=fossilfuel, pch=16)  # Display graph 

plot of chunk unnamed-chunk-1

par(pty="m")
## Now open pdf device and copy graph to it 
dev.copy(pdf, file="fossilcopy.pdf") 
## pdf 
##   3
dev.off()            # Close pdf device 
## pdf 
##   2
## Sec 15.2: Plotting characters, symbols, line types and colors 
##                           Font families 
##                               Colors 
##                  Contour and filled contour plots 

## Sec 15.3: Formatting and plotting of text and equations 
## ss 15.3.1: Symbolic substitution of symbols in an expression 
library(DAAG) 
## Replace the symbol tx by its value 
specnam <- "Acmena smithii" 
plot(wood ~ dbh, data=rainforest, subset=species==specnam) 
title(main=substitute(italic(tx) * ": " * "wood vs dbh", 
                      list(tx=specnam))) 

plot of chunk unnamed-chunk-1

## ss 15.3.2: Plotting expressions in parallel 
## Sample from bivariate normal (x, y) with correlation rho 
rho <- 0.75; x <- rnorm(40); y <- rnorm(40) + rho/sqrt(1-rho^2)*x 
## Now map all values of x onto 3 integer values 
x0 <- cut(x, breaks=c(-Inf,-1,1,Inf)) 
plot(unclass(x0), y, xaxt="n") 
axis(3, at=1:3, labels=expression("(-Inf,-1]", "(-1,1]", "(1, Inf]")) 

tiklab <- lapply(strsplit(levels(x0), "Inf"),   
              function(x)if(length(x)>1)substitute(x1*infinity*x2,   
                                   list(x1=x[1],x2=x[2])) else x) 
axis(1, at=1:3, labels=as.expression(tiklab)) 

plot of chunk unnamed-chunk-1

##                 Symbolic substitution in parallel 
 seps <- c(-12.5, -10.523, -8.523, -6.85, -6.6, -3.523, -1, 1, 2.5) 
 plot(range(10^seps), c(0, 1), axes=FALSE, xlab="Wavelength (m)", 
      ylab="", type="n", xaxs="i", yaxs="i", log="x") 
 log10ticks <- pretty(seps,6) 
 log10ticks <- log10ticks[log10ticks > seps[1]] 
 ticklabs <- lapply(round(log10ticks,2),  
                    function(x)substitute(10^a, list(a=x))) 
 axis(1, at=10^log10ticks, labels=as.expression(ticklabs)) 
 ## Use rect() to divide the axis a/c to types of radiation 
 len <- length(seps) 
 rect(xleft=10^seps[-len], ybottom=rep(0, len-1),    
      xright=10^seps[-1], ytop=rep(1, len-1), 
      border=c(NA, rep(1,len-3), NA), 
      col=c(rep("gray70",3), "white", rep("gray70",4))) 
 radtypes <- c("gamma ray","X-ray","ultraviolet", "visible", 
               "infra-red", "microwave", "radio, TV", "long-wave") 
 mid <- 0.5*(seps[-1]+seps[-len]) 
 text(10^mid, rep(0.08, len), radtypes, adj=c(0,0.5), srt=90) 

plot of chunk unnamed-chunk-1

##                  A further use for {substitute()} 
testfun <- function(x, y)deparse(substitute(x)) 
testfun(x = AnyValidName) 
## [1] "AnyValidName"
## Sec 15.4: Multiple graphs on a single graphics page 
par(fig = c(0, 1, 0.38, 1))  # xleft, xright, ybottom, ytop 
## Plot graph A 
par(fig = c(0, 1, 0, 0.38), new=TRUE) 
## Plot graph B 
par(fig = c(0, 1, 0, 1))     # Reset to default 

library(lattice) 
cuckoos.strip <- stripplot(species ~ length, xlab="", data=cuckoos) 
print(cuckoos.strip, position=c(0,0.5,1,1)) 
                       # x, y, x, y; left, bottom, right, top 
cuckoos.bw <- bwplot(species ~ length, xlab="", data=cuckoos) 
print(cuckoos.bw, position=c(0,0,1,0.5), newpage=FALSE) 
  # NB: Use newpage=FALSE in order to retain the first graph 

##          Base and trellis plots on the same graphics page 
par(new=FALSE)

plot of chunk unnamed-chunk-1

plot(0:1, 0:1, type="n", bty="n", axes=FALSE, xlab="", ylab="") 
mtext(side=3, line=3, "Lattice bwplot (i.e., boxplot)") 
cuckoos.bw <- bwplot(species~length, data=cuckoos)  
print(cuckoos.bw, newpage=FALSE) 

plot of chunk unnamed-chunk-1

##              Use of {layout()} -- base graphics only 

## Sec 15.5: Lattice Graphics, and the {grid} Package 
##                 Lattice Graphics vs Base Graphics 
## ss 15.5.1: Groups within data, and/or columns in parallel 
## Simple version of plot 
grogplot <- xyplot(Beer+Spirit+Wine ~ Year | Country, data=grog, 
                   aspect=1, outer=FALSE, auto.key=list(columns=3)) 
## Enhance, and print enhanced code 
update(grogplot, ylim=c(0,5.5), 
       xlab="", ylab="Amount consumed (per person)", 
       par.settings=simpleTheme(pch=c(1,3,4))) 

plot of chunk unnamed-chunk-1

## Footnote Code
## Update trellis object, then print 
frillyplot <- update(grogplot, xlab="", ylab="Amount consumed (per person)", 
                     ylim=c(0,5.5), par.settings=simpleTheme(pch=c(1,3,4))) 
print(frillyplot) 

plot of chunk unnamed-chunk-1

xyplot(Beer+Spirit+Wine ~ Year, groups=Country, outer=TRUE, aspect=1, 
       data=grog, auto.key=list(columns=2) ) 

plot of chunk unnamed-chunk-1

##               Keys -- {auto.key}, {key} and {legend} 
## ss 15.5.2: Lattice parameter settings 
##     Point, line and fill color settings, using {simpleTheme()} 
settings1 <- simpleTheme(pch = c(1,3,4), cex=1.5) 
settings2 <- simpleTheme(pch = c(1,3,4), alpha=0.75) 

grogplot0 <- xyplot(Beer+Spirit+Wine ~ Year | Country, outer=FALSE, 
                    aspect=1, data=grog, ylim=c(0,5.5)) 
grogplot1 <- update(grogplot0, par.settings=settings1) 
  # settings1 are then stored with the object grogplot1 
print(grogplot1) 

plot of chunk unnamed-chunk-1

trellis.par.set(settings2) 

##        Setting that are not available using {simpleTheme()} 
names(trellis.par.get()) 
##  [1] "grid.pars"         "fontsize"          "background"       
##  [4] "panel.background"  "clip"              "add.line"         
##  [7] "add.text"          "plot.polygon"      "box.dot"          
## [10] "box.rectangle"     "box.umbrella"      "dot.line"         
## [13] "dot.symbol"        "plot.line"         "plot.symbol"      
## [16] "reference.line"    "strip.background"  "strip.shingle"    
## [19] "strip.border"      "superpose.line"    "superpose.symbol" 
## [22] "superpose.polygon" "regions"           "shade.colors"     
## [25] "axis.line"         "axis.text"         "axis.components"  
## [28] "layout.heights"    "layout.widths"     "box.3d"           
## [31] "par.xlab.text"     "par.ylab.text"     "par.zlab.text"    
## [34] "par.main.text"     "par.sub.text"
trellis.device(color=FALSE) 
show.settings() 
trellis.device(color=TRUE) 
show.settings() 

trellis.par.set(list(fontsize = list(text = 7, points = 4))) 

##      Parameters that affect axes, tick marks, and axis labels 
## Create a simplified version of the graphics object 
jobsB.xyplot <- 
  xyplot(Ontario+Quebec+BC+Alberta+Prairies+Atlantic ~ Date,  
         data=jobs, type="b", layout=c(3,2), ylab="Number of jobs", 
         scales=list(y=list(relation="sliced", log=TRUE)),  
         outer=TRUE) 
## Update jobsB.xyplot, with various enhancements 
ylabpos <- exp(pretty(log(unlist(jobs[,-7])), 100)) 
ylabels <- paste(round(ylabpos),"\n(", log(ylabpos), ")", sep="")  
## Create a date object 'startofmonth'; use this instead of 'Date' 
atdates <- seq(from=95, by=0.5, length=5) 
datelabs <- format(seq(from=as.Date("1Jan1995", format="%d%b%Y"),   
                   by="6 month", length=5), "%b%y") 
update(jobsB.xyplot, xlab="", between=list(x=0.5, y=0.5), 
       scales=list(x=list(at=atdates, labels=datelabs), 
                   y=list(at=ylabpos, labels=ylabels), tck=0.6) ) 

##                  An example -- the {ais} dataset 
with(ais, table(sex,sport)) 
##    sport
## sex B_Ball Field Gym Netball Row Swim T_400m T_Sprnt Tennis W_Polo
##   f     13     7   4      23  22    9     11       4      7      0
##   m     12    12   0       0  15   13     18      11      4     17
basic1 <-  xyplot(hc ~ rcc | sex, groups=sport[drop=TRUE],  
                  data=subset(ais, sport %in% c("B_Ball", "Swim", 
                                                "Tennis")), 
                  ylab="Blood cell to plasma ratio (%)") 
parSet <- simpleTheme(pch = c(1:3), lty=1:3, lwd=1.5, 
                      col.line=c("gray40","black","black")) 
xaxlab <- expression("Red cell count (10"^{12}*"."*L^{-1}*")") 
basic2 <- update(basic1, par.settings=parSet,  
                 auto.key=list(columns=3, lines=TRUE), 
                 scales=list(tck=0.5), xlab=xaxlab) 
print(update(basic2, type=c("p", "r"))) 

## ss 15.5.3: Panel functions, strip functions, strip labels, and other annotation 
xyplot(species ~ length, xlab="", data=cuckoos) 
xyplot(species ~ length, xlab="", data=cuckoos,  
          panel=panel.xyplot) 

##        A panel function that fits and plots parallel lines 
update(basic2, 
  strip=strip.custom(factor.levels=c("Female","Male")), 
  # In place of level names c("f", "m"), use c("Female", "Male") 
  panel=function(x, y, groups, subscripts, ...){ 
     panel.superpose(x,y, groups=groups, 
                     subscripts=subscripts, ...) 
     ## Obtain fitted values for parallel line model 
     parallel.fitted <- fitted(lm(y ~ groups[subscripts] + x)) 
     panel.superpose(x, parallel.fitted, groups=groups, type="l", 
                     subscripts=subscripts, ...) 
     }) 

library(grid) 
stripplot(species ~ length, xlab="", data=cuckoos, 
          legend=list(top=list(fun=textGrob, 
                      args=list(label="Stripplot", x=0)))) 
  # Here, x=0 is equivalent to x=unit(0,"npc"); the range is (0,1) 

##                  Modification of the strip labels 
tau <- (0:5)/2.5; m <- length(tau); n <- 200; SD <- 2 
x0 <- rnorm(n, mean=12.5, sd=SD)  # Generate x-values 
df <- data.frame(sapply(tau, function(s)x0+rnorm(n, sd=2*s))) 
  # Columns after the first are x-values with addedd error 
df$y = 15+2.5*x0 
names(df) <- c(paste("X", tau, sep=""), "y") 
lab <- c(list("0"), 
         lapply(tau[-1], function(x)substitute(A*s[z], list(A=x)))) 
form <- formula(paste("y ~ ", paste(paste("X", tau, sep=""), 
                                  collapse="+"))) 
xyplot(form, data=df, outer=TRUE, 
  strip=strip.custom(strip.names=TRUE, var.name="SD(added err)", 
          sep=expression(" = "), factor.levels=as.expression(lab))) 

library(reshape) 
longdf <- melt(df, measure.vars=1:6, id.vars="y",  
               variable_name="tau") 
  # Columns of "measure" variables are stacked in a column that has 
  # the name "value": column 1, then 2, then 3, ... 
xyplot(y ~ value | tau, data=longdf, 
       strip=strip.custom(strip.names=TRUE, 
                var.name="SD(added err)", sep=expression(" = "), 
                factor.levels=as.expression(lab))) 

##                   Annotation using {textGrob()} 
library(grid) 
stripplot(species ~ length, xlab="", data=cuckoos, 
          legend=list(top=list(fun=textGrob, 
               args=list(label="Stripplot of cuckoo data", x=0)))) 
  # Here, x=0 is equivalent to x=unit(0,"npc"); the range is (0,1) 

## ss 15.5.4: Interaction with lattice (and other) plots -- the {playwith} package 
## install.packages("playwith", dependencies=TRUE) 

# library(DAAG) 
# library(playwith) 
# playwith(xyplot(age ~ distance, data=hotspots),  
#          labels=hotspots$name) 
# 
# gph <- xyplot(age ~ distance, data=hotspots) 
# playwith(update(gph), labels=hotspots$name) 

## ss 15.5.5: Interaction with lattice plots -- focus, interact,  unfocus 
##      Use of {panel.text()} to label points or add annotation 
# xyplot(Brainwt ~ Bodywt, data=primates) 
# trellis.focus("panel", row=1, column=1, clip.off=TRUE,  
#               highlight=FALSE)  # Non-interactive use 
# xyetc <- trellis.panelArgs() 
# panel.text( x=xyetc$x, y=xyetc$y, labels=row.names(primates), pos=3) 
# trellis.unfocus() 
# 
# stripplot(species ~ length, xlab="", data=cuckoos) 
# trellis.focus(name="toplevel", highlight=FALSE) 
# panel.text("Stripplot of cuckoo data", x=0.5, y=0.98)   
#   # x=0.05 translates to x=unit(0.05,"npc"); the range is (0,1) 
# trellis.unfocus() 

##      Use of {panel.identify()} to label points interactively 
## Use of xyplot(): data frame tinting (DAAG) 
library(lattice) 
# xyplot(it ~ csoa | sex, data=tinting) 
# trellis.focus("panel", column=1, row=1) 
# panel.identify(labels=as.character(tinting$target))     
# ## Now click near points as required. 
# ## Terminate by right clicking inside the panel. 
# ## Now interact with panel 2 
# trellis.focus("panel", row=1, column=2) 
# panel.identify(labels=as.character(tinting$target)) 
## Click, ..., and right click 

## ss 15.5.6: Overlaid plots with different scales 
library(latticeExtra) 
## Loading required package: RColorBrewer
library(DAAG) 
CO2smu <- with(edcCO2, supsmu(age, co2, span=.005)) 
tempsmu <- with(edcT, supsmu(Age, dT, span=.001)) 
co2_layer <- xyplot(y ~ I(-x), data=CO2smu, type="l", 
                    xlab="Years Before Present",  
                    ylab = expression(CO[2]*" (ppm)")) 
temp_layer <- xyplot(I(y+54.5) ~ I(-x), data=tempsmu, type="l", 
                     ylab=expression("Temperature ("*degree*C*")")) 
doubleYScale(co2_layer, temp_layer, add.ylab2 = TRUE) 


## Sec 15.6: An Implementation of Wilkinson's {Grammar of Graphics} 
##                        Australian rain data 
library(DAAG) 
library(ggplot2) 
## 
## Attaching package: 'ggplot2'
## 
## The following object is masked from 'package:latticeExtra':
## 
##     layer
## Default loess smooth, with SE bands added. 
quickplot(Year, mdbRain, data=bomregions, geom=c("point","smooth"), 
          span=0.5, xlab="", se=TRUE, ylab="Av. rainfall, M-D basin") 
## geom_smooth: method="auto" and size of largest group is <1000, so using loess. Use 'method = x' to change the smoothing method.
ggplot(bomregions, aes(x=Year, y=mdbRain)) +  
  geom_point() +                       # Specify a scatterplot 
  geom_smooth(span=0.5, se=TRUE) +     # Add a smooth curve 
  xlab("") +                           # Blank out x-axis label 
  ylab("Av. rainfall, M-D basin")      # Specify y-axis label 
## geom_smooth: method="auto" and size of largest group is <1000, so using loess. Use 'method = x' to change the smoothing method.
## NB: Note use of aes() to supply x- and y-axis variables 

## No legal code:
library(splines) 
library(quantreg) 
## Loading required package: SparseM
## 
## Attaching package: 'SparseM'
## 
## The following object is masked from 'package:base':
## 
##     backsolve
quickplot(Year, mdbRain, data=bomregions, geom=c("point", "quantile"), 
          formula = y ~ ns(x,5), quantiles=c(0.2,0.5,0.8) )  

##            Physical measurements of Australian athletes 
## Overlay scatterplots with boxplots and with density contours 
quickplot(wt, ht, data=ais, geom=c("boxplot", "point", "density2d"), 
          facets = . ~ sex)  

quickplot(wt, ht, xlab="Weight (kg)", ylab="Height (cm)", data=ais, 
          facets = . ~ sex) + 
  geom_boxplot(outlier.size=1.75, outlier.colour="gray",  
               color="gray") + 
  geom_point(shape=2, size=1) + 
  geom_density2d(color="gray")          

## Extract from ais the data for rowers and swimmers 
aisBS <- subset(ais, sport %in% c("Row","Swim")) 
aisBS$sport <- factor(aisBS$sport) 

## Single panel: distinguish sexes by colors; sports by symbols 
## Use quickplot() 
quickplot(wt, ht, data=aisBS, geom="point",  
          color=sex, shape=sport) 
## Follow a call to ggplot() with a call to geom_point() 
ggplot(aisBS) + geom_point(aes(wt, ht, color=sex, shape=sport)) 

## Two panels (sports); sexes have different colors and shapes 
quickplot(wt, ht, data=aisBS, geom="point", size=I(2.5), 
          color=sex, shape=sex, facets = . ~ sport) 
## Single panel: distinguish sexes by colors; sports by shape 
quickplot(wt, ht, data=aisBS, geom="point", 
          color=sex, shape=sport, size=I(2.5)) 
## Single panel: distinguish sexes by color; sports by size 
quickplot(wt, ht, data=aisBS, geom="point", color=sex, size=sport) 

##                   Aesthetic mappings vs settings 
##                 Available geometries and settings 
##                         Themes and updates 
## Set theme_bw() defaults: black gridlines & white background 
## Also set base text size to 8pt 
old <- theme_set(theme_bw(base_size=8)) 
## Reduce default size for geom_point(), from 2 to 1.5  
update_geom_defaults("point", aes(size=1.5)) 
theme_set(old)       # Restore the earlier settings 


## Sec 15.7: Dynamic Graphics -- the {rgl} and {rggobi} packages 
## The Rcmdr and rgl packages must be installed  
# library(Rcmdr)       # This makes scatter3d() available 
# with(nihills, scatter3d(x=log(dist), y=log(climb), z=log(time), 
#                         grid=FALSE, surface.col="gray", 
#                         point.col="black", axis.scales=FALSE)) 
# # with(nihills, identify3d(x=log(dist), y=log(climb), z=log(time), 
# #                          labels=row.names(nihills), col="gray40")) 
# ## NB: Use the middle or right mouse button to drag a rectangle  
# ## around any pony that is to be labeled. 
# 
# open3d(userMatrix= matrix(c(              
#        1, 0.000,  0.000,    0, 
#        0, 0.966, -0.259,    0, 
#        0, 0.259,  0.966,    0, 
#        0, 0.000,  0.000,    1), ncol=4, byrow=TRUE)) 
# par3d(cex=0.6)       # Optional; see help(par3d) for other parameters 


## Sec 15.8: Further Reading