How to read this document

This notebook is written to be read like a textbook chapter: every code chunk is followed (or preceded) by an explanation of what the code does, line by line, why we’re doing it, and what to look for in the output. If you are new to R, do not skim the code — read the explanation, then re-read the code slowly, then run the chunk and compare what you see against what the text told you to expect. That loop (read → predict → run → check) is how you actually learn a language.

A few R survival facts used constantly below, collected here once:

%>% (the pipe) takes the thing on its left and feeds it in as the first argument of the function on its right. x %>% f(y) is exactly f(x, y). Pipes let you read a data recipe top-to-bottom instead of inside-out: data %>% filter(...) %>% group_by(...) %>% summarise(...) reads as “take the data, keep some rows, form groups, compute summaries.”
A data frame is R’s table: columns have names and types (numeric, character, factor…), rows are observations. dplyr verbs (mutate, filter, group_by, summarise, left_join, …) each take a data frame and return a new data frame — they never modify the original in place unless you re-assign with <-.
A factor is R’s type for categorical variables. It stores the set of allowed values (“levels”) alongside the data. Models care deeply about factors: glm() automatically turns a factor with 10 levels into 9 dummy (0/1) variables, measured against a “reference level” — we will exploit that when we extract relativities.
function(x) { ... } defines your own function. The last expression evaluated is the return value (no return() needed, though allowed).
<- is assignment (“gets”). x <- 5 stores 5 under the name x.
Vectorization: most R operations work on whole vectors at once. pmin(v, 1) caps every element of v at 1 in one call — no loop needed. When we do write explicit for loops below, it is deliberate (the credibility algorithm is genuinely sequential).

What we are building, in one paragraph

An insurer needs to charge different premiums in different places, but a single postal code has far too little data to be priced on its own experience. So we (1) start from raw policy data and a raw shapefile (the file format for map polygons) and connect the two; (2) stabilize each postcode’s loss experience by blending it with its geographic neighbors using actuarial credibility theory (Jennings, CAS 2008); (3) group the ~1,100 postcodes into ~10 territories with k-means clustering, using both risk and location so the groups come out geographically coherent; and (4) fit generalized linear models (GLMs) on the policy-level data with territory as one rating factor among many, so each territory’s final price relativity is measured holding driver age, bonus-malus, vehicle characteristics, etc. constant. Steps 2-3 are “Part 1”, step 4 is “Part 2”, and step 1 — the unglamorous foundation everything sits on — is “Part 0”.

Unlike freMTPL2 (which only carries 22 pre-made administrative regions), this dataset gives you ~1,100 postcodes with coordinates — so the credibility problem, the contiguity problem, and the shapefile work are all real, not toy.

Loading packages

required_pkgs <- c(
  "dplyr", "tidyr", "ggplot2", "forcats", "stringr",
  "sf", "geosphere",
  "cluster", "factoextra",
  "broom", "scales", "readr", "tibble", "knitr", "patchwork"
)
missing_pkgs <- required_pkgs[!required_pkgs %in% rownames(installed.packages())]
if (length(missing_pkgs) > 0) install.packages(missing_pkgs)
invisible(lapply(required_pkgs, library, character.only = TRUE))

theme_set(theme_minimal(base_size = 11))
relativity_fill <- function(name = "Relativity") {
  scale_fill_gradient2(midpoint = 1, low = "#2166AC", mid = "white",
                        high = "#B2182B", name = name)
}

What this chunk does, piece by piece:

required_pkgs <- c(...) builds a character vector (a list of text strings) naming every package this notebook uses. c() is “combine” — the most basic R function, gluing values into a vector.
installed.packages() returns a matrix with one row per installed package; rownames() of it is therefore the vector of installed package names. The expression !required_pkgs %in% rownames(...) asks, for each required package, “is it not installed?” (%in% tests membership; ! negates). Indexing required_pkgs[...] with that TRUE/FALSE vector keeps only the missing ones — this TRUE/FALSE (“logical”) indexing pattern is everywhere in R.
if (length(missing_pkgs) > 0) install.packages(missing_pkgs) installs only what’s missing, so the notebook is runnable on a fresh machine but doesn’t waste time reinstalling on yours.
lapply(required_pkgs, library, character.only = TRUE) calls library(pkg) for each package name. lapply = “apply this function to each element of the list.” character.only = TRUE is needed because library() normally expects a bare name (library(dplyr)), not a string (library("dplyr")). invisible() just suppresses the printed output.
theme_set(theme_minimal(...)) sets a default visual theme for every ggplot2 chart in the document, so we don’t repeat it on each plot.
The chunk header says cache=FALSE, overriding the document-wide cache=TRUE from the setup chunk — and the reason is a trap every R Markdown user hits eventually. Caching means “if this chunk’s code hasn’t changed, skip running it and reuse the saved results.” But library() calls have no saved result — their entire effect is the side effect of attaching packages to the running session. A cached packages chunk gets skipped on re-knit, no packages get attached, and any non-cached chunk downstream fails with could not find function "..." for functions that worked fine yesterday. Rule: chunks whose purpose is a side effect (loading packages, setting options, connecting to databases) must never be cached.
relativity_fill is our first user-defined function. It wraps scale_fill_gradient2(), a ggplot color scale with three anchor colors: blue below the midpoint, white at it, red above. We set midpoint = 1 because a relativity of exactly 1.00 means “priced at the base level” — so on our final maps, blue will literally mean “cheaper than base” and red “more expensive than base.” Defining it once guarantees every relativity map in the document uses identical coloring.

Package roles (so the imports aren’t mysterious): dplyr/tidyr manipulate data frames; ggplot2 draws every chart; sf handles spatial data (shapefiles, polygons, projections); geosphere computes distances on the Earth’s surface; cluster and factoextra provide clustering algorithms and diagnostics; broom turns model output into tidy data frames; patchwork lets us stack two ggplots with /; the rest are small utilities (string handling, scales, CSV writing).

Part 0 - Foundation: policy data + shapefile from scratch

0.1 The portfolio: `CASdatasets::beMTPL97`

CASdatasets is a large collection of real insurance datasets maintained by Christophe Dutang. It is not on CRAN (R’s default package repository) because it’s too big, so install.packages("CASdatasets") alone would fail — you must point R at the maintainer’s own repository via the repos argument. This is a one-time installation; the chunk is marked eval=FALSE so knitting the document doesn’t re-run it every time.

(xts and zoo are time-series packages that CASdatasets depends on; installing them explicitly first avoids a dependency error. type = "source" tells R to build the package from source code rather than look for a pre-built binary, which this repository doesn’t provide.)

if (!requireNamespace("CASdatasets", quietly = TRUE)) {
  stop("CASdatasets is not installed.")
}

data("beMTPL97", package = "CASdatasets")

mtpl_raw <- beMTPL97

dplyr::glimpse(mtpl_raw)

## Rows: 163,212
## Columns: 18
## $ id       <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18…
## $ expo     <dbl> 1.00000000, 1.00000000, 1.00000000, 1.00000000, 0.04657534, 1…
## $ claim    <fct> 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ nclaims  <int> 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ amount   <dbl> 1618.0010, 0.0000, 0.0000, 0.0000, 155.9746, 0.0000, 155.9746…
## $ average  <dbl> 1618.0010, NaN, NaN, NaN, 155.9746, NaN, 155.9746, NaN, NaN, …
## $ coverage <fct> TPL, TPL+, TPL, TPL, TPL, TPL, TPL++, TPL, TPL++, TPL+, TPL++…
## $ ageph    <int> 50, 64, 60, 77, 28, 26, 26, 58, 59, 34, 39, 55, 77, 49, 63, 4…
## $ sex      <fct> male, female, male, male, female, male, male, female, male, m…
## $ bm       <int> 5, 5, 0, 0, 9, 11, 11, 11, 0, 7, 1, 11, 9, 8, 0, 1, 0, 10, 0,…
## $ power    <int> 77, 66, 70, 57, 70, 70, 55, 47, 98, 74, 85, 87, 104, 71, 104,…
## $ agec     <int> 12, 3, 10, 15, 7, 12, 8, 14, 3, 6, 2, 6, 6, 19, 7, 2, 7, 6, 1…
## $ fuel     <fct> gasoline, gasoline, diesel, gasoline, gasoline, gasoline, gas…
## $ use      <fct> private, private, private, private, private, private, private…
## $ fleet    <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0…
## $ postcode <int> 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1…
## $ long     <dbl> 4.355223, 4.355223, 4.355223, 4.355223, 4.355223, 4.355223, 4…
## $ lat      <dbl> 50.84539, 50.84539, 50.84539, 50.84539, 50.84539, 50.84539, 5…

Line by line:

requireNamespace("CASdatasets", quietly = TRUE) returns TRUE/FALSE for “is this package available?” without attaching it. Wrapping it in if (!...) stop(...) is defensive programming: if the package is missing we fail immediately with a clear message, instead of 30 chunks later with a confusing one. You will see this fail-fast pattern throughout the notebook — professional analysis code checks its own assumptions.
data("beMTPL97", package = "CASdatasets") loads a dataset that ships inside a package into your workspace. After this line an object named beMTPL97 exists.
mtpl_raw <- beMTPL97 copies it under our own name. Convention used here: _raw = untouched input; we never overwrite the raw object, so we can always compare “before cleaning” vs. “after.”
dplyr::glimpse() prints a compact sideways summary: every column, its type, and its first few values. The dplyr:: prefix explicitly says which package the function comes from — good habit when two loaded packages export functions with the same name. Always glimpse a dataset before touching it.

Data dictionary (one row per policy):

Field	Meaning
`id`	Policy identifier
`expo`	Exposure: fraction of 1997 the policy was in force
`claim` / `nclaims`	Claim indicator / number of claims filed
`amount`	Total claim amount (EUR) across the policy’s claims
`average`	Mean claim amount = `amount / nclaims`
`coverage`	`TPL`, `TPL+` (adds limited material damage), `TPL++`
`ageph`	Policyholder age
`sex`	Policyholder sex
`bm`	Belgian bonus-malus level, 0 (best) - 22 (worst)
`power`	Vehicle power (kW)
`agec`	Vehicle age (years)
`fuel`	`gasoline` / `diesel`
`use`	`private` / `work`
`fleet`	Fleet-policy indicator
`postcode`	4-digit Belgian postal code — our building block
`long`, `lat`	Centroid of the postcode’s municipality — our geography

Three of these deserve a moment, because the whole analysis stands on them:

Exposure (expo) is the actuarial unit of “how much insurance was provided.” A policy in force all year contributes 1.0 policy-years; one that started in July contributes ~0.5. Claim counts must always be judged per unit of exposure — 3 claims on 100 policy-years is very different from 3 claims on 10.
Frequency = claims / exposure. Severity = cost per claim. Pure premium = losses / exposure = frequency × severity. The pure premium is the break-even price of the risk before expenses and profit — it is the quantity territorial rating tries to estimate by location.
postcode + long/lat make this dataset special: we know where every policy lives, at a fine enough grain (~1,100 areas) that geography becomes a genuine modeling problem instead of a lookup.

Cleaning

Real data always needs cleaning, and it should happen once, at the top, visibly — not scattered through the analysis.

mtpl <- mtpl_raw %>%
  mutate(
    postcode = as.integer(as.character(postcode)),
    expo     = pmin(expo, 1),
    fleet    = factor(fleet)
  ) %>%
  filter(expo > 0)

# A tiny number of rows have claim amounts recorded with nclaims = 0 (or
# vice versa) - inspect, then drop the inconsistent ones rather than guess.
mtpl %>% filter((nclaims > 0) != (amount > 0)) %>% nrow()

## [1] 0

mtpl <- mtpl %>% filter((nclaims > 0) == (amount > 0))

nrow(mtpl)

## [1] 163212

sum(mtpl$expo)

## [1] 145216.8

n_distinct(mtpl$postcode)

## [1] 583

What each line does and why:

mutate() adds or overwrites columns. Inside it:
- as.integer(as.character(postcode)) — a classic R idiom. If postcode arrived as a factor, converting it straight to integer with as.integer() would return its internal level codes (1, 2, 3…) rather than the actual postcodes (1000, 1030…) — a famous R trap. Going through as.character() first gives the text “1000”, which then converts to the number 1000 correctly. We want postcodes numeric because the shapefile stores them numerically and joins must match types.
- pmin(expo, 1) — “parallel minimum”: element-by-element, take the smaller of expo and 1. A policy cannot be exposed more than a full year within a one-year dataset, so values above 1 are data errors; we cap rather than drop, keeping the policy but fixing the impossible value.
- factor(fleet) — fleet is coded 0/1, but it’s a category (fleet policy: yes/no), not a quantity. Declaring it a factor tells the GLM later to estimate a separate effect for fleet policies rather than treat “1” as “one unit more fleet than 0” (same fit for binary, but factor is the honest type and prints better output).
filter(expo > 0) keeps only rows where the condition is TRUE. Zero exposure means the policy was never actually at risk — and it would crash log(expo) in the GLM offset later (log of 0 is -Inf).
The consistency check is worth internalizing as a pattern. The claim count and the claim amount come from two different administrative systems, and they can disagree: a row with nclaims = 0 but positive amount (or the reverse) is internally contradictory. The expression (nclaims > 0) != (amount > 0) is TRUE exactly for those contradictory rows. We first count them (%>% nrow()) so the anomaly is visible in the knitted output — never silently delete data — and then keep only the consistent rows with the mirrored condition ==.
The three closing lines print our post-cleaning vital signs: number of policies, total policy-years of exposure, and n_distinct(postcode) — the number of distinct building blocks we’ll be mapping.

~1,100 distinct postcodes for ~160k policies — an average of ~145 policies per postcode, but very unevenly spread. That unevenness is the whole reason Part 1 exists.

