Summarize features by polygon/cell ID

Aggregates an sf point dataset into one row per cell. By default computes counts and means, but the aggregation function is configurable.

Usage

summarize_by_cell(
  assigned_points_sf,
  response_var = NULL,
  predictor_vars = NULL,
  id_col = "poly_id",
  agg_funs = list(mean = function(x) mean(x, na.rm = TRUE)),
  cells_sf = NULL,
  deff = 1,
  sac = NULL,
  deff_max_n = 500L,
  quiet = TRUE,
  conf_level = NULL,
  area = FALSE
)

Arguments

assigned_points_sf

An sf object with a cell identifier column.

response_var

Optional response column name for per-cell aggregation: a single character string. A name that is not a column, or a column that is not numeric, is skipped with a warning.

predictor_vars

Optional predictor column names for per-cell aggregation. Names that are not columns, and columns that are not numeric, are skipped with a warning.

id_col

Preferred name of the polygon/cell ID column.

agg_funs

Named list of aggregation functions. Default list(mean = \(x) mean(x, na.rm = TRUE)). Additional common options: median, sum, sd. A single function (agg_funs = median) or a character vector of function names (c("median", "sum")) is also accepted. A single function is named after the expression passed: a name gives that name (median, or f for a variable f holding a function), stats::median gives median, and any other expression gives agg1; pass a named list to choose the name. Anything else falls back to the default mean with a warning.

cells_sf

Optional polygon sf layer to join cell geometries onto the output. When supplied, the return value is an sf object with the polygon geometry from cells_sf, with one row per cell in cells_sf. Cells that no feature fell in are kept, with NA summaries. Duplicate ID values in cells_sf would multiply those rows, so they are reported with a warning. Its ID column is the first of id_col, "poly_id", "polygon_id", "id", "cell_id" and "grid_id" it carries, the list and order assign_features_to_polygons() reads the polygons' IDs from. A summarised ID that matches no cell is reported with a warning, since the join drops it with its points; a cells_sf with none of those columns, or one that is not an sf object, gives a warning and a plain data frame (an error with area = TRUE). When NULL (default), a plain data.frame/tibble is returned (previous behaviour).

deff

Design-effect adjustment for standard errors. One of:

1 (default): No adjustment; the classic IID standard error sd / sqrt(n), which assumes the observations within a cell are independent. No "deff_applied" attribute is attached, so a result carrying none was computed this way.
"variogram": Compute a per-cell design effect from a fitted variogram: for n points in a cell with correlation matrix R, the effective sample size of the mean is n^2 / sum(R), so deff = sum(R) / n. This generalises Kish (substituting a constant off-diagonal correlation recovers 1 + (n - 1) * rho exactly) but lets correlation decay with distance, which matters increasingly as cells get larger and Kish's single-rho assumption degrades. Supply the fit via sac, or it is estimated when response_var is given and 'gstat' is available. Exponential, spherical and Gaussian models are supported, with a nugget and with several structured components (each weighted by its partial sill) and with gstat's 2-D geometric anisotropy (vgm(..., anis = c(angle, ratio))), applied as gstat applies it. A model of any other family falls back to deff = 1 with a warning naming it, and so does a request with no usable model (none supplied and none could be estimated, or a rejected fit that could not be replaced by an estimate); see "Value" for how to detect a fallback. A model with no structured component (a pure nugget) implies that distinct observations are uncorrelated, so it is applied as a design effect of 1 in every cell, not treated as a fallback. Points with empty geometry count towards their cells' values but not towards the correlation, with a warning.
"kish": Estimate per-variable-type intra-class correlations (ICCs) from the grouped data using a one-way random-effects ANOVA decomposition (one ICC for the response variable and a separate ICC for the predictor variables), then apply Kish's formula per cell: deff_i = 1 + (n_i - 1) * rho. The ICC is the ANOVA (method-of-moments) estimator with Donner's n0 for unequal cell sizes, not the REML estimate a mixed model returns: on a single unbalanced draw the two can differ by 0.1–0.2 (one check with cell sizes 3 to 77 and a true ICC of 0.5 gave 0.33 against REML's 0.51), while balanced designs agree to about 0.01. Neither is wrong, so do not read the difference as a defect. When multiple columns are pooled for a single ICC estimate (e.g. several predictor variables), each column is z-scored before pooling so that variables with different scales contribute equally to the variance decomposition. The response-specific ICC is used for response SEs and the predictor-specific ICC for predictor SEs. The response's ICC is never applied to predictor columns or vice versa, and it is the response's ICC (not the predictors') that sets cell_weight whenever a response was given. Requires at least 2 cells with 2+ observations and at least 2 residual degrees of freedom (N - k >= 2); the ICC is taken as 0 (no correction for that variable type) otherwise, and likewise when the estimate itself comes out at or below 0. The "deff_applied" attribute is attached when either ICC is positive, so when both are 0 there is none.
A positive number: Applied as a uniform design effect to every cell, as sd * sqrt(deff / n), exactly sqrt(deff) times the naive SE. Use when you have an external estimate of the design effect. Anything that is not a single finite number >= 1 (including a value below 1, which would shrink the standard errors, and NA or Inf) is refused with a warning and replaced by 1.

