<< 目次へ
最終更新日 ： 2016.9.17
名義・順序データの解析
| 【注意】 本ページの初版（？）は2009年5月13日に作成されたものであり、VGAMパッケージを用いた解析例を説明しています。しかし、現在ではVGAMパッケージが拡張され（より良いものに変わり）p値の計算などデフォルトで出力される仕様になっています。したがって、本ページの内容はもはや不要なものである部分が多いです。ページ修繕を試みましたが、当時の環境を再現できず（R本体もパッケージも新しくなっています）、その手間を考えて本ページの修繕は諦めています。最近では参考になる書籍も豊富になってきているので、そういったものを参照するようにしてください。特にページ後半の 実行例2 以下の部分は不要なものです。 |
名義データの解析 -多項ロジットモデル-
2値データに対して用いたロジスティック回帰分析を拡張して、3つ以上のカテゴリをもつ名義データに対して用いることができる。このようなモデルを多項ロジスティックモデルといい、このモデルの従属変数は多項分布に従うことを仮定している。
| データの種類 | 分布 |
| 2値データ（2つのカテゴリを持つ名義データ） | 従属変数は二項分布に従うと仮定 |
| 3つ以上のカテゴリを持つ名義データ | 従属変数は多項分布に従うと仮定 |
多項ロジスティックモデルをRで実行するにはVGAMパッケージにあるvglm()という関数を用いる。このパッケージは標準では実装されていないので、新たにインストールする必要がある。
# VGAMパッケージのインストール
> install.packages("VGAM")
# VGAMの呼び出し
> library(VGAM)
ちなみに、このレベルの話題になると私自信も勉強しないとよく分からない。また、残念なことに日本語で書かれた参考図書やwebサイトはあまり多くないので、海外の（主に英語の）webサイトなりを自力で調べる必要性も出てくるだろう。
以下、とりあえずRでの実行例をあげておく。
実行例1 文献【1】p190-表8.1の解析例
このようなデータファイルを用意して、適当なディレクトリに保存する（ワークディレクトリに保存するのが良い。ワークディレクトリはgetwd()で確認することができる。）。
# ワークディレクトリの確認
# ここにTable81.datというデータファイルが保存されているとする
> getwd()
[1] "C:/Documents and Settings/Owner/My Documents/R/workdir"
# データの読み込み
> dat <- read.table("Table81.dat")
# 従属変数esaを用意する
> esa <- dat[,2]
> esa
[1] I I I I I F F O I O O I F F F F I I I O F F O F F F F I F F O F I F F F
[37] F F F I F F F O F I I I I F F F F I O I F I F
Levels: F I O
# 独立変数taicyoを用意する
> taicyo <- dat[,1]
> taicyo
[1] 1.24 1.45 1.63 1.78 1.98 2.36 2.79 3.68 1.30 1.45 1.65 1.78 2.03 2.39
[15] 2.84 3.71 1.30 1.47 1.65 1.78 2.03 2.41 3.25 3.89 1.32 1.47 1.65 1.80
[29] 2.16 2.44 3.28 1.32 1.50 1.65 1.80 2.26 2.46 3.33 1.40 1.52 1.68 1.85
[43] 2.31 2.56 3.56 1.42 1.55 1.70 1.88 2.31 2.67 3.58 1.42 1.60 1.73 1.93
[57] 2.36 2.72 3.66
# vglm()を用いてモデル解析
> res <- vglm(esa ~ taicyo, multinomial)
> summary(res) # 係数表の出力
Coefficients:
Value Std. Error t value
(Intercept):1 1.61773 1.30727 1.23749
(Intercept):2 5.69744 1.79374 3.17629
taicyo:1 -0.11011 0.51708 -0.21294
taicyo:2 -2.46545 0.89960 -2.74059
このデータセットはワニの胃袋にどのような餌が入っていたかを調べたものである。従属変数に用いたesa（餌）はF = 魚、I = 無脊椎動物、O = その他という3つのカテゴリをもっている。独立変数taicyo（体長）は連続型である。
要するに
餌 = 体長
というモデルを解析するということは、ワニの体長からどのような餌を食べているかを予測することができるだろう、ということである。
今回の例のように連続型である独立変数が1つのみある場合は、比較的、簡単に結果を表すことができる。以下のような3つの関数を定義しておくと非常に便利である。
体長がxのとき、F（魚）を餌としている確率を求める関数 ：
f1 <- function(x) exp(1.62 - 0.11 * x) / (1 + exp(1.62 - 0.11 * x) + exp(5.70 - 2.47 * x))
体長がxのとき、I（無脊椎動物）を餌としている確率を求める関数 ：
f2 <- function(x) exp(5.70 - 2.47 * x) / (1 + exp(1.62 - 0.11 * x) + exp(5.70 - 2.47 * x))
体長がxのとき、O（その他）を餌としている確率を求める関数 ：
f3 <- function(x) 1 / (1 + exp(1.62 - 0.11 * x) + exp(5.70 - 2.47 * x))
ここで、例えば体長が3mであるワニがF, I, Oそれぞれを摂食している確率を求めてみると次のようになる。
# Fを摂食している確率 > f1(3.0) [1] 0.7598129 # Iを摂食している確率 > f2(3.0) [1] 0.03782886 # Oを摂食している確率 > f3(3.0) [1] 0.2091543
F, I, Oそれぞれの確率を体長の関数として描かれる曲線は次のようになる。
# 魚を摂食する確率
plot(f1, xlim=c(0, 10), ylim=c(0, 1), lwd=2)
# 曲線を描くためのx軸のデータ
x2 <- seq(0, 20, 0.1)
# 無脊椎動物を摂食する確率
I.y <- f2(x2)
lines(x2, I.y, col="red", lwd=2)
# その他を摂食する確率
O.y <- f3(x2)
lines(x2, O.y, col="blue", lwd=2)
# 凡例の追加
legend(6, 0.9, # 凡例を書く座標
c("F", "I", "O"), # 凡例のラベル
col=c("black", "red", "blue"), # 凡例のラインの色
lwd=2 # ラインの太さ
)
簡単に結論を述べると、体長が0mから1.5m位までのワニは無脊椎動物を摂食している確率が高い。体長が3m前後のワニは魚を食べている確率が高いが、4mを越えるワニになると、魚を摂食する確率は徐々に減少していき、逆にその他の餌を食べている確率が増えていくようである。
もっとも、この曲線から判断できるのはせいぜい4mまでで、現実的にそれ以上を越えるワニというのは存在する可能性も低くなるだろうし、何を食べているかも予測するのが困難であると考えるのが妥当であろう。まして8mを越えるワニなど存在するかどうか灰色なものである。
実行例2 文献【1】p193-表8.3の解析例
# 以下のようなデータを用意する必要がある
# Rでは先頭の列を削除(ベースライン=0)として解析するので,
# levelsで適当に水準の順序を変更する必要がある。
# (ただし、どのような順序であっても最終的に得られる結果は同じである)
# Raceは人種を表すカテゴリカル変数
# 教科書ではWhite=1としているので,
# RではBlack=0, White=1という順で水準を指定する。
Race <- rep(c("White", "Black"), c(860, 131))
Race <- factor(Race, levels=c("Black", "White"))
# Sexは性別を表すカテゴリカル変数
Sex <- rep(c("Female", "Male", "Female", "Male"), c(494, 366, 88, 43))
Sex <- factor(Sex, levels=c("Male", "Female"))
# Beliveは死生観を表すカテゴリカル変数
Belive <- rep(c(
"Yes", "Dknow", "No",
"Yes", "Dknow", "No",
"Yes", "Dknow", "No",
"Yes", "Dknow", "No"), c(
371, 49, 74,
250, 45, 71,
64, 9, 15,
25, 5, 13))
Belive <- factor(Belive, levels=c("Yes", "Dknow", "No"))
> library(VGAM)
> result <- vglm(Belive ~ Sex + Race, multinomial)
> summary(result)
Coefficients:
Value Std. Error t value
(Intercept):1 0.88305 0.24264 3.63939
(Intercept):2 -0.75801 0.36135 -2.09769
SexFemale:1 0.41855 0.17125 2.44407
SexFemale:2 0.10506 0.24651 0.42621
RaceWhite:1 0.34177 0.23703 1.44190
RaceWhite:2 0.27098 0.35412 0.76520
# どんなデータを解析しているのか確認してみる。
# データフレームmydataを作成
> mydata <- data.frame(
+ R = Race, S = Sex, B = Belive)
# 膨大な行数になるので,コンソール画面上にではなくテキストファイルへ出力
# ワークディレクトリにTable8_3_dataframe.txtというファイルが生成される
> sink("Table8_3_dataframe.txt")
> mydata # このデータフレームは先ほどのテキストファイルに出力される
> sink() # 出力先をコンソール画面に戻す
# fitted()かpredict( , type="response")で死生観の予測確率を求めることができる
> a <- fitted(result)
> b <- predict(result, type="response")
# これまた膨大な量になるので適当なテキストファイルに書き出したほうが無難
> sink("a_out.txt")
> a
> sink()
> sink("b_out.txt")
> b
> sink()
*** *** *** *** ***
vglm()では、summary()によって得られる係数表にp値が表示されない。それなので、以下のような関数を用意しておくと便利である。
# --- ここから --- toridashi <- function(res.dat, n){ p <- length(res.dat$t.value) # 推定したパラメータ数 DF <- (n - p) # 自由度 t.value <- res.dat[,3] # t値を代入 p.value <- c() # p値を格納する変数 # p値を計算していく for(i in 1:p){ p.value[i] <- (1 - pt(abs(t.value[i]), DF)) * 2 } # 係数表にp値のデータベクトルを結合させて出力 res <- cbind(res.dat, p.value) cat("Coefficents: \n") res } # --- ここまで ---
使い方は少し面倒なのだが、まずsummary()によって得られる係数表（Coefficents: という部分）をまるごとコピーして適当なテキストエディタに貼り付ける。つづいて行ラベルと列ラベルの文字列をダブルクォーテーション" "でくくる（下のスクリーンショットのように）。それができたらvglm.txt（としてRのワークディレクトリに保存する。なお、ファイル名と保存先はべつにどこでもよい。これで前準備は完了である。
 |
