差分

このページの2つのバージョン間の差分を表示します。

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
mselab:2012:stat:week2:r2 [2012/12/11 10:33] – [実験の流れ] watalumselab:2012:stat:week2:r2 [不明な日付] (現在) – 外部編集 (不明な日付) 127.0.0.1
行 1: 行 1:
 +===== 統計工学 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&#データの説明|データの説明]]に、変数名を参考までに追記しました。
 +
 +
 === 概要 === === 概要 ===
  
行 5: 行 16:
   * データマイニング：保険のデータのV86の予測モデルの構築   * データマイニング：保険のデータのV86の予測モデルの構築
     * 回帰分析     * 回帰分析
-    * ロジスティック回帰 
     * 決定木     * 決定木
 +    * ロジスティック回帰 (オプション)
  
 である。3週間の流れがデータ分析(第1週)を行って、回帰分析等によるモデル構築(第2週)、そして第3週の他の手法も組み合わせたモデル構築と結論へと繋がるため、今週は[[http://kyoumu.office.uec.ac.jp/syllabus/2012/21/21_17123216.html|多変量解析]]との内容の重複を避けずにおいた。また同科目の履修を前提としていないため、回帰分析の学習も自習内容として含めてある。 である。3週間の流れがデータ分析(第1週)を行って、回帰分析等によるモデル構築(第2週)、そして第3週の他の手法も組み合わせたモデル構築と結論へと繋がるため、今週は[[http://kyoumu.office.uec.ac.jp/syllabus/2012/21/21_17123216.html|多変量解析]]との内容の重複を避けずにおいた。また同科目の履修を前提としていないため、回帰分析の学習も自習内容として含めてある。
  
   * 先週に続いてデータの把握、特にV86を中心に。   * 先週に続いてデータの把握、特にV86を中心に。
-  * Rコマンダーを用いた回帰分析の2つの課題(自習、練習に相当)+  * 回帰分析の2つの課題、同じ課題を決定木も(自習、練習に相当)
   * 解析データを用いた1つの課題(本番)   * 解析データを用いた1つの課題(本番)
  
 に取り組んで貰う。 に取り組んで貰う。
 +==== 実験の流れ ====
  
-=== 実験の流れ ===+時間内課題：回帰分析(のたぶん復習)と決定木分析 
 + 
 +  * 回帰分析 
 +    * 配付資料とこのページにあるコードを参考に、単回帰(説明変数1つ)と重回帰(説明変数2つ)の場合を自習して、回帰分析の結果の読み方について学ぶ。(時間内は少なくとも、分析結果やグラフをすべて得て、Wordに貼るなどしてお持ち帰りするところまで) 
 +    * 重回帰(説明変数2つ)の場合で、さらに変数を分類尺度に変換し、分類尺度の説明変数を含む場合の回帰分析の結果の読み方(特に推定値Estimateの意味や解釈)について考える。 
 +  * 決定木分析 
 +    * 回帰分析と同じデータを、決定木分析にかけてみて、どのようなモデルかを考える。 
 + 
 +課題：保険会社の顧客データのデータマイニング (時間内はできるところまで) 
 + 
 +  * まずは分類数の多い変数(V1とV5)の数を減らす。 
 +  * 回帰分析と決定木分析を行う。 
 +  * 得たモデルを考察する。 
 + 
 +本実験に関係するメモ
  
   - 配付資料とRの実行結果を見比べながら、lm関数を用いた回帰分析で出力される情報のどれが、配付資料のどれに対応するのかを把握する   - 配付資料とRの実行結果を見比べながら、lm関数を用いた回帰分析で出力される情報のどれが、配付資料のどれに対応するのかを把握する
行 32: 行 58:
     * てこ比: Leverage     * てこ比: Leverage
     * 標準化残差: Studentized Residuals     * 標準化残差: Studentized Residuals
-    * 変数増減法 
   - rpart関数を用いた決定木分析の出力を見ながら、決定木分析はどのようなモデルを構築するのかを検討する   - rpart関数を用いた決定木分析の出力を見ながら、決定木分析はどのようなモデルを構築するのかを検討する
   - 保険データの回帰分析、に取り組む (保険データの回帰分析)   - 保険データの回帰分析、に取り組む (保険データの回帰分析)
行 42: 行 67:
     - 以上を繰り返す     - 以上を繰り返す
   - 保険のデータの決定木分析、に取り組む (時間外課題)   - 保険のデータの決定木分析、に取り組む (時間外課題)
 +  - オプション課題
 +    - 回帰分析で変数選択を行う
 +      - 手動でP値を見ながらの変数減少法
 +      - AICを用いた変数増減法 stepAIC()
 +      - 決定木の学習パラメータの調整
 +
 +
 +==== 時間内課題１：回帰分析の自習 ====
 +=== 単回帰分析 ===
 +
 +配付資料の表4.1のデータの入力は次の通り。
 +<code>
 +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)
 +</code>
 +
 +回帰分析の実行にはlm関数を用いる。
 +<code>
 +lm(y~x, data=table.4.1)
 +</code>
 +「y~x」は、変数yを変数xで説明する回帰分析を行うための表現。これは
 +<jsmath>
 +Y = \beta_0+\beta_1 X+\epsilon
 +</jsmath>
 +という回帰式を推定せよ、という意味である。
 +
 +回帰分析の詳細を示すには、lm関数の実行結果を別の変数lm.4.1などに代入してから、詳細を出力させるsummary関数を用いる。
 +
 +<code>
 +lm.4.1 <- lm(y~x, data=table.4.1)
 +print(lm.4.1)
 +summary(lm.4.1)
 +</code>
 +
 +**問１：このsummary関数の出力結果を、配付資料と対比させて理解せよ。**
 +
 +回帰分析の結果を診断するために、幾つかのグラフを出力する機能がある。lm関数の実行結果をplot関数にかければ、そのようなグラフが4枚、順次出力される。
 +<code>
 +plot(test.lm)
 +</code>
 +
 +**問２：このplot関数の出力結果を、配付資料と対比させて理解せよ。**
 +
 +=== 重回帰分析 ===
 +
 +表5.1のデータは次のように入力する。
 +
 +<code>
 +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)
 +</code>
 +
 +回帰分析の実行には、説明変数が複数あっても(重回帰でも)lm関数を用いる。
 +<code>
 +lm(y~x.1+x.2, data=table.5.1)
 +</code>
 +今度は説明変数が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関数を用いる。
 +
 +<code>
 +lm.5.1 <- lm(y~x.1+x.2, data=table.5.1)
 +print(lm.5.1)
 +summary(lm.5.1)
 +</code>
 +
 +**問４：このplot関数の出力結果を、配付資料と対比させて理解せよ。**
 +
 +=== 重回帰応用(水準変数が説明変数に含まれる場合) ===
 +
 +上のデータを、広さを「広め(w)」「狭め(n)」とし、築年数も「新しめ(new)」「古め(old)」にする
 +<code>
 +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)
 +</code>
 +
 +**問５：このデータで、回帰分析を行い、グラフも作成して、先のモデルと何が異なるか、検討せよ。**
 +
 +<code>
 +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)
 +</code>
 +
 +==== 時間内課題２：決定木分析の自習 ====
 +
 +まず、次の一行を実行しておく。
 +
 +<code>
 +library(mvpart)
 +</code>
 +
 +上の実行例で「lm」とあるところをすべて「rpart」で置き換える。
 +
 +<code>
 +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)
 +</code>
 +
 +<code>
 +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)
 +</code>
 +
 +<code>
 +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)
 +</code>
 +
 +**問６：これらがどのようなモデルか、データと出力に照らして検討せよ。** グラフと画面出力の比較をまず行うとよい。また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|リンク先]]を参考に、ロジスティック回帰にも取り組んでみよ。(重回帰と決定木がとりあえずの目標)
 +
 +前の班との違い
 +
 +  * モデルの学習(推定、構築)は、学習用データのみでいい
 +  * 分類変数のグループ化はできるだけ、行った方がいい
 +  * モデルを得るだけでいい
 +
 +やらなくて良いこと。
 +
 +  * 検証用データでの検証
 +  * 訪問対象の明確化 (モデルのみでいいから)
  
 === データの説明 === === データの説明 ===
