Introduction to SSIMmap

Example data and provenance

Toronto_SSIM contains 279 Aggregate Dissemination Areas (ADAs) in Toronto, using 2021 Census geography. Its columns include ADAUID, the CIMD residential instability score (CIMD) and the social dimension of the Pampalon index (PP), and a derived commuting measure (Commute). The numeric dataset is unchanged from version 0.4.0. The former description as 2016 Census Tracts is incorrect. Commute is the unweighted mean of all stored DA percentages for all intersecting DAs, including boundary-only intersections and stored zeros. The nine zero-population DAs have stored zero percentages, which are included in the means. For positive-population DAs, the percentage is 100 times the Census count for commuting to a different CSD and CD within the province or territory of residence (characteristic 2601), divided by total DA Census population (characteristic 1), rounded half-up to two decimals. It is not a household or worker share, a pooled ADA population rate, or the share of all commuting outside the CSD. Sixteen source numerator cells are officially suppressed but stored as zero in the preserved source; their values cannot be independently verified from the public release. See the installed extdata/Toronto_SSIM_provenance.txt note.

The manuscript FWI examples use packed terra rasters for British Columbia at 2 km resolution:

• fwi_0816_bc: 16 August 2023.
• fwi_0818_bc: 18 August 2023.
• fwi_1111_bc: 11 November 2023, added in version 0.4.1 from the actual source observation.

The legacy fwi_1101_bc remains the unchanged 1 November 2023 observation for compatibility. It is not an alias for November 11. The new raster was cropped and masked to the BC boundary without resampling or transformation; all three manuscript input grids match. Source checksums and preparation details are in system.file("extdata", "fwi_1111_bc_provenance.txt", package = "SSIMmap").

data("Toronto_SSIM")
data("fwi_0816_bc")
data("fwi_0818_bc")
data("fwi_1111_bc")
aug16 <- terra::unwrap(fwi_0816_bc)
aug18 <- terra::unwrap(fwi_0818_bc)
nov11 <- terra::unwrap(fwi_1111_bc)
nrow(Toronto_SSIM)
## [1] 279
terra::compareGeom(aug16, aug18, nov11)
## [1] TRUE

Choose the comparison scale explicitly

The default transformation is "none" in all three functions from version 0.4.2. Alternatives are "normal_score", "percentile", and "minmax". The manuscript explicitly uses percentiles for the Toronto variables and "none" for FWI. A high SIM or SIV after percentile transformation describes local agreement in ranks, not agreement in original variable units. Inverse-normal scores can be negative; the polygon function’s component-range checks can reject some signed-input comparisons.

The dimensionless defaults are k1 = 0.01 and k2 = 0.03. The functions compute stabilizing constants after transformation. Raster calculations use the pooled full range on common non-missing support; polygon calculations use the pooled 99th-minus-1st percentile range, with the documented fallback behavior in the source.

Polygon bandwidth selection

ssim_bandwidth() evaluates integer neighbor counts from 12 to the user-supplied max_bandwidth. Bias and variance curves are separately rescaled to [0,1] across that candidate range. Each variable’s preferred count minimizes the absolute difference between those two curves. The midpoint option rounds the mean of the two preferred counts upward.

The manuscript uses a maximum of 200. The two midpoint bandwidths are 87 and 91. The square-root alternative is ceiling(sqrt(nrow(Toronto_SSIM))), or 17.

S1 <- ssim_bandwidth(Toronto_SSIM, "CIMD", "PP",
                     max_bandwidth = 200,
                     transform = "percentile", option = "midpoint")
S2 <- ssim_bandwidth(Toronto_SSIM, "CIMD", "Commute",
                     max_bandwidth = 200,
                     transform = "percentile", option = "midpoint")
c(deprivation = S1$bandwidth, commuting = S2$bandwidth,
  square_root = ceiling(sqrt(nrow(Toronto_SSIM))))
## deprivation   commuting square_root 
##          87          91          17
p1 <- S1$plot + labs(tag = "a)") +
  theme(plot.tag = element_text(face = "plain"))
p2 <- S2$plot + labs(tag = "b)") +
  theme(plot.tag = element_text(face = "plain"))
(p1 + p2) + plot_layout(guides = "collect")

Polygon SSIM

The descriptive local result is an ordinary sf object containing SSIM, SIM, SIV, and SIP.

local1 <- ssim_polygon(Toronto_SSIM, "CIMD", "PP", global = FALSE,
                       bandwidth = S1$bandwidth,
                       transform = "percentile", do_test = FALSE)
ggplot(local1) + geom_sf(aes(fill = SSIM), colour = NA) +
  scale_fill_distiller(palette = "PRGn", limits = c(-1, 1)) + theme_void()

For polygon inference, both variables are independently permuted over polygons. From version 0.4.3, the global means and local values of all four statistics use a two-sided test centered on the corresponding permutation mean. Global p-values are unadjusted. With fdr = TRUE, local p-values are BH-adjusted across polygons separately for each statistic. Results have _p, _q, and _sig suffixes, and significance uses q <= alpha, matching rasters. With fdr = FALSE, polygon _q columns equal the unadjusted p-values and significance uses p <= alpha. The existing p_value, q_value, and sig columns remain aliases for SSIM_p, SSIM_q, and SSIM_sig.

The longer inference examples are not executed during vignette building. Running them calculates new seeded results; it does not load significance results from the manuscript.

global1 <- ssim_polygon(Toronto_SSIM, "CIMD", "PP", global = TRUE,
                        bandwidth = S1$bandwidth, transform = "percentile",
                        do_test = TRUE, R = 1000, seed = 1)
