===== 統計工学 2週目 =====
==== はじめに ====
=== 連絡 2012.12.11 ===
* 自分で頑張って、とお願いした、回帰分析と決定木分析のコードを追記しました。
* このページを実験時間中に改訂しましたが、[[http://stat.inf.uec.ac.jp/dokuwiki/doku.php?id=mselab:2012:stat:week2:r2&#時間内課題１_回帰分析の自習|自習部分その1]]と[[http://stat.inf.uec.ac.jp/dokuwiki/doku.php?id=mselab:2012:stat:week2:r2&#時間内課題２_決定木分析の自習|自習部分その2]]の内容はそのままです。決定木に関して、少しグラフの出力などを、追加しました。
* TICデータの[[http://stat.inf.uec.ac.jp/dokuwiki/doku.php?id=mselab:2012:stat:week2:r2&#回帰分析|重回帰分析のコード]]と、[[http://stat.inf.uec.ac.jp/dokuwiki/doku.php?id=mselab:2012:stat:week2:r2&#決定木分析|決定木分析のコード]]を追記しました。[[http://stat.inf.uec.ac.jp/dokuwiki/doku.php?id=mselab:2012:stat:week2:r2&#分類コード_水準_因子変数_の扱い|分類コードの扱い]]はほぼ今朝のまま。[[http://stat.inf.uec.ac.jp/dokuwiki/doku.php?id=mselab:2012:stat:week2:r2&#回帰分析|重回帰分析のコード]]の中に変数選択がありますが、これは参考です。
* [[http://stat.inf.uec.ac.jp/dokuwiki/doku.php?id=mselab:2012:stat:week2:r2&#練習課題_参考|練習課題(参考)]]はあくまでも参考まで、です。[[http://stat.inf.uec.ac.jp/dokuwiki/doku.php?id=mselab:2012:stat:week2:r2&#参考1|参考]]も次週の内容を含んでいるので、参考まで、です。
* [[http://stat.inf.uec.ac.jp/dokuwiki/doku.php?id=mselab:2012:stat:week2:r2&#データの説明|データの説明]]に、変数名を参考までに追記しました。
=== 概要 ===
今週の実験の内容は
* データマイニング：保険のデータのV86の予測モデルの構築
* 回帰分析
* 決定木
* ロジスティック回帰 (オプション)
である。3週間の流れがデータ分析(第1週)を行って、回帰分析等によるモデル構築(第2週)、そして第3週の他の手法も組み合わせたモデル構築と結論へと繋がるため、今週は[[http://kyoumu.office.uec.ac.jp/syllabus/2012/21/21_17123216.html|多変量解析]]との内容の重複を避けずにおいた。また同科目の履修を前提としていないため、回帰分析の学習も自習内容として含めてある。
* 先週に続いてデータの把握、特にV86を中心に。
* 回帰分析の2つの課題、同じ課題を決定木も(自習、練習に相当)
* 解析データを用いた1つの課題(本番)
に取り組んで貰う。
==== 実験の流れ ====
時間内課題：回帰分析(のたぶん復習)と決定木分析
* 回帰分析
* 配付資料とこのページにあるコードを参考に、単回帰(説明変数1つ)と重回帰(説明変数2つ)の場合を自習して、回帰分析の結果の読み方について学ぶ。(時間内は少なくとも、分析結果やグラフをすべて得て、Wordに貼るなどしてお持ち帰りするところまで)
* 重回帰(説明変数2つ)の場合で、さらに変数を分類尺度に変換し、分類尺度の説明変数を含む場合の回帰分析の結果の読み方(特に推定値Estimateの意味や解釈)について考える。
* 決定木分析
* 回帰分析と同じデータを、決定木分析にかけてみて、どのようなモデルかを考える。
課題：保険会社の顧客データのデータマイニング (時間内はできるところまで)
* まずは分類数の多い変数(V1とV5)の数を減らす。
* 回帰分析と決定木分析を行う。
* 得たモデルを考察する。
本実験に関係するメモ
- 配付資料とRの実行結果を見比べながら、lm関数を用いた回帰分析で出力される情報のどれが、配付資料のどれに対応するのかを把握する
* 回帰係数の推定値: Estimate
* 切片: Intercept
* 寄与率: R-Squared
* 自由度調整済み寄与率: Adjusted R-squared
* F値: F-statistic
* 自由度: DF, degrees of freedom
* t値: t value
* P値: Pr(>|t|), p-value
* 標準誤差: Std. Error, standard error
* 残差: Residuals
* (回帰)係数: Coefficients
* てこ比: Leverage
* 標準化残差: Studentized Residuals
- rpart関数を用いた決定木分析の出力を見ながら、決定木分析はどのようなモデルを構築するのかを検討する
- 保険データの回帰分析、に取り組む (保険データの回帰分析)
- 回帰係数の推定
- 変数の加工(今回初)
- 分散分析によるモデルの有意性の検討や回帰係数の有意性の検討
- てこ比や標準化残差などの検討
- 変数の増減 (stepAIC関数を使う?)
- 以上を繰り返す
- 保険のデータの決定木分析、に取り組む (時間外課題)
- オプション課題
- 回帰分析で変数選択を行う
- 手動でP値を見ながらの変数減少法
- AICを用いた変数増減法 stepAIC()
- 決定木の学習パラメータの調整
==== 時間内課題１：回帰分析の自習 ====
=== 単回帰分析 ===
配付資料の表4.1のデータの入力は次の通り。
x <- c(2.2,4.1,5.5,1.9,3.4,2.6,4.2,3.7,4.9,3.2)
y <- c(71,81,86,72,77,73,80,81,85,74)
table.4.1 <- data.frame(x=x,y=y)
回帰分析の実行にはlm関数を用いる。
lm(y~x, data=table.4.1)
「y~x」は、変数yを変数xで説明する回帰分析を行うための表現。これは
Y = \beta_0+\beta_1 X+\epsilon
という回帰式を推定せよ、という意味である。
回帰分析の詳細を示すには、lm関数の実行結果を別の変数lm.4.1などに代入してから、詳細を出力させるsummary関数を用いる。
lm.4.1 <- lm(y~x, data=table.4.1)
print(lm.4.1)
summary(lm.4.1)
**問１：このsummary関数の出力結果を、配付資料と対比させて理解せよ。**
回帰分析の結果を診断するために、幾つかのグラフを出力する機能がある。lm関数の実行結果をplot関数にかければ、そのようなグラフが4枚、順次出力される。
plot(test.lm)
**問２：このplot関数の出力結果を、配付資料と対比させて理解せよ。**
=== 重回帰分析 ===
表5.1のデータは次のように入力する。
x.1 <- c(51,38,57,51,53,77,63,69,72,73)
x.2 <- c(16,4,16,11,4,22,5,5,2,1)
y <- c(3.0,3.2,3.3,3.9,4.4,4.5,4.5,5.4,5.4,6.0)
table.5.1 <- data.frame(x.1=x.1,x.2=x.2,y=y)
回帰分析の実行には、説明変数が複数あっても(重回帰でも)lm関数を用いる。
lm(y~x.1+x.2, data=table.5.1)
今度は説明変数が2つあるので、「y~x.1+x.2」となる。これは
Y = \beta_0+\beta_1 X_1 + \beta_2 X_2+\epsilon
という回帰式を推定せよ、という意味である。
回帰分析の詳細を示すには、lm関数の実行結果を別の変数lm.5.1などに代入してから、詳細を出力させるsummary関数を用いる。
lm.5.1 <- lm(y~x.1+x.2, data=table.5.1)
print(lm.5.1)
summary(lm.5.1)
**問４：このplot関数の出力結果を、配付資料と対比させて理解せよ。**
=== 重回帰応用(水準変数が説明変数に含まれる場合) ===
上のデータを、広さを「広め(w)」「狭め(n)」とし、築年数も「新しめ(new)」「古め(old)」にする
x.1 <- c("n","n","n","n","n","w","w","w","w","w")
x.2 <- c("old","new","old","old","new","old","new","new","new","new")
y <- c(3.0,3.2,3.3,3.9,4.4,4.5,4.5,5.4,5.4,6.0)
table.5.1.c <- data.frame(x.1=x.1,x.2=x.2,y=y)
**問５：このデータで、回帰分析を行い、グラフも作成して、先のモデルと何が異なるか、検討せよ。**
lm.5.1.c <- lm(y~x.1+x.2, data=table.5.1.c)
print(lm.5.1.c)
summary(lm.5.1.c)
plot(lm.5.1.c)
==== 時間内課題２：決定木分析の自習 ====
まず、次の一行を実行しておく。
library(mvpart)
上の実行例で「lm」とあるところをすべて「rpart」で置き換える。
rpart.4.1 <- rpart(y~x, data=table.4.1)
print(rpart.4.1)
summary(rpart.4.1)
plot(rpart.4.1)
text(rpart.4.1)
rpart.5.1 <- rpart(y~x.1+x.2, data=table.5.1)
print(rpart.5.1)
summary(rpart.5.1)
plot(rpart.5.1)
text(rpart.5.1)
rpart.5.1.c <- rpart(y~x.1+x.2, data=table.5.1.c)
print(rpart.5.1.c)
summary(rpart.5.1.c)
plot(rpart.5.1.c)
text(rpart.5.1.c)
**問６：これらがどのようなモデルか、データと出力に照らして検討せよ。** グラフと画面出力の比較をまず行うとよい。またlm.5.1、lm.5.1.c、rpart.5.1、rpart.5.1.cの間の違いは考察せよ。
==== 課題：保険会社の顧客データのデータマイニング ====
* TICデータは、V65からV86までが22種類の保険商品の契約件数である。今回はV86の契約に漕ぎ着けるためのモデル構築、を目的とする。
* 重回帰分析を行い、分析結果を考察せよ。
* 重回帰分析の結果を用いて、指定された保険商品の契約に関する予測式を構築せよ。
* 線形モデルと重回帰分析を用いて、このデータをデータマイニングすることについて、考察せよ。
* 余裕があったら、[[http://stat.inf.uec.ac.jp/dokuwiki/doku.php?id=mselab:2012:stat:week3|リンク先]]を参考に、ロジスティック回帰にも取り組んでみよ。(重回帰と決定木がとりあえずの目標)
前の班との違い
* モデルの学習(推定、構築)は、学習用データのみでいい
* 分類変数のグループ化はできるだけ、行った方がいい
* モデルを得るだけでいい
やらなくて良いこと。
* 検証用データでの検証
* 訪問対象の明確化 (モデルのみでいいから)
=== データの説明 ===
== TIC2000 ==
[[http://kdd.ics.uci.edu/databases/tic/tic.data.html|tic.data.txt]]からの要約。
* CoIL 2000 Challengeで用いられた保険会社の顧客に関するデータ。86個の変数は、契約状況(V44-V85)と社会人口統計学的な変数(V1-V43)を含んでいる。この調査は "Can you predict who would be interested in buying a caravan insurance policy and give an explanation why?" という問いに答えるように集められた。
* このデータはオランダのデータマイニング会社Sentinent Machine Researchから提供され、現実のビジネスの問題に基づいている。学習用データ(ticdata2000.txt)は5000レコードでcaravan insurance policyの契約の有無(V86)を含んでおり、検証用データ(ticeval2000.txt)は4000レコードで契約の有無(V86)は含んでいない。検証用データの正解は、CoIL 2000 Challengeの開催時には公開されていなかったが、現在はテストデータ(tictest2000.txt)として公開されている。
* V1-V43のうち、コード化が指定されていない変数はすべて、郵便番号の一桁目のエリアを指している。たとえばV30が9ならばその顧客は郵便番号が9で始まるエリアに家を借りていることを、V31が5ならば郵便番号が5のエリアに持ち家があることを意味する。職業、社会層などもすべて、該当するエリアの箇所が郵便番号の一桁目で埋まっている。
== 変数 ==
[[http://kdd.ics.uci.edu/databases/tic/dictionary.txt|dictionary.txt]]からの抜粋と要約、の日本語版。
^変数^分類^メモ^
|V1|顧客分類２|L0でコード化されている、数字の大きさに意味なし|
|V2|住居数|大きいほど住む箇所が多い|
|V3|世帯構成員数の平均|人数|
|V4|世帯構成員の平均年齢|L1でコード化されている、年齢|
|V5|顧客分類１|L2でコード化されている、数字の大きさに意味なし|
|V6-V9|宗教|L3でコード化されている、V6+V7+V8+V9は9から12の間。それぞれの宗教を信じる割合？|
|V10-V13|結婚|場所を表す変数, 例えばV10が0ならば無し？|
|V14-V15|世帯の大きさ|L3でコード化されている、なぜかV14+V15は10以下。割合？|
|V16-V18|教育水準|L3でコード化されている、なぜかV16+V17+V18はほぼ10、それぞれの年数？割合？|
|V19-V24|職業|L3でコード化されている、なぜかV19+V20+V21+V22+V23+V24は9から13の間|
|V25-V29|社会層|L3でコード化されている、なぜかV25+V26+V27+V28+V29は9から12の間|
|V30-V31|住居|L3でコード化されている、なぜかV30+V31は9か10|
|V32-V34|自動車|L3でコード化されている、なぜかV32+V33+V34は9から11の間|
|V35-V36|健康保険|L3でコード化されている、なぜかV35+V36は9か10|
|V37-V41|収入|L3でコード化されている、なぜかV37+V38+V39+V40+V41は9から13の間|
|V42|平均収入|L3でコード化されている|
|V43|購買力|L3でコード化されている、1から8の間。|
|V44-V64|各種保険支払い額|L4でコード化|
|V65-V85|各種保険契約件数|件数|
メモの確認用のコード。
table((tic.learn$V16+tic.learn$V17+tic.learn$V18))
table((tic.learn$V19+tic.learn$V20+tic.learn$V21+tic.learn$V22+tic.learn$V23+tic.learn$V24))
table((tic.learn$V25+tic.learn$V26+tic.learn$V27+tic.learn$V28+tic.learn$V29))
table(tic.learn$V30+tic.learn$V31)
table(tic.learn$V32+tic.learn$V33+tic.learn$V34)
table(tic.learn$V35+tic.learn$V36)
table(tic.learn$V37+tic.learn$V38+tic.learn$V39+tic.learn$V40+tic.learn$V41)
== 変数詳細 ==
^Nr^Name^Description Domain^
|V1|MOSTYPE|Customer Subtype see L0|
|V2|MAANTHUI|Number of houses 1 - 10|
|V3|MGEMOMV|Avg size household 1 - 6|
|V4|MGEMLEEF|Avg age see L1|
|V5|MOSHOOFD|Customer main type see L2|
|V6|MGODRK|Roman catholic see L3|
|V7|MGODPR|Protestant ...|
|V8|MGODOV|Other religion|
|V9|MGODGE|No religion (無宗教)|
|V10|MRELGE|Married (既婚)|
|V11|MRELSA|Living together (同居)|
|V12|MRELOV|Other relation (その他)|
|V13|MFALLEEN|Singles (独身)|
|V14|MFGEKIND|Household without children (子供のいない世帯)|
|V15|MFWEKIND|Household with children (子供のいる世帯)|
|V16|MOPLHOOG|High level education (高等教育)|
|V17|MOPLMIDD|Medium level education (中等教育)|
|V18|MOPLLAAG|Lower level education (初等教育)|
|V19|MBERHOOG|High status|
|V20|MBERZELF|Entrepreneur|
|V21|MBERBOER|Farmer (農業)|
|V22|MBERMIDD|Middle management (中間管理職)|
|V23|MBERARBG|Skilled labourers (熟練労働者)|
|V24|MBERARBO|Unskilled labourers (非熟練労働者)|
|V25|MSKA|Social class A|
|V26|MSKB1|Social class B1|
|V27|MSKB2|Social class B2|
|V28|MSKC|Social class C|
|V29|MSKD|Social class D|
|V30|MHHUUR|Rented house|
|V31|MHKOOP|Home owners|
|V32|MAUT1|1 car (保有車1台)|
|V33|MAUT2|2 cars (保有車2台)|
|V34|MAUT0|No car (保有車なし)|
|V35|MZFONDS|National Health Service|
|V36|MZPART|Private health insurance|
|V37|MINKM30|Income < 30.000|
|V38|MINK3045|Income (収入) 30-45.000|
|V39|MINK4575|Income (収入) 45-75.000|
|V40|MINK7512|Income (収入) 75-122.000|
|V41|MINK123M|Income (収入) >123.000|
|V42|MINKGEM|Average income (平均収入)|
|V43|MKOOPKLA|Purchasing power class|
|V44|PWAPART|Contribution (契約高) private third party insurance see L4|
|V45|PWABEDR|Contribution (契約高) third party insurance (firms) ...|
|V46|PWALAND|Contribution (契約高) third party insurane (agriculture)|
|V47|PPERSAUT|Contribution (契約高) car policies|
|V48|PBESAUT|Contribution (契約高) delivery van policies|
|V49|PMOTSCO|Contribution (契約高) motorcycle/scooter policies|
|V50|PVRAAUT|Contribution (契約高) lorry policies|
|V51|PAANHANG|Contribution (契約高) trailer policies|
|V52|PTRACTOR|Contribution (契約高) tractor policies|
|V53|PWERKT|Contribution (契約高) agricultural machines policies |
|V54|PBROM|Contribution (契約高) moped policies|
|V55|PLEVEN|Contribution (契約高) life insurances|
|V56|PPERSONG|Contribution (契約高) private accident insurance policies|
|V57|PGEZONG|Contribution (契約高) family accidents insurance policies|
|V58|PWAOREG|Contribution (契約高) disability insurance policies|
|V59|PBRAND|Contribution (契約高) fire policies|
|V60|PZEILPL|Contribution (契約高) surfboard policies|
|V61|PPLEZIER|Contribution (契約高) boat policies|
|V62|PFIETS|Contribution (契約高) bicycle policies|
|V63|PINBOED|Contribution (契約高) property insurance policies|
|V64|PBYSTAND|Contribution (契約高) social security insurance policies|
|V65|AWAPART|Number of (契約口数) private third party insurance 1 - 12|
|V66|AWABEDR|Number of (契約口数) third party insurance (firms) ...|
|V67|AWALAND|Number of (契約口数) third party insurane (agriculture)|
|V68|APERSAUT|Number of (契約口数) car policies|
|V69|ABESAUT|Number of (契約口数) delivery van policies|
|V70|AMOTSCO|Number of (契約口数) motorcycle/scooter policies|
|V71|AVRAAUT|Number of (契約口数) lorry policies|
|V72|AAANHANG|Number of (契約口数) trailer policies|
|V73|ATRACTOR|Number of (契約口数) tractor policies|
|V74|AWERKT|Number of (契約口数) agricultural machines policies|
|V75|ABROM|Number of (契約口数) moped policies|
|V76|ALEVEN|Number of (契約口数) life insurances|
|V77|APERSONG|Number of (契約口数) private accident insurance policies|
|V78|AGEZONG|Number of (契約口数) family accidents insurance policies|
|V79|AWAOREG|Number of (契約口数) disability insurance policies|
|V80|ABRAND|Number of (契約口数) fire policies|
|V81|AZEILPL|Number of (契約口数) surfboard policies|
|V82|APLEZIER|Number of (契約口数) boat policies|
|V83|AFIETS|Number of (契約口数) bicycle policies|
|V84|AINBOED|Number of (契約口数) property insurance policies|
|V85|ABYSTAND|Number of (契約口数) social security insurance policies|
|V86|CARAVAN|Number of (契約口数) mobile home policies 0 - 1|
== 各変数のコーディング ==
L0:分類を表す数字なので、大小関係に意味がなく、名義尺度である。そのままでは説明変数にならない。
^Value^Label^
|1|High Income, expensive child|
|2|Very Important Provincials|
|3|High status seniors|
|4|Affluent senior apartments|
|5|Mixed seniors|
|6|Career and childcare|
|7|Dinki's (double income no kids)|
|8|Middle class families|
|9|Modern, complete families|
|10|Stable family|
|11|Family starters|
|12|Affluent young families|
|13|Young all american family|
|14|Junior cosmopolitan|
|15|Senior cosmopolitans|
|16|Students in apartments|
|17|Fresh masters in the city|
|18|Single youth|
|19|Suburban youth|
|20|Etnically diverse|
|21|Young urban have-nots|
|22|Mixed apartment dwellers|
|23|Young and rising|
|24|Young, low educated |
|25|Young seniors in the city|
|26|Own home elderly|
|27|Seniors in apartments|
|28|Residential elderly|
|29|Porchless seniors: no front yard|
|30|Religious elderly singles|
|31|Low income catholics|
|32|Mixed seniors|
|33|Lower class large families|
|34|Large family, employed child|
|35|Village families|
|36|Couples with teens 'Married with children'|
|37|Mixed small town dwellers|
|38|Traditional families|
|39|Large religous families|
|40|Large family farms|
|41|Mixed rurals|
L1:大きさが年齢の順なので、そのまま説明変数に使える。
|1|20-30 years|
|2|30-40 years|
|3|40-50 years|
|4|50-60 years|
|5|60-70 years|
|6|70-80 years|
L2:数字は分類を表すだけなので、連続尺度でも順序尺度でもなく、名義尺度。そのままでは説明変数にならない。
|1|Successful hedonists|
|2|Driven Growers|
|3|Average Family|
|4|Career Loners|
|5|Living well|
|6|Cruising Seniors|
|7|Retired and Religeous|
|8|Family with grown ups|
|9|Conservative families|
|10|Farmers|
L3:順序尺度。このまま連続尺度の説明変数として用いる。
|0|0%|
|1|1 - 10%|
|2|11 - 23%|
|3|24 - 36%|
|4|37 - 49%|
|5|50 - 62%|
|6|63 - 75%|
|7|76 - 88%|
|8|89 - 99%|
|9|100%|
L4: 順序尺度。今回はこのまま連続尺度の変数として用いる。
|0|f 0|
|1|f 1 - 49|
|2|f 50 - 99|
|3|f 100 - 199|
|4|f 200 - 499|
|5|f 500 - 999|
|6|f 1000 - 4999|
|7|f 5000 - 9999|
|8|f 10.000 - 19.999|
|9|f 20.000 - ?|
== データのダウンロードと読み込み ==
演習で用いる保険データは、大学内からであれば、Rに次の命令を実行すれば読み込める。
Sys.setenv("http_proxy"="http://130.153.8.66:8080/")
tic.learn <- read.table("http://kdd.ics.uci.edu/databases/tic/ticdata2000.txt")
tic.eval <- read.table("http://kdd.ics.uci.edu/databases/tic/ticeval2000.txt")
tic.test <- read.table("http://kdd.ics.uci.edu/databases/tic/tictgts2000.txt")
tic.eval <- cbind(tic.eval, tic.test)
colnames(tic.eval)[86] <- "V86"
rm(tic.test)
学外もしくは自宅等であれば、最初の1行(Sys.setenvで始まる行)は不要。
== 参考 ==
kernlabパッケージに、加工済みのデータが入っていて、それを使うこともできる。
install.packages(c("kernlab"), dependencies=TRUE)
tic.learn <- ticdata[1:5822,]
tic.eval <- ticdata[5823:9822,]
=== 準備 ===
上のデータの読み込みを実行済みであれば、この課題ではMASSライブラリを使う可能性があるので、それを読み込んでおく。
library(MASS)
=== 回帰分析 ===
== 回帰分析をするだけ ==
V1からV85を説明変数にして、V86を予測するだけなら、次の命令をコピーして貼り付けたら、パラメータの最小二乗推定値は得られる。
lm(V86~V1 +V2 +V3 +V4 +V5 +V6 +V7 +V8 +V9 +V10
+V11+V12+V13+V14+V15+V16+V17+V18+V19+V20
+V21+V22+V23+V24+V25+V26+V27+V28+V29+V30
+V31+V32+V33+V34+V35+V36+V37+V38+V39+V40
+V41+V42+V43+V44+V45+V46+V47+V48+V49+V50
+V51+V52+V53+V54+V55+V56+V57+V58+V59+V60
+V61+V62+V63+V64+V65+V66+V67+V68+V69+V70
+V71+V72+V73+V74+V75+V76+V77+V78+V79+V80
+V81+V82+V83+V84+V85, data=tic.learn)
先ほどと同様に、と回帰分析を実行するには、次のようにすれば、各種統計量も診断用のグラフも出力される。
lm.learn.0 <- lm(V86~V1 +V2 +V3 +V4 +V5 +V6 +V7 +V8 +V9 +V10
+V11+V12+V13+V14+V15+V16+V17+V18+V19+V20
+V21+V22+V23+V24+V25+V26+V27+V28+V29+V30
+V31+V32+V33+V34+V35+V36+V37+V38+V39+V40
+V41+V42+V43+V44+V45+V46+V47+V48+V49+V50
+V51+V52+V53+V54+V55+V56+V57+V58+V59+V60
+V61+V62+V63+V64+V65+V66+V67+V68+V69+V70
+V71+V72+V73+V74+V75+V76+V77+V78+V79+V80
+V81+V82+V83+V84+V85, data=tic.learn)
print(lm.learn.0)
summary(lm.learn.0)
plot(lm.learn.0)
ただ、これではV1やV5の扱いを「数値」としてしまっているし、有意でない(*や**や***がついてない)変数も多くモデルに取り込まれてしまっている。
そこで、データマイニング手法を「線形モデル」に限るならば、モデルに対して次の2つの対応をとることが必要になる。
- 予測の役に立たない不要な変数を削る
- 本来は「分類コード」として扱うべき変数を、「数値データ」として扱っている現状を解消する
それぞれ順に説明する。
== 予測に不要な変数を削る ==
変数増減法という方法がある。上で得たlm.learn.0という「すべての変数を説明変数に加えた回帰モデル」に対して、
lm.learn.AIC <- stepAIC(lm.learn.0, direction="both")
との1行を実行すると、AICという基準に照らして、不要な変数を1つずつ取り除か加えるか、検討してくれて、一番不要な変数から順に出し入れすることで、最適なモデルに到達するよう頑張ってくれる。
=== 分類コード(水準, 因子変数)の扱い ===
V1とV5の中身は、前述のように分類コードである。前回の課題で、ヒストグラムを描こうとすると、エラーが発生したかもしれない。
Rでは、中身が分類コードや水準である変数のことを、因子変数(factor)と呼ぶ。
V1とV5を因子変数に変換するには、次の2行を実行して、V1とV5を「分類尺度」に変換した変数をそれぞれ、V1fおよびV5fと名付ける。
tic.learn$V1f <- as.factor(tic.learn$V1)
tic.learn$V5f <- as.factor(tic.learn$V5)
通常の数値変数(numeric)と因子変数(factor)の違いは、以下、幾つかのグラフを描いているので、1行ずつ貼って実行していくと、違いが分かる。
V1やV5のヒストグラムは、hist関数を用いて
hist(tic.learn$V1)
hist(tic.learn$V5)
などで描ける。一方、V1fやV5fのヒストグラムは
hist(tic.learn$V1f)
hist(tic.learn$V5f)
を実行するとエラーになる。これを得るには、table関数で集計してから、barplot関数で棒グラフを描く。
barplot(table(tic.learn$V1f))
barplot(table(tic.learn$V5f)
更にカテゴリごとの契約の有無の割合を見るには
barplot(table(tic.learn$V86,tic.learn$V1f))
barplot(table(tic.learn$V86,tic.learn$V5f))
とする。色の濃い部分が契約無、色の薄い部分が契約有である。
参考までに、因子変数を用いた回帰分析を実行するには、
lm.learn.1 <- lm(V86~V1f+V2 +V3 +V4 +V5f+V6 +V7 +V8 +V9 +V10
+V11+V12+V13+V14+V15+V16+V17+V18+V19+V20
+V21+V22+V23+V24+V25+V26+V27+V28+V29+V30
+V31+V32+V33+V34+V35+V36+V37+V38+V39+V40
+V41+V42+V43+V44+V45+V46+V47+V48+V49+V50
+V51+V52+V53+V54+V55+V56+V57+V58+V59+V60
+V61+V62+V63+V64+V65+V66+V67+V68+V69+V70
+V71+V72+V73+V74+V75+V76+V77+V78+V79+V80
+V81+V82+V83+V84+V85, data=tic.learn)
print(lm.learn.1)
summary(lm.learn.1)
plot(lm.learn.1)
とすればよい。ただ、これではV5fについての推定ができず、まだまだ問題があることがわかる。
== 分類変数のグループ化 ==
推定に問題が発生したり、カテゴリが多すぎることを解消するには、V1とV86のtable関数による集計結果や、V5とV86の集計結果を参考に、V1とV5について「幾つかのカテゴリをまとめてしまうグルーピング」という操作を行うことが1案である。
回帰分析でも決定木分析でも数千のサンプルで40以上の分類は、これだけ変数があると細かすぎる。
まず集計方法を再掲する。
library(mvpart)
table(tic.learn$V1f, tic.learn$V86)
table(tic.learn$V5f, tic.learn$V86)
これらをグラフに描くには
barchart(table(tic.learn$V1f, tic.learn$V86))
barchart(table(tic.learn$V5f, tic.learn$V86))
とする。
さて、V86の予測にV1のどのカテゴリが不要かを見るために、決定木分析を行ってみる。
決定木(実際には回帰木)を用いて、V86の増加に寄与しそうなカテゴリのみを残す。
目論見は次の2点。
* 契約にあまり寄与しないカテゴリには、コード0を割り当ててしまおう。
* 契約にしそうなカテゴリは残そう。
library(mvpart)
rpart(V86~V1f, data=tic.learn, control=rpart.control(cp=0))
剪定せずにそのまま出力した回帰樹は次の通り。
n= 5822
node), split, n, deviance, yval
* denotes terminal node
1) root 5822 327.1989000 0.05977327
2) V1f=2,4,5,7,9,10,11,13,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41 4880 220.0654000 0.04733607
4) V1f=4,5,9,15,16,17,18,19,21,22,23,24,25,26,27,28,29,30,31,35,40,41 2075 56.3788000 0.02795181
8) V1f=15,16,17,18,19,21,23,25,26,27,28,29,40,41 885 14.7457600 0.01694915
16) V1f=15,16,17,18,19,21,28,40 163 0.0000000 0.00000000 *
17) V1f=23,25,26,27,29,41 722 14.6883700 0.02077562
34) V1f=23 251 3.9362550 0.01593625 *
35) V1f=25,26,27,29,41 471 10.7431000 0.02335456
70) V1f=26,27 98 1.9591840 0.02040816
140) V1f=27 50 0.9800000 0.02000000 *
141) V1f=26 48 0.9791667 0.02083333 *
71) V1f=25,29,41 373 8.7828420 0.02412869
142) V1f=29 86 1.9534880 0.02325581 *
143) V1f=25,41 287 6.8292680 0.02439024 *
9) V1f=4,5,9,22,24,30,31,35 1190 41.4462200 0.03613445
18) V1f=24,30,31 503 14.5526800 0.02982107
36) V1f=24,31 385 10.6857100 0.02857143
72) V1f=24 180 4.8611110 0.02777778 *
73) V1f=31 205 5.8243900 0.02926829 *
37) V1f=30 118 3.8644070 0.03389831 *
19) V1f=4,5,9,22,35 687 26.8588100 0.04075691
38) V1f=4,35 266 9.6240600 0.03759398
76) V1f=35 214 7.7009350 0.03738318 *
77) V1f=4 52 1.9230770 0.03846154 *
39) V1f=5,9,22 421 17.2304000 0.04275534
78) V1f=22 98 3.8367350 0.04081633 *
79) V1f=5,9 323 13.3931900 0.04334365
158) V1f=9 278 11.4820100 0.04316547 *
159) V1f=5 45 1.9111110 0.04444444 *
5) V1f=2,7,10,11,13,20,32,33,34,36,37,38,39 2805 162.3301000 0.06167558
10) V1f=10,11,32,33,34,39 1779 94.3788600 0.05621135
20) V1f=34 182 8.5549450 0.04945055 *
21) V1f=10,11,32,33,39 1597 85.8146500 0.05698184
42) V1f=10 165 8.5090910 0.05454545 *
43) V1f=11,32,33,39 1432 77.3044700 0.05726257
86) V1f=32,33 951 50.9337500 0.05678233
172) V1f=32 141 7.5460990 0.05673759 *
173) V1f=33 810 43.3876500 0.05679012 *
87) V1f=11,39 481 26.3700600 0.05821206
174) V1f=39 328 17.8993900 0.05792683 *
175) V1f=11 153 8.4705880 0.05882353 *
11) V1f=2,7,13,20,36,37,38 1026 67.8060400 0.07115010
22) V1f=7,38 383 24.2349900 0.06788512
44) V1f=38 339 21.4395300 0.06784661 *
45) V1f=7 44 2.7954550 0.06818182 *
23) V1f=2,13,20,36,37 643 43.5645400 0.07309487
46) V1f=2,13,36 486 32.4794200 0.07201646
92) V1f=36 225 14.8622200 0.07111111 *
93) V1f=2,13 261 17.6168600 0.07279693
186) V1f=13 179 12.0558700 0.07262570 *
187) V1f=2 82 5.5609760 0.07317073 *
47) V1f=20,37 157 11.0828000 0.07643312
94) V1f=37 132 9.2424240 0.07575758 *
95) V1f=20 25 1.8400000 0.08000000 *
3) V1f=1,3,6,8,12 942 102.4682000 0.12420380
6) V1f=1,3,6 492 44.9187000 0.10162600
12) V1f=3,6 368 33.2798900 0.10054350
24) V1f=3 249 22.4899600 0.10040160 *
25) V1f=6 119 10.7899200 0.10084030 *
13) V1f=1 124 11.6371000 0.10483870 *
7) V1f=8,12 450 57.0244400 0.14888890
14) V1f=12 111 13.6936900 0.14414410 *
15) V1f=8 339 43.3274300 0.15044250 *
上の決定木の結果を見るなどして、下の「V1の値をV1grに代入するだけのコード」の右辺をいろいろ変えてみて、「V1の値の種類を減らす工夫」(グルーピング)をする。
これをメモ帳に貼り付けて幾つかを0にしたり、幾つかの変数に同じ値を割り振るなど、が実際の作業となる。
勿論、グルーピングの結果はレポートに記すこと。
tic.learn$V1gr <- 0
tic.learn$V1gr[tic.learn$V1==1] <- 1
tic.learn$V1gr[tic.learn$V1==2] <- 2
tic.learn$V1gr[tic.learn$V1==3] <- 3
tic.learn$V1gr[tic.learn$V1==4] <- 4
tic.learn$V1gr[tic.learn$V1==5] <- 5
tic.learn$V1gr[tic.learn$V1==6] <- 6
tic.learn$V1gr[tic.learn$V1==7] <- 7
tic.learn$V1gr[tic.learn$V1==8] <- 8
tic.learn$V1gr[tic.learn$V1==9] <- 9
tic.learn$V1gr[tic.learn$V1==10] <- 10
tic.learn$V1gr[tic.learn$V1==11] <- 11
tic.learn$V1gr[tic.learn$V1==12] <- 12
tic.learn$V1gr[tic.learn$V1==13] <- 13
#tic.learn$V1gr[tic.learn$V1==14] <- 14
tic.learn$V1gr[tic.learn$V1==15] <- 15
tic.learn$V1gr[tic.learn$V1==16] <- 16
tic.learn$V1gr[tic.learn$V1==17] <- 17
tic.learn$V1gr[tic.learn$V1==18] <- 18
tic.learn$V1gr[tic.learn$V1==19] <- 19
tic.learn$V1gr[tic.learn$V1==20] <- 10
tic.learn$V1gr[tic.learn$V1==21] <- 21
tic.learn$V1gr[tic.learn$V1==22] <- 22
tic.learn$V1gr[tic.learn$V1==23] <- 23
tic.learn$V1gr[tic.learn$V1==24] <- 24
tic.learn$V1gr[tic.learn$V1==25] <- 25
tic.learn$V1gr[tic.learn$V1==26] <- 26
tic.learn$V1gr[tic.learn$V1==27] <- 27
tic.learn$V1gr[tic.learn$V1==28] <- 28
tic.learn$V1gr[tic.learn$V1==29] <- 29
tic.learn$V1gr[tic.learn$V1==30] <- 30
tic.learn$V1gr[tic.learn$V1==31] <- 31
tic.learn$V1gr[tic.learn$V1==32] <- 32
tic.learn$V1gr[tic.learn$V1==33] <- 33
tic.learn$V1gr[tic.learn$V1==34] <- 34
tic.learn$V1gr[tic.learn$V1==35] <- 35
tic.learn$V1gr[tic.learn$V1==36] <- 36
tic.learn$V1gr[tic.learn$V1==37] <- 37
tic.learn$V1gr[tic.learn$V1==38] <- 38
tic.learn$V1gr[tic.learn$V1==39] <- 39
tic.learn$V1gr[tic.learn$V1==40] <- 40
tic.learn$V1gr[tic.learn$V1==41] <- 41
tic.learn$V1gr <- as.factor(tic.learn$V1gr)
V5についても同様。
rpart(V86~V5f, data=tic.learn, control=rpart.control(cp=0))
V5の値も幾つかのグループにまとめる。
tic.learn$V5gr <- 0
tic.learn$V5gr[tic.learn$V5==1] <- 1
tic.learn$V5gr[tic.learn$V5==2] <- 2
tic.learn$V5gr[tic.learn$V5==3] <- 3
tic.learn$V5gr[tic.learn$V5==4] <- 4
tic.learn$V5gr[tic.learn$V5==5] <- 5
tic.learn$V5gr[tic.learn$V5==6] <- 6
tic.learn$V5gr[tic.learn$V5==7] <- 7
tic.learn$V5gr[tic.learn$V5==8] <- 8
tic.learn$V5gr[tic.learn$V5==9] <- 9
tic.learn$V5gr[tic.learn$V5==10] <- 10
tic.learn$V5gr <- as.factor(tic.learn$V5gr)
参考までに、V1とV5で回帰樹を作ってみると・・・。
rpart(V86~V1f+V5f, data=tic.learn, control=rpart.control(cp=0))
保険データへの決定木分析の適用に際しては、[[http://stat.inf.uec.ac.jp/dokuwiki/doku.php?id=dmb:2011:q2|このページ]]や[[http://stat.inf.uec.ac.jp/dokuwiki/doku.php?id=mselab:2012:stat:week3|このページ]]が参考になるかもしれない。
ここまで行うと、回帰分析を実行する関数は、次のようになる。(V1grとV5grを定義していなければエラーが表示されるだけ)
lm.learn.2 <- lm(V86~V1gr+V2 +V3 +V4 +V5gr+V6 +V7 +V8 +V9 +V10
+V11+V12+V13+V14+V15+V16+V17+V18+V19+V20
+V21+V22+V23+V24+V25+V26+V27+V28+V29+V30
+V31+V32+V33+V34+V35+V36+V37+V38+V39+V40
+V41+V42+V43+V44+V45+V46+V47+V48+V49+V50
+V51+V52+V53+V54+V55+V56+V57+V58+V59+V60
+V61+V62+V63+V64+V65+V66+V67+V68+V69+V70
+V71+V72+V73+V74+V75+V76+V77+V78+V79+V80
+V81+V82+V83+V84+V85, data=tic.learn)
print(lm.learn.2)
summary(lm.learn.2)
plot(lm.learn.2)
=== 決定木分析 ===
== とりあえず決定木分析を実行してみる ==
rpart.learn.0 <- rpart(V86~V1f+V2 +V3 +V4 +V5f+V6 +V7 +V8 +V9 +V10
+V11+V12+V13+V14+V15+V16+V17+V18+V19+V20
+V21+V22+V23+V24+V25+V26+V27+V28+V29+V30
+V31+V32+V33+V34+V35+V36+V37+V38+V39+V40
+V41+V42+V43+V44+V45+V46+V47+V48+V49+V50
+V51+V52+V53+V54+V55+V56+V57+V58+V59+V60
+V61+V62+V63+V64+V65+V66+V67+V68+V69+V70
+V71+V72+V73+V74+V75+V76+V77+V78+V79+V80
+V81+V82+V83+V84+V85, data=tic.learn)
print(rpart.learn.0)
summary(rpart.learn.0)
plot(rpart.learn.0)
text(rpart.learn.0)
== 決定木分析をいろいろチューニングする (オプション課題) ==
rpart関数の分割を決める制御変数には
^rpart.controlの引数^意味^デフォルト値^
|minsplit|ノードの分割を試みる最小のレコード数|20|
|minbucket|終端ノードのレコード数の最小値|round(minsplit/3)|
|cp|決定木の複雑さを調整するパラメータ|0.05|
|maxdepth|決定木の最大の深さ|30|
などがある．他にも，
^rpart.controlの引数^意味^デフォルト値^
|maxcompete|the number of competitor splits retained in the output|4|
|maxsurrogate|the number of surrogate splits retained in the output|5|
|usesurrrogate|how to use surrogates in the splitting process|2|
|xval|number of cross-validations|10|
|surrogatestyle|controls the selection of a best surrogate|0|
などがあるが，このデータの分析では使わない．
例えば次のように、分割する際のノードの最小データ数minsplitを5、複雑さの最小値cpを0.001とすると、結構、複雑な決定木に成長する。
rpart.learn.1 <- rpart(V86~V1f+V2 +V3 +V4 +V5f+V6 +V7 +V8 +V9 +V10
+V11+V12+V13+V14+V15+V16+V17+V18+V19+V20
+V21+V22+V23+V24+V25+V26+V27+V28+V29+V30
+V31+V32+V33+V34+V35+V36+V37+V38+V39+V40
+V41+V42+V43+V44+V45+V46+V47+V48+V49+V50
+V51+V52+V53+V54+V55+V56+V57+V58+V59+V60
+V61+V62+V63+V64+V65+V66+V67+V68+V69+V70
+V71+V72+V73+V74+V75+V76+V77+V78+V79+V80
+V81+V82+V83+V84+V85, data=tic.learn,
control=c(minsplit=5, cp=0.001))
print(rpart.learn.1)
summary(rpart.learn.1)
plot(rpart.learn.1)
text(rpart.learn.1)
これはきっと、成長させすぎ・・・。複雑さの下限cpをさらに0.0001と小さくすると，より複雑な木が得られる。
グラフもノード数が増えるので，字が小さくなる。
==== レポート提出について ====
レポート提出要領:下記「XXXXXXX」は各自の学籍番号(半角文字)で置き換えること
^項目^指定^
|提出期限|実験実施の翌週の月曜日の午後7時0分まで (今回は12月17日の午後7時)|
|提出方法|電子メールに添付 (宛先は配付資料に記載)|
|ファイル形式|Wordファイル (LaTeXで作成する場合は、dvipdfmxでPDFに変換すること)|
|メールの件名|統計工学実験2レポート提出(XXXXXXX)|
|レポートファイルの名称|統計工学実験2_XXXXXXX.doc あるいは 統計工学実験2_XXXXXXX.docx|
|提出部数|レポートは各自1通ずつ。{{:mselab:report-header-2012.doc|レポートの表紙}}に、共同実験者の学籍番号と氏名を記すこと。|
==== 参考文献 ====
* 永田・棟近 (2001) [[http://www.saiensu.co.jp/?page=book_details&ISBN=ISBN978-4-7819-0980-6&YEAR=2001|多変量解析法入門]], サイエンス社
==== 練習課題(参考) ====
決定木について、理解が難しい人は、この部分の決定木の箇所だけでも実行して考えてみるといい。この箇所は今回の課題ではない。
=== タイタニック号 ===
タイタニック号の乗客の生死のデータがある。Rで
Titanic
と実行すると、表示される。
これを個別のレコードに展開し、更に救出の優先順位を高く設定された女性もしくは子供と、低めに設定された成人男性という変数も加える。
data("Titanic", package = "datasets")
titanic <- as.data.frame(Titanic)
titanic <- titanic[rep(1:nrow(titanic), titanic$Freq), 1:4]
names(titanic)[2] <- "Gender"
titanic <- transform(titanic, Treatment = factor(
Gender == "Female" | Age == "Child",
levels = c(FALSE, TRUE),
labels = c("Normal\n(Male&Adult)", "Preferential\n(Female|Child)")
))
さらに、数値に変換したデータも用意する。これは、線形学習機械とロジスティック学習機械のために用いる。
titanic.2 <- data.frame(Gender=(titanic$Gender=="Female")*1,
Age=(titanic$Age=="Child")*1,
Survived=(titanic$Survived=="Yes")*1,
Class=(titanic$Class=="Crew")*1+(titanic$Class=="3rd")*2+(titanic$Class=="2nd")*3+(titanic$Class=="1st")*4)
|変数|数値化情報|
|性別(Gender)|女性(Female)=1, 男性(Male)=0|
|年齢(Age)|子供(Child)=1, 大人(Adult)=0|
|客室等級(Class)|1等(1st)=4,2等(2nd)=3,3等(3rd)=2,乗組員(Crew)=1|
|生存(Survived)|生存(Yes)=1,死亡(No)=0|
=== 先週の課題と同等の課題 ===
先週の課題は、タイタニック号で、生存率が高くなる条件を求めよ、という問題と同等。
まず titanic というデータに含まれる変数の一覧を
names(titanic)
で取り出す。すると
[1] "Class" "Gender" "Age" "Survived" "Treatment"
のように5個の変数があることが分かる。
これを用いて、生存(Survived)についての集計を、客室等級(Class)、性別(Gender)、年齢層(Age)の組み合わせで行うには、次の1行を実行する。
「$」の前がデータ名、「$」の後ろに変数名(フィールド名)を付ける。
ftable(titanic$Gender,titanic$Age,titanic$Class, titanic$Survived)
すると
No Yes
Male Child 1st 0 5
2nd 0 11
3rd 35 13
Crew 0 0
Adult 1st 118 57
2nd 154 14
3rd 387 75
Crew 670 192
Female Child 1st 0 1
2nd 0 13
3rd 17 14
Crew 0 0
Adult 1st 4 140
2nd 13 80
3rd 89 76
Crew 3 20
と出力されるのを、各自確認せよ。
同じことは、数値化したデータを用いても得られる。
ftable(titanic.2$Gender,titanic.2$Age,titanic.2$Class, titanic.2$Survived)
上の出力と下の出力を見比べて、各変数の値の数値化を確認せよ。
0 1
0 0 1 670 192
2 387 75
3 154 14
4 118 57
1 1 0 0
2 35 13
3 0 11
4 0 5
1 0 1 3 20
2 89 76
3 13 80
4 4 140
1 1 0 0
2 17 14
3 0 13
4 0 1
=== 線形学習機械 ===
上のデータ(titanic.2)に対して
lm(Survived~., data=titanic.2)
を実行すると、全変数を用いた線形学習機械が最小二乗法により、学習される。学習結果のみなら、この一行で
Call:
lm(formula = Survived ~ ., data = titanic.2)
Coefficients:
(Intercept) Gender Age Class
0.11725 0.46493 0.09655 0.05029
と出力される。
この結果に、統計的推測の結果を付与するなら
summary(lm(Survived~., data=titanic.2))
を実行すればよい。summary()の出力は次のように得られる。
Call:
lm(formula = Survived ~ ., data = titanic.2)
Residuals (残差の分布):
Min 1Q Median 3Q Max
-0.7833 -0.2178 -0.1675 0.2207 0.8325
Coefficients (回帰係数):
Estimate Std. Error t value Pr(>|t|) 推定値 標準誤差 t値 P値
(Intercept) 0.117252 0.019113 6.135 1.01e-09 *** 切片
Gender 0.464930 0.023325 19.933 < 2e-16 ***
Age 0.096553 0.040856 2.363 0.0182 *
Class 0.050289 0.008989 5.594 2.49e-08 ***
---
Signif. codes (有意水準の表示法): 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error (残差の標準偏差): 0.4131 on 2197 degrees of freedom
Multiple R-squared (決定係数): 0.2209, Adjusted R-squared (自由度調整済み決定係数): 0.2198
F-statistic (F統計量、分散分析の結果): 207.7 on 3 and 2197 DF, p-value: < 2.2e-16
== 線形学習機械のチューニング ==
AICを用いた説明変数の選択を行うのであれば、
library(MASS)
とMASSパッケージを読み込んでから使える関数stepAIC()を用いて、
stepAIC(lm(Survived~., data=titanic.2))
を実行する。
Start: AIC=-3887.24
Survived ~ Gender + Age + Class
Df Sum of Sq RSS AIC
374.99 -3887.2
- Age 1 0.953 375.95 -3883.7
- Class 1 5.342 380.33 -3858.1
- Gender 1 67.815 442.81 -3523.4
Call:
lm(formula = Survived ~ Gender + Age + Class, data = titanic.2)
Coefficients:
(Intercept) Gender Age Class
0.11725 0.46493 0.09655 0.05029
これはAICを変数増減法で最適化する関数であり、最終的に採用するモデルを出力してくれる便利な関数である。
== 関係するグラフ ==
titanic.lm <- lm(Survived~., data=titanic.2)
par(mfrow=c(2,2))
plot(titanic.lm)
par(mfrow=c(1,1))
== 学習結果の解釈 ==
上の回帰分析から
Survived = 0.11725 + 0.46493 * Gender + 0.09655 * Age + 0.05029 * Class
という学習結果が得られた。Survivedを大きくする(=生存する)ために、Genderが0よりは1の方が良いことは、回帰係数 0.46493 を見れば分かる。
同様に Age も0よりは1の方が良く、Classも1よりは4の方が良い。
そしてどの変数も有意であることから、生存するためには、
* 男性よりは女性
* 大人よりは子供
* 客室等級は高めの部屋を選ぶ
のが好条件となることが分かる。
=== ロジスティック線形学習機械 ===
生存確率をpとして
log(p/1-p)
を与える学習機械が、ロジスティック回帰モデルである。Rではglm()という関数を用いて
glm(Survived~., data=titanic.2, family="binomial")
で、ロジスティック回帰モデルを最尤推定で学習させることができる。
結果は
Call: glm(formula = Survived ~ ., family = "binomial", data = titanic.2)
Coefficients:
(Intercept) Gender Age Class
-1.8622 2.0580 0.5115 0.2783
Degrees of Freedom: 2200 Total (i.e. Null); 2197 Residual
Null Deviance: 2769
Residual Deviance: 2299 AIC: 2307
と、lm()と似た出力を得る。
これもsummary()を加えて
summary(glm(Survived~., data=titanic.2, family="binomial"))
を実行すると、
Call:
glm(formula = Survived ~ ., family = "binomial", data = titanic.2)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.7597 -0.6926 -0.6109 0.7055 1.8818
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.86224 0.11615 -16.033 < 2e-16 ***
Gender 2.05802 0.12604 16.328 < 2e-16 ***
Age 0.51147 0.22292 2.294 0.0218 *
Class 0.27834 0.05047 5.515 3.48e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 2769.5 on 2200 degrees of freedom
Residual deviance: 2299.2 on 2197 degrees of freedom
AIC: 2307.2
Number of Fisher Scoring iterations: 4
と、種々の検定統計量が一緒に出力される。
== ロジスティック学習機械のチューニング ==
AICを用いた説明変数の選択を行うのであれば、
library(MASS)
とMASSパッケージを読み込んでから使える関数stepAIC()を用いて、
stepAIC(glm(Survived~., data=titanic.2, family="binomial"))
を実行する。
Start: AIC=2307.21
Survived ~ Gender + Age + Class
Df Deviance AIC
2299.2 2307.2
- Age 1 2304.4 2310.4
- Class 1 2329.1 2335.1
- Gender 1 2595.5 2601.5
Call: glm(formula = Survived ~ Gender + Age + Class, family = "binomial",
data = titanic.2)
Coefficients:
(Intercept) Gender Age Class
-1.8622 2.0580 0.5115 0.2783
Degrees of Freedom: 2200 Total (i.e. Null); 2197 Residual
Null Deviance: 2769
Residual Deviance: 2299 AIC: 2307
これはAICを変数増減法で最適化する関数であり、最終的に採用するモデルを出力してくれる便利な関数である。
== 関係するグラフ ==
titanic.glm <- glm(Survived~., data=titanic.2, family="binomial")
par(mfrow=c(2,2))
plot(titanic.glm)
par(mfrow=c(1,1))
== 学習結果の解釈 ==
上の回帰分析から
log(Pr[Survived]/(1-Pr[Survived])) = -1.86224 + 2.050802 * Gender + 0.51147 * Age + 0.27834 * Class
という学習結果が得られた。Pr[Survived](生存する確率)を大きくするために、Genderが0よりは1の方が良いことは、回帰係数 2.050802 を見れば分かる。
同様に Age も0よりは1の方が良く、Classも1よりは4の方が良い。
そしてどの変数も有意であることから、生存するためには、
* 男性よりは女性
* 大人よりは子供
* 客室等級は高めの部屋を選ぶ
のが好条件となることが分かる。
=== 決定木 ===
決定木は生存確率pの高低を際立たせるような、データの分割を表現するモデルである。Rではrpartパッケージの中のrpart()、もしくはmvpartパッケージの中のmvpart()という関数を用いて
library(mvpart)
rpart(Survived~.,data=titanic)
で、決定木を学習させることができる。
結果は
n= 2201
node), split, n, loss, yval, (yprob)
* denotes terminal node
1) root 2201 711 No (0.6769650 0.3230350)
2) Gender=Male 1731 367 No (0.7879838 0.2120162)
4) Age=Adult 1667 338 No (0.7972406 0.2027594) *
5) Age=Child 64 29 No (0.5468750 0.4531250)
10) Class=3rd 48 13 No (0.7291667 0.2708333) *
11) Class=1st,2nd 16 0 Yes (0.0000000 1.0000000) *
3) Gender=Female 470 126 Yes (0.2680851 0.7319149)
6) Class=3rd 196 90 No (0.5408163 0.4591837) *
7) Class=1st,2nd,Crew 274 20 Yes (0.0729927 0.9270073) *
と、これまでとは違ったものになる。各行はノードと呼ばれる分割単位を表し、
たとえば
1) root 2201 711 No (0.6769650 0.3230350)
は、
* データ全体(root)はレコード数が2201
* このノードの代表値(一番多い値)はNo
* 2201のうち711はNoでない
* YesとNoの比は0.6769650:0.3230350
を意味する。これを性別で分割すると
2) Gender=Male 1731 367 No (0.7879838 0.2120162)
3) Gender=Female 470 126 Yes (0.2680851 0.7319149)
を得た、とある。これは
* データ全体を性別で分割すると、生存確率がもっとも差が開く
* 男性の代表値はNo(死亡)、0.7879838は死亡確率を意味するので、生存確率は0.2120162
* 女性の代表値はYes(生存)、0.7319149は生存確率
と読み取ることができる。
生存確率の高いノードを探すと
7) Class=1st,2nd,Crew 274 20 Yes (0.0729927 0.9270073) *
が目に付く。これは、
- Gender=Female
- Class=1st,2nd,Crew
と条件がふたつついたノードであり、「女性かつ１等客室の乗客、２等客室の乗客、または乗組員」となる。
この属性を持つ人々は、生存確率が0.9270073と最も高い。
逆に最も生存確率が低いのは
10) Class=3rd 48 13 No (0.7291667 0.2708333) *
であり、このノードまでを
- Gender=Male
- Age=Adult
- Class=3rd
と辿れることから、「男性で成人で３等客室の乗客」は生存確率が0.2708333と最も低かったことが分かる。
これもsummary()を加えて
summary(rpart(Survived~.,data=titanic))
を実行すると、
Call:
rpart(formula = Survived ~ ., data = titanic)
n= 2201
CP nsplit rel error xerror xstd
1 0.30661041 0 1.0000000 1.0000000 0.03085662
2 0.02250352 1 0.6933896 0.7102672 0.02774463
3 0.01125176 2 0.6708861 0.6947961 0.02752965
4 0.01000000 4 0.6483826 0.6877637 0.02743006
Node number 1: 2201 observations, complexity param=0.3066104
predicted class=No expected loss=0.323035
class counts: 1490 711
probabilities: 0.677 0.323
left son=2 (1731 obs) right son=3 (470 obs)
Primary splits:
Gender splits as LR, improve=199.821600, (0 missing)
Treatment splits as LR, improve=198.792000, (0 missing)
Class splits as RRLL, improve= 69.684100, (0 missing)
Age splits as RL, improve= 9.165241, (0 missing)
Surrogate splits:
Treatment splits as LR, agree=0.971, adj=0.864, (0 split)
Node number 2: 1731 observations, complexity param=0.01125176
predicted class=No expected loss=0.2120162
class counts: 1364 367
probabilities: 0.788 0.212
left son=4 (1667 obs) right son=5 (64 obs)
Primary splits:
Age splits as RL, improve=7.726764, (0 missing)
Treatment splits as LR, improve=7.726764, (0 missing)
Class splits as RLLL, improve=7.046106, (0 missing)
Surrogate splits:
Treatment splits as LR, agree=1, adj=1, (0 split)
Node number 3: 470 observations, complexity param=0.02250352
predicted class=Yes expected loss=0.2680851
class counts: 126 344
probabilities: 0.268 0.732
left son=6 (196 obs) right son=7 (274 obs)
Primary splits:
Class splits as RRLR, improve=50.015320, (0 missing)
Age splits as LR, improve= 1.197586, (0 missing)
Surrogate splits:
Age splits as LR, agree=0.619, adj=0.087, (0 split)
Node number 4: 1667 observations
predicted class=No expected loss=0.2027594
class counts: 1329 338
probabilities: 0.797 0.203
Node number 5: 64 observations, complexity param=0.01125176
predicted class=No expected loss=0.453125
class counts: 35 29
probabilities: 0.547 0.453
left son=10 (48 obs) right son=11 (16 obs)
Primary splits:
Class splits as RRL-, improve=12.76042, (0 missing)
Node number 6: 196 observations
predicted class=No expected loss=0.4591837
class counts: 106 90
probabilities: 0.541 0.459
Node number 7: 274 observations
predicted class=Yes expected loss=0.0729927
class counts: 20 254
probabilities: 0.073 0.927
Node number 10: 48 observations
predicted class=No expected loss=0.2708333
class counts: 35 13
probabilities: 0.729 0.271
Node number 11: 16 observations
predicted class=Yes expected loss=0
class counts: 0 16
probabilities: 0.000 1.000
と、分割の経緯が一緒に出力される。
== 決定木のチューニング ==
上のsummary()の出力の中に「complexity param」という項目が見られる。
rpart()では、この値の下限を指定することで、生成する決定木の深さを選択できる。
次の3行による学習結果の違いを観察してみよ。
print(rpart(Survived~.,data=titanic, control=c(cp=0.5)))
print(rpart(Survived~.,data=titanic, control=c(cp=0.05)))
print(rpart(Survived~.,data=titanic, control=c(cp=0.005)))
== 関係するグラフ ==
titanic.rpart <- rpart(Survived~., data=titanic)
plot(titanic.rpart)
text(titanic.rpart)
== 学習結果の解釈 ==
上の決定木の結果から、生存確率が高い条件は
- 1等・2等客室の女性乗客、または女性の乗組員
- 1等・2等客室の子供乗客
であり、生存確率が特に低いのは
- 3等客室の男性乗客
であったことが分かる。
=== 練習課題についての課題 ===
* 3つの分析手法を適用せよ。(時間内課題)
* 3つの分析手法の結果をまとめ、比較検討せよ。(レポート課題)
* これら以外に、Rで2値判別を行う手法を探し、適用して、比較に加えてみよ。(課外課題)
=== サポート欄 ===
* data.2のyがひとつ足りなかったのを、追加しました。(1122a)
* tic.learnというデータ名をtic.leaanとミスタイプしていたのを修正しました。(1122a)
* インターネットに繋がらないパソコンを使っている人は、TAさんから次の2つのファイルを貰ってください。(1122a)
* このWikiページのPDFファイル
* ticdata2000.txt
* 回帰分析の結果から標準化残差とテコ比の散布図を描くとき、配布資料では残差を標準化するのに、「残差の平方和を残差の自由度で割ったもの」を誤差分散の推定値としますが、Rでは「残差の標本分散」を誤差分散の推定値としています。第5章の例題では、それぞれ「残差平方和/7」と「残差平方和/9」ですので、Rが描く標準化残差のグラフはすべて、配布資料よりも9/7だけ原点から拡大されることになります。(0211p)
* V1は使わないのがおすすめ。番号の順序に意味がなく、各コードごとの頻度を集計させると、次のようになるため。(0237p)
> table(tic.learn$V1)
1 2 3 4 5 6 7 8 9 10 11 12 13 15 16 17 18 19 20
124 82 249 52 45 119 44 339 278 165 153 111 179 5 16 9 19 3 25
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
15 98 251 180 82 48 50 25 86 118 205 141 810 182 214 225 132 339 328
40 41
71 205
* 保険商品ごとに難易度が異なります。V86が一番簡単。
|V|0|1|2|3|4|5|6|
|V75|5426|382|14| | | | |
|V76|5529|173|100|11|8|1| |
|V77|5791|31| | | | | |
|V78|5784|38| | | | | |
|V79|5799|19|4| | | | |
|V80|2666|3017|126|7|3|2|1|
|V81|5819|3| | | | | |
|V82|5789|31|2| | | | |
|V83|5675|111|34|2| | | |
|V84|5777|44|1| | | | |
==== 参考 ====
=== V86も因子変数にしてみると ==
V86まで因子変数に変えると、glm関数の挙動が少し変わるかも。
tic.learn$V86f <- as.factor(tic.learn$V86)
=== 訪問客リストを作成したい場合 ===
== 考えたルールに基づく対象限定 ==
各変数に閾値を設けてルールを生成したとする。
たとえば、「V47が5.5以上かつV44が1未満」または「V47が5.5以上かつV1が{1,3,6,8,12,20}のどれか」、というルールは
次のように記す。
(tic.learn$V47>5.5 & tic.learn$V44<1) | (tic.learn$V47>5.5 & (tic.learn$V1==1 |tic.learn$V1==3 | tic.learn$V1==6 | tic.learn$V1==8 | tic.learn$V1==12 | tic.learn$V1==20) )
「&」が「かつ(AND)」、「|」が「または(OR)」である。
このルールを検証用データに適用するには、
tic.learn.visit <- (tic.learn$V47>5.5 & tic.learn$V44<1) | (tic.learn$V47>5.5 & (tic.learn$V1==1 |tic.learn$V1==3 | tic.learn$V1==6 | tic.learn$V1==8 | tic.learn$V1==12 | tic.learn$V1==20) )
と、訪問するか否かを二値(TRUE, FALSE)で表すオブジェクトを生成する。
このモデルに予測に基づいた訪問の成果を検証するには、訪問対象のリストtic.visitと検証用データの正解V86のクロス集計を行えばよい。
table(tic.learn.visit)
FALSE TRUE
3029 971
table(tic.learn.visit, tic.learn$V86)
tic.learn.visit 0 1
FALSE 2878 151
TRUE 884 87
ここでは、訪問対象に884+87=971人を選定し、そのうちの87人が実際に契約してくれる人だったことになる。
契約率は87/971=8.96%。また誤判別率は
(884+151)/4000
で25.9%となる。
== モデルに基づく対象限定 ==
学習したモデルに基づいて、訪問対象を狭めるには、predict()という関数を用いて、訪問対象か否かというリストを作成する。
まず、設定まで調整したモデルを、学習用データ(tic.learn)から得る。
tic.rpart <- rpart(V86~., data=tic.learn, control=c(cp=0.005))
次に、このモデル(ここではtic.rpart)を検証用データ(tic.eval)に適用して、契約してくれるか否かの予測を行う。
この際、0.05という閾値も調整の必要がある。
tic.eval.visit <- predict(tic.rpart, newdata=tic.eval)[,2]>0.05
このモデルに予測に基づいた訪問の成果を検証するには、訪問対象のリストtic.visitと検証用データの正解V86のクロス集計を行えばよい。
V86の保険商品が分析対象の場合は、
table(tic.eval.visit)
tic.eval.visit
FALSE TRUE
2389 1611
table(tic.eval.visit, tic.eval$V86)
tic.eval.visit 0 1
FALSE 2310 79
TRUE 1452 159
ここでは、訪問対象に1452+159=1611人を選定し、そのうちの159人が実際に契約してくれる人だったことになる。契約率は159/1452=11.0%。
また誤判別率は
(79+1452)/4000
で38.275%となる。
=== TICデータでロジスティック回帰を行う場合のメモ ===
== 想定される困難 ===
次の1行を実行すると、かなり時間がかかってエラーになる。
tic.glm.step <- step(glm(V86~., family="binomial", data=tic.learn)
次の４行、いずれもエラーになる。変数間の関係が悪すぎるよう。変数の意味を考えて、追加しないといけないかも。
tic.glm <- glm(V86~V1+V2+V3+V4+V5+V6+V7+V8+ V10+
V11+V12+ V14+ V16+V17+ V19+V20+
V21+V22+V23+ V25+V26+V27+V28+ V30+
V33+V34+V35+ V37+V38+V39+V40+
V42+V43+V44+V45+V46+V47+V48+V49+V50+
V51+V52+V53+V54+V55+V56+V57+V58+V59+V60+
V61+V62+V63+V64+V65+V66+V67+V68+V69+V70+
V71+V72+V73+V74+V75+V76+V77+V78+V79+V80+
V81+V82+V83+V84+V85, family="binomial", data=tic.learn)
table(predict(tic.glm, newdata=tic.eval)>0.5)
tic.glm <- glm(V86~ V2+V3+V4+ V6+V7+V8+ V10+
V11+V12+ V14+ V16+V17+ V19+V20+
V21+V22+V23+ V25+V26+V27+V28+ V30+
V33+V34+V35+ V37+V38+V39+V40+
V42+V43+V44+V45+V46+V47+V48+V49+V50+
V51+V52+V53+V54+V55+V56+V57+V58+V59+V60+
V61+V62+V63+V64+V65+V66+V67+V68+V69+V70+
V71+V72+V73+V74+V75+V76+V77+V78+V79+V80+
V81+V82+V83+V84+V85, family="binomial", data=tic.learn)
tic.glm <- glm(V86~V44+V45+V46+V47+V48+V49+V50+
V51+V52+V53+V54+V55+V56+V57+V58+V59+V60+
V61+V62+V63+V64+V65+V66+V67+V68+V69+V70+
V71+V72+V73+V74+V75+V76+V77+V78+V79+V80+
V81+V82+V83+V84+V85, family="binomial", data=tic.learn)
tic.glm <- glm(V86~V44+V45+V46+V47+V48+V49+V50+
V51+V52+V53+V54+V55+V56+V57+V58+V59+V60+
V61+V62+V63+V64, family="binomial", data=tic.learn)