行 56: 行 228:
 [[http://kdd.ics.uci.edu/databases/tic/dictionary.txt|dictionary.txt]]からの抜粋と要約、の日本語版。 [[http://kdd.ics.uci.edu/databases/tic/dictionary.txt|dictionary.txt]]からの抜粋と要約、の日本語版。
  
-|変数|分類|メモ|+^変数^分類^メモ^
 |V1|顧客分類２|L0でコード化されている、数字の大きさに意味なし| |V1|顧客分類２|L0でコード化されている、数字の大きさに意味なし|
 |V2|住居数|大きいほど住む箇所が多い| |V2|住居数|大きいほど住む箇所が多い|
行 89: 行 261:
 </code> </code>
  
 +== 変数詳細 ==
 +
 +^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|
  
 == 各変数のコーディング == == 各変数のコーディング ==
行 94: 行 355:
 L0:分類を表す数字なので、大小関係に意味がなく、名義尺度である。そのままでは説明変数にならない。 L0:分類を表す数字なので、大小関係に意味がなく、名義尺度である。そのままでは説明変数にならない。
  
-|Value|Label|+^Value^Label^
 |1|High Income, expensive child| |1|High Income, expensive child|
 |2|Very Important Provincials| |2|Very Important Provincials|
行 184: 行 445:
 |8|f 10.000 - 19.999| |8|f 10.000 - 19.999|
 |9|f 20.000 - ?| |9|f 20.000 - ?|
 +
 +== データのダウンロードと読み込み ==
 +
 +演習で用いる保険データは、大学内からであれば、Rに次の命令を実行すれば読み込める。
 +<code>
 +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)
 +</code>
 +
 +学外もしくは自宅等であれば、最初の1行(Sys.setenvで始まる行)は不要。
  
 == 参考 == == 参考 ==
行 194: 行 470:
 tic.eval <- ticdata[5823:9822,] tic.eval <- ticdata[5823:9822,]
 </code> </code>
 +=== 準備 ===
  
-=== 今回の課題 === +上のデータの読み込みを実行済みであれば、この課題ではMASSライブラリを使う可能性があるので、それを読み込んでおく。
-== 概要 ==+
  
-  * 回帰分析について学ぶ。 
-  * TICデータは、V65からV86までが22種類の保険商品の契約件数である。今回はV86の契約に漕ぎ着けるためのモデル構築、を目的とする。 
-  * 重回帰分析を行い、分析結果を考察せよ。 
-  * 重回帰分析の結果を用いて、指定された保険商品の契約に関する予測式を構築せよ。 
-  * 線形モデルと重回帰分析を用いて、このデータをデータマイニングすることについて、考察せよ。 
- 
-=== 時間内課題１：回帰分析の自習 === 
-== 単回帰分析 == 
- 
-配付資料の表4.1のデータの入力は次の通り。 
 <code> <code>
-x <- c(2.2,4.1,5.5,1.9,3.4,2.6,4.2,3.7,4.9,3.2) +library(MASS)
-y <- c(71,81,86,72,77,73,80,81,85,74) +
-table.4.1 <- data.frame(x=x,y=y)+
 </code> </code>
  
-回帰分析の実行にはlm関数を用いる。+=== 回帰分析 === 
 +== 回帰分析をするだけ == 
 + 
 +V1からV85を説明変数にして、V86を予測するだけなら、次の命令をコピーして貼り付けたら、パラメータの最小二乗推定値は得られる。 
 <code> <code>
-lm(y~x, data=table.4.1)+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)
 </code> </code>
