差分

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

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
mselab:2012:stat:week2:r2 [2012/12/11 13:57] – [分類コード(水準変数)の扱い] watalumselab:2012:stat:week2:r2 [不明な日付] (現在) – 外部編集 (不明な日付) 127.0.0.1
行 1: 行 1:
 ===== 統計工学 2週目 ===== ===== 統計工学 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&#データの説明|データの説明]]に、変数名を参考までに追記しました。
 +
 +
 === 概要 === === 概要 ===
  
行 13: 行 22:
  
   * 先週に続いてデータの把握、特にV86を中心に。   * 先週に続いてデータの把握、特にV86を中心に。
-  * 回帰分析の2つの課題(自習、練習に相当)+  * 回帰分析の2つの課題、同じ課題を決定木も(自習、練習に相当)
   * 解析データを用いた1つの課題(本番)   * 解析データを用いた1つの課題(本番)
  
 に取り組んで貰う。 に取り組んで貰う。
- 
 ==== 実験の流れ ==== ==== 実験の流れ ====
  
行 103: 行 111:
 **問２：このplot関数の出力結果を、配付資料と対比させて理解せよ。** **問２：このplot関数の出力結果を、配付資料と対比させて理解せよ。**
  
-== 重回帰分析 ==+=== 重回帰分析 ===
  
 表5.1のデータは次のように入力する。 表5.1のデータは次のように入力する。
行 134: 行 142:
 **問４：このplot関数の出力結果を、配付資料と対比させて理解せよ。** **問４：このplot関数の出力結果を、配付資料と対比させて理解せよ。**
  
-== 重回帰応用(水準変数が説明変数に含まれる場合) ==+=== 重回帰応用(水準変数が説明変数に含まれる場合) ===
  
 上のデータを、広さを「広め(w)」「狭め(n)」とし、築年数も「新しめ(new)」「古め(old)」にする 上のデータを、広さを「広め(w)」「狭め(n)」とし、築年数も「新しめ(new)」「古め(old)」にする
行 153: 行 161:
 </code> </code>
  
-=== 時間内課題２：決定木分析の自習 ===+==== 時間内課題２：決定木分析の自習 ====
  
 まず、次の一行を実行しておく。 まず、次の一行を実行しておく。
行 167: 行 175:
 print(rpart.4.1) print(rpart.4.1)
 summary(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) rpart.5.1 <- rpart(y~x.1+x.2, data=table.5.1)
 print(rpart.5.1) print(rpart.5.1)
 summary(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) rpart.5.1.c <- rpart(y~x.1+x.2, data=table.5.1.c)
 print(rpart.5.1.c) print(rpart.5.1.c)
 summary(rpart.5.1.c) summary(rpart.5.1.c)
 +plot(rpart.5.1.c)
 +text(rpart.5.1.c)
 </code> </code>
  
-**問６：これらがどのようなモデルか、データと出力に照らして検討せよ。** 特にlm.5.1、lm.5.1.c、rpart.5.1、rpart.5.1.cの間の違いは考察せよ。 +**問６：これらがどのようなモデルか、データと出力に照らして検討せよ。** グラフと画面出力の比較をまず行うとよい。またlm.5.1、lm.5.1.c、rpart.5.1、rpart.5.1.cの間の違いは考察せよ。
 ==== 課題：保険会社の顧客データのデータマイニング ==== ==== 課題：保険会社の顧客データのデータマイニング ====
  
行 209: 行 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|住居数|大きいほど住む箇所が多い|
行 242: 行 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|
  
 == 各変数のコーディング == == 各変数のコーディング ==
行 247: 行 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|
行 362: 行 470:
 tic.eval <- ticdata[5823:9822,] tic.eval <- ticdata[5823:9822,]
 </code> </code>
- 
 === 準備 === === 準備 ===
  
行 371: 行 478:
 </code> </code>
  
