差分

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

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
mselab:2012:stat:week2:r2 [2012/12/11 10:42] – [概要] 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&#データの説明|データの説明]]に、変数名を参考までに追記しました。
 +
 +
 === 概要 === === 概要 ===
  
行 11: 行 22:
  
   * 先週に続いてデータの把握、特にV86を中心に。   * 先週に続いてデータの把握、特にV86を中心に。
-  * 回帰分析の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の契約に漕ぎ着けるためのモデル構築、を目的とする。 
-  * 重回帰分析を行い、分析結果を考察せよ。 
-  * 重回帰分析の結果を用いて、指定された保険商品の契約に関する予測式を構築せよ。 
-  * 線形モデルと重回帰分析を用いて、このデータをデータマイニングすることについて、考察せよ。 
-  * 余裕があったら、[[http://stat.inf.uec.ac.jp/dokuwiki/doku.php?id=mselab:2012:stat:week3|リンク先]]を参考に、ロジスティック回帰にも取り組んでみよ。(重回帰と決定木がとりあえずの目標) 
- 
-前の班との違い 
- 
-  * モデルの学習(推定、構築)は、学習用データのみでいい 
-  * 分類変数のグループ化はできるだけ、行った方がいい 
-  * モデルを得るだけでいい 
- 
-やらなくて良いこと。 
- 
-  * 検証用データでの検証 
-  * 訪問対象の明確化 (モデルのみでいいから) 
- 
-=== 時間内課題１：回帰分析の自習 === 
-== 単回帰分析 == 
- 
-配付資料の表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関数を用いる。 +=== 回帰分析 === 
-<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関数を用いる。+V1からV85を説明変数にして、V86を予測するだけなら、次の命令をコピーして貼り付けたら、パラメータの最小二乗推定値は得られる。
  
 <code> <code>
-lm.4.1 <- lm(y~x, data=table.4.1) +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)
 </code> </code>
  
-**問１：このsummary関数の出力結果を、配付資料と対比させて理解せよ。**+先ほどと同様に、と回帰分析を実行するには、次のようにすれば、各種統計量も診断用のグラフも出力される。
  
-回帰分析の結果を診断するために、幾つかのグラフを出力する機能がある。lm関数の実行結果をplot関数にかければ、そのようなグラフが4枚、順次出力される。 
 <code> <code>
-plot(test.lm)+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)
 </code> </code>
  
-**問２：このplot関数の出力結果を、配付資料と対比させて理解せよ。**+ただ、これではV1やV5の扱いを「数値」としてしまっているし、有意でない(*や**や***がついてない)変数も多くモデルに取り込まれてしまっている。 
 +そこで、データマイニング手法を「線形モデル」に限るならば、モデルに対して次の2つの対応をとることが必要になる。
  
-== 重回帰分析 ==+  - 予測の役に立たない不要な変数を削る 
 +  - 本来は「分類コード」として扱うべき変数を、「数値データ」として扱っている現状を解消する
  
-表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関数を用いる。 +変数増減法という方法がある。上で得たlm.learn.0という「すべての変数を説明変数に加えた回帰モデル」に対して、
-<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> <code>
-lm.5.1 <- lm(y~x.1+x.2, data=table.5.1) +lm.learn.AIC <- stepAIC(lm.learn.0, direction="both")
-print(lm.5.1) +
-summary(lm.5.1)+
 </code> </code>
  
-**問４：このplot関数の出力結果を、配付資料と対比させて理解せよ。**+との1行を実行すると、AICという基準に照らして、不要な変数を1つずつ取り除か加えるか、検討してくれて、一番不要な変数から順に出し入れすることで、最適なモデルに到達するよう頑張ってくれる。
  
-== 重回帰応用(水準変数が説明変数に含まれる場合) ==+=== 分類コード(水準, 因子変数)の扱い ===
  
