二项式数据的回归克里金插值

Question

我使用gstat来预测二项式数据，但是预测值超过了1或者低于0。有人知道如何解决这个问题吗？谢谢。

data(meuse)
data(meuse.grid)
coordinates(meuse) <- ~x+y
coordinates(meuse.grid) <- ~x+y
gridded(meuse.grid) <- TRUE

#glm model
glm.lime <- glm(lime~dist+ffreq, meuse, family=binomial(link="logit"))
summary(glm.lime)

#variogram of residuals
var <- variogram(lime~dist+ffreq, data=meuse)
fit.var <- fit.variogram(var, vgm(nugget=0.9, "Sph", range=sqrt(diff(meuse@bbox\[1,\])^2 + diff(meuse@bbox\[2,\])^2)/4, psill=var(glm.lime$residuals)))   
plot(var, fit.var, plot.nu=T)

#universal kriging
kri <- krige(lime~dist+ffreq, meuse, meuse.grid, fit.var)
spplot(kri[1])

- Geo-sp

Answer 1

一般来说，使用这种回归克里金方法并不能保证模型的有效性，因为趋势和残差的计算是分开的。关于您的代码，需要注意的是，您使用variogram来计算残差半变异函数，但variogram使用正常的线性模型来计算趋势和残差。您需要从glm中确定残差，然后基于此计算残差半变异函数。

你可以手动完成这个过程，或者查看GSIF包中的fit.gstatModel函数。您还可以查看geoRglm包中的binom.krige。 R-sig-geo上的该主题也可能会很有趣：

“从GLM中获取残差与使用指示变量有所不同。此外，根据获取的GLM残差类型可能会有所不同。运行GLM并通过半变异函数探索残差是我认为的例行实践，但这并不总是能告诉您整个故事。拟合广义线性地理统计模型（GLGM）可能更具决定性，因为您可以对模型参数进行推断并更客观地访问空间项的相关性。这是geoRglm一次性完成所有建模而不是通过两个步骤（如拟合没有相关性的模型，然后建模残差）的初衷。这带来了校准MCMC算法的额外负担。后来，spBayes出现了，并且似乎提出了一个更通用的框架，而geoRglm则更具体于单变量二项式和泊松模型。”

正如Roger所说，还可以尝试使用其他替代方案，例如GLMM或MCMCpack，但这肯定不是“开箱即用”的，代码需要适应空间目的。

- Paul Hiemstra