mselab:2012:stat:week2:r2

統計工学 2週目

はじめに

連絡 2012.12.11

自分で頑張って、とお願いした、回帰分析と決定木分析のコードを追記しました。
このページを実験時間中に改訂しましたが、自習部分その1と自習部分その2の内容はそのままです。決定木に関して、少しグラフの出力などを、追加しました。
TICデータの重回帰分析のコードと、決定木分析のコードを追記しました。分類コードの扱いはほぼ今朝のまま。重回帰分析のコードの中に変数選択がありますが、これは参考です。
練習課題(参考)はあくまでも参考まで、です。参考も次週の内容を含んでいるので、参考まで、です。
データの説明に、変数名を参考までに追記しました。

概要

今週の実験の内容は

データマイニング：保険のデータのV86の予測モデルの構築
- 回帰分析
- 決定木
- ロジスティック回帰 (オプション)

である。3週間の流れがデータ分析(第1週)を行って、回帰分析等によるモデル構築(第2週)、そして第3週の他の手法も組み合わせたモデル構築と結論へと繋がるため、今週は多変量解析との内容の重複を避けずにおいた。また同科目の履修を前提としていないため、回帰分析の学習も自習内容として含めてある。

先週に続いてデータの把握、特に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関数を用いた決定木分析の出力を見ながら、決定木分析はどのようなモデルを構築するのかを検討する
保険データの回帰分析、に取り組む (保険データの回帰分析)
1. 回帰係数の推定
2. 変数の加工(今回初)
3. 分散分析によるモデルの有意性の検討や回帰係数の有意性の検討
4. てこ比や標準化残差などの検討
5. 変数の増減 (stepAIC関数を使う?)
6. 以上を繰り返す
保険のデータの決定木分析、に取り組む (時間外課題)
オプション課題
1. 回帰分析で変数選択を行う
  1. 手動でP値を見ながらの変数減少法
  2. AICを用いた変数増減法 stepAIC()
  3. 決定木の学習パラメータの調整

時間内課題１：回帰分析の自習

単回帰分析

配付資料の表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で説明する回帰分析を行うための表現。これは <jsmath> Y = \beta_0+\beta_1 X+\epsilon </jsmath> という回帰式を推定せよ、という意味である。

回帰分析の詳細を示すには、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」となる。これは <jsmath> Y = \beta_0+\beta_1 X_1 + \beta_2 X_2+\epsilon </jsmath> という回帰式を推定せよ、という意味である。

回帰分析の詳細を示すには、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の契約に漕ぎ着けるためのモデル構築、を目的とする。
重回帰分析を行い、分析結果を考察せよ。
重回帰分析の結果を用いて、指定された保険商品の契約に関する予測式を構築せよ。
線形モデルと重回帰分析を用いて、このデータをデータマイニングすることについて、考察せよ。
余裕があったら、リンク先を参考に、ロジスティック回帰にも取り組んでみよ。(重回帰と決定木がとりあえずの目標)

前の班との違い

モデルの学習(推定、構築)は、学習用データのみでいい
分類変数のグループ化はできるだけ、行った方がいい
モデルを得るだけでいい

やらなくて良いこと。

検証用データでの検証
訪問対象の明確化 (モデルのみでいいから)

データの説明

TIC2000

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のエリアに持ち家があることを意味する。職業、社会層などもすべて、該当するエリアの箇所が郵便番号の一桁目で埋まっている。

変数

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))

保険データへの決定木分析の適用に際しては、このページやこのページが参考になるかもしれない。

ここまで行うと、回帰分析を実行する関数は、次のようになる。(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通ずつ。レポートの表紙に、共同実験者の学籍番号と氏名を記すこと。

参考文献

永田・棟近 (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
<none>                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
<none>        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)

目次

統計工学 2週目

はじめに

連絡 2012.12.11

概要

実験の流れ

時間内課題１：回帰分析の自習

単回帰分析

重回帰分析

重回帰応用(水準変数が説明変数に含まれる場合)

時間内課題２：決定木分析の自習

課題：保険会社の顧客データのデータマイニング

データの説明

TIC2000

変数

変数詳細

各変数のコーディング

データのダウンロードと読み込み

参考

準備

回帰分析

回帰分析をするだけ

予測に不要な変数を削る

分類コード(水準, 因子変数)の扱い

分類変数のグループ化

決定木分析

とりあえず決定木分析を実行してみる

決定木分析をいろいろチューニングする (オプション課題)

レポート提出について

参考文献

練習課題(参考)

タイタニック号

先週の課題と同等の課題

線形学習機械

線形学習機械のチューニング

関係するグラフ

学習結果の解釈

ロジスティック線形学習機械

ロジスティック学習機械のチューニング

関係するグラフ

学習結果の解釈

決定木

決定木のチューニング

関係するグラフ

学習結果の解釈

練習課題についての課題

サポート欄

参考

V86も因子変数にしてみると

訪問客リストを作成したい場合

考えたルールに基づく対象限定

モデルに基づく対象限定

TICデータでロジスティック回帰を行う場合のメモ

想定される困難