生存時間解析における「打ち切り」と「切断」

　生存時間解析には、打ち切りと切断という、少々紛らわしい用語が存在する。院生ゼミの時にちょっと関連知識を紹介した。その時の様子をツィートしたら、院生のところに「どう違うのか」と別の大学の院生から質問が来たらしい（すぐに質問するフットワークの軽さは見習うべきだけど。。。）。いい機会なので、理解している範囲でまとめてみる。

　結論から言うと、打ち切りありも切断ありも、観測データが不完全であることを指す点では同じだ。しかし、（少なくとも生存時間解析の文脈では）切断ありが観測対象の選び方の不完全さに関連するのに対し、打ち切りは観測自体の不完全さ（例えば、観測期間の短さ）に由来する。以下に例を挙げながら説明する。

1. 打ち切りについて

　打ち切りは、生態学関連でも出てくるし、RESTでも使うので私にとってが馴染み深い。

　（なんでもいいのだけど、あえてRESTっぽい状況を想定して）動物のある行動の持続時間を測りたいとしよう。行動開始から終了までを常に観察できれば理想的だが、必ずしもそうはいかない。例えば、動物を観察・記録するための装置が「ある時間が経過したら勝手に撮影を終了してしまう」仕様になっているかもしれないし、行動開始時点ではうっかりしていて気づかずに「途中からしか観察出来なかった』ということもあるかもしれない。

　これらのデータをそのまま使って行動持続時間を推定すると、過小評価してしまう。このような制約を伴ったデータのことを打ち切りデータと呼ぶ。しかし、打ち切りデータも「少なくとも〇〇より長い」ということを教えてくれるので、貴重なデータなので捨ててしまうのはもったいない。

2. 打ち切りの種類

　打ち切りには、左打ち切りと右打ち切り、両側打ち切り、そして間隔打ち切りがある。

　左打ち切りというのは、先程あげた例の2番目が該当する。すなわち、開始時点が不明な状況だ。右打ち切りとは、先ほどの1番目の例のように、終了時点が不明な状況を指す。そして、両側打ち切りとは、この両方が生じている場合を指す。

　最後の間隔打ち切りは、前提とする観測条件が少しだけ違う。例えば、動物がある行動を続けているかどうかを1分間隔で記録したとする。この場合、持続時間の最小値と最大値は分かるが、ピッタリこれだけの長さというデータを得ることは原理的にできない。このようなデータを間隔打ち切りという。野外データのほとんどは間隔打ち切りな気もするが、意外にそれとして扱われることは少ない気もする。

　確率分布を仮定した解析を行う場合、打ち切りデータの尤度は割と簡単に定式化できる。それについては、こちらでも説明しているので参考にしてほしい。リンク先でも述べている通り、理屈は簡単だが、実際の解析でRの関数を使おうとすると、打ち切り（観測なし）を1として指定するものや、観測あり（打ち切りなし）を1と表現するものがあったりして、結構面倒くさかったりする。こういうので躓くのは避けたいところだ。

3. 切断とは？

　では、切断とはどのような状況を指すのだろうか。こちらは、生態学ではあんまり出てこないが、医学の世界では頻出なようだ。

　ちょうどゼミでも使った「厨二病の発症年齢を推定する」という状況を考えてみよう。中学生を入学時点から追跡調査し、いつ発症したかを記録したとする。まあ、厨二病だと発症のタイミングが不明確だが、それが仮に分かるとしよう。また、後で図にしやすくするために、13歳の誕生日で中学に入学することにしておく。

　発症年齢を推定する上で、この調査方法には大きな制約があることは明らかだ。調査対象を中学生の厨二病患者に絞っているので、①小学生の時に発症した人と、②高校生になってから発症する人が対象から抜け落ちているからだ。

　図示すると、以下のようになる。①に該当するのデータを左切断データと呼び、②に該当するのが右切断データだ。切断という表現は、真の分布のが13歳と16歳の時点で「切れている」ことを表現したものだ。

