今週の実験の内容は

• 回帰分析を用いたデータマイニング

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

• Rコマンダーを用いた回帰分析の2つの課題(自習、練習に相当)
• 解析データを用いた1つの課題(本番)

に取り組んで貰う。

1. 配付資料とRコマンダーを照らし合わせながら、出力される情報のどれが配付資料のどれに対応するのかを把握する
   • 回帰係数の推定値: Estimate
     • 切片: Intercept
   • 寄与率: R-Squared
   • 自由度調整済み寄与率
   • てこ比
   • 標準化残差
   • 変数増減法
2. 保険データの回帰分析、に取り組む (保険データの回帰分析)
   1. 回帰係数の推定
   2. 分散分析によるモデルの有意性の検討や回帰係数の有意性の検討
   3. てこ比や標準化残差などの検討
   4. 変数の増減
   5. 以上を繰り返す
3. 回帰分析の結果に基づいて、訪問する顧客層を絞り込む (訪問ルールの作成)
4. 必要に応じて、保険データの回帰分析と訪問ルールの作成を繰り返す

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からの抜粋と要約、の日本語版。

メモの確認用のコード。

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)

L0:分類を表す数字なので、大小関係に意味がなく、名義尺度である。そのままでは説明変数にならない。

L1:大きさが年齢の順なので、そのまま説明変数に使える。

L2:数字は分類を表すだけなので、連続尺度でも順序尺度でもなく、名義尺度。そのままでは説明変数にならない。

L3:順序尺度。このまま連続尺度の説明変数として用いる。

L4: 順序尺度。今回はこのまま連続尺度の変数として用いる。

kernlabパッケージに、加工済みのデータが入っていて、それを使うこともできる。

install.packages(c("kernlab"), dependencies=TRUE)
tic.learn <- ticdata[1:5822,]
tic.eval <- ticdata[5823:9822,]

• TICデータは、V65からV85までが21種類の保険商品の契約件数である。実験ペアのどちらかの学籍番号の末尾1桁を用いて、分析対象とする保険商品を決めよ。

• 重回帰分析を行い、分析結果を考察せよ。
• 重回帰分析の結果を用いて、指定された保険商品の契約に関する予測式を構築せよ。
• 線形モデルと重回帰分析を用いて、このデータをデータマイニングすることについて、考察せよ。

この課題ではMASSライブラリのみ、使う可能性がある。

library(MASS)

1つ目のデータは、Rに次の命令を実行させておく。

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)
data.1 <- data.frame(x=x,y=y)
rm(x,y)

２つ目のデータは、Rに次の命令を実行させておく。

x1 <- c(51,38,57,51,53,77,63,69,72,73)
x2 <- 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)
data.2 <- data.frame(x1=x1,x2=x2,y=y)
rm(x1,x2,y)

演習で用いる保険データは、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)

あとはRコマンダーで、data.1、data.2、tic.learnそれぞれについて、回帰分析を進める。

library(Rcmdr)

Rコマンダーでの回帰分析の手順は、次の通り。

[統計量] -> [モデルへの適合]
[モデル] -> [モデルを要約]
[モデル] -> [逐次モデル選択]
[モデル] -> [部分モデル選択]
[モデル] -> [仮説検定] -> [分散分析]
[モデル] -> [グラフ] -> [基本的診断プロット]
[モデル] -> [グラフ] -> [影響プロット]

この手順で分析を進めながら、参考資料の解析ストーリーと対比させよ。 「モデル選択」は、添付の資料に従うなら、変数増減法だが、その他のことも考えてよいし、 必ずしも選択されたモデルが最適であることもないので、少し変えても構わない。

レポート提出要領:下記「XXXXXXX」は各自の学籍番号(半角文字)で置き換えること

• 永田・棟近 (2001) 多変量解析法入門, サイエンス社

• 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が一番簡単。

以下の6行は、実行しない方がいい場合もある。

tic.learn$V1 <- as.factor(tic.learn$V1)
tic.learn$V5 <- as.factor(tic.learn$V5)
tic.learn$V86 <- as.factor(tic.learn$V86)
tic.eval$V1 <- as.factor(tic.eval$V1)
tic.eval$V5 <- as.factor(tic.eval$V5)
tic.eval$V86 <- as.factor(tic.eval$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%となる。

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