-=== 回帰分析をするだけ ===+=== 回帰分析 === 
 +== 回帰分析をするだけ ==
  
 V1からV85を説明変数にして、V86を予測するだけなら、次の命令をコピーして貼り付けたら、パラメータの最小二乗推定値は得られる。 V1からV85を説明変数にして、V86を予測するだけなら、次の命令をコピーして貼り付けたら、パラメータの最小二乗推定値は得られる。
行 412: 行 520:
 それぞれ順に説明する。 それぞれ順に説明する。
  
-=== 予測に不要な変数を削る ===+== 予測に不要な変数を削る ==
  
 変数増減法という方法がある。上で得たlm.learn.0という「すべての変数を説明変数に加えた回帰モデル」に対して、 変数増減法という方法がある。上で得たlm.learn.0という「すべての変数を説明変数に加えた回帰モデル」に対して、
行 477: 行 585:
  
  
-=== 分類変数のグループ化 ===+== 分類変数のグループ化 ==
  
-説明変数の選択とは別に、ticデータのV1とV5は分類コード(水準、分類尺度)の変数なので、これらをまずは分類コードに変換する。 +推定に問題が発生したり、カテゴリが多すぎることを解消するには、V1とV86のtable関数による集計結果や、V5とV86の集計結果を参考に、V1とV5について「幾つかのカテゴリをまとめてしまうグルーピング」という操作を行うことが1案である。 
- +回帰分析でも決定木分析でも数千のサンプルで40以上の分類は、これだけ変数があると細かすぎる。
- +
-「分類変数のグループ化」について、少し手順を解説する。 +
- +
-  * 分類尺度の変数の、分類数を減らす (回帰分析でも決定木分析でも数千のサンプルで40以上の分類は細かすぎる)+
  
 +まず集計方法を再掲する。
  
 <code> <code>
行 493: 行 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>
行 582: 行 690:
 </code> </code>
  
-下記は、V1の値をV1grに代入するだけのコード。これで幾つかを0にしたり、幾つかの変数に同じ値を割り振るなどして、「V1grの値の種類を減らす」工夫をする。+上の決定木の結果を見るなどして、下の「V1の値をV1grに代入するだけのコード」の右辺をいろいろ変えてみて、「V1の値の種類を減らす工夫」(グルーピング)をする。 
 +これをメモ帳に貼り付けて幾つかを0にしたり、幾つかの変数に同じ値を割り振るなど、が実際の作業となる。 
 +勿論、グルーピングの結果はレポートに記すこと。
  
 <code> <code>
行 636: 行 746:
 </code> </code>
  
- V5の値も幾つかのグループにまとめる。+V5の値も幾つかのグループにまとめる。
  
 <code> <code>
行 660: 行 770:
  
 保険データへの決定木分析の適用に際しては、[[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と小さくすると，より複雑な木が得られる。
 +グラフもノード数が増えるので，字が小さくなる。
 +
  
 ==== レポート提出について ==== ==== レポート提出について ====
行 676: 行 864:
  
   * 永田・棟近 (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値判別を行う手法を探し、適用して、比較に加えてみよ。(課外課題)
  
 === サポート欄 === === サポート欄 ===
行 708: 行 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>
  
-あとはそのまま。 +=== 訪問客リストを作成したい場合 ===
 == 考えたルールに基づく対象限定 == == 考えたルールに基づく対象限定 ==
  
行 805: 行 1493:
 で38.275%となる。 で38.275%となる。
  
-== 想定される困難 ==+=== TICデータでロジスティック回帰を行う場合のメモ === 
 +== 想定される困難 ===
  
 次の1行を実行すると、かなり時間がかかってエラーになる。 次の1行を実行すると、かなり時間がかかってエラーになる。
行 852: 行 1541:
 V61+V62+V63+V64, family="binomial", data=tic.learn) V61+V62+V63+V64, family="binomial", data=tic.learn)
 </code> </code>
 +