カテゴリカルデータ解析読書メモ(第2章)

気付いたら買って本棚に飾ってあったので読んでいる。

カテゴリカルデータ解析 (Rで学ぶデータサイエンス 1)

カテゴリカルデータ解析 (Rで学ぶデータサイエンス 1)

準備

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)

f:id:Rion778:20160725003837p:plain
クロス集計をしたい場合は、条件を+で繋ぐ。

> 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) # margin 1:行和を1 2:列和を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), "%"))

f:id:Rion778:20160725004739p:plain
多変数間のクロス集計結果を図示するにはモザイクプロットが便利である。

> Ctable3 <- xtabs(~ Treatment +  Sex + Improved, data = Arthritis)
> mosaicplot(Ctable3, col = TRUE)

f:id:Rion778:20160725005102p:plain
なお、個票形式のデータであれば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)

f:id:Rion778:20160725005654p:plain

表形式→集計形式→個票

既に表形式になってしまっているデータを集計形式に戻すのは簡単で、単に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