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なんだってか？

共分散構造分析[R編]―構造方程式モデリング ついにでた。

ついに、待望の「共分散構造分析[R編]―構造方程式モデリング」が出版されました。AmosやOpenMxといった商用ソフトが容易に使える環境にない私には良い報せです。幾分天邪鬼な私は、`楽天ブックス`に注文していましたが、この金曜日に配送されてきました。手に取ると[Amos編]より幾分重い感じがしました。ページ数は1割増で、内容もぎっしりです。やはり'lavaan'を使用した内容になっていました。当分の間、本物の「共分散構造分析[R編]」を読みつつ勉強したいと思います。

[Amos編]の再現も、7章まではやってみようと思っていましたが、ここでこのプロジェクトも一区切りという感じです。答えがわかっていてやっているのだから気楽なものですが、それでも収束させることができなかったモデルも残ってしまいました。残念。気が向いたらまた試してみることにしよう。

生兵法は怪我のもと。ついにブログ閉鎖の危機?

先日、見慣れないサイトからトラフィックがあったので、リンクをたどってみると中澤先生のサイトでした。「Rによる統計解析の基礎」持ってます。「なんでちゃんとした先生のサイトにリンクされているんだ、変だな」と思って見てみると、関西大学の水本先生のツィートにつながっていて、そこに豊田先生が共分散構造分析[R編]を出版される！という情報がありました。ついにですか。

私がおもしろ半分で始めたブログもついに終りの日を迎えることになるのでしょうか。このようなパチモンの存在は、豊田先生も許されるはずもありません。いやむしろ、眼中にないか・・・。

水本先生は「これと比べてみるのが楽しみ」などと書かれておられましたが、「これ」がもし私のブログなら、そもそも比べる対象が間違っていますよ(μガンダムとザクＩを比べるのか！！「スペック違い過ぎ」とせがれが言ってました)。 私は、統計学なんぞきちんと勉強したこともない地方の一介のサラリーマンに過ぎません。ここは恥ずかしいのでブログを閉鎖するのが良いのか、パチモンとしても水本先生のお楽しみのために残しておくのが良いのか。どちらもちょっと悲しい。

今はっきりしていることは、・・・・・・・・・・・・・・・共分散構造分析[R編]が出版されたら必ず買います。
Amos編の例では収束すらさせられないモデルがあるのですよ、私なんか。どうしたら良いのかわかるのかな。

lavaanを使われているかもですね。それともスクラッチからプログラムしちゃうのでしょうか？

共分散構造分析[R編]につまづき、CATDAP1に逃避

SEMパッケージは収束条件が厳しい感じで、完全情報最尤推定法でつまずいたのをきっかけに、このブログもまたまた随分放置していました()。

　全く話が変わるのですが、随分以前、統計数理研究所のセミナーで坂本先生のcatdapの話を聞きました。先生は、データが集まったら、とりあえずcatdapで変数選択してみるみたいな紹介をされたように記憶しています。便利かもと思って、作ったのが下の関数。実は、香川大学経済学部教授でいらした故堀啓三先生がSPSSマクロで作成されていたcatdapを下敷きに、ほぼ盗作？
　いまや本家本元の統計数理研究所でcatdap1とcatdap2のスクリプトが公開されている(ここ)ので、下にあるような怪しげなスクリプトを使うメリットはあるはずないよね、でへっ。いや、絶対に使ってはいけないぞ。
　実は、人に見せる集計表を選択するとき、よく、このcatdapを使っているのだけれど、自分はAICを説明しきれないのでcatdapは日陰者あつかいです。変数選択は頭の中でひねり出したことになってます。

下記スクリプトはcatdap01というからには、質的変数にしか使えないのです。

catdap01(y.data, x.data, zero.cell = exp(-1), nvar = 2, printout = 10)

• x.data : 従属変数ベクトル
• y.data : 独立変数候補のデータフレーム
• nvar : 考慮する独立変数候補の数(4層まで)
• zero.cell : 0になったセルに入れる値
• printout : コンソールに出力するデータの数