-「y~x」は、変数yを変数xで説明する回帰分析を行うための表現。これは 
-<jsmath> 
-Y = \beta_0+\beta_1 X+\epsilon 
-</jsmath> 
-という回帰式を推定せよ、という意味である。 
  
-回帰分析の詳細を示すには、lm関数の実行結果を別の変数lm.4.1などに代入してから、詳細を出力させるsummary関数を用いる。+先ほどと同様に、と回帰分析を実行するには、次のようにすれば、各種統計量も診断用のグラフも出力される。
  
 <code> <code>
-lm.4.1 <- lm(y~x, data=table.4.1) +lm.learn.0 <- lm(V86~V1 +V2 +V3 +V4 +V5 +V6 +V7 +V8 +V9 +V10 
-print(lm.4.1) +                +V11+V12+V13+V14+V15+V16+V17+V18+V19+V20 
-summary(lm.4.1)+                +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)
 </code> </code>
  
-**問１：このsummary関数の出力結果を、配付資料と対比させて理解せよ。**+ただ、これではV1やV5の扱いを「数値」としてしまっているし、有意でない(*や**や***がついてない)変数も多くモデルに取り込まれてしまっている。 
 +そこで、データマイニング手法を「線形モデル」に限るならば、モデルに対して次の2つの対応をとることが必要になる。 
 + 
 +  - 予測の役に立たない不要な変数を削る 
 +  - 本来は「分類コード」として扱うべき変数を、「数値データ」として扱っている現状を解消する 
 + 
 +それぞれ順に説明する。 
 + 
 +== 予測に不要な変数を削る == 
 + 
 +変数増減法という方法がある。上で得たlm.learn.0という「すべての変数を説明変数に加えた回帰モデル」に対して、
  
-回帰分析の結果を診断するために、幾つかのグラフを出力する機能がある。lm関数の実行結果をplot関数にかければ、そのようなグラフが4枚、順次出力される。 
 <code> <code>
-plot(test.lm)+lm.learn.AIC <- stepAIC(lm.learn.0, direction="both")
 </code> </code>
  
-**問２：このplot関数の出力結果を、配付資料と対比させて理解せよ。**+との1行を実行すると、AICという基準に照らして、不要な変数を1つずつ取り除か加えるか、検討してくれて、一番不要な変数から順に出し入れすることで、最適なモデルに到達するよう頑張ってくれる。
  
-== 重回帰分析 ==+=== 分類コード(水準, 因子変数)の扱い ===
  
