Неизвестные сообщения об ошибках в R с использованием nls () с SSlogis () и предсказать ()

Я подгоняю логистическую модель (самозапуск; SSlogis) к данным о нескольких популяциях птиц с помощью nls (). Моя цель - подогнать ожидаемую функцию к данным (используя только часть каждого набора данных) и отобразить меру отклонения от ожидаемого значения на графике. Затем я хочу подобрать и построить график наблюдаемой функции (используя весь набор данных для каждой популяции), чтобы определить, укладывалась ли наблюдаемая динамика в отклонения от ожидаемого. Вот мой текущий код, написанный для этого:

    CE.mod = nls(CE.observed ~ SSlogis(t.CattleEgret, Asym, xmid, scal))

    with(collapse.data, plot(CE.time, CE.obs))

    CE.extrap = predict(CE.mod, data.frame(t.CattleEgret = CE.time))
    lines(CE.time, CE.extrap)

    CE.se.fit = sqrt(apply(attr(CE.extrap, "gradient"), 1, function(x) 
    sum(vcov(CE.mod)*outer(x,x))))

    matplot(CE.time, CE.extrap+outer(CE.se.fit, qnorm(c(0.5, 0.025, 0.975))),
    type = "l", lty = c(1,1,1), ylab = "Abundance (# per party hour)",
    xlab = "Time (year)", main = "Cattle Egret Collapse Analysis", 
    pch = 15, font.lab = 2, font.axis = 2, cex = 4, cex.lab = 1.5, 
    cex.axis = 2, cex.main = 2, frame.plot = FALSE, lwd = 4, 10)

    with(collapse.data, matpoints(CE.time, CE.obs, pch = 15, cex = 3))
    lines(CE.time, predict(nls(CE.obs ~ SSlogis(log(CE.time), 
    Asym, xmid, scal))), lty = 3, lwd = 4)

Где (из файла collapse.data):

    t.CattleEgret = c(1:20)
    CE.time = c(1:45)
    CE.obs = c(0.3061324, 0.0000100, 0.2361211, 0.5058240, 2.0685032, 2.1944544, 
               4.2689494, 4.9508297, 3.1334720, 3.6570752, 5.6753381, 10.9133183,
               5.4518257, 20.4166979, 15.9741054, 19.0970426, 13.7559959, 14.1358153, 
               15.9986416, 29.6762828, 10.3760667, 8.4284488, 6.1060359, 3.7099982, 
               3.3584060, 2.5981386, 2.5697082, 2.8091952, 5.5487979, 1.6505442,
               2.2696972, 2.1835692, 3.6747876, 4.8307886, 3.5019731, 2.8397137,
               1.8605288, 11.1848738, 2.6268683, 4.1215127, 2.3996210, 2.6569938, 
               2.1987387, 3.0267252, 2.4420927)
    CE.observed = c(0.3061324, 0.0000100, 0.2361211, 0.5058240, 2.0685032, 2.1944544, 
               4.2689494, 4.9508297, 3.1334720, 3.6570752, 5.6753381, 10.9133183,
               5.4518257, 20.4166979, 15.9741054, 19.0970426, 13.7559959, 14.1358153, 
               15.9986416, 29.6762828)

Этот код работает нормально и дает такую ​​цифру:

Анализ гибели цапли крупного рогатого скота

Если, однако, я удалю «log ()» из последней строки кода, чтобы написать следующее:

    lines(CE.time, predict(nls(CE.obs ~ SSlogis(CE.time, 
    Asym, xmid, scal))), lty = 3, lwd = 4),

