为什么我的 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 起作用，您必须将变量指定为 'raw' 值而不是 $:

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 所建议的，您还必须在 newdata.

中为 n 指定一个（或多个）值

此模型似乎与二项式回归相同：是否有理由不使用

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

最后，请注意您的第三远站点是一个异常值（尽管样本量小意味着它不会造成太大伤害）。

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

Why are the predicted values of my GLM cyclical?

statistics

r

graph

glm

mle