0.2 Download and read the real Belgian postcode shapefile

What is a shapefile? The near-universal file format for map geometry, invented by Esri in the 1990s. One “shapefile” is actually a bundle of files sharing a name: .shp (the polygon coordinates), .dbf (a table of attributes — one row per polygon, e.g. its postcode), .shx (an index), and .prj (the coordinate reference system — see below). They must travel together, which is why shapefiles are distributed zipped.

The one we use (npc96_region_Project1.shp, one polygon per Belgian postcode) is the exact file used in Henckaerts et al. (2018), publicly hosted in the PE-pricing-analytics workshop repo. The workflow below — download → unzip → st_read() → fix the CRS → standardize the join key — is the same one you’d use with any statistical-office shapefile (US Census ZCTAs, StatCan FSAs, …).

shp_dir <- "data/belgium_postcode_shp"
if (!dir.exists(shp_dir)) {
  dir.create(shp_dir, recursive = TRUE)
  zip_url  <- "https://github.com/katrienantonio/PE-pricing-analytics/raw/master/scripts/shape%20file%20Belgie%20postcodes/postcode.zip"
  zip_path <- file.path(shp_dir, "postcode.zip")
  download.file(zip_url, zip_path, mode = "wb")
  unzip(zip_path, exdir = shp_dir)
}

shp_file <- list.files(shp_dir, pattern = "\\.shp$", full.names = TRUE, recursive = TRUE)[1]
shp_file

## [1] "data/belgium_postcode_shp/npc96_region_Project1.shp"

belgium_shape <- st_read(shp_file, quiet = TRUE)
names(belgium_shape)

## [1] "POSTCODE"   "NAAM1"      "NAAM2"      "POSTCODELA" "Shape_Leng"
## [6] "Shape_Area" "geometry"

st_crs(belgium_shape)$Name

## [1] "BD72 / Belgian Lambert 72"

Line by line:

The whole download block sits inside if (!dir.exists(shp_dir)) — a caching pattern. First run: the folder doesn’t exist, so we create it (dir.create(..., recursive = TRUE) also creates the parent data/ folder if needed), download, and unzip. Every later run: the folder exists, the block is skipped, no network access. Cache anything you download; re-downloading on every knit is slow and fragile.
file.path() joins path pieces with the correct separator for your operating system — always use it instead of pasting “/” yourself.
download.file(..., mode = "wb"): wb = write binary. On Windows the default mode can corrupt binary files (zips, images) by translating line endings; mode = "wb" prevents that. Forgetting it is a classic Windows-specific bug.
unzip(zip_path, exdir = shp_dir) extracts the bundle into the folder.
list.files(..., pattern = "\\.shp$", full.names = TRUE, recursive = TRUE) searches the folder for the geometry file. The pattern is a regular expression: \\. is a literal dot (an unescaped . means “any character” in regex) and $ anchors the match to the end of the filename — so this matches files ending in .shp and nothing else. [1] takes the first match. Finding the file by pattern instead of hard-coding its name makes the code survive a renamed or restructured zip.
st_read() (from sf) reads the whole bundle into an sf object: an ordinary data frame — you can filter() it, mutate() it, join to it — with one special geometry column holding each row’s polygon. This “spatial data is just a data frame” design is why sf has become the standard: everything you know about data frames transfers directly.
names() shows the attribute columns that came from the .dbf; st_crs(...)$Name prints the coordinate reference system read from the .prj. Inspect both before using any shapefile — you need to know what the join key column is called and what units the coordinates are in.

Coordinate reference systems (CRS), the two-minute version. The Earth is a curved surface; a map is flat. A CRS is the rulebook for converting between the two. There are two families you’ll meet constantly:

Geographic CRS — positions as longitude/latitude degrees on the globe. The universal standard is WGS84, code EPSG:4326 (what GPS and Google Maps use). Our policy data’s long/lat are in this.
Projected CRS — positions in meters on a flat plane, using a projection tuned to a region. National mapping agencies use these; Belgium’s is Lambert 72, and that’s what this shapefile is stored in.

Every spatial layer carries its CRS, and layers must share one before you can combine or compare them. Mixing CRSs doesn’t error — it silently puts your data in the wrong place, which is worse.

# The projection on file is a Belgian national CRS (Lambert 72). Our policy
# data has plain WGS84 long/lat, so transform the polygons to match (EPSG 4326).
belgium_shape <- st_transform(belgium_shape, 4326)

# Identify the postcode attribute column (named POSTCODE in this file; the
# grep keeps this robust if the shapefile version differs) and standardize it.
pc_col <- grep("post", names(belgium_shape), ignore.case = TRUE, value = TRUE)[1]
pc_col

## [1] "POSTCODE"

belgium_shape <- belgium_shape %>%
  mutate(postcode = as.integer(as.character(.data[[pc_col]]))) %>%
  select(postcode, geometry)

nrow(belgium_shape)

## [1] 1146

Line by line:

st_transform(belgium_shape, 4326) mathematically re-projects every polygon vertex from Lambert 72 meters into WGS84 degrees. The 4326 is the EPSG code — a numeric catalogue ID for CRSs; memorize this one.
grep("post", names(...), ignore.case = TRUE, value = TRUE) searches the column names for any containing “post” (case-insensitive) and returns the matching names (value = TRUE; without it grep returns positions). Why bother? Different shapefile versions name this column POSTCODE, Postcode, nouveau_PO… Finding it by pattern makes the code robust. We print pc_col so the knitted document shows which column was found — transparency again.
.data[[pc_col]] is how you refer to a column whose name is stored in a variable inside dplyr verbs. Plain pc_col inside mutate() would be taken literally as a column named “pc_col”; the .data[[ ]] pronoun says “look up this string as a column name.”
The as.integer(as.character(...)) factor-safety idiom appears again — same reason as in cleaning: both sides of the upcoming join must be plain integers.
select(postcode, geometry) drops every other attribute column. The polygons and the key are all we need; carrying extra columns risks name collisions in later joins.

ggplot(belgium_shape) +
  geom_sf(fill = "grey95", color = "grey40", linewidth = 0.1) +
  labs(title = "Belgium: one polygon per postal code",
       subtitle = paste(nrow(belgium_shape), "building blocks")) +
  theme_void()

Your first map — and your first look at ggplot2’s grammar. A ggplot is built from layers glued with +: ggplot(data) starts the plot, geom_*() layers draw things, labs() adds titles, theme_*() controls appearance. geom_sf() is the geometry layer that understands sf objects — it reads the geometry column and draws each polygon, handling the projection automatically. Here every polygon gets the same fill (light grey) and thin border, because the point of this plot is purely to verify the shapefile: does it look like Belgium? Is it complete? Roughly the right number of polygons? theme_void() strips axes and gridlines — for maps, coordinate axes are usually clutter. Always draw the bare geometry before putting data on it — if the foundation is wrong, every later map is wrong.

0.3 Aggregate the portfolio to postcode level and join to the map

We now cross from policy-level (one row = one policy, ~160k rows) to postcode-level (one row = one postcode, ~1,100 rows). This group_by() %>% summarise() maneuver is the fundamental move of data analysis — learn it cold.

pc_stats <- mtpl %>%
  group_by(postcode) %>%
  summarise(
    long       = first(long),
    lat        = first(lat),
    exposure   = sum(expo),
    n_policies = n(),
    n_claims   = sum(nclaims),
    amount     = sum(amount),
    frequency  = n_claims / exposure,
    severity   = amount / pmax(n_claims, 1),
    pure_prem  = amount / exposure,
    .groups = "drop"
  )

summary(pc_stats$exposure)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##    1.099   73.942  141.123  249.085  277.929 4505.726

How to read it: group_by(postcode) invisibly partitions the 160k rows into ~1,100 groups. summarise() then collapses each group to a single row, computing each named expression within the group:

first(long), first(lat) — every policy in a postcode carries the same municipality centroid, so taking the first is just a way to carry the coordinates through the aggregation.
sum(expo), n(), sum(nclaims), sum(amount) — total policy-years, row count (n() is summarise’s “how many rows in this group”), total claims, total losses.
Then note something subtle and very useful: frequency = n_claims / exposure uses n_claims and exposure — columns created two lines above in the same summarise() call. dplyr allows this sequential reference, letting you build derived measures in one readable block.
pmax(n_claims, 1) in the severity denominator guards division: a postcode with zero claims would give 0/0 = NaN (“not a number”), which poisons everything downstream. pmax (parallel maximum) makes the denominator at least 1, so a no-claim postcode gets severity 0/1 = 0 instead of NaN. Cheap insurance against a whole class of bugs.
.groups = "drop" removes the grouping from the result. If you leave a data frame grouped, later verbs silently operate per-group — a common source of “why is my mutate doing that?!” confusion. Ungroup when done.

summary() on the exposure column prints min / quartiles / max. Read this output: the smallest postcodes have single-digit policy-years, the biggest have thousands. Raw pure premium on 5 policy-years is noise — one large claim triples it, zero claims makes it 0. Keep that spread in mind; it is the entire motivation for credibility smoothing.

map_pc <- belgium_shape %>% left_join(pc_stats, by = "postcode")

# QA the join, never assume it worked:
sum(is.na(map_pc$exposure))                                  # polygons with no policies (white on maps - fine)

## [1] 563

setdiff(pc_stats$postcode, belgium_shape$postcode)           # portfolio postcodes missing a polygon (should be few/none)

## integer(0)

Joins, properly. left_join(x, y, by = "postcode") keeps every row of x (the polygons) and attaches matching columns from y (the statistics) wherever postcode matches; polygons with no match get NA in the attached columns. Direction matters: we join stats onto the map (not the map onto the stats) because for drawing we want every polygon present, even ones where we sold no policies.

The two lines after the join are join quality-assurance, and they are not optional decoration — every join you ever write should be followed by a check like this, because a join that silently matches nothing produces a perfectly valid-looking data frame full of NA:

sum(is.na(map_pc$exposure)) — how many polygons got no data? (is.na() gives TRUE/FALSE per element; summing logicals counts the TRUEs — another idiom to memorize.) Some unmatched polygons are expected: postcodes where this insurer simply has no customers. They’ll render white on the maps.
setdiff(a, b) returns elements of a not in b — here, portfolio postcodes with no polygon. This direction is the dangerous one: those policies would vanish from every map. It should be empty or nearly so; if it were large, you’d investigate the key mismatch (type? leading zeros? old vs. new postcode vintages?) before proceeding.

p_expo <- ggplot(map_pc) +
  geom_sf(aes(fill = exposure), color = NA) +
  scale_fill_viridis_c(option = "mako", trans = "log10", na.value = "white",
                        name = "Exposure\n(log scale)") +
  labs(title = "Earned exposure by postcode") + theme_void(base_size = 10)

p_raw_pp <- ggplot(map_pc) +
  geom_sf(aes(fill = pmin(pure_prem, quantile(pc_stats$pure_prem, 0.95))), color = NA) +
  scale_fill_viridis_c(option = "magma", na.value = "white",
                        name = "Pure premium\n(capped at p95)") +
  labs(title = "RAW pure premium by postcode - mostly noise in thin postcodes") +
  theme_void(base_size = 10)

p_expo / p_raw_pp

Your first choropleths (a choropleth = a map whose regions are colored by a data value). New ideas in this chunk:

aes(fill = exposure) — aes() (“aesthetics”) maps data to visual properties. Fill-color-driven-by-a-column is what turns geometry into a statistical graphic. Contrast the bare map, where fill = "grey95" sat outside aes(): outside = a fixed constant, inside = data-driven. This inside/outside distinction is the single most common beginner ggplot confusion — get it straight now.
color = NA removes polygon borders — with 1,100 small polygons, borders overwhelm the fill colors.
scale_fill_viridis_c() applies a viridis palette (_c = continuous). Viridis palettes are the modern default: perceptually uniform (equal data steps look like equal color steps) and legible to colorblind readers.
trans = "log10" on the exposure map: exposure spans orders of magnitude, and on a linear scale Brussels and Antwerp would be the only visible features. A log scale spends color on ratios, revealing structure everywhere.
pmin(pure_prem, quantile(..., 0.95)) caps the fill at the 95th percentile. Without it, one or two extreme postcodes stretch the color scale so far that the rest of the country renders as one flat color. Capping the display (the data itself is untouched) is standard for heavy-tailed insurance quantities. quantile(x, 0.95) = value below which 95% of observations fall.
na.value = "white" — the polygons the join QA told us about render white (“no data”) rather than being colored as if they were zeros.
p_expo / p_raw_pp — the patchwork package overloads / to mean “stack vertically” (and | for side-by-side). We name plots (p_expo <- ggplot(...)) and compose them at the end.

Read the pure-premium map: it’s speckled — isolated extreme postcodes sitting next to ordinary ones, with no smooth pattern. That salt-and-pepper texture is the visual signature of a credibility problem: true geographic risk varies smoothly (neighboring villages share roads, weather, traffic), so the jagged variation you see is mostly sampling noise from thin data. Part 1 fixes exactly this.

0.4 Interlude - the centroid method: why we collapse polygons to points

Everything in Part 1 — the distance matrix, the credibility rings, the lat/long clustering features — quietly rests on one modeling decision: each postcode, a two-dimensional polygon, is represented by a single point. This section makes that decision explicit: what the point is, why we do it, and where it can mislead. It is worth slowing down for, because the centroid method is used constantly in spatial work and is almost never explained.

What is a centroid?

The centroid of a polygon is its geometric center of mass — the point where a cardboard cutout of the shape would balance on a pin. (It is not the average of the boundary vertices — a shape with many vertices crowded on one side would pull a vertex-average toward that side; the true area centroid integrates over the whole surface.) In sf, st_centroid() computes it for every polygon in one call.

