Как получить доверительные интервалы, не инвертируя особую матрицу Гессе в R? - PullRequest
Новая
       16

Как получить доверительные интервалы, не инвертируя особую матрицу Гессе в R?

0 голосов
/ 21 апреля 2010

Я студент, работающий над моделью эпидемиологии в 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, например, уже объявлено. Я сделал это просто для того, чтобы мои последующие вызовы функции не обязательно содержали ужасные векторы.

Затем я оптимизировал CV1, CV2, v1 и v2, используя mle2 (библиотека bbmle). Это также выглядит немного грубо и выглядит так: 

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, выглядят как графики оригинальной работы ). Я также прочитал соответствующие части руководства, а также книгу Болкера Экологические модели в R . Я также попробовал различные методы оптимизации, которые привели к увеличению времени выполнения, но с теми же ошибками. Метод «SANN» не завершил работу в течение часа, поэтому я не стал ждать результата.

В двух словах: мои доверительные интервалы плохие. Есть ли относительно простой способ исправить их в 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

Также возможно, что я делаю что-то еще совершенно неправильно, но мои результаты кажутся разумными, поэтому я не уловил это.

1 Ответ

6 голосов
/ 27 апреля 2010

Вы можете изменить параметризацию, чтобы ограничения всегда выполнялись.Перепишите вероятность как функцию от 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)))
 }
...