#======================================================
#	2 Ｒにおけるグラフィックスの初歩
#======================================================

#======================================================
#	2.1 伝統的なグラフィックス
#======================================================

#------------<基本的なグラフのタイプ>------------
x <- 1:10
y <- c(1, 3, 6, 4, 2, 5, 8, 6, 7, 10)

plot(x, y)   　          #ベーシックな散布図

#図2.1
op <- par(mfrow=c(2, 3))              #グラフを２行３列に配置
plot(x, y, type="l", main='"type=l"') #lines（線グラフ）
plot(x, y, type="h", main='"type=h"') #high density（高密度グラフ）
plot(x, y, type="b", main='"type=b"') #both（線と記号の双方によるグラフ）
plot(x, y, type="o", main='"type=o"') #overplotted（上記の変種）
plot(x, y, type="s", main='"type=s"') #steps（階段グラフ）
plot(x, y, type="S", main='"type=S"') #other steps（上記の変種）
par(op)                               #グラフ配列を元に戻す。

#------------<１次元データ>------------

#------------<ヒストグラム>------------
市区町村$高齢者人口比率 <- 市区町村$人口65歳以上/(市区町村$人口15歳未満 + 
	市区町村$人口15〜64歳 + 市区町村$人口65歳以上) * 100

#図2.2
hist(市区町村$高齢者人口比率, main="市区町村別高齢者人口比率（2010年）",
	xlab="高齢者人口比率", ylab="市区町村の数", col="steelblue")

colors() #指定可能な色の名前を確認

#色の指定を数字で行い、度数ではなく構成比で表示する。
#図2.3
hist(市区町村$高齢者人口比率, main="市区町村別高齢者人口比率（2010年）",
	xlab="高齢者人口比率", ylab="構成比", col=1:11, freq=F)

#------------<密度プロット>------------

#図2.4
plot(density(市区町村$高齢者人口比率), main="市区町村別高齢者人口比率（2010年）",
	xlab="高齢者人口比率", ylab="確率密度")
#ラグプロットを追加
rug(市区町村$高齢者人口比率)

#------------<ＱＱプロット>------------

#正規ＱＱプロット
#図2.5
qqnorm(市区町村$高齢者人口比率, main="市区町村別高齢者人口比率（2010年）")
#正規分布に従うとした場合の直線を描く
qqline(市区町村$高齢者人口比率)

#------------<円グラフ>-----------
#連続データを集計
bks <- c(5, seq(from=20, to=40, by=5), 60) #比率の区分
tab <- table(cut(市区町村$高齢者人口比率, breaks=bks))
tab

#名前を付ける。
names(tab) <- c("20％以下","20％超〜\n25％以下","25％超〜\n30％以下",
	"30％超〜\n35％以下","35％超〜\n40％以下","40％超") 

#反時計回り（極座標）
#図2.6
pie(tab)

#時計回り
#図2.7
pie(tab, clockwise=T)

#------------<点グラフ>------------

#図2.8
dotchart(unlist(tab), xlab="市区町村の数", main="高齢者人口比率")

#x軸の範囲を指定
#図2.9
dotchart(tab, xlab="市区町村の数", main="高齢者人口比率", xlim=c(0, 600))

#------------<棒グラフ>------------

#図2.10
op <- par(mfrow=c(3, 1)) #グラフィックス・パラメータの設定
#y軸のラベルを横書きにし（las=1）、テキスト（市区町村の数）を表示する。
x <- barplot(tab, las=1, main="(a)")
text(x, tab+70, tab)
#グラフ領域外にも、文字列が表示されるようにする（xpd=T）。
x <- barplot(tab, las=1, main="(b)")
text(x, tab+70, tab, xpd=T)
#あるいは、y軸の範囲を指定する。
x <- barplot(tab, las=1, ylim=c(0, 700), main="(c)")
text(x, tab+70, tab)
box() #グラフに四角形の枠を描く
par(op)　#グラフィックス・パラメータを元に戻す。

#------------<２次元データ>------------
#------------<散布図：連続変数×連続変数>------------

#図2.11
#x, yの順に変数を指定
plot(市区町村$人口密度, 市区町村$高齢者人口比率, log="xy",
	xlab="人口密度", ylab="高齢者人口比率")

#線形回帰などのモデル式と同じ記法
plot(高齢者人口比率 ~ 人口密度, data=市区町村, log="xy")