sac

Optional sac_range object from estimate_sac_range(), used when deff = "variogram". Supplying one avoids re-fitting the variogram and lets you inspect the fit the design effect is based on. A sac_range whose fit was rejected (its value is NA and it carries a rejected_reason attribute) carries no usable correlation function, so it is set aside rather than correcting by a shape that was not trusted enough to report a range: the variogram is then estimated as if no sac had been given, with a plain warning saying so, or, where that is not possible, deff falls back to 1 with the fallback warning, which names the rejection. A sac with no variogram_model attribute – a plain number or a units object, say – is a range without a correlation function, and is set aside the same way. A sac fitted to residuals (attr(sac, "detrended") TRUE) is used as given, with a warning when it corrects response columns (see "Design effects and variable types").

deff_max_n

Cells with more than this many points are subsampled before forming the n x n correlation matrix used by deff = "variogram". Default 500. It must be a single number of at least 2 when deff = "variogram"; anything else is an error.

quiet

Logical; suppress this function's progress message()s. It does not silence R warnings, nor the package's console log echo (see spatialkit_quiet for that). Default TRUE, unlike the tessellation functions, whose default is FALSE.

conf_level

Optional confidence level in (0, 1), such as 0.95. When given, every numeric response and predictor column also gets ..neff_*, ..df_*, ..ci_lo_* and ..ci_hi_* (see "Confidence intervals"). Default NULL: no interval columns, and the frame is exactly what it was before this argument existed.

area

Logical, default FALSE. With TRUE, and cells_sf supplied, the result gains cell_area (each cell's planar area in the squared units of cells_sf's CRS) and n_per_area (the count of rows in the cell over that area, a point density; a rate of anything else is that thing's agg_funs sum over cell_area). The request is refused with an error, not answered with a number, when the cells' CRS distorts areas across them by more than 1 percent, measured as the spread of planar-to-geodesic area ratios over the cells: a density is a comparison between cells and is meaningless where the map scale differs from one cell to the next. Inside a UTM zone the spread is under 0.3 percent and the request goes through; the conterminous United States forced into one zone (14 percent), or a few degrees of latitude in Web Mercator (4 percent at 48N), does not. ensure_projected() with purpose = "area" chooses an equal-area CRS for lon/lat input; build the cells in it. Lon/lat cells are measured geodesically instead: their cell_area is sf::st_area()'s area in square metres (on the sphere with s2, sf's default), not a planar area in squared degrees, and the distortion check passes by construction; with s2 switched off, sf needs the lwgeom package for that area and the request is refused without it. The measured spread is attached as attr(, "area_error").

Value

A tibble/data.frame (or sf if cells_sf given) with per-cell summaries: the ID column, n (rows in the cell), one column per agg_funs entry per variable, ..sd_* / ..se_* for every numeric response and predictor, ..neff_* / ..df_* / ..ci_lo_* / ..ci_hi_* for the same columns when conf_level is given, cell_weight, deff_applied when a design effect was requested (any deff other than 1), and cell_area / n_per_area when area = TRUE. An input column also called n is not allowed to shadow the count.

deff_applied is TRUE on every row when the requested correction was applied (exactly when the "deff_applied" attribute below is attached; under "kish", when any standard-error column was corrected) and FALSE when it fell back to the uncorrected standard errors: a refused deff, a "variogram" request with no usable model, or "kish" ICCs of 0 for every variable type (and NA on a cells_sf row no point fell in). Unlike the attribute it survives rbind() and dplyr::bind_rows() of many results. A fallback is also signalled by a warning of class "spatialkit_deff_fallback" (a Kish ICC of 0 is reported on attr(, "icc") instead), which tryCatch(spatialkit_deff_fallback = ) catches without matching the message; it is raised only when the standard errors really are the uncorrected ones.

