I'm stuck on which statistic to use with a spatial data set to resolve population concentrations in a large area, when I have only sampled a small area relative to that large area. Here's an example of my problem: somewhere there exists some hypothetical organism, the volcanobeast, that has a distribution in space, with a true concentration somewhere between 1 and 2 organisms per meter ^2, as follows.



v2 = (volcano - min(volcano))/max(volcano) * 4

v3 <- reshape2::melt(v2) %>% rename(Xpos = Var1, Ypos = Var2)

ggplot(v3, aes(x = Xpos, y = Ypos, fill = value)) + geom_raster() + scale_fill_viridis_c()

true volcanobeast abundances

Here each square is a kilometer plot of land and the value is the number of organisms per meter squared on that plot of land. Observe that the data are autocorrelated in space. Nearby squares have similar values to eachother.

Let's say I sample one meter square within each kilometer of land. I think, that if these organisms are randomly distributed, I will get values per square meter that follow a poisson distribution (please correct me if I am wrong).

v4 <- v3 %>% mutate(meas = map_dbl(value, ~rpois(1, lambda = .)))

ggplot(v4, aes(x = Xpos, y = Ypos, fill = meas)) + geom_raster() + scale_fill_viridis_c()

observed volcanobeast abundances

So here, we have the number of volcanobeasts sampled per square meter, on each 1 km plot of land. I am wondering if, using my knowledge that the data are autocorrelated and that I am drawing from a poisson distribution, I can get an estimate and confidence interval of the true concentration of volcanobeasts in each square kilometer.

Any suggestion on what might be the appropriate statistic for this? Thanks!

  • 2
    $\begingroup$ There are certainly methods for this, more than I can mention here! Check GAMs with smoothing splines for spatial data, you could check out a 2D smooth on Xpos and Ypos as a good place to start. The mgcv package is good for this $\endgroup$ – ASeaton Mar 4 at 15:25

Per @ASeaton's answer, I tried using a GAM for this and I got a pretty satisfactory answer.


gmod <- gam(meas ~ s(Xpos, Ypos), data = v4, family = "poisson")
gpred <- predict(gmod, type = "response", se.fit = TRUE) %>% as.data.frame()
v5 <- bind_cols(v4, gpred)
fit_plt <- ggplot(v5, aes(x = Xpos, y = Ypos, fill = fit)) + geom_raster() + scale_fill_viridis_c(limits = c(0,2))

fitted data

Also, modifying directions found below, I was able to come up with upper and lower confidence intervals for this estimate.


# calculate upper and lower bounds
fam <- family(gmod)
ilink <- fam$linkinv

ndata <- predict(gmod, type = "link", se.fit = TRUE) %>% as.tibble() %>% setNames(c("fit_link", "se_link")) %>%
  mutate(fit_resp = ilink(fit_link),
         right_upr = ilink(fit_link + (3 * se_link)),
         right_lwr = ilink(fit_link - (3 * se_link))) %>%

# display upper and lower bounds

upr_plt <- ggplot(ndata, aes(x = Xpos, y = Ypos, fill = right_upr)) + geom_raster() + scale_fill_viridis_c(limits = c(0, 2.1))

lwr_plt <- ggplot(ndata, aes(x = Xpos, y = Ypos, fill = right_lwr)) + geom_raster() + scale_fill_viridis_c(limits = c(0,2.1))

plot_grid(upr_plt, lwr_plt)



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.