-上のデータを、広さを「広め(w)」「狭め(n)」とし、築年数も「新しめ(new)」「古め(old)」にする +V1とV5の中身は、前述のように分類コードである。前回の課題で、ヒストグラムを描こうとすると、エラーが発生したかもしれない。 
-<code> +Rでは、中身が分類コードや水準である変数のことを、因子変数(factor)と呼ぶ。 
-x.1 <- c("n","n","n","n","n","w","w","w","w","w") +V1とV5を因子変数に変換するには、次の2行を実行して、V1とV5を「分類尺度」に変換した変数をそれぞれ、V1fおよびV5fと名付ける。
-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> <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) 
-rpart.5.1 <- rpart(y~x.1+x.2, data=table.5.1) 
-print(rpart.5.1) 
-summary(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) 
-</code> 
- 
-**問６：これらがどのようなモデルか、データと出力に照らして検討せよ。** 特にlm.5.1、lm.5.1.c、rpart.5.1、rpart.5.1.cの間の違いは考察せよ。 
- 
-=== レポート課題 === 
- 
-この課題ではMASSライブラリのみ、使う可能性がある。 
- 
-<code> 
-library(MASS) 
-</code> 
- 
-1つ目のデータは、Rに次の命令を実行させておく。 
-<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) 
-data.1 <- data.frame(x=x,y=y) 
-rm(x,y) 
-</code> 
- 
-２つ目のデータは、Rに次の命令を実行させておく。 
-<code> 
-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) 
-</code> 
- 
-演習で用いる保険データは、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> 
- 
-あとは、上の回帰分析と同様に行えばよいのだが、V1とV5が特殊なので、これらを例に「分類変数のグループ化」について、少し手順を解説する。 
- 
-=== 分類変数のグループ化 === 
- 
- 
-  * 分類尺度の変数の、分類数を減らす (回帰分析でも決定木分析でも数千のサンプルで40以上の分類は細かすぎる) 
- 
-V1とV5を「分類尺度」の変数に変換し、それぞれV1fおよびV5fと名付ける。 
- 
-<code> 
-tic.learn <- read.table("http://kdd.ics.uci.edu/databases/tic/ticdata2000.txt") 
 tic.learn$V1f <- as.factor(tic.learn$V1) tic.learn$V1f <- as.factor(tic.learn$V1)
 tic.learn$V5f <- as.factor(tic.learn$V5) tic.learn$V5f <- as.factor(tic.learn$V5)
 </code> </code>
  
-以下、幾つかのグラフを描いているので、1行ずつ貼って実行していくことを進める。+通常の数値変数(numeric)と因子変数(factor)の違いは、以下、幾つかのグラフを描いているので、1行ずつ貼って実行していくと、違いが分かる。
  
 V1やV5のヒストグラムは、hist関数を用いて V1やV5のヒストグラムは、hist関数を用いて
行 401: 行 564:
 </code> </code>
 とする。色の濃い部分が契約無、色の薄い部分が契約有である。 とする。色の濃い部分が契約無、色の薄い部分が契約有である。
 +
 +参考までに、因子変数を用いた回帰分析を実行するには、
 +<code>
 +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)
 +</code>
 +
 +とすればよい。ただ、これではV5fについての推定ができず、まだまだ問題があることがわかる。
 +
 +
 +
 +== 分類変数のグループ化 ==
 +
 +推定に問題が発生したり、カテゴリが多すぎることを解消するには、V1とV86のtable関数による集計結果や、V5とV86の集計結果を参考に、V1とV5について「幾つかのカテゴリをまとめてしまうグルーピング」という操作を行うことが1案である。
 +回帰分析でも決定木分析でも数千のサンプルで40以上の分類は、これだけ変数があると細かすぎる。
 +
 +まず集計方法を再掲する。
  
 <code> <code>
行 408: 行 598:
 </code> </code>
  
 +これらをグラフに描くには
 <code> <code>
-rpart(V86~V1f, data=tic.learn, control=rpart.control(cp=0)) +barchart(table(tic.learn$V1f, tic.learn$V86)) 
-<code>+barchart(table(tic.learn$V5f, tic.learn$V86)) 
 +</code> 
 +とする。
  