-表5.1のデータは次のように入力する。+V1とV5の中身は、前述のように分類コードである。前回の課題で、ヒストグラムを描こうとすると、エラーが発生したかもしれない。 
 +Rでは、中身が分類コードや水準である変数のことを、因子変数(factor)と呼ぶ。 
 +V1とV5を因子変数に変換するには、次の2行を実行して、V1とV5を「分類尺度」に変換した変数をそれぞれ、V1fおよびV5fと名付ける。
  
 <code> <code>
-x.1 <- c(51,38,57,51,53,77,63,69,72,73) +tic.learn$V1f <- as.factor(tic.learn$V1) 
-x.2 <- c(16,4,16,11,4,22,5,5,2,1) +tic.learn$V5f <- as.factor(tic.learn$V5)
-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)+
 </code> </code>
  
-回帰分析の実行には、説明変数が複数あっても(重回帰でも)lm関数を用いる。+通常の数値変数(numeric)と因子変数(factor)の違いは、以下、幾つかのグラフを描いているので、1行ずつ貼って実行していくと、違いが分かる。 
 + 
 +V1やV5のヒストグラムは、hist関数を用いて
 <code> <code>
-lm(y~x.1+x.2, data=table.5.1)+hist(tic.learn$V1) 
 +hist(tic.learn$V5) 
 +</code> 
 +などで描ける。一方、V1fやV5fのヒストグラムは 
 +<code> 
 +hist(tic.learn$V1f) 
 +hist(tic.learn$V5f) 
 +</code> 
 +を実行するとエラーになる。これを得るには、table関数で集計してから、barplot関数で棒グラフを描く。 
 +<code> 
 +barplot(table(tic.learn$V1f)) 
 +barplot(table(tic.learn$V5f)
 </code> </code>
-今度は説明変数が2つあるので、「y~x.1+x.2」となる。これは +更にカテゴリごとの契約の有無の割合を見るには 
-<jsmath> +<code> 
-Y = \beta_0+\beta_1 X_1 + \beta_2 X_2+\epsilon +barplot(table(tic.learn$V86,tic.learn$V1f)) 
-</jsmath> +barplot(table(tic.learn$V86,tic.learn$V5f)) 
-という回帰式を推定せよ、という意味である。 +</code> 
- +とする。色の濃い部分が契約無、色の薄い部分が契約有である。
-回帰分析の詳細を示すには、lm関数の実行結果を別の変数lm.5.1などに代入してから、詳細を出力させるsummary関数を用いる。+
  
 +参考までに、因子変数を用いた回帰分析を実行するには、
 <code> <code>
-lm.5.1 <- lm(y~x.1+x.2, data=table.5.1) +lm.learn.1 <- lm(V86~V1f+V2 +V3 +V4 +V5f+V6 +V7 +V8 +V9 +V10 
-print(lm.5.1) +                +V11+V12+V13+V14+V15+V16+V17+V18+V19+V20 
-summary(lm.5.1)+                +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)
 </code> </code>
  
-**問４：このplot関数の出力結果を、配付資料と対比させて理解せよ。**+とすればよい。ただ、これではV5fについての推定ができず、まだまだ問題があることがわかる。 
 + 
  
-== 重回帰応用(水準変数が説明変数に含まれる場合) ==+== 分類変数のグループ化 == 
 + 
 +推定に問題が発生したり、カテゴリが多すぎることを解消するには、V1とV86のtable関数による集計結果や、V5とV86の集計結果を参考に、V1とV5について「幾つかのカテゴリをまとめてしまうグルーピング」という操作を行うことが1案である。 
 +回帰分析でも決定木分析でも数千のサンプルで40以上の分類は、これだけ変数があると細かすぎる。 
 + 
 +まず集計方法を再掲する。
  
-上のデータを、広さを「広め(w)」「狭め(n)」とし、築年数も「新しめ(new)」「古め(old)」にする 
 <code> <code>
-x.1 <- c("n","n","n","n","n","w","w","w","w","w") +library(mvpart) 
-x.2 <- c("old","new","old","old","new","old","new","new","new","new") +table(tic.learn$V1f, tic.learn$V86) 
-y <- c(3.0,3.2,3.3,3.9,4.4,4.5,4.5,5.4,5.4,6.0) +table(tic.learn$V5f, tic.learn$V86)
-table.5.1.c <- data.frame(x.1=x.1,x.2=x.2,y=y)+
 </code> </code>
  
-**問５：このデータで、回帰分析を行い、グラフも作成して、先のモデルと何が異なるか、検討せよ。** +これらをグラフに描くには
 <code> <code>
-lm.5.1.c <- lm(y~x.1+x.2, data=table.5.1.c) +barchart(table(tic.learn$V1f, tic.learn$V86)) 
-print(lm.5.1.c) +barchart(table(tic.learn$V5f, tic.learn$V86))
-summary(lm.5.1.c) +
-plot(lm.5.1.c)+
 </code> </code>
 +とする。
  
-=== 時間内課題２：決定木分析の自習 ===+さて、V86の予測にV1のどのカテゴリが不要かを見るために、決定木分析を行ってみる。 
 +決定木(実際には回帰木)を用いて、V86の増加に寄与しそうなカテゴリのみを残す。 
 +目論見は次の2点。
  
-まず、次の一行を実行しておく。+  * 契約にあまり寄与しないカテゴリには、コード0を割り当ててしまおう。 
 +  * 契約にしそうなカテゴリは残そう。
  
 <code> <code>
 library(mvpart) library(mvpart)
 +rpart(V86~V1f, data=tic.learn, control=rpart.control(cp=0))
 </code> </code>
  
-上の実行例で「lm」とあるところをすべて「rpart」で置き換える。+剪定せずにそのまま出力した回帰樹は次の通り。
  
 <code> <code>
-rpart.4.1 <- rpart(y~x, data=table.4.1) +n= 5822  
-print(rpart.4.1) + 
-summary(rpart.4.1) +node), split, n, deviance, yval 
-rpart.5.1 <- rpart(y~x.1+x.2, data=table.5.1) +      * denotes terminal node 
-print(rpart.5.1) + 
-summary(rpart.5.1) +  1) root 5822 327.1989000 0.05977327   
-rpart.5.1.c <- rpart(y~x.1+x.2, data=table.5.1.c) +    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   
-print(rpart.5.1.c) +      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   
-summary(rpart.5.1.c)+        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 *
 </code> </code>
  
