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