[R program]Rグラフィックの勉強

グラフィックに定評のあるRなので、以下の本の勉強しました。

ggplot2がある今となっては、次の本をやった方が有意義かも。

どちらの本もコードと例は豊富なので勉強にはなります。

上の本で自分的に役に立ったのはこれ。

axTicks(1)      #でx軸の目盛りを取得する

テキストとかを追記したいときに目盛が分かるので便利。

だけどパッケージなんかによっては取得できないことも。

あとは矢印を描いたり、ヒストグラムに密度を付けたりとかいろいろ。

3章までやって力尽きた。。。

##---------１章　作図の紹介
#par(ask=)

plot(pressure)
text(150,600,"Pressure(mm Hg)\nversus\nTemperature(Celsius)")

#---------信頼区間のヒゲ
#------plotCI
#install.packages("gplots")
library("gplots")
example(plotCI)

#------errbar
#install.packages("Hmisc")
library("Hmisc")
example(errbar)

#---------作図出力
postscript(file="myplot.ps")
plot(pressure)
dev.off()

png(file="myplot.png")
plot(pressure)
dev.off()


##---------２章　traditional作図を素朴に使ってみる

example(barplot)

#---------プロットのタイプ
y <- rnorm(20)
plot(y, type="p")
plot(y, type="l")
plot(y, type="b")
plot(y, type="h")

#---------回帰診断
lm.SR <-lm(sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings) 
plot(lm.SR)

#---------クラスタリング
library("cluster")
subset <- sample(1:150, 20)				#ランダムサンプルっぽい
cS <- as.character(Sp <- iris$Species[subset])		#as.characterで文字列のままcSに格納する
cS[Sp == "setosa"] <- "S"
cS[Sp == "versicolor"] <- "V"
cS[Sp == "virginica"] <- "g"
ai <- agnes(iris[subset, 1:4])
plot(ai, labels = cS)

#---------その他
example(matplot)
example(stripchart)
example(curve)
example(stem)

example(labcurve)					#Hmiscパッケージ

#------箱ヒゲ図等
boxplot(decrease ~ treatment, data = OrchardSprays, log = "y", col="light grey")
boxplot(decrease ~ treatment, data = OrchardSprays, log = "y", col="light grey", boxwex=0.5)

barplot(VADeaths[1:2,], angle = c(45, 135), density = 40, col = "grey", names=c("RM", "RF", "UM", "UF"))
barplot(VADeaths[1:2,], angle = c(45, 135), density = 20, col = "grey", names=c("RM", "RF", "UM", "UF"), horiz=TRUE)

y <- rnorm(20)
plot(y, type="l", lwd=3)
plot(y, type="l", col="grey")
plot(y, type="l", lty="dashed")
plot(y, type="l", ylim=c(-4, 4))

#---------多変数のグラフ
example(persp)			#3次元サーフェス
example(contour)		#３番目の変数の値を表す等高線
example(filled.contour)
example(image)			#３番目の変数を表す正方形グリッド
example(symbols)		#３番目の変数をシンボルで表現
example(fourfoldplot)		#2*2*kの分割表
example(pairs)
example(stars)
example(mosaicplot)

#install.packages("scatterplot3d")
#install.packages("rgl")
#install.packages("Rggobi")

library("scatterplot3d")
library("rgl")
library("Rggobi")

example(play3d)			#３次元図を動かす

example("dotchart")
example("coplot")		#条件付きグラフを描く
example("sunflowerplot")	#離散データの重複を花弁で表す

#install.packages("hexbin")
library("hexbin")
example(hexbin)			#散布図を密度で表すから容量の節約になる

example(jitter)			#離散データをランダムにずらす
example(assocplot)		#Cohen-Friendlyの相関グラフ、2*2表各セルのχ２乗値の棒グラフ

example(dendrogram)		#樹形図を描いてその一部を抜き出すこともできる

#install.packages("rpart")
library("rpart")
example(rpart)

#install.packages("maptree")
library("maptree")
example(maptree)		#サンプルプログラムがネット上に。結果は謎

locator(n=2)			#作図デバイス上の座標を取得する
identify(x)			#デバイス上の観測値の識別番号を表示する(Escボタンで抜ける)


