Inconsistent results between custom function using terra and a simplified version to calculate global statistic of a raster

I am running an access poverty analysis on a raster of travel times. Using a slightly altered version of the Foster-Greer-Thorbecke index with exponent Alpha = 0, the following equality holds:

Here y is a real number indicating travel distance from pixel i to, say, the next hospital. z is a threshold above which a pixel is coded as access-poor. I is the number of all pixels in the raster, indexed by i from the smallest to the largest travel time, and i-hat is the first pixel in that order that has a travel time higher than z (i.e. it's the "marginally access-poor" pixel).

I tried to implement this in R using terra (Apologies for not offering a reprex here, but somehow the reprex package didn't handle the raster operations):

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

^{Created on 2024-11-05 with reprex v2.1.1.9000}

Setting the threshold to the global median should, by definition, yield 0.5 considering that half of the values in the raster fall above (below) the threshold. The last three lines of code above confirm that.

Can you please help me understand why my function fgt(), as I defined it above, gives the wrong result? In fact, with alpha = 0, it always returns 1. Why is that?

Update The answer was rather simple. Let's start by taking a closer look at the expression ((upper-z)/z)^alpha. It can be split up into the following three steps

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

Turns out the third line applies the exponent to all cells, whether they are NA or not. If alpha > 0 it gives the correct answer, i.e. that NA to the power of alpha is still NA. However, it seems that "anything to the power of zero equals 1" is hard-coded, and so it yields 1 for any base with exponent zero...including NAs. This is why I get N times 1, which divided by N gives 1.

The solution is to only apply the relevant part of the function to the non-missing cells, e.g. by using an if-else condition.

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

^{Created on 2024-11-05 with reprex v2.1.1.9000}