When a correction was actually applied, an attribute "deff_applied" is attached recording it: method plus icc_resp/icc_pred and deff for "kish" (deff is the primary variable's per-cell design effect, all 1 when only the predictor ICC was positive), deff/deff_rows/rbar/crs/max_n for "variogram" (deff is the design effect at the primary variable's non-missing count per cell, deff_rows at the cell's row count, which is the vector the log line summarises as a median and a max), and deff alone for a fixed number. When cells_sf is supplied, every per-cell vector in that attribute (deff, deff_rows and rbar alike) is realigned to the joined row order, so deff[i] and rbar[i] still describe row i; cells with no observations carry NA. No attribute is attached when no correction was applied: deff = 1, a deff = "kish" request whose ICCs are all 0, or a "variogram" request that could not be fitted. A deff = "kish" request always records the ICCs it estimated on an attribute "icc" (resp and pred, NA for a variable type with no numeric column), whether or not they were positive enough to apply, so a result with no "deff_applied" still says what the ICC came out as.

The ID column keeps its input type when cells_sf's ID column and the summarised IDs already have the same class. When the classes differ, both are coerced to character in order to join (logged as a warning), and the returned ID column is therefore character. Whole numbers are written out in full for that ("100000", never "1e+05"), so an integer and a double ID of the same cell still match.

Details

This is the third step of the package's pipeline, taking the labelled layer from assign_features_to_polygons() down to cell level. What distinguishes it from a plain dplyr::group_by() + summarise() is that it carries the uncertainty of each aggregate with it: alongside every mean it returns a within-cell standard deviation, a standard error and an observation count, and it can correct that standard error for within-cell spatial autocorrelation via deff. Reach for it whenever the cell-level values will be modelled or mapped, because a cell mean over 2 observations and one over 200 are not the same measurement and nothing downstream can tell them apart otherwise.

In addition to user-specified aggregation functions, this function always computes within-cell standard deviation (..sd_<var>) and standard error (..se_<var>) for every numeric response/predictor column, plus an n column (rows falling in the cell) and a cell_weight column. These columns let downstream models account for the fact that a cell with 2 observations carries more aggregation uncertainty than one with 200.

cell_weight is the effective sample size of the primary variable: the response when one was supplied, otherwise the first predictor. It counts that variable's non-missing rows, not all rows (a cell of 10 rows with 3 finite responses carries 3 observations' worth of information about the response, not 10), and it is divided by that cell's design effect when deff applied one. With deff = 1 and no missing values it equals n. Pass it as the weights argument of a downstream regression.

Spatial autocorrelation and standard-error bias

By default (deff = 1), the ..se_* columns are computed as sd / sqrt(n), which treats the observations within each cell as independent. When data are spatially autocorrelated (the common case for the spatial workflows this package supports), within-cell observations are typically positively correlated: they share the cell's departure from the population mean, so as an estimate of the population (grand) mean a cell mean has an effective sample size smaller than n. For that estimand the naive SE is anticonservative (too small), and a downstream weighted regression that uses cell_weight or the ..se_* columns for population-level inference will produce overconfident standard errors for cells with strong intra-cell correlation. For the cell's own mean the naive SE is the right one when the cell's points are spread through it, and the corrected SE is too wide; see "What the standard error estimates" before setting deff.

Setting deff = "kish" applies an approximate correction using Kish's design effect. Separate intra-class correlations (ICCs) are estimated for response and predictor variables via a one-way random-effects decomposition across all cells. Each variable type's ICC is used for its own SE adjustment, and each cell's effective sample size is reduced to n_i / (1 + (n_i - 1) * rho). This is a first-order correction that does not require a full spatial covariance model but does require enough cells and observations for a stable ICC estimate.

