1
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

ggplot2でfor構文を用いる(種ごとの作図を行う場合など)

Last updated at Posted at 2019-06-08

#ggplot2とfor構文を組み合わせて、一気に複数枚の作図をする

階層群集モデルの結果をCSVファイルに出力後、種ごとの結果を作図する際に使用したコードです。
前半は気合が入ってますが、後半に行くにつれてとてもコードが汚くなっています...

最終的には、こんな画像が種ごとに描けます。
INtSpecies0.png


library("ggplot2")                    ##ggplot2パッケージの読み込み
windowsFonts("YuGo" = windowsFont("游ゴシック"))    ##フォントの読み込み
library(gridExtra)
library(onemap)
library(png)
library(grid)

setwd("c:/Users/...") ##ディレクトリの指定
comm3 <- read.csv("結果ファイル.csv", header = F)                 ##MCMCの結果が入ったcsvファイル
comm4 <- read.csv("タイトルファイル.csv", header = F, stringsAsFactors=F)        ##一枚一枚の作図のタイトルを入力したcsvファイル

#新たな297行3列のマトリックスを作成
#MCMCの結果が3種のハイパーパラメータと33種のパラメーター9個で、297行あるため
#MCMCの出力結果から描写に使うのが、中央値・95%信用区間の上端と下端なので3列分
x <- matrix(NA,297,3)

x[,1] <- comm3[,4] # comm3の結果のうちの4列目のみ抜き出してxに付加
x[,2] <- comm3[,5] # comm3の結果のうちの5列目のみ抜き出してxに付加
x[,3] <- comm3[,6] # comm3の結果のうちの6列目のみ抜き出してxに付加

#新たな21行1列のマトリックスを作成
p <- matrix(NA,297,1)

p[,1] <- comm4[,1] 
# 別途作成していたタイトルファイルについても行列を作成、ここでは種名(F1. Eurasian Skylark)に該当する文が33種ぶん含まれている