Our dataset actually ships with a second kind of center: the long/lat columns give the municipality center — a curated, “where the town actually is” point, which can differ from the pure geometric centroid when a municipality’s populated core sits off-center in its administrative boundary (think of a town at the edge of a large forest polygon: the geometric centroid is in the trees; the municipality center is in the town). For risk work, a population-weighted or settlement point is better than the geometric centroid — the point should represent where the insured vehicles actually are. We use the dataset’s long/lat for all distance computations for exactly that reason; the shapefile’s geometric centroids below are computed only to illustrate the comparison.

Why collapse a polygon to a point at all?

Three reasons, in increasing order of importance:

“Distance between two areas” is otherwise ambiguous. Between two polygons there are infinitely many distances — nearest edges? farthest corners? average over all point pairs? Any choice is defensible and they disagree. Between two points there is exactly one distance. The ~1,100 × 1,100 distance matrix that powers the credibility rings is only well-defined because each postcode became a point.
Statistical methods eat numbers, not shapes. K-means requires each observation to be a vector of numbers. A polygon is not a number — but its centroid’s latitude and longitude are two perfectly good ones. That’s precisely how geography enters the clustering in §1.2.
Computation. Point-to-point distances for all pairs of 1,100 points is a fraction of a second. True polygon-to-polygon operations (shared borders, minimum edge distances) cost orders of magnitude more and need specialized data structures. The centroid method buys you 95% of the spatial signal for 1% of the effort.

Seeing it on the real map

# Planar geometry (centroids, buffers, areas) should be computed in a
# PROJECTED crs (meters), not in raw long/lat degrees - so: project to the
# Belgian national grid (EPSG:31370), take centroids there, transform back.
belgium_centroids <- belgium_shape %>%
  st_transform(31370) %>%
  st_centroid() %>%
  st_transform(4326)

p_country <- ggplot() +
  geom_sf(data = belgium_shape, fill = "grey95", color = "grey70", linewidth = 0.1) +
  geom_sf(data = belgium_centroids, size = 0.3, color = "#B2182B") +
  labs(title = "Every polygon collapsed to a single point",
       subtitle = "Red dots = geometric centroids of the postcode polygons") +
  theme_void(base_size = 10)

p_zoom <- ggplot() +
  geom_sf(data = belgium_shape, fill = "grey95", color = "grey60", linewidth = 0.2) +
  geom_sf(data = belgium_centroids, size = 1.8, color = "#B2182B") +
  geom_point(data = pc_stats, aes(long, lat), size = 1.8, shape = 21,
             fill = "#2166AC", color = "white", stroke = 0.3) +
  coord_sf(xlim = c(4.2, 4.6), ylim = c(50.75, 50.95)) +
  labs(title = "Zoom: the Brussels area",
       subtitle = "Red = shapefile geometric centroid; blue = the dataset's own municipality-centre long/lat") +
  theme_void(base_size = 10)

p_country / p_zoom

Notes on the code: the centroid computation is deliberately sandwiched between two st_transform() calls. Centroids, buffers, and areas are planar operations — they assume flat coordinates in consistent units. Raw long/lat degrees are neither flat nor consistent (see §1.1), so the robust habit is: project to a local metric CRS (Belgium’s is Lambert 72, EPSG:31370, in meters), do the geometry, transform back to 4326 for plotting alongside everything else. In the zoom panel, coord_sf(xlim, ylim) crops the view without filtering the data, and note the two different layer types coexisting: geom_sf() for the sf objects, plain geom_point(aes(long, lat)) for the ordinary data frame — ggplot overlays them happily as long as the coordinates are in the same CRS.

What to see in the zoom: red and blue points mostly sit close together — for small, compact urban postcodes the geometric centroid and the municipality center nearly coincide, which is why the method is safe here. Where they visibly separate, it’s the larger, odd-shaped polygons — a preview of the failure modes below.

Failure mode 1: the centroid can fall outside its own polygon

For convex shapes the centroid is always inside. For concave shapes — crescents, horseshoes, coastal municipalities wrapped around a bay — the balance point can land outside the shape entirely:

# A synthetic C-shaped "municipality"
horseshoe <- st_sfc(st_polygon(list(rbind(
  c(0, 0), c(4, 0), c(4, 1), c(1, 1), c(1, 3),
  c(4, 3), c(4, 4), c(0, 4), c(0, 0)
))))

centroid_pt <- st_centroid(horseshoe)
surface_pt  <- st_point_on_surface(horseshoe)

cc <- st_coordinates(centroid_pt)
sc <- st_coordinates(surface_pt)

ggplot() +
  geom_sf(data = horseshoe, fill = "grey90", color = "grey40") +
  geom_sf(data = centroid_pt, color = "#B2182B", size = 3) +
  geom_sf(data = surface_pt, color = "#2166AC", size = 3) +
  annotate("text", x = cc[1], y = cc[2] + 0.3, label = "st_centroid() - outside the shape!",
           color = "#B2182B", size = 3.4) +
  annotate("text", x = sc[1], y = sc[2] - 0.3, label = "st_point_on_surface()",
           color = "#2166AC", size = 3.4) +
  labs(title = "A concave region's centroid can fall outside the region") +
  theme_void(base_size = 11)

The code builds a polygon by hand — st_polygon() takes a matrix of boundary vertices (first row repeated last to close the ring), st_sfc() wraps it into a geometry column — which is a useful trick to know for constructing test shapes. st_point_on_surface() is the standard fix when you need a representative point guaranteed inside the polygon (e.g. for placing labels). For distance work the outside-the-shape centroid is usually tolerable; for point-in-polygon assignment it is a genuine bug waiting to happen.

Failure mode 2: centroid distance is not adjacency

The deeper limitation: two postcodes can share a border while their centroids are far apart, and can have nearby centroids while never touching.

rects <- data.frame(
  name = c("A", "B", "C", "D"),
  xmin = c(0, 0, 4, 6.5), xmax = c(1, 1, 6, 8.5),
  ymin = c(0, 6, 4, 4),   ymax = c(6, 12, 6, 6)
)
cent <- transform(rects, x = (xmin + xmax) / 2, y = (ymin + ymax) / 2)

ggplot() +
  geom_rect(data = rects, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
            fill = "grey92", color = "grey30") +
  annotate("segment", x = 0.5, y = 3, xend = 0.5, yend = 9,
           linetype = "dashed", color = "#B2182B") +
  annotate("segment", x = 5, y = 5, xend = 7.5, yend = 5,
           linetype = "dashed", color = "#B2182B") +
  geom_point(data = cent, aes(x, y), size = 3) +
  geom_text(data = cent, aes(x, y, label = name), vjust = -1) +
  annotate("text", x = 1.3, y = 6, hjust = 0, size = 3.3,
           label = "A and B SHARE a border,\nyet centroid distance = 6") +
  annotate("text", x = 6.25, y = 3.3, hjust = 0.5, size = 3.3,
           label = "C and D NEVER touch,\nyet centroid distance = 2.5") +
  coord_equal() +
  labs(title = "Centroid distance is not adjacency",
       subtitle = "Elongated shapes (river valleys, coastal strips) are where the two notions diverge most") +
  theme_void(base_size = 11)

Real-world shapes that behave like A and B: municipalities stretched along a river valley, a coastline, or a mountain ridge. A ring-based method sees their centroids as distant and may exclude a literal next-door neighbor from the 5 km ring, while including a compact non-neighbor across a provincial border. (This is also Jennings’ coastline caveat: near a coast or border, half of every ring is water or foreign territory, so rings capture systematically less data — one reason he suggests oval “bands” that follow a coastline as a refinement.)

Why the centroid method is nevertheless the right call here

Belgian postcode polygons are small and compact (typically a few km across) relative to our ring radii (5-50 km). When polygons are much smaller than the distances you measure between them, the centroid-vs-adjacency distortion becomes second-order: anything sharing a border is nearly always inside the 5 km ring anyway.
We use the municipality-center coordinates, which are population-representative — closer to “where the risks live” than any geometric abstraction.
The alternative — true polygon adjacency via spdep::poly2nb() and contiguity-constrained clustering (ClustGeo, skater) — is the principled upgrade, listed in “Where to take it next.” It costs significant additional machinery; the centroid method is the standard first implementation, and for a first pass through this material it is the right trade.

One further assumption worth naming out loud: by giving every policy its postcode’s single point, we assume all policies in a postcode live at the same place. Within-postcode variation (the safe cul-de-sac vs. the arterial road) is invisible to this entire analysis — territorial rating at postcode grain can never be sharper than the postcode itself.

Part 1 - Credibility smoothing, then clustering into territories

1.1 Jennings’ concentric-ring credibility smoothing

The theory, before the code. Classical credibility theory answers: “how much should I trust this group’s own experience versus a broader benchmark?” The workhorse formula is

\[\text{estimate} = z \cdot (\text{own experience}) + (1 - z) \cdot (\text{benchmark}), \qquad z \in [0, 1]\]

where the credibility weight $z$ grows with data volume. Jennings (CAS 2008) adapted this to geography with one elegant idea: the best benchmark for a postcode is not the national average — it’s the neighborhood (“principle of locality”: nearby places share risk drivers). His scheme:

The postcode’s own experience gets $z = \dfrac{\text{exposure}}{\text{exposure} + K}$ — the Bühlmann form: more exposure → more weight, with $K$ setting how much exposure earns 50% credibility.
The complement $(1-z)$ is not dumped on the national average but assigned to expanding concentric rings around the postcode (we use 5, 10, 15, 25, 50 km — Belgium is small; Jennings used miles). Each ring’s credibility is computed by the same formula on the ring’s cumulative data, minus what’s already been assigned — so credibilities sum to at most 1.
Whatever remains after the widest ring goes to the national average.
For frequency, the credibility formula is the claim-count-based $z = \min\left(1, \sqrt{n / K}\right)$ — the classical limited-fluctuation (“square root”) rule.

The result: a thick urban postcode keeps mostly its own experience; a thin rural one borrows mostly from its 10-25 km surroundings; and only truly data-poor corners fall back to the national number.

First we need to know how far apart all the postcodes are:

# ~1,100 x 1,100 great-circle distance matrix between postcode centroids - small enough to hold in memory
dist_km <- geosphere::distm(cbind(pc_stats$long, pc_stats$lat)) / 1000
dim(dist_km)

## [1] 583 583

cbind() binds vectors as columns, producing an n×2 matrix of (longitude, latitude) — the input format geosphere expects (longitude first!).
geosphere::distm() computes the great-circle distance — distance along the Earth’s curved surface — between every pair of points, returning an n×n matrix in meters; / 1000 converts to km. You cannot use plain Euclidean distance on raw long/lat degrees: one degree of longitude is ~111 km at the equator but only ~70 km at Belgian latitudes, so naive degree-distance is badly distorted.
Why precompute the full matrix? The smoothing algorithm will repeatedly ask “which postcodes are within r km of postcode i?” — with the matrix in memory, that’s one cheap row lookup (dist_km[i, ]). An n×n matrix of doubles at n≈1,100 is ~10 MB — trivially fine. (At n = 40,000 US ZIPs it would be ~13 GB — not fine; you’d switch to sf::st_is_within_distance or a spatial index. Always sanity-check memory before scaling up a matrix-based method.)
dim() prints the dimensions — confirm it’s square and matches nrow(pc_stats).

Aside: why “great-circle”? Degrees are not distances

It’s worth seeing the distortion that distm() corrects. A degree of latitude is ~111 km everywhere on Earth, but a degree of longitude shrinks as the meridians converge toward the poles — at Belgian latitude it’s only ~70 km:

# 1 degree NORTH from Brussels vs 1 degree EAST from Brussels vs 1 degree EAST at the equator
geosphere::distHaversine(c(4.35, 50.85), c(4.35, 51.85)) / 1000

## [1] 111.3195

geosphere::distHaversine(c(4.35, 50.85), c(5.35, 50.85)) / 1000

## [1] 70.28134

geosphere::distHaversine(c(30, 0),      c(31, 0))       / 1000

## [1] 111.3195

(distHaversine() is the single-pair version of distm(); points are c(longitude, latitude) — longitude first.) The same “1 degree” is 111 km going north but 70 km going east — a 60% stretch. Treating long/lat as flat x/y coordinates and using Euclidean distance would therefore make every east-west neighborhood look 60% farther away than an equally distant north-south one: credibility rings would become tall ellipses, silently borrowing too much from the north and south and too little from the east and west. Great-circle distance measures along the curved surface and is immune. Rule of habit: degrees are for locating things; meters/km from a proper distance function are for measuring things.

Seeing the rings on the map

Before running the algorithm, look at what “concentric rings around a postcode” literally means for one real postcode — this picture is the whole method:

# Pick a postcode with roughly median exposure - a "typical" thin-ish one
example_i  <- which.min(abs(pc_stats$exposure - median(pc_stats$exposure)))
example_pc <- pc_stats$postcode[example_i]

# Band every postcode by its centroid distance to the example postcode
ring_band <- cut(dist_km[example_i, ],
                 breaks = c(-Inf, 0, 5, 10, 15, 25, 50, Inf),
                 labels = c("itself", "0-5 km", "5-10 km", "10-15 km",
                            "15-25 km", "25-50 km", "beyond 50 km"))

ring_map <- belgium_shape %>%
  left_join(tibble::tibble(postcode = pc_stats$postcode, ring = ring_band),
            by = "postcode")

# Draw the literal circles: buffer the centre point in a metric CRS, transform back
centre_pt <- st_sfc(st_point(c(pc_stats$long[example_i], pc_stats$lat[example_i])),
                    crs = 4326)
circles <- st_sf(
  radius_km = c(5, 10, 15, 25, 50),
  geometry  = do.call(c, lapply(c(5, 10, 15, 25, 50) * 1000, function(r) {
    centre_pt %>% st_transform(31370) %>% st_buffer(r) %>% st_transform(4326)
  }))
)