　最初に「切断ありが観測対象の選び方の不完全さに関連するのに対し、打ち切りは観測自体の不完全さ（例えば、観測期間の短さ）に由来する」と言ったが、その意味が分かってもらえるはずだ。

4. より複雑な状況

　実際には、もっと複雑な状況が起こりうる。先ほどの厨二病の発症と同時に、治癒についても記録していたとする。そして、発症年齢、治癒年齢、治癒に要する期間を推定したいとする。この場合、様々な状況を考慮しないといけなくなる．

　全部で上の図の6通りが存在する．何を推定する場合に，どれが「打ち切り」や「切断」に相当するのかを考えてみると頭を整理できるはずだ．医学の世界では，多分，「発症のタイミングによって治癒までの期間がどうなるか」みたいな解析をする必要があるんだろうと思う．しかも，途中で治療法を変えたりしながら…これは結構大変な状況だけど，生態学では，幸いそこまで複雑な状況はまれだろう（多分）．

検出エリア内の滞在時間から移動速度を推定する

はじめに

　ちょっとしんどいことがあって，仕事にもあまり手がつかない感じなのだが，研究モチベーションの維持のためにも最近取り組んでいたことをまとめてみた．もう少し落ち着いたら後で振り返ろう．

　昔から「RESTで使う正三角形内での滞在時間から移動速度を出せないのか？」と聞かれることが多い．「やればできるだろうなあ」と思いながら，類似の問題（目的と言うより手法が類似）を，元筑波大のAさんが取り組んでいたはずなので手を出さないでいた（あの研究はどうなったんだろう？）．考えてみると，まあ，我々が使っている正三角形バージョンに限定してしまえば，とくに重なるところもないはず．ちょっとやってみることにした．

問いの確認

　取り組みたい問題は，以下の通り．まず，自動撮影カメラの前に設定した正三角形の検出エリアがある．調査の後に，動画を見ながら，動物がその検出エリア内にとどまった時間（あるいは，通過するのに要する時間．以下，滞在時間）を測定したとする（REST法を使えば，滞在時間と撮影頻度から動物の密度が推定できる）．この滞在時間には，動物の移動速度の情報が含まれていることは自明だろう．もし動物が，1. ランダムな場所から侵入しランダムな場所から脱出する，2. 検出エリア内では真っすぐ動くならば，滞在時間から動物の移動速度の推定も可能な気がする．では，具体的にはどうしたらいいのだろうか？

推定の手順

滞在時間の分解

　最初に考えなければならないのは，測定される滞在時間と速度の関係である．これは，難しい話でもなんでもなくて，\(滞在時間t = 距離y / 移動速度x\)になる．

　さて，距離yや移動速度xは，動物が通過するたびに確率的に変わる値，すなわち，確率変数である．今，移動速度xが従う確率分布をf(x)，距離yが従う確率分布をg(y)としよう．すると，滞在時間tは，2つの確率変数の商が生成する確率分布d(t)に従う確率変数と見なせる．一般に，確率変数yとxの「商」は，以下のような確率分布に従う（導出は省略．途中，ガンマ関数とか出てきてややこしいが，以下のような形に落ち着くはず¹）．

\[ d(y/x) = \int y f(x) g(y) dy \]

すなわち，f(x)とg(y)が決まれば，yを積分消去することで，滞在時間が従う確率分布をxの関数として定義できることになる．そこまでいけば，滞在時間データから速度が従う確率分布のパラメータを決めることも可能だろう（今回は，最尤法で推定する）．

移動速度の確率密度関数f(x)

　では，移動速度は，どのような確率分布に従うと想定できるだろうか？これについては，イギリスのRowcliffe大先生がパナマの熱帯雨林で詳細に調べた研究から，対数正規分布に従うことが多いということになっている．おそらく，動物は，時々立ち止まるので，裾野長い分布になるということだろう．対数正規分布は，その逆数も対数正規分布に従う（v ~ dlognorm(mealog, sdlog)の場合，1/v ~ dlognorm(-meanlong, sd)になる）という便利な性質があるので，これはうれしい（今回は使わない）．

