## 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
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)))
## 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))
## 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)
## 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(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)
## 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)))
## 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)
xyplot(Beer+Spirit+Wine ~ Year, groups=Country, outer=TRUE, aspect=1,
data=grog, auto.key=list(columns=2) )
## 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)
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