ggplot() +
  geom_sf(data = ring_map, aes(fill = ring), color = NA) +
  geom_sf(data = circles, fill = NA, color = "black", linewidth = 0.3, linetype = "dashed") +
  geom_sf(data = centre_pt, color = "black", size = 1.5) +
  scale_fill_brewer(palette = "YlOrRd", direction = -1, na.value = "white", name = "Ring") +
  labs(title = paste0("Concentric credibility rings around postcode ", example_pc),
       subtitle = "A postcode belongs to a ring if its CENTROID falls inside - the smoothing walks outward ring by ring") +
  theme_void(base_size = 10)

How the illustration is built:

The example postcode is chosen deterministically — the one whose exposure is closest to the portfolio median (which.min(abs(x - median(x))) is a handy “find the most typical row” idiom). A median-exposure postcode is the interesting case: credible enough to keep some weight itself, thin enough to genuinely need its neighbors.
cut() bands the example’s row of the distance matrix into the exact radii the algorithm will use. Note this reuses the same dist_km[i, ] lookup the smoothing function performs — the map is a faithful picture of what the code will do, not an artist’s impression.
The dashed circles are made with st_buffer() — “grow a shape outward by a distance.” Same projection discipline as the centroid interlude: buffer by meters in EPSG:31370, then transform back to 4326 for plotting. Buffering a long/lat point “by 5” would attempt a 5-degree buffer — a circle the size of France. do.call(c, lapply(...)) builds the five buffers in a loop and concatenates them into one geometry column; attaching it to a data frame of radii via st_sf() yields a tidy sf object of five circles.
scale_fill_brewer(palette = "YlOrRd", direction = -1) — a ColorBrewer sequential palette (yellow → red), reversed so the postcode itself is darkest: visual weight matches credibility weight, fading outward exactly as the algorithm’s influence does.

Read the picture like the algorithm: the black dot’s own experience gets z = expo/(expo + 5000) of the weight; the darkest surrounding band gets the next slice; each paler band gets a smaller slice; everything in the palest zone contributes only through the national average. Also notice the imperfection the picture exposes honestly — polygons straddling a dashed circle are assigned by where their centroid falls (a polygon can be mostly outside the 10 km circle yet count as “5-10 km”) — the centroid method’s approximation, §0.4, visible in the wild. And if your example postcode sits near the French or Dutch border, part of every ring is empty foreign territory: Jennings’ border caveat, also visible.

Now the main event. This chunk is long because it is written the way production actuarial code should be: validate inputs → clean → compute benchmarks → define functions → run → validate outputs. The internal comment banners mark each section; the walkthrough after the chunk follows the same order.

radii_km <- c(0, 5, 10, 15, 25, 50)

# Credibility constants
K_expo   <- 5000
K_claims <- 1082


# ------------------------------------------------------------
# Basic validation
# ------------------------------------------------------------

required_columns <- c(
  "postcode",
  "exposure",
  "n_claims",
  "amount",
  "pure_prem",
  "frequency"
)

missing_columns <- setdiff(required_columns, names(pc_stats))

if (length(missing_columns) > 0) {
  stop(
    "Missing columns in pc_stats: ",
    paste(missing_columns, collapse = ", ")
  )
}

if (!exists("dist_km")) {
  stop("The object dist_km does not exist.")
}

if (
  nrow(dist_km) != nrow(pc_stats) ||
  ncol(dist_km) != nrow(pc_stats)
) {
  stop(
    "dist_km must be a square matrix with one row and column ",
    "for every row of pc_stats."
  )
}


# ------------------------------------------------------------
# Clean potentially invalid values
# ------------------------------------------------------------

pc_stats$exposure <- as.numeric(pc_stats$exposure)
pc_stats$n_claims <- as.numeric(pc_stats$n_claims)
pc_stats$amount   <- as.numeric(pc_stats$amount)
pc_stats$pure_prem <- as.numeric(pc_stats$pure_prem)
pc_stats$frequency <- as.numeric(pc_stats$frequency)

pc_stats$exposure[
  !is.finite(pc_stats$exposure) |
  pc_stats$exposure < 0
] <- 0

pc_stats$n_claims[
  !is.finite(pc_stats$n_claims) |
  pc_stats$n_claims < 0
] <- 0


# ------------------------------------------------------------
# National benchmarks
# ------------------------------------------------------------

total_exposure <- sum(
  pc_stats$exposure,
  na.rm = TRUE
)

total_claims <- sum(
  pc_stats$n_claims,
  na.rm = TRUE
)

total_amount <- sum(
  pc_stats$amount,
  na.rm = TRUE
)

if (total_exposure <= 0) {
  stop("Total exposure must be greater than zero.")
}

national_pp <- total_amount / total_exposure

national_freq <- total_claims / total_exposure


# ------------------------------------------------------------
# Credibility functions
# ------------------------------------------------------------

cred_expo <- function(idx) {
  exposure_local <- sum(
    pc_stats$exposure[idx],
    na.rm = TRUE
  )

  exposure_local / (exposure_local + K_expo)
}


cred_claims <- function(idx) {
  claims_local <- sum(
    pc_stats$n_claims[idx],
    na.rm = TRUE
  )

  min(
    1,
    sqrt(claims_local / K_claims)
  )
}


# ------------------------------------------------------------
# Smoothing function
# ------------------------------------------------------------

smooth_postcode <- function(
  i,
  values,
  credibility_function,
  national_value
) {

  assigned_credibility <- 0
  smoothed_estimate <- 0

  for (radius in radii_km) {

    nearby_idx <- which(
      is.finite(dist_km[i, ]) &
      dist_km[i, ] <= radius
    )

    valid_idx <- nearby_idx[
      is.finite(values[nearby_idx]) &
      is.finite(pc_stats$exposure[nearby_idx]) &
      pc_stats$exposure[nearby_idx] > 0
    ]

    if (length(valid_idx) == 0) {
      next
    }

    cumulative_credibility <-
      credibility_function(valid_idx)

    incremental_credibility <- max(
      cumulative_credibility - assigned_credibility,
      0
    )

    incremental_credibility <- min(
      incremental_credibility,
      1 - assigned_credibility
    )

    if (incremental_credibility <= 0) {
      next
    }

    local_estimate <- stats::weighted.mean(
      x = values[valid_idx],
      w = pc_stats$exposure[valid_idx],
      na.rm = TRUE
    )

    smoothed_estimate <-
      smoothed_estimate +
      incremental_credibility * local_estimate

    assigned_credibility <-
      assigned_credibility +
      incremental_credibility

    if (assigned_credibility >= 1) {
      break
    }
  }

  remaining_credibility <- max(
    1 - assigned_credibility,
    0
  )

  smoothed_estimate +
    remaining_credibility * national_value
}


# ------------------------------------------------------------
# Calculate smoothed pure premiums
# ------------------------------------------------------------

pp_smooth_vec <- numeric(nrow(pc_stats))

for (i in seq_len(nrow(pc_stats))) {

  pp_smooth_vec[i] <- smooth_postcode(
    i = i,
    values = pc_stats$pure_prem,
    credibility_function = cred_expo,
    national_value = national_pp
  )
}


# ------------------------------------------------------------
# Calculate smoothed frequencies
# ------------------------------------------------------------

freq_smooth_vec <- numeric(nrow(pc_stats))

for (i in seq_len(nrow(pc_stats))) {

  freq_smooth_vec[i] <- smooth_postcode(
    i = i,
    values = pc_stats$frequency,
    credibility_function = cred_claims,
    national_value = national_freq
  )
}


# ------------------------------------------------------------
# Attach results to pc_stats
# ------------------------------------------------------------

pc_stats$pp_smooth <- pp_smooth_vec
pc_stats$freq_smooth <- freq_smooth_vec


# ------------------------------------------------------------
# Validate results
# ------------------------------------------------------------

stopifnot(
  length(pc_stats$pp_smooth) == nrow(pc_stats),
  length(pc_stats$freq_smooth) == nrow(pc_stats)
)

cat("Number of postcodes:", nrow(pc_stats), "\n")

## Number of postcodes: 583

cat("National pure premium:", national_pp, "\n")

## National pure premium: 182.2445

cat("National frequency:", national_freq, "\n")

## National frequency: 0.1392056

cat(
  "Missing smoothed pure premiums:",
  sum(!is.finite(pc_stats$pp_smooth)),
  "\n"
)

## Missing smoothed pure premiums: 0

cat(
  "Missing smoothed frequencies:",
  sum(!is.finite(pc_stats$freq_smooth)),
  "\n"
)

## Missing smoothed frequencies: 0

summary(
  pc_stats[
    c(
      "exposure",
      "n_claims",
      "pure_prem",
      "pp_smooth",
      "frequency",
      "freq_smooth"
    )
  ]
)

##     exposure           n_claims        pure_prem         pp_smooth    
##  Min.   :   1.099   Min.   :  0.00   Min.   :   0.00   Min.   :138.2  
##  1st Qu.:  73.942   1st Qu.:  8.00   1st Qu.:  91.89   1st Qu.:168.9  
##  Median : 141.123   Median : 18.00   Median : 144.01   Median :181.0  
##  Mean   : 249.085   Mean   : 34.67   Mean   : 180.09   Mean   :185.0  
##  3rd Qu.: 277.929   3rd Qu.: 35.50   3rd Qu.: 214.04   3rd Qu.:198.1  
##  Max.   :4505.726   Max.   :703.00   Max.   :1896.68   Max.   :286.1  
##    frequency       freq_smooth    
##  Min.   :0.0000   Min.   :0.1047  
##  1st Qu.:0.1001   1st Qu.:0.1233  
##  Median :0.1259   Median :0.1330  
##  Mean   :0.1281   Mean   :0.1385  
##  3rd Qu.:0.1506   3rd Qu.:0.1458  
##  Max.   :0.3724   Max.   :0.2720

Walkthrough, section by section:

Constants. radii_km <- c(0, 5, 10, 15, 25, 50) — the ring radii. The leading 0 is a small trick: “within 0 km of postcode i” is postcode i itself, so the first loop iteration handles the postcode’s own experience using the same code as the rings — no special case needed. K_expo = 5000 means a postcode needs 5,000 policy-years for 50% credibility ($z = 5000/(5000+5000)$) — so with our largest postcodes in the low thousands, no postcode fully stands alone, which is deliberate. K_claims = 1082 is Jennings’ own full-credibility standard for frequency. These are judgment parameters: change them and re-knit to see how the smoothed map hardens (small K) or softens (large K).

Basic validation. Three fail-fast checks before any computation: the needed columns exist in pc_stats (via setdiff of required vs. actual names); dist_km exists (exists() checks by name — catches “forgot to run the previous chunk”, the classic notebook error); and the matrix dimensions match the data (|| is logical OR). If anything’s wrong we stop() with a message that says exactly what — five seconds to diagnose instead of an hour chasing a cryptic subscript error 200 lines later.

Clean invalid values. as.numeric() coercions force numeric types, then indexed assignment like pc_stats$exposure[!is.finite(pc_stats$exposure) | pc_stats$exposure < 0] <- 0 replaces bad values in place. Read it inside-out: build a logical mask of offending elements (is.finite() is FALSE for NA, NaN, Inf; | is element-wise OR), index by the mask, assign 0 to exactly those slots. Find the offenders, replace the offenders — a core R editing pattern.

National benchmarks. Totals with na.rm = TRUE (“remove NAs before summing” — without it a single NA makes the whole sum NA), a guard that total exposure is positive (division ahead), then the national pure premium and frequency that serve as the final fallback layer.

Credibility functions. cred_expo(idx) and cred_claims(idx) each take a vector of row indices, pool the exposure (or claims) of those rows, and return the credibility of the pooled data — one is the Bühlmann formula, the other the square-root rule, exactly as in the theory box. Note they read pc_stats from the enclosing environment rather than receiving it as an argument — fine at notebook scale (and keeps call sites short), but know that’s what’s happening. The design point: the smoothing algorithm below takes a credibility function as a parameter, so one algorithm serves both measures. Functions-as-arguments is a superpower of R — same idea as passing kmeans to fviz_nbclust() later.

The smoothing function — the heart of the method. For one postcode i, it maintains two accumulators: assigned_credibility (how much of the total weight is spoken for) and smoothed_estimate (the weight × value sum so far). Then, for each radius in order:

which(is.finite(dist_km[i, ]) & dist_km[i, ] <= radius) — find the cumulative disc: all postcodes within radius km of i (which() converts the TRUE/FALSE mask to numeric indices). At radius 0 this is just i itself.
Filter to valid_idx: members with finite values and positive exposure — a member with no exposure contributes nothing and would distort the weighted mean. If nobody valid is in range, next skips to the following radius.
cumulative_credibility = credibility of the whole disc’s pooled data. The incremental credibility of this ring is the cumulative figure minus what’s already assigned — clamped below by 0 (max(..., 0)) and above by the unassigned remainder (min(..., 1 - assigned)), so total assigned weight can never exceed 1. This subtract-the-prefix step is exactly how Jennings makes disc-based credibilities behave like ring-based ones.
weighted.mean(values[valid_idx], w = exposure[valid_idx]) — the disc’s experience, exposure-weighted so a big postcode influences the local average more than a village. This is the “own experience” of the disc.
Accumulate: add incremental_credibility × local_estimate to the estimate; add the increment to assigned_credibility; break out early if we’ve reached full credibility.

After the loop, whatever credibility remains goes to national_value. Trace the two extremes in your head: a huge urban postcode gets most of its weight at radius 0 and barely touches the rings; a tiny rural one gets almost nothing at 0, accumulates weight as the discs grow, and may still lean on the national average at the end. One algorithm, sensible behavior at both ends — that’s the elegance of the ring scheme.