A design effect estimated from the data ("kish" or "variogram") comes with a second, small-sample correction. The within-cell correlation that inflates the variance of the mean to sigma^2 * deff / n also biases the within-cell sample variance downward: under exchangeable correlation rho (Kish's own assumption), with deff = 1 + (n - 1) * rho, E[s^2] = sigma^2 * (n - deff) / (n - 1), so s^2 understates sigma^2 by very nearly the factor by which deff inflates the mean's variance, and the two errors compound rather than cancel. The standard error is therefore s * sqrt(deff / n) * sqrt((n - 1) / (n - deff)), and NA where deff >= n (the cell then holds one observation's worth of information and s carries none about sigma). This is the package's own derivation, not taken from a reference. Measured 95% interval coverage at n = 30 over 20,000 replicates: 0.921, 0.844 and 0.628 at rho = 0.2, 0.5 and 0.8 with s * sqrt(deff / n) alone, against 0.948, 0.950 and 0.949 with the rescaling.

You may also pass a fixed numeric design effect (e.g. deff = 2) to uniformly inflate standard errors: an externally supplied constant is applied as sd * sqrt(deff / n), exactly sqrt(deff) times the naive SE in every cell. (The E[s^2] correction that the estimated design effects also apply is derived from within-cell correlation and would not be justified for a number the caller chose.)

Even with the Kish correction, the adjusted SE is an approximation. For rigorous inference under spatial dependence, consider fitting an explicit spatial covariance model (e.g. via fit_bayesian_spatial_model).

What the standard error estimates

The ..se_* columns are the standard error of the cell mean as an estimate of the population (grand) mean: the unconditional quantity, in which the cell's own realised deviation is part of the error. That is the right quantity when cells are treated as samples from a common population, and the design-effect correction is calibrated for it: measured 95% interval coverage of the grand mean is 0.95 with deff = "kish" (and 0.29 with the naive SE) on exchangeable within-cell correlation, and 0.93 with deff = "variogram" on a simulated Gaussian field.

It is not the standard error of the cell's own mean (the block average over that cell), which is what a cell-level map or a regression on cell values usually wants. For that quantity the naive sd / sqrt(n) is the better of the two on offer: measured coverage 0.95 under exchangeable within-cell correlation, against very nearly 1.00 for the design-effect-corrected SE, which is too wide by the factor sqrt(deff / (1 - rho)) (4.6 at 20 points a cell and rho = 0.5). That 0.95 is exact under the exchangeable model and holds under a spatial covariance model only when the cell's points are spread through the cell; with clustered sampling inside a cell it is anticonservative for the block average too (measured 0.58), and the right quantity there is a block-kriging variance, which this function does not compute. Use deff when the cell means feed a population-level inference; leave it at 1 when they are measurements of the cells themselves and the sampling within cells is reasonably uniform.

Design effects and variable types

deff = "kish" estimates a separate ICC for response and predictor variables and applies each to its own columns. deff = "variogram" fits or accepts one correlation function and applies it to every numeric column, because a variogram is a property of the field being modelled rather than of a variable type; a predictor whose spatial structure differs markedly from the response's will have its SE corrected by the response's correlation. The internally estimated variogram is fitted to the response itself, never to OLS residuals, whatever predictor_vars holds: the ..se_resp_* columns estimate the grand mean of the response, so the correlation to correct for is the response's own. (A residual variogram, whose correlation is that of the part the predictors do not explain, is weaker; using it here dropped grand-mean coverage from 0.93 to 0.51 the moment a predictor was listed.) Pass sac explicitly when you want a different variogram, such as a residual one from estimate_sac_range(..., predictor_vars = ). A sac whose attr(sac, "detrended") is TRUE is used as given, but when it corrects response columns a warning says that their standard errors are understated, as kriging_adequacy() warns about the same mismatch.

Confidence intervals

With conf_level set, every numeric response and predictor column gains four more columns: ..neff_*, that column's effective sample size in the cell (its non-missing count over its design effect, the per-column version of cell_weight); ..df_*, the degrees of freedom the interval uses; and ..ci_lo_* / ..ci_hi_*, a t interval for the cell mean as an estimate of the grand mean. The estimand is the same as ..se_*'s, so everything in "What the standard error estimates" applies to it, including that it is not an interval for the cell's own block average. The interval is mean +/- qt((1 + conf_level) / 2, df) * se, centred on the plain mean of the column's non-missing values whatever agg_funs computes, and is NA wherever the standard error is (a single observation; complete redundancy under deff). ..neff_* and ..df_* are NA where the column has one non-missing value or none, like the standard error; cell_weight still counts such a cell (1, or 1 / deff for a numeric deff), so the two differ there.

The degrees of freedom are not neff - 1. The interval's spread comes from the within-cell sample variance, and under exchangeable correlation (deff = "kish") that variance keeps its n - 1 degrees of freedom whatever the design effect. The design effect inflates the mean's variance and biases s^2, both of which the standard error already corrects, and the resulting pivot is exactly t on n - 1 df. Measured 95% coverage of the grand mean on the Kish path, 20 cells of 20 at an ICC of 0.2 / 0.6 / 0.9: 0.954 / 0.953 / 0.952 with n - 1, against 0.992 / 1.000 / 1.000 with neff - 1, which is not an interval so much as a refusal to say anything. The effective-sample-size degrees of freedom of Faes et al. (2009) belong to a mean estimated across many correlated units whose variance is estimated from the spread between them; they do not transfer to a single cell's mean with a variance estimated from inside it. ..df_* is therefore n - 1 at deff = 1, for a numeric deff and for "kish". Under deff = "variogram" correlation decays with distance and the within-cell variance loses degrees of freedom to it; ..df_* is then the Satterthwaite (1946) moment-matched df of s^2 under the fitted correlation matrix, a fraction of n - 1 that shrinks as the range grows. Measured on a Gaussian field with an exponential range of 150 and 400 on a 1000-unit domain, 16 cells of about 25 points: 0.960 and 0.958 with that df, against 0.931 and 0.918 with n - 1, and 0.34 and 0.19 for the naive deff = 1 interval. The interval is only as good as the design effect under it: a mis-specified variogram, or an ICC estimated from too few cells, moves the coverage with it.

References

Faes, C., Molenberghs, G., Aerts, M., Verbeke, G. and Kenward, M. G. (2009). The effective sample size and an alternative small-sample degrees-of-freedom method. The American Statistician, 63(4), 389–399. doi:10.1198/tast.2009.08196

Satterthwaite, F. E. (1946). An approximate distribution of estimates of variance components. Biometrics Bulletin, 2(6), 110–114. doi:10.2307/3002019

Examples

library(sf)
set.seed(1)
n <- 200
east  <- runif(n, 0, 100)
north <- runif(n, 0, 100)
# A response with spatial structure, so the within-cell ICC is not zero.
pts <- st_as_sf(
  data.frame(x = 5e5 + east, y = 5e6 + north,
             val = 5 + 0.02 * east + 0.02 * north + rnorm(n, sd = 0.5)),
  coords = c("x", "y"), crs = 32632
)
bnd <- st_sf(geometry = st_sfc(st_polygon(list(rbind(
  c(5e5, 5e6), c(5e5 + 100, 5e6), c(5e5 + 100, 5e6 + 100),
  c(5e5, 5e6 + 100), c(5e5, 5e6)
))), crs = 32632))
grid <- create_grid_polygons(bnd, target_cells = 9, type = "square")
assigned <- assign_features_to_polygons(pts, grid)

# IID standard errors (default) vs Kish design-effect adjustment.  The
# correction is large here, and meant to be: about 22 points per cell that
# all share the cell's part of the trend carry far fewer than 22
# independent pieces of information about the cell mean.
naive <- summarize_by_cell(assigned, response_var = "val")
kish  <- summarize_by_cell(assigned, response_var = "val", deff = "kish")
data.frame(n = naive$n,
           se_naive = naive$..se_resp_val,
           se_kish  = kish$..se_resp_val)
#>    n   se_naive   se_kish
#> 1 25 0.09706775 0.7210227
#> 2 28 0.11815220 0.9279034
#> 3 25 0.11573069 0.8596517
#> 4 21 0.11556510 0.7881134
#> 5 23 0.13192931 0.9407000
#> 6 21 0.12707224 0.8665880
#> 7 10 0.18215988 0.8673238
#> 8 27 0.09534639 0.7355266
#> 9 20 0.14402822 0.9590654
attr(kish, "deff_applied")   # method, icc_resp, icc_pred, per-cell deff
#> $method
#> [1] "kish"
#> 
#> $icc_resp
#> [1] 0.6842466
#> 
#> $icc_pred
#> [1] NA
#> 
#> $deff
#> [1] 17.42192 19.47466 17.42192 14.68493 16.05343 14.68493  7.15822 18.79041
#> [9] 14.00069
#> 

# A 95% interval for each cell mean as an estimate of the grand mean, on
# the Kish-corrected standard error and n - 1 degrees of freedom.
ci <- summarize_by_cell(assigned, response_var = "val", deff = "kish",
                        conf_level = 0.95)
ci[, c("poly_id", "n", "resp_mean_val", "..neff_resp_val", "..df_resp_val",
       "..ci_lo_resp_val", "..ci_hi_resp_val")]
#> # A tibble: 9 × 7
#>   poly_id     n resp_mean_val ..neff_resp_val ..df_resp_val ..ci_lo_resp_val
#>     <int> <int>         <dbl>           <dbl>         <dbl>            <dbl>
#> 1       1    25          5.63            1.43            24             4.14
#> 2       2    28          6.30            1.44            27             4.40
#> 3       3    25          6.93            1.43            24             5.16
#> 4       4    21          6.42            1.43            20             4.77
#> 5       5    23          7.11            1.43            22             5.15
#> 6       6    21          7.89            1.43            20             6.09
#> 7       7    10          7.16            1.40             9             5.19
#> 8       8    27          7.77            1.44            26             6.25
#> 9       9    20          8.16            1.43            19             6.15
#> # ℹ 1 more variable: ..ci_hi_resp_val <dbl>