移動速度の確率密度関数g(x)

　一方，距離の方も（既存の確率分布にはもちろん従わないが）「ある距離になる確率密度」を与える関数を定義することは不可能ではなさそう．最初，これが結構面倒くさそうに感じていたが，超優秀な東大の院生M君と話していたら，割とエレガントな導出方がありそうだったので，やってみることにする．

　正三角形を考える．動物が底辺と左斜辺の任意の点を通過したとする（下図の赤線がその軌跡）．

この時，動物の軌跡bの長さは，正弦定理から以下のようになる．

\[b^2 = a^2 + b^2 - 2 \times a \times c \times cos(pi/3)\]

このbの長さが，どのような分布になるのかをシミュレーションで確かめてみる．辺の長さは，重要ではないので，以下では１辺の長さ1の正三角形の場合を考えてみる．

n <- 100000
z <- rep(0, n)

for(i in 1:n){
  x <- runif(1, 0, 1) # 辺の任意の点
  y <- runif(1, 0, 1)
  z[i] <- sqrt(x^2 + y^2 - x*y) # bの長さ
}

h <- hist(z, freq = F, breaks = 100)
lines(density(z), col = "blue")

　結構面白い形になる．0.9のちょっと前くらいまでは，直線的に増加していき，そのあと，滑らかに減衰していく．

　このようになる理由は，次のように考えるとわかりやすい．長さの決まった棒があり，その棒の両端を正三角形の辺の上に重ねるという遊びを考える．棒の長さが短いと以下の図のように底辺にぎりぎり沿う線から垂線まで動かせる．棒は，左図の黒丸Dから黒丸Bまで動かせると考えるとわかりやすい．棒の置きやすさは，DBの長さに比例するので，棒が長いほど，うまい具合に重ねやすくなる．

　棒がさらに長くなると（三角形の高さ = sqrt(3)/2より長くなると），垂直の棒が置けなくなってしまう（右図の垂線）．この場合，辺ABとの接点である黒丸が動ける範囲は，AからFの範囲になる．

　以上のように考えると，相対的な棒の置きやすさ（=確率密度に比例する値）は，棒の長さをyとすると，

• 0 < b < sqrt(3)/2の範囲では，左図のDBの長さ．すなわち，\(g(y) = 2y/\sqrt{3}\times 定数\)

• sqrt(3)/2 <= b < 1の範囲では，右図のAF

ということになる．

　正弦定理を使えば，AFの長さも容易に求まる．最終的にAFは以下のような形になる（綺麗な式になったらうれしいなと思ったけど，そうでもなかった）．

\[ g(y) = (1 - \sqrt{4y^2 - 3})\times 定数 \]

なお，定数と書いているのは，式で算出されるのは，相対的な生じやすさで，確率密度に比例する値だから．まあ，推定上は無視していい値だが，以下では一応値を出している．

念のためRで確認

　Rで確認する．確率密度にするためには，面積を１にする必要があるので，最初に積分して面積を出して調整する．

# 規格化するための定数の算出（総面積を1にするための準備）
# yが0 ~ sqrt(3)/2の範囲
scale.s <- function(x){
  return(2*sqrt(3)/3*x)
}
integrate(scale.s, 0, sqrt(3)/2) # 念のため確認 sqrt(3)/4に一致

## 0.4330127 with absolute error < 4.8e-15

ss <- sqrt(3)/4

# yがsqrt(3)/2 ~ 1の範囲
scale.l <- function(x){
  return(1 - sqrt(4*x^2 - 3))
}
sl <- integrate(scale.l, sqrt(3)/2, 1)$value # これも頑張れば手計算可能

scale.value <- 1/(ss + sl) # これを掛ければ面積1に


# 距離分布の作成
p.dis.s <- function(x){ # y < sqrt(3)/2
  return(2*sqrt(3)/3*x*scale.value)
}
p.dis.l <- function(x){ # sqrt(3)/2 < y < 1
  #　4*((sqrt(3)/2)^2) - 3 = 0にならないのでroundしておく
  return((1 - sqrt(round(4*(x^2) - 3, 10)))*scale.value) 
}

