r - 为什么我的 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 更健壮。)(另外:将两个参数拟合到四个数据点理论上是可以的，但你不应该试图将结果推得太远......特别是， GLM/MLE 的许多默认结果是渐近的......)

实际上，在仔细检查 MLE 与 GLM 的对应关系时，我意识到默认方法(“BFGS”，出于历史原因)实际上并没有给出正确的答案(!)；切换到 method="Nelder-Mead"改善事情。添加 control=list(parscale=c(a=1,br=0.001))到参数列表，或缩放河流距离(例如从“1 m”到“100 m”或“1 km”作为单位)，也可以解决问题。

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")

最后，请注意，您的第三远站点是一个异常值(尽管样本量小意味着它不会受到太大影响)。