#---------traditional作図のカスタマイズ、低水準作図関数

rgb()				#RGB色空間＋アルファチャンネル透過値
hsv()				#HSV色空間
rgb2hsv()
convertColor()			#色空間の変換
hcl()				#CIELUV色空間
palette()			#現在のカラーパレットの確認

#install.packages("RColorBrewer")
library("RColorBrewer")
example(RColorBrewer)		#様々なカラーパレット
brewer.pal(11,"PRGn")		#紫から緑までの11個のパレット

#parでのlend,lljoinオプションで太い線の終わりと継ぎ目を指定
#xaxs="i"：軸のタイプ、bty="o"：グラフの周りを囲む
#par(mfrow=c(3, 2))：グラフ領域を３×２行列に分割

#---------複数のグラフの配置
layout(rbind(c(1,3), c(0,0), c(2,2)), heights=c(2, lcm(0.5), 1), respect=rbind(c(0, 0), c(0, 0), c(1, 0)))
layout.show(3)

#---------plot regionへの追加注記
#------異なる変数のプロット
x <- 1:10
y <- matrix(sort(rnorm(30)), ncol=3)
plot(x, y[, 1], ylim=range(y), ann=F, axes=F, type="l", col="grey")
box(col="grey")

points(x, y[, 1])
lines(x, y[, 2], col="grey")
points(x, y[, 2], pch=2)
lines(x, y[, 3], col="grey")
points(x, y[, 3], pch=3)

example(segments)

#------テキストの追記
x <- c(4, 5, 2, 1)
y <- x
plot(x, y, ann=F, axes=F, col="grey", pch=16)
points(3, 3, col="grey", pch=16)
box(col="grey")

text(x, y, c("bottom", "left", "top", "right"), pos=1:4)		#pos optionで文字が重ならない
text(3, 3, "overlay")



#------データ範囲を囲む破線と凸囲み
x <- rnorm(100)
y <- rnorm(100)
plot(x, y, ann=F, axes=F, col="grey")
box(col="grey")

rect(min(x), min(y), max(x), max(y), lty="dashed")
hull <- chull(x, y)		#凸部分集合?を選ぶ
polygon(x[hull], y[hull])



#------直線、矢印
x <- runif(20, 1, 10)		#１から１０の一様分布
y <- x + rnorm(20)		#xに誤差を付ける
plot(x, y, ann=F, axes=F, col="grey", pch=16)
box(col="grey")

lmfit <- lm(y ~ x)
abline(lmfit)
arrows(5, 8, 7, predict(lmfit, data.frame(x=7)), length=0.1)		#predictで予測モデルの当てはめを行う
text(5, 8, "Line of best fit", pos=2)


#------ヒストグラムに密度を追加（絨毯プロット)
y <- rnorm(10)
hist(y, main="", xlab="", ylab="", axes=F, border="grey", col="light grey")
box(col="grey")
rug(y, ticksize=0.02)


#------ マージンにテキストを描く
y1 <- rnorm(100)
y2 <- rnorm(100)
par(mfrow=c(2, 1), xpd=NA)	#クリッピング領域がデバイス全体になる

plot(y1, type="l", axes=F, xlab="", ylab="", main="")
box(col="grey")
mtext("Left end of margin", adj=0, side=3)
lines(x=c(20, 20, 40, 40), y=c(-7, max(y1), max(y1), -7), lwd=3, col="grey")	#xpd=NAにしているから下の四角と繋がる

plot(y2, type="l", axes=F, xlab="", ylab="", main="")
box(col="grey")
mtext("Right end of margin", adj=1, side=3)
mtext("Label below x=30", at=30, side=1)
lines(x=c(20, 20, 40, 40), y=c(7, min(y2), min(y2), 7), lwd=3, col="grey")


#------凡例
with(iris, plot(Sepal.Length, Sepal.Width, pch=as.numeric(Species), cex=1.2))
legend(6.1, 4.4, c("sentosa", "versicolor", "virginica"), cex=1.5, pch=1:3)

barplot(VADeaths[1:2,], angle=c(45,135), density=20, col="grey", names=c("RM", "RF", "UM", "UF"))
legend(0.4, 38, c("55-59", "50-54"), cex=1.5, angle=c(135, 45), density=20, fill="grey")


