graphics.txt

里面有三个PDF文件

1>Figure 1.1 A simple scatterplot 
plot(pressure)
text(150, 600, "Pressure (mm Hg)\nversus\nTemperature (Celsius)")

#################################################################
2>Figure 1.2 Some standard plots 

#
#  Comment:
# 
#  Examples of the use of standard high-level plotting functions.
# 
#  In each case, extra output is also added using low-level 
#  plotting functions.
#


par(mfrow=c(3, 2))

# Scatterplot
x <- c(0.5, 2, 4, 8, 12, 16)
y1 <- c(1, 1.3, 1.9, 3.4, 3.9, 4.8)
y2 <- c(4, .8, .5, .45, .4, .3)
par(las=1, mar=c(4, 4, 2, 4))
plot.new()
plot.window(range(x), c(0, 6))
lines(x, y1)
lines(x, y2)
points(x, y1, pch=16, cex=2)
points(x, y2, pch=21, bg="white", cex=2)
par(col="grey50", fg="grey50", col.axis="grey50")
axis(1, at=seq(0, 16, 4))
axis(2, at=seq(0, 6, 2))
axis(4, at=seq(0, 6, 2))
box(bty="u")
mtext("Travel Time (s)", side=1, line=2, cex=0.8)
mtext("Responses per Travel", side=2, line=2, las=0, cex=0.8)
mtext("Responses per Second", side=4, line=2, las=0, cex=0.8)
text(4, 5, "Bird 131")
par(mar=c(5.1, 4.1, 4.1, 2.1), col="black", fg="black", col.axis="black")

# Histogram
# Random data
Y <- rnorm(50)
# Make sure no Y exceed [-3.5, 3.5]
Y[Y < -3.5 | Y > 3.5] <- NA
x <- seq(-3.5, 3.5, .1)
dn <- dnorm(x)
par(mar=c(4.5, 4.1, 3.1, 0))
hist(Y, breaks=seq(-3.5, 3.5), ylim=c(0, 0.5), 
     col="grey80", freq=FALSE)
lines(x, dnorm(x), lwd=2)
par(mar=c(5.1, 4.1, 4.1, 2.1))

# Barplot
# Modified from example(barplot)
par(mar=c(2, 3.1, 2, 2.1))
midpts <- barplot(VADeaths, col=grey(0.5 + 1:5/12), 
                  names=rep("", 4))
mtext(sub(" ", "\n", colnames(VADeaths)),
      at=midpts, side=1, line=0.5, cex=0.5)
text(rep(midpts, each=5), apply(VADeaths, 2, cumsum) - VADeaths/2,
     VADeaths, col=rep(c("white", "black"), times=2:3, cex=0.8))
par(mar=c(5.1, 4.1, 4.1, 2.1))

# Boxplot
# Modified example(boxplot) - itself from suggestion by Roger Bivand
par(mar=c(3, 4.1, 2, 0))
     boxplot(len ~ dose, data = ToothGrowth,
             boxwex = 0.25, at = 1:3 - 0.2,
             subset= supp == "VC", col="grey90",
             xlab="",
             ylab="tooth length", ylim=c(0,35))
     mtext("Vitamin C dose (mg)", side=1, line=2.5, cex=0.8)
     boxplot(len ~ dose, data = ToothGrowth, add = TRUE,
             boxwex = 0.25, at = 1:3 + 0.2,
             subset= supp == "OJ", col="grey70")
     legend(1.5, 9, c("Ascorbic acid", "Orange juice"), bty="n",
            fill = c("grey90", "grey70"))
par(mar=c(5.1, 4.1, 4.1, 2.1))

# Persp
# Almost exactly example(persp)
    x <- seq(-10, 10, length= 30)
     y <- x
     f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
     z <- outer(x, y, f)
     z[is.na(z)] <- 1
# 0.5 to include z axis label
par(mar=c(0, 0.5, 0, 0), lwd=0.1)
     persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "grey80")
