Rいろは：グラフ作成ggplot2編

###### ggplot2とは？ #################
# 標準仕様でもそれらしい科学的グラフを作ってくれるRパッケージ
# （plot関数の場合、相当なアレンジをしないと使い物にならない）
# 色分けや複合グラフ、重ね合わせなどに強い、応用の利きやすいメリットがある
# さらにcowplotパッケージを適用して、より学術用に適したグラフ作成を紹介します

# 項や足し算形式でグラフ要素を追加する形態になっている。例えば、散布図、棒グラフのようなグラフ形式や、回帰直線など。
# 項の指定方法などが（plot関数と違って）英単語を元にしているものが多く連想しやすい

# 使い方の早見表も参考になります

install.packages("ggplot2", dep=T) # 未インストールの場合
install.packages("cowplot", dep=T) # 未インストールの場合
install.packages("reshape", dep=T) # 未インストールの場合
require(ggplot2) # 呼び出し

###### まず定番のdata"iris"を用いて、originalのggplot2で描いてみます #########

data(iris) # iris呼び出し
head(iris) # irisの構成をチェック
  #Sepal.Length Sepal.Width Petal.Length Petal.Width Species
#1          5.1         3.5          1.4         0.2  setosa
#2          4.9         3.0          1.4         0.2  setosa
#3          4.7         3.2          1.3         0.2  setosa
#4          4.6         3.1          1.5         0.2  setosa
#5          5.0         3.6          1.4         0.2  setosa
#6          5.4         3.9          1.7         0.4  setosa



##### 散布図 ####################
ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width, fill=Species)) + geom_point(pch=21, size=2)

# ggplot(data=iris, # データ元の指定
#  aes(x=Sepal.Length, y=Sepal.Width, # x, yに用いる変数の指定（aes()内で変数を指定）
#      fill=Species)) # Speciesによって異なる色を使用
# + geom_point( # + geom_xxx: xxxの部分でグラフ種類を指定
#  pch=21,  # pch=21: 塗りと枠に異なる色を適用可能（21~25） # 全種に適用する場合、aes外に記述
#  size=2) # プロットサイズ（少し大きくしました）

# 以下でも同様のグラフになります（fill=Speciesを後半に記述）
ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width)) + geom_point(aes(fill=Species), pch=21, size=2)

# このままでは背景がうるさく、学術用に不向き
require(cowplot) # cowplotパッケージを呼び出し、同じようにグラフを作成します
theme_set(theme_cowplot()) # このコードを通します（すると以降、cowplotスタイルに変更される）

ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width, fill=Species)) + geom_point(pch=21, size=2)
# 背景が白地のスッキリしたものに変わりました。文字のサイズもちょうどいいです。



##### 箱ひげ図 ####################
ggplot(data=iris, aes(x=Species, y=Sepal.Width, fill=Species)) + geom_boxplot()



##### 棒グラフ＋SD ####################
# サンプルデータの準備
MEAN1 <- with(iris, tapply(Sepal.Width, Species, mean)) # 各種の平均（"対象データ", "グループ分け", "関数名"）
SD1 <- with(iris, tapply(Sepal.Width, Species, sd)) # 各グループの標準偏差
iris2 <- data.frame(Species=names(MEAN1), me=MEAN1, sd=SD1)
# グラフの作成
ggplot(data=iris2, aes(x=Species, y=me, fill=Species)) + geom_bar(stat="identity") + geom_linerange(aes(ymax=me + sd, ymin=me - sd))



##### 2要因の棒グラフ + SD ####################
# サンプルデータの準備
iris3 <- rbind( # 1.2倍したサンプルデータを追加
data.frame(Species=names(MEAN1), me=MEAN1, sd=SD1, treat="C"),
data.frame(Species=names(MEAN1), me=1.2*MEAN1, sd=1.2*SD1, treat="T") )
# グラフの作成
ggplot(data=iris3, aes(x=Species, y=me, fill=treat)) + geom_bar(stat="identity", position="dodge") + geom_linerange(aes(ymax=me + sd, ymin=me - sd), position=position_dodge(width=0.9))
# position, dodge: 接したグラフを作成するためのkeyword



##### ヒストグラム ####################
ggplot(iris, aes(x=Sepal.Width, fill=Species)) + geom_histogram()