Running it. pp_smooth_vec <- numeric(nrow(pc_stats)) pre-allocates a zero vector, then an explicit for loop fills it one postcode at a time (seq_len(n) is the safe way to write 1:n — it correctly produces an empty sequence if n were 0, where 1:0 disastrously gives c(1, 0)). Pre-allocate-then-fill matters: growing a vector inside a loop with c() re-copies it every iteration and turns a linear algorithm quadratic. The same loop runs twice — pure premium with cred_expo against national_pp, frequency with cred_claims against national_freq — and the results are attached as new pc_stats columns.

Validate results. stopifnot() asserts the outputs have the right length; the cat() lines print a mini-report (postcode count, national benchmarks, count of non-finite outputs — should be 0); and the final summary() puts raw and smoothed columns side by side. Check two things in that summary: the smoothed columns’ min-max range should be much narrower than the raw ones (extremes pulled toward local means — that’s smoothing working), and the means should be in the same ballpark (smoothing redistributes, it shouldn’t systematically shift the level).

Now look at what smoothing did:

map_pc <- belgium_shape %>% left_join(pc_stats, by = "postcode")
pp_cap <- quantile(pc_stats$pure_prem, 0.95)

p_before <- ggplot(map_pc) +
  geom_sf(aes(fill = pmin(pure_prem, pp_cap)), color = NA) +
  scale_fill_viridis_c(option = "magma", na.value = "white", limits = c(0, pp_cap),
                        name = "Pure premium") +
  labs(title = "Before: raw postcode pure premium") + theme_void(base_size = 10)

p_after <- ggplot(map_pc) +
  geom_sf(aes(fill = pmin(pp_smooth, pp_cap)), color = NA) +
  scale_fill_viridis_c(option = "magma", na.value = "white", limits = c(0, pp_cap),
                        name = "Pure premium") +
  labs(title = "After: credibility-smoothed pure premium") + theme_void(base_size = 10)

p_before / p_after

Two details make this comparison honest, and both are transferable habits: the join is re-run (map_pc is rebuilt) because pc_stats gained new columns since the last join — a stale merged copy is a classic notebook bug; and both maps share the same cap and, crucially, the same limits = c(0, pp_cap) on the color scale, so identical colors mean identical values across the two panels. Without shared limits, each map would auto-scale to its own range and the visual comparison would be a lie.

What you should see: the “before” map is static noise; the “after” map shows structure — elevated urban corridors (Brussels, Antwerp, Liège) against the quieter rural south — because averaging over rings cancels independent sampling noise but preserves shared geographic signal. This smoothed surface, not the raw one, is what we cluster.

1.2 Clustering features: Jennings’ four variables

Where we’re heading. K-means clustering treats each postcode as a point in a feature space and groups nearby points. Jennings’ feature choice — and ours — is four numbers per postcode:

smoothed pure premium, 2. smoothed frequency (the risk signal — kept separate because two postcodes can reach the same pure premium through different frequency/severity mixes, and that difference is worth preserving), 3. latitude, 4. longitude (included purely to make territories come out geographically contiguous — k-means knows nothing about maps unless you feed it coordinates).

Why standardize? K-means groups by Euclidean distance, and raw scales would rig the outcome: pure premium lives in the hundreds, latitude spans ~2 degrees — unstandardized, pure premium is the distance and geography is invisible. scale(x) converts each variable to a z-score ($z = (x - \bar{x}) / s$): mean 0, standard deviation 1, so all four variables enter with equal voice. (And note the deliberate lever: after standardizing, multiplying a column by 2 gives it twice the influence — that’s how you’d tune risk-vs-geography dominance.)

# ------------------------------------------------------------
# Build clustering features
# ------------------------------------------------------------

required_columns <- c(
  "postcode",
  "pp_smooth",
  "freq_smooth",
  "lat",
  "long"
)

missing_columns <- setdiff(
  required_columns,
  names(pc_stats)
)

if (length(missing_columns) > 0) {
  stop(
    "The following required columns are missing from pc_stats: ",
    paste(missing_columns, collapse = ", ")
  )
}


# ------------------------------------------------------------
# Ensure clustering variables are numeric
# ------------------------------------------------------------

pc_stats$postcode    <- as.integer(as.character(pc_stats$postcode))
pc_stats$pp_smooth   <- as.numeric(pc_stats$pp_smooth)
pc_stats$freq_smooth <- as.numeric(pc_stats$freq_smooth)
pc_stats$lat         <- as.numeric(pc_stats$lat)
pc_stats$long        <- as.numeric(pc_stats$long)


# ------------------------------------------------------------
# Keep only complete and finite observations
# ------------------------------------------------------------

valid_rows <- (
  is.finite(pc_stats$pp_smooth) &
  is.finite(pc_stats$freq_smooth) &
  is.finite(pc_stats$lat) &
  is.finite(pc_stats$long) &
  !is.na(pc_stats$postcode)
)

cluster_input <- data.frame(
  postcode = pc_stats$postcode[valid_rows],
  z_pp = as.numeric(scale(pc_stats$pp_smooth[valid_rows])),
  z_freq = as.numeric(scale(pc_stats$freq_smooth[valid_rows])),
  z_lat = as.numeric(scale(pc_stats$lat[valid_rows])),
  z_long = as.numeric(scale(pc_stats$long[valid_rows]))
)


# ------------------------------------------------------------
# Create numeric matrix for K-means
# ------------------------------------------------------------

cluster_matrix <- as.matrix(
  cluster_input[
    c(
      "z_pp",
      "z_freq",
      "z_lat",
      "z_long"
    )
  ]
)

storage.mode(cluster_matrix) <- "double"

rownames(cluster_matrix) <- as.character(
  cluster_input$postcode
)


# ------------------------------------------------------------
# Validate matrix
# ------------------------------------------------------------

if (any(!is.finite(cluster_matrix))) {
  stop("cluster_matrix contains NA, NaN or infinite values.")
}

if (nrow(cluster_matrix) < 2) {
  stop("Not enough valid postcodes to perform clustering.")
}

if (any(apply(cluster_matrix, 2, sd) == 0)) {
  stop(
    "At least one clustering variable has zero variation ",
    "and therefore cannot be standardized."
  )
}


# ------------------------------------------------------------
# Display checks
# ------------------------------------------------------------

cat("Postcodes available for clustering:", nrow(cluster_matrix), "\n")

## Postcodes available for clustering: 583

cat("Number of clustering variables:", ncol(cluster_matrix), "\n")

## Number of clustering variables: 4

head(cluster_input)

dim(cluster_matrix)

## [1] 583   4

summary(cluster_matrix)

##       z_pp             z_freq            z_lat             z_long        
##  Min.   :-2.0395   Min.   :-1.3820   Min.   :-3.1912   Min.   :-2.35720  
##  1st Qu.:-0.7031   1st Qu.:-0.6235   1st Qu.:-0.5691   1st Qu.:-0.70678  
##  Median :-0.1776   Median :-0.2254   Median : 0.1682   Median :-0.03037  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000  
##  3rd Qu.: 0.5694   3rd Qu.: 0.2984   3rd Qu.: 0.7314   3rd Qu.: 0.81107  
##  Max.   : 4.4035   Max.   : 5.4544   Max.   : 1.8733   Max.   : 2.11720

Walkthrough:

Build features / validation — same fail-fast column check as before (this chunk depends on pp_smooth/freq_smooth created in §1.1, so the check catches out-of-order execution).
Numeric coercion — belt-and-braces type forcing, including the factor-safe postcode idiom a third time.
valid_rows — a logical mask that is TRUE only where all four features are finite and the postcode is known (& chains element-wise ANDs). Filtering the mask once and building every column from [valid_rows] guarantees all columns stay aligned row-for-row — filtering each column separately could silently misalign them.
scale() and a subtlety: scale() returns a 1-column matrix with centering/scaling attributes attached, not a plain vector — a common surprise. Wrapping in as.numeric(...) strips that down to the plain numeric vector we want inside a data frame.
The matrix: kmeans() wants a pure numeric matrix, not a data frame, so we as.matrix() the four z-columns (postcode is deliberately excluded — it’s an identifier, not a feature; clustering on ID numbers would be meaningless). storage.mode(...) <- "double" forces floating-point storage, and rownames(...) <- postcode labels each row so cluster output remains traceable back to a postcode.
Validate matrix — three assertions: no non-finite cells (kmeans() refuses NAs with an unhelpful error; ours is helpful); at least 2 rows (can’t partition one point); and no zero-variance column (apply(m, 2, sd) applies sd down each column — dimension 2 means columns, 1 would mean rows; a constant column would have produced NaN z-scores).
Display checks — head() shows the first rows (eyeball that z-values look like z-values, roughly -3..+3); summary(cluster_matrix) should show near-0 means for every column. If a mean isn’t ~0, standardization went wrong — catch it here, not after clustering.

1.3 Choosing k

K-means requires you to choose the number of clusters k. There is no single right answer — there are diagnostics, and then there is judgment. We look at three diagnostics.

factoextra::fviz_nbclust(cluster_matrix, kmeans, method = "wss", k.max = 20, nstart = 10) +
  labs(title = "Elbow (within-cluster SS on the 4 clustering variables)")

factoextra::fviz_nbclust(cluster_matrix, kmeans, method = "silhouette", k.max = 20, nstart = 10) +
  labs(title = "Average silhouette")

fviz_nbclust() runs kmeans() at every k from 1 (or 2) to k.max and plots a quality measure per k. (Note it takes the function kmeans as an argument — functions are ordinary objects in R, passable like any value.)

Elbow (method = "wss"): within-cluster sum of squares — total squared distance of points to their cluster centers — always decreases as k grows (more clusters = tighter clusters, trivially). You look for the elbow: the k where the curve stops dropping steeply and starts flattening, i.e. where extra clusters stop buying much tightness.
Silhouette: for each point, compares its average distance to its own cluster vs. the nearest other cluster, giving a score in [-1, 1] (near 1 = clearly in the right cluster; near 0 = on a boundary). The plot shows the average per k — higher is better, and the peak suggests a k.

The elbow tends to be forgiving (a range of k look fine); the silhouette is more opinionated. Read both, expect them to disagree slightly, and treat them as bracketing a reasonable range, not delivering a verdict.

The third diagnostic is Jennings’ own, and it speaks the business’s language rather than the statistician’s: what fraction of the total (exposure-weighted) pure-premium variance remains within territories? If territories explained geography perfectly, within-territory variance would be tiny. Crucially this is computed on the raw, unstandardized, unsmoothed pure premium — the actual money quantity a filing cares about — not on the standardized clustering features.

homogeneity_pct <- function(cl) {
  df <- pc_stats %>% mutate(cl = cl)
  overall <- weighted.mean(df$pure_prem, df$exposure)
  within <- df %>%
    group_by(cl) %>%
    mutate(m = weighted.mean(pure_prem, exposure)) %>%
    ungroup() %>%
    summarise(w = sum(exposure * (pure_prem - m)^2)) %>% pull(w)
  total <- sum(df$exposure * (df$pure_prem - overall)^2)
  within / total
}

k_range <- 2:20
set.seed(20260719)
homog_curve <- sapply(k_range, function(k) {
  homogeneity_pct(kmeans(cluster_matrix, centers = k, nstart = 10)$cluster)
})

ggplot(tibble(k = k_range, homogeneity = homog_curve), aes(k, homogeneity)) +
  geom_line() + geom_point() +
  scale_y_continuous(labels = scales::percent) +
  labs(title = "Jennings' homogeneity metric vs. k",
       subtitle = "Within-territory share of (exposure-weighted) pure-premium variance - lower is better",
       y = "within / total variance")

How the function works: given a vector of cluster labels cl (one per postcode), attach it as a column; compute the national exposure-weighted mean pure premium; then — the key move — group_by(cl) %>% mutate(m = ...) computes each territory’s weighted mean and, because it’s mutate (not summarise), writes that group mean onto every row of the group, so the next line can take (pure_prem - m)^2 row by row. Sum those squared deviations, exposure-weighted, and divide by the total variance about the national mean. (pull(w) extracts the single number from the 1×1 result frame; ungroup() before summarising avoids accidental per-group summaries.) This is a decomposition-of-variance: within/total, where lower = territories capture more of the geographic cost variation.

sapply(k_range, function(k) ...) maps an anonymous function (defined inline, never named) over k = 2..20 and simplifies the results to a vector — the functional-programming sibling of a for loop, preferred when each iteration is independent. set.seed() fixes the random number generator so k-means’ random initializations — and therefore this whole document — reproduce exactly on every knit. Any analysis with randomness that might be reviewed, audited, or filed needs seeds.

Expected shape: steep drop, then flattening — with diminishing gains well before k = 20. Jennings made exactly this plot (his Exhibit 1) to justify his territory count to regulators.

k_chosen <- if (!is.na(params$k_territories)) params$k_territories else 10
k_chosen

## [1] 10

set.seed(20260719)
km <- kmeans(cluster_matrix, centers = k_chosen, nstart = 50)

pc_stats$Territory <- factor(paste0("T", sprintf("%02d", km$cluster)))

The first line reads the document parameter from the YAML header: in R, if/else is an expression that returns a value, so this is a one-line conditional assignment. Knit with rmarkdown::render("...", params = list(k_territories = 14)) to rerun the entire analysis at a different k without touching code — that’s what parameterized reports are for.
kmeans(cluster_matrix, centers = k_chosen, nstart = 50) runs the algorithm (place k centers → assign each point to nearest center → recompute centers → repeat to convergence). nstart = 50 matters: k-means only finds a local optimum that depends on the random starting centers, so it runs 50 different random starts and keeps the best (lowest within-SS). Jennings documented that different starts materially changed which postcodes landed where — many restarts is the standard fix. Cheap insurance at this data size.
The labeling line, inside-out: km$cluster is a vector of cluster numbers (1..k, one per row, in cluster_matrix row order — which is pc_stats row order by construction); sprintf("%02d", ...) formats with a leading zero (“01”, “02”… so “T02” sorts before “T10” — plain paste would sort “T10” before “T2”); paste0("T", ...) prefixes; factor() declares the result categorical, ready for the GLM.

