demo()の内容
demo()
demo(graphics)
##
##
## demo(graphics)
## ---- ~~~~~~~~
##
## > # Copyright (C) 1997-2009 The R Core Team
## >
## > require(datasets)
##
## > require(grDevices); require(graphics)
##
## > ## Here is some code which illustrates some of the differences between
## > ## R and S graphics capabilities. Note that colors are generally specified
## > ## by a character string name (taken from the X11 rgb.txt file) and that line
## > ## textures are given similarly. The parameter "bg" sets the background
## > ## parameter for the plot and there is also an "fg" parameter which sets
## > ## the foreground color.
## >
## >
## > x <- stats::rnorm(50)
##
## > opar <- par(bg = "white")
##
## > plot(x, ann = FALSE, type = "n")
##
## > abline(h = 0, col = gray(.90))
##
## > lines(x, col = "green4", lty = "dotted")
##
## > points(x, bg = "limegreen", pch = 21)
##
## > title(main = "Simple Use of Color In a Plot",
## + xlab = "Just a Whisper of a Label",
## + col.main = "blue", col.lab = gray(.8),
## + cex.main = 1.2, cex.lab = 1.0, font.main = 4, font.lab = 3)
##
## > ## A little color wheel. This code just plots equally spaced hues in
## > ## a pie chart. If you have a cheap SVGA monitor (like me) you will
## > ## probably find that numerically equispaced does not mean visually
## > ## equispaced. On my display at home, these colors tend to cluster at
## > ## the RGB primaries. On the other hand on the SGI Indy at work the
## > ## effect is near perfect.
## >
## > par(bg = "gray")
##
## > pie(rep(1,24), col = rainbow(24), radius = 0.9)
##
## > title(main = "A Sample Color Wheel", cex.main = 1.4, font.main = 3)
##
## > title(xlab = "(Use this as a test of monitor linearity)",
## + cex.lab = 0.8, font.lab = 3)
##
## > ## We have already confessed to having these. This is just showing off X11
## > ## color names (and the example (from the postscript manual) is pretty "cute".
## >
## > 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 Cream")
##
## > pie(pie.sales,
## + col = c("purple","violetred1","green3","cornsilk","cyan","white"))
##
## > title(main = "January Pie Sales", cex.main = 1.8, font.main = 1)
##
## > title(xlab = "(Don't try this at home kids)", cex.lab = 0.8, font.lab = 3)
##
## > ## Boxplots: I couldn't resist the capability for filling the "box".
## > ## The use of color seems like a useful addition, it focuses attention
## > ## on the central bulk of the data.
## >
## > par(bg="cornsilk")
##
## > n <- 10
##
## > g <- gl(n, 100, n*100)
##
## > x <- rnorm(n*100) + sqrt(as.numeric(g))
##
## > boxplot(split(x,g), col="lavender", notch=TRUE)
##
## > title(main="Notched Boxplots", xlab="Group", font.main=4, font.lab=1)
##
## > ## An example showing how to fill between curves.
## >
## > par(bg="white")
##
## > n <- 100
##
## > x <- c(0,cumsum(rnorm(n)))
##
## > y <- c(0,cumsum(rnorm(n)))
##
## > xx <- c(0:n, n:0)
##
## > yy <- c(x, rev(y))
##
## > plot(xx, yy, type="n", xlab="Time", ylab="Distance")
##
## > polygon(xx, yy, col="gray")
##
## > title("Distance Between Brownian Motions")
##
## > ## Colored plot margins, axis labels and titles. You do need to be
## > ## careful with these kinds of effects. It's easy to go completely
## > ## over the top and you can end up with your lunch all over the keyboard.
## > ## On the other hand, my market research clients love it.
## >
## > x <- c(0.00, 0.40, 0.86, 0.85, 0.69, 0.48, 0.54, 1.09, 1.11, 1.73, 2.05, 2.02)
##
## > par(bg="lightgray")
##
## > plot(x, type="n", axes=FALSE, ann=FALSE)
##
## > usr <- par("usr")
##
## > rect(usr[1], usr[3], usr[2], usr[4], col="cornsilk", border="black")
##
## > lines(x, col="blue")
##
## > points(x, pch=21, bg="lightcyan", cex=1.25)
##
## > axis(2, col.axis="blue", las=1)
##
## > axis(1, at=1:12, lab=month.abb, col.axis="blue")
##
## > box()
##
## > title(main= "The Level of Interest in R", font.main=4, col.main="red")
##
## > title(xlab= "1996", col.lab="red")
##
## > ## A filled histogram, showing how to change the font used for the
## > ## main title without changing the other annotation.
## >
## > par(bg="cornsilk")
##
## > x <- rnorm(1000)
##
## > hist(x, xlim=range(-4, 4, x), col="lavender", main="")
##
## > title(main="1000 Normal Random Variates", font.main=3)
##
## > ## A scatterplot matrix
## > ## The good old Iris data (yet again)
## >
## > pairs(iris[1:4], main="Edgar Anderson's Iris Data", font.main=4, pch=19)
##
## > pairs(iris[1:4], main="Edgar Anderson's Iris Data", pch=21,
## + bg = c("red", "green3", "blue")[unclass(iris$Species)])
##
## > ## Contour plotting
## > ## This produces a topographic map of one of Auckland's many volcanic "peaks".
## >
## > x <- 10*1:nrow(volcano)
##
## > y <- 10*1:ncol(volcano)
##
## > lev <- pretty(range(volcano), 10)
##
## > par(bg = "lightcyan")
##
## > pin <- par("pin")
##
## > xdelta <- diff(range(x))
##
## > ydelta <- diff(range(y))
##
## > xscale <- pin[1]/xdelta
##
## > yscale <- pin[2]/ydelta
##
## > scale <- min(xscale, yscale)
##
## > xadd <- 0.5*(pin[1]/scale - xdelta)
##
## > yadd <- 0.5*(pin[2]/scale - ydelta)
##
## > plot(numeric(0), numeric(0),
## + xlim = range(x)+c(-1,1)*xadd, ylim = range(y)+c(-1,1)*yadd,
## + type = "n", ann = FALSE)
##
## > usr <- par("usr")
##
## > rect(usr[1], usr[3], usr[2], usr[4], col="green3")
##
## > contour(x, y, volcano, levels = lev, col="yellow", lty="solid", add=TRUE)
##
## > box()
##
## > title("A Topographic Map of Maunga Whau", font= 4)
##
## > title(xlab = "Meters North", ylab = "Meters West", font= 3)
##
## > mtext("10 Meter Contour Spacing", side=3, line=0.35, outer=FALSE,
## + at = mean(par("usr")[1:2]), cex=0.7, font=3)
##
## > ## Conditioning plots
## >
## > par(bg="cornsilk")
##
## > coplot(lat ~ long | depth, data = quakes, pch = 21, bg = "green3")
##
## > par(opar)
demo(image)
##
##
## demo(image)
## ---- ~~~~~
##
## > # Copyright (C) 1997-2009 The R Core Team
## >
## > require(datasets)
##
## > require(grDevices); require(graphics)
##
## > x <- 10*(1:nrow(volcano)); x.at <- seq(100, 800, by=100)
##
## > y <- 10*(1:ncol(volcano)); y.at <- seq(100, 600, by=100)
##
## > # Using Terrain Colors
## >
## > image(x, y, volcano, col=terrain.colors(100),axes=FALSE)
##
## > contour(x, y, volcano, levels=seq(90, 200, by=5), add=TRUE, col="brown")
##
## > axis(1, at=x.at)
##
## > axis(2, at=y.at)
##
## > box()
##
## > title(main="Maunga Whau Volcano", sub = "col=terrain.colors(100)", font.main=4)
##
## > # Using Heat Colors
## >
## > image(x, y, volcano, col=heat.colors(100), axes=FALSE)
##
## > contour(x, y, volcano, levels=seq(90, 200, by=5), add=TRUE, col="brown")
##
## > axis(1, at=x.at)
##
## > axis(2, at=y.at)
##
## > box()
##
## > title(main="Maunga Whau Volcano", sub = "col=heat.colors(100)", font.main=4)
##
## > # Using Gray Scale
## >
## > image(x, y, volcano, col=gray(100:200/200), axes=FALSE)
##
## > contour(x, y, volcano, levels=seq(90, 200, by=5), add=TRUE, col="black")
##
## > axis(1, at=x.at)
##
## > axis(2, at=y.at)
##
## > box()
##
## > title(main="Maunga Whau Volcano \n col=gray(100:200/200)", font.main=4)
##
## > ## Filled Contours are even nicer sometimes :
## > example(filled.contour)
##
## flld.c> require(grDevices) # for colours
##
## flld.c> filled.contour(volcano, color = terrain.colors, asp = 1) # simple
##
## flld.c> x <- 10*1:nrow(volcano)
##
## flld.c> y <- 10*1:ncol(volcano)
##
## flld.c> filled.contour(x, y, volcano, color = terrain.colors,
## flld.c+ plot.title = title(main = "The Topography of Maunga Whau",
## flld.c+ xlab = "Meters North", ylab = "Meters West"),
## flld.c+ plot.axes = { axis(1, seq(100, 800, by = 100))
## flld.c+ axis(2, seq(100, 600, by = 100)) },
## flld.c+ key.title = title(main = "Height\n(meters)"),
## flld.c+ key.axes = axis(4, seq(90, 190, by = 10))) # maybe also asp = 1
##
## flld.c> mtext(paste("filled.contour(.) from", R.version.string),
## flld.c+ side = 1, line = 4, adj = 1, cex = .66)
##
## flld.c> # Annotating a filled contour plot
## flld.c> a <- expand.grid(1:20, 1:20)
##
## flld.c> b <- matrix(a[,1] + a[,2], 20)
##
## flld.c> filled.contour(x = 1:20, y = 1:20, z = b,
## flld.c+ plot.axes = { axis(1); axis(2); points(10, 10) })
##
## flld.c> ## Persian Rug Art:
## flld.c> x <- y <- seq(-4*pi, 4*pi, len = 27)
##
## flld.c> r <- sqrt(outer(x^2, y^2, "+"))
##
## flld.c> filled.contour(cos(r^2)*exp(-r/(2*pi)), axes = FALSE)
##
## flld.c> ## rather, the key *should* be labeled:
## flld.c> filled.contour(cos(r^2)*exp(-r/(2*pi)), frame.plot = FALSE,
## flld.c+ plot.axes = {})
demo(persp)
##
##
## demo(persp)
## ---- ~~~~~
##
## > ### Demos for persp() plots -- things not in example(persp)
## > ### -------------------------
## >
## > require(datasets)
##
## > require(grDevices); require(graphics)
##
## > ## (1) The Obligatory Mathematical surface.
## > ## Rotated sinc function.
## >
## > x <- seq(-10, 10, length.out = 50)
##
## > y <- x
##
## > rotsinc <- function(x,y)
## + {
## + sinc <- function(x) { y <- sin(x)/x ; y[is.na(y)] <- 1; y }
## + 10 * sinc( sqrt(x^2+y^2) )
## + }
##
## > sinc.exp <- expression(z == Sinc(sqrt(x^2 + y^2)))
##
## > z <- outer(x, y, rotsinc)
##
## > oldpar <- par(bg = "white")
##
## > persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue")
##
## > title(sub=".")## work around persp+plotmath bug
##
## > title(main = sinc.exp)
##
## > persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
## + ltheta = 120, shade = 0.75, ticktype = "detailed",
## + xlab = "X", ylab = "Y", zlab = "Z")
##
## > title(sub=".")## work around persp+plotmath bug
##
## > title(main = sinc.exp)
##
## > ## (2) Visualizing a simple DEM model
## >
## > z <- 2 * volcano # Exaggerate the relief
##
## > x <- 10 * (1:nrow(z)) # 10 meter spacing (S to N)
##
## > y <- 10 * (1:ncol(z)) # 10 meter spacing (E to W)
##
## > persp(x, y, z, theta = 120, phi = 15, scale = FALSE, axes = FALSE)
##
## > ## (3) Now something more complex
## > ## We border the surface, to make it more "slice like"
## > ## and color the top and sides of the surface differently.
## >
## > z0 <- min(z) - 20
##
## > z <- rbind(z0, cbind(z0, z, z0), z0)
##
## > x <- c(min(x) - 1e-10, x, max(x) + 1e-10)
##
## > y <- c(min(y) - 1e-10, y, max(y) + 1e-10)
##
## > fill <- matrix("green3", nrow = nrow(z)-1, ncol = ncol(z)-1)
##
## > fill[ , i2 <- c(1,ncol(fill))] <- "gray"
##
## > fill[i1 <- c(1,nrow(fill)) , ] <- "gray"
##
## > par(bg = "lightblue")
##
## > persp(x, y, z, theta = 120, phi = 15, col = fill, scale = FALSE, axes = FALSE)
##
## > title(main = "Maunga Whau\nOne of 50 Volcanoes in the Auckland Region.",
## + font.main = 4)
##
## > par(bg = "slategray")
##
## > persp(x, y, z, theta = 135, phi = 30, col = fill, scale = FALSE,
## + ltheta = -120, lphi = 15, shade = 0.65, axes = FALSE)
##
## > ## Don't draw the grid lines : border = NA
## > persp(x, y, z, theta = 135, phi = 30, col = "green3", scale = FALSE,
## + ltheta = -120, shade = 0.75, border = NA, box = FALSE)
##
## > ## `color gradient in the soil' :
## > fcol <- fill ; fcol[] <- terrain.colors(nrow(fcol))
##
## > persp(x, y, z, theta = 135, phi = 30, col = fcol, scale = FALSE,
## + ltheta = -120, shade = 0.3, border = NA, box = FALSE)
##
## > ## `image like' colors on top :
## > fcol <- fill
##
## > zi <- volcano[ -1,-1] + volcano[ -1,-61] +
## + volcano[-87,-1] + volcano[-87,-61] ## / 4
##
## > fcol[-i1,-i2] <-
## + terrain.colors(20)[cut(zi,
## + stats::quantile(zi, seq(0,1, length.out = 21)),
## + include.lowest = TRUE)]
##
## > persp(x, y, 2*z, theta = 110, phi = 40, col = fcol, scale = FALSE,
## + ltheta = -120, shade = 0.4, border = NA, box = FALSE)
##
## > ## reset par():
## > par(oldpar)
demo(plotmath)
##
##
## demo(plotmath)
## ---- ~~~~~~~~
##
## > # Copyright (C) 2002-2009 The R Core Team
## >
## > require(datasets)
##
## > require(grDevices); require(graphics)
##
## > ## --- "math annotation" in plots :
## >
## > ######
## > # create tables of mathematical annotation functionality
## > ######
## > make.table <- function(nr, nc) {
## + savepar <- par(mar=rep(0, 4), pty="s")
## + plot(c(0, nc*2 + 1), c(0, -(nr + 1)),
## + type="n", xlab="", ylab="", axes=FALSE)
## + savepar
## + }
##
## > get.r <- function(i, nr) {
## + i %% nr + 1
## + }
##
## > get.c <- function(i, nr) {
## + i %/% nr + 1
## + }
##
## > draw.title.cell <- function(title, i, nr) {
## + r <- get.r(i, nr)
## + c <- get.c(i, nr)
## + text(2*c - .5, -r, title)
## + rect((2*(c - 1) + .5), -(r - .5), (2*c + .5), -(r + .5))
## + }
##
## > draw.plotmath.cell <- function(expr, i, nr, string = NULL) {
## + r <- get.r(i, nr)
## + c <- get.c(i, nr)
## + if (is.null(string)) {
## + string <- deparse(expr)
## + string <- substr(string, 12, nchar(string) - 1)
## + }
## + text((2*(c - 1) + 1), -r, string, col="grey")
## + text((2*c), -r, expr, adj=c(.5,.5))
## + rect((2*(c - 1) + .5), -(r - .5), (2*c + .5), -(r + .5), border="grey")
## + }
##
## > nr <- 20
##
## > nc <- 2
##
## > oldpar <- make.table(nr, nc)
##
## > i <- 0
##
## > draw.title.cell("Arithmetic Operators", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x + y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x - y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x * y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x / y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %+-% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %/% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %*% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %.% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(-x), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(+x), i, nr); i <- i + 1
##
## > draw.title.cell("Sub/Superscripts", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x[i]), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x^2), i, nr); i <- i + 1
##
## > draw.title.cell("Juxtaposition", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x * y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(paste(x, y, z)), i, nr); i <- i + 1
##
## > draw.title.cell("Radicals", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(sqrt(x)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(sqrt(x, y)), i, nr); i <- i + 1
##
## > draw.title.cell("Lists", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(list(x, y, z)), i, nr); i <- i + 1
##
## > draw.title.cell("Relations", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x == y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x != y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x < y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x <= y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x > y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x >= y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %~~% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %=~% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %==% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %prop% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %~% y), i, nr); i <- i + 1
##
## > draw.title.cell("Typeface", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(plain(x)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(italic(x)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(bold(x)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(bolditalic(x)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(underline(x)), i, nr); i <- i + 1
##
## > # Need fewer, wider columns for ellipsis ...
## > nr <- 20
##
## > nc <- 2
##
## > make.table(nr, nc)
## $mar
## [1] 0 0 0 0
##
## $pty
## [1] "s"
##
##
## > i <- 0
##
## > draw.title.cell("Ellipsis", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(list(x[1], ..., x[n])), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x[1] + ... + x[n]), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(list(x[1], cdots, x[n])), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x[1] + ldots + x[n]), i, nr); i <- i + 1
##
## > draw.title.cell("Set Relations", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %subset% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %subseteq% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %supset% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %supseteq% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %notsubset% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %in% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %notin% y), i, nr); i <- i + 1
##
## > draw.title.cell("Accents", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(hat(x)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(tilde(x)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(ring(x)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(bar(xy)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(widehat(xy)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(widetilde(xy)), i, nr); i <- i + 1
##
## > draw.title.cell("Arrows", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %<->% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %->% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %<-% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %up% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %down% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %<=>% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %=>% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %<=% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %dblup% y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x %dbldown% y), i, nr); i <- i + 1
##
## > draw.title.cell("Symbolic Names", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(Alpha - Omega), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(alpha - omega), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(phi1 + sigma1), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(Upsilon1), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(infinity), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(32 * degree), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(60 * minute), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(30 * second), i, nr); i <- i + 1
##
## > # Need even fewer, wider columns for typeface and style ...
## > nr <- 20
##
## > nc <- 1
##
## > make.table(nr, nc)
## $mar
## [1] 0 0 0 0
##
## $pty
## [1] "s"
##
##
## > i <- 0
##
## > draw.title.cell("Style", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(displaystyle(x)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(textstyle(x)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(scriptstyle(x)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(scriptscriptstyle(x)), i, nr); i <- i + 1
##
## > draw.title.cell("Spacing", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x ~~ y), i, nr); i <- i + 1
##
## > # Need fewer, taller rows for fractions ...
## > # cheat a bit to save pages
## > par(new = TRUE)
##
## > nr <- 10
##
## > nc <- 1
##
## > make.table(nr, nc)
## $mar
## [1] 0 0 0 0
##
## $pty
## [1] "s"
##
##
## > i <- 4
##
## > draw.plotmath.cell(expression(x + phantom(0) + y), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x + over(1, phantom(0))), i, nr); i <- i + 1
##
## > draw.title.cell("Fractions", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(frac(x, y)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(over(x, y)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(atop(x, y)), i, nr); i <- i + 1
##
## > # Need fewer, taller rows and fewer, wider columns for big operators ...
## > nr <- 10
##
## > nc <- 1
##
## > make.table(nr, nc)
## $mar
## [1] 0 0 0 0
##
## $pty
## [1] "s"
##
##
## > i <- 0
##
## > draw.title.cell("Big Operators", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(sum(x[i], i=1, n)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(prod(plain(P)(X == x), x)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(integral(f(x) * dx, a, b)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(union(A[i], i==1, n)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(intersect(A[i], i==1, n)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(lim(f(x), x %->% 0)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(min(g(x), x >= 0)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(inf(S)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(sup(S)), i, nr); i <- i + 1
##
## > make.table(nr, nc)
## $mar
## [1] 0 0 0 0
##
## $pty
## [1] "s"
##
##
## > i <- 0
##
## > draw.title.cell("Grouping", i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression((x + y)*z), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x^y + z), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(x^(y + z)), i, nr); i <- i + 1
##
## > # have to do this one by hand
## > draw.plotmath.cell(expression(x^{y + z}), i, nr, string="x^{y + z}"); i <- i + 1
##
## > draw.plotmath.cell(expression(group("(", list(a, b), "]")), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(bgroup("(", atop(x, y), ")")), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(group(lceil, x, rceil)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(group(lfloor, x, rfloor)), i, nr); i <- i + 1
##
## > draw.plotmath.cell(expression(group("|", x, "|")), i, nr); i <- i + 1
##
## > par(oldpar)
demo(colors)
##
##
## demo(colors)
## ---- ~~~~~~
##
## > ### ----------- Show (almost) all named colors ---------------------
## >
## > ## 1) with traditional 'graphics' package:
## > showCols1 <- function(bg = "gray", cex = 0.75, srt = 30) {
## + m <- ceiling(sqrt(n <- length(cl <- colors())))
## + length(cl) <- m*m; cm <- matrix(cl, m)
## + ##
## + require("graphics")
## + op <- par(mar=rep(0,4), ann=FALSE, bg = bg); on.exit(par(op))
## + plot(1:m,1:m, type="n", axes=FALSE)
## + text(col(cm), rev(row(cm)), cm, col = cl, cex=cex, srt=srt)
## + }
##
## > showCols1()
##
## > ## 2) with 'grid' package:
## > showCols2 <- function(bg = "grey", cex = 0.75, rot = 30) {
## + m <- ceiling(sqrt(n <- length(cl <- colors())))
## + length(cl) <- m*m; cm <- matrix(cl, m)
## + ##
## + require("grid")
## + grid.newpage(); vp <- viewport(w = .92, h = .92)
## + grid.rect(gp=gpar(fill=bg))
## + grid.text(cm, x = col(cm)/m, y = rev(row(cm))/m, rot = rot,
## + vp=vp, gp=gpar(cex = cex, col = cm))
## + }
##
## > showCols2()
## Loading required package: grid
##
## > showCols2(bg = "gray33")
##
## > ###
## >
## > ##' @title Comparing Colors
## > ##' @param col
## > ##' @param nrow
## > ##' @param ncol
## > ##' @param txt.col
## > ##' @return the grid layout, invisibly
## > ##' @author Marius Hofert, originally
## > plotCol <- function(col, nrow=1, ncol=ceiling(length(col) / nrow),
## + txt.col="black") {
## + stopifnot(nrow >= 1, ncol >= 1)
## + if(length(col) > nrow*ncol)
## + warning("some colors will not be shown")
## + require(grid)
## + grid.newpage()
## + gl <- grid.layout(nrow, ncol)
## + pushViewport(viewport(layout=gl))
## + ic <- 1
## + for(i in 1:nrow) {
## + for(j in 1:ncol) {
## + pushViewport(viewport(layout.pos.row=i, layout.pos.col=j))
## + grid.rect(gp= gpar(fill=col[ic]))
## + grid.text(col[ic], gp=gpar(col=txt.col))
## + upViewport()
## + ic <- ic+1
## + }
## + }
## + upViewport()
## + invisible(gl)
## + }
##
## > ## A Chocolate Bar of colors:
## > plotCol(c("#CC8C3C", paste0("chocolate", 2:4),
## + paste0("darkorange", c("",1:2)), paste0("darkgoldenrod", 1:2),
## + "orange", "orange1", "sandybrown", "tan1", "tan2"),
## + nrow=2)
##
## > ##' Find close R colors() to a given color {original by Marius Hofert)
## > ##' using Euclidean norm in (HSV / RGB / ...) color space
## > nearRcolor <- function(rgb, cSpace = c("hsv", "rgb255", "Luv", "Lab"),
## + dist = switch(cSpace, "hsv" = 0.10, "rgb255" = 30,
## + "Luv" = 15, "Lab" = 12))
## + {
## + if(is.character(rgb)) rgb <- col2rgb(rgb)
## + stopifnot(length(rgb <- as.vector(rgb)) == 3)
## + Rcol <- col2rgb(.cc <- colors())
## + uniqC <- !duplicated(t(Rcol)) # gray9 == grey9 (etc)
## + Rcol <- Rcol[, uniqC] ; .cc <- .cc[uniqC]
## + cSpace <- match.arg(cSpace)
## + convRGB2 <- function(Rgb, to)
## + t(convertColor(t(Rgb), from="sRGB", to=to, scale.in=255))
## + ## the transformation, rgb{0..255} --> cSpace :
## + TransF <- switch(cSpace,
## + "rgb255" = identity,
## + "hsv" = rgb2hsv,
## + "Luv" = function(RGB) convRGB2(RGB, "Luv"),
## + "Lab" = function(RGB) convRGB2(RGB, "Lab"))
## + d <- sqrt(colSums((TransF(Rcol) - as.vector(TransF(rgb)))^2))
## + iS <- sort.list(d[near <- d <= dist])# sorted: closest first
## + setNames(.cc[near][iS], format(d[near][iS], digits=3))
## + }
##
## > nearRcolor(col2rgb("tan2"), "rgb")
## 0.0 21.1 25.8 29.5
## "tan2" "tan1" "sandybrown" "sienna1"
##
## > nearRcolor(col2rgb("tan2"), "hsv")
## 0.0000 0.0410 0.0618 0.0638 0.0667
## "tan2" "sienna2" "coral2" "tomato2" "tan1"
## 0.0766 0.0778 0.0900 0.0912 0.0918
## "coral" "sienna1" "sandybrown" "coral1" "tomato"
##
## > nearRcolor(col2rgb("tan2"), "Luv")
## 0.00 7.42 7.48 12.41 13.69
## "tan2" "tan1" "sandybrown" "orange3" "orange2"
##
## > nearRcolor(col2rgb("tan2"), "Lab")
## 0.00 5.56 8.08 11.31
## "tan2" "tan1" "sandybrown" "peru"
##
## > nearRcolor("#334455")
## 0.0867
## "darkslategray"
##
## > ## Now, consider choosing a color by looking in the
## > ## neighborhood of one you know :
## >
## > plotCol(nearRcolor("deepskyblue", "rgb", dist=50))
##
## > plotCol(nearRcolor("deepskyblue", dist=.1))
##
## > plotCol(nearRcolor("tomato", "rgb", dist= 50), nrow=3)
##
## > plotCol(nearRcolor("tomato", "hsv", dist=.12), nrow=3)
##
## > plotCol(nearRcolor("tomato", "Luv", dist= 25), nrow=3)
##
## > plotCol(nearRcolor("tomato", "Lab", dist= 18), nrow=3)
demo(hclColors)
##
##
## demo(hclColors)
## ---- ~~~~~~~~~
##
## > ### ------ hcl() explorations
## >
## > hcl.wheel <-
## + function(chroma = 35, lums = 0:100, hues = 1:360, asp = 1,
## + p.cex = 0.6, do.label = FALSE, rev.lum = FALSE,
## + fixup = TRUE)
## + {
## + ## Purpose: show chroma "sections" of hcl() color space; see ?hcl
## + ## ----------------------------------------------------------------------
## + ## Arguments: chroma: can be vector -> multiple plots are done,
## + ## lums, hues, fixup : all corresponding to hcl()'s args
## + ## rev.lum: logical indicating if luminance
## + ## should go from outer to inner
## + ## ----------------------------------------------------------------------
## + ## Author: Martin Maechler, Date: 24 Jun 2005
## +
## + require("graphics")
## + stopifnot(is.numeric(lums), lums >= 0, lums <= 100,
## + is.numeric(hues), hues >= 0, hues <= 360,
## + is.numeric(chroma), chroma >= 0, (nch <- length(chroma)) >= 1)
## + if(is.unsorted(hues)) hues <- sort(hues)
## + if(nch > 1) {
## + op <- par(mfrow= n2mfrow(nch), mar = c(0,0,0,0), xpd = TRUE)
## + on.exit(par(op))
## + }
## + for(i.c in 1:nch) {
## + plot(-1:1,-1:1, type="n", axes = FALSE, xlab="",ylab="", asp = asp)
## + ## main = sprintf("hcl(h = <angle>, c = %g)", chroma[i.c]),
## + text(0.4, 0.99, paste("chroma =", format(chroma[i.c])),
## + adj = 0, font = 4)
## + l.s <- (if(rev.lum) rev(lums) else lums) / max(lums) # <= 1
## + for(ang in hues) { # could do all this using outer() instead of for()...
## + a. <- ang * pi/180
## + z.a <- exp(1i * a.)
## + cols <- hcl(ang, c = chroma[i.c], l = lums, fixup = fixup)
## + points(l.s * z.a, pch = 16, col = cols, cex = p.cex)
## + ##if(do."text") : draw the 0,45,90,... angle "lines"
## + if(do.label)
## + text(z.a*1.05, labels = ang, col = cols[length(cols)/2],
## + srt = ang)
## + }
## + if(!fixup) ## show the outline
## + lines(exp(1i * hues * pi/180))
## + }
## + invisible()
## + }
##
## > ## and now a few interesting calls :
## >
## > hcl.wheel() # and watch it redraw when you fiddle with the graphic window
##
## > hcl.wheel(rev.lum= TRUE) # ditto
##
## > hcl.wheel(do.lab = TRUE) # ditto
##
## > ## Now watch:
## > hcl.wheel(ch = c(25,35,45,55))
##
## > hcl.wheel(ch = seq(10, 90, by = 10), p.cex = 0.4)
##
## > hcl.wheel(ch = seq(10, 90, by = 10), p.cex = 0.3, fixup = FALSE)
##
## > hcl.wheel(ch = seq(10, 90, by = 10), p.cex = 0.3, rev.lum = TRUE)
##
## > if(dev.interactive()) # new "graphics window" -- to compare with previous :
## + dev.new()
##
## > hcl.wheel(ch = seq(10, 90, by = 10), p.cex = 0.3, rev.lum = TRUE, fixup=FALSE)
demo(nlm)
##
##
## demo(nlm)
## ---- ~~~
##
## > # Copyright (C) 1997-2009 The R Core Team
## >
## > ### Helical Valley Function
## > ### Page 362 Dennis + Schnabel
## >
## > require(stats); require(graphics)
##
## > theta <- function(x1,x2) (atan(x2/x1) + (if(x1 <= 0) pi else 0))/ (2*pi)
##
## > ## but this is easier :
## > theta <- function(x1,x2) atan2(x2, x1)/(2*pi)
##
## > f <- function(x) {
## + f1 <- 10*(x[3] - 10*theta(x[1],x[2]))
## + f2 <- 10*(sqrt(x[1]^2+x[2]^2)-1)
## + f3 <- x[3]
## + return(f1^2+f2^2+f3^2)
## + }
##
## > ## explore surface {at x3 = 0}
## > x <- seq(-1, 2, length.out=50)
##
## > y <- seq(-1, 1, length.out=50)
##
## > z <- apply(as.matrix(expand.grid(x, y)), 1, function(x) f(c(x, 0)))
##
## > contour(x, y, matrix(log10(z), 50, 50))
##
## > str(nlm.f <- nlm(f, c(-1,0,0), hessian = TRUE))
## List of 6
## $ minimum : num 1.24e-14
## $ estimate : num [1:3] 1.00 3.07e-09 -6.06e-09
## $ gradient : num [1:3] -3.76e-07 3.49e-06 -2.20e-06
## $ hessian : num [1:3, 1:3] 2.00e+02 -4.07e-02 9.77e-07 -4.07e-02 5.07e+02 ...
## $ code : int 2
## $ iterations: int 27
##
## > points(rbind(nlm.f$estim[1:2]), col = "red", pch = 20)
##
## > ### the Rosenbrock banana valley function
## >
## > fR <- function(x)
## + {
## + x1 <- x[1]; x2 <- x[2]
## + 100*(x2 - x1*x1)^2 + (1-x1)^2
## + }
##
## > ## explore surface
## > fx <- function(x)
## + { ## `vectorized' version of fR()
## + x1 <- x[,1]; x2 <- x[,2]
## + 100*(x2 - x1*x1)^2 + (1-x1)^2
## + }
##
## > x <- seq(-2, 2, length.out=100)
##
## > y <- seq(-0.5, 1.5, length.out=100)
##
## > z <- fx(expand.grid(x, y))
##
## > op <- par(mfrow = c(2,1), mar = 0.1 + c(3,3,0,0))
##
## > contour(x, y, matrix(log10(z), length(x)))
##
## > str(nlm.f2 <- nlm(fR, c(-1.2, 1), hessian = TRUE))
## List of 6
## $ minimum : num 3.97e-12
## $ estimate : num [1:2] 1 1
## $ gradient : num [1:2] -6.54e-07 3.34e-07
## $ hessian : num [1:2, 1:2] 802 -400 -400 200
## $ code : int 1
## $ iterations: int 23
##
## > points(rbind(nlm.f2$estim[1:2]), col = "red", pch = 20)
##
## > ## Zoom in :
## > rect(0.9, 0.9, 1.1, 1.1, border = "orange", lwd = 2)
##
## > x <- y <- seq(0.9, 1.1, length.out=100)
##
## > z <- fx(expand.grid(x, y))
##
## > contour(x, y, matrix(log10(z), length(x)))
##
## > mtext("zoomed in");box(col = "orange")
##
## > points(rbind(nlm.f2$estim[1:2]), col = "red", pch = 20)
##
## > par(op)
##
## > fg <- function(x)
## + {
## + gr <- function(x1, x2) {
## + c(-400*x1*(x2 - x1*x1)-2*(1-x1), 200*(x2 - x1*x1))
## + }
## + x1 <- x[1]; x2 <- x[2]
## + res<- 100*(x2 - x1*x1)^2 + (1-x1)^2
## + attr(res, "gradient") <- gr(x1, x2)
## + return(res)
## + }
##
## > nlm(fg, c(-1.2, 1), hessian = TRUE)
## $minimum
## [1] 1.182096e-20
##
## $estimate
## [1] 1 1
##
## $gradient
## [1] 2.583521e-09 -1.201128e-09
##
## $hessian
## [,1] [,2]
## [1,] 802.24 -400.02
## [2,] -400.02 200.00
##
## $code
## [1] 1
##
## $iterations
## [1] 24
##
##
## > ## or use deriv to find the derivatives
## >
## > fd <- deriv(~ 100*(x2 - x1*x1)^2 + (1-x1)^2, c("x1", "x2"))
##
## > fdd <- function(x1, x2) {}
##
## > body(fdd) <- fd
##
## > nlm(function(x) fdd(x[1], x[2]), c(-1.2,1), hessian = TRUE)
## $minimum
## [1] 1.182096e-20
##
## $estimate
## [1] 1 1
##
## $gradient
## [1] 2.583521e-09 -1.201128e-09
##
## $hessian
## [,1] [,2]
## [1,] 802.24 -400.02
## [2,] -400.02 200.00
##
## $code
## [1] 1
##
## $iterations
## [1] 24
##
##
## > fgh <- function(x)
## + {
## + gr <- function(x1, x2)
## + c(-400*x1*(x2 - x1*x1) - 2*(1-x1), 200*(x2 - x1*x1))
## + h <- function(x1, x2) {
## + a11 <- 2 - 400*x2 + 1200*x1*x1
## + a21 <- -400*x1
## + matrix(c(a11, a21, a21, 200), 2, 2)
## + }
## + x1 <- x[1]; x2 <- x[2]
## + res<- 100*(x2 - x1*x1)^2 + (1-x1)^2
## + attr(res, "gradient") <- gr(x1, x2)
## + attr(res, "hessian") <- h(x1, x2)
## + return(res)
## + }
##
## > nlm(fgh, c(-1.2,1), hessian = TRUE)
## $minimum
## [1] 2.829175
##
## $estimate
## [1] -0.6786981 0.4711891
##
## $gradient
## [1] -0.4911201 2.1115987
##
## $hessian
## [,1] [,2]
## [1,] 366.1188 271.4593
## [2,] 271.4593 200.0000
##
## $code
## [1] 4
##
## $iterations
## [1] 100
demo(glm.vr)
##
##
## demo(glm.vr)
## ---- ~~~~~~
##
## > # Copyright (C) 1997-2009 The R Core Team
## >
## > #### -*- R -*-
## > require(stats)
##
## > Fr <- c(68,42,42,30, 37,52,24,43,
## + 66,50,33,23, 47,55,23,47,
## + 63,53,29,27, 57,49,19,29)
##
## > Temp <- gl(2, 2, 24, labels = c("Low", "High"))
##
## > Soft <- gl(3, 8, 24, labels = c("Hard","Medium","Soft"))
##
## > M.user <- gl(2, 4, 24, labels = c("N", "Y"))
##
## > Brand <- gl(2, 1, 24, labels = c("X", "M"))
##
## > detg <- data.frame(Fr,Temp, Soft,M.user, Brand)
##
## > detg.m0 <- glm(Fr ~ M.user*Temp*Soft + Brand, family = poisson, data = detg)
##
## > summary(detg.m0)
##
## Call:
## glm(formula = Fr ~ M.user * Temp * Soft + Brand, family = poisson,
## data = detg)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.20876 -0.99190 -0.00126 0.93542 1.97601
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.01524 0.10034 40.018 < 2e-16 ***
## M.userY -0.21184 0.14257 -1.486 0.13731
## TempHigh -0.42381 0.15159 -2.796 0.00518 **
## SoftMedium 0.05311 0.13308 0.399 0.68984
## SoftSoft 0.05311 0.13308 0.399 0.68984
## BrandM -0.01587 0.06300 -0.252 0.80106
## M.userY:TempHigh 0.13987 0.22168 0.631 0.52806
## M.userY:SoftMedium 0.08323 0.19685 0.423 0.67245
## M.userY:SoftSoft 0.12169 0.19591 0.621 0.53449
## TempHigh:SoftMedium -0.30442 0.22239 -1.369 0.17104
## TempHigh:SoftSoft -0.30442 0.22239 -1.369 0.17104
## M.userY:TempHigh:SoftMedium 0.21189 0.31577 0.671 0.50220
## M.userY:TempHigh:SoftSoft -0.20387 0.32540 -0.627 0.53098
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 118.627 on 23 degrees of freedom
## Residual deviance: 32.826 on 11 degrees of freedom
## AIC: 191.24
##
## Number of Fisher Scoring iterations: 4
##
##
## > detg.mod <- glm(terms(Fr ~ M.user*Temp*Soft + Brand*M.user*Temp,
## + keep.order = TRUE),
## + family = poisson, data = detg)
##
## > summary(detg.mod)
##
## Call:
## glm(formula = terms(Fr ~ M.user * Temp * Soft + Brand * M.user *
## Temp, keep.order = TRUE), family = poisson, data = detg)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.91365 -0.35585 0.00253 0.33027 0.92146
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.14887 0.10603 39.128 < 2e-16 ***
## M.userY -0.40521 0.16188 -2.503 0.01231 *
## TempHigh -0.44275 0.17121 -2.586 0.00971 **
## M.userY:TempHigh -0.12692 0.26257 -0.483 0.62883
## SoftMedium 0.05311 0.13308 0.399 0.68984
## SoftSoft 0.05311 0.13308 0.399 0.68984
## M.userY:SoftMedium 0.08323 0.19685 0.423 0.67245
## M.userY:SoftSoft 0.12169 0.19591 0.621 0.53449
## TempHigh:SoftMedium -0.30442 0.22239 -1.369 0.17104
## TempHigh:SoftSoft -0.30442 0.22239 -1.369 0.17104
## M.userY:TempHigh:SoftMedium 0.21189 0.31577 0.671 0.50220
## M.userY:TempHigh:SoftSoft -0.20387 0.32540 -0.627 0.53098
## BrandM -0.30647 0.10942 -2.801 0.00510 **
## M.userY:BrandM 0.40757 0.15961 2.554 0.01066 *
## TempHigh:BrandM 0.04411 0.18463 0.239 0.81119
## M.userY:TempHigh:BrandM 0.44427 0.26673 1.666 0.09579 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 118.627 on 23 degrees of freedom
## Residual deviance: 5.656 on 8 degrees of freedom
## AIC: 170.07
##
## Number of Fisher Scoring iterations: 4
##
##
## > summary(detg.mod, correlation = TRUE, symbolic.cor = TRUE)
##
## Call:
## glm(formula = terms(Fr ~ M.user * Temp * Soft + Brand * M.user *
## Temp, keep.order = TRUE), family = poisson, data = detg)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.91365 -0.35585 0.00253 0.33027 0.92146
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.14887 0.10603 39.128 < 2e-16 ***
## M.userY -0.40521 0.16188 -2.503 0.01231 *
## TempHigh -0.44275 0.17121 -2.586 0.00971 **
## M.userY:TempHigh -0.12692 0.26257 -0.483 0.62883
## SoftMedium 0.05311 0.13308 0.399 0.68984
## SoftSoft 0.05311 0.13308 0.399 0.68984
## M.userY:SoftMedium 0.08323 0.19685 0.423 0.67245
## M.userY:SoftSoft 0.12169 0.19591 0.621 0.53449
## TempHigh:SoftMedium -0.30442 0.22239 -1.369 0.17104
## TempHigh:SoftSoft -0.30442 0.22239 -1.369 0.17104
## M.userY:TempHigh:SoftMedium 0.21189 0.31577 0.671 0.50220
## M.userY:TempHigh:SoftSoft -0.20387 0.32540 -0.627 0.53098
## BrandM -0.30647 0.10942 -2.801 0.00510 **
## M.userY:BrandM 0.40757 0.15961 2.554 0.01066 *
## TempHigh:BrandM 0.04411 0.18463 0.239 0.81119
## M.userY:TempHigh:BrandM 0.44427 0.26673 1.666 0.09579 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 118.627 on 23 degrees of freedom
## Residual deviance: 5.656 on 8 degrees of freedom
## AIC: 170.07
##
## Number of Fisher Scoring iterations: 4
##
## Correlation of Coefficients:
##
## (Intercept) 1
## M.userY , 1
## TempHigh , . 1
## M.userY:TempHigh . , , 1
## SoftMedium , . . 1
## SoftSoft , . . . 1
## M.userY:SoftMedium . , . , . 1
## M.userY:SoftSoft . , . . , . 1
## TempHigh:SoftMedium . , . . . . 1
## TempHigh:SoftSoft . , . . . . . 1
## M.userY:TempHigh:SoftMedium . . . . , . , . 1
## M.userY:TempHigh:SoftSoft . . . . . , . , . 1
## BrandM . 1
## M.userY:BrandM . , 1
## TempHigh:BrandM . . . . 1
## M.userY:TempHigh:BrandM . . . . , 1
## attr(,"legend")
## [1] 0 ' ' 0.3 '.' 0.6 ',' 0.8 '+' 0.9 '*' 0.95 'B' 1
##
##
## > anova(detg.m0, detg.mod)
## Analysis of Deviance Table
##
## Model 1: Fr ~ M.user * Temp * Soft + Brand
## Model 2: Fr ~ M.user * Temp * Soft + Brand * M.user * Temp
## Resid. Df Resid. Dev Df Deviance
## 1 11 32.826
## 2 8 5.656 3 27.17
demo(lm.glm)
##
##
## demo(lm.glm)
## ---- ~~~~~~
##
## > ### Examples from: "An Introduction to Statistical Modelling"
## > ### By Annette Dobson
## > ###
## > ### == with some additions ==
## >
## > # Copyright (C) 1997-2008 The R Core Team
## >
## > require(stats); require(graphics)
##
## > ## Plant Weight Data (Page 9)
## > ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
##
## > trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
##
## > group <- gl(2,10, labels=c("Ctl","Trt"))
##
## > weight <- c(ctl,trt)
##
## > anova (lm(weight~group))
## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## group 1 0.6882 0.68820 1.4191 0.249
## Residuals 18 8.7293 0.48496
##
## > summary(lm(weight~group -1))
##
## Call:
## lm(formula = weight ~ group - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0710 -0.4938 0.0685 0.2462 1.3690
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## groupCtl 5.0320 0.2202 22.85 9.55e-15 ***
## groupTrt 4.6610 0.2202 21.16 3.62e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6964 on 18 degrees of freedom
## Multiple R-squared: 0.9818, Adjusted R-squared: 0.9798
## F-statistic: 485.1 on 2 and 18 DF, p-value: < 2.2e-16
##
##
## > ## Birth Weight Data (Page 14)
## > age <- c(40, 38, 40, 35, 36, 37, 41, 40, 37, 38, 40, 38,
## + 40, 36, 40, 38, 42, 39, 40, 37, 36, 38, 39, 40)
##
## > birthw <- c(2968, 2795, 3163, 2925, 2625, 2847, 3292, 3473, 2628, 3176,
## + 3421, 2975, 3317, 2729, 2935, 2754, 3210, 2817, 3126, 2539,
## + 2412, 2991, 2875, 3231)
##
## > sex <- gl(2,12, labels=c("M","F"))
##
## > plot(age, birthw, col=as.numeric(sex), pch=3*as.numeric(sex),
## + main="Dobson's Birth Weight Data")
##
## > lines(lowess(age[sex=='M'], birthw[sex=='M']), col=1)
##
## > lines(lowess(age[sex=='F'], birthw[sex=='F']), col=2)
##
## > legend("topleft", levels(sex), col=1:2, pch=3*(1:2), lty=1, bty="n")
##
## > summary(l1 <- lm(birthw ~ sex + age), correlation=TRUE)
##
## Call:
## lm(formula = birthw ~ sex + age)
##
## Residuals:
## Min 1Q Median 3Q Max
## -257.49 -125.28 -58.44 169.00 303.98
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1610.28 786.08 -2.049 0.0532 .
## sexF -163.04 72.81 -2.239 0.0361 *
## age 120.89 20.46 5.908 7.28e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 177.1 on 21 degrees of freedom
## Multiple R-squared: 0.64, Adjusted R-squared: 0.6057
## F-statistic: 18.67 on 2 and 21 DF, p-value: 2.194e-05
##
## Correlation of Coefficients:
## (Intercept) sexF
## sexF 0.07
## age -1.00 -0.12
##
##
## > summary(l0 <- lm(birthw ~ sex + age -1), correlation=TRUE)
##
## Call:
## lm(formula = birthw ~ sex + age - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -257.49 -125.28 -58.44 169.00 303.98
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## sexM -1610.28 786.08 -2.049 0.0532 .
## sexF -1773.32 794.59 -2.232 0.0367 *
## age 120.89 20.46 5.908 7.28e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 177.1 on 21 degrees of freedom
## Multiple R-squared: 0.9969, Adjusted R-squared: 0.9965
## F-statistic: 2258 on 3 and 21 DF, p-value: < 2.2e-16
##
## Correlation of Coefficients:
## sexM sexF
## sexF 1.00
## age -1.00 -1.00
##
##
## > anova(l1,l0)
## Analysis of Variance Table
##
## Model 1: birthw ~ sex + age
## Model 2: birthw ~ sex + age - 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 21 658771
## 2 21 658771 0 1.5134e-09
##
## > summary(li <- lm(birthw ~ sex + sex:age -1), correlation=TRUE)
##
## Call:
## lm(formula = birthw ~ sex + sex:age - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -246.69 -138.11 -39.13 176.57 274.28
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## sexM -1268.67 1114.64 -1.138 0.268492
## sexF -2141.67 1163.60 -1.841 0.080574 .
## sexM:age 111.98 29.05 3.855 0.000986 ***
## sexF:age 130.40 30.00 4.347 0.000313 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 180.6 on 20 degrees of freedom
## Multiple R-squared: 0.9969, Adjusted R-squared: 0.9963
## F-statistic: 1629 on 4 and 20 DF, p-value: < 2.2e-16
##
## Correlation of Coefficients:
## sexM sexF sexM:age
## sexF 0.00
## sexM:age -1.00 0.00
## sexF:age 0.00 -1.00 0.00
##
##
## > anova(li,l0)
## Analysis of Variance Table
##
## Model 1: birthw ~ sex + sex:age - 1
## Model 2: birthw ~ sex + age - 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 20 652425
## 2 21 658771 -1 -6346.2 0.1945 0.6639
##
## > summary(zi <- glm(birthw ~ sex + age, family=gaussian()))
##
## Call:
## glm(formula = birthw ~ sex + age, family = gaussian())
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -257.49 -125.28 -58.44 169.00 303.98
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1610.28 786.08 -2.049 0.0532 .
## sexF -163.04 72.81 -2.239 0.0361 *
## age 120.89 20.46 5.908 7.28e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 31370.04)
##
## Null deviance: 1829873 on 23 degrees of freedom
## Residual deviance: 658771 on 21 degrees of freedom
## AIC: 321.39
##
## Number of Fisher Scoring iterations: 2
##
##
## > summary(z0 <- glm(birthw ~ sex + age - 1, family=gaussian()))
##
## Call:
## glm(formula = birthw ~ sex + age - 1, family = gaussian())
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -257.49 -125.28 -58.44 169.00 303.98
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## sexM -1610.28 786.08 -2.049 0.0532 .
## sexF -1773.32 794.59 -2.232 0.0367 *
## age 120.89 20.46 5.908 7.28e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 31370.04)
##
## Null deviance: 213198964 on 24 degrees of freedom
## Residual deviance: 658771 on 21 degrees of freedom
## AIC: 321.39
##
## Number of Fisher Scoring iterations: 2
##
##
## > anova(zi, z0)
## Analysis of Deviance Table
##
## Model 1: birthw ~ sex + age
## Model 2: birthw ~ sex + age - 1
## Resid. Df Resid. Dev Df Deviance
## 1 21 658771
## 2 21 658771 0 -1.1642e-10
##
## > summary(z.o4 <- update(z0, subset = -4))
##
## Call:
## glm(formula = birthw ~ sex + age - 1, family = gaussian(), subset = -4)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -253.86 -129.46 -53.46 165.04 251.14
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## sexM -2318.03 801.57 -2.892 0.00902 **
## sexF -2455.44 803.79 -3.055 0.00625 **
## age 138.50 20.71 6.688 1.65e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 26925.39)
##
## Null deviance: 204643339 on 23 degrees of freedom
## Residual deviance: 538508 on 20 degrees of freedom
## AIC: 304.68
##
## Number of Fisher Scoring iterations: 2
##
##
## > summary(zz <- update(z0, birthw ~ sex+age-1 + sex:age))
##
## Call:
## glm(formula = birthw ~ sex + age + sex:age - 1, family = gaussian())
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -246.69 -138.11 -39.13 176.57 274.28
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## sexM -1268.67 1114.64 -1.138 0.268492
## sexF -2141.67 1163.60 -1.841 0.080574 .
## age 111.98 29.05 3.855 0.000986 ***
## sexF:age 18.42 41.76 0.441 0.663893
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 32621.23)
##
## Null deviance: 213198964 on 24 degrees of freedom
## Residual deviance: 652425 on 20 degrees of freedom
## AIC: 323.16
##
## Number of Fisher Scoring iterations: 2
##
##
## > anova(z0,zz)
## Analysis of Deviance Table
##
## Model 1: birthw ~ sex + age - 1
## Model 2: birthw ~ sex + age + sex:age - 1
## Resid. Df Resid. Dev Df Deviance
## 1 21 658771
## 2 20 652425 1 6346.2
##
## > ## Poisson Regression Data (Page 42)
## > x <- c(-1,-1,0,0,0,0,1,1,1)
##
## > y <- c(2,3,6,7,8,9,10,12,15)
##
## > summary(glm(y~x, family=poisson(link="identity")))
##
## Call:
## glm(formula = y ~ x, family = poisson(link = "identity"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.7019 -0.3377 -0.1105 0.2958 0.7184
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.4516 0.8841 8.428 < 2e-16 ***
## x 4.9353 1.0892 4.531 5.86e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 18.4206 on 8 degrees of freedom
## Residual deviance: 1.8947 on 7 degrees of freedom
## AIC: 40.008
##
## Number of Fisher Scoring iterations: 3
##
##
## > ## Calorie Data (Page 45)
## > calorie <- data.frame(
## + carb = c(33,40,37,27,30,43,34,48,30,38,
## + 50,51,30,36,41,42,46,24,35,37),
## + age = c(33,47,49,35,46,52,62,23,32,42,
## + 31,61,63,40,50,64,56,61,48,28),
## + wgt = c(100, 92,135,144,140,101, 95,101, 98,105,
## + 108, 85,130,127,109,107,117,100,118,102),
## + prot = c(14,15,18,12,15,15,14,17,15,14,
## + 17,19,19,20,15,16,18,13,18,14))
##
## > summary(lmcal <- lm(carb~age+wgt+prot, data= calorie))
##
## Call:
## lm(formula = carb ~ age + wgt + prot, data = calorie)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.3424 -4.8203 0.9897 3.8553 7.9087
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.96006 13.07128 2.828 0.01213 *
## age -0.11368 0.10933 -1.040 0.31389
## wgt -0.22802 0.08329 -2.738 0.01460 *
## prot 1.95771 0.63489 3.084 0.00712 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.956 on 16 degrees of freedom
## Multiple R-squared: 0.4805, Adjusted R-squared: 0.3831
## F-statistic: 4.934 on 3 and 16 DF, p-value: 0.01297
##
##
## > ## Extended Plant Data (Page 59)
## > ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
##
## > trtA <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
##
## > trtB <- c(6.31,5.12,5.54,5.50,5.37,5.29,4.92,6.15,5.80,5.26)
##
## > group <- gl(3, length(ctl), labels=c("Ctl","A","B"))
##
## > weight <- c(ctl,trtA,trtB)
##
## > anova(lmwg <- lm(weight~group))
## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## group 2 3.7663 1.8832 4.8461 0.01591 *
## Residuals 27 10.4921 0.3886
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## > summary(lmwg)
##
## Call:
## lm(formula = weight ~ group)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0710 -0.4180 -0.0060 0.2627 1.3690
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.0320 0.1971 25.527 <2e-16 ***
## groupA -0.3710 0.2788 -1.331 0.1944
## groupB 0.4940 0.2788 1.772 0.0877 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6234 on 27 degrees of freedom
## Multiple R-squared: 0.2641, Adjusted R-squared: 0.2096
## F-statistic: 4.846 on 2 and 27 DF, p-value: 0.01591
##
##
## > coef(lmwg)
## (Intercept) groupA groupB
## 5.032 -0.371 0.494
##
## > coef(summary(lmwg))#- incl. std.err, t- and P- values.
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.032 0.1971284 25.526514 1.936575e-20
## groupA -0.371 0.2787816 -1.330791 1.943879e-01
## groupB 0.494 0.2787816 1.771996 8.768168e-02
##
## > ## Fictitious Anova Data (Page 64)
## > y <- c(6.8,6.6,5.3,6.1,7.5,7.4,7.2,6.5,7.8,9.1,8.8,9.1)
##
## > a <- gl(3,4)
##
## > b <- gl(2,2, length(a))
##
## > anova(z <- lm(y~a*b))
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## a 2 12.7400 6.3700 25.8243 0.001127 **
## b 1 0.4033 0.4033 1.6351 0.248225
## a:b 2 1.2067 0.6033 2.4459 0.167164
## Residuals 6 1.4800 0.2467
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## > ## Achievement Scores (Page 70)
## > y <- c(6,4,5,3,4,3,6, 8,9,7,9,8,5,7, 6,7,7,7,8,5,7)
##
## > x <- c(3,1,3,1,2,1,4, 4,5,5,4,3,1,2, 3,2,2,3,4,1,4)
##
## > m <- gl(3,7)
##
## > anova(z <- lm(y~x+m))
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x 1 36.575 36.575 60.355 5.428e-07 ***
## m 2 16.932 8.466 13.970 0.0002579 ***
## Residuals 17 10.302 0.606
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## > ## Beetle Data (Page 78)
## > dose <- c(1.6907, 1.7242, 1.7552, 1.7842, 1.8113, 1.8369, 1.861, 1.8839)
##
## > x <- c( 6, 13, 18, 28, 52, 53, 61, 60)
##
## > n <- c(59, 60, 62, 56, 63, 59, 62, 60)
##
## > dead <- cbind(x, n-x)
##
## > summary( glm(dead ~ dose, family=binomial(link=logit)))
##
## Call:
## glm(formula = dead ~ dose, family = binomial(link = logit))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5941 -0.3944 0.8329 1.2592 1.5940
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -60.717 5.181 -11.72 <2e-16 ***
## dose 34.270 2.912 11.77 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 284.202 on 7 degrees of freedom
## Residual deviance: 11.232 on 6 degrees of freedom
## AIC: 41.43
##
## Number of Fisher Scoring iterations: 4
##
##
## > summary( glm(dead ~ dose, family=binomial(link=probit)))
##
## Call:
## glm(formula = dead ~ dose, family = binomial(link = probit))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5714 -0.4703 0.7501 1.0632 1.3449
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -34.935 2.648 -13.19 <2e-16 ***
## dose 19.728 1.487 13.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 284.20 on 7 degrees of freedom
## Residual deviance: 10.12 on 6 degrees of freedom
## AIC: 40.318
##
## Number of Fisher Scoring iterations: 4
##
##
## > summary(z <- glm(dead ~ dose, family=binomial(link=cloglog)))
##
## Call:
## glm(formula = dead ~ dose, family = binomial(link = cloglog))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.80329 -0.55135 0.03089 0.38315 1.28883
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -39.572 3.240 -12.21 <2e-16 ***
## dose 22.041 1.799 12.25 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 284.2024 on 7 degrees of freedom
## Residual deviance: 3.4464 on 6 degrees of freedom
## AIC: 33.644
##
## Number of Fisher Scoring iterations: 4
##
##
## > anova(z, update(z, dead ~ dose -1))
## Analysis of Deviance Table
##
## Model 1: dead ~ dose
## Model 2: dead ~ dose - 1
## Resid. Df Resid. Dev Df Deviance
## 1 6 3.446
## 2 7 285.222 -1 -281.78
##
## > ## Anther Data (Page 84)
## > ## Note that the proportions below are not exactly
## > ## in accord with the sample sizes quoted below.
## > ## In particular, the value 0.555 does not seem sensible.
## > ## [MM: huh? round(round(n*p)/n, 3) looks almost exactly like "p" !]
## > n <- c(102, 99, 108, 76, 81, 90)
##
## > p <- c(0.539,0.525,0.528,0.724,0.617,0.555)
##
## > x <- round(n*p)
##
## > ## x <- n*p
## > y <- cbind(x,n-x)
##
## > f <- rep(c(40,150,350),2)
##
## > (g <- gl(2,3))
## [1] 1 1 1 2 2 2
## Levels: 1 2
##
## > summary(glm(y ~ g*f, family=binomial(link="logit")))
##
## Call:
## glm(formula = y ~ g * f, family = binomial(link = "logit"))
##
## Deviance Residuals:
## 1 2 3 4 5 6
## 0.08269 -0.12998 0.04414 0.42320 -0.60082 0.19522
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.1456719 0.1975451 0.737 0.4609
## g2 0.7963143 0.3125046 2.548 0.0108 *
## f -0.0001227 0.0008782 -0.140 0.8889
## g2:f -0.0020493 0.0013483 -1.520 0.1285
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 10.45197 on 5 degrees of freedom
## Residual deviance: 0.60387 on 2 degrees of freedom
## AIC: 38.172
##
## Number of Fisher Scoring iterations: 3
##
##
## > summary(glm(y ~ g + f, family=binomial(link="logit")))
##
## Call:
## glm(formula = y ~ g + f, family = binomial(link = "logit"))
##
## Deviance Residuals:
## 1 2 3 4 5 6
## -0.5507 -0.2781 0.7973 1.1558 -0.3688 -0.6584
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.306643 0.167629 1.829 0.0674 .
## g2 0.405554 0.174560 2.323 0.0202 *
## f -0.000997 0.000665 -1.499 0.1338
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 10.4520 on 5 degrees of freedom
## Residual deviance: 2.9218 on 3 degrees of freedom
## AIC: 38.49
##
## Number of Fisher Scoring iterations: 3
##
##
## > ## The "final model"
## > summary(glm.p84 <- glm(y~g, family=binomial(link="logit")))
##
## Call:
## glm(formula = y ~ g, family = binomial(link = "logit"))
##
## Deviance Residuals:
## 1 2 3 4 5 6
## 0.17150 -0.10947 -0.06177 1.77208 -0.19040 -1.39686
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.1231 0.1140 1.080 0.2801
## g2 0.3985 0.1741 2.289 0.0221 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 10.452 on 5 degrees of freedom
## Residual deviance: 5.173 on 4 degrees of freedom
## AIC: 38.741
##
## Number of Fisher Scoring iterations: 3
##
##
## > op <- par(mfrow = c(2,2), oma = c(0,0,1,0))
##
## > plot(glm.p84) # well ?
##
## > par(op)
##
## > ## Tumour Data (Page 92)
## > counts <- c(22,2,10,16,54,115,19,33,73,11,17,28)
##
## > type <- gl(4,3,12,labels=c("freckle","superficial","nodular","indeterminate"))
##
## > site <- gl(3,1,12,labels=c("head/neck","trunk","extremities"))
##
## > data.frame(counts,type,site)
## counts type site
## 1 22 freckle head/neck
## 2 2 freckle trunk
## 3 10 freckle extremities
## 4 16 superficial head/neck
## 5 54 superficial trunk
## 6 115 superficial extremities
## 7 19 nodular head/neck
## 8 33 nodular trunk
## 9 73 nodular extremities
## 10 11 indeterminate head/neck
## 11 17 indeterminate trunk
## 12 28 indeterminate extremities
##
## > summary(z <- glm(counts ~ type + site, family=poisson()))
##
## Call:
## glm(formula = counts ~ type + site, family = poisson())
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0453 -1.0741 0.1297 0.5857 5.1354
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.7544 0.2040 8.600 < 2e-16 ***
## typesuperficial 1.6940 0.1866 9.079 < 2e-16 ***
## typenodular 1.3020 0.1934 6.731 1.68e-11 ***
## typeindeterminate 0.4990 0.2174 2.295 0.02173 *
## sitetrunk 0.4439 0.1554 2.857 0.00427 **
## siteextremities 1.2010 0.1383 8.683 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 295.203 on 11 degrees of freedom
## Residual deviance: 51.795 on 6 degrees of freedom
## AIC: 122.91
##
## Number of Fisher Scoring iterations: 5
##
##
## > ## Randomized Controlled Trial (Page 93)
## > counts <- c(18,17,15, 20,10,20, 25,13,12)
##
## > outcome <- gl(3, 1, length(counts))
##
## > treatment <- gl(3, 3)
##
## > summary(z <- glm(counts ~ outcome + treatment, family=poisson()))
##
## Call:
## glm(formula = counts ~ outcome + treatment, family = poisson())
##
## Deviance Residuals:
## 1 2 3 4 5 6 7
## -0.67125 0.96272 -0.16965 -0.21999 -0.95552 1.04939 0.84715
## 8 9
## -0.09167 -0.96656
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.045e+00 1.709e-01 17.815 <2e-16 ***
## outcome2 -4.543e-01 2.022e-01 -2.247 0.0246 *
## outcome3 -2.930e-01 1.927e-01 -1.520 0.1285
## treatment2 8.717e-16 2.000e-01 0.000 1.0000
## treatment3 4.557e-16 2.000e-01 0.000 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 10.5814 on 8 degrees of freedom
## Residual deviance: 5.1291 on 4 degrees of freedom
## AIC: 56.761
##
## Number of Fisher Scoring iterations: 4
##
##
## > ## Peptic Ulcers and Blood Groups
## > counts <- c(579, 4219, 911, 4578, 246, 3775, 361, 4532, 291, 5261, 396, 6598)
##
## > group <- gl(2, 1, 12, labels=c("cases","controls"))
##
## > blood <- gl(2, 2, 12, labels=c("A","O"))
##
## > city <- gl(3, 4, 12, labels=c("London","Manchester","Newcastle"))
##
## > cbind(group, blood, city, counts) # gives internal codes for the factors
## group blood city counts
## [1,] 1 1 1 579
## [2,] 2 1 1 4219
## [3,] 1 2 1 911
## [4,] 2 2 1 4578
## [5,] 1 1 2 246
## [6,] 2 1 2 3775
## [7,] 1 2 2 361
## [8,] 2 2 2 4532
## [9,] 1 1 3 291
## [10,] 2 1 3 5261
## [11,] 1 2 3 396
## [12,] 2 2 3 6598
##
## > summary(z1 <- glm(counts ~ group*(city + blood), family=poisson()))
##
## Call:
## glm(formula = counts ~ group * (city + blood), family = poisson())
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## -0.7520 3.0183 0.6099 -2.8137 0.1713 -0.4339 -0.1405 0.3977
## 9 10 11 12
## 0.9318 -2.2691 -0.7742 2.0648
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.39239 0.03476 183.92 < 2e-16 ***
## groupcontrols 1.90813 0.03691 51.69 < 2e-16 ***
## cityManchester -0.89800 0.04815 -18.65 < 2e-16 ***
## cityNewcastle -0.77420 0.04612 -16.79 < 2e-16 ***
## bloodO 0.40187 0.03867 10.39 < 2e-16 ***
## groupcontrols:cityManchester 0.84069 0.05052 16.64 < 2e-16 ***
## groupcontrols:cityNewcastle 1.07287 0.04822 22.25 < 2e-16 ***
## groupcontrols:bloodO -0.23208 0.04043 -5.74 9.46e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 26717.157 on 11 degrees of freedom
## Residual deviance: 29.241 on 4 degrees of freedom
## AIC: 154.32
##
## Number of Fisher Scoring iterations: 3
##
##
## > summary(z2 <- glm(counts ~ group*city + blood, family=poisson()),
## + correlation = TRUE)
##
## Call:
## glm(formula = counts ~ group * city + blood, family = poisson())
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## -3.7688 3.7168 3.2813 -3.4418 -1.7675 0.2387 1.5565 -0.2174
## 9 10 11 12
## -1.1458 -1.4687 1.0218 1.3275
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.51395 0.02663 244.60 <2e-16 ***
## groupcontrols 1.77563 0.02801 63.38 <2e-16 ***
## cityManchester -0.89800 0.04815 -18.65 <2e-16 ***
## cityNewcastle -0.77420 0.04612 -16.79 <2e-16 ***
## bloodO 0.18988 0.01128 16.84 <2e-16 ***
## groupcontrols:cityManchester 0.84069 0.05052 16.64 <2e-16 ***
## groupcontrols:cityNewcastle 1.07287 0.04822 22.25 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 26717.157 on 11 degrees of freedom
## Residual deviance: 62.558 on 5 degrees of freedom
## AIC: 185.63
##
## Number of Fisher Scoring iterations: 4
##
## Correlation of Coefficients:
## (Intercept) groupcontrols cityManchester
## groupcontrols -0.90
## cityManchester -0.52 0.50
## cityNewcastle -0.55 0.52 0.30
## bloodO -0.23 0.00 0.00
## groupcontrols:cityManchester 0.50 -0.55 -0.95
## groupcontrols:cityNewcastle 0.52 -0.58 -0.29
## cityNewcastle bloodO
## groupcontrols
## cityManchester
## cityNewcastle
## bloodO 0.00
## groupcontrols:cityManchester -0.29 0.00
## groupcontrols:cityNewcastle -0.96 0.00
## groupcontrols:cityManchester
## groupcontrols
## cityManchester
## cityNewcastle
## bloodO
## groupcontrols:cityManchester
## groupcontrols:cityNewcastle 0.32
##
##
## > anova(z2, z1, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: counts ~ group * city + blood
## Model 2: counts ~ group * (city + blood)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 5 62.558
## 2 4 29.241 1 33.318 7.827e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
demo(smooth)
##
##
## demo(smooth)
## ---- ~~~~~~
##
## > ### This used to be in example(smooth) before we had package-specific demos
## > # Copyright (C) 1997-2009 The R Core Team
## >
## > require(stats); require(graphics); require(datasets)
##
## > op <- par(mfrow = c(1,1))
##
## > ## The help(smooth) examples:
## > example(smooth, package="stats")
##
## smooth> require(graphics)
##
## smooth> ## see also demo(smooth) !
## smooth>
## smooth> x1 <- c(4, 1, 3, 6, 6, 4, 1, 6, 2, 4, 2) # very artificial
##
## smooth> (x3R <- smooth(x1, "3R")) # 2 iterations of "3"
## 3R Tukey smoother resulting from smooth(x = x1, kind = "3R")
## used 2 iterations
## [1] 3 3 3 6 6 4 4 4 2 2 2
##
## smooth> smooth(x3R, kind = "S")
## S Tukey smoother resulting from smooth(x = x3R, kind = "S")
## changed
## [1] 3 3 3 3 4 4 4 4 2 2 2
##
## smooth> sm.3RS <- function(x, ...)
## smooth+ smooth(smooth(x, "3R", ...), "S", ...)
##
## smooth> y <- c(1, 1, 19:1)
##
## smooth> plot(y, main = "misbehaviour of \"3RSR\"", col.main = 3)
##
## smooth> lines(sm.3RS(y))
##
## smooth> lines(smooth(y))
##
## smooth> lines(smooth(y, "3RSR"), col = 3, lwd = 2) # the horror
##
## smooth> x <- c(8:10, 10, 0, 0, 9, 9)
##
## smooth> plot(x, main = "breakdown of 3R and S and hence 3RSS")
##
## smooth> matlines(cbind(smooth(x, "3R"), smooth(x, "S"), smooth(x, "3RSS"), smooth(x)))
##
## smooth> presidents[is.na(presidents)] <- 0 # silly
##
## smooth> summary(sm3 <- smooth(presidents, "3R"))
## 3R Tukey smoother resulting from
## smooth(x = presidents, kind = "3R") ; n = 120
## used 4 iterations
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 44.0 57.0 54.2 71.0 82.0
##
## smooth> summary(sm2 <- smooth(presidents,"3RSS"))
## 3RSS Tukey smoother resulting from
## smooth(x = presidents, kind = "3RSS") ; n = 120
## used 5 iterations
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 44.00 57.00 55.45 69.00 82.00
##
## smooth> summary(sm <- smooth(presidents))
## 3RS3R Tukey smoother resulting from
## smooth(x = presidents) ; n = 120
## used 7 iterations
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 24.00 44.00 57.00 55.88 69.00 82.00
##
## smooth> all.equal(c(sm2), c(smooth(smooth(sm3, "S"), "S"))) # 3RSS === 3R S S
## [1] TRUE
##
## smooth> all.equal(c(sm), c(smooth(smooth(sm3, "S"), "3R"))) # 3RS3R === 3R S 3R
## [1] TRUE
##
## smooth> plot(presidents, main = "smooth(presidents0, *) : 3R and default 3RS3R")
##
## smooth> lines(sm3, col = 3, lwd = 1.5)
##
## smooth> lines(sm, col = 2, lwd = 1.25)
##
## > ## Didactical investigation:
## >
## > showSmooth <- function(x, leg.x = 1, leg.y = max(x)) {
## + ss <- cbind(x, "3c" = smooth(x, "3", end="copy"),
## + "3" = smooth(x, "3"),
## + "3Rc" = smooth(x, "3R", end="copy"),
## + "3R" = smooth(x, "3R"),
## + sm = smooth(x))
## + k <- ncol(ss) - 1
## + n <- length(x)
## + slwd <- c(1,1,4,1,3,2)
## + slty <- c(0, 2:(k+1))
## + matplot(ss, main = "Tukey Smoothers", ylab = "y ; sm(y)",
## + type= c("p",rep("l",k)), pch= par("pch"), lwd= slwd, lty= slty)
## + legend(leg.x, leg.y,
## + c("Data", "3 (copy)", "3 (Tukey)",
## + "3R (copy)", "3R (Tukey)", "smooth()"),
## + pch= c(par("pch"),rep(-1,k)), col=1:(k+1), lwd= slwd, lty= slty)
## + ss
## + }
##
## > ## 4 simple didactical examples, showing different steps in smooth():
## >
## > for(x in list(c(4, 6, 2, 2, 6, 3, 6, 6, 5, 2),
## + c(3, 2, 1, 4, 5, 1, 3, 2, 4, 5, 2),
## + c(2, 4, 2, 6, 1, 1, 2, 6, 3, 1, 6),
## + x1))
## + print(t(showSmooth(x)))
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## x 4 6 2 2 6 3 6 6 5 2
## 3c 4 4 2 2 3 6 6 6 5 2
## 3 4 4 2 2 3 6 6 6 5 3
## 3Rc 4 4 2 2 3 6 6 6 5 2
## 3R 4 4 2 2 3 6 6 6 5 3
## sm 4 4 4 3 3 6 6 6 5 3
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## x 3 2 1 4 5 1 3 2 4 5 2
## 3c 3 2 2 4 4 3 2 3 4 4 2
## 3 2 2 2 4 4 3 2 3 4 4 4
## 3Rc 3 2 2 4 4 3 3 3 4 4 2
## 3R 2 2 2 4 4 3 3 3 4 4 4
## sm 2 2 2 2 3 3 3 3 4 4 4
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## x 2 4 2 6 1 1 2 6 3 1 6
## 3c 2 2 4 2 1 1 2 3 3 3 6
## 3 2 2 4 2 1 1 2 3 3 3 3
## 3Rc 2 2 2 2 1 1 2 3 3 3 6
## 3R 2 2 2 2 1 1 2 3 3 3 3
## sm 2 2 2 2 2 2 2 3 3 3 3
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## x 4 1 3 6 6 4 1 6 2 4 2
## 3c 4 3 3 6 6 4 4 2 4 2 2
## 3 3 3 3 6 6 4 4 2 4 2 2
## 3Rc 4 3 3 6 6 4 4 4 2 2 2
## 3R 3 3 3 6 6 4 4 4 2 2 2
## sm 3 3 3 3 4 4 4 4 2 2 2
##
## > par(op)