Part 1: Outlier Detection
Why Detect Outliers?
An outlier is an observation that differs substantially from the expected pattern. In phenological data, outliers may indicate:
- Recording errors that should be excluded from
analysis
- A typo: DOY 250 instead of DOY 150 (roughly 3 months difference)
- Wrong year entered: observation assigned to incorrect season
- Unusual weather events worth investigating
- Very warm winter causing early flowering
- Late frost delaying spring phenology
- Second flowering or other abnormal phenological
events
- Plants flowering again in autumn/winter after spring flowering
- Increasingly observed with climate change
The challenge: You can’t simply delete all outliers. Some are errors, but others are scientifically valuable observations of unusual events.
Statistical Methods for Outlier Detection
The pep_flag_outliers() function provides four
statistical methods for identifying potential outliers. Each method has
different assumptions and strengths:
Method 1: 30-Day Rule (Simple Threshold)
The simplest approach: flag any observation more than 30 days from the median for that species/phase combination.
# Flag outliers using the default 30-day rule
outliers <- pep_flag_outliers(
pep = flowering,
method = "30day",
by = c("species", "phase_id")
)
print(outliers)
#> Phenological Outlier Detection Results
#> =============================================
#> Method: 30day
#> Threshold: 30 days
#> Grouping: species, phase_id
#>
#> Total observations: 1,310,094
#> Outliers flagged: 26573 (2.03%)
#> - Early outliers: 15011
#> - Late outliers: 11562
#>
#> Deviation summary for outliers:
#> Min: -167.0 days
#> Max: 223.0 days
#> Mean |deviation|: 38.4 daysHow it works:
Distribution of flowering dates
← Early Late →
----|----[=======|=======]----|-----
^ ^ ^
-30 days median +30 days
Anything outside the brackets is flagged as an outlierUnderstanding the Output
The function adds new columns to your original data:
| Column | What It Contains | How to Use It |
|---|---|---|
is_outlier |
TRUE/FALSE flag | Filter data: data[is_outlier == FALSE] |
deviation |
Days from expected | Prioritize investigation: larger = more extreme |
expected_doy |
Reference DOY (median) | Understand what “normal” looks like |
Method 2: MAD (Median Absolute Deviation) - Recommended
The MAD method is robust to the presence of multiple outliers. This is important because if your data contains many outliers, they can distort mean-based methods.
# MAD method: flag if > 3 MAD from median (robust)
outliers_mad <- pep_flag_outliers(flowering, method = "mad", threshold = 3)Why MAD is recommended:
| Scenario | Mean-based method | MAD method |
|---|---|---|
| 1-2 outliers | Works reasonably | Works well |
| Many outliers | Outliers inflate SD, may miss more | Stays robust |
| Asymmetric data | Assumes symmetry | Handles asymmetry |
Method 3: IQR (Interquartile Range) - Also Robust
The IQR method uses quartiles and is familiar from boxplot “whiskers”:
# IQR method: flag if outside 1.5 * IQR (standard boxplot rule)
outliers_iqr <- pep_flag_outliers(flowering, method = "iqr", threshold = 1.5)How it works: - Q1 = 25th percentile, Q3 = 75th percentile - IQR = Q3 - Q1 - Outlier if: value < Q1 - 1.5×IQR OR value > Q3 + 1.5×IQR
Method 4: Z-Score - Sensitive but Less Robust
The z-score method assumes normally distributed data:
# Z-score method: flag if |z| > 3 (assumes normal distribution)
outliers_zscore <- pep_flag_outliers(flowering, method = "zscore", threshold = 3)Caution: The z-score method is sensitive to existing outliers. If your data has many outliers, they will inflate the standard deviation, making the z-score threshold less effective at detecting them.
Method 5: GAM Residual - Model-based, Covariate-aware
The four methods above are all univariate within a group — they compare each observation to the central tendency of its own station-species-phase combination. That misses three realistic outlier patterns:
- Year-effect outliers. In a year where most of central Europe bloomed 8 days earlier, a station that flowered on its 30-year median looks “late” only because the reference ignores the shared weather signal.
- Covariate-inconsistent reports. A lowland station that reports a DOY typical of a mountain station is invisible to a within-station detector because the station’s own history may be short.
- Trend-inconsistent reports. A station whose DOY is constant while every other station in the network has a trend is anomalous only in relation to the trend.
method = "gam_residual" tackles all three. For each
group (typically species × phase), it fits a GAM of DOY on year,
altitude, latitude, and a station random intercept (customisable via
formula =) and flags observations whose residual —
standardised by a robust MAD scale — exceeds the threshold:
# Model-based detection across the whole species x phase group.
outliers_gam <- pep_flag_outliers(
apple_flowering,
by = c("genus", "species", "phase_id"),
method = "gam_residual",
threshold = 3.5
)The deviation column now holds GAM residuals in days
(not deviations from a within-station median), and
expected_doy holds the fitted values — so the station-year
whose residual is 20 days above the fit is flagged whether or not its
station’s history made that value look “normal”.
Groups with fewer observations than min_n_per_group
(default 50) or where the GAM fails to converge silently fall back to
the 30-day rule.
Method 6: Mahalanobis - Multivariate, Robust (MCD)
All methods above flag observations one at a time. A station-year could have an entirely reasonable BBCH 60 date and an entirely reasonable BBCH 65 date, but only one day apart — biologically impossible, because first flowering precedes full flowering by ~10 days.
method = "mahalanobis" catches this by treating each
station-year as a vector across phases and computing the
robust Mahalanobis distance of that vector from the
centre of the rest of the network. “Robust” matters: the classical
Mahalanobis distance uses the sample mean and covariance, which are
themselves pulled by outliers (masking). We use the Minimum
Covariance Determinant estimator (Rousseeuw 1984,
robustbase::covMcd()) which gives a high-breakdown robust
centre and scatter, so the outliers stay outside the robust
ellipsoid instead of contaminating it:
# Multivariate detection on the phase profile per station-year.
# by = c("genus", "species") pools across phases so the covariance
# is estimated from the joint (BBCH 60, BBCH 65, BBCH 89) distribution.
outliers_mahal <- pep_flag_outliers(
apple,
by = c("genus", "species"),
method = "mahalanobis"
)The default threshold is sqrt(qchisq(0.975, df = p))
where p is the number of phases — the standard chi-square
containment cut-off under multivariate normality. Flagged station-years
get all their phase rows marked; deviation stores the
robust MD (not a day-count, so it is unitless).
Choosing a Method
| Method | Best For | Threshold Meaning |
|---|---|---|
30day |
Simple, interpretable threshold | Absolute days from median |
mad |
Most situations (recommended for single-phase) | Number of MADs from median |
iqr |
Familiar boxplot-style detection | IQR multiplier |
zscore |
Clean data, normal distribution | Standard deviations |
gam_residual |
Year / elevation / station covariates available | Robust z on GAM residuals |
mahalanobis |
Multiple phases per station-year | sqrt(chi^2_{0.975, p}) default |
Summary Statistics
summary(outliers)
#> Phenological Outlier Summary
#> =============================================
#>
#> Total observations: 1,310,094
#> Total outliers: 26573 (2.03%)
#>
#> Groups with most outliers:
#> species phase_id n_obs n_outliers pct_outliers mean_dev
#> <fctr> <int> <int> <int> <num> <num>
#> 1: Solanum tuberosum 60 114245 6420 5.62 12.0
#> 2: Prunus domestica 60 114641 2958 2.58 10.5
#> 3: Prunus domestica 65 77393 2146 2.77 10.8
#> 4: Malus domestica 60 205403 1957 0.95 9.1
#> 5: Malus domestica 65 175706 1770 1.01 9.1
#> 6: Prunus persica 60 48074 1515 3.15 11.0
#> 7: Prunus cerasus 60 136588 1227 0.90 9.0
#> 8: Prunus cerasus 65 112597 1042 0.93 9.0
#> 9: Prunus armeniaca 60 22692 943 4.16 12.1
#> 10: Pyrus communis 60 50335 940 1.87 9.5
#>
#> Deviation distribution (all data):
#> 1% 5% 25% 50% 75% 95% 99%
#> -31 -21 -8 0 8 20 30Comparing Outlier Rates Across Species
Different species may have different outlier rates due to observation difficulty or data quality:
# Check outliers for grapevine
vine_flowering <- pep[species == "Vitis vinifera" & phase_id %in% c(60, 65)]
if (nrow(vine_flowering) > 0) {
outliers_vine <- pep_flag_outliers(
pep = vine_flowering,
method = "30day",
by = c("species", "phase_id")
)
cat("Outlier comparison:\n")
cat("All flowering species: ", round(100 * mean(outliers$is_outlier), 2), "%\n")
cat("Grapevine only: ", round(100 * mean(outliers_vine$is_outlier), 2), "%\n")
} else {
cat("No grapevine flowering data available for comparison.\n")
}
#> Outlier comparison:
#> All flowering species: 2.03 %
#> Grapevine only: 3.68 %Visualizing Outliers
Visual inspection is essential for understanding outlier patterns.
The pep_plot_outliers() function provides several
visualization types — four for the univariate methods plus two targeted
at the model-based and multivariate methods:
1. Overview Plot
Shows the big picture: how many outliers, which species/phases are affected:
# Overview of outlier patterns
pep_plot_outliers(outliers, type = "overview")
#> Warning: Removed 16 rows containing non-finite outside the scale range
#> (`stat_bin()`).
What to look for: - Which species have the most outliers? (Might indicate data quality issues) - Are outliers concentrated in certain phases? (Might indicate definition problems) - What’s the overall outlier rate? (Typical: 1-5% for well-curated data)
2. Seasonal Distribution
When in the year do outliers occur?
# When do outliers occur in the year?
pep_plot_outliers(outliers, type = "seasonal")
Interpretation guide:
| Outlier Timing | Likely Explanation |
|---|---|
| Very early (DOY < 50) | Possible data errors (flowering in January/February unlikely for most species) |
| Somewhat early | Warm winters or Mediterranean locations |
| Somewhat late | Cool springs or high-altitude stations |
| Very late (DOY > 250) | Possible second flowering events - investigate! |
3. Detailed Context
See outliers alongside normal observations:
# See outliers in context of all observations
pep_plot_outliers(outliers, type = "detail", n_top = 15)
This plot shows: - Full distribution of observations
(gray) - Flagged outliers (highlighted) - The n_top
parameter controls how many species/phases to show
4. Geographic Distribution
Where are outliers located?
# Where are outliers located?
pep_plot_outliers(outliers, type = "map")
What to look for: - Clustered outliers in one region? Might indicate local data quality issues - Outliers at network edges? Might be at environmental limits - Random scatter? Suggests individual observation errors
5. Model Diagnostic (for gam_residual /
mahalanobis)
type = "diagnostic" produces a paper-ready 4-panel
figure tailored to the detection method. For gam_residual
the panels are classical GAM diagnostics — residuals vs fitted, Q-Q
against Normal, |residual| vs altitude, and a per-station
worst-case-residual map. For mahalanobis the panels switch
automatically to MD-specific diagnostics: sorted MD with the chi-square
threshold line, Q-Q of MD² against χ²_p, mean/max MD over time
(exposes network-wide anomalous epochs), and the per-station worst-case
MD map.
# For gam_residual: classical model diagnostics
outliers_gam <- pep_flag_outliers(
apple_flowering,
by = c("genus", "species", "phase_id"),
method = "gam_residual"
)
pep_plot_outliers(outliers_gam, type = "diagnostic")
GAM-residual diagnostic on a synthetic fixture with one injected contamination at station 1, year 2000. Panel A: residuals vs fitted; the injected point appears at residual ≈ +50 days. Panel B: Q-Q against N(0,1); the right-tail departure corresponds to the outlier. Panel C: |residual| vs altitude, showing the model has no altitude bias away from the outlier. Panel D: per-station worst-case residual map.
# For mahalanobis: MD-specific diagnostics
outliers_mahal <- pep_flag_outliers(
apple,
by = c("genus", "species"),
method = "mahalanobis"
)
pep_plot_outliers(outliers_mahal, type = "diagnostic")
Mahalanobis-specific diagnostic on the same synthetic fixture. Panel A: sorted robust Mahalanobis distances across all station-years, with the \(\sqrt{\chi^2_{0.975,\,3}} \approx 3.06\) threshold as a red dashed line; flagged points are red. Panel B: Q-Q of observed \(\mathrm{MD}^2\) against \(\chi^2_3\) — upper-tail departure identifies the outlier mass. Panel C: mean (blue) and max (red dotted) MD per year — the spike at year 2000 coincides with the injected anomaly. Panel D: per-station worst-case MD map.
What to look for:
- Residuals vs fitted (GAM): a level trend with no structure is ideal; a funnel or curve signals model misspecification, not outliers.
- Q-Q (GAM): heavy-tailed deviations in the corners are the real outlier signal — light tails mean the threshold is picking up noise.
- Sorted MD with threshold (Mahalanobis): the steep rise at the right-hand end is where flagged station-years sit; a long flat “knee” across the threshold means the threshold is arguably too loose or too strict.
- Mean / max MD over time (Mahalanobis): spikes in the max line identify years in which many stations reported biologically inconsistent phase sequences — useful for correlating with known data-collection or instrumental events.
6. Phase-profile Plot (primary Mahalanobis figure)
The MD detector flags station-years whose shape across
phases is atypical, even when each marginal DOY looks fine.
type = "profile" plots the phase profile as parallel
coordinates — one line per station-year, flagged ones in red, with the
robust median profile overlaid in black as a reference:
# Phase-profile parallel coordinates (primary Mahalanobis figure).
pep_plot_outliers(outliers_mahal, type = "profile")
Phase-profile parallel-coordinates plot. Each line is one station-year; the grey envelope shows the natural variability of the network, the red lines are station-years flagged by the robust Mahalanobis detector, and the thick black line is the MCD robust per-phase median profile. Red lines that cross the phases at atypical angles — not just at unusual absolute DOYs — are the classic Mahalanobis flags.
How to read it:
- Every grey line is a “normal-looking” station-year; together they form the envelope of the robust ellipsoid.
- Red lines are flagged. Look where they bend oddly: a station-year whose BBCH 60 → 65 slope is much flatter or steeper than the reference is the classic MD flag.
- The black reference profile is the robust (MCD) per-phase median — the centre of the ellipsoid. Anything strongly parallel to this is safe; anything that crosses multiple phases at unusual angles is suspect.
