===== データマイニング (2012年度) =====
大学院講義。

==== 演習の準備 2012.06.14 ====
[[http://stat.inf.uec.ac.jp/library/dm.2011/counting-processes-20100708.pdf|このファイル]](昨年の課題)に収められたコードがすべて動くような環境を整えて下さい。
具体的には

  * R言語 (PerlやCと似た文法でMATLABのような機能を持つオブジェクト指向言語、処理系はインタプリタでコンパイラは無し)
  * R言語のsurvivalパッケージ (最近はバンドルされてるはず)

の2つです。

=== はじめに ===
演習を行うのにコーディング量が最も少なくて済む環境として、R言語を用います。
ただし、使って貰うのにR言語をすべて理解する必要はなく、こちらが提供したコードを少し改変する範囲で演習課題が終わるような工夫をしています。
演習を進める上で、コードの改変では済まないけどやってみたいこと、があったら、教員に相談してください。

Rは、1976年に発表され、1980年頃に外部リリースが始まり、1988年にオブジェクト指向が取り入れられた、AT&Tのベル研究所で開発されたS言語というプログラミング言語の派生言語(文法や関数はほぼ共通だが、実装方針が異なるので、クローンとも言える)です。
S言語の実装のひとつなら、Rと呼ぶべきですが、世間に倣ってR言語と呼びます。
SもRも、線形代数などの基本的な数値計算には、世界中で古くから用いられて現在ではとても枯れている[[http://www.netlib.org/|Netlib]]に含まれているライブラリを用いています。
例えば[[http://www.netlib.org/blas/|BLAS]]は1979年、[[http://www.netlib.org/lapack|LAPACK]]は1992年まで、遡ることができますが、いずれもFortran 77で書かれています。
既存のライブラリが広く普及しているのに、新たに実装することで、既知のバグを新たに引き起こしたり、メンテナンスの労力を分散させるのを避けるためです。

R言語の今日の時点での最新版は、2012/03/30にリリースされた2.15.0です。この言語に関して、[[http://developer.r-project.org/|開発者向けのウェブサイト]]に、[[http://developer.r-project.org/devel-guidelines.txt|ガイドライン]]が掲載されています。そこを読むと、
<code>
DO fix simple bugs
DO NOT fix bugs that require extensive modification
DO NOT fix exotic bugs that haven't bugged anyone
DO make small enhancements if they are badly needed
DO NOT let one user decide what is badly needed
DO fix configuration for broken platforms
DO NOT break functioning platforms in the process

DO NOT fix things that are not broken
DO NOT restructure code
DO NOT add experimental code
DO NOT modify the API (unless absolutely sure it is buggy)
DO NOT change defaults without a *very* good reason

DO clarify documentation
ONLY add or modify examples if needed for clarification
DO NOT reword messages

DO NOT modify regression tests (except if they were buggy)
Be very careful in adding regression tests and consider using them
  in R-devel only at first.
</code>
と書かれています。メジャーリビジョンを変更するようなドラスティックな設計変更や大規模な改修は行わないこと、実験的なコードは含めないこと、不具合の修正と機能の追加をしつつ定期的にマイナーリビジョンを更新する、などの開発方針が読み取れます。よって、使用するパッケージが要求しているバージョン以降であれば、多少古くても大丈夫です。これまでのところ、全ユーザに対して、特定のバージョン以降を用いるように、という推奨は出されていません。

まずはR言語をインストールします。
インストールは通常、[[http://cran.r-project.org/|CRAN]]からダウンロードして行います。
Rの国際的な開発組織が運営しているウェブサイトです。ここからのダウンロードは、とりあえず信用しても良いです。
このサイトは、[[http://essrc.hyogo-u.ac.jp/cran/mirrors.html|世界中にミラーサイトを持っており]]、一番近いサイトにアクセスすることが推奨されています。
以下では、この大学から一番ダウンロード速度が速い、[[http://cran.ism.ac.jp/|統計数理研究所のミラーサイト]]にリンクを貼りました。

  * WindowsかMac OX Sなら、まず[[http://cran.ism.ac.jp/|CRAN(のミラーサイト)]]にアクセスして下さい。
    - Windowsなら、"Windows"の文字をクリックし、[[http://cran.ism.ac.jp/bin/windows/|次に現れるページ]]で"base"をクリックすると、"Download R X.Y.Z for Windows"という文字が見つかるので、それをクリックしてダウンロードし、実行する。
    - Mac OS Xなら、"Mac OS X"の文字をクリックし、[[http://cran.ism.ac.jp/bin/macosx/|次に現れるページ]]で "R-X.Y.Z.pkg (latest version)"という文字が見つかるので、それをクリックしてダウンロードし、実行する。
  * Linuxなら、次の2行のどちらかが、インストールしてくれます。sudoが必要かも。<code>
yum install R R-*
apt-get install r-base
</code>

次に、Rを起動して、必要なパッケージを追加します。今回は次の一つだけ。

  * survivalパッケージをインストールするには、次の一行を実行します。<code>
install.packages(pkgs=c("survival"), repos='http://cran.ism.ac.jp/')
</code>

この後、[[http://stat.inf.uec.ac.jp/library/dm.2011/counting-processes-20100708.pdf|このファイル]]にあるコードをすべて実行してみてください。問題無く実行でき、赤い文字などでエラー(Error)やワーニング(Warning)が表示されなければ、成功です。

もしワーニングやエラーが表示されたら、[[mailto:watalu.yamamoto_at_gmail.com|山本]]まで、表示されたエラーメッセージなど、Rの出力画面の情報と共に、連絡を下さい。(ただし、メールアドレスの_at_は"@"に取り替えて、送信してください。)

=== 既知の不具合 ===

== インターネットに接続できないのが原因 ==

この大学には、プライベートIPアドレスを持つPCが多くあります。それらには、NATルータやファイアウォールを介してインターネットに接続されているものと、NATルータなどは介さずにProxyサーバと呼ばれるサーバを経由してインターネットには接続せずに間接的に利用しているものがあります。以下は、この後者で発生する不具合への対応です。


プライベートIPアドレスを持ち、NATルータを介していないコンピュータからは、インターネットに直接アクセスすることはできず、PROXYを通す必要があります。
Windowsの場合は、Rを起動する際に
<code>
C:\Program Files\R\R-X.Y.Z\bin\Rgui.exe  (Xはバージョン、Yはメジャーリビジョン、Zはマイナーリビジョン)
</code>
となっているスタートメニュー内のショートカットの代わりに
<code>
"C:\Program Files\R\R-X.Y.Z\bin\Rgui.exe" --internet2
</code>
というショートカットを新たに作ると、
<code>
〔コントロールパネル〕→〔ネットワークとインターネット接続〕→〔インターネットオプション〕→〔接続〕→〔LANの設定〕→〔プロキシサーバー〕
</code>
で設定したプロキシを使用できます。

その他の環境については、[[http://appl.stat.inf.uec.ac.jp/doku.php?id=r:how_to:internet_proxy|r;how_to:internet_proxy]]が参考になるかもしれません。

== データファイルの保存方法が原因 ==

ダウンロードしたrats.datは、改行コードがUNIX用になっているため、一度、ブラウザで開いた上で、すべてをコピーし、メモ帳などに貼り付けてから、rats.datとして保存しないと、Rで読み込む際にエラーとなることがあります。

== PDFが原因 ==

LaTeXで書いたテキストをdvipdfmxでPDFに変換すると、なぜか「~」の記号のフォントが、置き換えられてしまいます。見た目はほとんど変わりませんが、半角の「~」では無くなっています。改行コードが余計につくこともあります。
テキストのコードをそのままRに貼り付けると、プログラムとして認識できない文字が入力されたことになりますので、「~」は半角文字で打ち直すようにして下さい。
あるいはPDFファイルではなく、このページからコピー＆ペーストして下さい。

課題1について、問題なく動作するコードを下に記します。

<code>
dimnames(rats)[[2]] <- c("id","start","stop","status","enum","trt")
rats$rtrunc <- rep(NA, nrow(rats))
rats$tstatus <- ifelse(rats$start == 0, 1, 2)
rats[1:12, c("id","start","stop","status","rtrunc","tstatus","enum","trt")]

NPfit <- coxph(Surv(start,stop,status)~1, data=rats, subset=(trt==0),method="breslow")
KM <- survfit(NPfit, type="aalen")
NA.MF <- data.frame(time=c(0,KM$time), na=-log(c(1,KM$surv)))
plot(NA.MF, type="l")

NPfit <- coxph(Surv(start,stop,status)~1, data=rats, subset=(trt==1),method="breslow")
KM <- survfit(NPfit, type="aalen")
NA.MF <- data.frame(time=c(0,KM$time), na=-log(c(1,KM$surv)))
lines(NA.MF, type="l")

NPfit <- coxph(Surv(start, stop, status)~trt, data=rats, method="breslow")
summary(NPfit)
KM <- survfit(NPfit, type="aalen")
NA.MF <- data.frame(time=c(0,KM$time), na=-log(c(1,KM$surv)))
plot(NA.MF)

NPfit <- coxph(Surv(start, stop, status)~trt+cluster(id), data=rats, method="breslow")
summary(NPfit)
KM <- survfit(NPfit, type="aalen")
NA.MF <- data.frame(time=c(0,KM$time), na=-log(c(1,KM$surv)))
</code>

== ライブラリをロードしていないのが原因 ==

<code>
library(survival)
</code>

==== 演習課題1 2012.06.14 ====
|出題|2012.06.14|
|〆切|2012.06.28|

レポートの提出について

|提出先|mailto:watalu.yamamoto_at_gmail.com (_at_を＠に交換してください)|
|提出方法|PDFファイルに変換してメールに添付して送付してください|
|提出確認|たまにまとめて返信します|



=== データの説明 ===
今回は Bladder Cancer Recurrences dataと呼ばれる、個人の再来院データを用いる。
観測されているのは、当初時点の、腫瘍の数(number)、最も大きい腫瘍の大きさ(size)、治療の種類(rx)、
それから履歴として、再来院(event)、及び間隔(start, stop)、である。

|変数|説明|データ|
|number|当初の腫瘍の数|整数|
|size|当初の最も大きい腫瘍の大きさ|整数|
|rx|治療の種類|整数 (0=プラセボ、1=pyrodixine or thiotepa)|
|event|事象の種類|整数 (0=打ち切り、1=来院)|
|start|区間の開始時点|整数(単位は月)|
|stop|区間の終了時点|整数(単位は月)|
|id|患者番号|整数|

オリジナルデータには、再来院時の腫瘍の数や、大きさなど、時間依存の変数も含まれているが、今回は上のデータのみから、
病気の進行ではなく、再来院の頻度についての分析を行う。

=== 準備 ===

Rで下記のコードを実行し、データ(bladder2)を少し変換したデータ(bladder.data)を作成する。

<code>
library(survival)
data(bladder)
bladder.data <- bladder2
bladder.data$rx[bladder.data$rx==1] <- 0
bladder.data$rx[bladder.data$rx==2] <- 1
</code>

最後に使うので、MASSパッケージも読み込んでおく。
<code>
library(MASS)
</code>

=== 治療群毎の強度関数の推定 ===

まず治療群(Thiotepaを用いている患者群)の計数過程の再来院確率関数を推定し、グラフに描く。
そしてコントロール群(プラセボを用いている患者群)の計数過程の生存関数を推定し、先のグラフに追記する。
ついでに両側95%信頼区間も描いた。
ここではKaplan-Meier推定量という方法を用いているが、詳細は省略する。

<code>
par(mfrow=c(2,1))
NPfit.1 <- coxph(Surv(start,stop,event)~1, data=bladder.data, subset=(rx==1),method="breslow")
KM.1 <- survfit(NPfit.1, conf.int=.95, type="aalen")
plot(KM.1)
NPfit.2 <- coxph(Surv(start,stop,event)~rx, data=bladder.data, subset=(rx==0),method="breslow")
KM.2 <- survfit(NPfit.2, conf.int=.95, type="aalen")
lines(KM.2, lty=3)
plot(KM.2)
lines(KM.1, lty=3)
</code>

95%信頼区間を重ね書きできないので、Treatment群の生存関数と信頼区間を描いてPlacebo群の生存関数を上書きしたものと、
Placebo群の生存関数と信頼区間を描いてTreatment群の生存関数を上書きしたもの、の2枚を用意した。

次に、上記と同等だが、それぞれの累積強度関数を推定してみる。
ここではNelson-Aalen推定量という推定方法を用いているが、詳細は省略する。

<code>
par(mfrow=c(1,1))
NPfit <- coxph(Surv(start,stop,event)~1, data=bladder.data, subset=(rx==1),method="breslow")
KM <- survfit(NPfit, conf.int=.95, type="aalen")
NA.MF <- data.frame(time=c(0,KM$time), na=-log(c(1,KM$surv)))
plot(NA.MF, type="l", lty=1,ylim=c(0,2.4), xlim=c(0,50))

NPfit <- coxph(Surv(start,stop,event)~rx, data=bladder.data, subset=(rx==0),method="breslow")
KM <- survfit(NPfit, conf.int=.95, type="aalen")
NA.MF <- data.frame(time=c(0,KM$time), na=-log(c(1,KM$surv)))
lines(NA.MF, type="l", lty=2)
</code>

== 問1 ==
これら2枚のグラフから、何が読み取れるか？

=== 強度関数の回帰分析 ===

上と同じ事を、比例ハザードモデルで推定する。

<code>
NPfit <- coxph(Surv(start, stop, event)~factor(rx), data=bladder.data, method="breslow")
summary(NPfit)
KM <- survfit(NPfit, type="aalen")
NA.MF <- data.frame(time=c(0,KM$time), na=-log(c(1,KM$surv)))
lines(NA.MF, type="l", lty=4)
</code>

出力は次のようになる。

<code>
Call:
coxph(formula = Surv(start, stop, event) ~ factor(rx), data = bladder.data, method = "breslow")

  n= 178 

                  coef       exp(coef) se(coef)          z Pr(>|z|)  
factor(rx)1 -0.3655    0.6939   0.1976 -1.85   0.0644 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

            exp(coef) exp(-coef) lower .95 upper .95
factor(rx)1    0.6939      1.441    0.4711     1.022

Rsquare= 0.02   (max possible= 0.994 )
Likelihood ratio test= 3.52  on 1 df,   p=0.0605
Wald test            = 3.42  on 1 df,   p=0.06438
Score (logrank) test = 3.46  on 1 df,   p=0.06293
</code>

それぞれ、直訳すると、次の通り。

<code>
Call:
coxph(formula = Surv(start, stop, event) ~ factor(rx), data = bladder.data, method = "breslow")

  n= 178 

            　　推定値  exp(推定値)   標準誤差           t検定のp値 (5%以下なら有意)  
factor(rx)1 -0.3655    0.6939     0.1976 -1.85   0.0644 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

                 exp(coef) exp(-coef) 下側95%   上側95%
factor(rx)1    0.6939      1.441    0.4711     1.022

重相関係数の平方= 0.02   (max possible= 0.994 )  (今回は無視)
尤度比検定統計量= 3.52  on 1 df,   p=0.0605   (有意になると、モデルを統計的に支持する)
ワルド検定統計量= 3.42  on 1 df,   p=0.06438   (有意になると、モデルを統計的に支持する)
ログランク検定統計量= 3.46  on 1 df,   p=0.06293   (有意になると、モデルを統計的に支持する)
</code>

== 問2 ==
以下のモデルから出発し、t検定のp値を眺めて、5%以上の変数を外したりしながら、各種検定統計量が有意になるようなモデルを探して下さい。
<code>
NPfit <- coxph(Surv(start, stop, event)~factor(rx)*number*size+cluster(id), 
               data=bladder.data, method="breslow")
summary(NPfit)
scatter.smooth(residuals(NPfit))
abline(h=0,lty=3,col=2)
</code>

表示されるグラフは、残差プロットです。
本来はこのグラフも眺めながら、分析を進めていくものですが、今回はモデルを得ながら、このグラフも眺めてみよ。
注目するのは、点のばらつき具合。

上のモデルはフルモデルと呼ばれ、すべての変数の効果、すべての交互作用の効果を含んでいるため、下記のように書くこともできる。

<code>
NPfit <- coxph(Surv(start, stop, event)~factor(rx)+number+size
               +factor(rx):number+factor(rx):size+number:size+factor(rx):number:size+cluster(id), 
               data=bladder.data, method="breslow")
summary(NPfit)
scatter.smooth(residuals(NPfit))
abline(h=0,lty=3,col=2)
</code>

こちらから出発して、"+変数"を削る方が作業は簡単でしょう。
なお推定結果にcluster(id)の係数は現れません。
"+cluster(id)"は、各変数の標準誤差を、モデルが正しくはないかもしれない、ということを考慮しながら推定する、サンドイッチ推定量という方法を用いて推定するために必要な項なので、
これだけは削らないこと。

また、上のモデル当てはめと同時に、先の強度関数の推定値に、ベースライン強度関数の推定値を重ねていくと、モデルを変えてもほとんど変化しないことが見て取れる。
これは、回帰係数とベースライン強度関数を別々に推定できている、Cox比例ハザードモデルの特徴でもある。

=== 強度関数の回帰分析 ===
AICというモデル選択基準は、これを最小にすると、データの予測精度の意味でのモデルの当てはまりが一番良い、という解釈を持つ。
上では、手動でフルモデルから変数を削ってもらったが、このAIC基準を用いて、自動的にモデル選択をさせることもできる。


<code>
NPfit <- coxph(Surv(start, stop, event)~factor(rx)*number*size+cluster(id), data=bladder.data, method="breslow")
stepAIC(NPfit)
</code>

== 問3 ==

上で得た、AICで最適なモデルに+cluster(id)を加え直し、問2で得たモデルと、ベースラインハザード関数、回帰係数、残差プロットなど、比較してみよ。

=== まとめ ===

このデータの再来院に関して、以上の分析から、何が言えるか、まとめよ。

==== 演習課題2 2012.07.26 ====
=== 日程 ===
|出題|2012.07.26|
|〆切|2012.08.10|


=== 今回の解析の流れ ===

  - 解析対象はSoftware Debugging Data (Dalal and McIntosh, 1994; Cook and Lawless, 2007, Sec. 1.2.2 and Appendix C)
  - まずはデータを、テスト時間 vs 累積発見不具合数、テスト時間 vs 累積発見不具合数の一階差分、テスト時間 vs 累積発見不具合数の二階差分、テスト時間 vs 累積改修行数、テスト時間 vs 累積改修行数の一階差分、など図示しながら、吟味していき、観測期間が変曲点を通り過ぎている（観測期間の終わりの方で累積のグラフの傾きが減少傾向にある）ことを確認する
  - 飽和する成長曲線モデルを当てはめる、今回は観測系列が1系列のみなので、非定常ポアソン過程を用いる (ロバスト分散などは考えないし、検定なども行わず、強度関数のモデルを工夫するだけ、という宣言)
  - 当てはめた累積強度関数とデータを重ねて描く、マルチンゲール残差、あるいはポアソン残差などの残差プロットを描く、などしつつAICも計算しておく (これらの残差プロットを、Cook and Lawlessでは「特にモデルに瑕疵があるとは思えない」としていることを、記しておく)
  - 上のモデルをベースに、モデルを少しずつ複雑にしてみて、どのようなモデルが適当か、調べていく

最初に当てはめる強度関数は

<jsmath>\lambda\left(t|H\left(t\right)\right)=\alpha\theta e^{-\theta t}+\beta\theta e^{-\theta t} \int_0^t \frac{dC\left(u\right)}{e^{-\theta u}}</jsmath>

である。この初項は飽和する曲線で、第二項は直近の改修が大きな影響を、古い改修ほど小さな影響を与えるように指数平滑化したような関数であり、改修行数を計数過程の強度関数に反映させるよう提案されたモデルである。

これは、期間 <jsm>[t_{j-1}, t_{j})</jsm> 内の改修がすべて <jsm> t_{j} </jsm> に追加されたとして、この期間内の強度関数は

<jsmath>\lambda\left(t|H\left(t\right)\right)=\theta e^{-\theta t}\left\{\alpha+\beta\sum_{l=1}^{j-1} c_l e^{\theta t_l}\right\}</jsmath>

と一定になる、と書き直すことができる。すると、この期間内の期待発見件数は

<jsmath>E\left[N_j|H_j\right] = \int_{t_{j-1}}^{t_j} \lambda\left(t|H\left(t\right)\right)dt = \left(e^{-\theta t_{j-1}}-e^{-\theta t_j}\right)\left(\alpha + \beta \sum_{l=0}^{j-1} c_l e^{\theta t_l}\right)</jsmath>

となる。これを <jsm>\mu_j</jsm> とおくと、ポアソン過程の尤度関数が

<jsmath>
L\left(\theta, \alpha, \beta\right)=\prod_{j=1}^k e^{-\mu_j}{\mu_j}^N_j
</jsmath>

であることから、

<jsmath>
\log L = \sum_{j=1}^k \left(-\mu_j+N_j \log \mu_j\right)
</jsmath>

を未知母数 <jsm> \theta, \alpha, \beta </jsm> について最大化すれば、良いことがわかる。その他、モデル選択に必要なAICは

<jsmath>
AIC\left(\hat{\theta}, \hat{\alpha}, \hat{\beta}\right)=-2\log L\left(\hat{\theta}, \hat{\alpha}, \hat{\beta}\right)+2k
</jsmath>

ただし<jsm>k</jsm>は母数の数である。残差は

<jsmath>
N_j-\hat{\mu}_j=N_j-\left(e^{-\hat{\theta} t_{j-1}}-e^{-\hat{\theta} t_j}\right)\left(\hat{\alpha} + \hat{\beta} \sum_{l=0}^{j-1} c_l e^{\hat{\theta} t_l}\right)
</jsmath>

もしくは、

<jsmath>
\left(N_j-\hat{\mu}_j\right)/\sqrt{\hat{\mu}_j}
</jsmath>

今回の解析では、Cook and Lawless (2007)では、上の強度関数をこのデータに対して当てはめているが、これを少しだけ精緻化することを試みてもらう。

下記では、

<jsmath>
\lambda\left(t|H\left(t\right)\right)=\theta e^{-\theta t}\left\{\alpha+\beta\sum_{l=1}^{d} c_l e^{\theta t_l}+\beta_2\sum_{l=d+1}^{j-1} c_l e^{\theta t_l}\right\}
</jsmath>

と、改修履歴に依存する部分を <jsm>l=1,\ldots,d</jsm> と <jsm>l=d+1,\ldots,j-1</jsm> のように前半と後半に分け、それぞれに母数 <jsm>\beta, \beta_2</jsm> と付与してみる例を与えてある。<jsm>d</jsm>の与え方は、現時点から300時間遡った時点まで、と最初から300時間、の二種類の与え方は例示している。

これを更に拡張するなど、できたら試みて欲しいとは思っている。


=== 準備 ===

データファイルは http://stat.inf.uec.ac.jp/library/dm.2011/software.txt に置いてある。そこからダウンロードしても良いし、下記のコードを実行してRに直接読み込ませることもできる。
<code>
software.debugging <- read.table("http://stat.inf.uec.ac.jp/library/dm.2011/software.txt",
                           sep=",", header=T)
</code>

先週、紹介した通り、このデータは3つのフィールドから成る。

|t|Nt|Ct|
|テスト時間(人時間)|累積不具合発見数|累積改修行数|
|0.0|0|0|
|4.8|0|16012|
|6.0|0|16012|
|以下略|||

今回の演習では、尤度関数を自分で書き、それをRのoptim()という最適化の関数を用いて最大化し、パラメータを推定する。
その途中で、平滑化のグラフを作成するのに、gamパッケージを用いるので、最初にインストールしておく。

<code>
install.packages(pkgs=c("gam"), 
                 repos='http://cran.md.tsukuba.ac.jp/',
                 dependencies=TRUE)
library(gam)
</code>

=== テスト用コード ===

下記のコードを走らせて、エラーが出なければ、

  * データの読み込み
  * 必要なライブラリのインストール

が完了していて、準備は整っていることが確認できる。

<code>
plot(software.debugging$t, software.debugging$Nt,
     type="l",
     xlab="Staff Days", ylab="Cumulative Faults")
library(gam)
software.debugging.gam <- gam(Nt~-1+bs(t),
                            data=software.debugging)
plot(software.debugging.gam,
     xlab="Staff Days", ylab="Cumulative Faults")
points(software.debugging$t, software.debugging$Nt,
       type="l")
n <- dim(software.debugging)[1]
software.debugging.diff <- data.frame(t=software.debugging$t[2:n], 
                                      t.diff=software.debugging$t[2:n]-software.debugging$t[1:(n-1)],
                                      Nt.diff=software.debugging$Nt[2:n]-software.debugging$Nt[1:(n-1)],
                                      Ct.diff=software.debugging$Ct[2:n]-software.debugging$Ct[1:(n-1)])
software.debugging.diff$dCt <- software.debugging.diff$Ct.diff/software.debugging.diff$t.diff
software.debugging.diff$dNt <- software.debugging.diff$Nt.diff/software.debugging.diff$t.diff
software.debugging.diff <- rbind(c(0,0,0,0,0,0),software.debugging.diff)
lambda.t <- function(theta, alpha, beta, j, data) {
  diff.1 <- (exp(-theta*data$t[j-1])-exp(-theta*data$t[j]))
  diff.2 <- (alpha+beta*sum(data$Ct.diff[1:(j-1)]*exp(theta*data$t[1:(j-1)])))
  return(diff.1*diff.2)
}
neg.log.lik <- function(x) {
  theta <- x[1]
  alpha <- x[2]
  beta  <- x[3]
  J <- dim(software.debugging.diff)[1]
  log.lik.temp <- 0
  for( j in c(2:J) ) {
    lambda.j <- lambda.t(theta,alpha,beta,j, software.debugging.diff)
    log.lik.temp <- log.lik.temp - lambda.j
    log.lik.temp <- log.lik.temp + software.debugging.diff$Nt.diff[j]*log(lambda.j)
  }
  return(-log.lik.temp)
}
fitted <- optim(c(0.001, 19, 0.005), neg.log.lik)
print(fitted)
</code>
=== まずは不具合発見数の成長を眺める ===

累積不具合発見数のグラフは次の一行で描ける。

<code>
plot(software.debugging$t, software.debugging$Nt,
     type="l",
     xlab="Staff Days", ylab="Cumulative Faults")
</code>

特にモデルは用いていないが、平滑化した曲線も付与してみる。

<code>
library(gam)
software.debugging.gam <- gam(Nt~-1+bs(t),
                            data=software.debugging)
plot(software.debugging.gam,
     xlab="Staff Days", ylab="Cumulative Faults")
points(software.debugging$t, software.debugging$Nt,
       type="l")
</code>

飽和しているように見えるのは、平滑化の特性なので、ここでは無視する。

次に、不具合発見数、改修行数とも、時点毎の差分を取ってみる。

<code>
n <- dim(software.debugging)[1]
software.debugging.diff <- data.frame(t=software.debugging$t[2:n], 
                                    t.diff=software.debugging$t[2:n]-software.debugging$t[1:(n-1)],
                                    Nt.diff=software.debugging$Nt[2:n]-software.debugging$Nt[1:(n-1)],
                                    Ct.diff=software.debugging$Ct[2:n]-software.debugging$Ct[1:(n-1)])
software.debugging.diff$dCt <- software.debugging.diff$Ct.diff/software.debugging.diff$t.diff
software.debugging.diff$dNt <- software.debugging.diff$Nt.diff/software.debugging.diff$t.diff
</code>

時点間の不具合発見数のグラフを描く。

<code>
software.debugging.gam <- gam(Nt.diff~-1+bs(t),
                            data=software.debugging.diff)
plot(software.debugging.gam,
     xlab="Staff Days", ylab="1st Order Diff of Cumulative Faults")
points(software.debugging.diff$t,software.debugging.diff$Nt.diff,
     pch=20)
</code>

一階差分を平滑化してみると、減少しているように見える。

二階差分も確認してみる。
<code>
n <- dim(software.debugging.diff)[1]
software.debugging.diff.2 <- data.frame(t=software.debugging.diff$t[2:n], 
                                    t.diff=software.debugging.diff$t[2:n]-software.debugging.diff$t[1:(n-1)],
                                    Nt.diff.2=software.debugging.diff$Nt[2:n]-software.debugging.diff$Nt[1:(n-1)],
                                    Ct.diff.2=software.debugging.diff$Ct[2:n]-software.debugging.diff$Ct[1:(n-1)])
software.debugging.diff.2$dCt <- software.debugging.diff.2$Ct.diff.2/software.debugging.diff.2$t
software.debugging.diff.2$dNt <- software.debugging.diff.2$Nt.diff.2/software.debugging.diff.2$t
</code>

グラフを描いてみる。

<code>
software.debugging.gam <- gam(Nt.diff.2~-1+bs(t),
                            data=software.debugging.diff.2)
plot(software.debugging.gam,
     xlab="Staff Days", ylab="2nd Order Diff of Cumulative Faults")
points(software.debugging.diff.2$t,software.debugging.diff.2$Nt.diff.2,
     pch=20)
</code>

平滑化しないと、こんな感じ。

<code>
plot(software.debugging.diff.2$t,software.debugging.diff.2$Nt.diff.2,
     pch=20)
</code>

フラットかもしれない。

時間あたりに直してみる。

<code>
software.debugging.gam <- gam(dNt~-1+bs(t),
                            data=software.debugging.diff.2)
plot(software.debugging.gam,
     xlab="Staff Days", ylab="2nd Order Diff of Cumulative Faults per Days")
points(software.debugging.diff.2$t,software.debugging.diff.2$dNt,
     pch=20)
</code>

二階差分の非負値を青、負値を赤でプロットしてみると、徐々に負に向かっているようにも見える。

<code>
plot(software.debugging.diff.2$t[software.debugging.diff.2$dNt>=0],
     software.debugging.diff.2$dNt[software.debugging.diff.2$dNt>=0],
     xlim=c(min(software.debugging.diff.2$t),
            max(software.debugging.diff.2$t)),
     ylim=c(min(software.debugging.diff.2$dNt),
            max(software.debugging.diff.2$dNt)),
     pch=20,col="blue")
points(software.debugging.diff.2$t[software.debugging.diff.2$dNt<0],
       software.debugging.diff.2$dNt[software.debugging.diff.2$dNt<0],
       pch=20,col="red")
</code>

変曲点があるようにも見えるので、今回は飽和すると思って、分析を続ける。

=== 次に改修行数の成長を眺める ===

時点間の改修行数のグラフも描いてみる。

<code>
plot(software.debugging.diff$t,software.debugging.diff$Ct.diff,
     pch=20)
</code>

更に、改修行数を改修に要した期間で割ったグラフも描いてみる。

<code>
plot(software.debugging.diff$t,software.debugging.diff$dCt,
     pch=20)
</code>

いずれにせよ、初期に大規模な改修が行われているが、徐々に改修行数は減ってきているのが見て取れるはず。

<code>
software.debugging.gam <- gam(dCt~-1+bs(t), data=software.debugging.diff)
plot(software.debugging.gam,
     xlab="Staff Days", ylab="Correction Lines Per Staff Day")
points(software.debugging.diff$t, software.debugging.diff$dCt, pch=20)
</code>


=== モデルを当てはめる ===

今回は、用いるモデルがRには用意されていないので、

  * モデルの強度関数と、強度関数を含む対数尤度関数を、Rの関数として記述する。
  * 対数尤度関数を、optim()を用いて最大化する

という手順で、パラメータの推定値を得る。

  * このデータの場合、データ観測系列は1本なので、ロバストな推定方法は使わない。
  * 先週の説明の通り、飽和しかけているので、曲線を強度関数に用いる (累積データに「飽和する曲線」を当てはめる場合、変曲点後まで観測できているか否かで、当てはまりがとても異なるので注意)

データの先頭に1行、0を入れておく。

<code>
software.debugging.diff <- rbind(c(0,0,0,0,0,0),software.debugging.diff)
</code>

強度関数は，次の通り。

<code>
lambda.t <- function(theta, alpha, beta, j, data) {
  diff.1 <- (exp(-theta*data$t[j-1])-exp(-theta*data$t[j]))
  diff.2 <- (alpha+beta*sum(data$Ct.diff[1:(j-1)]*exp(theta*data$t[1:(j-1)])))
  return(diff.1*diff.2)
}
#lambda.t(0.001704,19.1745,0.003045,10,software.debugging.diff)
</code>

これを用いた対数尤度関数を、次のように定義する。
<code>
neg.log.lik <- function(x) {
  theta <- x[1]
  alpha <- x[2]
  beta  <- x[3]
  J <- dim(software.debugging.diff)[1]
  log.lik.temp <- 0
  for( j in c(2:J) ) {
    lambda.j <- lambda.t(theta,alpha,beta,j, software.debugging.diff)
    log.lik.temp <- log.lik.temp - lambda.j
    log.lik.temp <- log.lik.temp + software.debugging.diff$Nt.diff[j]*log(lambda.j)
  }
  return(-log.lik.temp)
}
</code>

対数尤度の最大化は、上の関数の最小化と等しい。

Rには最小化の関数が幾つかある。nlminb(), optim(), nlm()など。

これを最小化するのは、結構、困難であるが、例えば初期値を

<code>
fitted <- optim(c(0.001, 19, 0.005), neg.log.lik)
print(fitted)
</code>

のように与えると、それなりの値に収束する。optim()の最初の引数は、最適化の初期値である。これは Cook and Lawless (2007)掲載の推定値を参考に設定した。
以下で、モデルを大幅に変更すると、初期値の探索もしなければならなくなるので、気をつけられたし。

== 補足 ==
optimという関数は、目的関数と初期値を与えると最適化してくれる関数、で、幾つかの有名な最適化アルゴリズムを利用できて便利だが、探索範囲を制限することができない。

今回の演習で用いているのはNelder-Meed法という、Simplex法に似たアルゴリズムで、初期値から出発して、探索していく。その際に、大きな歩幅で探索すると、warningが表示されるように計算きない点に到達する。そのため、今回のような「計算可能な探索範囲」が狭い最適化問題では、目的関数を評価できない旨のwarningが沢山、ログに残る。ただし、初期値の近傍の局所最適解に到達している場合には、その解が表示されるので、warningがある旨のメッセージではなく、最終的な結果が出力されているか否かで、最適化が行えたかどうか、判断して欲しい。


=== 強度関数を変えてみる ===

ここでは、強度関数の変え方の一例を示す。

  * 強度関数を変えるにはlambda.t()だけでなく、neg.log.lik()も変える必要が出てくることもある。今回は、パラメータを追加したので、neg.log.lik()の中で x[4] の扱いを追加し、4変数関数としなければならない。
  * さらに閾値を変えてみるとして、betaとbeta.2の切り替え点を300から500に変えるなど、してみている。

過去の修正の影響がもっと軽い可能性はないかと、betaを2種類にしてみる。
まずは300日以内の影響をbetaとし、300日以降の影響をbeta.2とするように、強度関数を変更する。

<code>
lambda.t <- function(theta, alpha, beta, beta.2, j, data) {
  diff.1 <- (exp(-theta*data$t[j-1])-exp(-theta*data$t[j]))
  data.t <- data$t[1:(j-1)]
  beta.s <- data.t
  beta.s[(data$t[j]-data.t)<=300] <- beta
  beta.s[(data$t[j]-data.t)>300] <- beta.2
  diff.2 <- (alpha+sum(beta.s*data$Ct.diff[1:(j-1)]*exp(theta*data$t[1:(j-1)])))
  return(diff.1*diff.2)
}
lambda.t(0.001704,19.1745,0.003045,0.003045,10,software.debugging.diff)
</code>

強度関数のパラメータがひとつ増えたので、対数尤度関数の引数xの次元もひとつ増える。
そのため、一行追加する。

<code>
neg.log.lik <- function(x) {
  theta <- x[1]
  alpha <- x[2]
  beta  <- x[3]
  beta.2 <- x[4]
  J <- dim(software.debugging.diff)[1]
  log.lik.temp <- 0
  for( j in c(2:J) ) {
    lambda.j <- lambda.t(theta,alpha,beta,beta.2,j, software.debugging.diff)
    log.lik.temp <- log.lik.temp - lambda.j
    log.lik.temp <- log.lik.temp + software.debugging.diff$Nt.diff[j]*log(lambda.j)
  }
  return(-log.lik.temp)
}
</code>

直近の改修よりも、以前の改修の方が、不具合の発生に大きく影響を与える、という母数の推定値が得られる。
これは、初期に多くの改修を行ったことと、無関係ではないかもしれない。

以下はbetaとbeta.2の切り替え時点を300から500に変更してみる場合の強度関数である。

<code>
lambda.t <- function(theta, alpha, beta, beta.2, j, data) {
  diff.1 <- (exp(-theta*data$t[j-1])-exp(-theta*data$t[j]))
  data.t <- data$t[1:(j-1)]
  beta.s <- data.t
  beta.s[(data$t[j]-data.t)<=500] <- beta
  beta.s[(data$t[j]-data.t)>500] <- beta.2
  diff.2 <- (alpha+sum(beta.s*data$Ct.diff[1:(j-1)]*exp(theta*data$t[1:(j-1)])))
  return(diff.1*diff.2)
}
</code>

このような変更手段だけでも、かなりのバリエーションにはなるので、いろいろ試してみて欲しい。

=== モデルを選ぶ ===

AICは、モデルの当てはまりの良さを評価する指標のひとつである。

  * AICは予測精度の近似として導出されている。
  * AICを基準にモデルを選ぶ、とはこの値が最小となるように、モデルを選ぶことである。
  * AICを最小化するモデルが、必ずしも正解とは限らない。

今回は、次のように計算される量をAICとして用いる。

<code>
fitted <- optim(c(0.001,19,0.005,0.005),neg.log.lik)
print(fitted$value*2+2*length(fitted$par))
</code>

これをなるべく小さくするようなモデルを用いるのが望ましいが、対象についての知識・情報がある場合には、それに基づいて、最小ではないモデルを選択することも少なくない。

=== 残差をプロットする ===

残差とは、推定したモデルと推定に用いたデータとの差を表す量である。モデルで説明されなかった、残りの部分、とも言える。

  * 残差のばらつきが満足できる小ささであること (予測に使えそうな精度であること)
  * 残差が何らかの傾向を持たないこと (他に考慮に入れるべき共変量がないこと)

は、モデルを推定した後で確認しなければならない。

時点jと時点j-1の間の強度は
<code>
lambda.t(fitted$par[1], fitted$par[2], fitted$par[3], fitted$par[4],j,software.debugging.diff)
</code>
で求められる。テストデータに付与するには、
<code>
fitted.int <- function(theta, alpha, beta, beta.2, data) {
  lambda <- software.debugging.diff$t
  lambda <- 0
  J <- dim(software.debugging.diff)[1]
  for( j in c(2:J) ) {
    lambda[j] <- lambda.t(theta, alpha, beta, beta.2, j, data)
  }
  Lambda <- cumsum(lambda)
  return(list(time=data$t,lambda=lambda, Lambda=Lambda, Nt=cumsum(data$Nt.diff)))
}
</code>

を用いて、

<code>
fitted.diff <- fitted.int(fitted$par[1], fitted$par[2], fitted$par[3], fitted$par[4],
                   software.debugging.diff)
</code>

とする。グラフに描くには、

<code>
plot(fitted.diff$t, fitted.diff$Nt, type="l")
lines(fitted.diff$t, fitted.diff$Lambda, type="l", lty=2)
</code>

とする。残差のプロットは

<code>
plot(fitted.diff$t, fitted.diff$Nt-fitted.diff$Lambda, type="p", pch=20)
</code>

標準化残差のプロットは

<code>
plot(fitted.diff$t, (fitted.diff$Nt-fitted.diff$Lambda)/sqrt(fitted.diff$Lambda), type="p", pch=20)
</code>

で得られる。

=== 更に強度関数を変えてみる ===

<code>
lambda.t <- function(theta, alpha, beta, beta.2, j, data) {
  diff.1 <- (exp(-theta*data$t[j-1])-exp(-theta*data$t[j]))
  data.t <- data$t[1:(j-1)]
  beta.s <- data.t
  beta.s[(data$t[j])<=300] <- beta
  beta.s[(data$t[j])>300] <- beta.2
  diff.2 <- (alpha+sum(beta.s*data$Ct.diff[1:(j-1)]*exp(theta*data$t[1:(j-1)])))
  return(diff.1*diff.2)
}
</code>

今のところ、これが一番良いモデル。
他にないか、頑張ってみて。