如何根据权重从空间网格中采样？

Question

我正在运行一些模拟，并且想要生成空间自相关数据。可能是以完全错误的方式来解决这个问题......

首先，我生成了一个观测数据集，其中变量

success

随着

age

和

country

的变化而变化。因为它很大，所以它看起来像这样：

library(tidyverse)
library(sf)
library(sfdep)

#read in the dataset of simulated data
sim_data <- read_rds("sim_data.rds")

#summarise what the data look like
sim_data %>% 
  head(20)
#>    countries      age prevalence success       age_z
#> 1     Malawi 3.946353 0.16114073       1  0.94643938
#> 2     Malawi 1.538666 0.10983733       0 -1.46124749
#> 3     Malawi 3.627969 0.07027045       0  0.62805528
#> 4     Malawi 3.820259 0.15889337       0  0.82034525
#> 5     Malawi 2.830967 0.13286989       0 -0.16894678
#> 6     Malawi 3.876449 0.07469908       0  0.87653512
#> 7     Malawi 4.738689 0.17526246       0  1.73877511
#> 8     Malawi 2.021715 0.11844667       0 -0.97819858
#> 9     Malawi 2.849171 0.05639000       1 -0.15074259
#> 10    Malawi 4.760058 0.17564332       0  1.76014421
#> 11    Malawi 4.912906 0.16997609       0  1.91299183
#> 12    Malawi 1.469949 0.03180826       0 -1.52996444
#> 13    Malawi 2.899988 0.14249146       1 -0.09992556
#> 14    Malawi 3.241331 0.14018377       1  0.24141710
#> 15    Malawi 4.616126 0.08788227       0  1.61621167
#> 16    Malawi 1.554841 0.11851702       0 -1.44507321
#> 17    Malawi 4.955567 0.17073644       0  1.95565303
#> 18    Malawi 4.786673 0.09092192       0  1.78675905
#> 19    Malawi 1.329750 0.11450525       0 -1.67016365
#> 20    Malawi 3.056847 0.13689573       0  0.05693326

^{创建于 2023-08-04，使用 reprex v2.0.2}

现在我创建一个空间网格，并生成一些空间自相关数据（得到了很多帮助：https://josiahparry.com/posts/csr.html）

#make a spatial grid
grid <- st_make_grid(cellsize = c(1, 1), n = 10, offset = c(0, 0)) |> 
  as_tibble() |> 
  st_as_sf() |> 
  mutate(
    id = row_number(),
    nb = st_contiguity(geometry),
    wt = st_weights(nb)
    )

#generate some spatially correlated data
nb <- grid[["nb"]]
wt <- grid[["wt"]]

x <-  geostan::sim_sar(w = wt_as_matrix(nb, wt), rho = 0.78)

#function to simulate spatial weights by permutation
#example taken from here: https://josiahparry.com/posts/csr.html
permute_nb <- function(nb) {
  # first get the cardinality
  cards <- st_cardinalties(nb)
  # instantiate empty list to fill
  nb_perm <- vector(mode = "list", length = length(nb))
  
  # instantiate an index
  index <- seq_along(nb)
  
  # iterate through and full nb_perm
  for (i in index) {
    # remove i from the index, then sample and assign
    nb_perm[[i]] <- sample(index[-i], cards[i])
  }
  
  structure(nb_perm, class = "nb")
}

#we only simulate 1 weighted grid here for the example
nsim = 1 

sim_spatial <- replicate(
  nsim, {
    nb_perm <- permute_nb(nb)
    st_lag(x, nb_perm, st_weights(nb_perm))
  }
)

#tidy up the grid
sim_spatial <- sim_spatial %>%
  as.data.frame() %>%
  rename_all(funs(paste0("sim_", str_remove(.,"V")))) %>%
  mutate(id = row_number()) %>%
  left_join(grid)

sim_spatial %>%
  ggplot() +
  geom_sf(aes(fill=sim_1, geometry=geometry)) +
  scale_fill_viridis_c()

^{创建于 2023-08-04，使用 reprex v2.0.2}

现在，我想将

sim_data

中的每个观测值随机分配到

spatial_data

中的单元格，其中

1

值为

success

的观测值更有可能分配给具有更高

sim_1

值的单元格.

谢谢！

Answer 1

很难知道您希望如何将

sim_1

列转换为权重（因为其中一些是负数），但现在让我们将它们标准化为 0-1 范围。我们需要跟踪已分配的方格，因此我们也将

id

列拉出到它自己的变量中。

weights <- with(sim_spatial, (sim_1 - min(sim_1))/diff(range(sim_1)))
squares <- sim_spatial$id

现在我们可以使用

sample

作为权重，使用

 误导性命名的 

weights 参数来

prob

所有方块，理想情况下应该将其称为“权重”。我们只会对数据中

success == 1

的数据执行此操作

success1 <- sample(squares, sum(sim_data$success), prob = weights)

现在我们删除“用完”的方块，并使用反转权重再次采样，对于以下情况：

success == 0

：

weights       <- weights[-success1]
squares       <- squares[-success1]
weights       <- 1 - weights
success0      <- sample(squares, sum(sim_data$success == 0), prob = weights)

现在我们将索引放入

sim_data

中的一个新列中，我们可以使用它来

left_join

将结果返回到我们的空间数据。

sim_data$id <- 0
sim_data$id[sim_data$success == 1] <- success1
sim_data$id[sim_data$success == 0] <- success0
sim_spatial <- left_join(sim_spatial, sim_data, by = 'id')

如果我们绘制此图，我们会发现我们更有可能在具有高

sim_1

值的方格中取得成功：

sim_spatial %>%
  ggplot() +
  geom_sf(aes(fill = sim_1, geometry = geometry)) +
  geom_sf(aes(colour = factor(success), geometry = geometry), 
          data = . %>% filter(!is.na(success)), linewidth = 2, fill = NA) +
  scale_fill_viridis_c() +
  scale_color_manual('Success', values = c('cyan', 'red'))

但请注意，因为有很多方格并且其中大多数具有中等权重，所以随机选择很可能会选择中等权重的方格。您可以更改从

sim_1

生成权重的方式 - 例如，使用

10^sim_1

表示

success == 1

，使用

10^(-sim_1)

表示

success == 0

可提供更令人信服的分布：