iiinova — where your digital self keeps evolving

ercbk.bsky.social @ercbk.bsky.social · 17/06/2026

I decided to go w/a thin plate spline for the coordinates instead. The Catboost has a kind of a nonlinear, contour-looking surface now. The local kriging model has a zonal pattern which is due to kNN methodology and the data's spatial distribution (I think). krigeST was the better scoring model. fin

$full code - https://ercbk.github.io/Data-Science-Notebook/qmd/geospatial-spat-temp.html#sec-geo-sptemp-krig-ex Only showing code for the catboost row. Both are similar. p_cat_pred <- ggplot() + stars::geom_stars( data = star_preds_cat, aes(x = x, y = y, fill = predicted) ) + geom_sf( data = sf_germany_line, color = "grey40", linewidth = 0.3 ) + geom_sf( data = sf_obs_sample, color = "white", size = 0.8, alpha = 0.7 ) + geom_point( data = sf_obs_sample_nl |> st_drop_geometry(), aes(x = x, y = y, fill = pm10), shape = 21, size = 3, color = "white", stroke = 0.1 ) + facet_wrap(~as.Date(time), nrow = 1) + scico::scale_fill_scico( palette = "glasgow", limits = legend_limits, direction = -1, name = "PM10\n(µg/m³)", na.value = "transparent" ) + coord_sf(lims_method = "geometry_bbox") + labs(title = "Catboost") + theme_charcoal_map() theme_patch <- theme( plot.title = element_text( family = "Anybody", size = 25, color = "white", face = "bold" ), plot.subtitle = element_text( family = "Anybody", size = 16, color = "white" ), plot.background = element_rect(fill = "#191a1a", color = NA) ) p_krige_pred / p_cat_pred + plot_layout(guides = "collect") + plot_annotation( title = "PM10 Predictions", subtitle = "Station location color represents observed PM10 level", theme = theme_patch )$

020

ercbk.bsky.social @ercbk.bsky.social · 17/06/2026

With a variogram model, the kriging model can be fit. I used mlr3spatiotempcv to also choose a ML model along with tuning a local kriging model. First time around I gave the ML models calendar features and the raw projected coordinates. The result is this paint brush effect on the prediction surface

A map of Germany that looks like a paint brush was dipped in orange pain and lightly run across vertically and horizontally over a light blue surface.

210

ercbk.bsky.social @ercbk.bsky.social · 17/06/2026

Here's a video of 3-D surface interactivity. I added javascript so that the subplots are synced, but it's not a smooth sync.

110

ercbk.bsky.social @ercbk.bsky.social · 17/06/2026

The last plot, expands the heat map to 3-D. Semi-variance now has it's own axis, and the surface is still colored according to the semi-variance values. I find the static version of this plot the most difficult to judge the models. Plotly's interactivity helps with this. It's a bear to code though.

source - https://cran.r-project.org/web/packages/gstat/vignettes/st.pdf
full code - https://ercbk.github.io/Data-Science-Notebook/qmd/geospatial-spat-temp.html#sec-geo-sptemp-krig-var-ex

plot(vario_emp,
best_mod_fits_all,
map = FALSE,
wireframe = TRUE,
all = TRUE)

110

ercbk.bsky.social @ercbk.bsky.social · 17/06/2026

The heatmap has time on the y-axis and distance on the x-axis with cells colored by semi-variance. The noise of the empirical variogram (sample, observed) makes this more difficult to judge the models, but I'd still have to go w/the Sum-Metric (or Simple Sum Metric). (btw lowest SSE is Sum-Metric)

$source - https://cran.r-project.org/web/packages/gstat/vignettes/st.pdf full code - https://ercbk.github.io/Data-Science-Notebook/qmd/geospatial-spat-temp.html#sec-geo-sptemp-krig-var-ex # observed noisy, metric abrupt temporally, prod/sep gradual in space and time ggplot(tib_plot, aes(x = factor(floor(avgDist)), y = timelag, fill = gamma)) + geom_tile() + scale_fill_gradientn( colors = rev(scico::scico(100, palette = "glasgow")), na.value = "white" ) + scale_x_discrete( # skips every other value to prevent overlap of the labels breaks = function(x) x[seq(1, length(x), by = 2)] ) + facet_wrap(~ model, nrow = 2) + labs( title = "Spatio-Temporal Semi-Variance Surfaces", subtitle = "Model comparison between the best fit of each type of model", x = "distance (km)", y = "time lag (days)", fill = "Semi-\nVariance" ) + theme_notebook( strip.background = element_blank(), panel.grid = element_blank(), legend.key.height = unit(1.95, "cm"), )$

110

ercbk.bsky.social @ercbk.bsky.social · 17/06/2026

The empirical variogram is calculated based on the data. Then a variogram model is chosen according to how well it approximates the empirical variogram. {gstat} provides 3 diagnostic plots. In the multi-line, the spread represents temporal autocorr and y-axis — spatial autocorr. Sum Metric is best.

110

ercbk.bsky.social @ercbk.bsky.social · 17/06/2026

