気付いたら買って本棚に飾ってあったので読んでいる。
準備
install.packages("vcd")
library(vcd)
data("Arthritis")
Arthritisは腱鞘炎に関する臨床試験のデータで、処置と経過、性別と年齢のデータが個票として入っている。
個票形式のデータの集計とプロット
xtabs()は表形式の集計を行う。集計の仕方を式の形で与えることができる。
> xtabs(~ Improved, data = Arthritis)
Improved
None Some Marked
42 14 28
このままbarplot()に与えれば棒グラフが作成できる。
IMP2 <- IMP[order(IMP, decreasing = TRUE)]
barplot(IMP2)
クロス集計をしたい場合は、条件を+で繋ぐ。
> xtabs(~ Treatment + Improved, data = Arthritis)
Improved
Treatment None Some Marked
Placebo 29 7 7
Treated 13 7 21
3変数以上は配列になる。
> xtabs(~ Treatment + Improved + Sex, data = Arthritis)
, , Sex = Female
Improved
Treatment None Some Marked
Placebo 19 7 6
Treated 6 5 16
, , Sex = Male
Improved
Treatment None Some Marked
Placebo 10 0 1
Treated 7 2 5
I()を用いると、特定の変数の組み合わせをまとめて行や列に指定できる。表形式に抑えたい場合に便利だ。
> xtabs(~ I(Sex:Treatment) + Improved, data = Arthritis)
Improved
I(Sex:Treatment) None Some Marked
Female:Placebo 19 7 6
Female:Treated 6 5 16
Male:Placebo 10 0 1
Male:Treated 7 2 5
集計結果を比率で表示したい場合はprop.table()を用いる。margin=引数によって行和または列和を1にするように制御できる。
> prop.table(Ctable, margin = 1)
Improved
Treatment None Some Marked
Placebo 0.6744186 0.1627907 0.1627907
Treated 0.3170732 0.1707317 0.5121951
prop.table()で比率に変換した表をbarplot()に渡すことで、帯グラフが作成できる。
> Ctable.p <- prop.table(Ctable, margin = 1)
> barplot(t(Ctable.p), horiz = T, axes = F,
+ xlab = "Improved(None Some Marked)")
> axis(1, seq(0, 1, by = 0.2), paste(seq(0, 100, by = 20), "%"))
多変数間のクロス集計結果を図示するにはモザイクプロットが便利である。
> Ctable3 <- xtabs(~ Treatment + Sex + Improved, data = Arthritis)
> mosaicplot(Ctable3, col = TRUE)
なお、個票形式のデータであればxtabsにより集計する必要はなく、直接mosaicplot()に渡せばよい。
mosaicplot(~ Treatment + Sex + Improved, data = Arthritis,
col = TRUE)
集計データの扱い
集計データとは、変数の組み合わせとその度数をまとめた形式のデータである。
> data("DanishWelfare")
> head(DanishWelfare)
Freq Alcohol Income Status Urban
1 1 <1 0-50 Widow Copenhagen
2 4 <1 0-50 Widow SubCopenhagen
3 1 <1 0-50 Widow LargeCity
4 8 <1 0-50 Widow City
5 6 <1 0-50 Widow Country
6 14 <1 0-50 Married Copenhagen
集計データをxtabs()により集計する際には、式の左辺に度数を指定する必要がある。
> xtabs(Freq ~ Alcohol, data = DanishWelfare)
Alcohol
<1 1-2 >2
2339 2083 722
モザイクプロットを作成したい場合、xtabs()による集計を省略することはできないので注意する。
Ctable4 <- xtabs(Freq ~ Status + Alcohol, data = DanishWelfare)
mosaicplot(Ctable4, col = TRUE)
表形式→集計形式→個票
既に表形式になってしまっているデータを集計形式に戻すのは簡単で、単にdata.frame()に渡せばよい。
> dtable <- matrix(c(2, 18, 3, 17, 9, 11, 15, 5, 18, 2), nr = 2)
> dimnames(dtable) <- list(effect = list("あり", "なし"),
+ dose = list("1", "2", "3", "4", "5"))
> class(dtable) <- "table"
> dose_effect <- data.frame(dtable)
> head(dose_effect)
effect dose Freq
1 あり 1 2
2 なし 1 18
3 あり 2 3
4 なし 2 17
5 あり 3 9
6 なし 3 11
集計形式から個票形式に戻すにはapply()などを駆使する。
> DW <- data.frame(lapply(DanishWelfare,
+ function(i) rep(i, DanishWelfare[,"Freq"])))
> DW2 <- DW[-1]
> head(DW2)
Alcohol Income Status Urban
1 <1 0-50 Widow Copenhagen
2 <1 0-50 Widow SubCopenhagen
3 <1 0-50 Widow SubCopenhagen
4 <1 0-50 Widow SubCopenhagen
5 <1 0-50 Widow SubCopenhagen
6 <1 0-50 Widow LargeCity
カテゴリーの編集
カテゴリー変数を編集したい場合は、直接書き換えてしまうことができる。
> levels(DanishWelfare$Income)
[1] "0-50" "50-100" "100-150" ">150"
> levels(DanishWelfare$Income) <- c("0-100", "0-100",
+ "100<", "100<")
> xtabs(Freq ~ Income, data = DanishWelfare)
Income
0-100 100<
2042 3102
量的変数を大きさで区切ってカテゴリー変数にしたい場合は、cut()が便利だ。
> Arthritis$CatAge <- cut(Arthritis$Age,
+ breaks = 2:8 * 10,
+ labels = 2:7 * 10,
+ right = TRUE)
> head(Arthritis)
ID Treatment Sex Age Improved CatAge
1 57 Treated Male 27 Some 20
2 46 Treated Male 29 None 20
3 77 Treated Male 30 None 20
4 17 Treated Male 32 Marked 30
5 36 Treated Male 46 Marked 40
6 23 Treated Male 58 Marked 50