- V1の値を幾つかのグループにまとめる+さて、V86の予測にV1のどのカテゴリが不要かを見るために、決定木分析を行ってみる。 
 +決定木(実際には回帰木)を用いて、V86の増加に寄与しそうなカテゴリのみを残す。 
 +目論見は次の2点。
  
   * 契約にあまり寄与しないカテゴリには、コード0を割り当ててしまおう。   * 契約にあまり寄与しないカテゴリには、コード0を割り当ててしまおう。
   * 契約にしそうなカテゴリは残そう。   * 契約にしそうなカテゴリは残そう。
- 
-決定木(実際には回帰木)を用いて、V86の増加に寄与しそうなカテゴリのみを残す。 
  
 <code> <code>
行 497: 行 690:
 </code> </code>
  
-下記は、V1の値をV1grに代入するだけのコード。これで幾つかを0にしたり、幾つかの変数に同じ値を割り振るなどして、「V1grの値の種類を減らす」工夫をする。+上の決定木の結果を見るなどして、下の「V1の値をV1grに代入するだけのコード」の右辺をいろいろ変えてみて、「V1の値の種類を減らす工夫」(グルーピング)をする。 
 +これをメモ帳に貼り付けて幾つかを0にしたり、幾つかの変数に同じ値を割り振るなど、が実際の作業となる。 
 +勿論、グルーピングの結果はレポートに記すこと。
  
 <code> <code>
行 551: 行 746:
 </code> </code>
  
- V5の値も幾つかのグループにまとめる。+V5の値も幾つかのグループにまとめる。
  
 <code> <code>
行 576: 行 771:
 保険データへの決定木分析の適用に際しては、[[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|このページ]]が参考になるかもしれない。 保険データへの決定木分析の適用に際しては、[[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> 
 +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) 
 +</code> 
 + 
 +=== 決定木分析 === 
 +== とりあえず決定木分析を実行してみる == 
 + 
 +<code> 
 +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> 
 + 
 +== 決定木分析をいろいろチューニングする (オプション課題) == 
 + 
 +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> 
 +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> 
 + 
 +これはきっと、成長させすぎ・・・。複雑さの下限cpをさらに0.0001と小さくすると，より複雑な木が得られる。 
 +グラフもノード数が増えるので，字が小さくなる。 
 + 
 + 
 +==== レポート提出について ====
  
 レポート提出要領:下記「XXXXXXX」は各自の学籍番号(半角文字)で置き換えること レポート提出要領:下記「XXXXXXX」は各自の学籍番号(半角文字)で置き換えること
  
 ^項目^指定^ ^項目^指定^
-|提出期限|実験実施の翌週の月曜日の午後7時0分まで|+|提出期限|実験実施の翌週の月曜日の午後7時0分まで (今回は12月17日の午後7時)|
 |提出方法|電子メールに添付 (宛先は配付資料に記載)| |提出方法|電子メールに添付 (宛先は配付資料に記載)|
 |ファイル形式|Wordファイル (LaTeXで作成する場合は、dvipdfmxでPDFに変換すること)| |ファイル形式|Wordファイル (LaTeXで作成する場合は、dvipdfmxでPDFに変換すること)|
行 588: 行 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値判別を行う手法を探し、適用して、比較に加えてみよ。(課外課題)
  
 === サポート欄 === === サポート欄 ===
行 623: 行 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>
  
-あとはそのまま。 +=== 訪問客リストを作成したい場合 ===
 == 考えたルールに基づく対象限定 == == 考えたルールに基づく対象限定 ==
  
行 720: 行 1493:
 で38.275%となる。 で38.275%となる。
  
-== 想定される困難 ==+=== TICデータでロジスティック回帰を行う場合のメモ === 
 +== 想定される困難 ===
  
 次の1行を実行すると、かなり時間がかかってエラーになる。 次の1行を実行すると、かなり時間がかかってエラーになる。
行 767: 行 1541:
 V61+V62+V63+V64, family="binomial", data=tic.learn) V61+V62+V63+V64, family="binomial", data=tic.learn)
 </code> </code>
 +