What is the best way to estimate an error on stacked data? I have 400 pixel maps (which are different realisations of the same underlying system) which, individually, don't show much. If I stack them, however, I get a statistical signal; the maps tend to have slightly lower values in the centre of the map, so a stack will show this dip quite clearly. I evaluate the map by dividing it into 10 linear radial bins starting from the centre, within which I calculate the mean of the pixel values. I want to estimate an "error" on this stacked map.
My method of choice to do so is "delete 1 jackknife resampling", but I am a bit confused how to interpret the variances I get. What I have done so far is to stack 399 maps, always leaving out a different map, so that we get 400 slightly different stacked maps. I then draw a radial profile on all 400 of these maps, so I get a distribution of 400 mean values for each bin; for simplicity, let's only consider the mean values of the first bin for now, which I call $\mu_{i}^{1}$. Calculating the average value of those is then straightforward as $\bar{\mu^{1}} = 1/N \sum_{i=1}^N \mu_{i}^{1}$, but I saw different definitions of the variance and the difference is not clear to me.
$\text{var} (\theta_{\text{jack}}) = \frac{N-1}{N} \sum_{i=1}^{N} (\theta_{(i)} - \theta_{(.)})^2$
is one version, where we can identify $\theta_{(i)}$ with $\mu_{i}^{1}$, and $\theta_{(.)}$ with $\bar{\mu^{1}}$. Another version lists the same equation, but with a factor of $\frac{1}{N-1}$ in front instead. I don't really understand where this difference comes from and how the interpretation differs between them. Generally speaking, the latter definition (obviously) produces far lower estimates for the error. If I plot a "total" radial profile this way, taking the $\bar{\mu^{(i)}}$ as the datapoints and the square root of either of these two definitions of the variance as the errorbars, the latter, smaller errorbars look very "small" compared to the differences between the individual $\bar{\mu^{(i)}}$, which is why I (and the people I am working with) are taking issue with it, however, I still have not found a satisfying answer as to why the definition is not appropriate.
The major question, then: is the current approach sensible for the stated goal, and which of the two variances is more appropriate as an estimate for the "error" of the total stacked map?
Also, what is the difference if I first draw the radial profiles and then use jackknife on them, i.e. I create 400 radial profiles on the individual maps, and then use a jackknife resampling approach on those profiles? Is that even sensible?
 A: I have not so much experience with jackknifing, but e.g. a BCa-boostrap confidence interval is easy enough to construct for each point - if you really need a method for non-normal data (perhaps assuming normally distributed noise would just work fine and you could just calculate a standard error & CI from the multiple values). This is very straightforward, if the different maps are realizations of the same underlying state just with some randomness on top.
Let's for simplicity use data like this (49 realizations/measurements on a 100 by 100 grid - btw. 49 is just because a 7 by 7 plot looks nicer than what I can do for 50):

Then we get an mean like that:

It's easiest to show a confidence interval across a cut along one axis, e.g. at some fixed y-value:

This also illustrates a major downside of this: there's no smoothing of the pattern.
If we want to smooth this out, we could use some model with spatial correlation such as a GAM model. Obviously, this involves more assumptions (e.g. here I assume normally distributed errors, particular correlation structures across space and a very simple single random effect across the realizations) and a choice on how much smoothing one wants (one could of course determine that using cross-validation etc., but it certainly will take longer than the simple mean + bootstrapping or jackknifing).
The model fit I get for the example above looks like this:

A cross-section with CI then looks like this:

Besides smoothing the function out, the other advantage is that by using information from surrounding points we reduce uncertainty (of course assuming the model is correct) and get narrow confidence intervals. However, with the extent of smoothing chosen here, we see that we actually oversmoothed the trough of the curve versus the true pattern we simulated with.
R code is below:
library(tidyverse)
library(boot)
library(mgcv)

# Function to get BCa confidence interval
get_bca_bootci = function(values){
  meanFunc <- function(x,i){mean(x[i])}
  bootci <- boot.ci(boot(data=values$observed,statistic=meanFunc,R=199))
  return( setNames(data.frame(t(bootci$bca[4:5])), c("lcl", "ucl")) )
}

# Simulate some data
set.seed(1234)
simdata = expand_grid(x=1:100, y=1:100, repeatno=1:49) %>%
  mutate(latent = sqrt((x-50)^2+(y-50)^2),
         noise = rnorm(n=n(), mean=0, sd=10),
         observed = latent + noise)

# 7 by 7 grid plot
simdata %>%
  ggplot(aes(x=x,y=y,fill=observed)) +
  theme_bw(base_size=9) +
  theme(strip.text = element_text(margin = margin(0,0,0,0, "cm"))) +
  geom_tile() +
  facet_wrap(~repeatno)

# Get mean and bootstrap CIs
mean_and_boot = simdata %>%
  group_by(x,y) %>%
  summarize(mean=mean(observed)) %>%
  left_join( simdata %>%
               group_by(x,y) %>%
               do(get_bca_bootci(.)) %>%
               unnest(lcl, ucl),
             by=c("x","y"))

# Plot the means
mean_and_boot %>%
  ungroup() %>%
  ggplot(aes(x=x,y=y,fill=mean,col=mean)) +
  theme_bw(base_size=18) +
  geom_tile()

# Plot across y=50 cross-section with CI
mean_and_boot %>%
  ungroup() %>%
  filter(y==50) %>%
  ggplot(aes(x=x, y=mean, ymin=lcl, ymax=ucl)) +
  geom_ribbon(alpha=0.4) +
  theme_bw(base_size=18) +
  geom_line() +
  ylab("Mean at y=50")

# Fit GAM model
# Maybe 100 knots is a bit much and with even more data may make the fitting take too long
gamfit = gam(observed ~ s(repeatno, bs='re') + s(x, y, bs='gp', knots=100, m=2), 
             data=simdata,
             method = 'REML')

# Obtain model predictions
modelpreds = as_tibble(t(predict(gamfit, 
               newdata = expand_grid(x=1:100, y=1:100), 
               exclude="s(repeatno)", newdata.guaranteed=T,
               se.fit=T))) %>%
  unnest(c(fit, se.fit)) %>%
  mutate(x=expand_grid(mergevar=1, x=1:100, y=1:100)$x,
     y=expand_grid(mergevar=1, x=1:100, y=1:100)$y)

# Plot model predictions
modelpreds %>%
  ggplot(aes(x=x,y=y,fill=fit)) +
  theme_bw(base_size=18) +
  geom_tile()

# Plot one cross-section (same one as above) with CI
modelpreds %>%
  filter(y==50) %>%
  ggplot(aes(x=x, y=fit, ymin=fit-se.fit*qnorm(0.975), ymax=fit+se.fit*qnorm(0.975))) +
  geom_ribbon(alpha=0.4) +
  theme_bw(base_size=18) +
  geom_line() +
  ylab("Model fit for GAM at y=50")