par(mar=c(5.1, 4.1, 4.1, 2.1), lwd=1)

# Piechart
# Example 4 from help(pie)
par(mar=c(0, 2, 1, 2), xpd=FALSE, cex=0.5)
     pie.sales <- c(0.12, 0.3, 0.26, 0.16, 0.04, 0.12)
     names(pie.sales) <- c("Blueberry", "Cherry",
         "Apple", "Boston Cream", "Other", "Vanilla")
     pie(pie.sales, col = gray(seq(0.4,1.0,length=6)))
#################################################################
3>Figure 1.3  A customized scatterplot 
#
# Comment:
#
# A sophisticated example of adding further output to a basic plot.
# 
# Most of the functions defined are just for calculating values
# relevant to the data analysis.  
# 
# The function plotPars() is the one of interest for seeing how
# the drawing of the plot is done.
#


params <- function(N, breaks, p=seq(0.001, 1, length=100)) {
  list(N=N, T=1/breaks, p=p, q=1-p)
}

pdfcomp <- function(comp, params) {
  n <- params$T
  p <- params$p
  q <- params$q
  y <- round(comp/n)
  choose(n, comp)*p^comp*q^(n-comp) / (1 - q^n)
}

# Expected num sherds (for a vessel) [=completeness]
expcomp <- function(params) {
  params$T*params$p/(1-params$q^params$T)
}

# Variance of num sherds (for a vessel)
varcomp <- function(params) {
  n <- params$T
  p <- params$p
  q <- params$q
  # From Johnson & Kotz
  (n*p*q / (1 - q^n)) - (n^2*p^2*q^n / (1 - q^n)^2)
  # n^2 times Thomas Yee's formula
  # n^2*((p*(1 + p*(n - 1)) / (n*(1 - q^n))) - (p^2 / (1 - q^n)^2))
}

# Expected value of completeness (for a sample of vessels)
expmeancomp <- function(params) {
  expcomp(params)
}

# Variance of completeness (for a sample of vessels)
# Use the expected number of vessels in sample as denominator
varmeancomp <- function(params) {
  varcomp(params)/(numvess(params))
}

numvess <- function(params) {
  params$N*(1-params$q^params$T)
}

ecomp <- function(p, T, comp) {
  q <- 1 - p
  T*p/(1 - q^T) - comp
}

estN <- function(comp, broke, n) {
  T <- 1/broke
  n / (1 - (1 - uniroot(ecomp, c(0.00001, 1), T=T, comp=comp)$root)^T)
}

nvessscale <- function(params, xlim, ylim, new=TRUE) {
  if (new)
    par(new=TRUE)
  plot(0:1, c(1, params$N), type="n", axes=!new, ann=FALSE,
       xlim=xlim, ylim=ylim)
}

compscale <- function(params, xlim, ylim, new=TRUE) {
  if (new)
    par(new=TRUE)
  plot(0:1, c(1, params$T), type="n", axes=!new, ann=FALSE,
       xlim=xlim, ylim=ylim)
}

lowerCI <- function(p, N, breaks, lb) {
  params <- params(N, breaks, p)
  expmeancomp(params) - 2*sqrt(varmeancomp(params)) - lb
}

upperCI <- function(p, N, breaks, lb) {
  params <- params(N, breaks, p)
  expmeancomp(params) + 2*sqrt(varmeancomp(params)) - lb
}

critP <- function(comp, params) {
  c(uniroot(lowerCI, c(0.00001, 1), N=params$N,
            breaks=1/params$T, lb=max(comp))$root,
    if (upperCI(0.00001, params$N, 1/params$T, min(comp)) > 0) 0
    else uniroot(upperCI, c(0.00001, 1), N=params$N,
                 breaks=1/params$T, lb=min(comp))$root)
}

