Почему прогнозируемые значения моего GLM цикличны?

Я написал модель биномиальной регрессии для прогнозирования распространенности магматических камней v на месте археологических раскопок на основе близости к реке river_dist, но когда я использую функцию predict(), я получаю странные циклические результаты вместо кривой I. ожидал. Для справки мои данные:

    v   n river_dist
1 102 256       1040
2   1  11        720
3  19  24        475
4  12  15        611

Что мне подходит к этой модели:

library(bbmle)
m_r <- mle2(ig$v ~ dbinom(size=ig$n, prob = 1/(1+exp(-(a + br * river_dist)))),
    start = list(a = 0, br = 0), data = ig)

Это дает коэффициент, который при обратном преобразовании предполагает снижение вероятности появления изверженных камней на метр от реки примерно на 0,4% (br = 0,996):

exp(coef(m_r))

Это все хорошо. Но когда я пытаюсь предсказать новые значения, я получаю странную цикличность значений:

newdat <- data.frame(river_dist=seq(min(ig$river_dist), max(ig$river_dist),len=100))
newdat$v <- predict(m_r, newdata=newdat, type="response")
plot(v~river_dist, data=ig, col="red4")
lines(v ~ river_dist, newdat, col="green4", lwd=2)

Пример прогнозируемых значений:

   river_dist          v
1     475.0000 216.855114
2     480.7071   9.285536
3     486.4141  20.187424
4     492.1212  12.571487
5     497.8283 213.762248
6     503.5354   9.150584
7     509.2424  19.888471
8     514.9495  12.381805
9     520.6566 210.476312
10    526.3636   9.007289
11    532.0707  19.571218
12    537.7778  12.180629

Почему значения циклически меняются вверх и вниз, создавая сумасшедшие всплески на графике?


r graph statistics glm mle
person Lauren Pratt    schedule 25.05.2020    source источник
comment
кажется, что прогнозы циклически повторяются с шагом в четыре, когда вы передаете четыре строки данных - так что, возможно, он циклически повторяет n. Для прогнозов вы можете сделать plogis(tcrossprod(coef(m_r), cbind(1, newdat$river_dist))), но это не ответит на ваш вопрос.   -  person user20650    schedule 25.05.2020
comment
... поэтому попробуйте использовать newdat$n = 1 . ps вам не нужно использовать ig$, так как вы используете data=, т.е. используйте mle2(v ~ dbinom(size=n, ...   -  person user20650    schedule 25.05.2020
comment
@ user20650, не хотел воровать это. (Для комментариев, которые являются такими ясными, вы также можете опубликовать ответ...)   -  person Ben Bolker    schedule 25.05.2020

Ответы (1)


arrow_upward
2
arrow_downward

Чтобы newdata работал, вы должны указать переменные как «сырые» значения, а не с $:

library(bbmle)
m_r <- mle2(v ~ dbinom(size=n, prob = 1/(1+exp(-(a + br * river_dist)))),
    start = list(a = 0, br = 0), data = ig)

На этом этапе, как предлагает @user20650, вам также необходимо указать значение (или значения) для n в newdata.

Эта модель кажется идентичной биномиальной регрессии: есть ли причина не использовать

glm(cbind(v,n-v) ~ river_dist, data=ig, family=binomial)

? (bbmle:mle2 является более общим, но glm гораздо более надежным.) (Кроме того: подгонка двух параметров к четырем точкам данных теоретически хороша, но вы не должны пытаться зайти слишком далеко... в частности, многие значения по умолчанию результаты GLM/MLE асимптотичны...)

На самом деле, перепроверив соответствие MLE и GLM, я понял, что метод по умолчанию ("BFGS" по историческим причинам) на самом деле не дает правильного ответа (!); переключение на method="Nelder-Mead" улучшает ситуацию. Добавление control=list(parscale=c(a=1,br=0.001)) к списку аргументов, или масштабирование расстояния реки (например, переход от «1 м» к «100 м» или «1 км» в качестве единицы измерения) также решит проблему.

m_r <- mle2(v ~ dbinom(size=n,
        prob = 1/(1+exp(-(a + br * river_dist)))),
            start = list(a = 0, br = 0), data = ig,
            method="Nelder-Mead")
pframe <- data.frame(river_dist=seq(500,1000,length=51),n=1)
pframe$prop <- predict(m_r, newdata=pframe, type="response")
CIs <- lapply(seq(nrow(ig)),
              function(i) prop.test(ig[i,"v"],ig[i,"n"])$conf.int)
ig2 <- data.frame(ig,setNames(as.data.frame(do.call(rbind,CIs)),
              c("lwr","upr")))
library(ggplot2); theme_set(theme_bw())
ggplot(ig2,aes(river_dist,v/n))+
    geom_point(aes(size=n)) +
    geom_linerange(aes(ymin=lwr,ymax=upr)) +
    geom_smooth(method="glm",
                method.args=list(family=binomial),
              aes(weight=n))+
    geom_line(data=pframe,aes(y=prop),colour="red")

введите здесь описание изображения

Наконец, обратите внимание, что ваш третий по дальности сайт является исключением (хотя небольшой размер выборки означает, что это не сильно повредит).

person Ben Bolker    schedule 25.05.2020
comment
как вы выбираете значения parscale, пожалуйста? Это делается после подгонки модели один раз, а затем масштабирования при втором прогоне или? - person user20650; 25.05.2020
comment
вы можете запустить и масштабировать или просто установить parscale примерно на 1/<typical value of predictor> - person Ben Bolker; 25.05.2020
comment
Вау, это потрясающе! Благодарю вас! Я подозреваю, что мой первоначальный инструктор по статистике R научил нас mle2, поэтому мы были вынуждены явно думать о наших начальных значениях; здорово знать о функции glm на будущее. И этот график великолепен! Ты спасатель. - person Lauren Pratt; 26.05.2020
comment
Дополнительное примечание: есть ли преимущество в том, чтобы делать прогноз на 500-1000 м, а не на реальный минимум и максимум в моих данных? (475 и 1040 м?) - person Lauren Pratt; 26.05.2020
comment
дело вкуса/предпочтения, я думаю. - person Ben Bolker; 27.05.2020