##### 時系列データのグラフ ####################
data(airquality) # data"airquality"を使用
head(airquality) # データの中身をチェック
  #Ozone Solar.R Wind Temp Month Day
#1    41     190  7.4   67     5   1
#2    36     118  8.0   72     5   2
#3    12     149 12.6   74     5   3
#4    18     313 11.5   62     5   4
#5    NA      NA 14.3   56     5   5
#6    28      NA 14.9   66     5   6

# 時間形式に変換
airquality$Time <- as.POSIXct(strptime(with(airquality, paste(2018, Month, Day, sep="-")), "%Y-%m-%d")) # (2020.03.03少しコードを修正)
ggplot(airquality, aes(x=Time, y=Ozone)) + geom_point(pch=21, cex=2, fill="purple") + geom_line()



##### 時系列データのグラフ（複数系列） ####################
require(reshape) # 要インストール
airquality2 <- melt( # データ構成を変換
airquality[,c("Ozone","Temp","Time")], id.var="Time")
head(airquality2)
        #Time variable value
#1 2018-05-01    Ozone    41
#2 2018-05-02    Ozone    36
#3 2018-05-03    Ozone    12
#4 2018-05-04    Ozone    18
#5 2018-05-05    Ozone    NA
#6 2018-05-06    Ozone    28

# 2パネル構成
base <- ggplot(airquality2, aes(x=Time, y=value)) # ごちゃごちゃしてきたので分割します
poi <- geom_point(pch=21, cex=2, aes(fill=variable)) 
lin <- geom_line(aes(colour=variable)) 
base + poi + lin + facet_wrap(~ variable, scale="free_y", ncol=1) # variableでパネルを分割




###### 散布図＋回帰直線 ##################
summary(iris) # 値の幅をチェック
  #Sepal.Length    Sepal.Width     Petal.Length    Petal.Width          Species
 #Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100   setosa    :50
 #1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300   versicolor:50
 #Median :5.800   Median :3.000   Median :4.350   Median :1.300   virginica :50
 #Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199                
 #3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800                
 #Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500                

md1 <- lm(Sepal.Length ~ Sepal.Width, iris[iris$Species=="setosa",]) # 1種に限定
coef(md1)
#(Intercept) Sepal.Width
  #6.5262226  -0.2233611

# 散布図 + 回帰直線
x.interval <- seq(2, 5, 0.1) # x(Sepal.Width)の最小〜最大を0.1刻みで細かく分割
est.Y <- predict(md1, newdata=data.frame(Sepal.Width=x.interval))
fitted <- data.frame(Sepal.Width=x.interval, Sepal.Length=est.Y)
# fitted <- data.frame(Sepal.Width=x.interval, Sepal.Length=coef(md1)[1] + coef(md1)[2]*x.interval) # 上と同じ
base2 <- ggplot(iris[iris$Species=="setosa",], aes(x=Sepal.Width, y=Sepal.Length)) 
poi2 <- geom_point(fill="turquoise", pch=21, size=2) 
lin2 <- geom_line(data=fitted, colour="turquoise", lwd=1) # ここは"data="を明記する必要あり
base2 + poi2 + lin2



##### 散布図 + 回帰直線 + 信頼区間 ##################
est.Y2 <- predict(md1, newdata=data.frame(Sepal.Width=x.interval), interval="confidence")
head(est.Y2)
       #fit      lwr      upr
#1 4.019981 3.753101 4.286860
#2 4.089030 3.839589 4.338470
#3 4.158079 3.925981 4.390177
#4 4.227128 4.012251 4.442004
#5 4.296177 4.098369 4.493984
#6 4.365226 4.184291 4.546160

fitted2 <- data.frame(Sepal.Width=x.interval, Sepal.Length=est.Y2[,1], upperCI=est.Y2[,3], lowerCI=est.Y2[,2])
band <- geom_ribbon(data=fitted2, aes(ymin=lowerCI, ymax=upperCI), alpha=0.5, fill="turquoise", linetype="blank") # alpha: 半透明(0~1)
base2 + poi2 + band + lin2 # 重ねたい順で