#------軸
x <- 1:2
y <- runif(2, 0, 100)
par(mar=c(4, 4, 2, 4))
plot(x, y, type="n", xlim=c(0.5, 2.5), ylim=c(-10, 110), axes=F, ann=F)

axis(2, at=seq(0, 100, 20))
mtext("Temperature (Centigrade)", side=2, line=3)

axis(1, at=1:2, labels=c("Treatment 1", "Treatment 2"))
axis(4, at=seq(0, 100, 20), labels=seq(0, 100, 20)*9/5 + 32)

mtext("Temperature (Fahrenheit)", side=4, line=3)
box()

segments(x, 0, x, 100, lwd=20, col="dark grey")
segments(x, 0, x, 100, lwd=16, col="white")
segments(x, 0, x, y, lwd=16, col="light grey")

#axTicks(1)でx軸の目盛りを取得する


#------数式
#expression(bar(x) == sum(frac(x[i],n), i==1, n))
#expression(hat(beta) == (X^t * X)^{-1} * X^t *y)
#expression(z[i] == sqrt(x[i]^2 + y[i]^2))


#------座標系（奇麗だけど使い方が分からないな）
plot(0:1, 0:1, type="n", axes=F, ann=F)
usr <- par("usr")
pin <- par("pin")
xcm <- diff(usr[1:2])/(pin[1]*2.54)
ycm <- diff(usr[3:4])/(pin[2]*2.54)

par(xpd=NA)
rect(0 + 0.2*xcm, 0 - 0.2*ycm, 1 + 0.2*xcm, 1- 0.2*ycm, col="grey", border=NA)

rect(0, 0, 1, 1, col="white")
segments(seq(1, 8, 0.1)*xcm, 0, seq(1, 8, 0.1)*xcm, c(rep(c(0.5, rep(0.25, 4), 0.35, rep(0.25, 4)), 7), 0.5)*ycm)
text(1:8*xcm, 0.6*ycm, 0:7, adj=c(0.5, 0))
text(8.2*xcm, 0.6*ycm, "cm", adj=c(0, 0))


#------重ね合わせ
drunkenness <- ts(c(3875, 4846, 5128, 5773, 7327, 6688, 5582, 3473, 3186, rep(NA, 51)), start=1912, end=1971)

#---方法１
par(mar=c(5, 6, 2, 4))
plot(drunkenness, lwd=3, col="grey", ann=F, las=2)
mtext("Drunkenness\nRelated Arrests", side=2, line=3.5)
par(new=T)
plot(nhtemp, ann=F, axes=F)
mtext("Temperature (F)", side=4, line=3)
title("Using par(new=T)")
axis(4)

#---方法２
par(mar=c(5, 6, 2, 4))
plot(drunkenness, lwd=3, col="grey", ann=F, las=2)
mtext("Drunkenness\nRelated Arrests", side=2, line=3.5)
usr <- par("usr")
par(usr=c(usr[1:2], 47.6, 54.9))
lines(nhtemp)
mtext("Temperature (F)", side=4, line=3)
title("Using par(usr=...)")
axis(4)

#---高水準関数の追記方法
with(trees, 
	{
	plot(Height, Volume, pch=3, xlab="Height (ft)", ylab=expression(paste("Volume ", (ft^3))))
	symbols(Height, Volume, circles=Girth/12, fg="grey", inches=F, add=T)
	})
		
		
#---複雑なグラグ
xx <- c(1:50)
yy <- rnorm(50)
n <- 50
hline <- 0

plot(yy ~ xx, type="n", axes=F, ann=F)
polygon(c(xx[1], xx, xx[n]), c(min(yy), yy, min(yy)), col="grey", border=NA)

usr <- par("usr")
rect(usr[1], usr[3], usr[2], hline, col="white", border=NA)

lines(xx, yy)

abline(h=hline, col="grey")
box()
axis(1)
axis(2)


#---ビットマップ
#install.packages("pixmap")
library("pixmap")
example(addlogo)

