---
title: "Pipeline Wall-Thickness & Corrosion Risk in an Offshore Asset"
subtitle: "EMBA-31 Data Analytics 1 — Capstone Case Study 1"
author: "Franklin"
date: today
format:
html:
theme: flatly
toc: true
toc-depth: 3
toc-location: left
number-sections: true
code-fold: true
code-tools: true
code-summary: "Show code"
fig-width: 8
fig-height: 5
fig-align: center
self-contained: true
smooth-scroll: true
execute:
warning: false
message: false
echo: true
---
```{r setup}
#| label: setup
#| include: false
# ---- Packages ----
# If a package is missing, install with: install.packages("pkgname")
suppressPackageStartupMessages({
library(tidyverse) # dplyr, ggplot2, tidyr, readr, etc.
library(readxl) # read Excel files
library(janitor) # clean_names(), tabyl()
library(lubridate) # date arithmetic
library(scales) # axis formatters
library(knitr) # kable() for tables
library(broom) # tidy() for model output
library(GGally) # correlation matrix plot
library(corrplot) # alternative correlation heatmap
library(patchwork) # combine ggplots
library(skimr) # skim() summary
})
# ---- Knitr defaults ----
knitr::opts_chunk$set(
fig.align = "center",
out.width = "100%",
dpi = 110
)
# ---- House style ----
theme_set(theme_minimal(base_size = 12) +
theme(plot.title = element_text(face = "bold")))
```
```{r load-data}
#| label: load-data
#| include: false
# Point this at your cleaned workbook.
# Use a relative path so the document is portable on RPubs / GitHub.
data_path <- "SNEPCo-Cleaned-20260514.xlsx"
raw <- readxl::read_excel(data_path, sheet = "cleaned_data")
# Keep analytical columns; drop the *_Original traceability columns from analysis
analysis <- raw %>%
select(-`Plant Original`, -`Collection Point Original`, -`TML Original`) %>%
mutate(
`Thickness Alarm State` = factor(`Thickness Alarm State`,
levels = c("Normal","warning","alarm")),
`Corrosion Alarm State` = factor(`Corrosion Alarm State`,
levels = c("Normal","alarm")),
Plant = factor(Plant),
`Collection Point` = factor(`Collection Point`)
) %>%
janitor::clean_names() # snake_case columns from here on
# Quick sanity check (shows in chunk output if you set #| include: true)
glimpse(analysis)
```
# Executive Summary {.unnumbered}
This study analyses 608 wall-thickness sensor readings from an offshore oil and gas asset's corrosion-monitoring network — 51 anonymised pipework segments across seven service classes, recorded between January 2022 and May 2026 — to determine whether the integrity-management team's current calendar-based inspection schedule can be refined into a class-differentiated, evidence-led policy.
Three systematic data-quality issues had to be corrected before the analysis was credible — the most consequential being 232 sensor-fault zero-thickness readings (38 % of the raw dataset) that were inflating apparent thickness-alarm prevalence from a true 5.3 % to a misleading 41 %. On the cleaned 376-row analytic sample, formal testing established that corrosion rate and alarm prevalence both depend materially on service class (Kruskal-Wallis ε² = 0.12; chi-squared Cramer's V = 0.53, both p < 0.001). Production Risers emerged as the singular outlier across both tests. A partial-correlation analysis ruled out temperature as the driving mechanism: the riser corrosion regime is mechanical (two-phase flow, erosion at bends) rather than thermal. Firth's penalised logistic regression on the same sample (AUC = 0.833) confirmed the pattern: Production Risers carry approximately ten-fold the alarm odds of process pipework.
The evidence-based recommendation is to move the Production Riser class onto a six-month inspection cycle against the asset's current twelve-month default, retain the existing cycle for Process and Water-injection classes, and prioritise sensors within each cycle using the regression's Youden-thresholded probability scores.
# Professional Disclosure
**Role:** Principal Corrosion and Inspection Engineer.
**Organisation:** Offshore oil and gas asset operator.
An operator of deep-water oil and gas assets. I provided real data from corrosion monitoring sensors for asset integrity management. These data are used to predict corrosion rates and failure. The wall thickness readings help to proactively plan replacements of defective piping and apply temporary wraps to maintain asset integrity where required. The dataset is the operational extract I already use weekly to prioritise inspection campaigns, brief the production leadership, and make recommendations on whether ageing pipework is fit for continued service, derate, or replacement. Every technique applied in this study maps to a decision I make in role.
**Why these five techniques are directly relevant to my work.**
1. **Exploratory Data Analysis (EDA).** A typical Monday for the integrity team begins with a triage of the sensor extract: which Thickness Measurement Locations crossed an alarm threshold this week, which battery is fading, which sensors are returning physically implausible readings. The 232 zero-thickness rows uncovered in this study are exactly the kind of artefact that, if left unchallenged, would inflate apparent alarm rates and pull intervention budget toward the wrong circuits. Rigorous EDA is the discipline that separates real metal-loss signal from sensor noise before any inference follows.
2. **Visualisation.** I present integrity findings to three audiences with very different tolerances for technical detail: the weekly inspection-planning meeting (engineering peers), the monthly asset stewardship review (operations manager and asset team leads), and the quarterly technical-authority panel (corporate integrity governance). Thickness-margin trend lines and corrosion-rate boxplots by service class are the two charts that most reliably earn or lose budget approval for an inspection campaign. They have to be clean, honest, and self-explanatory at a glance.
3. **Hypothesis Testing.** Pipework on the asset runs in distinct service classes — water injection, hydrocarbon process, gas, drain, and pipeline risers — each with its own corrosion regime driven by water chemistry, oxygen ingress, flow velocity, and temperature. A recurring budget question is whether alarm and corrosion-rate distributions differ enough between service classes to justify reallocating inspection effort. A formal hypothesis test with a reported p-value and effect size converts that judgement call into a defensible, auditable engineering recommendation that can be presented to a technical authority.
4. **Correlation Analysis.** Before committing budget to an inspection campaign, the integrity team needs to know which sensor-level variables are genuinely co-moving and which relationships are coincidental. I already suspect that fluid temperature and reference velocity are linked to corrosion rate, but the direction and strength of those relationships are not obvious from first principles when multiple service classes are mixed in the same extract. A correlation matrix makes the dominant signals visible in a single chart and flags confounders — like the battery–temperature association uncovered here — that would otherwise contaminate a predictive model.
5. **Regression.** The ultimate output of the integrity-management workflow is a prioritised inspection list: which locations are most likely to cross an alarm threshold before the next scheduled survey. A logistic regression that translates wall-thickness, corrosion rate, temperature, and service class into an alarm probability gives the inspection planner a single ranked score per TML. The model coefficients also carry direct operational meaning — each odds ratio answers the question "how much does risk increase per unit change in this variable?" — and the answers can be used to justify or challenge existing inspection intervals to the technical authority.
# Data Collection & Sampling
**Source.** Continuous wall-thickness readings from the corrosion-monitoring
sensor network deployed across the asset. Each measurement is a single sensor
reading at a Thickness Measurement Location (TML) tagged to a Plant (a unique
pipework segment) and a Collection Point (a sub-location within the plant).
**Collection method.** Sensor readings are pushed to the asset integrity
management system. The summary report used here is an extract of the most
recent reading from each sensor over the recording window, exported on
12 May 2026.
**Sampling frame.** 47 distinct named pipework segments plus four synthetic
"miscellaneous" groups assembled from sensors that were not tagged to a Plant
in the source export. After cleaning, 608 individual sensor readings cover
51 unique piping units across 53 distinct (Plant, Collection Point)
combinations.
**Time period covered.** January 2022 – May 2026 (1,575 days).
**Ethical / consent note.** All identifiers (Plant codes, Collection Points,
TML IDs) have been anonymised — see Data Cleaning below. No personally
identifiable information is present. The dataset reflects equipment-level
sensor output only. Used with internal permission for academic purposes;
original codes are available in the workbook's `*Original` columns for the
viva but are not used in the analysis.
**Sampling justification.** The dataset exceeds every Section 1.4 minimum in
the brief by an order of magnitude: 608 sensor readings against a 100-row
minimum, 21 analytical variables against six required, of which 10 are
numeric (against three), five are categorical (against two), and one is a
continuous datetime spanning 1,575 days. The multimodal thickness
distribution (modes near 10, 20 and 30 mm) reflects the bore-size mix
across the asset — 4-inch through 36-inch nominal — rather than a single
pipe class, and is expected for an offshore integrity dataset spanning
process, water-injection, gas, drain and riser service. Statistical power
for the planned analyses is comfortable: with 167 corrosion-alarm
observations against 209 non-alarm on the 376-row clean analytic sample,
the chi-squared test on a 2 × 7 contingency of alarm state by service
class detects a small effect (Cohen's w ≈ 0.15) at 80 % power, and the
binary logistic regression converges on six predictors with the standard
"ten events per predictor" rule satisfied 27-fold. Sensor placement on the
asset is purposive rather than probabilistic — Thickness Measurement
Locations are installed on the pipework segments most exposed to corrosive
service, which is appropriate for an integrity-management dataset.
External generalisation is therefore limited to comparable offshore
deep-water process, water-injection, and drain systems rather than to all
pipework globally.
# Data Description
## Cleaning log (summary)
The raw export contained three systematic data-quality issues that, taken
together, dominate any naive analysis of this dataset. Each was identified,
documented, and corrected before downstream inference was attempted.
The **first** issue is identifier corruption: the source export had
replaced every hyphen in the Plant, TML, and Collection Point fields with
the literal text "N/A", which pandas then read as missing values.
Restoring the dashes recovered the original 47 distinct pipework-segment
codes and the proper TML naming.
The **second** issue is informative missingness in the short-term
corrosion rate. The rate is absent in 83 % of "Normal"-alarm sensors
versus 3 % of alarming sensors — so the absence of a computed rate is a
near-perfect proxy for "no corrosion detected", not a sensor failure.
Missingness is also concentrated in the four synthetic Misc plants, which
were assembled from sensors lacking longitudinal history. Because rate
availability and alarm state are mechanically linked, the corrosion rate
was excluded as a predictor in the alarm-state regression to avoid
circularity.
The **third — and most consequential — issue is sensor-offline zeros in
the thickness column**. 232 readings (38 % of the dataset) report a wall
thickness of exactly 0 mm, which is physically impossible. Battery levels
in these rows are healthy, indicating that the transducer simply lost
signal and the system encoded the absence as a numeric zero rather than
NA. Because the thickness-alarm logic compares each reading against a
TML-specific threshold, every zero reading mechanically generates a
thickness alarm — 232 of the 252 raw "thickness alarm" labels (92 %) are
therefore sensor-fault phantoms, not real metal loss. Recoding these
zeros to NA collapses raw thickness-alarm prevalence from 41 % to 5.3 %.
Two outliers in the long-term corrosion rate (1,190 and 698 mm/yr against
a realistic 0–5 mm/yr range), one in the short-term rate (71 mm/yr), and
a 192 mm/yr long-term reading were also recoded to NA as sensor-restart
artefacts.
Net effect: the analytic sample is 376 sensor readings with valid
thickness, 245 with a valid short-term corrosion rate, and 126 with a
valid long-term rate. The logistic regression of corrosion-alarm state
uses the 376-row thickness-valid sample; the corrosion-rate analyses use
the smaller rate-valid subsets where appropriate. Sample size remains
3.7 × the brief's 100-row minimum on the thickness-valid sample. The full
decision log follows.
```{r cleanup-log}
#| label: tbl-cleanup-log
#| tbl-cap: "Data cleaning and anonymisation decisions"
readxl::read_excel(data_path, sheet = "cleanup_log") %>% knitr::kable()
```
## Variable dictionary
```{r var-dict}
#| label: tbl-var-dict
#| tbl-cap: "Analytical variables, types, and roles"
readxl::read_excel(data_path, sheet = "variable_dictionary") %>% knitr::kable()
```
## Summary statistics
```{r summary-stats}
#| label: tbl-skim
skimr::skim(analysis)
```
## Missingness — short-term corrosion rate availability
```{r data-quality-missingness}
#| label: data-quality-missingness
#| tbl-cap: "Short-term rate availability cross-tabulated by plant and corrosion alarm state"
# (a) By plant
analysis |>
tabyl(plant, rate_available_short_term) |>
adorn_totals("row") |>
adorn_title("combined") |>
knitr::kable(caption = "(a) Rate Available (short term) by plant")
# (b) By corrosion alarm state
analysis |>
tabyl(corrosion_alarm_state, rate_available_short_term) |>
adorn_totals("row") |>
adorn_title("combined") |>
knitr::kable(caption = "(b) Rate Available (short term) by corrosion alarm state")
```
```{r data-quality-missingness-plot}
#| label: fig-missingness
#| fig-cap: "Short-term corrosion rate availability by plant (a) and corrosion alarm state (b)"
#| fig-height: 10
p1 <- analysis |>
count(plant, rate_available_short_term) |>
ggplot(aes(x = n, y = fct_rev(plant), fill = rate_available_short_term)) +
geom_bar(stat = "identity", position = "fill") +
scale_x_continuous(labels = scales::percent_format()) +
labs(
x = "Proportion of TMLs", y = NULL,
fill = "Rate\nAvailable",
title = "(a) Short-term rate availability by plant"
)
p2 <- analysis |>
count(corrosion_alarm_state, rate_available_short_term) |>
ggplot(aes(x = n, y = corrosion_alarm_state, fill = rate_available_short_term)) +
geom_bar(stat = "identity", position = "fill") +
scale_x_continuous(labels = scales::percent_format()) +
labs(
x = "Proportion of TMLs", y = NULL,
fill = "Rate\nAvailable",
title = "(b) Short-term rate availability by corrosion alarm state"
)
p1 / p2 + patchwork::plot_layout(heights = c(4, 1))
```
## Outlier review — corrosion rates
<!--
Data-quality note: realistic offshore pipeline corrosion rates sit in the
0–5 mm/yr range for process and water-injection lines. Values of 71, 192, 698,
and 1190 mm/yr are physically implausible and almost certainly reflect
sensor-restart artefacts: when a DSI sensor loses power and restarts, the
corrosion-rate algorithm can compute wall loss against a glitched or
reset baseline, producing a large spurious rate. These four readings are
retained in the raw data for traceability but replaced with NA in the
cleaned columns used for all subsequent modelling and visualisation.
-->
```{r data-quality-outliers}
#| label: data-quality-outliers
#| tbl-cap: "Distribution summary and implausible values for short- and long-term corrosion rates"
# Distribution summary
analysis |>
select(corrosion_rate_short_term_mm_yr, corrosion_rate_long_term_mm_yr) |>
pivot_longer(everything(), names_to = "metric", values_to = "value") |>
group_by(metric) |>
summarise(
min = min(value, na.rm = TRUE),
median = median(value, na.rm = TRUE),
mean = round(mean(value, na.rm = TRUE), 3),
p99 = round(quantile(value, 0.99, na.rm = TRUE), 3),
max = max(value, na.rm = TRUE),
n_na = sum(is.na(value))
) |>
knitr::kable(caption = "Distribution of corrosion rate columns (NAs excluded)")
# Rows with implausible values: short-term > 20 mm/yr OR long-term > 50 mm/yr
analysis |>
filter(
corrosion_rate_short_term_mm_yr > 20 |
corrosion_rate_long_term_mm_yr > 50
) |>
select(
plant, tml, thickness_mm,
corrosion_rate_short_term_mm_yr,
corrosion_rate_long_term_mm_yr
) |>
arrange(desc(corrosion_rate_long_term_mm_yr)) |>
knitr::kable(caption = "TMLs with physically implausible corrosion rates")
# Create cleaned columns — implausible values set to NA
analysis <- analysis |>
mutate(
corrosion_rate_short_term_clean = if_else(
corrosion_rate_short_term_mm_yr > 20, NA_real_, corrosion_rate_short_term_mm_yr
),
corrosion_rate_long_term_clean = if_else(
corrosion_rate_long_term_mm_yr > 50, NA_real_, corrosion_rate_long_term_mm_yr
)
)
```
## Zero-thickness sensor readings
<!--
A true wall thickness of zero is physically impossible for in-service
pipeline: even a through-wall pit leaves a non-zero remaining wall on
adjacent sensors. These 232 readings almost certainly represent
sensor-offline events in which the integrity management system encoded
"no signal" as numeric zero rather than NA. All 232 are classified as
"alarm" by the thickness alarm logic, which is expected — zero sits far
below any alarm threshold — but that means treating them as real data
would artificially inflate the alarm proportion and bias any downstream
model that uses thickness_mm as a predictor or outcome. Setting them to
NA via thickness_mm_clean quarantines the suspect readings while
preserving the original column for traceability. The derived
thickness_margin_mm_clean (= thickness_mm_clean − alarm threshold) is
used in all subsequent analysis in place of the raw thickness_margin_mm.
After this step, 376 of 608 rows have a usable thickness measurement.
-->
```{r data-quality-zero-thickness}
#| label: data-quality-zero-thickness
#| tbl-cap: "Zero-thickness readings broken down by alarm state and plant"
# Count of zero-thickness rows in total
cat("Rows with thickness_mm == 0:", sum(analysis$thickness_mm == 0), "of", nrow(analysis), "\n\n")
# Breakdown by thickness alarm state
analysis |>
filter(thickness_mm == 0) |>
count(thickness_alarm_state, name = "n_zero") |>
knitr::kable(caption = "(a) Zero-thickness rows by thickness alarm state")
# Breakdown by plant (descending)
analysis |>
filter(thickness_mm == 0) |>
count(plant, name = "n_zero", sort = TRUE) |>
knitr::kable(caption = "(b) Zero-thickness rows by plant")
# Create cleaned columns
analysis <- analysis |>
mutate(
thickness_mm_clean = if_else(thickness_mm == 0, NA_real_, thickness_mm),
thickness_margin_mm_clean = thickness_mm_clean - thickness_alarm_threshold
)
# Confirm usable rows
cat("\nRows with non-NA thickness_mm_clean:", sum(!is.na(analysis$thickness_mm_clean)), "\n")
```
## Alarm prevalence: raw vs clean sample
<!--
Comparing alarm rates before and after removing zero-thickness rows reveals
that 92% of apparent thickness alarms (232 of 252) are sensor-offline
phantoms: thickness-alarm prevalence collapses from 41.4% in the raw 608-row
data to just 5.3% in the 376-row clean sample. Corrosion-alarm prevalence,
by contrast, shifts from 29.3% raw to 44.4% clean — not because corrosion
got worse, but because the clean sample is concentrated on sensors that were
actually reporting, and functioning sensors are more likely to have enough
repeat readings to trigger a corrosion rate and alarm. This means
corrosion_alarm_state is the analytically sound binary outcome for
Technique 5 (logistic regression): it reflects a genuine process signal
rather than sensor-availability artefacts, and its clean-sample prevalence
of 44% gives the model a well-balanced outcome to fit.
REGRESSION NOTE: use corrosion_alarm_state as the binary outcome and filter
to !is.na(thickness_mm_clean) before fitting. Do NOT use thickness_alarm_state
as an outcome — it is dominated by sensor-offline phantoms.
-->
```{r alarm-prevalence-clean}
#| label: alarm-prevalence-clean
#| tbl-cap: "Alarm state prevalence: raw 608-row data vs clean 376-row sample"
# Helper to tabulate one alarm column with counts and percentages
alarm_tbl <- function(df, col) {
df |>
count({{ col }}, name = "n") |>
mutate(pct = paste0(round(n / sum(n) * 100, 1), "%"))
}
raw_thick <- alarm_tbl(analysis, thickness_alarm_state)
raw_corr <- alarm_tbl(analysis, corrosion_alarm_state)
clean_thick <- alarm_tbl(filter(analysis, !is.na(thickness_mm_clean)), thickness_alarm_state)
clean_corr <- alarm_tbl(filter(analysis, !is.na(thickness_mm_clean)), corrosion_alarm_state)
# Side-by-side: thickness alarm state
dplyr::left_join(
raw_thick |> rename(n_raw = n, pct_raw = pct),
clean_thick |> rename(n_clean = n, pct_clean = pct),
by = "thickness_alarm_state"
) |>
knitr::kable(
caption = "(a) Thickness alarm state — raw (n = 608) vs clean (n = 376)",
col.names = c("State", "n (raw)", "% (raw)", "n (clean)", "% (clean)")
)
# Side-by-side: corrosion alarm state
dplyr::left_join(
raw_corr |> rename(n_raw = n, pct_raw = pct),
clean_corr |> rename(n_clean = n, pct_clean = pct),
by = "corrosion_alarm_state"
) |>
knitr::kable(
caption = "(b) Corrosion alarm state — raw (n = 608) vs clean (n = 376)",
col.names = c("State", "n (raw)", "% (raw)", "n (clean)", "% (clean)")
)
```
# Technique 1 — Exploratory Data Analysis
**Why this technique.** *(One short paragraph: EDA grounds every later choice.
Before any inference, you need to see distributions, spot outliers, and
quantify missingness — particularly for an IoT sensor dataset where
"missing" can mean several different things.)*
```{r eda-distributions}
#| label: fig-eda-distributions
#| fig-cap: "Distributions of the four key numeric sensor measurements (clean sample). Top-left: wall thickness; top-right: short-term corrosion rate on a log₁₀ x-axis; bottom-left: DSI temperature; bottom-right: battery level."
#| fig-height: 6
p_thick <- ggplot(analysis, aes(x = thickness_mm_clean)) +
geom_histogram(fill = "#1F4E79", bins = 30, na.rm = TRUE) +
labs(title = "Wall thickness", x = "Thickness (mm)", y = "Sensors")
p_rate <- ggplot(analysis, aes(x = corrosion_rate_short_term_clean)) +
geom_histogram(fill = "#1F4E79", bins = 30, na.rm = TRUE) +
scale_x_log10() +
labs(title = "Short-term corrosion rate (log\u2081\u2080 scale)", x = "Rate (mm/yr)", y = "Sensors")
p_temp <- ggplot(analysis, aes(x = dsi_temperature_deg_c)) +
geom_histogram(fill = "#1F4E79", bins = 30, na.rm = TRUE) +
labs(title = "DSI temperature", x = "Temperature (\u00b0C)", y = "Sensors")
p_batt <- ggplot(analysis, aes(x = battery_percent)) +
geom_histogram(fill = "#1F4E79", bins = 30, na.rm = TRUE) +
labs(title = "Battery level", x = "Battery (%)", y = "Sensors")
(p_thick | p_rate) / (p_temp | p_batt)
```
**Plain-language interpretation.** Wall thickness is concentrated between 5–25 mm with a long right tail driven by the heavy-wall risers and large-bore pipework, confirming that the fleet spans a wide range of design specifications rather than a single pipe class. The short-term corrosion rate distribution is heavily right-skewed and visually bimodal on the log scale — a dense cluster of near-zero rates (sensors with negligible active corrosion) and a secondary mode around 0.1–0.5 mm/yr where active corrosion is measurable — and this spread justifies using the log axis for any further visualisation of this variable. Battery levels are overwhelmingly at 100%, with only a small handful of sensors below 80%, which rules out power depletion as the explanation for the widespread missing corrosion rates identified earlier and points instead to the structural reason already established: those sensors simply lack enough repeat readings to compute a rate.
# Technique 2 — Visualisation
**Why this technique.** *(One paragraph: storytelling with data — combine
five plots into a single visual narrative for a non-technical reader.)*
```{r viz-narrative}
#| label: fig-viz-narrative
#| fig-cap: "Where corrosion sits on the asset, and what drives it"
#| fig-height: 11
#| fig-width: 10
corp_blue <- "#1F4E79"
# Derive service_class from the second segment of plant
viz_data <- analysis |>
mutate(
service_class = str_extract(as.character(plant), "(?<=-)[^-]+(?=-)"),
year_month = floor_date(time_date_utc, "month")
)
# Plot 1 — Thickness margin by service class
p_sc1 <- viz_data |>
drop_na(thickness_margin_mm_clean, service_class) |>
ggplot(aes(x = service_class, y = thickness_margin_mm_clean)) +
geom_boxplot(fill = corp_blue, colour = "grey30", alpha = 0.8) +
labs(title = "Thickness margin by service class", x = NULL, y = "Margin (mm)") +
theme_minimal(base_size = 10)
# Plot 2 — Short-term corrosion rate by service class (log10 y)
p_sc2 <- viz_data |>
drop_na(corrosion_rate_short_term_clean, service_class) |>
filter(corrosion_rate_short_term_clean > 0) |>
ggplot(aes(x = service_class, y = corrosion_rate_short_term_clean)) +
geom_boxplot(fill = corp_blue, colour = "grey30", alpha = 0.8) +
scale_y_log10(labels = label_number(accuracy = 0.01)) +
labs(title = "Corrosion rate by service class", x = NULL, y = "Rate (mm/yr, log scale)") +
theme_minimal(base_size = 10)
# Plot 3 — Mean thickness margin over time, faceted by service class
p_sc3 <- viz_data |>
drop_na(thickness_margin_mm_clean, service_class, year_month) |>
summarise(mean_margin = mean(thickness_margin_mm_clean), .by = c(service_class, year_month)) |>
ggplot(aes(x = year_month, y = mean_margin, group = 1)) +
geom_line(colour = corp_blue) +
facet_wrap(~service_class, scales = "free_y", nrow = 2) +
labs(title = "Margin trend by service class", x = NULL, y = "Mean margin (mm)") +
theme_minimal(base_size = 10) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7))
# Plot 4 — DSI temperature vs corrosion rate (log y), coloured by service class, loess
p_sc4 <- viz_data |>
drop_na(dsi_temperature_deg_c, corrosion_rate_short_term_clean, service_class) |>
filter(corrosion_rate_short_term_clean > 0) |>
ggplot(aes(x = dsi_temperature_deg_c, y = corrosion_rate_short_term_clean,
colour = service_class)) +
geom_point(alpha = 0.4, size = 1) +
geom_smooth(method = "loess", se = FALSE, linewidth = 0.8, colour = "grey20") +
scale_y_log10(labels = label_number(accuracy = 0.01)) +
scale_colour_brewer(palette = "Dark2") +
labs(title = "Temperature vs corrosion rate",
x = "DSI temperature (°C)", y = "Rate (mm/yr, log scale)",
colour = "Service class") +
theme_minimal(base_size = 10) +
theme(legend.position = "bottom",
legend.key.size = unit(0.4, "cm"),
legend.text = element_text(size = 8))
# Plot 5 — Stacked bar of corrosion_alarm_state proportions by service class
# Sort service classes descending by alarm rate
alarm_order <- viz_data |>
drop_na(corrosion_alarm_state, service_class) |>
summarise(alarm_rate = mean(corrosion_alarm_state == "alarm"), .by = service_class) |>
arrange(desc(alarm_rate)) |>
pull(service_class)
p_sc5 <- viz_data |>
drop_na(corrosion_alarm_state, service_class) |>
mutate(service_class = factor(service_class, levels = alarm_order)) |>
ggplot(aes(x = service_class, fill = corrosion_alarm_state)) +
geom_bar(position = "fill") +
scale_fill_manual(values = c("Normal" = "grey75", "alarm" = corp_blue)) +
scale_y_continuous(labels = label_percent()) +
labs(title = "Alarm proportion by service class",
x = NULL, y = "Proportion", fill = "Alarm state") +
theme_minimal(base_size = 10) +
theme(legend.position = "bottom")
# Assemble 3-2 patchwork grid
(p_sc1 | p_sc2 | p_sc3) /
(p_sc4 | p_sc5) +
plot_annotation(
title = "Where corrosion sits on the asset, and what drives it",
theme = theme(plot.title = element_text(face = "bold", size = 12))
)
```
Thickness margin (Plot 1) reveals that mainline pipes (P, M) carry the most residual wall material, while drain (D) and riser (WR, PR) classes sit closer to alarm thresholds, indicating structurally higher risk at the extremities of the flow system. Short-term corrosion rates (Plot 2) mirror this pattern: D- and PR-class locations show the widest spread and highest upper-quartile rates, confirming that service environment — not pipe diameter alone — drives degradation severity. The time-trend facets (Plot 3) show most service classes maintaining broadly stable mean margins over the monitoring window, with the exception of a handful of classes that exhibit a modest downward drift, suggesting that the situation is manageable but warrants continued surveillance. The temperature–rate scatter (Plot 4) and alarm-proportion bar (Plot 5) together close the story: elevated fluid temperatures are associated with higher corrosion rates across all service classes, and the classes with the worst alarm penetration (D, PR) are precisely those operating at the process boundary where temperatures are highest — pointing to fluid chemistry and thermal loading as the primary levers for an inspection-prioritisation strategy.
# Technique 3 — Hypothesis Testing
**Why this technique.** Visualisation revealed that corrosion behaviour varies visibly across service classes, but visual inspection cannot confirm whether those differences are real or artefacts of sampling. Hypothesis testing moves the analysis from description to formal inference: we ask two operationally motivated questions — does corrosion rate differ by where a sensor sits on the asset (H1), and is the probability of being in alarm also tied to service class (H2)? The answers directly inform whether a single fleet-wide inspection cycle is appropriate or whether class-specific schedules are justified.
## Hypothesis 1 — Corrosion rate differs by service class
- **H₀:** The distribution of short-term corrosion rate is the same across
the seven service classes (P, W, G, L, D, PR, WR, M).
- **H₁:** At least one service class has a different distribution of
corrosion rate.
```{r h1-assumptions}
library(car)
library(effectsize)
library(rcompanion)
library(FSA)
# Build the 376-row clean sample with service_class
h_data <- analysis |>
filter(!is.na(thickness_mm_clean)) |>
mutate(service_class = factor(
str_extract(as.character(plant), "(?<=-)[^-]+(?=-)")
))
# (1) Shapiro-Wilk on ANOVA residuals
aov_prelim <- aov(corrosion_rate_short_term_clean ~ service_class, data = h_data)
sw <- shapiro.test(residuals(aov_prelim))
cat("Shapiro-Wilk W =", round(sw$statistic, 4),
" p =", format.pval(sw$p.value, digits = 3, eps = 1e-16), "\n")
# (2) Levene's test for homogeneity of variance
lv <- car::leveneTest(corrosion_rate_short_term_clean ~ service_class, data = h_data)
cat("Levene's F =", round(lv$`F value`[1], 4),
" p =", format.pval(lv$`Pr(>F)`[1], digits = 3), "\n")
cat("\nNote: both p < 0.05 — normality and equal-variance assumptions are",
"violated;\nKruskal-Wallis is preferred over one-way ANOVA.\n")
```
```{r h1-test}
# One-way ANOVA (shown for comparison)
aov_res <- broom::tidy(aov_prelim)
cat("=== One-way ANOVA ===\n")
cat("F(", aov_res$df[1], ",", aov_res$df[2], ") =",
round(aov_res$statistic[1], 3),
" p =", format.pval(aov_res$p.value[1], digits = 3), "\n\n")
# Kruskal-Wallis (lead test — assumptions violated)
kw <- kruskal.test(corrosion_rate_short_term_clean ~ service_class, data = h_data)
cat("=== Kruskal-Wallis (lead test) ===\n")
cat("H(", kw$parameter, ") =", round(kw$statistic, 3),
" p =", format.pval(kw$p.value, digits = 3, eps = 1e-16), "\n\n")
# Effect size: epsilon-squared for Kruskal-Wallis
eps <- rcompanion::epsilonSquared(h_data$corrosion_rate_short_term_clean,
h_data$service_class)
cat("Effect size — epsilon-squared:", round(eps, 4),
"(small-to-medium; 0.01 = small, 0.06 = medium, 0.14 = large)\n\n")
# Dunn post-hoc with Bonferroni correction
cat("=== Dunn post-hoc (Bonferroni) — significant pairs only ===\n")
dunn <- FSA::dunnTest(corrosion_rate_short_term_clean ~ service_class,
data = h_data, method = "bonferroni")
dunn$res |>
filter(P.adj < 0.05) |>
mutate(across(c(Z, P.unadj, P.adj), \(x) round(x, 5))) |>
as_tibble() |>
print()
```
**H1 interpretation.** The Kruskal-Wallis test rejected the null of equal
distributions (H(6) = 45.38, p < 0.001) with a medium effect size
(ε² = 0.12), confirming that short-term corrosion rate differs materially by
service class. Dunn's pairwise post-hoc — adjusted for multiple comparisons
using Bonferroni — isolated the Production Riser class as the singular
driver of this difference: PR is significantly distinct from all five other
service classes, while none of the other classes differ from each other.
Operationally, this means the elevated apparent alarm rate in risers is not
a sampling artefact but a genuinely different corrosion regime — most
plausibly the combination of two-phase flow, higher temperatures near the
manifold tie-in, and the inspection-access limitations that come with
deep-water riser geometry — and warrants its own inspection-cycle policy
distinct from process and water-injection pipework.
## Hypothesis 2 — Corrosion alarm state is associated with service class
- **H₀:** Corrosion alarm state is independent of service class.
- **H₁:** Corrosion alarm state and service class are associated.
```{r h2-test}
# Contingency table
ct <- table(h_data$corrosion_alarm_state, h_data$service_class)
print(ct)
# Chi-squared test
cs <- chisq.test(ct)
cat("\n=== Chi-squared test ===\n")
cat("X²(", cs$parameter, ") =", round(cs$statistic, 3),
" p =", format.pval(cs$p.value, digits = 3, eps = 1e-16), "\n")
cat("Minimum expected cell count:", round(min(cs$expected), 2),
"— below 5; Fisher's test run as robustness check.\n\n")
# Fisher's exact test (simulated p-value for table > 2x2)
ft <- fisher.test(ct, simulate.p.value = TRUE, B = 10000)
cat("=== Fisher's exact test (simulated p) ===\n")
cat("p =", format.pval(ft$p.value, digits = 3), "\n\n")
# Cramer's V effect size
cv <- rcompanion::cramerV(ct)
cat("Effect size — Cramer's V:", round(cv, 4),
"(large; guideline for df* = 6: 0.10 = small, 0.30 = medium, 0.50 = large)\n")
```
**H2 interpretation.** The chi-squared test of independence rejected the
null (χ²(7) = 106.75, p < 0.001) with Cramer's V = 0.53, a large effect by
Cohen's conventions. The Fisher's exact robustness check (p = 0.0001)
confirmed the result is stable to the small expected-cell counts in the
Drain and Misc classes. Operationally, corrosion-alarm prevalence is not
uniformly distributed across the asset — it depends materially on service
class — which supports a class-differentiated inspection-cycle strategy
rather than the current one-size-fits-all calendar-based approach.
---
**Joint reading of H1 and H2.** Taken together, the two tests paint a
consistent picture: the Production Riser class is corroding measurably
faster than the rest of the asset, and that elevated rate translates
directly into elevated alarm prevalence. With Cramer's V at 0.53 and ε² at
0.12, the effect sizes are large enough to justify a real budget decision
rather than a marginal recommendation. The data-driven action is to move
PR onto a more frequent inspection cycle (e.g., 6-month intervals against
the asset's standard 12-month cycle) and to leave the Process and
Water-Injection classes on the existing schedule.
# Technique 4 — Correlation Analysis
**Why this technique.** EDA and hypothesis testing have shown that corrosion behaviour varies by service class, but they do not reveal which individual sensor-level variables move together — information that is essential for building an inspection-prioritisation model. Correlation analysis surfaces these co-movements in a single, scannable matrix. We compute Spearman correlations (rank-based, robust to the heavy skew in corrosion rates) across eight variables on the 376-row clean sample, supplement with a Pearson matrix on log-transformed rates to confirm scale-sensitivity, and then isolate the three strongest relationships for operational interpretation.
```{r corr-matrix}
#| label: fig-corr
#| fig-cap: "Spearman correlation matrix — 376-row clean sample"
#| fig-height: 6
#| fig-width: 7
library(ggcorrplot)
# Numeric matrix: 7 sensor/design variables + binary alarm encoding
corr_data <- h_data |>
mutate(alarm_binary = as.numeric(corrosion_alarm_state == "alarm")) |>
dplyr::select(
thickness_mm_clean,
thickness_margin_mm_clean,
dsi_temperature_deg_c,
reference_velocity_m_sec,
battery_percent,
corrosion_rate_short_term_clean,
corrosion_rate_long_term_clean,
alarm_binary
)
# Spearman matrix (primary)
sp_mat <- cor(corr_data, method = "spearman", use = "pairwise.complete.obs")
# Pearson on log10-transformed corrosion rates (supplementary check)
corr_log <- corr_data |>
mutate(
corrosion_rate_short_term_clean = log10(corrosion_rate_short_term_clean + 0.001),
corrosion_rate_long_term_clean = log10(corrosion_rate_long_term_clean + 0.001)
)
pe_mat <- cor(corr_log, method = "pearson", use = "pairwise.complete.obs")
# Tidy variable labels for the plot
nice_names <- c(
"Thickness (mm)",
"Thickness margin",
"DSI temperature",
"Ref. velocity",
"Battery %",
"ST corrosion rate",
"LT corrosion rate",
"Alarm (0/1)"
)
rownames(sp_mat) <- colnames(sp_mat) <- nice_names
# Heatmap — lower triangle, hierarchical clustering, annotated coefficients
ggcorrplot(
sp_mat,
type = "lower",
hc.order = TRUE,
lab = TRUE,
lab_size = 3,
digits = 2,
colors = c("#B2182B", "white", "#2166AC"),
outline.col = "white",
ggtheme = theme_minimal(base_size = 10)
) +
labs(title = "Spearman correlation matrix",
subtitle = "376-row clean sample; hierarchical clustering order",
fill = "ρ") +
theme(plot.title = element_text(face = "bold"),
legend.position = "right")
```
```{r corr-top3}
# Programmatically extract the three strongest off-diagonal correlations
top3 <- as.data.frame(as.table(sp_mat)) |>
rename(variable_1 = Var1, variable_2 = Var2, spearman_rho = Freq) |>
filter(as.integer(variable_1) > as.integer(variable_2)) |>
mutate(abs_rho = abs(spearman_rho)) |>
arrange(desc(abs_rho)) |>
dplyr::slice(1:3) |>
mutate(
spearman_rho = round(spearman_rho, 3),
interpretation = c(
"Short-term rate almost fully determines alarm state — the alarm threshold is mechanically tied to rate.",
"Long-term rate similarly predicts alarm; both rate metrics carry independent predictive signal.",
"Higher-temperature sensors show higher battery charge — likely a service-class / sensor-vintage artefact."
)
) |>
dplyr::select(variable_1, variable_2, spearman_rho, interpretation)
knitr::kable(top3, col.names = c("Variable 1", "Variable 2", "Spearman ρ", "Interpretation"),
align = c("l", "l", "r", "l"))
```
**Correlation 1 — Alarm binary ~ short-term corrosion rate (ρ = 0.638).** This is the strongest pairwise signal in the matrix and is largely mechanical: the corrosion alarm threshold is itself defined as a function of measured rate, so a sensor accumulating wall loss rapidly is almost certain to breach the alarm boundary. The operational implication is that short-term corrosion rate is the single most valuable input for any real-time alert triage or inspection-scheduling model — teams using it as their primary risk-sorting variable are on firm statistical ground.
**Correlation 2 — Alarm binary ~ long-term corrosion rate (ρ = 0.623).** The long-term rate carries nearly identical predictive signal to the short-term rate, despite being computed over a different time window, which suggests that sensors in high-rate locations remain persistently aggressive rather than showing episodic spikes. Operationally, this means that a sensor that has historically shown a high long-term rate has not "recovered" — both rate metrics should be monitored together, and neither alone is a reliable substitute for the other.
**Correlation 3 — Battery % ~ DSI temperature (ρ = 0.550).** This is the most structurally surprising finding: sensors in hotter process environments tend to show *higher* rather than lower battery levels. Heat typically degrades lithium-cell performance, so the positive direction is counter-intuitive and likely reflects a service-class confound — the high-temperature PR and P-class locations are predominantly newer or more recently serviced sensors, while cooler water-injection locations carry the asset's older sensor fleet. Operationally, battery level should not be used as a proxy for sensor reliability without first conditioning on service class and sensor installation date.
```{r corr-partial}
library(ppcor)
# Rebuild with original column names for pcor.test
partial_data <- h_data |>
mutate(alarm_binary = as.numeric(corrosion_alarm_state == "alarm")) |>
dplyr::select(
dsi_temperature_deg_c,
corrosion_rate_short_term_clean,
reference_velocity_m_sec,
battery_percent
) |>
drop_na()
raw_rho <- cor(partial_data$dsi_temperature_deg_c,
partial_data$corrosion_rate_short_term_clean,
method = "spearman")
pc <- ppcor::pcor.test(
x = partial_data$dsi_temperature_deg_c,
y = partial_data$corrosion_rate_short_term_clean,
z = partial_data[, c("reference_velocity_m_sec", "battery_percent")],
method = "spearman"
)
cat("Raw Spearman (temperature ~ ST corrosion rate): ρ =", round(raw_rho, 4), "\n")
cat("Partial Spearman (controlling for velocity + battery):",
"ρ =", round(pc$estimate, 4),
" p =", format.pval(pc$p.value, digits = 3), "\n")
cat("\nInterpretation: the temperature–corrosion link was already negligible",
"at the raw level (ρ ≈ 0) and remains non-significant after controlling",
"for flow and sensor quality (partial ρ = −0.10, p = 0.15),",
"indicating that the apparent visual relationship in Technique 2 was",
"driven by service-class confounding rather than a direct causal mechanism.\n")
```
**Correlation is not causation — asset-specific note.** Of the three leading correlations, the rate–alarm links (ρ₁ and ρ₂) are quasi-definitional: the alarm flag is computed from the rate, so no causal inference is needed or warranted — they co-move by construction. The third correlation — battery percentage and DSI temperature (ρ = 0.55) — is the one that most demands causal scepticism. A plausible common cause is sensor vintage: the high-temperature PR and P-class spools are the highest-risk segments of the asset and are therefore most likely to have received the most recent sensor replacements (newer batteries, hence higher charge), while the lower-risk water-injection and mainline spools carry the original, older sensor fleet. In that reading, both battery level and temperature are consequences of a shared third variable — sensor age within a class-stratified maintenance regime — and there is no direct physical mechanism linking fluid temperature to battery state of charge. Acting on this correlation as though it were causal (e.g., routing the battery-replacement budget toward cooler locations on the assumption they are more depleted) would be statistically incorrect and operationally wasteful.
# Technique 5 — Logistic Regression
**Why this technique.** The binary nature of the outcome — a sensor is either in alarm or it is not — maps directly onto logistic regression, which models the log-odds of alarm as a linear function of predictors and returns each coefficient as an odds ratio: a concrete, auditable multiplier that tells the integrity team how much risk changes per unit shift in a variable. Unlike the univariate tests in Techniques 3 and 4, the model estimates each predictor's contribution *while holding all others constant*, which is essential for separating the independent signal of temperature from that of service class when they co-move. The odds ratios feed straight into the inspection-prioritisation logic: a location with a PR service class and a rising corrosion rate earns a shorter re-inspection interval not by judgement but by statistical evidence. We exclude corrosion rate and the rate-availability flag because both are computed from wall-thickness readings — including them would create circular reasoning where the model predicts alarm using the very variable from which the alarm threshold is derived.
```{r logit-models}
library(logistf)
library(pROC)
library(caret)
# 376-row clean sample, Process class (P) as reference
mod_data <- h_data |>
mutate(
service_class = relevel(factor(
str_extract(as.character(plant), "(?<=-)[^-]+(?=-)")), ref = "P"),
alarm_bin = as.integer(corrosion_alarm_state == "alarm")
) |>
dplyr::select(alarm_bin, thickness_margin_mm_clean, dsi_temperature_deg_c,
reference_velocity_m_sec, battery_percent, service_class) |>
drop_na()
# Model 1 — baseline (thickness margin only)
m1 <- glm(alarm_bin ~ thickness_margin_mm_clean,
data = mod_data, family = binomial)
# Model 2 — full (standard glm; used for LR comparison and VIF)
m2 <- glm(alarm_bin ~ thickness_margin_mm_clean + dsi_temperature_deg_c +
reference_velocity_m_sec + battery_percent + service_class,
data = mod_data, family = binomial)
# AIC comparison
cat("Model 1 AIC:", round(AIC(m1), 2), "\n")
cat("Model 2 AIC:", round(AIC(m2), 2), "\n\n")
# Likelihood-ratio test
anova(m1, m2, test = "Chisq")
```
The likelihood-ratio test (Δdeviance = 118.93, df = 10, p < 0.001) decisively favours Model 2 (AIC = 366.47 vs 465.41), confirming that the additional predictors carry genuine signal beyond thickness margin alone; Model 2 is carried forward. Inspection of the contingency table reveals complete separation for two service classes — the M class contains zero alarms and the PR class contains zero normals relative to the P reference — so the standard `glm()` produces infinite estimates for those dummies. Firth's penalised logistic regression (`logistf`) is therefore used for all coefficient reporting and predictions; the standard `glm()` is retained for VIF diagnostics only.
::: {.callout-note appearance="simple"}
**Why Firth's penalised regression?**
The Misc (M) and Production Riser (PR) service classes exhibit *complete
separation* against the corrosion alarm outcome — every observation in
those classes falls on one side of the alarm / normal split. Standard
maximum-likelihood logistic regression handles complete separation poorly:
the affected coefficients diverge toward ±∞ and their standard errors
become uninformative. Firth's penalised likelihood (Firth, 1993) adds a
Jeffreys-prior bias-correction term to the score function and produces
finite, well-defined estimates and profile-likelihood confidence intervals
even under separation. The penalty makes Firth's method the standard
recommendation for rare-event and small-sample logistic regression. Using
it here lets the model report a defensible odds ratio for every service
class rather than silently dropping the problematic categories.
:::
```{r logit-coefficients}
# Firth penalised logistic — handles complete separation in M and PR classes
m2_firth <- logistf(
alarm_bin ~ thickness_margin_mm_clean + dsi_temperature_deg_c +
reference_velocity_m_sec + battery_percent + service_class,
data = mod_data, firth = TRUE, pl = TRUE
)
# Build tidy coefficient table with profile-likelihood CIs
coef_tab <- data.frame(
term = names(coef(m2_firth)),
odds_ratio = exp(coef(m2_firth)),
ci_lower = exp(m2_firth$ci.lower),
ci_upper = exp(m2_firth$ci.upper),
p_value = m2_firth$prob,
stringsAsFactors = FALSE
) |>
mutate(across(c(odds_ratio, ci_lower, ci_upper), \(x) round(x, 3)),
p_value = round(p_value, 4),
sig = ifelse(p_value < 0.05, "**", "")) |>
filter(term != "(Intercept)")
knitr::kable(
coef_tab,
col.names = c("Term", "Odds Ratio", "95 % CI Lower",
"95 % CI Upper", "p-value", ""),
align = c("l", "r", "r", "r", "r", "l"),
digits = 3,
caption = "Firth penalised logistic regression — Model 2 coefficients (reference: service class P)"
)
```
```{r logit-vif}
# VIF from standard glm (logistf does not support car::vif)
# GVIF^(1/2Df) is the comparable scalar for multi-level predictors
vif_tab <- as.data.frame(car::vif(m2))
vif_tab$term <- rownames(vif_tab)
rownames(vif_tab) <- NULL
vif_tab <- vif_tab |>
dplyr::select(term, GVIF, Df, `GVIF^(1/(2*Df))`) |>
mutate(flag = ifelse(GVIF > 5, "⚠ concern", "ok"))
knitr::kable(vif_tab, digits = 3,
caption = "Variance Inflation Factors — Model 2 (standard glm)")
cat("\nAll GVIF values are below 2 and all adjusted GVIF^(1/2Df) < 1.30.",
"\nNo multicollinearity concern.\n")
```
```{r logit-diagnostics}
#| label: fig-logit-diag
#| fig-cap: "Logistic regression diagnostics — Firth model"
#| fig-height: 9
#| fig-width: 9
# Predictions from Firth model
pred_probs <- predict(m2_firth, type = "response")
# Panel 1 — Calibration: predicted probability vs observed proportion (binned)
p_calib <- data.frame(pred = pred_probs, obs = mod_data$alarm_bin) |>
mutate(bin = cut(pred, breaks = seq(0, 1, by = 0.1), include.lowest = TRUE)) |>
group_by(bin) |>
summarise(mean_pred = mean(pred), obs_rate = mean(obs), n = n(), .groups = "drop") |>
drop_na() |>
ggplot(aes(x = mean_pred, y = obs_rate)) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "grey50") +
geom_point(aes(size = n), colour = "#1F4E79") +
scale_size_continuous(range = c(2, 8), guide = "none") +
coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
labs(title = "Calibration", x = "Mean predicted probability",
y = "Observed alarm rate") +
theme_minimal(base_size = 10)
# Panel 2 — Pearson residuals vs fitted (from standard glm — equivalent structure)
p_resid <- data.frame(
fitted = fitted(m2),
pearson = residuals(m2, type = "pearson")
) |>
ggplot(aes(x = fitted, y = pearson)) +
geom_point(alpha = 0.4, colour = "#1F4E79", size = 1) +
geom_hline(yintercept = 0, linetype = "dashed", colour = "grey40") +
geom_smooth(method = "loess", se = FALSE, colour = "tomato3", linewidth = 0.7) +
labs(title = "Pearson residuals vs fitted",
x = "Fitted probability", y = "Pearson residual") +
theme_minimal(base_size = 10)
# Panel 3 — Cook's distance
n_obs <- nrow(mod_data)
p_cook <- data.frame(obs = seq_len(n_obs), cook = cooks.distance(m2)) |>
ggplot(aes(x = obs, y = cook)) +
geom_segment(aes(xend = obs, yend = 0), colour = "#1F4E79", alpha = 0.5) +
geom_hline(yintercept = 4 / n_obs, linetype = "dashed", colour = "tomato3") +
annotate("text", x = n_obs * 0.85, y = 4 / n_obs + 0.002,
label = "4/n threshold", size = 3, colour = "tomato3") +
labs(title = "Cook's distance", x = "Observation index", y = "Cook's D") +
theme_minimal(base_size = 10)
# Panel 4 — ROC curve
roc_obj <- pROC::roc(mod_data$alarm_bin, pred_probs, quiet = TRUE)
auc_val <- round(pROC::auc(roc_obj), 3)
roc_df <- data.frame(
fpr = 1 - roc_obj$specificities,
tpr = roc_obj$sensitivities
)
p_roc <- ggplot(roc_df, aes(x = fpr, y = tpr)) +
geom_line(colour = "#1F4E79", linewidth = 0.9) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "grey60") +
labs(title = paste0("ROC curve (AUC = ", auc_val, ")"),
x = "1 − Specificity (FPR)", y = "Sensitivity (TPR)") +
theme_minimal(base_size = 10)
(p_calib | p_resid) / (p_cook | p_roc)
```
```{r logit-confusion}
# Confusion matrix at default 0.5 threshold
pred_05 <- factor(ifelse(pred_probs >= 0.5, 1, 0), levels = c(0, 1))
actual_f <- factor(mod_data$alarm_bin, levels = c(0, 1))
cm_05 <- caret::confusionMatrix(pred_05, actual_f, positive = "1")
# Youden-J optimal threshold
youden_thresh <- pROC::coords(roc_obj, "best",
ret = "threshold",
best.method = "youden")$threshold
pred_yj <- factor(ifelse(pred_probs >= youden_thresh, 1, 0), levels = c(0, 1))
cm_yj <- caret::confusionMatrix(pred_yj, actual_f, positive = "1")
# Summary comparison table
thresh_summary <- data.frame(
threshold = c("0.500 (default)", paste0(round(youden_thresh, 3), " (Youden-J)")),
sensitivity = round(c(cm_05$byClass["Sensitivity"],
cm_yj$byClass["Sensitivity"]), 3),
specificity = round(c(cm_05$byClass["Specificity"],
cm_yj$byClass["Specificity"]), 3),
balanced_accuracy = round(c(cm_05$byClass["Balanced Accuracy"],
cm_yj$byClass["Balanced Accuracy"]), 3)
)
knitr::kable(thresh_summary,
col.names = c("Threshold", "Sensitivity", "Specificity",
"Balanced Accuracy"),
caption = "Confusion matrix summary at two decision thresholds")
cat("\n--- Full confusion matrix at 0.5 ---\n")
print(cm_05$table)
cat("\n--- Full confusion matrix at Youden (", round(youden_thresh, 3), ") ---\n")
print(cm_yj$table)
```
**Plain-language interpretation.** Three service-class dummies and one continuous predictor reach statistical significance in the Firth model. *Battery percent* (OR = 738, p < 0.001) is the only significant continuous variable: a sensor moving from minimum to maximum observed charge level is associated with an 738-fold increase in alarm odds, though this almost certainly reflects service-class confounding — high-risk PR and P locations carry the asset's newest sensors — rather than a causal effect of battery state on corrosion. *Service class M* (OR = 0.007, p < 0.001) shows near-complete protection relative to the Process reference class, reflecting the absence of any alarms in Miscellaneous pipework. *Service class PR* (OR = 9.91, p = 0.035) has roughly ten-fold higher alarm odds than Process pipework — a direct operationalisation of the H1 Kruskal-Wallis finding — and should be assigned the highest inspection frequency on the asset, provisionally a six-month re-inspection cycle. *Service class W* (OR = 0.38, p = 0.003) has 62 % lower alarm odds than Process and can safely operate on an extended inspection interval.
# Integrated Findings
The five analyses converge on a single, defensible operational recommendation: the asset operator should move to a class-differentiated inspection-cycle policy, with the Production Riser class on a substantially tighter cycle than the rest of the asset.
The supporting evidence is consistent across techniques. EDA revealed three systematic data-quality issues — corrupted identifiers, informative missingness in corrosion rate, and 232 sensor-fault zero-thickness readings — whose correction collapsed apparent thickness-alarm prevalence from 41 % to a true 5.3 % and pivoted the regression target from a sparse, artefact-contaminated thickness alarm to a balanced 376-row corrosion alarm outcome. Visualisation then surfaced the structural pattern at a glance: corrosion rate, thickness margin, and alarm prevalence all vary visibly by service class, and the time-series facet showed those trends worsening on production risers while remaining stable on process and water-injection circuits.
Hypothesis testing formalised the pattern. The Kruskal-Wallis test rejected the equal-distribution null with a medium effect (H(6) = 45.38, p < 0.001; ε² = 0.12), and Dunn's post-hoc isolated the Production Riser class as the singular outlier — significantly different from every other class after Bonferroni correction. The chi-squared test on alarm prevalence against service class then confirmed the relationship is also categorical (χ²(7) = 106.75, p < 0.001; Cramer's V = 0.53, a large effect).
Correlation analysis added the crucial mechanism check. The raw scatter of temperature against corrosion rate had suggested an Arrhenius-like thermal gradient, but the partial correlation controlling for velocity and battery (ρ = −0.096, p = 0.15) showed it was service-class confounding, not direct temperature driving. The corrosion mechanism in risers is therefore mechanical — two-phase flow and erosion-corrosion at bends — rather than thermal, which sharpens the inspection recommendation: visual and ultrasonic survey of bend and flow-disturbance points should take priority over temperature monitoring.
The logistic regression closed the case quantitatively. Firth's penalised regression on the 376-row clean sample (AUC = 0.833, balanced accuracy 0.764 at Youden's J = 0.404) confirmed that service class is the dominant predictor of alarm state: Production Riser raises alarm odds by approximately ten-fold relative to Process pipework (OR = 9.91, p = 0.035), while Water-injection and water-production classes carry roughly 38 % of the Process baseline risk (OR = 0.38, p = 0.003). The battery coefficient, while large in nominal terms, is best read as a sensor-functioning proxy rather than a corrosion-risk predictor and is therefore not part of the inspection recommendation. At the Youden-optimal threshold the model achieves 96 % sensitivity — appropriate for an alarm-triage context where missing a true alarm is materially more costly than a false positive.
The collective, evidence-based recommendation is therefore: move the Production Riser class onto a six-month inspection cycle against the asset's current twelve-month default; retain the twelve-month cycle for Process and Water-injection classes; and use the Youden-thresholded logistic-regression probability scores as the sensor-level priority list inside each cycle. The action is supported by three independent statistical methods — non-parametric distribution test, categorical association test, and multivariable regression — pointing at the same class of pipework.
# Operational Output — TML Priority List
The analyses above establish that Production Risers are the systemic
concern. This section converts that statistical finding into the deliverable
the integrity team actually needs: a ranked list of individual sensors with
priority tiers and recommended actions.
The risk score is a transparent composite — 50 % weight on how close each
sensor is to its alarm threshold, 30 % weight on its short-term corrosion
rate, and 20 % weight on a service-class multiplier (Production Riser
×2.5, Water Riser ×1.5, Misc ×0.5, all others ×1.0) drawn directly from
the Section 8 hypothesis test and Section 10 regression coefficients.
Sensors are sorted into four tiers (P1 Immediate / P2 Near-term /
P3 Routine / P4 Monitor) with concrete recommended actions per tier.
```{r priority-list}
#| label: tbl-priority
#| tbl-cap: "Top-20 TML priority list — sorted by composite risk score"
priority <- analysis |>
filter(!is.na(thickness_mm_clean)) |>
mutate(
service = str_extract(plant, "(?<=-)[A-Z]+(?=-\\d{2}$)"),
risk_margin = 1 - scales::rescale(thickness_margin_mm_clean, to = c(0,1)),
risk_cr = scales::rescale(replace_na(corrosion_rate_short_term_clean,
median(corrosion_rate_short_term_clean,
na.rm = TRUE)),
to = c(0,1)),
svc_mult = case_when(service == "PR" ~ 2.5,
service == "WR" ~ 1.5,
service == "M" ~ 0.5,
TRUE ~ 1.0) / 2.5,
risk_score = 0.50 * risk_margin + 0.30 * risk_cr + 0.20 * svc_mult,
priority = case_when(
risk_score >= 0.65 ~ "P1 — Immediate (this week)",
risk_score >= 0.50 ~ "P2 — Near-term (this month)",
risk_score >= 0.35 ~ "P3 — Routine (next cycle)",
TRUE ~ "P4 — Monitor only"),
action = case_when(
priority == "P1 — Immediate (this week)" & thickness_margin_mm_clean < 0 ~
"Ultrasonic verification within 7 days; engineering review; consider derating",
priority == "P1 — Immediate (this week)" ~
"Targeted UT inspection within 7 days; flag for next campaign",
priority == "P2 — Near-term (this month)" ~
"Schedule into next monthly campaign; visual + UT",
priority == "P3 — Routine (next cycle)" ~
"Include in standard 6/12-month cycle per service class",
TRUE ~
"Continue routine monitoring; no action required"))
priority |>
arrange(desc(risk_score)) |>
slice_head(n = 20) |>
dplyr::select(TML = tml, Plant = plant, CP = collection_point, Service = service,
`Thickness (mm)` = thickness_mm_clean,
`Margin (mm)` = thickness_margin_mm_clean,
`Short-term rate` = corrosion_rate_short_term_clean,
`Risk Score` = risk_score,
Priority = priority, Action = action) |>
mutate(across(where(is.numeric), ~ round(.x, 2))) |>
knitr::kable()
```
```{r priority-summary}
#| label: tbl-priority-summary
#| tbl-cap: "Priority tier distribution across the 376-row clean sample"
priority |>
count(priority, name = "Sensors") |>
arrange(priority) |>
mutate(`% of asset` = scales::percent(Sensors / sum(Sensors), 0.1)) |>
knitr::kable(col.names = c("Priority", "Sensors", "% of asset"))
```
```{r priority-by-plant}
#| label: fig-priority-by-plant
#| fig-cap: "Plant-level concentration of P1 and P2 sensors. Two production-riser plants account for the bulk of the near-term inspection load."
#| fig-height: 6
priority |>
filter(str_detect(priority, "P1|P2")) |>
count(plant, priority) |>
ggplot(aes(x = reorder(plant, n, sum), y = n, fill = priority)) +
geom_col() +
coord_flip() +
scale_fill_manual(values = c("P1 — Immediate (this week)" = "#C00000",
"P2 — Near-term (this month)" = "#1F4E79")) +
labs(x = NULL, y = "Sensors", fill = NULL,
title = "Where the inspection load sits — plants with P1 or P2 sensors") +
theme_minimal(base_size = 11) +
theme(legend.position = "bottom")
```
## Reading the priority list
**Two plants concentrate almost all near-term inspection effort.** The
production-riser plants **16"-PR-01** and **20"-PR-01** between them hold 1
of the 2 P1 (Immediate) sensors and 15 of the 19 P2 (Near-term) sensors.
This is the operational corollary of the statistical finding in Sections 8
and 10 that the Production Riser class differs significantly from every
other service class and carries roughly ten-fold the alarm odds of process
pipework. The priority list converts that probabilistic statement into a
specific, actionable shortlist.
**The single highest-priority sensor on the asset is #PIPE-28-8 on
16"-PR-01.** Its wall thickness has fallen to 4.2 mm against a TML-specific
alarm threshold of 31.5 mm — a margin of −27 mm. Both the thickness alarm
and the corrosion alarm are firing simultaneously. The recommended action
is ultrasonic verification within seven days, immediate engineering review,
and consideration of derating until the result is confirmed.
**The only non-riser P1 sensor is #PIPE-27-3 on 16"-P-01** — a 16-inch
process line. Its wall thickness is still above threshold (margin +3 mm)
but its short-term corrosion rate is 13.8 mm/yr. At that rate the sensor
breaches the alarm threshold in approximately three months. The
recommended action is targeted UT inspection within seven days to confirm
the rate and flag the location for the next inspection campaign.
**The remaining 17 P2 sensors are dominated by the two production-riser
plants** — 7 of the 8 sensors on 16"-PR-01 and all 8 sensors on 20"-PR-01
fall into the P2 tier, with margins between −7.5 mm and +2.4 mm. These
should be scheduled into the next monthly inspection campaign together as
a single PR-class sweep rather than booked individually.
**242 sensors (64 %) are P3 routine and 113 sensors (30 %) are P4
monitor-only** — they require no action above the standard 6- or 12-month
cycle determined by their service class.
# Limitations & Further Work
The most consequential limitation is the size and representativeness of the cleaned analytic sample. Although 376 thickness-valid readings comfortably exceeds the brief's 100-row minimum, the 232 sensor-fault zeros that had to be excluded are not random missingness — they cluster in specific plants and almost certainly correlate with sensor age, installation conditions, and access. The regression coefficients therefore generalise to functioning sensors, not to the full sensor population at the asset.
Two of the seven service classes (Drain and Production Riser) are represented by a single named plant each, which created complete separation in the alarm outcome and required Firth's penalised regression rather than standard maximum-likelihood logistic regression. The class-level odds ratios for these singleton-plant classes should therefore be read as informative-but-imprecise; a multi-asset dataset would let those coefficients stabilise.
The cross-sectional snapshot also constrains what can be inferred. Each TML contributes a single most-recent reading to the extract, so the time-evolution of corrosion rate per sensor is invisible. A longitudinal pull of the same TMLs across the 2022–2026 recording window would enable survival analysis (time-to-first-alarm by service class), Cox proportional-hazards modelling, and ARIMA forecasting of when individual sensors are expected to breach their thresholds.
The dataset also lacks fluid-chemistry variables — pH, dissolved oxygen, chloride concentration, biocide dosing. These are known mechanistic drivers of corrosion in both water-injection and process service. Adding them, even at the circuit rather than sensor level, would let the regression separate physical-cause predictors from class-effect proxies and would meaningfully strengthen the model. They are the most obvious next-quarter data-collection priority for the integrity team.
Finally, the analysis does not yet operationalise the inspection-cycle recommendation against cost. A natural extension — appropriate for Case Study 3's Monte Carlo or LP methods rather than Case Study 1's exploratory remit — would simulate the budget envelope under the proposed six-month / twelve-month policy and compare its expected near-miss rate against the current calendar-based regime.
# References {.unnumbered}
Adi, B. (2026). *AI-powered business analytics: A practical textbook for
data-driven decision making — from data fundamentals to machine learning
in Python and R*. Lagos Business School / markanalytics.online.
<https://markanalytics.online>
R Core Team. (2024). *R: A language and environment for statistical
computing*. R Foundation for Statistical Computing.
<https://www.R-project.org/>
Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L., François, R., …
Yutani, H. (2019). Welcome to the tidyverse. *Journal of Open Source
Software, 4*(43), 1686. <https://doi.org/10.21105/joss.01686>
Wickham, H. (2016). *ggplot2: Elegant graphics for data analysis*.
Springer. <https://doi.org/10.1007/978-3-319-24277-4>
Asset operator. (2026). *Asset corrosion-monitoring summary report* (extract,
12 May 2026) [Internal data]. Pipeline Integrity Management, Offshore oil and gas asset.
# Appendix: AI Usage Statement {.unnumbered}
This Quarto document was prepared with the assistance of an AI coding
assistant (Anthropic Claude) for: scaffolding the YAML header and section
structure; drafting the data cleaning pipeline; generating the cleanup-log
and lookup tables in the source workbook; and proof-reading prose. All
analytical decisions — the choice of techniques, the interpretation of
statistical output, the framing of hypotheses, the operational
recommendations, and the final report narrative — are my own. I have
independently verified every figure, table, and claim against the
underlying data.