anncomp <- function(params, comp, xlim, ylim, cylim) {
  cp <- critP(comp, params)
  nv <- numvess(params(params$N, 1/params$T, cp))
  nvessscale(params, xlim, ylim)
  polygon(c(cp[2], cp[2], 0, 0, cp[1], cp[1]),
          c(0, nv[2], nv[2], nv[1], nv[1], 0),
          col="grey90", border=NA)
  text(0, nv[1], paste(round(nv[1]),
                       " (", round(100*nv[1]/params$N), "%)", sep=""),
       adj=c(0, 0), col="grey")
  text(0, nv[2], paste(round(nv[2]), 
                       " (", round(100*nv[2]/params$N), "%)", sep=""),
       adj=c(0, 1), col="grey")
  compscale(params, xlim, cylim)
  segments(1, min(comp), cp[2], comp, col="grey")
  segments(1, max(comp), cp[1], comp, col="grey")
  text(1, comp, paste(comp, collapse="-"), adj=c(1, 0), col="grey")
}

plotPars <- function(params, comp, xlim=NULL, ylim=NULL) {
  mean <- expmeancomp(params)
  var <- 2*sqrt(varmeancomp(params))
  lb <- mean - var
  ub <- mean + var
  par(mar=c(5, 4, 4, 4))
  if (is.null(ylim))
    cylim <- ylim
  else
    cylim <- c(1 + ((ylim[1] - 1)/(params$N - 1))*(params$T - 1),
               1 + ((ylim[2] - 1)/(params$N - 1))*(params$T - 1))
  nvessscale(params, xlim, ylim, new=FALSE)
  compscale(params, xlim, cylim)
  polygon(c(params$p, rev(params$p)), c(lb, rev(ub)),
          col="grey90", border=NA)
  anncomp(params, comp, xlim, ylim, cylim)
  nvessscale(params, xlim, ylim)
  mtext("Number of Vessels", side=2, line=3)
  mtext("Sampling Fraction", side=1, line=3)
  lines(params$p, numvess(params))
  par(new=TRUE)
  compscale(params, xlim, cylim)
  mtext("Completeness", side=4, line=3)
  axis(4)
  lines(params$p, mean, lty="dashed")
  lines(params$p, lb, lty="dotted")
  lines(params$p, ub, lty="dotted")
  mtext(paste("N = ", round(params$N),
              "     brokenness = ", round(1/params$T, 3), sep=""),
        side=3, line=2)
}

par(cex=0.8, mar=c(3, 3, 3, 3))
p6 <- params(estN(1.2, 0.5, 200), 0.5)
plotPars(p6, 1.2)
nvessscale(p6, NULL, NULL)
pcrit <- 1 - (1 - 200/estN(1.2, 0.5, 200))^(1/p6$T)
lines(c(0, pcrit), c(200, 200))
lines(c(pcrit, pcrit), c(200, 0))
##################################################################
4>Figure 1.4 
A Trellis dotplot 
#
# Comment:
#
# A slightly modified version of Figure 1.1 from 
# Cleveland's book "Visualizing Data"
#


library(lattice)
trellis.par.set(theme = canonical.theme("postscript", col=FALSE))
trellis.par.set(list(fontsize=list(text=6),
	             par.xlab.text=list(cex=1.5),
                     add.text=list(cex=1.5),
                     superpose.symbol=list(cex=.5)))
key <- simpleKey(levels(barley$year), space = "right")
key$text$cex <- 1.5
print(
     dotplot(variety ~ yield | site, data = barley, groups = year,
             key = key,
             xlab = "Barley Yield (bushels/acre) ",
             aspect=0.5, layout = c(1,6), ylab=NULL)
)

###################################################################
5>Figure 1.5 
A map of New Zealand produced using R 
#
# Comment:
#
# A bit of mucking around is required to get the second (whole-world)
# map positioned correctly;  this provides an example of calling a 
# plotting function to perform calculations but do no drawing (see the
# second call to the map() function).
#
# Makes use of the "maps" and "mapproj" packages to draw the maps.
#