moonPhase <- function(x, y, phase, size=.07) {
  # phase 1: first quarter
  #       2: full
  #       3: last quarter
  #       4: new
  # size is in inches
  n <- 17
  angle <- seq(0, 2*pi, length=n)
  xx <- x + cos(angle)*xinch(size)
  yy <- y + sin(angle)*yinch(size)
  if (phase == 4)
    fill <- "black"
  else
    fill <- "white"
  polygon(xx, yy, col=fill)
  if (phase == 1)
    polygon(xx[(n/4):(n*3/4) + 1],
            yy[(n/4):(n*3/4) + 1],
            col="black")
  if (phase == 3)
    polygon(xx[c(1:(n/4 + 1), (n*3/4 + 1):n)],
            yy[c(1:(n/4 + 1), (n*3/4 + 1):n)],
            col="black")
}

# Data from Land Information New Zealand
# http://hydro.linz.govt.nz
# + 1 for daylight saving
hours <- c(18, 18, 19, 20, 21, 22, 23,
           0, 1, 2, 3, 4, 5, 5,
           6, 7, 8, 9, 10, 11, 12,
           13, 13, 14, 15, 15, 16) + 1
hours[7] <- 0 # 23 + 1 = 0 on a 24-hour clock
mins <- c(9, 57, 49, 46, 48, 54, 59,
          59, 52, 41, 28, 14, 1, 52,
          36, 43, 41, 39, 38, 35, 26,
          10, 49, 26, 2, 39, 16)
lowTideDate <- ISOdatetime(2005, 2, c(1:6,8:28),
                           hours, mins, 0)
lowTideHour <- ISOdatetime(2005, 2, 1,
                           hours, mins, 0)
phases <- ISOdatetime(2005, 2, c(2, 9, 16, 24),
                      c(19, 10, 12, 16) + 1,
                      c(28, 30, 16, 55), 0)
mainHours <- ISOdatetime(2005, 2, 1,
                         c(0, 4, 8, 12, 16, 20, 23), 
                         c(rep(0, 6), 59), 
                         c(rep(0, 6), 59))

library(pixmap)
# Original image from NASA
# http://grin.hq.nasa.gov/ABSTRACTS/GPN-2000-000473.html
moon <- read.pnm(system.file(file.path("Images", 
                                       "GPN-2000-000473halfsize.pnm"), 
                             package="RGraphics"))
par(pty="s", xaxs="i", yaxs="i", cex=.7)
plot.new()
addlogo(moon, 0:1, 0:1, asp=1)
par(new=TRUE, xaxs="r", yaxs="r", las=1)
plot(lowTideDate, lowTideHour, type="n",
     ylim=range(mainHours), axes=FALSE, ann=FALSE)
# dashed reference lines
midday <- ISOdatetime(2005, 2, 1, 12, 0, 0)
abline(h=midday, v=phases,
       col="white", lty="dashed")
# grey "repeat" of tide info to show gradient
lines(lowTideDate[6:7], 
      c(ISOdatetime(2005, 1, 31,
                    hours[6], mins[6], 0),
        lowTideHour[7]),
      lwd=2, col="grey50")
points(lowTideDate[6:7], 
       c(ISOdatetime(2005, 1, 31,
                     hours[6], mins[6], 0),
         lowTideHour[7]),
       pch=16, cex=.7, col="grey50")
for (subset in list(1:6, 7:27)) {
  lines(lowTideDate[subset], lowTideHour[subset],
        lwd=2, col="white")
  points(lowTideDate[subset], lowTideHour[subset],
         pch=16, cex=.7, col="white")
}
box()
axis.POSIXct(1, lowTideDate)
axis.POSIXct(2, at=mainHours, format="%H:%M")
mtext("Time of Low Tide (NZDT)", side=2, line=4, las=0)
mtext("Auckland, New Zealand 2005", side=1, line=3)
axis(3, at=phases, labels=FALSE)
par(xpd=NA)
ymax <- par("usr")[4]
for (i in 1:4)
  moonPhase(phases[i], ymax + yinch(.2), c(3, 4, 1, 2)[i])
mtext("Phases of the Moon", side=3, line=3)


