现在遇到几个问题，求助各位。目前，这几个问题以注释/注意的形式放在正文中
- nlme 包的函数
Variogram() 拟合结果中，n.pairs 总对数不是 157*156/2 = 12246 而是 156*155/2 = 12090，不知道是哪个地方有问题？函数 Variogram() 帮助文档显示 n*(n-1)/2 对的。
- 我将之前的简单克里金插值（无权重约束）换成普通克里金插值（带权重约束），公式是从《Statistics for Spatial Data》(Cressiee 1993) 第 122 页公式 3.2.16 仿写的。全书的数学公式符号与 Cressiee 的书不同，我看了书籍材料之后，手推了一遍，手推仿写的可能有问题，也可能是我的实现代码有问题，理由是普通克里金插值预测方差图看着有问题，见本书图 19.17 的上子图。不知道能否帮我更正？
我先自报家门，回答两个大家可能关心的问题。也想和大家探讨，我这写作的思路是否合理？
- 其实，有很多 R 包可以用来介绍空间数据分析的内容，为什么选用 nlme 包和 spaMM 包？
空间随机效应是随机效应的一种，本章想将所有内容纳入广义线性混合效应模型框架下，参数的估计和预测都是基于频率派的。nlme 包可以求解带空间随机效应的线性混合效应模型，nlme 包高亮突出广义最小二乘和限制极大似然法，它们是很基础和核心的内容，因此，本章特意选择 nlme 包，但是对空间数据，nlme 包不支持预测，这部分就要手推公式、写代码实现普通克里金插值预测的过程。最后，与 spaMM 包实现的空间广义线性混合效应模型比较。spaMM 包基于似然函数推断的（属于频率派的），限制极大似然法估计方差-协方差矩阵的参数，nlme 包也是。nlme 包是 R 内置的，spaMM 包计算很高效，依赖都相对比较少且轻，安装方便。
- 从变差函数拟合、网格划分、克里金插值预测、画图，10-20 行代码就可以了，举例如下，为什么书中写了那么多？
因不想做调包侠，想搞通一条道路。计算过程都被 R 包函数封装了，而且不同的 R 包大多是不同的人开发的，各有自己的一套理论视角，以 nlme 包作为参照，有的参数估计结果与之完全一样，有的需要做一些转化才能获得一样的结果，有的方法不同只能获得近似的结果，更别提还有贝叶斯派方法。
library(gstat)
# 变差模型结构
vario <- variogram(object = log(counts / time) ~ 1, data = rongelap_sf)
# 变差模型的协方差结构和参数初值
fit_vgm <- vgm(psill = 1, model = "Exp", range = 200, nugget = 1)
fit_vgm
# 拟合变差
fit_vario <- fit.variogram(object = vario, model = fit_vgm)
fit_vario
library(abind)
library(stars)
rongelap_coastline_sf <- st_as_sf(rongelap_coastline, coords = c("cX", "cY"), dim = "XY")
rongelap_coastline_sfp <- st_cast(st_combine(st_geometry(rongelap_coastline_sf)), "POLYGON")
rongelap_coastline_buffer <- st_buffer(rongelap_coastline_sfp, dist = 50)
rongelap_grid_stars <- st_bbox(rongelap_coastline_buffer) |>
st_as_stars(nx = 150, ny = 75) |>
st_crop(rongelap_coastline_sfp)
rongelap_grid_stars
# 预测
fit_krige_pred <- krige(
formula = log(counts / time) ~ 1, locations = rongelap_sf,
newdata = rongelap_grid_stars, model = fit_vario
)
fit_krige_pred
library(ggplot2)
# 预测值
ggplot() +
geom_stars(
data = fit_krige_pred, aes(fill = exp(var1.pred)),
na.action = na.omit
) +
scale_fill_viridis_c(
option = "C", breaks = 0:12,
guide = guide_colourbar(
barwidth = 1, barheight = 15
)
) +
theme_bw() +
labs(x = "横坐标", y = "纵坐标", fill = "辐射强度")
# 预测方差
ggplot() +
geom_stars(
data = fit_krige_pred, aes(fill = var1.var),
na.action = na.omit
) +
scale_fill_viridis_c(
option = "C", breaks = 0.1 * 0:12 / 4,
guide = guide_colourbar(
barwidth = 1, barheight = 15
)
) +
theme_bw() +
labs(x = "横坐标", y = "纵坐标", fill = "预测方差")
