=== 概要 === 今週の実験の内容は * 回帰分析を用いたデータマイニング である。3週間の流れがデータ分析(第1週)を行って、回帰分析によるモデル構築(第2週)、そして第3週の他の手法を用いたモデル構築へと繋がるため、今週は[[http://kyoumu.office.uec.ac.jp/syllabus/2012/21/21_17123216.html|多変量解析]]との内容の重複を避けずにおいた。また同科目の履修を前提としていないため、回帰分析の学習も自習内容として含めてある。 * Rコマンダーを用いた回帰分析の2つの課題(自習、練習に相当) * 解析データを用いた1つの課題(本番) に取り組んで貰う。 === 実験の流れ === - 配付資料とRコマンダーを照らし合わせながら、出力される情報のどれが配付資料のどれに対応するのかを把握する * 回帰係数の推定値: Estimate * 切片: Intercept * 寄与率: R-Squared * 自由度調整済み寄与率 * てこ比 * 標準化残差 * 変数増減法 - 保険データの回帰分析、に取り組む (保険データの回帰分析) - 回帰係数の推定 - 分散分析によるモデルの有意性の検討や回帰係数の有意性の検討 - てこ比や標準化残差などの検討 - 変数の増減 - 以上を繰り返す - 回帰分析の結果に基づいて、訪問する顧客層を絞り込む (訪問ルールの作成) - 必要に応じて、保険データの回帰分析と訪問ルールの作成を繰り返す === データの説明 === == 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) == 各変数のコーディング == 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 - ?| == 参考 == kernlabパッケージに、加工済みのデータが入っていて、それを使うこともできる。 install.packages(c("kernlab"), dependencies=TRUE) tic.learn <- ticdata[1:5822,] tic.eval <- ticdata[5823:9822,] === 今回の課題 === == 概要 == * TICデータは、V65からV85までが21種類の保険商品の契約件数である。実験ペアのどちらかの学籍番号の末尾1桁を用いて、分析対象とする保険商品を決めよ。 |学籍番号の末尾1桁|解析する保険商品| |0|V75| |1|V76| |2|V77| |3|V78| |4|V79| |5|V80| |6|V81| |7|V82| |8|V83| |9|V84| * 重回帰分析を行い、分析結果を考察せよ。 * 重回帰分析の結果を用いて、指定された保険商品の契約に関する予測式を構築せよ。 * 線形モデルと重回帰分析を用いて、このデータをデータマイニングすることについて、考察せよ。 === 実験準備 === この課題では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コマンダーで回帰分析をする際に用いるメニュー === Rコマンダーでの回帰分析の手順は、次の通り。 [統計量] -> [モデルへの適合] [モデル] -> [モデルを要約] [モデル] -> [逐次モデル選択] [モデル] -> [部分モデル選択] [モデル] -> [仮説検定] -> [分散分析] [モデル] -> [グラフ] -> [基本的診断プロット] [モデル] -> [グラフ] -> [影響プロット] この手順で分析を進めながら、参考資料の解析ストーリーと対比させよ。 「モデル選択」は、添付の資料に従うなら、変数増減法だが、その他のことも考えてよいし、 必ずしも選択されたモデルが最適であることもないので、少し変えても構わない。 === レポート === レポート提出要領:下記「XXXXXXX」は各自の学籍番号(半角文字)で置き換えること ^項目^指定^ |提出期限|実験実施の翌週の火曜日の午前10時30分まで| |提出方法|電子メールに添付 (宛先は配付資料に記載)| |ファイル形式|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|多変量解析法入門]], サイエンス社 === サポート欄 === * 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| | | | | === 参考 === == 少し加工する == 以下の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)