-**問６：これらがどのようなモデルか、データと出力に照らして検討せよ。** 特にlm.5.1、lm.5.1.c、rpart.5.1、rpart.5.1.cの間の違いは考察せよ。+上の決定木の結果を見るなどして、下の「V1の値をV1grに代入するだけのコード」の右辺をいろいろ変えてみて、「V1の値の種類を減らす工夫」(グルーピング)をする。 
 +これをメモ帳に貼り付けて幾つかを0にしたり、幾つかの変数に同じ値を割り振るなど、が実際の作業となる。 
 +勿論、グルーピングの結果はレポートに記すこと。
  
-=== レポート課題 ===+<code> 
 +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) 
 +</code>
  
-この課題ではMASSライブラリのみ、使う可能性がある。+V5についても同様。
  
 <code> <code>
-library(MASS)+rpart(V86~V5f, data=tic.learn, control=rpart.control(cp=0))
 </code> </code>
  
-1つ目のデータは、Rに次の命令を実行させておく。+V5の値も幾つかのグループにまとめる。 
 <code> <code>
-x <- c(2.2,4.1,5.5,1.9,3.4,2.6,4.2,3.7,4.9,3.2) +tic.learn$V5gr <- 0 
-y <- c(71,81,86,72,77,73,80,81,85,74) +tic.learn$V5gr[tic.learn$V5==1] <- 1 
-data.1 <- data.frame(x=x,y=y) +tic.learn$V5gr[tic.learn$V5==2] <- 2 
-rm(x,y)+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)
 </code> </code>
  
-２つ目のデータは、Rに次の命令を実行させておく。+参考までに、V1とV5で回帰樹を作ってみると・・・。 
 <code> <code>
-x1 <- c(51,38,57,51,53,77,63,69,72,73) +rpart(V86~V1f+V5f, data=tic.learn, control=rpart.control(cp=0))
-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)+
 </code> </code>
  
-演習で用いる保険データは、Rに次の命令を実行させておく。+保険データへの決定木分析の適用に際しては、[[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を定義していなければエラーが表示されるだけ) 
 <code> <code>
-Sys.setenv("http_proxy"="http://130.153.8.66:8080/") +lm.learn.2 <- lm(V86~V1gr+V2 +V3 +V4 +V5gr+V6 +V7 +V8 +V9 +V10 
-tic.learn <- read.table("http://kdd.ics.uci.edu/databases/tic/ticdata2000.txt") +                +V11+V12+V13+V14+V15+V16+V17+V18+V19+V20 
-tic.eval <- read.table("http://kdd.ics.uci.edu/databases/tic/ticeval2000.txt") +                +V21+V22+V23+V24+V25+V26+V27+V28+V29+V30 
-tic.test <- read.table("http://kdd.ics.uci.edu/databases/tic/tictgts2000.txt") +                +V31+V32+V33+V34+V35+V36+V37+V38+V39+V40 
-tic.eval <- cbind(tic.eval, tic.test) +                +V41+V42+V43+V44+V45+V46+V47+V48+V49+V50 
-colnames(tic.eval)[86] <- "V86" +                +V51+V52+V53+V54+V55+V56+V57+V58+V59+V60 
-rm(tic.test)+                +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)
 </code> </code>
  
-あとはRコマンダーで、data.1、data.2、tic.learnそれぞれについて、回帰分析を進める。+=== 決定木分析 === 
 +== とりあえず決定木分析を実行してみる ==
  
 <code> <code>
-library(Rcmdr)+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)
 </code> </code>
  
