Почему прогнозируемые значения моей GLM цикличны? - PullRequest
Новая
       30

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

2 голосов
/ 25 мая 2020

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

    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

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

1 Ответ

2 голосов
/ 25 мая 2020

Для того, чтобы 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 гораздо более надежным.) (Также: подгонка двух параметров к четырем точкам данных теоретически нормально, но вы не должны пытаться получить sh результаты слишком далеко ... в в частности, многие результаты по умолчанию из GLM / MLE являются асимптотическими c ...)

На самом деле, дважды проверяя соответствие соответствия 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")

enter image description here

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

...