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))
}
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()
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 件のコメント:
コメントを投稿