Rで混合モデル (R Advent Calendar 2013 16日目)
R Advent Calendar 2013の16日目です．
「Rと同位体分析: siarパッケージの紹介」の内容を文章にして，補足情報を載せました．
背景
生態学の分野では，ある生物がさまざまな食資源をそれぞれどのくらいの割合で摂取しているかを知ることが，重要な研究課題のひとつになっています．この研究課題に対して，ここ数十年のあいだに，同位体分析と混合モデル (mixing model) によるアプローチがよく行なわれるようになってきました．同位体分析と混合モデル
同位体分析についてすごーく簡単に述べます．ある生物の体組織の同位体比を測れば，その生物の摂取してきた食物の内容がおおまかにわかります．同位体比というのは，質量数の異なる同一元素の存在比のことで，炭素の場合は，具体例を出すと，以下のようになります．アマモとアオサのみを食べているガチョウがいたとして (図1)，それらの食資源とガチョウ体組織の炭素同位体比 (値) を測って連立方程式を解けば，食物全体に占めるアオサ，アマモの相対的な寄与率を求めることができます．
| 同位体比の測定結果 | δ13C (‰) |
|---|---|
| ガチョウ (同位体分別補正済み) | -13.2 |
| アマモ | -11.2 |
| アオサ | -14.1 |
図1 ガチョウとその食資源の炭素同位体比
δ13C{ガチョウ} = a δ13C{アマモ} + b δ13C{アオサ}
a + b = 1
同位体比測定結果をもとににこの連立方程式を解くと，以下のようになります．
アマモ寄与率 a = 0.31
アオサ寄与率 b = 0.69
用いる同位体の指標を2つに増やせば，3つの食資源の寄与率を計算できます．先ほどの計算に，新たに窒素同位体比 () の結果と食資源の草を加えて，アオサ・アマモ・草のみを食べているガチョウを考えてみます (図2)．
| 同位体比の測定結果 | δ13C (‰) | δ15N (‰) |
|---|---|---|
| ガチョウ (同位体分別補正済み) | -13.2 | 6.8 |
| アマモ | -11.2 | 6.5 |
| アオサ | -14.1 | 9.8 |
| 草 | -30.9 | 4.4 |
図2 ガチョウとその食資源の炭素・窒素同位体比
δ13C{ガチョウ} = a δ13C{アマモ} + b δ13C{アオサ} + c δ13C{草}
δ15N{ガチョウ} = a δ15N{アマモ} + b δ15N{アオサ} + c δ15N{草}
a + b + c = 1
連立方程式を解くと，以下のようになります．
アマモ寄与率 a = 0.81
アオサ寄与率 b = 0.13
草寄与率 c = 0.06
混合モデルとは，こうした，食資源と消費者の同位体比を元に，それぞれの食資源の寄与率を計算するモデルのことです．どうしてこんな簡単な計算にモデルが必要かというと，以下3つの大きな問題があるためです．
- 利用可能な同位体比の指標 (せいぜい2-4種類) に対して，摂取割合を推定したい食資源の数が多い場合，連立方程式を解くことができなくなる．
- 生物現象のため，食資源や消費者体組織の値にばらつきがあり，平均値で代表させた値の計算では不正確なことがある．
- 上の例では示していないが，食資源の元素濃度や，消費者体内での同位体の扱いの違いなど，さらに考慮しなければならない変数がある．
こうした問題を解決するため，モンテカルロ法による解析，総あたり計算，濃度依存混合モデルなど，さまざまな解決案が提示されてきました．このあたりの事情は，Phillips (2012) をご参照ください．
ベイズ推定を利用した混合モデル"siar"
そんななか，2010年に，ベイズ推定を利用した強力なモデルが，Parnell et al. (2010) によって提案されました．このモデルは，RパッケージとしてCRANに公開され，Google Scholarで調べてみると，論文の引用数はすでに320を超えています (2013年12月16日現在: 同位体生態学の論文では多いほうだと思います)．このパッケージsiar (Stable Isotope Analysis in R) の紹介をするのが本エントリの目的です (やっとRの話が出てきました！)．ユーザインタフェースがつくりこまれており，ふだんあまりRに触れない方でも，けっこう簡単に使えてしまいます．
特徴は以下の通りです．
- (用いる同位体比の指標数 + 1) 以上の食資源からの寄与率であっても計算できる．
- ベイズ推定のフレームワークを採用しているので，値のばらつきをモデルに組み込み，各食資源の寄与率を事後分布として推定できる．事前分布も設定可能．
- 食資源の元素濃度を組み込んだモデルや，消費者体内での同位体の扱いの変動などにも対応可能．
モデルの中身については，元の論文 (Parnell et al. 2010) をご覧ください．仕組みは非常に簡単なので，同位体生態学以外の分野にも応用が効きそうです．
デモ
以下に，様子を載せておきます．# まずインストール
install.packages("siar")
# パッケージの読み込み
library(siar)
# メインの関数を開いてみましょう
siarmenu()そうすると，以下のような選択式のメニューがあらわれ，あとはデータを読み込んで，数字を選んでいくだけで，ベイズ推定を組み込んだ混合モデルを利用したひと通りの解析ができます．
------------------------------- Welcome to Stable Isotope Analysis in R version 4.1 Author: Andrew Parnell, University College Dublin Please report bugs to: Andrew.Parnell@ucd.ie ------------------------------- Useful: Press 0 at a prompt to return to the main menu or Esc to exit. The available options are: 1: Load in some data (run this first) 2: Run SIAR for a single group 3: Run SIAR for multiple groups 4: Plot single group proportions 5: Matrix plot of proportions 6: Plot of proportions by source 7: Save parameter output to a file 8: Summary information and convergence diagnostics 9: Demo (for first time users) 10: Exit Selection:
はじめは9のデモを選んでみると，一連の流れがよくわかります．上記に示したガチョウの食性データはこのデモから引用したものです．
関数群
自分好みに変数や設定を変えたりしたい場合は，ヘルプと，以下のガイドが参考になります．パッケージsiarの作者らが作成した文書で，環境中に読み込まれる関数などについて説明しています．これらの関数を使えば，いちいち，関数siarmenuによるユーザインタフェースを通して解析する必要がなくなり，通常ようなコマンドラインからの操作ができます．SIAR V4 (Stable Isotope Analysis in R) An Ecologist's Guide (PDF)
たとえば，デモで見たのと同じ操作が以下のコマンドで可能です．
# データ読み込み data(geese1demo) data(sourcesdemo) data(correctionsdemo) data(concdepdemo) # MCMCを走らせる siarresult <- siarmcmcdirichletv4( data = geese1demo, sources = sourcesdemo, corrections = correctionsdemo, concdep = concdepdemo) # 生データのプロット (図3) siarplotdata(siarresult) # 混合モデル解析結果のプロット (図4) siarhistograms(siarresult)
図3 ガチョウとその食資源同位体比のプロット
図4 各食資源の寄与率推定結果
まとめ
混合モデルは，同位体生態学の研究以外にも応用可能です．「何か異なるものが混ざり合ってひとつのものができていて，元の混合比率を知りたい」という課題に対して，ここで紹介したパッケージsiarは強力な解法を提供してくれます．なにか機会があったら，ぜひ使ってみてください．引用文献
いずれもオープンアクセスです．Parnell AC, Inger R, Bearhop S, and Jackson AL. 2010. Source partitioning using stable isotopes: coping with too much variation. PLoS One 5:e9672.
Phillips DL. 2012. Converting isotope values to diet composition: the use of mixing models. J Mammal 93:342-352.
統計ソフトの言及数推移
学術文献における統計解析ソフトの言及数を，最近の約10年間で年ごとに調べてみました．
方法
Google Scholarを利用して，以下の検索クエリを含む学術文献数 (概算) を，2000年から2012年までの期間で取得した．得られた数値は，R，Bioconductor，SPSS，Excel にそれぞれ言及した学術文献の出版数とみなした．- "R Development Core Team"
- "Bioconductor"
- "SPSS"
- "Excel"
結果
結果は表と図に示す通り (図は対数表示)．RとBioconductorの言及論文数が増加しており，SPSSとExcelはほぼ横ばいまたは微減であった．表 各統計ソフト言及学術文献数の年推移
| 年 | R | Bioconductor | SPSS | Excel |
|---|---|---|---|---|
| 2012 | 17500 | 4060 | 72800 | 88500 |
| 2011 | 13400 | 3380 | 66500 | 85500 |
| 2010 | 9480 | 2900 | 91300 | 121000 |
| 2009 | 7090 | 2410 | 126000 | 151000 |
| 2008 | 5280 | 2090 | 143000 | 165000 |
| 2007 | 3570 | 1560 | 161000 | 171000 |
| 2006 | 2210 | 1190 | 166000 | 175000 |
| 2005 | 1210 | 971 | 167000 | 159000 |
| 2004 | 354 | 385 | 149000 | 154000 |
| 2003 | 91 | 162 | 125000 | 130000 |
| 2002 | 42 | 36 | 96000 | 107000 |
| 2001 | 24 | 17 | 71100 | 85000 |
| 2000 | 7 | 8 | 43000 | 63400 |
図 各統計ソフト言及学術文献数の年推移
考察
結果より，RやBioconductorを利用した，もしくは視野に入れた学術文献の絶対数は，ここ10年間程度で着実に大きく増加してきていることが強く示唆される．ただし，以下に示すような留意点もある．
- Google Scholarにおけるヒット件数が言及学術論文数とイコールになるとは限らない．特にRの検索クエリはRを「正しく」引用するときに用いる一連の語句だが，他3者はその統計ソフトウェアの名前自体をクエリとしているため，Rの言及数は過小評価されている可能性がある．
- 学術文献全体の出版数が不明であるため，相対的な出版数については増減の判断がつかない．
参考
Google Scholar本エントリはKashiwa.R#7で発表した内容を元にしています (Kashiwa.R#7特設サイト)．
Kashiwa.Rについて: #7へのイントロダクション (Kahiwa.R#7発表スライド)
irisの正体 (R Advent Calendar 2012 6日目)
R Advent Calendar 2012の6日目です．
Rのデータセットでもっとも有名なのが，irisではないでしょうか．FisherもしくはAndersonのアヤメの計測データですね．しかし，有名なわりには，このデータセットの正体はそれほどよく知られていないように思います (私もよく知らずに使っていました)．
そこで，このエントリでは，どんな研究でirisデータセットが報告されたのか，元の論文2報について，その概要を述べてみたいと思います ※1．
アヤメ
Iris属 (アヤメ属) はアヤメ科に含まれる属のうちのひとつで，世界の温帯に150種が知られているようです (アヤメ属 - Wikipedia)．irisデータセットに含まれるIris setosa，Iris versicolor，Iris virginicaは，このアヤメ属に含まれる植物のうちの3種です．どんな論文で報告されたの？
米国の植物学者Edgar Andersonによって1936年に発表された研究のなかで，もともとの計測が行なわれています (Anderson, 1936 ※2)．数値データが報告されたのは，英国の統計学者Ronald Aylmer Fisherの1936年の論文です (Fisher, 1936 ※3)．Andersonの研究
アヤメ属の3種 (Iris setosa，Iris versicolor，Iris virginica) の種分化について，形態，遺伝，生態，進化の考察を行なった大論文です．論旨をごく簡単にまとめると，以下になります．
- setosaに比べて，versicolorとvirginicaの形態はよく似ている．
- setosaとの違いは，virsinicaのほうが大きく，versicolorで小さい．
- versicolorは，setosaとvirginicaの雑種に由来しているかもしれない．
以下に，論文中のいろんなトピックを述べてみました．しかし，生物学の専門用語なんかもでてきてしまったので，難しいようでしたら読み飛ばしてください．
形態
大きさや形の特徴を比べています．全体的な形態はよく似ているが，sepal (がく片) とpetal (花弁) の長さ・幅が，特に種間で異なる特徴である，と述べられています ※4．分布
以下のようになっているそうです．- setosaは亜北極帯に，virginicaはもっと南に分布し，生息域は重なっていない．
- setosaとversicolorの分布域は重なるが，交配は稀．
- virsicolarとvirginicaの分布域が重なるところでは交配が確認されている．
核相
染色体数は，setosaで38本，virginicaで70-72本，vergicolarで108本であり，virsinicaは二倍体で，versicolorはsetosaとvirsinicaの複二倍体ではないかと考察しています．まあつまり，{setosaの38 + virginicaの70 = versicolorの108} ということを言いたいわけですね．
種分化のプロセス
アヤメ属はもともと雑種形成がさかんな種で，コロニーをつくるその生息様式と関連して，氷期 (最終氷期？) の前か間氷期くらいの昔に，種分化が起こったのではないか，と考察しているようです．sepal (がく片) とpetal (花弁) ※4
これらは，3種の分類に利用できる形態的な特徴として，論文中でも中心的なデータとして提示されています．まず，sepalとpetalの図です．右からsetosa，versicolor，virginicaで，sepalは左側，petalは右側です．
irisのデータは，これらの長さと幅を測った結果だったわけですね ※5．
以下のように，Andersonは，長さ・幅に対応する長方形を作成し，sepalとpetalでの差分を示すことで，種間の違いを視覚的に提示しています．
ついでに
謝辞には，S Wright，JBS Haldane，RA Fisherといった集団遺伝学の大家が名を連ねています．この研究が行なわれた頃，Andersonは英国の研究所の特別研究員に招かれ，HaldaneやFisherと共同研究をしていたようです (Edgar Anderson - Wikipedia)．こんな経緯があり，Fisher (1963) でAndersonの計測データが使われることになったりしているのでしょうか．Fisherの研究多変数の線形判別分析の理論を提示したあと，Andersonの計測したアヤメsepal・petalの長さ・幅データにその方法を適用しています．
判別分析に関しては，Rと判別分析 (フリーソフトによるデータ解析・マイニング 第17回) などが参考になるかと思います．
結果
- setosaに比べて，versicolorとvirginicaはよく似ており，判別がときに困難である．
- setosaとの違いは，virsinicaのほうが大きく，versicolorで小さい．
ついでに
発表された雑誌が"Annals of Eugenics" (『優生学紀要』) だったりするあたり，時代を感じますね…．まとめ
Anderson (1936) の趣旨は，以下の2点でした．- 形態的に，Iris versicolorとIris virginicaは似ている一方，Iris setosaとは異なっており，違いはIris virginicaのほうでより大きい．
- 形態，分布，遺伝，生態，進化に関する証拠より，Iris versicolorは，Iris setosaとIris virginicaの雑種 (複二倍体) に由来するかもしれない．
Fisher (1936) は，Anderson (1936) で計測されたデータを利用して，多変数の線形判別分析を行ないました．
ということで，本エントリの内容は以上なのですが，あまりにRと関係ない話題なので，ちょっと心が痛んできました．ですので，申し訳程度にコードを載せておきます．どうぞ参考になさってください．
?iris
注
※1 そうは言ったものの，植物系統分類学は学部の授業で習ったくらいで，論文もざっと読んだ程度ですし，最新の分子生物学がアヤメの種分化に関する知見をどのくらいアップデートさせたかも知りませんので，もしお詳しい方いらっしゃいましたら，つっこみやコメントをくださいましたら幸いです．※2 Anderson E. 1936. The Species Problem in Iris. Annals of the Missouri Botanical Garden 23:457-509.
http://biostor.org/reference/11559
(CC BY-NCで公開されています．フォントがなんだか優雅な論文です)
※3 Fisher RA. 1936. The use of multiple measurements in taxonomic problems. Annals of Eugenics 7:179-188.
http://onlinelibrary.wiley.com/doi/10.1111/j.1469-1809.1936.tb02137.x/abstract
※4 アヤメの場合，正確には「がく片」「花弁」はそれぞれ「外花被片」「内花被片」と呼ぶそうですが (アヤメ - 国立科学博物館)，ここではややこしくなるので単純に「がく片」「花弁」と呼びました．
※5 実際の花ではこんなふうになっているようです → ヒオウギアヤメ（檜扇文目） - 群馬大 青木先生のサイト
重複のため消しました
以下のエントリと記事がかぶっていたので，消しました．(2012.12.06)
http://d.hatena.ne.jp/tsutatsutatsuta/20121206/1354737461
Rで密度推定
ヒストグラムとカーネル密度推定についてすごく簡単にまとめました．Kashiwa.R#4で発表する内容を文章にしたものです．
使用するデータセット
まずここでは，faithful データセットの，eruption を例に用います．イエローストーン国立公園にあるOld Faithful間欠泉の噴出時間 (分) だそうです．
# データセットの確認 head(faithful) # 噴出時間データのみを抜き出す fe <- faithful$eruption
まず度数分布を見てみると以下のようになります．
stem(fe)
The decimal point is 1 digit(s) to the left of the | 16 | 070355555588 18 | 000022233333335577777777888822335777888 20 | 00002223378800035778 22 | 0002335578023578 24 | 00228 26 | 23 28 | 080 30 | 7 32 | 2337 34 | 250077 36 | 0000823577 38 | 2333335582225577 40 | 0000003357788888002233555577778 42 | 03335555778800233333555577778 44 | 02222335557780000000023333357778888 46 | 0000233357700000023578 48 | 00000022335800333 50 | 0370 >
ヒストグラム
これをヒストグラム (相対度数) で表すのですが，区切り幅や端点の選択によって，形が大きく変化してしまいます．# Sturges (1926) の方法にしたがった区切り幅 (デフォルト) hist(fe, freq = FALSE)
# Scott (1992) の方法にしたがった区切り幅 # KernSmoothパッケージのインストールが必要です library(KernSmooth) fe.h <- dpih(fe) fe.bins <- seq(min(fe), max(fe) + fe.h, by = fe.h) hist(fe, freq = FALSE, breaks = fe.bins)
# 区切り幅の1/2だけ端点をずらしてみる fe.bins.2 <- seq(min(fe) - fe.h / 2, max(fe) + fe.h * 3/2, by = fe.h) hist(fe, freq = FALSE, breaks = fe.bins.2)
困りました…
なにが問題かというと，階級内に含まれるデータの分布を考慮していない点です．
たとえば，c(0.2, 0.6, 0.8) という3つの標本があったとき，示した図はどれも，データを表すヒストグラムになり得てしまうわけです．
data <- c(0.2, 0.7, 0.9) oldpar <- par(mfrow = c(1, 3)) hist(data, breaks = seq(0, 1, 1/2)) points(data, rep(-0, 3), pch = 16) hist(data, breaks = seq(0, 1, 1/3)) points(data, rep(-0, 3), pch = 16) hist(data, breaks = seq(0, 1, 1/6)) points(data, rep(-0, 3), pch = 16)
カーネル密度推定
カーネル密度推定では，この問題を回避して，境界の設定に関係なく分布を推定できるようにするための方法です，ひとつひとつの標本に，カーネル関数から決定された分布を与えて，それを積み上げて全体の分布とするようなイメージです．カーネル密度推定 > 定義 - Wikipedia
カーネル密度推定 > 直感的説明 - Wikipedia
Rでは，density()関数を使えば，簡単にカーネル密度推定ができます．
fe.d <- density(fe) hist(fe, freq = FALSE, breaks = fe.bins) lines(fe.d)
ただし，個々の標本に与える分布の形 (カーネル関数) とその裾野 (バンド幅) をどう設定するかによって，推定の結果が変わってきます．デフォルトでは，カーネル関数に正規分布の確率密度関数，バンド幅に bw.nrd0 関数で計算された値が使用されますが，たとえば以下のように，バンド幅を変えると，滑らかさ (平滑化の程度) が変化します．
fe.d <- density(fe) fe.d.1 <- density(fe, bw = 0.1) fe.d.2 <- density(fe, bw = 1) hist(fe, freq = FALSE, breaks = fe.bins) lines(fe.d) lines(fe.d.1, col = "red") lines(fe.d.2, col = "blue")
と，このようなわけで，絶対的な解はないので気をつけましょうね，というところで終わってしまいます．はい．
ヒストグラムの平滑化
まずヒストグラムを作成して，それをノンパラメトリック回帰を利用して平滑化する，という方法もあるそうです．『Rによるノンパラメトリック回帰の入門講義』 (竹澤邦夫. 2009. メタブレーン) という本に詳しく載っています．
ちょっと時間がないので，これは別の機会に．
カテゴリカルデータの解析 (その2)
その1からのつづきです．
変数間の関係性を調べる
集計したカテゴリカルデータを解析する方法を紹介します．独立性の検定
「変数のあいだに関連性があるか」を調べる検定です．先に紹介したArthritisで，処置と改善度合いのあいだに関連性があるかを調べてみましょう．この場合には，カイ二乗検定 (データ数が少ない場合にはフィッシャーの直接確率検定のほうがよい) を使います．
帰無仮説「処置と改善度合いに関連性はない」を，カイ二乗検定で検討してみます．
> chisq.test(arthritis.imp.tre) Pearson's Chi-squared test data: arthritis.imp.tre X-squared = 13.055, df = 2, p-value = 0.001463
結果，p値 < 0.0015で，対立仮説「処置と改善度合いには関連性がある」をとることになります．
カイ二乗検定の用途には，他に，適合度の検定があり，「観察値と理論値が同じか」を調べることができます．
カイ二乗検定の原理については，以下などがわかりやすいかもしれません．
『カテゴリカルデータ解析 (Rで学ぶデータサイエンス1)』には，以下に示したように，ほかにもいろいろな場合に使う検定法が記されていますので，参考にしてみてください．
- マンテル検定: カテゴリー間の順序関係をスコアに変換/li>
- ウィルコクソン順位和検定: スコアではなく順位の状態で検定
- コクラン・アーミテージ検定: ロジスティック回帰分析のように順序カテゴリーを対象とする
- クラスカル・ワリス検定: 3群以上の比較に拡張したウィルコクソン検定
- マクネマー検定: 変数間に対応関係がある場合に使う
回帰分析
ひとつの変数を目的変数として，目的変数を，その他の変数で説明するかたちのモデルを用いて分析を行なう手法です．これについては以前まとめました．
ほかにも，目的変数が個数データのときに用いるポアソン回帰モデルや，目的変数を特に定めないで変数の関係性を調べるときに用いる対数線形モデルがありますが，ここでは大胆に省略します:(
対応分析
カテゴリカル変数間の関係をうまく表そうとする方法．複数変数のカテゴリーに，相関が高くなるような数値を割り当てて，関係を分析します．たとえば，DanishWelfareのデータで，婚姻状態と居住地域の関係を表してみましょう．
> danish.sta.urb <- xtabs(Freq ~ Status + Urban, data = DanishWelfare)
> round(prop.table(danish.sta.urb, margin = 2), 3)
Urban
Status Copenhagen SubCopenhagen LargeCity City Country
Widow 0.210 0.148 0.093 0.121 0.064
Married 0.522 0.700 0.705 0.705 0.773
Unmarried 0.268 0.151 0.202 0.173 0.163
> chisq.test(danish.sta.urb)
Pearson's Chi-squared test
data: danish.sta.urb
X-squared = 158.1145, df = 8, p-value < 2.2e-16独立性に関するカイ二乗検定を行なうと，p値 < 0.001で，これらふたつの変数間に何らかの関係性がありそうです．
対応分析では，それぞれのカテゴリーに数値を割り当ててみて，変数どうしの相関が高くなる最適な数値の組み合わせを探索します．
最適化問題はcorresp()関数で解くことができます．それぞれのカテゴリーを2次元の値で表現するように指定 (nf = 2) して，結果のベクトルを2次元平面上に表してみましょう．
> danish.corresp <- corresp(danish.sta.urb, nf = 2) > plot(danish.corresp)
割り当てられた数値が近いほど関連が近いと考えられ，平面上でもカテゴリー名が近くに表示されます．Copenhagen以外の地域は比較的まとまっており，これらは同じような傾向を示していることがみてとれます．また，この4地点の近くにMarriedが表示されており，これらの地域では婚姻者の割合が大きこともわかります．
このようにして，変数のカテゴリー間の対応関係を調べる方法が対応分析です．
決定木
説明変数を利用して各個体を段階的に分類していき，どのカテゴリーに属するかを決定する方法です．フィッシャーのアヤメのデータを，決定木を使って解析してみます．mvpartまたはrpartパッケージを使います．
> install.packages("mvpart")
> library(mvpart)
> (iris.part <- rpart(Species ~ ., data = iris))
n= 150
node), split, n, loss, yval, (yprob)
* denotes terminal node
1) root 150 100 setosa (0.33333333 0.33333333 0.33333333)
2) Petal.Length< 2.45 50 0 setosa (1.00000000 0.00000000 0.00000000) *
3) Petal.Length>=2.45 100 50 versicolor (0.00000000 0.50000000 0.50000000)
6) Petal.Width< 1.75 54 5 versicolor (0.00000000 0.90740741 0.09259259)
12) Petal.Length< 4.95 48 1 versicolor (0.00000000 0.97916667 0.02083333) *
13) Petal.Length>=4.95 6 2 virginica (0.00000000 0.33333333 0.66666667) *
7) Petal.Width>=1.75 46 1 virginica (0.00000000 0.02173913 0.97826087) *結果の意味は以下のとおりになります．
最初のnは，対象としたデータの行数を表しています．
次の行は，出力値の意味を示しています．nodeは分岐する点の順番，splitはデータの分け方，nはそのグループに用いた行数，lossは予想が外れた数，yvalはグループの予想値，(yprob)はそれぞれのカテゴリーの割合 (アルファベット順に配置) を表すそうです．*は，分岐の終末点を表します．
このままではわかりづらいので図にしましょう．
> plot(iris.part) > text(iris.part, use.n = TRUE, all = TRUE)
ノード1) のPetal.Lengthによる分類で，setosaと{versicolor, virginica}が100%分類できていることがわかります．
後者は，ノード3) のPetal.Widthによる分類でさらに分けられていますが，versicolorに分類されてしまったvirginicaが5/50個体，virginicaに分類されてしまったversicolorが1/50個体いることがわかります．
ノード 6)では，Petal.Lengthによる分類で，versicolor 49個体とvirginica 5個体を分けています．
irisデータの場合，説明変数はカテゴリーではなく量的変数です．そうした場合，さまざまな可能性が探索されながら，もっとも良い値が各変数の分割点に採用されるようになっています．
説明変数がカテゴリーとして表されている場合でも決定木の方法は使えます．
また，predict()関数を使った予測や，rpartOrdinalパッケージを使った順序カテゴリカルデータに対する方法もありますが，ここでは省略します:( 詳細は『カテゴリカルデータ解析 (Rで学ぶデータサイエンス 1)』をご覧ください．
みせかけの相関
最後に，みせかけの相関の問題について．母集団全体と，いずれかの変数に着目して分割した集団とで，相関関係が異なっている場合があるので気をつけましょう，というものです．
例によって，ここでは詳しく述べませんが，以下などを参考にしてみてください．
シンプソンのパラドックス (Wikipedia)
シンプソンのパラドックスの図解 (社会学者の研究メモ)
ふたつの変数間に相関が見られたとしても，両方に関連する別な因子がないかを常に意識することが大切ですね．
参考
Kashiwa.R#3で発表する内容を文章にしたものです．内容のほぼすべては，『カテゴリカルデータ解析 (Rで学ぶデータサイエンス 1)』 藤井良宜 2010年 共立出版 を大いに参考にしています．
カテゴリカルデータの解析 (その1)
カテゴリカルデータの解析
「アンケート」などに代表されるカテゴリカルデータの解析法についてざっとまとめました．Kashiwa.R#3で発表する内容を文章にしたものです．
内容のほぼすべては，『カテゴリカルデータ解析 (Rで学ぶデータサイエンス 1)』 藤井良宜 2010年 共立出版 を大いに参考にしています．
要約
- xtab()で集計
- 変数間の関係性を調べるいろいろな方法がある
- みかけ上の相関に注意
カテゴリカルデータとは？
連続変数ではなくて，グループで分類されるような変数をもつデータです．「アンケート」が良い例．グループには，順序グループ (例: なし，軽度，重度) と名義グループ (例: 男性，女性，その他) の2種類があります．
ここでは，データの項目のひとつを「変数」と呼び (例: 症状，性別)，項目に含まれる分類を「カテゴリー」と呼びます (例: {なし，軽度，重度}，{男性，女性，その他})．
データをR上に載せるまで
省略します…．read.table()関数を使ったり，直にデータフレームを打ち込んだりしてみてください．
個票データと集計データ
個票データ
個票データはひとつひとつの個体情報のリストです．たとえば，vcdパッケージに含まれるArthritis．これは，関節炎に関する臨床試験のデータで，84人の患者さんひとりひとりの情報が行ごとに記録されているようです．
> install.packages("vcd")
> library(vcd)
> head(Arthritis)
ID Treatment Sex Age Improved
1 57 Treated Male 27 Some
2 46 Treated Male 29 None
3 77 Treated Male 30 None
4 17 Treated Male 32 Marked
5 36 Treated Male 46 Marked
6 23 Treated Male 58 Marked
個票データの場合，summary()でそれぞれの変数の概要を知ることができます．
> summary(Arthritis)
ID Treatment Sex Age Improved
Min. : 1.00 Placebo:43 Female:59 Min. :23.00 None :42
1st Qu.:21.75 Treated:41 Male :25 1st Qu.:46.00 Some :14
Median :42.50 Median :57.00 Marked:28
Mean :42.50 Mean :53.36
3rd Qu.:63.25 3rd Qu.:63.00
Max. :84.00 Max. :74.00
集計データ
集計データは，ひとつひとつの個体情報の集計結果です．同じくvcdパッケージに含まれる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
集計データの場合は，単にsummary()を適用しても意味がありません．
どちらにしてもxtab()
それではどうするかというと，個票データでも集計データでも，まずはxtab()を適用します．xtab()を使うことで，ひとつ以上の変数について，各カテゴリーの組み合わせごとの度数を表す表形式のデータ (クロス表) が得られます．> # Arthritisで，処置 (Treatment) と病状がどれだけ改善したか (Improved) の
> # カテゴリーの組み合わせごとの頻度．
> (arthritis.imp.tre <- xtabs(~ Improved + Treatment, data = Arthritis))
Treatment
Improved Placebo Treated
None 29 13
Some 7 7
Marked 7 21
> # DanishWalfareで，アルコール消費量 (Alcohol) と収入 (Income) の
> # カテゴリーの組み合わせごとの頻度．
> (danish.alc.inc <- xtabs(Freq ~ Alcohol + Income, data = DanishWelfare))
Income
Alcohol 0-50 50-100 100-150 >150
<1 382 748 273 936
1-2 150 567 437 929
>2 34 161 144 383個票データではデータのカウント数を指定しないのに対し，集計データでは頻度 (Freq) を指定しているのに注意してください．
なんとも簡単に集計作業ができてしまいます:)
度数でなく割合を見たいときにはprop.table()関数を使います．
marginで，何に対する割合にするかを指定します．
> round(prop.table(arthritis.imp.tre, margin = 2), 2)
Treatment
Improved Placebo Treated
None 0.67 0.32
Some 0.16 0.17
Marked 0.16 0.51
> # アルコール摂取量 (Alcohol) ごとのそれぞれの収入状況 (Income) の割合
> round(prop.table(danish.alc.inc, margin = 1), 2)
Income
Alcohol 0-50 50-100 100-150 >150
<1 0.16 0.32 0.12 0.40
1-2 0.07 0.27 0.21 0.45
>2 0.05 0.22 0.20 0.53
xtabs()を用いた操作をマスターすれば，カテゴリカルデータを解析する作業が非常に楽になります．
table()関数でも同じことができるようですが，ここでは省略してしまいます:(
その2 につづく
