0

私は最尤法を使用して、R で疫学モデルに取り組んでいる学生です。負の対数尤度関数を作成しました。ちょっと大げさですが、こんな感じです。

NLLdiff = function(v1, CV1, v2, CV2, st1 = (czI01 - czV01), st2 = (czI02 - czV02), st01 = czI01, st02 = czI02, tt1 = czT01, tt2 = czT02) { 
    prob1 = (1 + v1 * CV1 * tt1)^(-1/CV1)
    prob2 = ( 1 + v2 * CV2 * tt2)^(-1/CV2) 
    -(sum(dbinom(st1, st01, prob1, log = T)) + sum(dbinom(st2, st02, prob2, log = T)))
 }

最初の行がひどく見える理由は、そこに入力されるデータのほとんどがそこにあるからです。czI01、たとえば、すでに宣言されています。私がこれを行ったのは、後で関数を呼び出すときに、ひどいベクトルが含まれている必要がないようにするためです。

次に、mle2 (ライブラリ bbmle) を使用して、CV1、CV2、v1、および v2 用に最適化しました。これも少し見栄えが悪く、次のようになります。

ml.cz.diff = mle2 (NLLdiff, start=list(v1 = vguess, CV1 = cguess, v2 = vguess, CV2 = cguess), method="L-BFGS-B", lower = 0.0001)

さて、ここまでは問題なく動作しています。ml.cz.diff は、データに適切に適合するプロットに変換できる値を提供します。また、いくつかの異なるモデルがあり、それらを比較するために AICc 値を取得できます。ただし、v1、CV1、v2、および CV2 の周りの信頼区間を取得しようとすると、問題が発生します。基本的に、私は CV1 に負の境界を取得しますが、これは実際には生物学モデルの平方数といくつかの警告を表すため、不可能です。

信頼区間を取得するより良い方法はありますか? それとも、本当に、ここで意味のある信頼区間を取得する方法はありますか?

私が見ているのは、偶然にも、私のヘッセ行列が最適化空間のいくつかの値に対して特異であるということです。しかし、私は 4 つの変数を最適化していて、プログラミングの知識があまりないため、ヘシアンに依存しない適切な最適化方法を思いつくことができません。私は問題をグーグルで検索しました-それは私のモデルが悪いことを示唆していましたが、私のモデルが本当にひどいものではないことを示唆する以前に行われたいくつかの作業を再構築しています(ml.cz.diffを使用して作成したプロットは、元の作業のプロットのように見えます) ）。マニュアルの関連部分とボルカーの本Ecological Models in Rも読みました. さまざまな最適化方法も試しましたが、実行時間は長くなりましたが、同じエラーが発生しました。「SANN」メソッドは 1 時間以内に実行が終了しなかったため、結果が表示されるのを待ちませんでした。

一言で言えば、私の信頼区間は悪いです。Rでそれらを修正する比較的簡単な方法はありますか?

私のベクトルは次のとおりです。

czT01 = c(5, 5, 5, 5, 5, 5, 5, 25, 25, 25, 25, 25, 25, 25, 50, 50, 50, 50, 50, 50, 50)
czT02 = c(5, 5, 5, 5, 5, 10, 10, 10, 10, 10, 25, 25, 25, 25, 25, 50, 50, 50, 50, 50, 75, 75, 75, 75, 75)
czI01 = c(25, 24, 22, 22, 26, 23, 25, 25, 25, 23, 25, 18, 21, 24, 22, 23, 25, 23, 25, 25, 25)
czI02 = c(13, 16, 5, 18, 16, 13, 17, 22, 13, 15, 15, 22, 12, 12, 13, 13, 11, 19, 21, 13, 21, 18, 16, 15, 11)
czV01 = c(1, 4, 5, 5, 2, 3, 4, 11, 8, 1, 11, 12, 10, 16, 5, 15, 18, 12, 23, 13, 22)
czV02 = c(0, 3, 1, 5, 1, 6, 3, 4, 7, 12, 2, 8, 8, 5, 3, 6, 4, 6, 11, 5, 11, 1, 13, 9, 7)

そして、私は次のように推測します：

v = -log((c(czI01, czI02) - c(czV01, czV02))/c(czI01, czI02))/c(czT01, czT02)
vguess = mean(v)
cguess = var(v)/vguess^2

私が完全に間違ったことをしている可能性もありますが、私の結果は合理的であるように見えるので、私はそれを捕まえていません.

4

1 に答える 1

6

制約が常に満たされるようにパラメーター化を変更できます。尤度を ln(CV1) と ln(CV2) の関数として書き直すと、CV1 と CV2 が厳密に正のままであることを確認できます。

NLLdiff_2 = function(v1, lnCV1, v2, lnCV2, st1 = (czI01 - czV01), st2 = (czI02 - czV02), st01 = czI01, st02 = czI02, tt1 = czT01, tt2 = czT02) { 
prob1 = (1 + v1 * exp(lnCV1) * tt1)^(-1/exp(lnCV1))
prob2 = ( 1 + v2 * exp(lnCV2) * tt2)^(-1/exp(lnCV2)) 
-(sum(dbinom(st1, st01, prob1, log = T)) + sum(dbinom(st2, st02, prob2, log = T)))
 }
于 2010-04-27T08:38:10.887 に答える