Ошибка: начальное значение в 'vmmin' не является конечным не в mle2(), а в confint()

Я знаю, что в Интернете полно вопросов (и ответов) об ошибке 'initial value in vmmim is not finite' при попытке подогнать параметры для объекта mle2. У меня нет этой ошибки при создании моего объекта mle2, но я получаю эту ошибку при попытке найти 95% CI для параметра из объекта mle2.

Вот воспроизводимый пример:

Вот данные:

d = structure(list(SST_1YR = c(11.6, 11.7, 11.9, 12, 12.1, 12.2, 
12.3, 12.4, 12.5, 12.6, 12.7, 12.8, 12.9, 13, 13.1, 13.2, 13.3, 
13.4, 13.5, 13.6, 13.7, 13.8, 13.9, 14, 14.2, 14.3, 14.4, 14.5, 
14.6, 14.7, 14.8, 14.9, 15, 15.1, 15.2, 15.3, 15.4, 15.5, 15.6, 
15.7, 15.8, 15.9, 16, 16.2, 16.3, 16.5, 16.6, 16.7, 16.9, 17, 
17.1, 17.2, 17.3, 17.4, 17.5, 17.6, 17.7, 17.8, 17.9), DML = structure(c(84.5, 
71, 114.75, 90.9473684210526, 31.7631578947368, 92.5, 80.4, 98.7021276595745, 
70.8, 66.8382352941177, 70.2553191489362, 98.1111111111111, 86.5241379310345, 
59.7209302325581, 38.7692307692308, 78.2028985507246, 86.3503649635037, 
69.1161290322581, 61.9122807017544, 60.1212121212121, 98.5490196078431, 
94.3145161290323, 76.5643564356436, 39.4230769230769, 98.42, 
95.6129032258064, 65.9673202614379, 39, 64.0576923076923, 42.4166666666667, 
59.6989247311828, 62.8039215686275, 74.5263157894737, 50.8888888888889, 
64.35, 40.5, 53.7466666666667, 42, 49.5, 23.8888888888889, 39.6170212765957, 
74.8947368421053, 42.8518518518519, 40.0344827586207, 53, 39.3333333333333, 
24.1333333333333, 30, 39.4880952380952, 94.4883720930233, 69.1428571428571, 
33.7179487179487, 26.1538461538462, 37.8965517241379, 38.4117647058824, 
44.2727272727273, 68.3157894736842, 37.3, 43.4444444444444), .Dim = 59L, .Dimnames = list(
    c("11.6", "11.7", "11.9", "12", "12.1", "12.2", "12.3", "12.4", 
    "12.5", "12.6", "12.7", "12.8", "12.9", "13", "13.1", "13.2", 
    "13.3", "13.4", "13.5", "13.6", "13.7", "13.8", "13.9", "14", 
    "14.2", "14.3", "14.4", "14.5", "14.6", "14.7", "14.8", "14.9", 
    "15", "15.1", "15.2", "15.3", "15.4", "15.5", "15.6", "15.7", 
    "15.8", "15.9", "16", "16.2", "16.3", "16.5", "16.6", "16.7", 
    "16.9", "17", "17.1", "17.2", "17.3", "17.4", "17.5", "17.6", 
    "17.7", "17.8", "17.9")))), .Names = c("SST_1YR", 
"DML"), row.names = c(NA, -59L), class = "data.frame")

Вот создание объекта mle2 (без предупреждений...)

m = mle2(DML~dgamma(scale=(a+b*SST_1YR)/sh, shape=sh), start=list(a=170, b=-7.4, sh=10), data=d)

И здесь я получаю NA и предупреждение vmmin для нижней границы параметра b:

confint(m)

Я пытался изменить начальные значения, но ничего из того, что я пробовал, не помогло. Я создал другие модели с теми же данными, но с другим распределением и без ошибок. Может ли кто-нибудь помочь мне понять, что вызывает эту ошибку?

Использование пакета bbmle-1.0.17


r mle
person CephBirk    schedule 12.12.2014    source источник

Ответы (1)


arrow_upward
3
arrow_downward

Здесь есть несколько вещей, которые можно попробовать. Сначала посмотрите на данные (всегда хорошая идея):

library("ggplot2"); theme_set(theme_bw())
ggplot(d,aes(SST_1YR,DML)) + geom_point()+
    geom_smooth(method="glm",family=Gamma(link="identity"))+
        geom_smooth(method="lm",colour="red",fill="red")

Обратите внимание, что в этом случае гамма-регрессия выглядит почти идентично обычной линейной регрессии (т. е. параметр формы большой). Кроме того, распределение значений x далеко от начала координат — это может привести к числовым проблемам.

library("bbmle")
m <- mle2(DML~dgamma(scale=(a+b*SST_1YR)/sh, shape=sh),
          start=list(a=170, b=-7.4, sh=10), data=d)
confint(m)

Подтверждает проблему:

##        2.5 %     97.5 %
## a  132.05952 203.192159
## b         NA  -4.407289
## sh   6.83566  13.933383

Я думал, что установка parscale может помочь, но похоже, что проблема усугубляется, а не улучшается:

m2 <- update(m,control=list(parscale=c(a=170,b=8,sh=10)))
confint(m2)
##       2.5 %     97.5 %
## a        NA 203.153230
## b        NA  -4.407281
## sh 6.835659  13.933383

Помогает ли центрирование переменной-предиктора? scale(x,scale=FALSE) центрирует, но не масштабирует x ... (использование SST_1YR-mean(SST_1YR) может быть понятнее, тогда у нас не будет трех плавающих scale в выражении...

m3 <- mle2(DML~dgamma(scale=(a+b*scale(SST_1YR,scale=FALSE))/sh, shape=sh),
          start=list(a=170, b=-7.4, sh=10), data=d)

confint(m3)
##       2.5 %    97.5 %
## a  56.462610 66.754118
## b  -9.421521 -4.407262
## sh  6.835662 13.933384

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

Оказывается, вы также можете установить эту модель через

glm(DML~SST_1YR,family=Gamma(link="identity"),data=d)

хотя confint() снова терпит неудачу довольно загадочно (Error in y/mu: non-conformable arrays).

Некоторые другие вещи, которые я пробовал, но не работали особенно хорошо (включены сюда только для полноты):

  1. попытайтесь предотвратить отрицательную линейную регрессию:
mle2(DML~dgamma(scale=pmin((a+b*SST_1YR)/sh,1e-5),
                      shape=sh),
          start=list(a=170, b=-7.4, sh=10), data=d)
  1. используйте оштрафованную форму dgamma, чтобы возвращать плохие вероятности, а не NA, когда x<0:
dgamma_pen <- function(x,...,log=FALSE) {
   r <- if (x<0) (-100) else dgamma(x,...,log=TRUE)
   if (log) r else exp(r)
}

m4 <- mle2(DML~dgamma_pen(scale=pmin((a+b*SST_1YR)/sh,1e-5),
                    shape=sh),
         start=list(a=170, b=-7.4, sh=10), data=d)
person Ben Bolker    schedule 13.12.2014
comment
Замечательно! Я бы проголосовал дважды, если бы мог! Благодарю вас! - person CephBirk; 14.12.2014
comment
Я нашел этот ответ в поиске описания связей идентичности с гамма-ошибками, что показалось мне действительно странной комбинацией. Не сталкивается ли он с трудностями при действительно заниженных оценках средств? - person IRTFM; 09.01.2015