##### 散布図 + 回帰直線（複数グループ） ##################
x.interval3 <- seq(1,5,0.1)
new.dat <- data.frame( # 種数分、繰り返す
Species=rep(unique(iris$Species),each=length(x.interval3)), 
Sepal.Width=rep(x.interval3, length(unique(iris$Species))))
md3 <- lm(Sepal.Length ~ Sepal.Width*Species, iris) # 全種を対象、Speciesによる違い（ANCOVA）
est.Y3 <- predict(md3, newdata=new.dat)
fitted3 <- data.frame(Species=new.dat$Species, Sepal.Width=new.dat$Sepal.Width, Sepal.Length=est.Y3) # Species列を追加
base3 <- ggplot(iris, aes(x=Sepal.Width, y=Sepal.Length)) 
poi3 <- geom_point(aes(fill=Species), pch=21, size=2) 
lin3 <- geom_line(data=fitted3, aes(colour=Species), lwd=1) # ここは"data="を明記する必要あり
base3 + poi3 + lin3




##### 複数タイプのグラフの組み合わせ（cowplot利用） ##################
g1 <- ggplot(data=iris, aes(x=Sepal.Length, y=Sepal.Width)) + geom_point(aes(fill=Species), pch=21, size=2) 
g2 <- ggplot(data=iris, aes(x=Species, y=Petal.Width)) + geom_boxplot(aes(fill=Species)) # 箱ひげ図
+ theme(axis.text.x = element_text(angle=90, vjust=0.5)) # label方向

g12 <- plot_grid(g1, g2, labels="auto", align="h") 
# labels="auto"：小文字でラベリング（大文字は"AUTO"）
# 直接入力の例：labels=c("a","b")
# align="h"：図の端を揃える（"h", "v", "hv"）

#### グラフの調整（例：凡例を統一、サイズ調整）
g1n <- g1 + theme(legend.position="none") # 一方の凡例消去
bd <- panel_border(colour=1, size=1) # 外枠追加
g1n2 <- plot_grid(g1n + bd, g2+ bd, labels="auto", 
   rel_widths=c(1, 1.4), # 横幅を1:1.4に変更
   scale=0.9, # パネル間の間隔を広げる
   vjust=0.2) # ラベル位置（a, b）の調整
 
 
 
#### ファイルの保存 ##############################
save_plot("directry/filename.pdf", g1n2) # cowplotの関数、お任せで適度なサイズになる。指定も可能（2020.03.03: ggsaveからsave_plotへ変更）

#### 水平線・垂線 ##############################
abl <- geom_hline(yintercept=0, linetype="dotted") # 水平 （"dotted": 点線）
abl2 <- geom_vline(xintercept=0, linetype="dotted") # 垂直


#### 軸とラベルの調整 ############################
x.ax <- scale_x_continuous(breaks=c(1:5), limits=c(1, 5), expand=c(0.001, 0.001), name="Sepal width (mm)")
# ex.
base3 + poi3 + lin3 + x.ax


#### 色の指定 ##############################
# グレースケール
set.fillG <- scale_fill_grey(start=0.5, end=1) # グレー 〜 白の場合（塗りつぶし色）
# "fill"の部分を"colour"にすると線の色の指定（"shape": 形状、

set.fillC <- scale_fill_manual(values=c("turquoise", "salmon", "blue")) # 塗りつぶし色
set.colC <- scale_colour_manual(values=c("turquoise", "salmon", "blue")) # 線の色
# ex.
base3 + poi3 + lin3 + set.fillC + set.colC

#### cf. 色コードの選択 ##########################
"見やすい"色の組み合わせ例の参照サイト（ColorBrewer）
fill="#2ca25f" のようにしてカラーコードで色の指定が可能



#### 90° 回転（グラフを横向きにする）#######################
 + coord_flip()


#### 凡例の調整 ##############################
+ theme(legend.position="none") # 凡例消し
+ theme(legend.position="bottom") # ラベルの右揃え（cf. vjust)

# 長すぎるラベルを改行する（ \n が使える）
scale_x_continuous(name = "Length\n (cm)")

# 上付き、ギリシャ文字など
expression(m^2) 
expression(Delta) # Δ



#### cowplotの調整 ##############################
+ panel_border(colour=1, size=1) # defaultの線・文字色がなぜかグレーなので黒に変更
+ theme_cowplot(font_size = 12) # 文字が大きく感じる場合はサイズ指定

# 日本語の利用（ただしアウトライン化されてしまう模様 ）
ggsave(..., family="Japan1GothicBBB",...)