library(maps)
par(mar=rep(0, 4))
map("nz", fill=TRUE, col="grey80")
points(174.75, -36.87, pch=16, cex=2)
arrows(172, -36.87, 174, -36.87, lwd=3)
text(172, -36.87, "Auckland  ", adj=1, cex=2)
# mini world map as guide
maplocs <- map(projection="sp_mercator", wrap=TRUE, lwd=0.1, 
               col="grey", ylim=c(-60, 75),
               interior=FALSE, orientation=c(90, 180, 0), add=TRUE,
               plot=FALSE)
xrange <- range(maplocs$x, na.rm=TRUE)
yrange <- range(maplocs$y, na.rm=TRUE)
aspect <- abs(diff(yrange))/abs(diff(xrange))
# customised to 6.5 by 4.5 figure size
par(fig=c(0.99 - 0.5, 0.99, 0.01, 0.01 + 0.5*aspect*4.5/6.5), 
    mar=rep(0, 4), new=TRUE)
plot.new()
plot.window(xlim=xrange,
            ylim=yrange)
map(projection="sp_mercator", wrap=TRUE, lwd=0.1, ylim=c(-60, 75),
    interior=FALSE, orientation=c(90, 180, 0), add=TRUE)
symbols(-.13, -0.8, circles=1, inches=0.1, add=TRUE)