#------------<箱ひげ図：カテゴリー×連続変数>------------

#図2.12
boxplot(高齢者人口比率 ~ 市区町村別, 
	ylab="高齢者人口比率", data=市区町村)

#------------<モザイクプロット：カテゴリー×カテゴリー>------------
#カテゴリカル変数を新しく作る。
市区町村$第１次産業比率 <- 市区町村$第１次産業/(市区町村$第１次産業 + 
        市区町村$第２次産業 + 市区町村$第３次産業) * 100
(quant <- quantile(市区町村$第１次産業比率)) #四分位点


市区町村$第１次産業クラス <- ordered(0, 
        levels=c("低い", "やや低い", "やや高い", "高い"))
市区町村$第１次産業クラス[市区町村$第１次産業比率 < quant[2]] <- "低い"
市区町村$第１次産業クラス[市区町村$第１次産業比率 >= quant[2] &
        市区町村$第１次産業比率 < quant[3]] <- "やや低い"
市区町村$第１次産業クラス[市区町村$第１次産業比率 >= quant[3] &
        市区町村$第１次産業比率 < quant[4]] <- "やや高い"
市区町村$第１次産業クラス[市区町村$第１次産業比率 >= quant[4]] <- "高い"

#クロス集計
(tab2 <- table(市区町村$市区町村別, 市区町村$第１次産業クラス))

#図2.13
mosaicplot(tab2, main="市区町村別×第１次産業クラス", las=1, 
        col="springgreen")

#------------<３次元データ>------------
#------------<等高線グラフ>------------

#図2.14
rx <- range(x <- 10*1:nrow(volcano)) #xの値の範囲
ry <- range(y <- 10*1:ncol(volcano)) #yの実際の値の範囲
ry <- ry + c(-1, 1) * (diff(rx) - diff(ry))/2 #yの範囲に上下の余白部分を加える。
tcol <- terrain.colors(12) #色の設定
#空のプロットでプロット領域を確保
plot(x = 0, y = 0, type = "n", xlim = rx, ylim = ry, xlab = "", ylab = "")
u <- par("usr") #ユーザー領域 c(x1, x2, y1, y2) を得る。
#ユーザー領域全体に色を塗り、領域の境界の色を赤とする。
rect(u[1], u[3], u[2], u[4], col = tcol[8], border = "red")
#等高線を描く。
contour(x, y, volcano, col = tcol[2], lty = "solid", add = TRUE,
     vfont = c("sans serif", "plain"))
title("A Topographic Map of Maunga Whau", font = 4) #メインタイトル

#------------<イメージ>------------

#図2.15
x <- 10*(1:nrow(volcano))　#xの格子点
y <- 10*(1:ncol(volcano))　#yの格子点
#イメージの描画
image(x, y, volcano, col = terrain.colors(100), axes = FALSE)
#等高線の追加
contour(x, y, volcano, levels = seq(90, 200, by = 5),
             add = TRUE, col = "peru")
axis(1, at = seq(100, 800, by = 100))　#x軸
axis(2, at = seq(100, 600, by = 100))　#y軸
box()
title(main = "Maunga Whau Volcano", font.main = 4) #メインタイトル

#------------<鳥瞰図>------------

demo(persp)

#------------<多次元データ>------------
#------------<散布図行列>------------

#図2.16
pairs(~ log(人口密度) + log(高齢者人口比率) + log(第１次産業比率), 
        data=市区町村, panel=panel.smooth, main="市区町村(2010年)",
        lwd=3)

#======================================================
#	2.2 トレリス・グラフィックス
#======================================================

library(lattice) #パッケージのロード

#------------<１次元データ>------------
#------------<ヒストグラム>------------

#図2.17
histogram(~ 高齢者人口比率, data=市区町村, 
        main="高齢者人口比率（2010年）",
        xlab="高齢者人口比率", ylab="市区町村の数", col="steelblue")

#------------<密度プロット>------------

#図2.18
densityplot(~ 高齢者人口比率, data=市区町村, 
        main="高齢者人口比率（2010年）",
        xlab="高齢者人口比率", ylab="確率密度")

#------------<ＱＱプロット>------------