上述したようにテキストファイルを用意することができたら、データファイルをread.table()で読み込む。
# データファイルの読み込み
> num <- read.table("vglm.txt")
# 中身を確認してみる
> num
Value.Std. Error t.value
(Intercept):1 0.88305 0.24264 3.63939
(Intercept):2 -0.75801 0.36135 -2.09769
SexFemale:1 0.41855 0.17125 2.44407
SexFemale:2 0.10506 0.24651 0.42621
RaceWhite:1 0.34177 0.23703 1.44190
RaceWhite:2 0.27098 0.35412 0.76520
# まっ、これは特に気にしなくてよい^^;
> names(num)
[1] "Value.Std." "Error" "t.value"
> rownames(num)
[1] "(Intercept):1" "(Intercept):2" "SexFemale:1" "SexFemale:2" "RaceWhite:1"
[6] "RaceWhite:2"
それから、自作した関数toridashi()を使う。toridashi()の第1引数は係数表が格納された変数（この場合はnumである）、第2引数はサンプルサイズである。サンプルサイズはlength()で確認すればよい。
# 目的変数Beliveをlength()で調べればサンプルサイズが分かるでしょ？
> length(Belive)
[1] 991
# 第1引数: 係数表が格納されているデータフレーム
# 第2引数: サンプルサイズ
> toridashi(num, 991)
Coefficents:
Value.Std. Error t.value p.value
(Intercept):1 0.88305 0.24264 3.63939 0.0002874843
(Intercept):2 -0.75801 0.36135 -2.09769 0.0361868772
SexFemale:1 0.41855 0.17125 2.44407 0.0146971181
SexFemale:2 0.10506 0.24651 0.42621 0.6700479525
RaceWhite:1 0.34177 0.23703 1.44190 0.1496484449
RaceWhite:2 0.27098 0.35412 0.76520 0.4443357576
*** *** *** *** ***
説明変数であるRace(Levels: Black, White)とSex（Levels: Male, Female）の各水準の組み合わせごとの予測確率を求める。これはすでに記したように、関数fitted()を用いれば求められる。
> dat.eins <- data.frame(R = Race, S = Sex) > a <- fitted(result) > dat.zwei <- data.frame(dat.eins, a) > names(dat.zwei) [1] "R" "S" "Yes" "Dknow" "No"
ここに、dat.zweiには4つの変数が存在している。しかし、これらの中身を表示させると991行にもわたる大量のデータが表示されてしまうので、Rのコンソール画面に表示させるのは控えた方がよい。かといって、sink()を利用してテキストファイルに出力しても、結局はごちゃごちゃして見にくいのには変わりがない。そこで、こんな風にすると比較的コンパクトにまとめられる。
# 白人かつ女性の予測確率
> apply(dat.zwei[dat.zwei$R=="White" & dat.zwei$S=="Female", c(3, 4, 5)], 2, mean)
Yes Dknow No
0.75456083 0.09956287 0.14587630
# 白人かつ男性
> apply(dat.zwei[dat.zwei$R=="White" & dat.zwei$S=="Male", c(3, 4, 5)], 2, mean)
Yes Dknow No
0.6782703 0.1224479 0.1992817
# 黒人かつ女性
> apply(dat.zwei[dat.zwei$R=="Black" & dat.zwei$S=="Female", c(3, 4, 5)], 2, mean)
Yes Dknow No
0.7073517 0.1001812 0.1924671
# 黒人かつ男性
> apply(dat.zwei[dat.zwei$R=="Black" & dat.zwei$S=="Male", c(3, 4, 5)], 2, mean)
Yes Dknow No
0.6221640 0.1205594 0.2572766
表にまとめてみると分かりやすい。これは文献【1】のp295の表8.5と同じものであることが確認できるだろう（小数点を切り捨てているので若干異なるが）。
Race Sex 信じる 分からない 信じない 白人 女性 0.75 0.09 0.14 白人 男性 0.67 0.12 0.19 黒人 女性 0.70 0.10 0.19 黒人 男性 0.62 0.12 0.25
どうやら、上の2つの問題をクリアできたようだ。
順序データの解析 -累積ロジットモデル-
vglm(Y ~ X, cumlative)とすることでできそうな感じ。
http://www.stat.cmu.edu/~brian/720-2004/R/2004-11-16-notes.txt
---
最後に従属変数が順序データである場合の解析法を紹介する。例として文献【2】のp168-表8.1のデータを用いて分析してみよう。まずはRでの解析例を以下に提示するが、従属変数である順序データはカテゴリカル型であることに注意してもらいたい。
# まずは分割表をもとにデータフレームを作成する
# 性別のデータ
Sex <- rep(c("female", "male"), c(150, 150))
# 年齢のデータ
Age <- rep(c("18-23", "24-40", ">40",
"18-23", "24-40", ">40"),
c(45, 45, 60,
65, 44, 41))
# 反応のデータ(順序データ)
Reaction <- rep(c(
"unimportant", "normal", "important",
"unimportant", "normal", "important",
"unimportant", "normal", "important",
"unimportant", "normal", "important",
"unimportant", "normal", "important",
"unimportant", "normal", "important"),
c(26, 12, 7,
9, 21, 15,
5, 14, 41,
40, 17, 8,
17, 15, 12,
8, 15, 18))
# カテゴリカル型であることを認識させる
# 関数factor()を用いて、各変数のカテゴリの優先順位を指定する
# 順序データである場合は引数にordered=TRUEを追加する
Sex <- factor(Sex, levels=c("female", "male"))
Age <- factor(Age, levels=c("18-23", "24-40", ">40"), ordered=TRUE)
Reaction <- factor(Reaction, levels=c("unimportant", "normal", "important"), ordered=TRUE)
# データフレームの作成
mydata <- data.frame(SEX = Sex, AGE = Age, REACTION = Reaction)
# 作成したデータフレームをテキストファイルに出力する
# 行数が多いので、Rのコンソール画面で表示させるのは避けた方がよい
sink("Table_8_1.txt")
mydata
sink()
# vglm()を使ってモデル解析をする
# familyの指定に注意すること
> result2 <- vglm(REACTION ~ SEX + AGE, family=cumulative(parallel=T), data=mydata)
> summary(result2)
Coefficients:
Value Std. Error t value
(Intercept):1 0.043533 0.23030 0.18903
(Intercept):2 1.654972 0.25359 6.52617
SEXmale 0.576223 0.22611 2.54844
AGE24-40 -1.147094 0.27727 -4.13712
AGE>40 -2.232452 0.29041 -7.68728
Log-likelihood: -290.6478 on 595 degrees of freedom
Rの出力結果について少し補足をしておこう。Coefficients（係数表）において、文献【2】のp177-表8.4に記されている係数とは符号が逆転しているが、結果としては同じものである（教科書と同じ結果を望むのであれば、Rによって出力された係数の符号をすべて逆転させればよいだけのことである）。
また、教科書でいうところの最大対数尤度にあたるものが"Log-likelihood"となっており、この値は-190.6478といったように、教科書の値と一致しているのが確認できるだろう。
ところで、「性別が女性で、年齢が18-23歳」である場合の重要度（"unimportant", "normal", "important"）の確率は、関数fitted()を利用することによって簡単に得られる。ただし、これも出力される結果の行数が多い（300行）ので、sink()を使ってテキストファイルに出力した方がよい。見やすい結果の表示の仕方はこのページの一番下を参照。
> sink("ore.txt")
> fitted(result2)
> sink()
このようにすることによって、Rのワークディレクトリ（関数getwd()で確認できる）にore.txtというテキストファイルが作成されているはずである。実際に確認してみると良いのだが、最初に作成したTable8_1.txtと並べて見ると分かりやすい。例えば、下の2つのスクリーンショット（Windowsではプリントスクリーン機能を利用できる。ノートパソコンなら[Fn] + [PrtSc]でスクリーンショットを撮ることができる。）
 SEX = female, AGE = 18-23である場合 ：
