2014年5月11日日曜日

ggplotでモザイクプロット


Rにはmosaicplot()があって、何の問題もなくモザイクプロットが描けるのに、ggplotでモザイクプロットを描けないかという人もいるよ。なんでggplotが良いのかさっぱりわからないという話もあるのは承知の上で、私もまたggplotを使ってモザイクプロットを描いてみました。軟弱なので、標準化残差の計算などは当然chisq.test()に頼ることになるのです。標準化残差の区切り方も手抜き。
でも、質的変数ならば、タイルとタイルの間は隙間を開けておきたいなあなどと、こだわったりなんかもするのです。

ggmosaic <- function(data , x, y, weight = 1, stdres = FALSE, 
    gap = c(2,2), correct = TRUE){
  require(ggplot2)
  gap <- rep(gap,2)[1:2]/100
  data <- cbind(data, 1)
  tbl <- xtabs( data[,deparse(substitute(weight))] ~ data[,deparse(substitute(x))] +  
    data[,deparse(substitute(y))])

  colnames(tbl)[is.na(colnames(tbl))] <- "<NA>"
  rownames(tbl)[is.na(rownames(tbl))] <- "<NA>"
  ptbl <- prop.table(tbl,1)                               #行の和が1になるような比率に変換
  chi2test <- chisq.test(tbl, correct = correct)
  dat <- data.frame(
    xname = factor(rep(rownames(tbl), ncol(tbl)), levels = rownames(tbl)),        #x軸のラベル
    xmin = rep(c(0, cumsum(rowSums(tbl)))[1:nrow(tbl)], ncol(tbl) ),              #xの最小値
    xmax = rep(cumsum(rowSums(tbl)), ncol(tbl) ),                                 #xの最大値
    yname = factor(rep(colnames(tbl), each = nrow(tbl)), levels = colnames(tbl)), #yのラベル
    ymin = as.data.frame.table(t(apply(cbind(0, ptbl)[,1:ncol(tbl)], 1, cumsum)))[,3],           #yの最小値
    ymax = as.data.frame.table(t(apply(ptbl, 1, cumsum)))[,3],                                   #yの最大値
    stdres = cut(as.data.frame.table(chi2test$stdres)[,3], breaks = c(-Inf,seq(-4,4,2),Inf)))    #残差の段階
  dat$xmin <- dat$xmin + c(0, cumsum(rep(sum(tbl) * gap[1], (nrow(tbl) - 1))))    #xのギャップ挿入
  dat$xmax <- dat$xmax + c(0, cumsum(rep(sum(tbl) * gap[1], (nrow(tbl) - 1))))    #xのギャップを挿入
  dat$ymin <- dat$ymin + rep((0:(ncol(tbl)-1))/ (ncol(tbl)-1), each = nrow(tbl)) * gap[2]  #yのギャップ挿入
  dat$ymax <- dat$ymax + rep((0:(ncol(tbl)-1))/ (ncol(tbl)-1), each = nrow(tbl)) * gap[2]  #yのギャップ挿入
 
  fig <- ggplot(dat) +
    geom_rect(aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax, fill = yname),
              color = 1)  +
    scale_fill_manual(name="", values = terrain.colors(ncol(tbl)),
                      guide = guide_legend(reverse = TRUE)) +
    scale_x_continuous(breaks = rowSums(dat[1:nrow(tbl),2:3])/2, labels= rownames(tbl)) +
    scale_y_continuous(breaks = tapply((dat$ymin + dat$ymax)/2, INDEX = dat$yname, FUN = mean),
    labels = colnames(tbl)) +
    labs(x="",y="",z="")
  if(stdres){   
    fig <- fig +
      geom_point(aes(x= (xmin + xmax) / 2, y = (ymin + ymax) / 2, color = stdres), shape = 15, size = 4) +
      scale_color_brewer(name="標準化残差", palette = "RdBu", drop = FALSE,
      guide = guide_legend(reverse = TRUE)) +
      geom_point(aes(x= (xmin + xmax) / 2, y = (ymin + ymax) / 2), color = "black", shape = 7, size = 4)
  }
  print(fig)
  print(chi2test)
  invisible(list(figure = fig, test = chi2test))
}

といったものを書いてみました。

##例えば
DF <- as.data.frame(UCBAdmissions)
ggmosaic(data = DF, x = Dept, y = Admit, weight = Freq, stdres = TRUE)
##横にしてみる
DF2 <- cbind(DF, GD = paste(DF$Gender,DF$Dept))
ggfig <- ggmosaic(data = DF2, x = GD, y = Admit, weight = Freq, stdres = FALSE)
ggfig$figure + coord_flip()

次のようなモザイクプロットが描けるのではないかと思います。で、なんでggplotなんだってか?

 

0 件のコメント:

コメントを投稿