Hot spot analysis extends the interpretation of local measures by classifying the cluster centers. The class depends on the statistic value and which quadrant of the Moran plot the observation lands. Interpretation exs: High - cluster of high values Low-High - low value center, high value neighbors

source - https://r-spatial.org/book/15-Measures.html#local-gearys-c_i
full code - https://ercbk.github.io/Data-Science-Notebook/qmd/association-spatial.html#sec-assoc-spat-local-ex

tmap_options_reset()
tmap_theme_notebook()

notebook_colors <- unname(swatches::read_ase(here::here("palettes/Forest Floor.ase")))

tm_shape(clust_mcg_iturn_pres15) +
tm_polygons(col = "#4A3728") +
tm_fill(
fill = "clust_type",
fill.scale = tm_scale(
values = c(notebook_colors[[5]], notebook_colors[[6]], notebook_colors[[7]]),
value.na = notebook_colors[[2]],
label.na = "Not Determined"
),
fill.legend = tm_legend(
"Hotspot Type",
frame = FALSE,
position = tm_pos_on_top(
pos.v = 0.4,
pos.h = 0.66
)
)
) +
tm_facets_wrap("method_type", ncol = 2) +
tm_title(
"High and Low Voter Turnout Clusters in Poland (2015)",
size = 2,
fontface = "bold",
# get title flush with left side to line up with subtitle
# every number except 0.05 is default
padding = c(0.25, 0.15, 0.25, 0.25)
) +
tm_title(
"Comparing Local Moran's I, Local Getis-Ord G, and Local Geary's C methods",
size = 1.25
) +
tm_layout(
panel.labels = c(
"MoranI-Permutation",
"Getis-OrdG-Permutation",
"GearyC-Permutation"
)
)

110

ercbk.bsky.social @ercbk.bsky.social · 17/06/2026

For the local Moran's I test, units with I statistics significantly below or above the expected value indicate negative or positive autocorrelation respectively. p-values should be adjusted. Using the 0.005 threshold, p-values in darker blue shades indicate potential cluster centers of voter turnout

source - https://r-spatial.org/book/15-Measures.html#local-morans-i_i
full code - https://ercbk.github.io/Data-Science-Notebook/qmd/association-spatial.html#sec-assoc-spat-local-ex

pv_brks <- c(0, 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1)

tmap_options_reset()
tmap_theme_charcoal()

tm_shape(pv_iturn_pres15) +
tm_polygons(col = "gray40") +
tm_fill(
fill = "p_value",
fill.scale = tm_scale(breaks = pv_brks, values = "scico.davos"),
fill.legend = tm_legend(
"Adjusted P-Values (BH)",
frame = FALSE,
position = tm_pos_on_top(
pos.v = 0.4,
pos.h = 0.66
)
)
) +
tm_facets_wrap("pval_type", ncol = 2) +
tm_title(
"Local Moran's I P-Values",
size = 2,
fontface = "bold"
) +
tm_title("With regards to Poland's voter turnout (2015)", size = 1.25) +
tm_layout(
panel.labels = c(
"Analytical-Conditional",
"Permutation-Std.Deviate",
"Permutation-Rank"
)
)

110

ercbk.bsky.social @ercbk.bsky.social · 17/06/2026

As mentioned before, spatial autocorrelation tests use spatial weights created from neighbor graphs. Most measures have global and local versions. One, Moran's I, has an extension called a Moran Plot. Among other things, it labels locations with strong influence on the measure.

source - https://r-spatial.org/book/15-Measures.html#local-measures
full code - https://ercbk.github.io/Data-Science-Notebook/qmd/association-spatial.html#sec-assoc-spat-global-ex

data(pol_pres15, package = "spDataLarge")
pol_pres15 <- pol_pres15 |>
janitor::clean_names()

iturnout <- pol_pres15 |>
st_drop_geometry() |>
pull(i_turnout)

lw_w_pol <- nb_poly_pol |> nb2listw(style = "W")

infl_turnout_pol <- moran.plot(
iturnout,
listw = lw_w_pol,
labels = pol_pres15$TERYT,
cex = 1,
pch = ".",
xlab = "I round turnout",
ylab = "lagged turnout"
)

120

ercbk.bsky.social @ercbk.bsky.social · 17/06/2026

This heatmap (tidyplots) and table (tinytable) summarize the previous map but for all nodes. For the first lag (first row) in the heatmap, the brightest yellow cell says around 16 census tracts have 6 First Order neighbors. Connectedness peaks around the 3rd lag. More explicitly shown in the table.