k = 10 is a deliberate, defensible default here, not a law: the elbow and silhouette flatten well before 20, the Jennings curve shows shrinking marginal homogeneity gains, and ~10 territories is a size a rate manual can actually file. Override via params$k_territories and re-knit to see the sensitivity.

territory_summary <- pc_stats %>%
  group_by(Territory) %>%
  summarise(
    n_postcodes = n(),
    exposure    = sum(exposure),
    n_claims    = sum(n_claims),
    frequency   = sum(n_claims) / sum(exposure),
    pure_prem   = sum(amount) / sum(exposure),
    .groups = "drop"
  ) %>%
  arrange(desc(pure_prem))

territory_summary

The same group_by %>% summarise pattern as §0.3, now grouping by territory. One conceptual point hides in the frequency line: sum(n_claims) / sum(exposure) — the ratio of sums, not the average of the postcode-level frequencies (mean(frequency)). The ratio of sums correctly weights each postcode by its exposure; averaging ratios would let a 5-policy village count as much as Antwerp. This ratio-of-sums vs. mean-of-ratios distinction recurs in every aggregation an actuary does — be deliberate about it every time. arrange(desc(pure_prem)) sorts riskiest-first. Check in the output: territory pure premiums should spread apart (that’s discrimination between territories — the point of the exercise), and no territory should have trivially small exposure (a too-thin territory won’t support a stable GLM coefficient in Part 2).

1.4 The territory map

map_pc <- belgium_shape %>% left_join(pc_stats, by = "postcode")

ggplot(map_pc) +
  geom_sf(aes(fill = Territory), color = NA) +
  scale_fill_viridis_d(option = "turbo", na.value = "white") +
  labs(title = paste0("Rating territories from k-means (k = ", k_chosen, ")"),
       subtitle = "Postcode polygons colored by territory; lat/long in the feature set keeps them (mostly) contiguous") +
  theme_void(base_size = 11)

The join is refreshed once more (pc_stats gained Territory), and the fill scale switches to scale_fill_viridis_d — _d for discrete, because Territory is a factor, not a number (using a continuous scale on a factor errors out; this trips beginners constantly). paste0() in the title splices the chosen k into the text so the figure self-documents.

Inspect this map critically — this is where the manual adjustment stage of real territory work begins (GUIDE.md §3.5). K-means never saw the polygons, only centroid coordinates, so mostly-contiguous is the best it can promise: a stray postcode of one territory marooned inside another usually gets manually reassigned for contiguity; boundaries that split a city in half get reviewed with the sales/market side. The k-means output is the starting point, not the filed product.

dir.create("data", showWarnings = FALSE)
saveRDS(pc_stats %>% select(postcode, Territory), "data/06_postcode_territory_lookup.rds")

saveRDS() writes a single R object to disk exactly as-is — types, factor levels and all (readRDS() restores it perfectly; contrast CSV, which flattens everything to text). We save only the two columns that constitute the Stage-1 deliverable: the postcode → territory lookup table. showWarnings = FALSE silences the “directory already exists” warning on re-runs.

Part 2 - Full GLM: rating variables + territorial pricing

Why a second stage at all? Stage 1 looked at geography univariately — it compared places on their loss experience alone. But suppose territory T03’s high pure premium is partly because young drivers concentrate there. If we price T03 off its raw average and the rating plan also has a driver-age variable, young T03 drivers get surcharged twice — once through territory, once through age. The fix: put territory into a multivariate model alongside every other rating variable, and let the model attribute each euro of cost to its proper cause. A GLM coefficient for T03 then answers the right question: how much more expensive is T03 for the same driver with the same car and the same bonus-malus?

GLMs in three sentences. A generalized linear model extends linear regression: a linear score $\eta = \beta_0 + \beta_1 x_1 + \dots$ is connected to the mean of the response through a link function, and the response’s randomness is described by a family (Poisson for counts, Gamma for positive skewed amounts…). With a log link, $\log(\mu) = \eta$, so effects that add on the linear scale multiply on the real scale: $\mu = e^{\beta_0} e^{\beta_1 x_1} \cdots$. That multiplicative structure is exactly how insurance rating works — base rate × age factor × territory factor × … — which is why log-link GLMs are the industry-standard pricing model, and why $e^{\beta}$ (not $\beta$ itself) will be our “relativity.”

2.1 Rating variables at policy level

rating_data <- mtpl %>%
  left_join(pc_stats %>% select(postcode, Territory), by = "postcode") %>%
  mutate(
    AgeBand   = cut(ageph, breaks = c(17, 24, 29, 39, 49, 59, 69, Inf),
                     labels = c("18-24", "25-29", "30-39", "40-49", "50-59", "60-69", "70+")),
    PowerBand = cut(power, breaks = c(-Inf, 40, 55, 75, Inf),
                     labels = c("<=40kW", "41-55kW", "56-75kW", ">75kW")),
    CarAgeBand = cut(agec, breaks = c(-Inf, 2, 6, 11, Inf),
                      labels = c("0-2", "3-6", "7-11", "12+")),
    BmBand    = cut(bm, breaks = c(-Inf, 0, 3, 7, 11, Inf),
                     labels = c("0", "1-3", "4-7", "8-11", "12+"))
  )
stopifnot(!any(is.na(rating_data$Territory)))

# Base level = highest-exposure level, for every factor (relativity 1.00 base
# that's meaningful to communicate, instead of R's alphabetical default).
base_level <- function(data, var) {
  data %>% count(.data[[var]], wt = expo, sort = TRUE) %>% slice(1) %>% pull(1) %>% as.character()
}
for (v in c("Territory", "AgeBand", "PowerBand", "CarAgeBand", "BmBand",
            "fuel", "sex", "coverage", "use", "fleet")) {
  rating_data[[v]] <- relevel(factor(rating_data[[v]]), ref = base_level(rating_data, v))
}
base_territory <- levels(rating_data$Territory)[1]
base_territory

## [1] "T01"

Three things happen here:

1. The join that connects the stages. left_join attaches each policy’s territory via its postcode — Stage 1’s output becoming a Stage-2 input, at policy level. We deliberately do not aggregate to territory level: the GLM needs individual policies so it can separate the territory effect from the age/vehicle/bonus-malus effects, which is the entire point. The stopifnot(!any(is.na(...))) assertion right after verifies no policy missed a territory — assert your join worked; don’t assume.

2. Banding with cut(). cut(x, breaks, labels) converts a continuous variable into a factor of intervals: breaks = c(17, 24, 29, ...) makes bins (17,24], (24,29], … — intervals include their right endpoint by default, so age 24 lands in “18-24”, 25 in “25-29”. -Inf and Inf as end breaks catch everything below/above. Why band at all instead of fitting age as a straight line? Risk is rarely linear (the age-risk curve is U-shaped: young worst, middle-aged best, elderly rising again); bands let each range have its own level with zero functional-form assumptions, they match how rate manuals are actually published, and their coefficients are directly readable. The cost is losing within-band variation — choosing band edges is genuine actuarial judgment (these follow conventional market groupings).

3. Reference levels, chosen deliberately. For a factor, glm() estimates coefficients for every level except the first — the reference level, absorbed into the intercept, against which all other levels are measured. R’s default reference is whatever sorts first alphabetically — statistically harmless but terrible for communication (“all relativities are relative to… whichever band spells first?”). The custom base_level() function instead returns the highest-exposure level: count(.data[[var]], wt = expo, sort = TRUE) tallies exposure (not rows — wt = makes it a weighted count) by level and sorts descending; slice(1) takes the top row; pull(1) extracts the level. The for loop then applies relevel(factor(...), ref = ...) to all ten rating factors — relevel() moves the chosen level to first position. (Note rating_data[[v]] with double brackets: programmatic column access by name-in-a-variable, the base-R cousin of .data[[ ]].) Now every relativity in the output reads as “versus our most common risk,” which is exactly what a rate filing wants.

2.2 Frequency GLM (Poisson, log link, exposure offset)

Model: claim counts are non-negative integers, so the natural family is Poisson; with a log link the rating structure is multiplicative. One more ingredient is essential: policies have different exposures, and a policy insured for half a year should be expected to have half the claims. Formally we want to model the rate $\lambda_j$ per policy-year with $\mathbb{E}[N_j] = \text{expo}_j \cdot \lambda_j$; taking logs,

\[\log \mathbb{E}[N_j] = \underbrace{\log(\text{expo}_j)}_{\text{offset}} + \beta_0 + \beta_1 x_{1j} + \dots\]

An offset is a term in the linear predictor with coefficient fixed at 1 — not estimated. Passing offset = log(expo) bakes exposure into the model so the $\beta$s describe pure per-policy-year rates. (This is the standard actuarial frequency setup; forgetting the offset silently biases everything toward long-tenured policies.)

freq_glm <- glm(
  nclaims ~ Territory + AgeBand + BmBand + PowerBand + CarAgeBand +
            fuel + sex + coverage + use + fleet,
  offset = log(expo),
  family = poisson(link = "log"),
  data   = rating_data
)

summary(freq_glm)

## 
## Call:
## glm(formula = nclaims ~ Territory + AgeBand + BmBand + PowerBand + 
##     CarAgeBand + fuel + sex + coverage + use + fleet, family = poisson(link = "log"), 
##     data = rating_data, offset = log(expo))
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -2.2480545  0.0289226 -77.727  < 2e-16 ***
## TerritoryT02     -0.0864535  0.0270094  -3.201 0.001370 ** 
## TerritoryT03     -0.1288743  0.0310311  -4.153 3.28e-05 ***
## TerritoryT04      0.2498061  0.0306648   8.146 3.75e-16 ***
## TerritoryT05     -0.2027457  0.0528240  -3.838 0.000124 ***
## TerritoryT06      0.1235797  0.0252916   4.886 1.03e-06 ***
## TerritoryT07      0.1943171  0.0275712   7.048 1.82e-12 ***
## TerritoryT08      0.0485781  0.0258305   1.881 0.060020 .  
## TerritoryT09     -0.0137335  0.0333398  -0.412 0.680395    
## TerritoryT10      0.5562656  0.0296331  18.772  < 2e-16 ***
## AgeBand18-24      0.3094321  0.0348480   8.879  < 2e-16 ***
## AgeBand25-29      0.1577005  0.0256400   6.151 7.72e-10 ***
## AgeBand30-39      0.0006798  0.0206264   0.033 0.973709    
## AgeBand50-59     -0.0999728  0.0227536  -4.394 1.11e-05 ***
## AgeBand60-69     -0.2535277  0.0265953  -9.533  < 2e-16 ***
## AgeBand70+       -0.2138736  0.0316757  -6.752 1.46e-11 ***
## BmBand1-3         0.1709150  0.0201223   8.494  < 2e-16 ***
## BmBand4-7         0.3337320  0.0210870  15.826  < 2e-16 ***
## BmBand8-11        0.5751579  0.0228425  25.179  < 2e-16 ***
## BmBand12+         0.7500647  0.0355803  21.081  < 2e-16 ***
## PowerBand<=40kW  -0.1153583  0.0196129  -5.882 4.06e-09 ***
## PowerBand56-75kW  0.0183490  0.0186681   0.983 0.325652    
## PowerBand>75kW    0.1392857  0.0226342   6.154 7.57e-10 ***
## CarAgeBand0-2     0.0087604  0.0249727   0.351 0.725738    
## CarAgeBand3-6    -0.0563601  0.0177621  -3.173 0.001508 ** 
## CarAgeBand12+    -0.0243375  0.0210794  -1.155 0.248269    
## fueldiesel        0.1904527  0.0158365  12.026  < 2e-16 ***
## sexfemale         0.0224949  0.0163603   1.375 0.169141    
## coverageTPL+     -0.1156024  0.0177294  -6.520 7.01e-11 ***
## coverageTPL++    -0.1181105  0.0247963  -4.763 1.91e-06 ***
## usework          -0.0651480  0.0336170  -1.938 0.052629 .  
## fleet1           -0.1388834  0.0435748  -3.187 0.001436 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 89880  on 163211  degrees of freedom
## Residual deviance: 86712  on 163180  degrees of freedom
## AIC: 124468
## 
## Number of Fisher Scoring iterations: 6

dispersion <- sum(residuals(freq_glm, type = "pearson")^2) / df.residual(freq_glm)
dispersion   # ~1 = Poisson variance assumption OK; >>1 = refit quasipoisson for honest SEs

## [1] 1.16245

Reading the call: the formula nclaims ~ Territory + AgeBand + ... reads “model nclaims as a function of these predictors” — ~ separates response from predictors, + adds terms. Because every predictor is a factor, R automatically expands each into dummy variables against its reference level: this one formula line generates ~40 estimated coefficients.

Reading summary(freq_glm): the coefficient table has one row per non-reference level. Estimate is $\beta$ on the log scale — a Territory row with estimate 0.25 means that territory’s frequency is $e^{0.25} \approx 1.28$, i.e. 28% above the base territory, all else held equal (that italicized clause is what Part 2 buys us over Part 1). Negative = below base. Std. Error measures estimation precision; Pr(>|z|) is the p-value for “is this level really different from the reference?”. Don’t read p-values mechanically — a level can be statistically significant yet commercially trivial, or vice versa on thin data.