-=== Rコマンダーで回帰分析をする際に用いるメニュー ===+== 決定木分析をいろいろチューニングする (オプション課題) ==
  
-Rコマンダーでの回帰分析の手順は、次の通り。+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とすると、結構、複雑な決定木に成長する。
 <code> <code>
-[統計量] -> [モデルへの適合] +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)
 </code> </code>
  
-この手順で分析を進めながら、参考資料の解析ストーリーと対比させよ。 +これはきっと、成長させすぎ・・・。複雑さの下限cpをさらに0.0001と小さくすると，より複雑な木が得られる。 
-「モデル選択」は、添付の資料に従うなら、変数増減法だが、その他のことも考えてよいし、 +グラフもノード数が増えるので，字が小さくなる。 
-必ずしも選択されたモデルが最適であることもないので、少し変えても構わない。 + 
-=== レポート ===+ 
 +==== レポート提出について ====
  
 レポート提出要領:下記「XXXXXXX」は各自の学籍番号(半角文字)で置き換えること レポート提出要領:下記「XXXXXXX」は各自の学籍番号(半角文字)で置き換えること
  
 ^項目^指定^ ^項目^指定^
-|提出期限|実験実施の翌週の火曜日の午前10時30分まで|+|提出期限|実験実施の翌週の月曜日の午後7時0分まで (今回は12月17日の午後7時)|
 |提出方法|電子メールに添付 (宛先は配付資料に記載)| |提出方法|電子メールに添付 (宛先は配付資料に記載)|
 |ファイル形式|Wordファイル (LaTeXで作成する場合は、dvipdfmxでPDFに変換すること)| |ファイル形式|Wordファイル (LaTeXで作成する場合は、dvipdfmxでPDFに変換すること)|
行 386: 行 861:
 |提出部数|レポートは各自1通ずつ。{{:mselab:report-header-2012.doc|レポートの表紙}}に、共同実験者の学籍番号と氏名を記すこと。| |提出部数|レポートは各自1通ずつ。{{:mselab:report-header-2012.doc|レポートの表紙}}に、共同実験者の学籍番号と氏名を記すこと。|
  