|
 SEX = male, AGE = 24-40である場合 ：
|
cum.dat <- data.frame(S = Sex, A = Age)
cum.res <- fitted(result2)
matome <- data.frame(cum.dat, cum.res)
sink("matome.txt")
matome
sink()
attach(matome)
levels(A)
levels(S)
mean(matome[matome$S=="female" & matome$A=="18-23", c(3, 4, 5)])
mean(matome[matome$S=="female" & matome$A=="24-40", c(3, 4, 5)])
mean(matome[matome$S=="female" & matome$A==">40", c(3, 4, 5)])
mean(matome[matome$S=="male" & matome$A=="18-23", c(3, 4, 5)])
mean(matome[matome$S=="male" & matome$A=="24-40", c(3, 4, 5)])
mean(matome[matome$S=="male" & matome$A==">40", c(3, 4, 5)])
p値を得るには先ほどやったようにデータファイルを用意して以下のようにする。
> dat <- read.table("cum.txt")
> toridashi(dat, n)
Coefficents:
Value Std..Error t.value p.value
(Intercept):1 -1.082982 0.18032 -6.00601 5.584156e-09
(Intercept):2 0.528456 0.17138 3.08354 2.239125e-03
SEXmale 0.576223 0.22611 2.54844 1.132732e-02
AGE.L -1.578582 0.20535 -7.68728 2.249312e-13
AGE.Q 0.025203 0.19504 0.12922 8.972716e-01
問題点
多項ロジットモデルにおいて挙げた2つの問題点に加え、今回新たにもう1つの問題点が見つかった。
→ parallel=TRUEにすると比例オッズモデル（proportional odds model）として扱われます。それに対してparallel=FALSEとすると累積ロジットモデル（cumulative logit model）として扱われます。これでまた新たに問題が1つ解決。しかも、何やらよさそうな資料を見つけたので、さらにvglm()を使いこなせそう。しばらくしたら、また新たな解説を加える。
参考文献・参考ウェブ資料
【1】 Alan Agresti著, 渡邉ほか訳「カテゴリカルデータ解析入門」サイエンティスト社（2005）
【2】 Annette J. Dobson著, 田中ほか訳「一般化線形モデル入門」共立出版（2008）
【3】 http://support.sas.com/techsup/technote/ts722g.pdf SASのサポートページ（PDFファイル）
【4】 http://mjin.doshisha.ac.jp/R/200808_61.pdf 同志社大学の金先生の資料（PDFファイル）
【5】 http://www.ats.ucla.edu/stat/R/dae/mlogit.htm Rのvglm()の使用例
【6】 http://www.ats.ucla.edu/stat/sas/output/SAS_mlogit_output.htm SASで多項ロジットモデルを扱った例