# 一応，シミュレーション結果と一致することを確認
n <- 100000
z <- rep(0, n)

for(i in 1:n){
  x <- runif(1, 0, 1)
  y <- runif(1, 0, 1)
  z[i] <- sqrt(x^2 + y^2 - x*y)
}

h <- hist(z, freq = F, breaks = 100)
x0 <- seq(0, sqrt(3)/2, by = 0.001)
y0 <- p.dis.s(x0)
lines(x0, y0, col = "red", lwd = 2)

x1 <- seq(sqrt(3)/2, 1, by = 0.001)
y1 <- p.dis.l(x1)
lines(x1, y1, col = "red", lwd = 2)

赤の線が式によって得られたモノ．シミュレーションとぴったりと一致したので，大丈夫そう．

尤度関数の定義

　では，相対的な生じやすさが分かったので，最初の式に戻って，滞在時間の確率密度関数を定義することにする．積分が出てくるが，多分陽には解けないので，Rで近似計算する．関数の中でfor関数を回すのは結構気持ち悪いけど，この辺りは慣れの問題．あと，最尤推定の際に必要なので，対数尤度にマイナスをつけた形で定義している．

# マイナス対数尤度関数の作成
minus.loglf <- function(stay){
  lik <- numeric(0)
  return(function(para){
    for(i in 1:length(stay)){
      int.fs <- function(d){
        d * dlnorm(d/stay[i], para[1], para[2]) * p.dis.s(d)
      }
      int.fl <- function(d){
        d * dlnorm(d/stay[i], para[1], para[2]) * p.dis.l(d)
      }
      bon.d <- sqrt(3)/2
      lik[i] <- integrate(int.fs, lower = 0, upper = bon.d)$value + integrate(int.fl, lower = bon.d, upper = 1)$value
    }
    loglik <- -sum(log(lik))
    return(loglik)})
}

最尤推定

　optim関数を使って，最尤推定してみる．

# データ準備
nsim <- 300 # 多めに300個の滞在時間を作成
meanlog <- 1 # これが正解
sdlog <- 0.2 # これが正解
v <- rlnorm(nsim, meanlog = meanlog, sdlog = sdlog) # 速度が対数正規分布に従う

# 移動距離の作成
dist <- rep(0, nsim)
for(i in 1:nsim){
  x <- runif(1, 0, 1)
  y <- runif(1, 0, 1)
  dist[i] <- sqrt(x^2 + y^2 - x*y) # 距離の乱数
}
hist(dist)

stay <- dist/v # 速度と距離から滞在時間を作成
hist(stay,  breaks = 10)

# 距離分布の作成
p.dis.s <- function(x){ # d < sqrt(3)/2
  return(2*sqrt(3)/3*x)
}
p.dis.l <- function(x){ # sqrt(3)/2 < d < 1
   #　4*((sqrt(3)/2)^2) - 3 = 0にならないのでroundしておく
  return((1 - sqrt(round(4*(x^2) - 3, 10))))
}
res <- optim(par = c(0.8, 0.15), fn = minus.loglf(stay)) # parは初期値
res

## $par
## [1] 1.032135 0.201001
## 
## $value
## [1] 899.061
## 
## $counts
## function gradient 
##       55       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

hist(v, freq = F, ylim = c(0,5), breaks = 30)
curve(dlnorm(x, meanlog, sdlog), add=TRUE, col = 'blue') # 正解
curve(dlnorm(x, res$par[1], res$par[2]), add = TRUE, col = 'green') # 推定線

　青が真の確率分布，青が推定されたパラメータによるもの．無事に一致した．

　ということで，無事，滞在時間から移動速度のパラメータを決めることが出来た．実際の動物は，直線的に移動するとは限らないので，実用的な価値は分からないが，とりあえずできてうれしい．

---

1. 「はず」というのは，色んなサイト見ながら最後は自分でやったから．割と間違った（少なくとも，推定が上手くいかない）式が書かれていることもあるので要注意．↩︎