-=== 参考文献 ===+==== 参考文献 ====
  
   * 永田・棟近 (2001) [[http://www.saiensu.co.jp/?page=book_details&ISBN=ISBN978-4-7819-0980-6&YEAR=2001|多変量解析法入門]], サイエンス社   * 永田・棟近 (2001) [[http://www.saiensu.co.jp/?page=book_details&ISBN=ISBN978-4-7819-0980-6&YEAR=2001|多変量解析法入門]], サイエンス社
 +
 +
 +==== 練習課題(参考) ====
 +
 +決定木について、理解が難しい人は、この部分の決定木の箇所だけでも実行して考えてみるといい。この箇所は今回の課題ではない。
 +
 +=== タイタニック号 ===
 +
 +タイタニック号の乗客の生死のデータがある。Rで
 +<code>
 +Titanic
 +</code>
 +と実行すると、表示される。
 +
 +これを個別のレコードに展開し、更に救出の優先順位を高く設定された女性もしくは子供と、低めに設定された成人男性という変数も加える。
 +<code>
 +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)")
 +))
 +</code>
 +
 +さらに、数値に変換したデータも用意する。これは、線形学習機械とロジスティック学習機械のために用いる。
 +<code>
 +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)
 +</code>
 +
 +|変数|数値化情報|
 +|性別(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 というデータに含まれる変数の一覧を
 +<code>
 +names(titanic)
 +</code>
 +で取り出す。すると
 +<code>
 +[1] "Class"     "Gender"    "Age"       "Survived"  "Treatment"
 +</code>
 +のように5個の変数があることが分かる。
 +
 +これを用いて、生存(Survived)についての集計を、客室等級(Class)、性別(Gender)、年齢層(Age)の組み合わせで行うには、次の1行を実行する。
 +「$」の前がデータ名、「$」の後ろに変数名(フィールド名)を付ける。
 +<code>
 +ftable(titanic$Gender,titanic$Age,titanic$Class, titanic$Survived)
 +</code>
 +すると
 +<code>
 +                    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
 +</code>
 +と出力されるのを、各自確認せよ。
 +
 +同じことは、数値化したデータを用いても得られる。
 +<code>
 +ftable(titanic.2$Gender,titanic.2$Age,titanic.2$Class, titanic.2$Survived)
 +</code>
 +
 +上の出力と下の出力を見比べて、各変数の値の数値化を確認せよ。
 +<code>
 +         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
 +</code>
 +
 +=== 線形学習機械 ===
 +
 +上のデータ(titanic.2)に対して
 +<code>
 +lm(Survived~., data=titanic.2)
 +</code>
 +を実行すると、全変数を用いた線形学習機械が最小二乗法により、学習される。学習結果のみなら、この一行で
 +<code>
 +Call:
 +lm(formula = Survived ~ ., data = titanic.2)
 +
 +Coefficients:
 +(Intercept)       Gender          Age        Class  
 +    0.11725      0.46493      0.09655      0.05029  
 +</code>
 +と出力される。
 +
 +この結果に、統計的推測の結果を付与するなら
 +<code>
 +summary(lm(Survived~., data=titanic.2))
 +</code>
 +を実行すればよい。summary()の出力は次のように得られる。
 +<code>
 +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 
 +</code>
 +
 +== 線形学習機械のチューニング ==
 +
 +AICを用いた説明変数の選択を行うのであれば、
 +<code>
 +library(MASS)
 +</code>
 +とMASSパッケージを読み込んでから使える関数stepAIC()を用いて、
 +<code>
 +stepAIC(lm(Survived~., data=titanic.2))
 +</code>
 +を実行する。
 +<code>
 +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  
 +</code>
 +これはAICを変数増減法で最適化する関数であり、最終的に採用するモデルを出力してくれる便利な関数である。
 +
 +== 関係するグラフ ==
 +<code>
 +titanic.lm <- lm(Survived~., data=titanic.2)
 +par(mfrow=c(2,2))
 +plot(titanic.lm)
 +par(mfrow=c(1,1))
 +</code>
 +
 +== 学習結果の解釈 ==
 +
 +上の回帰分析から
 +<code>
 +Survived = 0.11725 + 0.46493 * Gender + 0.09655 * Age + 0.05029 * Class
 +</code>
 +という学習結果が得られた。Survivedを大きくする(=生存する)ために、Genderが0よりは1の方が良いことは、回帰係数 0.46493 を見れば分かる。
 +同様に Age も0よりは1の方が良く、Classも1よりは4の方が良い。
 +そしてどの変数も有意であることから、生存するためには、
 +
 +  * 男性よりは女性
 +  * 大人よりは子供
 +  * 客室等級は高めの部屋を選ぶ
 +
 +のが好条件となることが分かる。
 +
 +=== ロジスティック線形学習機械 ===
 +
 +生存確率をpとして
 +<code>
 +log(p/1-p)
 +</code>
 +を与える学習機械が、ロジスティック回帰モデルである。Rではglm()という関数を用いて
 +<code>
 +glm(Survived~., data=titanic.2, family="binomial")
 +</code>
 +で、ロジスティック回帰モデルを最尤推定で学習させることができる。
 +結果は
 +<code>
 +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 
 +</code>
 +と、lm()と似た出力を得る。
 +これもsummary()を加えて
 +<code>
 +summary(glm(Survived~., data=titanic.2, family="binomial"))
 +</code>
 +を実行すると、
 +<code>
 +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
 +</code>
 +と、種々の検定統計量が一緒に出力される。
 +
 +== ロジスティック学習機械のチューニング ==
 +
 +AICを用いた説明変数の選択を行うのであれば、
 +<code>
 +library(MASS)
 +</code>
 +とMASSパッケージを読み込んでから使える関数stepAIC()を用いて、
 +<code>
 +stepAIC(glm(Survived~., data=titanic.2, family="binomial"))
 +</code>
 +を実行する。
 +<code>
 +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 
 +</code>
 +これはAICを変数増減法で最適化する関数であり、最終的に採用するモデルを出力してくれる便利な関数である。
 +
 +== 関係するグラフ ==
 +
 +<code>
 +titanic.glm <- glm(Survived~., data=titanic.2, family="binomial")
 +par(mfrow=c(2,2))
 +plot(titanic.glm)
 +par(mfrow=c(1,1))
 +</code>
 +
 +== 学習結果の解釈 ==
 +
 +上の回帰分析から
 +<code>
 +log(Pr[Survived]/(1-Pr[Survived])) = -1.86224 + 2.050802 * Gender + 0.51147 * Age + 0.27834 * Class
 +</code>
 +という学習結果が得られた。Pr[Survived](生存する確率)を大きくするために、Genderが0よりは1の方が良いことは、回帰係数 2.050802 を見れば分かる。
 +同様に Age も0よりは1の方が良く、Classも1よりは4の方が良い。
 +そしてどの変数も有意であることから、生存するためには、
 +
 +  * 男性よりは女性
 +  * 大人よりは子供
 +  * 客室等級は高めの部屋を選ぶ
 +
 +のが好条件となることが分かる。
 +
 +=== 決定木 ===
 +
 +決定木は生存確率pの高低を際立たせるような、データの分割を表現するモデルである。Rではrpartパッケージの中のrpart()、もしくはmvpartパッケージの中のmvpart()という関数を用いて
 +<code>
 +library(mvpart)
 +rpart(Survived~.,data=titanic)
 +</code>
 +で、決定木を学習させることができる。
 +結果は
 +<code>
 +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) *
 +</code>
 +と、これまでとは違ったものになる。各行はノードと呼ばれる分割単位を表し、
 +たとえば
 +<code>
 + 1) root 2201 711 No (0.6769650 0.3230350)  
 +</code>
 +は、
 +  * データ全体(root)はレコード数が2201
 +  * このノードの代表値(一番多い値)はNo
 +  * 2201のうち711はNoでない
 +  * YesとNoの比は0.6769650:0.3230350
 +を意味する。これを性別で分割すると
 +<code>
 +   2) Gender=Male 1731 367 No (0.7879838 0.2120162)  
 +   3) Gender=Female 470 126 Yes (0.2680851 0.7319149)  
 +</code>
 +を得た、とある。これは
 +  * データ全体を性別で分割すると、生存確率がもっとも差が開く
 +  * 男性の代表値はNo(死亡)、0.7879838は死亡確率を意味するので、生存確率は0.2120162
 +  * 女性の代表値はYes(生存)、0.7319149は生存確率
 +と読み取ることができる。
 +生存確率の高いノードを探すと
 +<code>
 +     7) Class=1st,2nd,Crew 274  20 Yes (0.0729927 0.9270073) *
 +</code>
 +が目に付く。これは、
 +  - Gender=Female
 +  - Class=1st,2nd,Crew
 +と条件がふたつついたノードであり、「女性かつ１等客室の乗客、２等客室の乗客、または乗組員」となる。
 +この属性を持つ人々は、生存確率が0.9270073と最も高い。
 +逆に最も生存確率が低いのは
 +<code>
 +      10) Class=3rd 48  13 No (0.7291667 0.2708333) *
 +</code>
 +であり、このノードまでを
 +  - Gender=Male
 +  - Age=Adult
 +  - Class=3rd
 +と辿れることから、「男性で成人で３等客室の乗客」は生存確率が0.2708333と最も低かったことが分かる。
 +
 +これもsummary()を加えて
 +<code>
 +summary(rpart(Survived~.,data=titanic))
 +</code>
 +を実行すると、
 +<code>
 +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 
 +</code>
 +と、分割の経緯が一緒に出力される。
 +== 決定木のチューニング ==
 +
 +上のsummary()の出力の中に「complexity param」という項目が見られる。
 +rpart()では、この値の下限を指定することで、生成する決定木の深さを選択できる。
 +
 +次の3行による学習結果の違いを観察してみよ。
 +<code>
 +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)))
 +</code>
 +
 +
 +== 関係するグラフ ==
 +<code>
 +titanic.rpart <- rpart(Survived~., data=titanic)
 +plot(titanic.rpart)
 +text(titanic.rpart)
 +</code>
 +
 +== 学習結果の解釈 ==
 +
 +上の決定木の結果から、生存確率が高い条件は
 +  - 1等・2等客室の女性乗客、または女性の乗組員
 +  - 1等・2等客室の子供乗客
 +であり、生存確率が特に低いのは
 +  - 3等客室の男性乗客
 +であったことが分かる。
 +
 +=== 練習課題についての課題 ===
 +
 +  * 3つの分析手法を適用せよ。(時間内課題)
 +  * 3つの分析手法の結果をまとめ、比較検討せよ。(レポート課題)
 +  * これら以外に、Rで2値判別を行う手法を探し、適用して、比較に加えてみよ。(課外課題)
  
 === サポート欄 === === サポート欄 ===
行 421: 行 1402:
 |V84|5777|44|1| | | | | |V84|5777|44|1| | | | |
  
-=== 参考 === +==== 参考 ==== 
-== 少し加工する ==+=== V86も因子変数にしてみると ==
  
-以下の6行は、実行しない方がいい場合もある。+V86まで因子変数に変えると、glm関数の挙動が少し変わるかも。
  
 <code> <code>
-tic.learn$V1 <- as.factor(tic.learn$V1) +tic.learn$V86f <- as.factor(tic.learn$V86)
-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)+
 </code> </code>
  
-あとはそのまま。 +=== 訪問客リストを作成したい場合 ===
 == 考えたルールに基づく対象限定 == == 考えたルールに基づく対象限定 ==
  
行 518: 行 1493:
 で38.275%となる。 で38.275%となる。
  
-== 想定される困難 ==+=== TICデータでロジスティック回帰を行う場合のメモ === 
 +== 想定される困難 ===
  
 次の1行を実行すると、かなり時間がかかってエラーになる。 次の1行を実行すると、かなり時間がかかってエラーになる。
行 565: 行 1541:
 V61+V62+V63+V64, family="binomial", data=tic.learn) V61+V62+V63+V64, family="binomial", data=tic.learn)
 </code> </code>
 +