#正規ＱＱプロット
#図2.19
qqmath(~ 高齢者人口比率, data=市区町村, 
	main="高齢者人口比率（2010年）",
       panel = function(x, ...) {
       	panel.qqmathline(x, ...)
       	panel.qqmath(x, ...)
       })

#指数分布に従うとした場合のＱＱプロット
#図2.20
qqmath(~ 高齢者人口比率, data=市区町村, 
        main="高齢者人口比率（2010年）\n指数分布との比較",
        distribution = qexp,	　#指数分布
        panel = function(x, ...) {
               panel.qqmathline(x, ...)
               panel.qqmath(x, ...)
       })

#------------<点グラフ>------------

#連続データを集計
bks <- c(5, seq(from=20, to=40, by=5), 60) #比率の区分
tab <- table(cut(市区町村$高齢者人口比率, breaks=bks))
tab

#名前を付ける。
names(tab) <- c("20％以下","20％超〜\n25％以下","25％超〜\n30％以下",
	"30％超〜\n35％以下","35％超〜\n40％以下","40％超")

#図2.21
dotplot(tab, xlab="市区町村の数", main="高齢者人口比率", 
        xlim=c(0, 600))

#------------<棒グラフ>------------

#図2.22
barchart(tab, xlab="市区町村の数", main="高齢者人口比率")

#関数revを使用して順序を反転
barchart(rev(tab), xlab="市区町村の数", main="高齢者人口比率")

#------------<２次元データ>------------
#------------<散布図：連続変数×連続変数>------------

#図2.23
xyplot(高齢者人口比率 ~ 人口密度, data=市区町村, 
        xlab="人口密度", ylab="高齢者人口比率", 
        scales = list(x = list(log = 10, equispaced.log=F), 
        y = list(log = 10, equispaced.log=F)))

#同じ
xyplot(高齢者人口比率 ~ 人口密度, data=市区町村, 
        xlab="人口密度", ylab="高齢者人口比率", 
        scales = list(log=10, equispaced.log=F))
#------------<箱ひげ図：カテゴリー×連続変数>------------

#図2.24
bwplot(高齢者人口比率 ~ 市区町村別, data=市区町村)

#------------<３次元データ>------------

datax <- rep(seq(from=-1.5, to=1.5, length=50), times=50)
datay <- rep(seq(from=-1.5, to=1.5, length=50), times = rep(50, times=50))
dataz <- exp(-(datax^2 + datay^2 + datax*datay))
gauss <- data.frame(datax, datay, dataz)

#------------<等高線グラフ>------------

#図2.25
contourplot(dataz ~ datax*datay, data=gauss, aspect=1, 
	at = seq(from=0.1, to=0.9, by=0.2))

#------------<イメージ>------------

#図2.26
levelplot(dataz ~ datax*datay, data=gauss, aspect=1, cuts=6)

#------------<鳥瞰図>------------

#図2.27
wireframe(dataz ~ datax*datay, data=gauss, screen = list(z=45, x=-60, y=0))

#------------<３次元散布図>------------

#図2.28
cloud(log(第１次産業比率) ~ log(人口密度) + log(高齢者人口比率), data=市区町村,
	screen = list(z=-30, x=-60, y=0), 
		main="市区町村(2010年)")

#パッケージscatterplot3dによるプロット
library(scatterplot3d)    #パッケージscatterplot3dのロード
scatterplot3d(log(市区町村$第１次産業比率) ~ log(市区町村$人口密度) + log(市区町村$高齢者人口比率), data=市区町村,
	screen = list(z=-30, x=-60, y=0), 
		main="市区町村(2010年)")
scatterplot3d(datax, datay, dataz)

#------------<動的グラフ>------------

library(rgl)　#パッケージrglのロード
plot3d(datax, datay, dataz, col=rainbow(1000))

#------------<多次元データ>------------
#------------<散布図行列>------------

#図2.29
splom(~ data.frame(log(人口密度), log(高齢者人口比率), log(第１次産業比率)), 
       data=市区町村,
	main="市区町村(2010年)",
	panel=function(...) {
		panel.splom(...)
		panel.loess(col="red", lwd=3, ...)
	})

#------------<因子による条件付け>------------

#図2.30
xyplot(log(高齢者人口比率)~log(人口密度)|市区町村別,
 	data=市区町村, layout=c(1, 3),
	as.table=T, main="市区町村(2010年)",
	panel = function(...) {
		panel.xyplot(...)  
		panel.loess(..., col="red", lwd=3)
 	})

