使用 terra 的自定义函数与计算栅格全局统计数据的简化版本之间的结果不一致

Question

我正在对旅行时间栅格进行访问贫困分析。使用指数 Alpha = 0 的 Foster-Greer-Thorbecke 指数 的稍微修改版本，以下等式成立：

这里 y 是一个实数，表示从像素 i 到下一家医院的行进距离。 z 是一个阈值，高于该阈值的像素将被编码为访问不良。 I 是栅格中所有像素的数量，按 i 从最小到最大行进时间索引，i-hat 是按该顺序的第一个行进时间高于 z 的像素（即，它是“边际访问-可怜的”像素）。

我尝试使用 terra 在 R 中实现这一点（很抱歉在这里没有提供 reprex，但不知何故 reprex 包没有处理光栅操作）：

library(terra)
#> Warning: package 'terra' was built under R version 4.3.3
#> terra 1.7.81

# make example raster
r <- rast(ncol=10, nrow=10)
set.seed(0)
values(r) <- runif(ncell(r))

# mid point of raster values
median_val <- global(r, median)$global

fgt <- function(alpha, z, raster){
  
  require(terra)
  
  # get number of non-missing pixels
  I <- global(raster, "notNA")
  
  # recode to only retain pixels above poverty threshold
  upper <- ifel(raster > z, raster, NA)
  
  # calculate poverty gap ratio for pixels above poverty threshold
  gap <- identity(((upper-z)/z)^alpha)
  
  # sum up poverty gap ratios across all poor pixels
  sumsign <- global(x = gap, fun = "sum", na.rm = T)
  
  # divide by number of all (non-missing) pixels
  sumsign/I 
}


# calculate FGT index (alpha = 0) for median-value threshold
fgt(0, median_val, r)
#>       sum
#> lyr.1   1


# alternative way to compute head count ratio
poor <- global(ifel(r > median_val, r, NA), "notNA")

all <- global(r, "notNA")

poor/all
#>       notNA
#> lyr.1   0.5

^{创建于 2024 年 11 月 5 日，使用 reprex v2.1.1.9000}

根据定义，考虑到栅格中一半的值高于（低于）阈值，将阈值设置为全局中位数应产生 0.5。上面的最后三行代码证实了这一点。

您能帮我理解为什么我的函数

fgt()

（正如我上面定义的那样）给出了错误的结果吗？事实上，当 alpha = 0 时，它总是返回 1。这是为什么？

Answer 1

更新 答案相当简单。让我们首先仔细看看表达式

((upper-z)/z)^alpha

。可以分为以下三个步骤

a <- upper-z # this works as expected
b <- a/z # this also works as expected
c <- b^alpha # this is where the issue lies.

结果第三行将指数应用于所有单元格，无论它们是否为

NA

。如果 alpha > 0，则给出正确答案，即 NA 的 alpha 次方仍然是 NA。然而，似乎“任何零次方都等于 1”是硬编码的，因此对于任何指数为零的基数，它都会产生 1...包括 NA。这就是为什么我得到 N 乘以 1，除以 N 得到 1。

解决方案是仅将函数的相关部分应用于非缺失单元格，例如通过使用 if-else 条件。

library(terra)
#> Warning: package 'terra' was built under R version 4.3.3
#> terra 1.7.81

# make example raster
r <- rast(ncol=10, nrow=10)
set.seed(0)
values(r) <- runif(ncell(r))

# mid point of raster values
median_val <- global(r, median)$global

fgt <- function(alpha, z, raster){
  
  require(terra)
  
  # get number of non-missing pixels
  I <- global(raster, "notNA")
  
  # recode to only retain pixels above poverty threshold
  upper <- ifel(raster > z, raster, NA)
  
  # UPDATED: if-else statement solves the issue!!!
  gap <- ifel(is.na(upper), NA, ((upper-z)/z)^alpha)
  
  # sum up poverty gap ratios across all poor pixels
  sumsign <- global(x = gap, fun = "sum", na.rm = T)
  
  # divide by number of all (non-missing) pixels
  sumsign/I 
}


# calculate FGT index (alpha = 0) for 50 minute threshold
fgt(0, median_val, r)
#>       sum
#> lyr.1 0.5

^{创建于 2024 年 11 月 5 日，使用 reprex v2.1.1.9000}