Rで散布図とヒストグラム - mecobalamin’s diary

irisのデータから散布図とヒストグラムを作る
前回も使ったけどirisのデータセットって？

irisの正体を調べた方がいた
irisの正体 (R Advent Calendar 2012 6日目) - どんな鳥も
"アヤメ属の3種 (Iris setosa，Iris versicolor，Iris virginica) の"、
"sepal (がく片) とpetal (花弁) の長さ(length)・幅(width)"の測定結果
ということだそう。

前回は種ごとに長さと幅の平均値を計算してグラフにしたが
今回は散布図とヒストグラムを書いて分布を調べる

結果はこんな感じ
f:id:mecobalamin:20180912113934p:plain

setosa(緑)は他の2種にくらべてがく片の幅が大きいが
それ以外は小さく分布も離れている
versicolor(赤)とvirginica(青)はがく片のサイズは似ているが
花弁のサイズで区別がつきそう
virginicaは外れ値なのかもう1グループありそうにも見える

コードはこちら
MakingSDGraph()からGraphSD()とGraphHist()を呼び出している
irisのデータセットの読み込み方は前回と同じ

Data.Iris <- iris
# setosa, versicolor, virginica
Species.List <- as.character(unique(Data.Iris$Species))
# Sepal.Length, Sepal.Width, Petal.Length, Petal.Width, Species
Col.Data <- colnames(Data.Iris)
# グラフの作成
MakingSDGraph(Data.Iris, Species.List, Col.Data[5], Col.Data)

MakingSDGraph()はpngファイルの書き込む用意をして
GraphSD()とGraphHsit()に渡すデータセットを整形する

MakingSDGraph <- function(data.all, sp.list, col.name, col.list){

  # データごとの色の定義
  data.color <- data.all[col.name]

  # グラフを書き込むファイルの用意, 5x5の領域を確保する
  png.name = sprintf("SD-Hist.png")
  png(png.name, width = 1200, height = 1200)
  par(mfcol = c(4, 4), lwd = 2, mar = c(5, 5, 1, 1), ps = 12)

  n <- c(); m <- c()
  for (i in 1:4){
    for (j in 1:4){
      par1 = as.character(col.list[i]); par2 = as.character(col.list[j])

      # 各列を変数に代入
      xT <- data.all[[par1]]; yT <- data.all[[par2]]
      xyT <- as.data.frame(cbind(xT, yT, data.color))

      # 最小値、最大値を求める
      minimum <- c(min(xT), min(yT)); maximum <- c(max(xT), max(yT))
      minimum.sign <- sign(minimum); maximum.sign <- sign(maximum)
      minimum <- floor(minimum * minimum.sign) * minimum.sign - 1
      maximum <- ceiling(maximum * maximum.sign) * maximum.sign + 1

      if(par1 != par2){
        GraphSD(par1, par2, minimum, maximum, xyT, sp.list)
      }
      else{
        GraphHist(par1, minimum, maximum, xyT, sp.list)
      }
    }
  }
  dev.off()
}

GraphSD()で散布図を作成する

GraphSD <- function(p1, p2, min, max, d, s){
    plot(0, 0, type = "n",
    xlim = c(min[1], max[1]),
    ylim = c(min[2], max[2]),
    xlab = p1,
    ylab = p2,
    cex.lab = 2,
    cex.axis = 2
  )
  # 点を描く
  d[,3] <- gsub(s[1], "#9bbb5980", d[,3])
  d[,3] <- gsub(s[2], "#c0504d80", d[,3])
  d[,3] <- gsub(s[3], "#4f81bd80", d[,3])
  points(d[,1], d[,2], col = d[,3], pch = 16, cex = 2.5)
}

GraphHist()でヒストグラムを作成する
ビン幅はkで定義しているがちょうどよいところは何度か試す
スタージェスの公式で求めた値そのままよりも3分の1にしたほうが良さげ

GraphHist <- function(p1, min, max, x, s){
  Lower <- min[1] - 1 ; Higher <- max[1] + 1
  k <- seq(Lower, Higher, (Higher - Lower) / (1 + log2(nrow(x))) / 3)

  # リガンドごとにデータを分ける
  sub_xr <- subset(x, (x[,3] == s[2]))
  sub_xb <- subset(x, (x[,3] == s[3]))
  sub_xg <- subset(x, (x[,3] == s[1]))

  # 分布の最大値を求める
  hr <- hist(as.numeric(sub_xr[,1]), breaks = k, plot = FALSE)
  hb <- hist(as.numeric(sub_xb[,1]), breaks = k, plot = FALSE)
  hg <- hist(as.numeric(sub_xg[,1]), breaks = k, plot = FALSE)
  max.hist <- max(hr$count, hb$count, hg$count)

  # ヒストグラム赤
  hist(as.numeric(sub_xr[,1]), col = "#c0504d60",
    breaks = k, axes = F,
    xlim = c(Lower, Higher),
    ylim = c(0, max.hist),
    main = "", xlab = p1,
    cex.lab = 2,
  )

  # ヒストグラム青
  hist(as.numeric(sub_xb[,1]), col = "#4f81bd60",
    breaks = k, axes = F,
    xlim = c(Lower, Higher),
    main = "", xlab = p1,
    add = TRUE
  )

  #  ヒストグラム緑
  hist(as.numeric(sub_xg[,1]), col = "#9bbb5960",
    breaks = k, axes = F,
    xlim = c(Lower, Higher),
    main = "", xlab = p1,
    add = TRUE
  )

  Axis(side = 2, cex.axis = 2)
  Axis(side = 1, at = (Lower + 1):(Higher - 1),
    cex.axis = 2,
    lty = NULL)
}