在 ggplot 中绘制经验和拟合半变异函数

Question

我正在使用半变异函数研究数据中的空间自相关。我的数据：

Response <- c(21L, 36L, 30L, 29L, 30L, 45L, 100L, 0L, 0L, 0L, 0L, 0L, 59L, 
18L, 24L, 23L, 26L, 29L, 23L, 21L, 14L, 30L, 43L, 14L, 8L, 0L, 
0L, 0L, 0L, 0L, 23L, 38L, 20L, 28L, 45L, 21L, 46L, 23L, 6L, 4L, 
0L, 0L, 0L, 0L, 0L, 17L, 10L, 41L, 24L, 31L, 16L, 23L, 31L, 6L, 
2L, 0L, 0L, 0L, 0L, 0L, 8L, 20L, 18L, 18L, 40L, 9L, 1L, 25L, 
4L, 34L, 0L, 0L, 0L, 0L, 0L, 39L, 8L, 7L, 22L, 16L, 18L, 23L, 
11L, 25L, 28L, 0L, 0L, 0L, 0L, 0L, 3L, 22L, 11L, 9L, 123L, 50L, 
12L, 1L, 46L, 1L, 4L, 1L, 2L, 0L, 37L)

Covar1 <- structure(c(1L, 3L, 1L, 1L, 3L, 3L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 
3L, 3L, 1L, 3L, 1L, 1L, 3L, 3L, 1L, 1L, 3L, 3L, 2L, 2L, 2L, 2L, 
2L, 1L, 3L, 1L, 1L, 3L, 3L, 1L, 1L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 
1L, 3L, 1L, 1L, 3L, 3L, 1L, 1L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 1L, 
3L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 3L, 
3L, 3L, 3L, 1L, 1L, 3L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 1L, 
3L, 1L, 1L, 3L, 3L, 1L, 2L, 2L, 2L, 2L, 2L, 1L), .Label = c("A", 
"B", "C"), class = "factor")

Covar2 <- structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L), .Label = c("1", 
"2", "3", "4", "5", "6", "7"), class = "factor")

df <- data.frame(Response, Covar1, Covar2)

我运行一个简单的模型，使用 gstat 从残差和空间坐标制作经验半变异函数和拟合半变异函数，并绘制它们：

mod1 <- glm(Response ~ Covar1 * Covar2, data = df)

geo <- as.data.frame(resid(mod1))

geo$x <- c(34.59481, 34.60548, 34.59825, 34.59039, 34.56546, 34.56749, 
34.5964, 34.40986, 34.40083, 34.39536, 34.41291, 34.40512, 34.36381, 
34.35335102, 34.32548, 34.59481, 34.60548, 34.59825, 34.59039, 
34.56749, 34.56546, 34.5964, 34.36381, 34.35335102, 34.32548, 
34.41291, 34.40986, 34.40512, 34.40083, 34.39536, 34.59481, 34.60548, 
34.59825, 34.59039, 34.56749, 34.56546, 34.5964, 34.36381, 34.35335102, 
34.32548, 34.41291, 34.40986, 34.40512, 34.40083, 34.39536, 34.59481, 
34.60548, 34.59825, 34.59039, 34.56749, 34.56546, 34.5964, 34.36381, 
34.35335102, 34.32548, 34.41291, 34.40986, 34.40512, 34.40083, 
34.39536, 34.59481, 34.60548, 34.59825, 34.59039, 34.36381, 34.35335102, 
34.32548, 34.56749, 34.56546, 34.5964, 34.41291, 34.40986, 34.40512, 
34.40083, 34.39536, 34.36381, 34.35335102, 34.32548, 34.56546, 
34.56749, 34.5964, 34.59481, 34.60548, 34.59825, 34.59039, 34.41291, 
34.40986, 34.40512, 34.40083, 34.39536, 34.32548, 34.35335102, 
34.59481, 34.60548, 34.59039, 34.59825, 34.56546, 34.56749, 34.5964, 
34.41291, 34.40986, 34.40512, 34.40083, 34.39536, 34.36381)