#======================================================
#	2.3 ggplot2
#======================================================
library(ggplot2) #パッケージのロード 
#------------<１次元データ>------------
#------------<ヒストグラム>------------

#図2.31
qplot(高齢者人口比率, data=市区町村, main="高齢者人口比率（2010年）",
      xlab="高齢者人口比率", ylab="市区町村の数", xlim=c(5, 60), 
      binwidth=5, fill=I("steelblue"))

#------------<密度プロット>------------

#図2.32
p <- qplot(高齢者人口比率, data=市区町村, 
        main="高齢者人口比率（2010年）",
        xlab="高齢者人口比率", ylab="市区町村の数",geom="density", 
        fill=I("steelblue"))
p + geom_rug() #ラグプロットの追加

#------------<ＱＱプロット>------------

#図2.33
p <- qplot(sample=高齢者人口比率, data=市区町村, 
	main="高齢者人口比率（2010年）", stat="qq")
#参照線の作成
y <- quantile(市区町村$高齢者人口比率, c(0.25, 0.75))
x <- qnorm(c(0.25, 0.75))
slope <- diff(y)/diff(x)
int <- y[1L] - slope * x[1L]
p + geom_abline(slope=slope, intercept=int)

#------------<円グラフ>------------
#図2.34
bks <- c(5, seq(from=20, to=40, by=5), 60) #比率の区分
tab <- table(cut(市区町村$高齢者人口比率, breaks=bks))
names(tab) <- c("20％以下","20％超〜\n25％以下","25％超〜\n30％以下",
	"30％超〜\n35％以下","35％超〜\n40％以下","40％超") #名前を付ける。
df <- data.frame(x = names(tab), y = as.numeric(tab)) #データフレームに変換
y.breaks <- cumsum(df$y) - df$y/2 #ラベルの位置

p <- ggplot(df, aes(x=1, y=y, fill=x))
p <- p + geom_bar(stat="identity", color="black") #棒グラフ
p <- p + coord_polar(theta="y")                   #極座標
p <- p + theme(axis.ticks=element_blank(),        #不要な要素を削除
        axis.title=element_blank(),
        axis.text.y=element_blank(),
        legend.position="none")
p + scale_y_continuous(breaks=y.breaks, labels=df$x) #ラベルを付ける。

#------------<点グラフ>------------

#図2.35
qplot(y, names(tab), data=df, main="高齢者人口比率（2010年）",
      xlab="市区町村の数", ylab="高齢者人口比率", xlim=c(0, 600), geom="point")

#------------<棒グラフ>------------

#図2.36
p <- qplot(names(tab), y, data=df, main="高齢者人口比率（2010年）",
      xlab="高齢者人口比率", ylab="市区町村の数", ylim=c(0, 600), 
      fill=I("steelblue"), geom="bar", stat="identity")
p + coord_flip() #x軸とy軸を入れ替える。

#------------<２次元データ>------------
#------------<散布図>------------

#図2.37
p <- qplot(人口密度, 高齢者人口比率, data = 市区町村, log="xy",
	geom = "point", colour=I("steelblue"),
	main="人口密度と高齢者人口比率（対数尺度）")
p + stat_smooth(method="loess", size=2, colour="tomato")

#------------<箱ひげ図>------------

#図2.38
qplot(市区町村別, 高齢者人口比率, data=市区町村, geom="boxplot")

#------------<３次元データ>------------
#------------<等高線グラフ>------------

#図2.39
library(reshape2) #関数meltを使用するため
volcano3d <- melt(volcano)
names(volcano3d) <- c("x", "y", "z")

p <- ggplot(volcano3d, aes(x, y, z=z))
p + stat_contour(binwidth=2)

#------------<イメージ>------------

#図2.40
p <- qplot(x, y, data=volcano3d, fill=z, geom="tile")
p + stat_contour(binwidth=2, aes(x, y, z=z))

#------------<多次元データ>------------
#------------<散布図行列>------------

#図2.41
library(GGally) #関数ggpairsを使用するため

data <- data.frame(人口密度=log(市区町村$人口密度), 
         高齢者人口比率=log(市区町村$高齢者人口比率),
         第１次産業比率=log(市区町村$第１次産業比率))
ggpairs(data)