The dispersion check is the standard model-adequacy diagnostic for Poisson models, which make a strong hidden assumption: variance = mean. If reality is noisier than that (“overdispersion” — common in insurance), coefficients remain fine but standard errors are understated and p-values overconfident. The statistic — sum of squared Pearson residuals (residuals standardized by their expected spread) over residual degrees of freedom — should be ≈ 1 under the Poisson assumption. Materially above 1? Refit with family = quasipoisson(link = "log"): identical point estimates, honestly inflated standard errors.

2.3 Severity GLM (Gamma, log link, weighted by claim count)

Model: given a claim happened, how big? Claim sizes are positive, right-skewed (many moderate, few huge) — the classic family is Gamma, whose constant coefficient-of-variation assumption (spread proportional to mean) suits cost data far better than the normal’s constant-variance assumption. Log link again, for multiplicative relativities.

Two structural differences from the frequency model:

Data: only claim-having policies enter (severity is undefined otherwise) — that’s the filter(nclaims > 0, average > 0). Expect roughly an eighth of the portfolio.
weights = nclaims, not an offset: the response average is the policy’s mean claim amount. A mean of 3 claims is a more reliable observation than a single claim — same expected value, one-third the variance. GLM weights encode exactly that reliability (variance $\propto 1/w$). Remember the pairing as a rule: counts → offset; averages → weights. Mixing them up is one of the most common real errors in actuarial GLM work.

sev_data <- rating_data %>% filter(nclaims > 0, average > 0)
nrow(sev_data)

## [1] 18276

sev_glm <- glm(
  average ~ Territory + AgeBand + BmBand + PowerBand + CarAgeBand +
            fuel + sex + coverage + use + fleet,
  family  = Gamma(link = "log"),
  weights = nclaims,
  data    = sev_data
)

summary(sev_glm)

## 
## Call:
## glm(formula = average ~ Territory + AgeBand + BmBand + PowerBand + 
##     CarAgeBand + fuel + sex + coverage + use + fleet, family = Gamma(link = "log"), 
##     data = sev_data, weights = nclaims)
## 
## Coefficients:
##                   Estimate Std. Error t value    Pr(>|t|)    
## (Intercept)       7.051653   0.076578  92.085     < 2e-16 ***
## TerritoryT02      0.197379   0.071253   2.770     0.00561 ** 
## TerritoryT03      0.002567   0.081805   0.031     0.97497    
## TerritoryT04      0.086609   0.080830   1.071     0.28396    
## TerritoryT05      0.099850   0.139288   0.717     0.47347    
## TerritoryT06      0.099310   0.066727   1.488     0.13669    
## TerritoryT07     -0.062262   0.072696  -0.856     0.39175    
## TerritoryT08      0.152789   0.068105   2.243     0.02488 *  
## TerritoryT09      0.437799   0.087842   4.984 0.000000629 ***
## TerritoryT10     -0.011935   0.078274  -0.152     0.87881    
## AgeBand18-24      0.088180   0.092011   0.958     0.33789    
## AgeBand25-29      0.044714   0.067472   0.663     0.50753    
## AgeBand30-39     -0.093345   0.054329  -1.718     0.08579 .  
## AgeBand50-59     -0.144031   0.059991  -2.401     0.01636 *  
## AgeBand60-69     -0.080866   0.070036  -1.155     0.24825    
## AgeBand70+        0.110460   0.083554   1.322     0.18618    
## BmBand1-3         0.114524   0.052988   2.161     0.03068 *  
## BmBand4-7         0.098181   0.055448   1.771     0.07663 .  
## BmBand8-11        0.133710   0.060346   2.216     0.02672 *  
## BmBand12+         0.075682   0.093297   0.811     0.41726    
## PowerBand<=40kW   0.040784   0.051580   0.791     0.42913    
## PowerBand56-75kW  0.033597   0.049346   0.681     0.49598    
## PowerBand>75kW    0.029919   0.059798   0.500     0.61685    
## CarAgeBand0-2     0.061921   0.066537   0.931     0.35206    
## CarAgeBand3-6    -0.028350   0.046619  -0.608     0.54312    
## CarAgeBand12+    -0.006057   0.055677  -0.109     0.91338    
## fueldiesel       -0.019786   0.041877  -0.472     0.63659    
## sexfemale        -0.034470   0.042806  -0.805     0.42067    
## coverageTPL+     -0.196280   0.046879  -4.187 0.000028406 ***
## coverageTPL++     0.167038   0.065637   2.545     0.01094 *  
## usework           0.005119   0.088284   0.058     0.95376    
## fleet1           -0.115750   0.114989  -1.007     0.31413    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 6.944682)
## 
##     Null deviance: 41611  on 18275  degrees of freedom
## Residual deviance: 40725  on 18244  degrees of freedom
## AIC: 326397
## 
## Number of Fisher Scoring iterations: 8

Only ~1 policy in 8 has a claim, so the severity model runs on far less data than the frequency model — expect several rating variables that were highly significant for frequency to wash out here, and much wider intervals on the territorial severity relativities. That is normal (and actuarially meaningful: where you drive strongly influences whether you crash, but much less how expensive the crash is), and it’s why many filed plans vary frequency relativities by territory but keep severity flatter.

2.4 Territorial relativities, net of the rest of the plan

Time to convert coefficients into the deliverable: a relativity table. Because both models use a log link, a level’s multiplicative factor is $e^{\beta}$, and — since expected pure premium = frequency × severity — a territory’s pure-premium relativity is its frequency relativity times its severity relativity.

tidy_rels <- function(fit, prefix) {
  broom::tidy(fit, conf.int = TRUE) %>%
    dplyr::filter(startsWith(term, prefix)) %>%
    dplyr::transmute(
      level      = stringr::str_remove(term, stringr::fixed(prefix)),
      relativity = exp(estimate),
      ci_low     = exp(conf.low),
      ci_high    = exp(conf.high)
    )
}

territory_relativities <- tibble::tibble(
  Territory = levels(rating_data$Territory)
) %>%
  dplyr::left_join(
    tidy_rels(freq_glm, "Territory") %>%
      dplyr::rename(
        rel_freq    = relativity,
        rel_freq_lo = ci_low,
        rel_freq_hi = ci_high
      ),
    by = c("Territory" = "level")
  ) %>%
  dplyr::left_join(
    tidy_rels(sev_glm, "Territory") %>%
      dplyr::rename(
        rel_sev    = relativity,
        rel_sev_lo = ci_low,
        rel_sev_hi = ci_high
      ),
    by = c("Territory" = "level")
  ) %>%
  dplyr::mutate(
    # The reference territory has no estimated coefficient,
    # so its relativity and confidence limits are 1.
    dplyr::across(
      dplyr::starts_with("rel_"),
      ~ tidyr::replace_na(.x, 1)
    ),
    rel_pure_premium = rel_freq * rel_sev
  ) %>%
  dplyr::arrange(dplyr::desc(rel_pure_premium))

territory_relativities

The helper tidy_rels(): broom::tidy(fit, conf.int = TRUE) converts a fitted model into a plain data frame — one row per coefficient, with columns term, estimate, conf.low, conf.high — which we can then manipulate with ordinary dplyr instead of picking at summary() output. (That models-become-data-frames idea is broom’s whole reason to exist.) startsWith(term, prefix) keeps only rows whose term begins with “Territory” — R names factor coefficients as factor-name + level, e.g. TerritoryT03. transmute() (mutate-then-keep-only-these) strips the prefix from the term (fixed() says “match this as literal text, not as a regular expression”) and exponentiates estimate and both confidence limits — valid because $x \mapsto e^x$ is monotonic, so transforming the endpoints transforms the interval.

The assembly: we start from tibble(Territory = levels(...)) — the complete list of territories — and left-join the frequency and severity relativities onto it. Why? The reference territory has no coefficient row in either model (it’s the baseline!), so building from coefficients alone would silently drop it from the table. Built this way, the base territory appears with NAs — and the across(starts_with("rel_"), ~ replace_na(.x, 1)) line fills every relativity-column NA with exactly 1, which is by definition the base’s relativity. (across() applies a function to all matching columns at once; ~ replace_na(.x, 1) is a compact lambda — .x stands for “the column being processed.”) Then rel_pure_premium = rel_freq * rel_sev combines the two models, and we sort most-expensive-first. (The pervasive dplyr:: prefixes are explicitness, not necessity — useful defense against function-name collisions, e.g. stats::filter vs dplyr::filter.)

ggplot(territory_relativities,
       aes(x = reorder(Territory, rel_pure_premium), y = rel_pure_premium)) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  geom_pointrange(aes(ymin = rel_freq_lo * rel_sev, ymax = rel_freq_hi * rel_sev),
                  color = "#B2182B") +
  coord_flip() +
  labs(title = "Territorial pure-premium relativities (GLM-adjusted)",
       subtitle = paste0("Base = ", base_territory, " at 1.00; interval reflects frequency uncertainty"),
       x = NULL, y = "Relativity")

New plotting pieces: reorder(Territory, rel_pure_premium) sorts the axis by value instead of alphabetically (sorted dot plots are far easier to read); geom_hline(yintercept = 1, ...) draws the reference line — everything above it is surcharged territory, below it discounted; geom_pointrange() draws a dot (estimate) with a whisker (interval); and coord_flip() swaps axes so territory names sit readably on the left. The interval shown propagates the frequency CI only (holding severity at its point estimate — a pragmatic simplification, disclosed in the subtitle; a full interval for a product of two estimates needs the delta method or bootstrapping). Read it like an underwriter: whose interval is clear of 1.00 (confidently different from base)? Whose straddles it (weak evidence)? How wide is the overall spread (how much is geography worth in this portfolio)?

2.5 Stage-1 vs. Stage-2: what the GLM changed

The payoff comparison: territory relativities computed the naive way (territory pure premium ÷ national pure premium — what you’d file if you skipped the GLM) against the GLM-adjusted ones.

overall_pp <- with(pc_stats, sum(amount) / sum(exposure))

compare_tbl <- territory_summary %>%
  transmute(Territory, n_postcodes, exposure,
            stage1_relativity = pure_prem / overall_pp) %>%
  left_join(territory_relativities %>%
              select(Territory, stage2_glm_relativity = rel_pure_premium),
            by = "Territory") %>%
  mutate(gap = stage2_glm_relativity - stage1_relativity) %>%
  arrange(desc(stage1_relativity))

compare_tbl

(with(pc_stats, expr) evaluates the expression with the data frame’s columns in scope — a tidy way to avoid writing pc_stats$ three times. Note again ratio-of-sums for the overall pure premium. Inside select(), stage2_glm_relativity = rel_pure_premium renames while selecting.)

How to read the table: where the two columns diverge, the gap is the part of the raw territorial cost difference that the GLM re-attributed to other rating variables — driver-age mix, bonus-malus distribution, fleet share… — rather than to geography itself. That is precisely the double-counting Stage 2 exists to remove. Typically the GLM relativities are compressed toward 1.00 relative to the univariate ones: part of what looked like “expensive place” was really “expensive drivers who live there.” If some territory moves a lot, that’s not a bug — it’s the model telling you its raw experience was heavily confounded with its risk mix.

2.6 The territorial pricing map

map_pc <- map_pc %>%
  left_join(territory_relativities %>% select(Territory, rel_pure_premium), by = "Territory")

ggplot(map_pc) +
  geom_sf(aes(fill = rel_pure_premium), color = NA) +
  relativity_fill("Pure-premium\nrelativity") +
  labs(title = "Final territorial relativities on the postcode map",
       subtitle = paste0("GLM-adjusted, base territory ", base_territory, " = 1.00")) +
  theme_void(base_size = 11)

The summary map of the whole pipeline. This join is by Territory (each postcode inherits its territory’s relativity), and here our relativity_fill() helper from the packages chunk finally earns its keep: diverging blue-white-red centered at 1.00, so color encodes meaning — blue areas priced below base, red above, white ≈ base. A sequential palette (viridis) would be wrong here: it has no natural “neutral” point, and 1.00 is a neutral point. Choosing sequential vs. diverging scales by whether your quantity has a meaningful midpoint is a small decision that separates good maps from confusing ones. (\n in the legend name is a line break.)

The same map, interactive

A static map of 1,100 polygons has a fundamental limitation: you can see the pattern but you cannot interrogate it — there is no way to know which postcode that suspicious dark-red sliver is, or what its exposure was. An interactive map fixes exactly that: hover any polygon and its identity, territory, relativity, and supporting data appear; zoom into a city and inspect boundary decisions polygon by polygon. This is the right tool for the manual adjustment stage from §1.4 — finding a postcode marooned inside the wrong territory takes seconds by hovering, versus cross-referencing a static map against a lookup table.

The package is leaflet, an R wrapper around the Leaflet JavaScript mapping library (the one most websites’ maps are built on). It produces an htmlwidget — a self-contained interactive component embedded in the knitted HTML. Two consequences to know up front: it only works in HTML output (in a PDF you’d get a static snapshot at best), and htmlwidgets don’t cooperate reliably with knitr’s chunk cache, so this chunk sets cache=FALSE explicitly.

if (!requireNamespace("leaflet", quietly = TRUE)) install.packages("leaflet")
library(leaflet)

# Simplify polygon boundaries before handing them to the browser: full-detail
# polygons make a sluggish widget. Simplification is a planar operation, so -
# same discipline as ever - do it in meters (EPSG:31370), then return to
# 4326, which leaflet REQUIRES (it only understands WGS84 long/lat).
map_interactive <- map_pc %>%
  st_transform(31370) %>%
  st_simplify(dTolerance = 100) %>%
  st_transform(4326)

# Rate-manual-style relativity bands, centred on 1.00
rel_bins <- c(0, 0.80, 0.90, 0.97, 1.03, 1.10, 1.25, Inf)
rel_pal  <- colorBin(
  palette  = "RdBu",
  domain   = map_interactive$rel_pure_premium,
  bins     = rel_bins,
  reverse  = TRUE,      # RdBu's default puts red at the LOW end; we want red = expensive
  na.color = "#f0f0f0"
)

