为什么我的GLM的预测值是周期性的？

Question

我写了一个二项式回归模型来预测火成岩石的普遍性。v，在一个考古遗址，根据靠近河流的程度。river_dist但当我使用predict()函数时，我得到的是奇怪的周期性结果，而不是我所期望的曲线。作为参考，我的数据。

    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

为什么这些值会这样循环上升和下降 当绘制图表时产生疯狂的尖峰？

Answer 1

为了让 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 是更稳健的。" (另外：对四个数据点拟合两个参数在理论上是可以的，但不要把结果推得太远......特别是，GLMMLE的很多默认结果都是渐变的。) (另：对四个数据点拟合两个参数，理论上是可以的，但不要试图把结果推得太远......特别是GLMMLE的很多默认结果都是渐变的......)

其实，在仔细检查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")

最后，请注意，你的第三个最远地点是一个离群值（尽管样本量小，这并不影响）。