==== レポート課題 (Quiz) #1 (2013.06.27, Due:2013.07.10) ====

1対象あたり高々1回しかイベントが発生しない状況での、マルチンゲール中心極限定理の観察。

=== 課題 Quizes ===

Quizes

  - Observe how martingale central limit theorems work.
  - Guess the relationships among means, variances, sample sizes (n.target), and simulation times (n.simulation).

課題

  - マルチンゲール中心極限定理が作用している様子を観察(確認)しなさい。
  - 背景となる分布についての関心は横においておき、分布ごとに、計数過程の平均、分散、対象数 n.target、シミュレーション回数 n.simulation の関係を調べなさい。

=== 準備 (Preparation) ===

Instruction
  - Download and install R from [[http://cran.r-project.org/|CRAN]]. ([[http://cran.ism.ac.jp/|The mirror site at ISM (Tachikawa)]] is  highly recommended for your choice.)
  - Run R either by clicking the icon of R or by typing "R" on your console.
  - Copy the following code.
  - Paste those lines to R.

手順
  - Rというプログラミング言語(インタプリタ)のインストーラを[[http://cran.r-project.org/|CRAN]]からダウンロードしてインストールする (この大学からは、[[http://cran.ism.ac.jp/|統計数理研究所のミラーサイト]]を利用するのが、速くてとてもお勧め)
  - Rのアイコンをクリックするか、コンソールでRと入力して、Rを起動する
  - 下の囲みの中をコードをすべてコピーする
  - 起動したRに貼り付けて実行させておく

<code>
# A function to plot the sample path of the given event history.
plot.counting.process <- function(event.history,
                                  origin=c(0,0),
                                  printing=FALSE,
                                  add=FALSE) {
  n.events <- length(event.history)
  if(printing==TRUE) print(length(event.history))

  time.init <- origin[1]
  time.history <- append(time.init,
                         time.init)
  time.history <- append(time.history,
                         sort(append(event.history, event.history)))
  if(printing==TRUE) print(length(time.history))

  count.init <- origin[2]
  count.history <- append(count.init,
                          sort(append(c(1:n.events),
                                      c(1:n.events)
                                      )
                               )-1+count.init)
  count.history <- append(count.history,
                          n.events+count.init)
  if(printing==TRUE) print(length(count.history))
  if(add==TRUE) {
    lines(time.history, count.history)
  } else {
    plot(time.history, count.history, type="l")
  }
}

# A function to plot the sample path of
# the martingale calculated from the event history
# and the compensator.
plot.counting.process.martingale <- function(event.history,
                                             compensator,
                                             origin=c(0,0),
                                             printing=FALSE,
                                             add=FALSE) {
  n.events <- length(event.history)
  if(printing==TRUE) print(length(event.history))

  time.init <- origin[1]
  time.history <- append(time.init,
                         time.init)
  time.history <- append(time.history,
                         sort(append(event.history, event.history)))
  if(printing==TRUE) print(length(time.history))

  count.init <- origin[2]
  count.history <- append(count.init,
                          sort(append(c(1:n.events),
                                      c(1:n.events)
                                      )
                               )-1+count.init)
  count.history <- append(count.history,
                          n.events+count.init)
  martingale <- count.history - compensator(time.history)
  if(printing==TRUE) print(length(count.history))
  if(add==TRUE) {
    lines(time.history, martingale, col="grey")
  } else {
    plot(time.history, martingale, type="l")
  }
}

# A function to count the sample path of
# the martingale calculated from the event history
# and the compensator.
count.counting.process <- function(event.history,
                                   time.to.count) {
  n.time <- length(time.to.count)
  if( n.time > 1 ) {
    c.events <- NULL
    for( i in c(1:n.time) ) {
      c.events <- append(c.events,
                         sum(event.history<=time.to.count[i]))
    }
  } else {
    c.events <- sum(event.history<=time.to.count)
  }
  return(c.events)
}

count.martingale <- function(event.history,
                             compensator,
                             time.to.count) {
  n.time <- length(time.to.count)
  if( n.time > 1 ) {
    c.events <- NULL
    for( i in c(1:n.time) ) {
      c.events <- append(c.events,
                         sum(event.history<=time.to.count[i]))
    }
  } else {
    c.events <- sum(event.history<=time.to.count)
  }
  m.events <- c.events - compensator(time.to.count)
  return(m.events)
}
</code>

以上のコードは次の4つの関数を含む。

  - plot.counting.process: 計数過程のサンプルパスの描画
  - plot.counting.process.martingale: 計数過程から作ったマルチンゲールのサンプルパス (Doob-Meyer分解でN=A+Mとあるところを、M=N-AとしてMのプロット)
  - count.counting.process: 特定の時点までの累積発生回数の数え上げ
  - count.counting.process.martingale: 未完

=== シミュレーションデータの生成コード ===

シミュレーションデータを生成する前に、対象数をまず設定する。

<code>
# Generate "an event history"
# The number of targets, each has an event
# during its lifetime.
n.target <- 10
n.target <- 100
n.target <- 1000
</code>

今回の課題では「事象発生を支配する確率法則」を3つ提供する。

  - ワイブル分布
  - 指数分布
  - ガンマ分布

パラメータは法則ごとに異なる変数で保持しているが、データの生成はgenerator.x()、compensatorの関数はcompensator.x()と、名前を共通にしてある。
これはRの関数渡し(関数の引数に関数を持たせることができる)の機能を使っているため。
シミュレーションの実施の際には、どの法則に基づいたシミュレーションを行っているか、意識して確認しながら進める必要がある。

<code>
# I offer three underlying models for the counting process.
# We restrict ourselves to the cases where each target has
# an event once at most. 

# Weibull distribution:
#       with shape parameter 2 and scale parameter 1/gamma(1+1/2)
#       so as to we have E[X] = 1
wb.shape <- 2
wb.scale <- 1/gamma(1+1/2)
generator.x <- function(n) {
  return(rweibull(n=n.target, shape=wb.shape, scale=wb.scale))
}
compensator.x <- function(x) {
  return(n.target * pweibull(x, shape=wb.shape, scale=wb.scale))
}

# Exponential distribution:
#       with rate 1
#       so as to we have E[X] = 1
exp.rate <- 1
generator.x <- function(n) {
  return(rexp(n=n.target, rate=exp.rate))
}
compensator.x <- function(x) {
  return(n.target * pexp(x, rate=exp.rate))
}

# Gamma distribution:
#       with shape parameter 2 and scale parameter 1/2
#       so as to we have E[X] = 1
gamma.shape <- 2
gamma.scale <- 1/2
generator.x <- function(n) {
  return(rgamma(n=n.target, shape=gamma.shape, scale=gamma.scale))
}
compensator.x <- function(x) {
  return(n.target * pgamma(x, shape=gamma.shape, scale=gamma.scale))
}
</code>

=== シミュレーションの実施 ===

<code>
# Simulations
# The number of repetitions
n.simulation <- 10
n.simulation <- 100
n.simulation <- 1000

# Generate event histories with "generator.x" function
event.history.data <- NULL
for( i in c(1:n.simulation) ) {
  event.history.data <- rbind(event.history.data,
                              generator.x(n.target))
}

# You can check the sample mean and the sample standard deviation
# by running the following two lines.
mean(event.history.data)
sqrt(var(c(event.history.data)))

# Plot event histories and the compensator
plot(c(0, max(event.history.data)), c(0, n.target),
     xlab="event time",
     ylab="event counts", 
     type="n")
for( i in c(1:dim(event.history.data)[1]) ) {
  plot.counting.process(event.history.data[i,],
                        printing=FALSE, add=TRUE)
}
plot(compensator.x, xlim=c(0, n.target), add=TRUE)
</code>

マルチンゲールをプロットしてみる。ついでに特定の時点でのイベントの累積発生回数も集計させる。
たとえば4時点 0.5, 1.0, 1.5, 2 での集計を行うには、
<code>
time.list <- c(0.5, 1, 1.5, 2)
</code>
としてから、
<code>
# Plot martingales
plot(c(0, max(event.history.data)), c(-10, 10),
     xlab="event time",
     ylab="martingales", 
     type="n")
cp.clt <- NULL
for( i in c(1:dim(event.history.data)[1]) ) {
  plot.counting.process.martingale(event.history.data[i,],
                                   compensator.x,
                                   printing=FALSE,
                                   add=TRUE)
  cp.clt <- rbind(cp.clt,
                  count.counting.process(event.history.data[i,],
                                         time.list))
}
lines(c(0, max(event.history.data)), c(0,0))
colnames(cp.clt) <- time.list
</code>
を実行する。これで表示されるグラフは縦軸が(-10, 10)の範囲で固定されているので、「c(-10, 10)」という箇所を変更して調整してほしい。
これを例えば「c(-5,5)」とすると、(-5, 5)の範囲の描画になる。

<code>
par(mfrow=c(ceiling(sqrt(length(time.list))),
            ceiling(sqrt(length(time.list)))) )
par(cex=0.5)
for( i in c(1: length(time.list)) ) {
  hist(cp.clt[,i])
}
par(mfrow=c(1,1), cex=1.0)
</code>

<code>
pairs(cp.clt)
</code>

その他、

<code>
summary(cp.clt)
apply(cp.clt,2,mean)
sqrt(apply(cp.clt,2,var))
</code>

なども参考になるかもしれない。

=== データの出力 ===

データの処理に、RではなくExcelもしくは他の言語が便利だと思う人は、下記のようにCSVファイルに出力して、そちらで集計や作図を行ってもよい。 (If you are familier with spread sheet software like Excel or any other data analysis environments, you can export the simulation data as CSV files.)

<code>
# If you want to export
write.csv(event.history.data, "event-history-data.csv")
write.csv(cp.clt, "counts-at-1-2-3.csv")
</code>

=== レポート提出 Report Submission ===

電子メールにて、[[mailto:watalu_at_inf.uec.ac.jp|山本]]まで送付すること。メールアドレスには適切な修正が必要だが、たぶん推測できるはず。

(Please send your report file to [[mailto:watalu_at_at_inf|me]]. You need to correct my email address.  But I believe you can guess it.)

==== レポート課題 (Quiz) #2 (2013.07.25, Due:2013.08.08) ====

再発事象データの分析を２ケース分。

=== 課題 Quizes ===

Quizes

  - Analyse a data set of recurrent events using estimating function approach. (coxph)
  - Analyze a data set of software tests by fitting intensity functions. 

課題

  - 比例ハザードモデルを用いて、Bladder Cancer Recurrences Dataを解析してみなさい。
  - Software Debugging Dataに、過去の履歴に依存するような強度関数を当てはめて解析してみなさい。

=== Bladder Cancer Recurrences Data ===

ここでは 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で得たモデルと、ベースラインハザード関数、回帰係数、残差プロットなど、比較してみよ。

=== まとめ ===

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

=== Software Debugging Data ===

  - 解析対象は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>

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

=== レポート提出 Report Submission ===

電子メールにて、[[mailto:watalu_at_inf.uec.ac.jp|山本]]まで送付すること。メールアドレスには適切な修正が必要だが、たぶん推測できるはず。

(Please send your report file to [[mailto:watalu_at_at_inf.uec.ac.jp|me]]. You need to correct my email address.  But I believe you can guess it.)