# One little HTML info card per polygon, shown on hover
labels <- sprintf(
  "<b>Postcode %s</b><br/>Territory: %s<br/>Relativity: %.3f<br/>Exposure: %s policy-years<br/>Raw pure premium: %s",
  map_interactive$postcode,
  as.character(map_interactive$Territory),
  map_interactive$rel_pure_premium,
  formatC(map_interactive$exposure, format = "d", big.mark = ","),
  formatC(map_interactive$pure_prem, format = "d", big.mark = ",")
) %>% lapply(htmltools::HTML)

leaflet(map_interactive) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    fillColor        = ~rel_pal(rel_pure_premium),
    fillOpacity      = 0.75,
    color            = "white",
    weight           = 0.4,
    highlightOptions = highlightOptions(weight = 2, color = "#333333",
                                        bringToFront = TRUE),
    label            = labels,
    labelOptions     = labelOptions(textsize = "12px")
  ) %>%
  addLegend(
    pal      = rel_pal,
    values   = ~rel_pure_premium,
    title    = "Pure-premium<br/>relativity",
    position = "bottomright"
  )

Walkthrough — leaflet’s grammar, and how it differs from ggplot’s:

st_simplify(dTolerance = 100) removes boundary vertices that deviate less than 100 m from a straightened line. The browser has to hold and redraw every vertex on every pan and zoom, so this is the difference between a snappy widget and a frozen tab. At Belgium’s scale a 100 m tolerance is visually undetectable. Note the now-familiar project → operate → transform-back sandwich; leaflet will silently misplace (or refuse) anything not in EPSG:4326, which is why the transform back is not optional.
colorBin() is leaflet’s palette factory: give it a palette name, the data range (domain), and cut points (bins), and it returns a function that maps numbers to colors — functions returning functions again. We bin rather than use a continuous gradient (colorNumeric()) for two reasons: bins match how rate manuals actually express relativities (discount bands / surcharge bands), and a binned legend is far easier to read against a hover value. The bins are deliberately asymmetric-looking but centred on 1.00 — a “neutral” band (0.97-1.03) flanked by widening discount and surcharge bands — and reverse = TRUE flips RdBu so red = expensive, matching our static maps. na.color handles the no-portfolio polygons in light grey, consistent with the white of the static maps.
The labels are built with sprintf(), which fills %s/%.3f placeholders from vectors — one formatted string per polygon in one vectorized call. formatC(..., big.mark = ",") inserts thousands separators. The strings contain HTML tags (<b>, <br/>), and lapply(htmltools::HTML) marks each one as trusted HTML so leaflet renders the tags instead of printing them literally.
The pipeline: where ggplot composes layers with +, leaflet chains layer functions with the familiar %>%: start from leaflet(data), add a basemap (addProviderTiles — CartoDB.Positron is a deliberately pale basemap that stays out of your data’s way), add the polygons, add the legend. Inside leaflet layers, ~ means “look this up in the data” — ~rel_pal(rel_pure_premium) plays exactly the role aes(fill = ...) plays in ggplot: it distinguishes data-driven values from fixed constants like weight = 0.4.
highlightOptions darkens a polygon’s outline as the mouse touches it (bringToFront stops the highlight hiding under neighbours) — a small thing that makes 1,100 tiny polygons feel navigable.

How to actually use it: zoom to Brussels or Antwerp and hover across a territory boundary — you can watch the relativity step change between adjacent postcodes and judge whether the boundary sits somewhere defensible. Then look for single off-color polygons inside an otherwise uniform region and hover them: low exposure plus an extreme raw pure premium is the signature of a postcode whose assignment you might override in the manual-adjustment pass. That inspection loop — spot, hover, judge — is precisely what this widget exists for, and it is how a territory review meeting with underwriters actually runs.

2.7 Diagnostics: actual vs. predicted by territory

A model used for pricing must be checked where it will be used: does predicted loss match actual loss by territory?

rating_data <- rating_data %>%
  mutate(
    fitted_claims = predict(freq_glm, type = "response"),                       # includes the exposure offset
    fitted_sev    = predict(sev_glm, newdata = rating_data, type = "response"),
    fitted_loss   = fitted_claims * fitted_sev
  )

avp <- rating_data %>%
  group_by(Territory) %>%
  summarise(
    exposure     = sum(expo),
    actual_pp    = sum(amount) / sum(expo),
    predicted_pp = sum(fitted_loss) / sum(expo),
    .groups = "drop"
  ) %>%
  mutate(avp_ratio = actual_pp / predicted_pp)

avp

predict() mechanics — three details that matter:

type = "response" returns predictions on the natural scale (expected claim counts, expected euro severities). The default, type = "link", returns the linear predictor — log scale. Forgetting this argument and treating log-scale output as real quantities is a rite-of-passage bug; consider yourself warned.
The frequency prediction includes each policy’s exposure offset — so it’s the expected claim count for that policy over its actual exposure period, not a per-year rate. That’s exactly what we want for summing to territory totals.
newdata = rating_data for the severity model: sev_glm was fit on claim-having policies only, but pricing needs an expected severity for every policy — “if this policy had a claim, how big would we expect it to be?” Supplying newdata scores the model on all rows. (The frequency call omits newdata because freq_glm was fit on all of rating_data already.)

fitted_loss = fitted_claims * fitted_sev is the frequency-severity pricing identity at policy level. The aggregation then computes actual and predicted pure premium per territory (ratio-of-sums, once more) and their ratio — the A/P ratio, the standard actuarial calibration diagnostic: 1.00 = perfectly calibrated, above 1 = model underpredicts (undercharges), below 1 = overpredicts.

ggplot(avp, aes(x = reorder(Territory, -exposure), y = avp_ratio)) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  geom_col(fill = "#4477AA") +
  coord_cartesian(ylim = c(0.5, 1.5)) +
  labs(title = "Actual / predicted pure premium by territory (in-sample)",
       subtitle = "Ordered by exposure, largest first; should hug 1.00",
       x = NULL, y = "A / P")

reorder(Territory, -exposure) puts the biggest territories first (the minus reverses the sort) — deviations in high-exposure territories matter far more than in small ones. One subtle choice: coord_cartesian(ylim = c(0.5, 1.5)) zooms without dropping data — unlike scale_y_continuous(limits = ...), which would remove out-of-range bars entirely (and with geom_col, silently distort the plot). Zoom with coord_cartesian; filter with scale_* limits only when you mean to.

Reading it honestly: bars hugging 1.00 confirm territory-level calibration — but this is in-sample (evaluated on the training data), so it’s a necessary check, not proof of predictive power. A territory stuck away from 1.00 even in-sample signals internal heterogeneity or a missing interaction (e.g. territory × age). For a real filing, refit everything on a training window and check A/P out-of-time on a later period.

Wrap-up: deliverables

final_rating_table <- territory_relativities %>%
  select(Territory, rel_freq, rel_sev, rel_pure_premium) %>%
  left_join(territory_summary %>% select(Territory, n_postcodes, exposure), by = "Territory") %>%
  arrange(desc(exposure))

final_rating_table

readr::write_csv(final_rating_table, "data/06_territory_relativities_beMTPL.csv")
# plus data/06_postcode_territory_lookup.rds saved in Part 1 - together these
# two files ARE the territorial rating plan: postcode -> territory -> relativity.

The final table joins the relativities to territory size (a relativity should always be read next to the exposure supporting it) and goes out as CSV via write_csv() — CSV rather than RDS because this file’s audience is people and other systems (pricing engines, filing exhibits, spreadsheets), whereas the RDS lookup’s audience is R. Two small files are the entire operational product of this notebook: to price a policy, look up its postcode → territory, then territory → relativity, and multiply into the rest of the rating algorithm.

What you should now be able to do

Read any shapefile into R, diagnose and fix its CRS, and join tabular data to it safely (with join QA).
Aggregate policy data to a geographic building block and compute exposure, frequency, severity, and pure premium correctly (ratio-of-sums!).
Implement ring-based credibility smoothing — and explain to a reviewer why raw small-area experience must not be used directly.
Standardize features, choose k with multiple diagnostics, run k-means reproducibly (set.seed, nstart), and map the result.
Fit and read frequency (Poisson + offset) and severity (Gamma + weights) GLMs, set meaningful reference levels, extract relativities with broom, and check calibration with A/P ratios.

Where to take it next

Real Jennings variant: cluster on smoothed severity as a fifth feature; or swap the discrete rings in §1.1 for an inverse-distance / sigmoid taper (GUIDE.md §3.2 caveat) and compare territory boundaries.
Polygon adjacency instead of centroid distance: with real polygons you can use true shared-border contiguity — spdep::poly2nb(belgium_shape) — both for ring construction and for enforcing contiguity as a hard constraint (e.g. ClustGeo, or skater in spdep).
GAM smooth surface as the modern alternative: Henckaerts et al. (2018) fit mgcv::gam(nclaims ~ s(long, lat, bs = "tp") + ..., offset = log(expo)) and then bin the fitted spatial effect into territories — clustering the model output instead of the raw experience. Same data, same shapefile, directly comparable to the k-means territories built here.
Random-effect territories: lme4::glmer(nclaims ~ ... + (1 | Territory)) credibility-shrinks thin territories automatically (Yao et al. 2023).
US/Canada transfer: replace §0.2 with tigris::zctas(cb = TRUE) or the StatCan FSA boundary file, and postcode with ZCTA/FSA — every subsequent chunk runs unchanged.

References

Henckaerts, R., Antonio, K., Clijsters, M., Verbelen, R., “A data driven binning strategy for the construction of insurance tariff classes,” Scandinavian Actuarial Journal, 2018 — the paper this dataset and shapefile come from.
Antonio, K., PE-pricing-analytics workshop repo (github.com/katrienantonio/PE-pricing-analytics) — source of the Belgian postcode shapefile and the spatial-maps workflow in §0.2.
Dutang, C., Charpentier, A., CASdatasets — ?beMTPL97.
Jennings, P.J., “Using Cluster Analysis to Define Geographical Rating Territories,” CAS 2008 DPP — the ring-credibility + k-means method; full treatment in GUIDE.md.
Goldburd, M., Khare, A., Tevet, D., “GLMs for Insurance Rating,” CAS Monograph No. 5, 2016.
Yao, J., et al., “Estimating Territory Risk Relativity Using GLMMs and Fuzzy C-Means Clustering,” Risks, 2023.

Territorial Ratemaking from the Foundation: Belgian MTPL at Postal-Code Level

A teaching reference: real data (CASdatasets::beMTPL97) + a real shapefile - build the map, credibility-smooth, cluster territories, then a full GLM with territorial relativities

How to read this document

What we are building, in one paragraph

Loading packages

Part 0 - Foundation: policy data + shapefile from scratch

0.1 The portfolio: `CASdatasets::beMTPL97`

Cleaning

0.2 Download and read the real Belgian postcode shapefile

0.3 Aggregate the portfolio to postcode level and join to the map

0.4 Interlude - the centroid method: why we collapse polygons to points

What is a centroid?

Why collapse a polygon to a point at all?

Seeing it on the real map

Failure mode 1: the centroid can fall outside its own polygon

Failure mode 2: centroid distance is not adjacency

Why the centroid method is nevertheless the right call here

Part 1 - Credibility smoothing, then clustering into territories

1.1 Jennings’ concentric-ring credibility smoothing

Aside: why “great-circle”? Degrees are not distances

Seeing the rings on the map

1.2 Clustering features: Jennings’ four variables

1.3 Choosing k

1.4 The territory map

Part 2 - Full GLM: rating variables + territorial pricing

2.1 Rating variables at policy level

2.2 Frequency GLM (Poisson, log link, exposure offset)

2.3 Severity GLM (Gamma, log link, weighted by claim count)

2.4 Territorial relativities, net of the rest of the plan

2.5 Stage-1 vs. Stage-2: what the GLM changed

2.6 The territorial pricing map

The same map, interactive

2.7 Diagnostics: actual vs. predicted by territory

Wrap-up: deliverables

What you should now be able to do

Where to take it next

References

Territorial Ratemaking from the Foundation: Belgian MTPL at Postal-Code Level

A teaching reference: real data (CASdatasets::beMTPL97) + a real shapefile - build the map, credibility-smooth, cluster territories, then a full GLM with territorial relativities

How to read this document

What we are building, in one paragraph

Loading packages

Part 0 - Foundation: policy data + shapefile from scratch

0.1 The portfolio: CASdatasets::beMTPL97

Cleaning

0.2 Download and read the real Belgian postcode shapefile

0.3 Aggregate the portfolio to postcode level and join to the map

0.4 Interlude - the centroid method: why we collapse polygons to points

What is a centroid?

Why collapse a polygon to a point at all?

Seeing it on the real map

Failure mode 1: the centroid can fall outside its own polygon

Failure mode 2: centroid distance is not adjacency

Why the centroid method is nevertheless the right call here

Part 1 - Credibility smoothing, then clustering into territories

1.1 Jennings’ concentric-ring credibility smoothing

Aside: why “great-circle”? Degrees are not distances

Seeing the rings on the map

1.2 Clustering features: Jennings’ four variables

1.3 Choosing k

1.4 The territory map

Part 2 - Full GLM: rating variables + territorial pricing

2.1 Rating variables at policy level

2.2 Frequency GLM (Poisson, log link, exposure offset)

2.3 Severity GLM (Gamma, log link, weighted by claim count)

2.4 Territorial relativities, net of the rest of the plan

2.5 Stage-1 vs. Stage-2: what the GLM changed

2.6 The territorial pricing map

The same map, interactive

2.7 Diagnostics: actual vs. predicted by territory

Wrap-up: deliverables

What you should now be able to do

Where to take it next

References

0.1 The portfolio: `CASdatasets::beMTPL97`