Линия не строится, и я получаю эту ошибку:

    Error in nls(y ~ 1/(1 + exp((xmid - x)/scal)), data = xy, start = list(xmid = 
    aux[1L],  : step factor 0.000488281 reduced below 'minFactor' of 0.000976562

который я не могу изменить, даже если поиграю с nls.controls и изменю значение minFactor. Я также получаю это сообщение об ошибке после начальной строки, определяющей мод (часть ##. Mod) для некоторых групп населения.

Кроме того, для некоторых групп я получаю сообщение об ошибке после последней строки кода, в котором сообщается об этом:

    Error in qr.solve(QR.B, cc) : singular matrix 'a' in solve

Я не могу придумать никакой рационализации для естественного преобразования данных в журнал, и мне остается предположить, что я просто изменил данные (в данном случае произвольно зарегистрировал их) таким образом, чтобы разрешить функции predic () и SSlogis () функции, чтобы функционировать должным образом, но я не знаю почему. Мне не удалось найти подходящих ответов ни на одном форуме по такой проблеме. Любая помощь будет принята с благодарностью.

* Обновление: я попытался реализовать функцию nlsLM в соответствии с рекомендациями Roland (см. ниже). Это действительно очищает часть кода с запутанным использованием log ():

    lines(CE.time, predict(nlsLM(CE.obs ~ Asym/(1 + exp((xmid - CE.time)/scal)), start 
    = list(Asym = max(CE.obs), xmid = popsizetime[1], scal = 1), control = 
    nls.lm.control(maxiter = 1000))

Однако для других групп я сталкиваюсь с тем же сообщением об ошибке, что и выше, в начальной спецификации модели:

    ChMa.mod = nls(ChMa.observed ~ SSlogis(t.ChestnutMannikin, Asym, xmid, scal))

    Error in nls(y ~ 1/(1 + exp((xmid - x)/scal)), data = xy, start = list(xmid = 
    aux[1L],  : step factor 0.000488281 reduced below 'minFactor' of 0.000976562

Переключено на:

    ChMa.mod = nlsLM(ChMa.observed ~ Asym/(1 + exp((xmid - t.ChestnutMannikin)/
scal)), start = list(Asym = max(ChMa.obs), xmid = popsizetime[2], 
scal = 1), control = nls.lm.control(maxiter = 1000))

Где

    ChMa.observed = c(4.02785074, 0.33847154, 0.99029776, 2.86516540, 0.59588068, 
    0.01334333, 2.07693362, 0.62485994, 3.48979515, 3.67785202, 20.84180181)
    t.ChestnutMannikin = c(1:11)
    popsizetime[2] = 11

Хотя этот переключатель позволяет избежать сообщения об ошибке, nlsLM оценивает функцию, но не оценивает градиент. Без оценки градиента я не могу использовать код se.fit и, следовательно, не могу получить оценку дисперсии для построения графика.


person Nigel Stackhouse    schedule 15.01.2013    source источник
comment
В анализе выживаемости практика log () - ging «время» позволяет получить модели ускоренного времени отказа. Интересно, наткнулись ли вы на случай, когда это подходящая стратегия?   -  person IRTFM    schedule 15.01.2013
comment
Это интересная мысль, но если я аналогичным образом зарегистрирую другой компонент «времени», t.CattleEgret, процесс снова завершится ошибкой. Кроме того, я могу избежать регистрации «времени», если я буду регистрировать компонент «изобилие» (CE.obs).   -  person Nigel Stackhouse    schedule 15.01.2013
comment
Старайтесь не использовать функцию самозапуска, используя другие начальные значения, и используйте nlsLM из пакета minpack.lm. Я попробовал, и мне показалось, что это сработало, но у меня нет времени проверить, верен ли результат, поэтому я не публикую ответ.   -  person Roland    schedule 15.01.2013
comment
Я начал играть с этим на rpubs.com/bbolker/3652. Я надеюсь показать, как использовать AD Model Builder / R2admb для решения этой проблемы.   -  person Ben Bolker    schedule 20.01.2013
comment
Между прочим, разница здесь не между nls и nlsLM; nls отлично работает с другим набором начальных условий.   -  person Ben Bolker    schedule 21.01.2013


Ответы (1)


Я нашел ответ на свои проблемы: мне нужно добавить компонент моей модели, который генерирует градиент для функции, которую я регрессирую с помощью nlsLM.

    log.model = function(t.RedventedBulbul, Asym, xmid, scal) {
            numericDeriv(quote(Asym/(1 + exp((xmid - t.RedventedBulbul)/scal))),
            c("Asym", "xmid", "scal"), parent.frame())
    }
person Nigel Stackhouse    schedule 17.01.2013