最近のPCはパワフルなので、変数候補が10くらいあっても3層ぐらいなら、茶を飲んで待てるくらい(2014年現在)で終了すると思います。

catdap01(Titanic[,4], Titanic[,-4], nvar = 1)だと、統計数理研究所のcatdapのヘルプにあるcatdap1と同じになり、
catdap01( (Titanic[,4], Titanic[,-4], nvar = 3)だと、統計数理研究所のcatdap2を使って、
catdap2(Titanic[,4], 2, Titanic[,-4])と同じ値が出力されるはず。

nvar=1　のとき、表示に問題があるので修正(2014/03/20)

catdap01 <- function (y.data, x.data, zero.cell = exp(-1), nvar = 2, printout = 10)
{
    if (!is.factor(y.data)) {
        y.data <- as.factor(y.data)
    }
    printout <- trunc(printout)
    if (printout < 1) {
        printout <- 1
    }
    nvar <- min(nvar, 4)
    nexp <- ncol(x.data)
    aicbox <- c(0)
    varbox <- list("independent")
    ct <- list()
    varlist <- NULL
    cat("\n\ncombining variables...\n")
    varlist <- combn(1:nexp, 1, simplify = FALSE)
    if (nvar > 1) {
        varlist <- c(varlist, combn(1:nexp, 2, simplify = FALSE))
    }
    if (nvar > 2) {
        varlist <- c(varlist, combn(1:nexp, 3, simplify = FALSE))
    }
    if (nvar > 3) {
        varlist <- c(varlist, combn(1:nexp, 4, simplify = FALSE))
    }
    cat("tabulating ", length(varlist), "combinations....\n")
    fun1 <- function(x, zero.cell) {
        x[x == 0] <- zero.cell
        return(x)
    }
    calcaic <- function(x) {
        if (is.vector(x))
            x <- as.matrix(x)
        rn <- ncol(x)
        cn <- nrow(x)
        csum <- colSums(x)
        rsum <- rowSums(x)
        total <- sum(x)
        aic.temp <- sum(x * log((total * x)/(rsum %*% t(csum))))
        aic.temp <- (-2) * aic.temp + 2 * (rn - 1) * (cn - 1)
        return(aic.temp)
    }
    rvarlist <- lapply(varlist, function(x) rev(x) + 1)
    rvarlist <- lapply(rvarlist, function(x) c(1, x))
    all.var <- data.frame(y.data, x.data)
    temp.table <- NULL
    temp.table <- lapply(rvarlist, function(x) data.frame(xtabs(~.,
        data = all.var[, x])))
    temp.table <- lapply(temp.table, function(x) matrix(x$Freq,
        length(x$Freq)/length(levels(y.data)), length(levels(y.data)),
        byrow = TRUE))
    temp.table <- lapply(temp.table, function(x) x[rowSums(x) >  0, ])
    temp.table <- lapply(temp.table, function(x) fun1(x, zero.cell))
    cat("calculating AICs...\n\n")
    aicbox <- sapply(temp.table, function(x) calcaic(x))
    aicorder <- order(aicbox)
    aicbox <- aicbox[aicorder]
    varlist <- varlist[aicorder]
    vars <- sapply(varlist, function(x) c(colnames(x.data)[rev(x)],
        "", "", "")[1:nvar])
    if(is.vector(vars){ vars <- as.matrix(vars) }else{ vars <- t(vars) }
    output.mat <- data.frame(round(aicbox, 3), vars)
    rownames(output.mat) <- NULL
    colnames(output.mat) <- c("AIC", paste("explanatory variables",
        1:nvar))
    cat("\n<<< catdap result>>>\n\n")
    print(output.mat[1:min(printout, 100), ])
    result <- list(output = output.mat, aics = aicbox, explvar = as.list(varlist))
    return(invisible(result))
}