大学院講義。

このファイル(昨年の課題)に収められたコードがすべて動くような環境を整えて下さい。

まずはRをダウンロードしてインストールします。CRANは、Rの開発組織が運営しているウェブサイトで、ここからのダウンロードは、とりあえず信用しても良いです。

• WindowsかMac OX Sなら、まずCRANにアクセスして下さい。
  1. Windowsなら、“Windows”の文字をクリックし、次に現れるページで“base”をクリックすると、“Download R X.Y.Z for Windows”という文字が見つかるので、それをクリックしてダウンロードし、実行する。
  2. Mac OS Xなら、“Mac OS X”の文字をクリックし、次に現れるページで “R-X.Y.Z.pkg (latest version)“という文字が見つかるので、それをクリックしてダウンロードし、実行する。
• Linuxなら、次の2行のどちらかが、インストールしてくれます。sudoが必要かも。

  yum install R R-*
  apt-get install r-base

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

• survivalパッケージをインストールするには、次の一行を実行します。

  install.packages(pkgs=c("survival"), repos='http://cran.md.tsukuba.ac.jp/')

この後、このファイルにあるコードをすべて実行してみてください。問題無く実行でき、赤い文字などでエラー(Error)やワーニング(Warning)が表示されなければ、成功です。

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

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

プライベートIPアドレスを持ち、NATルータを介していないコンピュータからは、インターネットに直接アクセスすることはできず、PROXYを通す必要があります。 Windowsの場合は、Rを起動する際に

C:\Program Files\R\R-X.Y.Z\bin\Rgui.exe  (Xはバージョン、Yはメジャーリビジョン、Zはマイナーリビジョン)

となっているスタートメニュー内のショートカットの代わりに

"C:\Program Files\R\R-X.Y.Z\bin\Rgui.exe" --internet2

というショートカットを新たに作ると、

〔コントロールパネル〕→〔ネットワークとインターネット接続〕→〔インターネットオプション〕→〔接続〕→〔LANの設定〕→〔プロキシサーバー〕

で設定したプロキシを使用できます。

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

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

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)))

library(survival)

2011.06.16に、今日より少し詳しい説明をしますので、コードを実行した出力だけは手元に持っていてもらえると助かります。

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

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

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

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

最後に使うので、MASSパッケージも読み込んでおく。

library(MASS)

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

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)

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

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

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)

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

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

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)

出力は次のようになる。

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

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

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   (有意になると、モデルを統計的に支持する)

以下のモデルから出発し、t検定のp値を眺めて、5%以上の変数を外したりしながら、各種検定統計量が有意になるようなモデルを探して下さい。

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)

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

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

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)

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

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

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

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

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

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

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

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

<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時間、の二種類の与え方は例示している。

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

データの読み込み

software.debugging <- read.table("http://stat.inf.uec.ac.jp/library/dm.2011/software.txt",
                           sep=",", header=T)

このデータは3つのフィールドを持つ。

今回は二つのパッケージを用いるので、インストールしておく。

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

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

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.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)

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

二階差分も確認してみる。

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

グラフを描いてみる。

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)

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

plot(software.debugging.diff.2$t,software.debugging.diff.2$Nt.diff.2,
     pch=20)

フラットかもしれない。

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

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)

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

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")

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

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

plot(software.debugging.diff$t,software.debugging.diff$Ct.diff,
     pch=20)

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

plot(software.debugging.diff$t,software.debugging.diff$dCt,
     pch=20)

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

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)

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

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

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)
}
#lambda.t(0.001704,19.1745,0.003045,10,software.debugging.diff)

これを用いた対数尤度関数を、次のように定義する。

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)
}

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

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

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

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

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

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

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

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)
}
#lambda.t(0.001704,19.1745,0.003045,10,software.debugging.diff)

これを用いた対数尤度関数を、次のように定義する。

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)
}

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

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

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

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

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

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

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

(質問を受けたので補足)

過去の修正の影響がもっと軽い可能性はないかと、betaを2種類にしてみる。

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)

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)
}

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

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)
}

AICは、モデルの予測精度を評価する指標として導出されており、次のように計算できる。

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

時点jと時点j-1の間の強度は

lambda.t(fitted$par[1], fitted$par[2], fitted$par[3], fitted$par[4],j,software.debugging.diff)

で求められる。テストデータに付与するには、

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)))
}

を用いて、

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

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

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

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

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

標準化残差のプロットは

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

で得られる。

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)
}

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

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)