global1$p_global
local1_test <- ssim_polygon(Toronto_SSIM, "CIMD", "PP", global = FALSE,
                            bandwidth = S1$bandwidth, transform = "percentile",
                            do_test = TRUE, R = 1000, fdr = TRUE, seed = 1)
head(local1_test[, c("SSIM", "SSIM_p", "SSIM_q", "SSIM_sig",
                    "SIM_p", "SIV_p", "SIP_p")])

Original-unit FWI comparisons

w is the half-window size: w = 1, 2, and 3 give 3x3, 5x5, and 7x7 windows. The first step restricts the inputs to their common non-missing cells. Windows containing missing cells are excluded, so valid spatial support can change with window size or the date pair.

summer <- ssim_raster(aug16, aug18, w = 1, transform = "none")
fall <- ssim_raster(aug16, nov11, w = 1, transform = "none")
summer$summary
##   metric      mean        min       max         sd
## 1   SSIM 0.7686929 -0.1976864 0.9999991 0.29695831
## 2    SIM 0.7888708  0.0034073 1.0000000 0.29035795
## 3    SIV 0.9603730  0.1852894 1.0000000 0.08892729
## 4    SIP 0.9915390 -0.2499358 1.0000000 0.02908319
fall$summary
##   metric      mean           min max         sd
## 1   SSIM 0.2259786 -0.2139426534   1 0.33614608
## 2    SIM 0.2339652  0.0001376689   1 0.33838033
## 3    SIV 0.7107652  0.0308876725   1 0.28760524
## 4    SIP 0.9940730 -0.2700576644   1 0.03425971

These original-unit 3x3 comparisons have mean SSIM approximately 0.7687 for August 16 versus August 18 and 0.2260 for August 16 versus November 11. Using November 1, omitting transform = "none", or changing the window will produce different results.

summer_local <- ssim_raster(aug16, aug18, global = FALSE,
                            w = 1, transform = "none")
fall_local <- ssim_raster(aug16, nov11, global = FALSE,
                          w = 1, transform = "none")
par(mfrow = c(1, 2))
plot(summer_local[["SSIM"]], range = c(-1, 1), main = "16 vs 18 August")
plot(fall_local[["SSIM"]], range = c(-1, 1), main = "16 August vs 11 November")

Raster inference holds the first map fixed and permutes the second map over common non-missing cells. From version 0.4.2, both global and local tests use the two-sided p-value (1 + sum(abs(null - mean(null)) >= abs(observed - mean(null)))) / (R + 1), matching the permutation-mean-centered definition used for polygons. Departures above or below the null center can be significant; significance does not by itself mean high similarity. Local centers are estimated separately for each cell and statistic. Local inference uses two passes through the same permutations to keep memory use independent of the number of permutations, while retaining the random-number state of one pass. With fdr = TRUE, the four global p-values are BH-adjusted together; local p-values are adjusted across valid cells separately for each statistic.

Local significance requires both do_test = TRUE and local_test = TRUE. The returned layers use suffixes: SSIM_p, SIM_p, SIV_p, SIP_p, their _q counterparts when FDR adjustment is requested, and _sig flags. Raster significance uses q <= alpha when fdr = TRUE.

fall_test <- ssim_raster(aug16, nov11, global = FALSE, w = 1,
                         transform = "none", do_test = TRUE,
                         local_test = TRUE, R = 1000, fdr = TRUE,
                         seed = 202610281)
plot(fall_test[[c("SSIM_p", "SIM_p", "SIV_p", "SIP_p")]])

The analytical inputs are installed as package datasets. The reproduction directory contains fixed seeds, scripts and reference results for Tables 1 and 2, the Figure 6 bandwidth selection, and the map values in Figures 3–7. Full reproduction includes 1,000 permutations for each comparison and verifies local significance against the fixed-seed Figure 7 results (279 and 124 ADAs with q < 0.01). It can take substantial time and is not run during vignette building. Use a new writable output directory for each run:

source(system.file("reproduction", "reproduce_manuscript.R", package="SSIMmap"))
reproduce_manuscript("ssim-observed", full=FALSE)
reproduce_manuscript("ssim-full", full=TRUE)

See reproduction/README.txt for parallel case execution and the distinction between reproducing the analysis from the packaged inputs and reconstructing every original external Toronto source record. The portable map-rendering helper is in system.file("figures", package="SSIMmap"). Static workflow illustrations (Figures 1 and 2) and their export helper are in system.file("workflows", package="SSIMmap"); they are explanatory diagrams, not computed analysis results.

list.files(system.file("reproduction", package = "SSIMmap"))
## [1] "analysis_manifest.json" "expected"               "historical_0.4.1"      
## [4] "README.txt"             "reproduce_fwi_table1.R" "reproduce_manuscript.R"
## [7] "Table1_raw_final.csv"

References

• Wang, Z., Bovik, A. C., Sheikh, H. R., and Simoncelli, E. P. (2004). Image quality assessment: from error visibility to structural similarity. IEEE Transactions on Image Processing, 13(4), 600–612. https://doi.org/10.1109/TIP.2003.819861
• Brunsdon, C., Fotheringham, A. S., and Charlton, M. (2002). Geographically weighted summary statistics. Computers, Environment and Urban Systems, 26(6), 501–524. https://doi.org/10.1016/S0198-9715(01)00009-6
• Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate. Journal of the Royal Statistical Society: Series B, 57(1), 289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x