R では質的変数、すなわち順序尺度、名義尺度を表現する方法として「因子 (factor)」が用意されています。モデル構築にあたっては、順序尺度は「順序付き因子 (ordered factor)」、名義尺度は「順序無し因子 (unordered factor)」として正しくデータを表現しておく必要がありますので、よく理解しておきましょう。
因子とは何か？
因子 (factor) は、質的変数を R で表現する方法である。質的変数には、順序尺度と名義尺度があるが、R ではそれぞれに対して順序付き因子と順序無し因子が対応する。因子ではない数値ベクトルは、量的変数 (間隔尺度や比例尺度) とみなされる。R では順序無し因子、順序付き因子、因子ではない数値ベクトルは、モデル構築において異なる取り扱い (各尺度として適切な取り扱い) がされるようになっている。
因子の使用例
オーストラリアの会計士30人の本拠地所在州と収入のデータを分析するとする。state は州名略称の文字列ベクトル、incomes は収入の数値ベクトルとする。
> state <- c("tas", "sa", "qld", "nsw", "nsw", "nt", "wa", "wa", "qld", "vic", "nsw", "vic", "qld", "qld", "sa",
+ "tas", "sa", "nt", "wa", "vic", "qld", "nsw", "nsw", "wa", "sa", "act", "nsw", "vic", "vic", "act")
> incomes <- c(60, 49, 40, 61, 64, 60, 59, 54, 62, 69, 70, 42, 56, 61, 61,
+ 61, 58, 51, 48, 65, 49, 49, 41, 48, 52, 46, 59, 46, 58, 43)
> state
[1] "tas" "sa" "qld" "nsw" "nsw" "nt" "wa" "wa" "qld" "vic" "nsw" "vic" "qld" "qld" "sa" "tas" "sa" "nt" "wa" "vic"
[21] "qld" "nsw" "nsw" "wa" "sa" "act" "nsw" "vic" "vic" "act"
> incomes
[1] 60 49 40 61 64 60 59 54 62 69 70 42 56 61 61 61 58 51 48 65 49 49 41 48 52 46 59 46 58 43
ちなみに略称と州名の関係は以下の通り。州と準州合わせて8つある。
| act | the Australian Capital Territory |
| nsw | New South Wales |
| nt | the Northern Territory (準州) |
| qld | Queensland |
| sa | South Australia |
| tas | Tasmania |
| vic | Victoria |
| wa | Western Australia |
因子の作成
州名は名義尺度なので、state から順序無し因子 statef を作る。
> statef <- factor(state) > statef [1] tas sa qld nsw nsw nt wa wa qld vic nsw vic qld qld sa tas sa nt wa vic qld nsw nsw wa sa act nsw vic vic act Levels: act nsw nt qld sa tas vic wa
statef を表示してみると、要素が "" で括られずに表示された。これは文字列ベクトルではなくなったためである。Levels 欄には「水準」が表示されている。因子の水準とは、その因子の取り得る値だ。水準はアルファベット順にソートされているように見える。
因子の水準の取得
> levels(statef) [1] "act" "nsw" "nt" "qld" "sa" "tas" "vic" "wa" > nlevels(statef) [1] 8
levels() で水準の一覧が、nlevels() で水準の数が取れる。
因子の性質
因子の性質を把握するため、いろいろ表示してみる。
> summary(statef) act nsw nt qld sa tas vic wa 2 6 2 5 4 2 5 4
summary() を使うと、各水準の度数が表示された。
> mode(statef) [1] "numeric"
モードは numeric になっており、実体は数値ベクトルのようだ。
> attributes(statef) $levels [1] "act" "nsw" "nt" "qld" "sa" "tas" "vic" "wa" $class [1] "factor"
属性には levels に水準を持っており、class は factor だ。
> unclass(statef) [1] 6 5 4 2 2 3 8 8 4 7 2 7 4 4 5 6 5 3 8 7 4 2 2 8 5 1 2 7 7 1 attr(,"levels") [1] "act" "nsw" "nt" "qld" "sa" "tas" "vic" "wa"
unclass() してみると、実体は水準の添字を要素として持つ数値ベクトルになっているようだ。
つまり、因子は factor クラスのオブジェクトで、水準に辞書順に数値を振り、要素を水準の添字で置き換えた数値ベクトルということのようだ。
因子を使ったデータ分析
因子を活用して、収入の分析を行ってみよう。まずは復習がてら、全標本の標本平均を出してみる。
> mean(incomes) [1] 54.73333
続いて、州毎の標本平均を出してみる。tapply() という関数を用いる。
> incmeans <- tapply(incomes, statef, mean)
> incmeans
act nsw nt qld sa tas vic wa
44.50000 57.33333 55.50000 53.60000 55.00000 60.50000 56.00000 52.25000
tapply(X, INDEX, FUN) は、ベクトル X を因子 INDEX を用いて水準別に分解し、それぞれに関数 FUN を適用して、その結果を水準名で名前付けした配列で返すもの。incomes と statef は同じ要素数なので、incomes に因子 statef を適用できる。適用することで、incomes のどの要素が statef のどの水準に対応するかが決定される。それを用いて水準別に incomes の要素を分解し、水準ごとに1つの数値ベクトルを作る (8州あるので、計8個)。その8州ごとに分解されたベクトルをそれぞれ mean() に渡すことで、各州の標本平均が算出される。結果は配列に格納される。
州ごとの標本平均に対する標準誤差を計算するとする。標準誤差を計算する組込み関数がないので、独自に定義する (stderr() は標準エラー出力オブジェクトだ)。
> stderror <- function(x) sqrt(var(x)/length(x))
> incstderror <- tapply(incomes, statef, stderror)
> incstderror
act nsw nt qld sa tas vic wa
1.500000 4.310195 4.500000 4.106093 2.738613 0.500000 5.244044 2.657536
対象がベクトルのときは tapply() を用いるが、配列・行列の場合は apply() を用いる。また、tapply() の第2引数は factor でなくとも as.factor() で強制変換することで因子になりうるものであれば、指定できる。今回の例では state はその条件を満たしているので、tapply(incomes, state, mean) とすることも可能だった。
母集団の平均の95%信頼区間を求めてみる。母集団の平均が $\mu$、標本平均が $\mu_0$、標準偏差が $s$、標本数が $n$ のとき、$t = (\mu_0-\mu)/(\frac{s}{\sqrt{n}})$ が自由度 $n-1$ の $t$-分布になることを利用すると、$\mu$ の信頼区間は $\mu_0 \pm t s \sqrt{n}$ として求まる。
まずは全標本に対して。
95% にあたる t の値を求める。
> t <- qt(1-(1-0.95)/2, length(incomes)-1) > t [1] 2.04523
qt() は $t$-分布の Quantile function。
> w <- t * sqrt(var(incomes)) / sqrt(length(incomes)) > w [1] 3.117875
信頼区間の幅 (片側分) を w として計算。
> r <- c(mean(incomes) - w, mean(incomes) + w) > r [1] 51.61546 57.85121
平均 $\pm$ 幅が求める信頼区間。
これを州ごとに行ってみる。数値ベクトルを引数に取ってその信頼区間を返す関数を定義しておいて、それを tapply() で適用すれば良い。
> f <- function(x) {t <- qt(1-(1-0.95)/2, length(x)-1); w <- t * sqrt(var(x)) / sqrt(length(x)); c(mean(x)-w, mean(x)+w)}
> incrange <- tapply(incomes, statef, f)
> incrange
$act
[1] 25.44069 63.55931
$nsw
[1] 46.25363 68.41304
$nt
[1] -1.677921 112.677921
$qld
[1] 42.19966 65.00034
$sa
[1] 46.28451 63.71549
$tas
[1] 54.1469 66.8531
$vic
[1] 41.4402 70.5598
$wa
[1] 43.79253 60.70747
複数のカテゴリを持つ場合にも tapply() は使用できる。性別の因子を付加してみる。
> sex <- factor(c('m','f','m','f','m','f','m','f','m','f','m','f','m','f','m','f','m','f','m','f','m','f','m','f','m','f','m','f','m','f'))
> sex
[1] m f m f m f m f m f m f m f m f m f m f m f m f m f m f m f
Levels: f m
tapply() に複数の因子を適用するには、因子を list にまとめれば良い。
> tapply(incomes, list(statef, sex), mean)
f m
act 44.5 NA
nsw 55.0 58.50
nt 55.5 NA
qld 61.0 51.75
sa 49.0 57.00
tas 61.0 60.00
vic 55.5 58.00
wa 51.0 53.50
この場合、結果は2次元配列で返ってくる。因子の指定の仕方によっては結果が不調和配列になっている可能性がある。
※ R 入門には不調和配列 (Ragged array) と書かれている。通常 Ragged array と言えば、列数の異なる行が集まっている配列のことだが、R の配列はそうはなり得ないように思われる。ここで Ragged array と言っているのは、因子の組み合わせに該当する要素が無いために結果に NA が含まれることを言っているのだと思う。
順序付き因子
順序付き因子を作るには factor() のかわりに ordered() を用いる。順序付き因子はクラス ordered でクラス factor も継承している。