$source - https://r-spatial.github.io/spdep/articles/nb.html#higher-order-neighbours full code - https://ercbk.github.io/Data-Science-Notebook/qmd/geospatial-spatial-weights.html#sec-geo-swgt-diag plot_col <- function(d, color, ...) { ggplot(d, aes(x = tot_neigh, y = tot_nodes)) + geom_col(color = color, fill = color) + # rm axes w/o using theme elts scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(expand = c(0, 0)) + theme_void() } tt(dat_nb_wide) |> style_tt(background = "#FFFDF9FF") |> # headers style_tt(i = 0, bold = TRUE, color = "#000E54") |> # row names style_tt(j = 1, bold = TRUE, color = "#000E54") |> # only color nonzero values style_tt(i = dat_nb_wide > 0, background = "#0095A8FF", color = "white") |> # add column chart plot_tt(j = 26, fun = plot_col, data = lst_tbl_nb, color = "#0095A8FF", height = 2)$

130

ercbk.bsky.social @ercbk.bsky.social · 17/06/2026

Spatial lags of neighbors are the number of steps away from a node. The neighbors at a particular spatial lag are that lag’s order of neighbors (e.g. 4th lag has 4th order neighbors). I was having trouble visualizing exactly what what that meant, so I thought it was finally time to try gganimate.

120

ercbk.bsky.social @ercbk.bsky.social · 17/06/2026

First up is neighbor algorithms and spatial weights which represent spatial autocorrelation using relationships between nodes on a graph. Econometric regression and autocorrelation tests make use of these. This is a map of NY census tracts illustrating k-NN where all nodes have at least k neighbors.

$source - https://r-spatial.github.io/spdep/articles/nb_sf.html#k-nearest-neighbours full code - https://ercbk.github.io/Data-Science-Notebook/qmd/geospatial-spatial-weights.html#sec-geo-swgt-nbalg-distb nb2tbl <- function(nb_list, coords, color = "black") { nb_len <- length(nb_list) if (length(color) == 1) { nb_col <- rep(color, nb_len) } else { nb_col <- color } # flattens neighbor vectors into 1 vector cardnb <- card(nb_list) # source node indexes i <- rep(1:nb_len, cardnb) # neighbor node indexes j <- unlist(nb_list) nb_tib <- tibble::tibble( # segment coords x_i = st_coordinates(coords)[, 1, drop = TRUE][i], x_j = st_coordinates(coords)[, 1, drop = TRUE][j], y_i = st_coordinates(coords)[, 2, drop = TRUE][i], y_j = st_coordinates(coords)[, 2, drop = TRUE][j], # segment colors seg_col = nb_col[i] ) return(nb_tib) } ny_nb_ll <- knn2nb(knearneigh(ny_ct_ll, k = 4, longlat = TRUE)) # Convert neighbors list to tibble tbl_nb <- nb2tbl(ny_nb_ll, ny_ct_ll) ggplot() + geom_sf(data = ny_sf_ll, fill = "#FFFDF9FF") + geom_sf(data = ny_ct_ll, color = "#000E54") + geom_segment(aes(x = x_i, y = y_i, xend = x_j, yend = y_j), data = tbl_nb, color = "#000E54") + ggtitle( "Central New York", subtitle = "k-NN with k = 4" )$

130

ercbk.bsky.social @ercbk.bsky.social · 27/05/2026

Thanks for the package. Facetted plots are a nightmare to create in plotly. The interactivity and rotation (janky) are nice for eda, but I might've scrapped this one if I could've created a static version version this easily.

030

ercbk.bsky.social @ercbk.bsky.social · 08/12/2025

To make it flow from left to right, you have to interchange x, y, xend, yend and add negative signs.

library(ggdag); library(ggplot2)

notebook_colors <- unname(swatches::read_ase(here::here("palettes/Forest Floor.ase")))

da_dag <- dagify(
Y ~ A + B + C,
D ~ A + B,
C ~ D,
outcome = "Y",
exposure = "A"
) |>
tidy_dagitty(layout = "sugiyama")

da_dag |>
ggplot(aes(
x = -.data$y,
y = -.data$x,
xend = -.data$yend,
yend = -.data$xend)
) +
geom_dag_node(color = notebook_colors[[6]]) +
geom_dag_edges() +
geom_dag_text() +
theme_dag(
panel.background = element_rect(fill='#FFFDF9FF',
color = NA),
panel.border = element_blank(),
panel.grid.minor = element_blank(),
plot.background = element_rect(fill='#FFFDF9FF',
color = NA),
legend.background = element_rect(color = NA,
fill='#FFFDF9FF'),
legend.box.background = element_rect(fill='#FFFDF9FF',
color = NA)
)

020

ercbk.bsky.social @ercbk.bsky.social · 08/12/2025

Btw if anyone is looking to do this in {ggdag}, it looks like layout = "sugiyama" gets it done.

110

ercbk.bsky.social @ercbk.bsky.social · 26/09/2025

This is what I get from gtsummary::tbl_regression(nom_mod, intercept = TRUE, exponentiate = TRUE). Are you not setting intercept = TRUE?

000

ercbk.bsky.social @ercbk.bsky.social · 28/08/2025

If it's just for quick exploration, the duckdb ui provides some charts, duckdb.org/2025/03/12/d... Other UIs that can connect duckdb allow you to export charts sql-workbench.com/docs/queryin... www.timestored.com/qstudio/help...

DuckDB UI showing a table with bar and density plots in a side panel

031