geo$y <- c(-2.18762, -2.18308, -2.16174, -2.16018, -2.14787, -2.15296, 
-2.12863, -2.14325, -2.14552, -2.1454, -2.13926, -2.14652, -2.12463, 
-2.121925978, -2.10213, -2.18762, -2.18308, -2.16174, -2.16018, 
-2.15296, -2.14787, -2.12863, -2.12463, -2.121925978, -2.10213, 
-2.13926, -2.14325, -2.14652, -2.14552, -2.1454, -2.18762, -2.18308, 
-2.16174, -2.16018, -2.15296, -2.14787, -2.12863, -2.12463, -2.121925978, 
-2.10213, -2.13926, -2.14325, -2.14652, -2.14552, -2.1454, -2.18762, 
-2.18308, -2.16174, -2.16018, -2.15296, -2.14787, -2.12863, -2.12463, 
-2.121925978, -2.10213, -2.13926, -2.14325, -2.14652, -2.14552, 
-2.1454, -2.18762, -2.18308, -2.16174, -2.16018, -2.12463, -2.121925978, 
-2.10213, -2.15296, -2.14787, -2.12863, -2.13926, -2.14325, -2.14652, 
-2.14552, -2.1454, -2.12463, -2.121925978, -2.10213, -2.14787, 
-2.15296, -2.12863, -2.18762, -2.18308, -2.16174, -2.16018, -2.13926, 
-2.14325, -2.14652, -2.14552, -2.1454, -2.10213, -2.121925978, 
-2.18762, -2.18308, -2.16018, -2.16174, -2.14787, -2.15296, -2.12863, 
-2.13926, -2.14325, -2.14652, -2.14552, -2.1454, -2.12463)

library(sp)     
names(geo) <- c("resids", "x", "y")
coordinates(geo) <- ~ x + y
proj4string(geo) <- CRS("+proj=longlat +datum=WGS84")

library(gstat)
var1 <- variogram(resids ~ x + y, data = geo)
v.fit1 = fit.variogram(var1, vgm(50, "Exp", 2, 50))
plot(var1, v.fit1)

情节是 class 'trellis'，它不接受标准基础 R 图形的参数，所以我想使用 ggplot 来创建我的图形。我可以绘制我的经验变异函数（仅限点）：

ggplot(var1, aes(x=dist,y=gamma)) +
geom_point()

但是我在绘制拟合模型（直线）时遇到了问题。任何帮助将不胜感激。

Answer 1

您可以通过包 gstat 中的 variogramLine() 为给定的变差函数模型生成半方差值。然后，您可以使用 geom_line() 中的那些来绘制那些拟合值。

我只在 var1 中获取最大距离的值，因此两个数据集的范围将相同。

preds = variogramLine(v.fit1, maxdist = max(var1$dist))
head(preds)
          dist    gamma
1 1.037174e-05 67.13427
2 5.212964e-02 71.24628
3 1.042489e-01 75.23463
4 1.563682e-01 79.10305
5 2.084874e-01 82.85514
6 2.606067e-01 86.49440

现在使用这个新数据集向图中添加一个 geom_line() 图层。 x 和 y 变量的名称与 var1 中的相同，因此您不需要映射任何新的美学。

ggplot(var1, aes(x = dist, y = gamma)) +
     geom_point() +
     geom_line(data = preds)

Answer 2

在 aosmith 非常有见地的回答之后添加一些有用的东西。您可能想在经验变异函数中添加比例点。我设法使用以下代码做到了这一点：

plot_variogram <- function(v, m) {
  preds = variogramLine(m, maxdist = max(v$dist))
  ggplot() + 
    geom_point(data = v, aes(x = dist, y = gamma, size=np)) +
    geom_line(data = preds, aes(x = dist, y = gamma))
}

v <- variogram(zinc ~ 1, meuse)
m <- fit.variogram(v, vgm(c("Exp", "Sph")))
plot_variogram(v, m)

您可以查看此 post 以获得更多主题：Variography with gstat and ggplot2

在 ggplot 中绘制经验和拟合半变异函数

Plotting empirical and fitted semivariogram in ggplot

r

ggplot2

plot

spatial

trellis