##################################################################
6>Figure 1.6 
Some polar-coordinate plots 
"polar.plot" <- 
function (r, theta, theta.zero = 0, theta.clw = FALSE, method = 1, 
    rlabel.axis = 0, dir = 8, rlimits = NULL, grid.circle.pos = NULL, 
    grid.lwd = 1, grid.col = "black", points.pch = 20, points.cex = 1, 
    lp.col = "black", lines.lwd = 1, lines.lty = 1, polygon.col = NA, 
    polygon.bottom = TRUE, overlay = NULL, pi2.lab = TRUE, text.lab = NULL, 
    num.lab = NULL, rlabel.method = 1, rlabel.pos = 3, rlabel.cex = 1, 
    rlabel.col = "black", tlabel.offset = 0.1, tlabel.cex = 1.5, 
    tlabel.col = "black", main = NULL, sub = NULL) 
{ 
# r: (vector of) radial data. 
# theta: (vector of) angular data (in radians). 
# theta.zero: angular direction on plot of theta = 0 (in radians). 
# theta.clw: clockwise orientation of theta values (default = FALSE). 
# 
# method: (plotting of (r,theta)-data): 
# 1: points (default) 
# 2: line 
# 3: polygon 
# 
# rlabel.axis: angular direction on the plot of radial label axis (in radians). 
# dir: number of radial grid lines (default=8). 
# rlimts: Interval for radial axis as a numeric vector: c(lower,upper). Interval will be extended by the default use of pretty()-function. (default = NULL). 
# grid.circle.pos: radial axis position of grid circles as numeric vector of minimum length 2. Overrides the default positioning of grid circles by pretty()-function. (default = NULL). 
# grid.lwd. grid line width. 
# grid.col: grid line color. 
# 
# points.pch: points plotting symbol. 
# point.cex: character expansion factor for points. 
# lp.col: color of points (method 1) or lines (method 2 and method 3). In method 3, set lp.col=0 for polygons without border. 
# lines.lwd: line width for plotting methods 2 and 3 (default = 1). 
# lines.lty: line type (default = 1). 
# polygon.col: color of polygon (defalut = NA). 
# polygon.bottom: polygon to the back i.e. behind the grid (default = TRUE). 
# 
# overlay: NULL (default), no overlay 
# 1, overlay data on existing plot 
# 2, overlay data, grid and labels on existing plot. 
# 
# pi2.lab: angular labels in radians (0, pi/2, pi, 3*pi/2) (default). 
# text.lab: angular axis labels from a character vector c("N","E","S","W") (default = NULL). 
# num.lab: numeric angular axis labels in interval [0;num.lab[ (default = NULL). Number of labels: dir. 
# 
# rlabel.method (plotting of radial axis labels): 
# 0: no radial labels. 
# 1: labels at pretty radial distances (default). 
# 2: exclude label at radial distace 0. 
# 3: exclude label at maximum radial distance. 
# 4: exclude radial labels at distance 0 and at maximum radial distance. 
# rlabel.pos: text position of radial axis labels (NULL,1,2,3,4). 
# rlabel.cex: cex for radial axis labels. 
# rlabel.col: color of the radial labels. 
# 
# tlabel.offset: radial offset for angular axis labels in fraction of maximum radial value (default = 0.1). 
# tlabel.cex: cex for angular axis labels. 
# tlabel.col: angular labels color. 
# 
# main: plot main title. 
# sub: plot sub title. 
   
    fit.rad <- function(x, twop = 2 * pi) { 
        for (i in 1:length(x)) { 
            while (x[i] < 0) x[i] <- x[i] + twop 
            while (x[i] >= twop) x[i] <- x[i] - twop 
        } 
        return(x) 
    } 
    if (is.null(rlimits)) 
        rpretty <- pretty(range(abs(r), 0, na.rm = TRUE)) 
    if (is.numeric(rlimits) & length(rlimits) == 2) 
        rpretty <- pretty(range(abs(rlimits[1]), abs(rlimits[2]))) 
    if (is.numeric(grid.circle.pos) & length(grid.circle.pos) > 
        1) 
        rpretty <- grid.circle.pos 
    lab.dist <- max(rpretty) 
    if (!is.null(text.lab) || is.numeric(num.lab) || pi2.lab) { 
        lab.dist <- lab.dist * (tlabel.offset + 1) 
    } 
    if (is.null(overlay)) { 
        plot.new() 
        ps <- max(lab.dist, max(rpretty)) 
        plot.window(xlim = c(-ps, ps), ylim = c(-ps, ps), asp = 1) 
        title(main = main, sub = sub) 
    } 
    drawgrid <- function() { 
        if (dir > 0) { 
            rDir <- seq(0, 2 * pi, length = dir + 1)[-(dir + 
                1)] 
            segments(0, 0, max(rpretty) * cos(rDir), max(rpretty) * 
                sin(rDir), col = grid.col, lwd = grid.lwd) 
        } 
        grid <- seq(0, 2 * pi, length = 360/4 + 1) 
        for (rad in rpretty) { 
            if (rad > 0) 
                lines(rad * cos(grid), rad * sin(grid), col = grid.col, 
                  lwd = grid.lwd) 
        } 
        if (rlabel.method != 0) { 
            if (rlabel.method == 1) 
                radLabels <- 1:length(rpretty) 
            if (rlabel.method == 2) 
                radLabels <- 2:length(rpretty) 
            if (rlabel.method == 3) 
                radLabels <- 1:(length(rpretty) - 1) 
            if (rlabel.method == 4) { 
                if (length(rpretty) > 2) 
                  radLabels <- 2:(length(rpretty) - 1) 
                else radLabels <- NULL 
            } 
            if (!is.null(radLabels)) { 
                xpos <- rpretty[radLabels] * cos(rlabel.axis) 
                ypos <- rpretty[radLabels] * sin(rlabel.axis) 
                text(xpos, ypos, rpretty[radLabels], cex = rlabel.cex, 
                  pos = rlabel.pos, col = rlabel.col) 
            } 
        } 
        if (!is.numeric(num.lab)) {
            t.lab <- ""
            if (pi2.lab & !is.character(text.lab)) 
                t.lab <- expression(0, pi/2, pi, 3 * pi/2) 
            if (!pi2.lab & is.character(text.lab)) 
                t.lab <- text.lab 
            labDir <- seq(0, 2 * pi, length = length(t.lab) + 
                1)[-(length(t.lab) + 1)] 
            labDir <- fit.rad(theta.zero + (!theta.clw) * labDir - 
                (theta.clw) * labDir) 
            text(lab.dist * cos(labDir), lab.dist * sin(labDir), 
                t.lab, cex = tlabel.cex, col = tlabel.col)