#---軸の取得
#棒グラフ
y <- sample(1:10)		#ランダム並べ替え
midpts <- barplot(y, col=" light grey")
width <- diff(midpts[1:2])/4
left <- rep(midpts, y - 1) - width
right <- rep(midpts, y - 1) + width
heights <- unlist(apply(matrix(y, ncol=10), 
                        2, seq))[-cumsum(y)]
segments(left, heights, right, heights,
         col="white")

#箱ヒゲ図
with(ToothGrowth, 
     {
       boxplot(len ~ supp, border="grey", 
               col="light grey", boxwex=0.5)
       points(jitter(rep(1:2, each=30), 0.5), 
              unlist(split(len, supp)),
              cex=0.5, pch=16)
     })


#---地図
#install.packages("maps")
library(maps)
coplot(lat ~ long | depth, data = quakes, number=4, 
       panel=function(x, y, ...) {
         usr <- par("usr")
         rect(usr[1], usr[3], usr[2], usr[4], col="white")
         map("world2", regions=c("New Zealand", "Fiji"),
             add=TRUE, lwd=0.1, fill=TRUE, col="grey")
         text(180, -13, "Fiji", adj=1, cex=0.7)
         text(170, -35, "NZ", cex=0.7)
         points(x, y, pch=".")
       })


#---３次元
#install.packages("RGraphics")
library(RGraphics)
par(mar=rep(0, 4), lwd=0.1)
z <- 2 * volcano        
x <- 10 * (1:nrow(z))   
y <- 10 * (1:ncol(z))   
trans <- persp(x, y, z, theta = 135, phi = 30, 
               scale = FALSE, ltheta = -120, 
               # shade=0.5, border=NA, 
               box = FALSE)
box(col="grey", lwd=1)

trans3d <- function(x,y,z,pmat) {
  tmat <- cbind(x,y,z,1)%*% pmat
  tmat[,1:2] / tmat[,4]
}		#persp関数が作った２次元変換行列を使って変換、４次元行列なのはなぜ？

summit <- trans3d(x[20], y[31], max(z), trans)
points(summit[1], summit[2], pch=16)
summitlabel <- trans3d(x[20], y[31], max(z) + 50, trans)
text(summitlabel[1], summitlabel[2], "Summit")

drawRoad <- function(x, y, z, trans) {
  road <- trans3d(x, y, z, trans)
  lines(road[,1], road[,2], lwd=5)
  lines(road[,1], road[,2], lwd=3, col="grey")
}
with(volcano.summitRoad,
     drawRoad(srx, sry, srz, trans))
with(volcano.upDownRoad,
     {
       clipudx <- udx
       clipudx[udx < 230 & udy < 300 | 
               udx < 150 & udy > 300] <- NA
       drawRoad(clipudx, udy, udz, trans)
     })
with(volcano.accessRoad,
     drawRoad(arx, ary, arz, trans))


#---新規のグラフを作る
plot.new()
plot.window(range(pressure$temperature), range(pressure$pressure))
plot.xy(pressure, type="p")
box()
axis(1)
axis(2)

plot(pressure)

#複雑な棒グラフ
groups <- c("cows", "sheep", "horses", 
            "elephants", "giraffes")
males <- sample(1:10, 5)
females <- sample(1:10, 5)

par(mar=c(0.5, 5, 0.5, 1))

plot.new()
plot.window(xlim=c(-10, 10), ylim=c(-1.5, 5.5))

ticks <- seq(-10, 10, 5)
y <- 1:5
h <- 0.2

lines(rep(0, 2), c(-1.5, 5.5), col="grey")
segments(-10, y, 10, y, lty="dotted")
rect(-males, y-h, 0, y+h, col="dark grey")
rect(0, y-h, females, y+h, col="light grey")
mtext(groups, at=y, adj=1, side=2, las=2)
par(cex.axis=0.5, mex=0.5)
axis(1, at=ticks, labels=abs(ticks), pos=0)		#pos=0でplot region内にx軸を描く

tw <- 1.5*strwidth("females")				#凡例の長さを指定
rect(-tw, -1-h, 0, -1+h, col="dark grey")
rect(0, -1-h, tw, -1+h, col="light grey")
text(0, -1, "males", pos=2)
text(0, -1, "females", pos=4)				#xもyも同じ点から左右に描いている

box("inner", col="grey")