####################################################################
##################plotting each species#############################
####################################################################
value <-  0:32
for (z in 0:32){

res <- paste("plot",z, sep="")
    eval(call("<-", as.name(z),value))
sum <- data.frame(sx = c(5,4,3,2,1),
                  sy = c(exp(x[9*z+1,1]),exp(x[9*z+2,1]),
                       exp(x[9*z+3,1]),exp(x[9*z+4,1]),exp(x[9*z+5,1])))
####データを指定する数列を作成、ここではパラメータW, A, P, C, Rの部分を指定


sum2 <- data.frame(sx = c(1,2,3,4),
                   sy = c(x[9*z+9,1],x[9*z+8,1],x[9*z+7,1],x[9*z+6,1]))                
####データを指定する数列を作成、ここではパラメータTemp, Prec, Temp:A, Prec:Aの部分を指定

####以下、コード汚い
c1 <-exp(x[9*z+1,2])                                             ####Lower 95% CI of W
c2 <-exp(x[9*z+2,2])                                             ####Lower 95% CI of A
c3 <-exp(x[9*z+3,2])                                             ####Lower 95% CI of P
c4 <-exp(x[9*z+4,2])                                             ####Lower 95% CI of C
c5 <-exp(x[9*z+5,2])                                             ####Lower 95% CI of R 
c6 <-exp(x[9*z+1,3])                                             ####Upper 95% CI of W
c7 <-exp(x[9*z+2,3])                                             ####Upper 95% CI of A
c8 <-exp(x[9*z+3,3])                                             ####Upper 95% CI of P
c9 <-exp(x[9*z+4,3])                                             ####Upper 95% CI of C
c10<-exp(x[9*z+5,3])                                             ####Upper 95% CI of R

tem1 <-x[9*z+6,2]                                                ####Lower 95% CI of Temp
tem2 <-x[9*z+6,3]                                                ####Upper 95% CI of Temp
pre1 <-x[9*z+7,2]                                                ####Lower 95% CI of Prec
pre2 <-x[9*z+7,3]                                                ####Upper 95% CI of Prec
tema1 <-x[9*z+8,2]                                               ####Lower 95% CI of Temp
tema2 <-x[9*z+8,3]                                               ####Upper 95% CI of Temp
prea1 <-x[9*z+9,2]                                               ####Lower 95% CI of Prec
prea2 <-x[9*z+9,3]                                               ####Upper 95% CI of Prec
name  <-p[9*z+1,1]                                               ####name


####ここから作図、まず上半分
lnduse <- ggplot(comm3) +
       theme_bw(base_size=11, base_family='') +
       xlim(0,5) +  
       scale_x_reverse() +
       scale_x_continuous(breaks=c(1,2,3,4,5),
                          labels=c("R","C","P","A","W"))+ 
       xlab("") + ylab("Abundance (/3 ha)")+ 
       geom_errorbar(mapping = aes(x = 1, ymin=c5, ymax =c10),
                        size = 1, width = 0.3, color = "black") + 
       geom_errorbar(mapping = aes(x = 2, ymin=c4 , ymax =c9),
                        size = 1, width = 0.3, color = "black") +
       geom_errorbar(mapping = aes(x = 3, ymin=c3 , ymax =c8),
                        size = 1, width = 0.3, color = "black") +
       geom_errorbar(mapping = aes(x = 4, ymin=c2 , ymax =c7),
                        size = 1, width = 0.3, color = "black") +
       geom_errorbar(mapping = aes(x = 5, ymin=c1 , ymax =c6),
                        size = 1, width = 0.3, color = "black") +
       geom_point(data = sum, aes(x = sx, y = sy), 
                  size = 3, color = "black") +
       labs(title = name) +
       theme(plot.title = element_text(hjust = 0.5),
             text = element_text(size = 12, family = "YuGo"))+
       theme(axis.title.x = element_text(size=12, family = "YuGo"), 
             axis.title.y = element_text(size=8, family = "YuGo")) +
       theme(axis.text.x = element_text(size=12, family = "YuGo"),
             axis.text.y = element_text(size=14, family = "YuGo")) +
       guides(fill = FALSE) +
       coord_flip()
       

####ここ以下も作図、下半分に相当
coeff  <- ggplot(comm3) +
       theme_bw(base_size=11, base_family='') +
       xlim(0,4.2) +  
       scale_x_reverse() +
       scale_x_continuous(breaks=c(1,2,3,4),
                          labels=c("Prec:A", "Temp:A","Prec ","Temp "))+ 
       xlab("") + ylab("Coefficients")+ 
       geom_errorbar(mapping = aes(x = 1, ymin=prea1, ymax =prea2),
                        size = 1, width = 0.3, color = "black") + 
       geom_errorbar(mapping = aes(x = 2, ymin=tema1, ymax =tema2),
                        size = 1, width = 0.3, color = "black") +
       geom_errorbar(mapping = aes(x = 3, ymin=pre1, ymax =pre2),
                        size = 1, width = 0.3, color = "black") + 
       geom_errorbar(mapping = aes(x = 4, ymin=tem1, ymax =tem2),
                        size = 1, width = 0.3, color = "black") + 
       geom_segment(aes(x=0.6, y=0, xend=4.2, yend=0), linetype = 2,size=0.7, color = "black") +
       geom_point(data = sum2, aes(x = sx, y = sy), 
                  size = 3, color = "black") +
       theme(plot.title = element_text(hjust = 0.5),
             text = element_text(size = 14, family = "YuGo"))+
       theme(axis.title.x = element_text(size=12, family = "YuGo"), 
             axis.title.y = element_text(size=8, family = "YuGo")) +
       theme(axis.text.x = element_text(size=12, family = "YuGo"),
             axis.text.y = element_text(size=14, family = "YuGo")) +
       guides(fill = FALSE) +
       coord_flip()

res <- ggdraw() +
          draw_plot(lnduse, x = 0.12, y = 0.43, width = 0.78, height = 0.52) +
          draw_plot(coeff,  x = 0, y = 0, width = 0.9, height = 0.40)
res
ggsave(plot = res, file = paste("RES190526/INtSpecies", z, ".png", sep=""),
       width = 4.3, height = 4.8, dpi = 400)}
1
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
1
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?