diff --git a/DESCRIPTION b/DESCRIPTION
index a8d6d50b..8b8ec2b9 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,7 +2,7 @@ Type: Package
 Package: GTFShift
 Title: Explore and Analyse General Transit Feed Specification (GTFS) Files
     with a Focus on Urban Mobility
-Version: 1.0.0
+Version: 1.0.0.9000
 Authors@R: c(
     person(c("Gonçalo", "F."), "Matos", , "goncaloafmatos@tecnico.pt", role = c("aut", "cre"),
            comment = c(ORCID = "0009-0001-3489-1732")),
diff --git a/NAMESPACE b/NAMESPACE
index 1cb466b9..038bb21e 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -12,6 +12,7 @@ export(get_network_extension)
 export(get_prioritisation_stats)
 export(get_route_frequency_hourly)
 export(get_stop_frequency_hourly)
+export(get_trip_speed_profile)
 export(get_way_frequency_hourly)
 export(load_feed)
 export(multiline_to_sorted_linestring)
@@ -55,7 +56,12 @@ importFrom(jsonlite,
 importFrom(purrr,map_dfr)
 importFrom(rlang,.data)
 importFrom(stats,
+  IQR,
+  median,
+  quantile,
+  sd,
   setNames,
+  var,
   weighted.mean
 )
 importFrom(stringr,
diff --git a/NEWS.md b/NEWS.md
new file mode 100644
index 00000000..349d2b1a
--- /dev/null
+++ b/NEWS.md
@@ -0,0 +1,12 @@
+# GTFShift 1.0.0.900
+
+(Dev version, not released on CRAN yet)
+
+- add `get_trip_speed_profile` to compute commercial speed and other speed related metrics at trip level
+- extended `rt_extend_prioritisation` output with `speed_p15` and `speed_p85`
+
+
+# GTFShift 1.0.0
+
+- Stable version
+- CRAN release
\ No newline at end of file
diff --git a/R/get_trip_speed_profile.R b/R/get_trip_speed_profile.R
new file mode 100644
index 00000000..11350b97
--- /dev/null
+++ b/R/get_trip_speed_profile.R
@@ -0,0 +1,318 @@
+#' Get trip speed profile from GTFS-RT speed estimates
+#'
+#' Computes aggregated speed metrics for trips or groups from the output of
+#' \code{GTFShift::rt_average_speed()}. Aggregates by trip, route, and day
+#' (customizable via the \code{by} parameter), calculating commercial speed
+#' (distance between first and last updates along geometry divided by elapsed time),
+#' alternative commercial speed (considering the 2nd and penultimate updates to avoid
+#' terminal wait biases), and measures of centrality and spread for speed observations.
+#'
+#' @param rt_speed data.frame or sf data.frame. The result of
+#'   \code{GTFShift::rt_average_speed()}. Must contain at least the columns
+#'   \code{distance_along_geometry} and \code{timestamp}.
+#' @param by Character vector (Default \code{c("trip_id", "route_id", "day")}).
+#'   Columns to aggregate by. If \code{"day"} is included in \code{by} but not present
+#'   in \code{rt_speed}, it is automatically derived from the \code{timestamp} column.
+#'   Set to \code{NULL} or \code{character(0)} to compute metrics across the entire dataset.
+#' @param dist_col Character (Default \code{"distance_along_geometry"}). Column name
+#'   present in \code{rt_speed} representing cumulative distance along geometry (in meters).
+#' @param speed_col Character (Default \code{"speed_kmh"}). Column name present in
+#'   \code{rt_speed} representing estimated speed between consecutive updates (in km/h).
+#' @param time_col Character (Default \code{"timestamp"}). Column name present in
+#'   \code{rt_speed} representing update timestamps (numeric epoch, POSIXct, or Date).
+#'
+#' @details
+#' For each group defined by \code{by} (by default, each unique combination of
+#' \code{trip_id}, \code{route_id}, and \code{day}), observations are ordered
+#' chronologically by \code{timestamp}.
+#'
+#' Let \eqn{\{(d_i, t_i)\}_{i=1}^n} denote the ordered sequence of updates, where
+#' \eqn{d_i} is the distance along geometry (meters) and \eqn{t_i} is the timestamp (seconds).
+#'
+#' \strong{Commercial speed} is calculated as the total distance traveled between the
+#' first and last updates divided by the elapsed time:
+#' \deqn{v_{\mathrm{commercial}} = \frac{|d_n - d_1|}{1000} \div \frac{t_n - t_1}{3600}}
+#' If \eqn{n < 2} or \eqn{t_n \le t_1}, \code{commercial_speed} is \code{NA}.
+#'
+#' \strong{Alternative commercial speed} (\code{commercial_speed_alt}) uses the 2nd and
+#' penultimate (\eqn{n-1}) observations to eliminate potential dwell times or layovers at
+#' the terminal stops:
+#' \deqn{v_{\mathrm{commercial\_alt}} = \frac{|d_{n-1} - d_2|}{1000} \div \frac{t_{n-1} - t_2}{3600}}
+#' If \eqn{n < 4} or \eqn{t_{n-1} \le t_2}, \code{commercial_speed_alt} is \code{NA}.
+#'
+#' If \code{by} does not include \code{trip_id} (e.g., aggregating at route or day level),
+#' \code{commercial_speed} and \code{commercial_speed_alt} are computed per trip and averaged
+#' across trips within each group.
+#'
+#' \strong{Measures of centrality and spread} are calculated from all valid (non-NA, finite)
+#' speed observations in \code{speed_kmh} for each group:
+#' \describe{
+#'   \item{timestamp_min}{Earliest timestamp in the trip/group.}
+#'   \item{timestamp_max}{Latest timestamp in the trip/group.}
+#'   \item{commercial_speed}{Commercial speed between first and last update (km/h).}
+#'   \item{commercial_speed_alt}{Alternative commercial speed between 2nd and penultimate update (km/h).}
+#'   \item{speed_avg}{Arithmetic mean speed (km/h).}
+#'   \item{speed_median}{Median speed (km/h).}
+#'   \item{speed_sd}{Standard deviation of speeds (km/h).}
+#'   \item{speed_var}{Variance of speeds (km/h)^2.}
+#'   \item{speed_min}{Minimum observed speed (km/h).}
+#'   \item{speed_max}{Maximum observed speed (km/h).}
+#'   \item{speed_p15}{15th percentile speed (km/h).}
+#'   \item{speed_p25}{25th percentile (1st quartile) speed (km/h).}
+#'   \item{speed_p75}{75th percentile (3rd quartile) speed (km/h).}
+#'   \item{speed_p85}{85th percentile speed (km/h).}
+#'   \item{speed_iqr}{Interquartile range (speed_p75 - speed_p25) (km/h).}
+#'   \item{speed_count}{Count of valid speed observations.}
+#'   \item{n_updates}{Total number of position updates in the group.}
+#' }
+#'
+#' In addition, an attribute \code{"global_summary"} is attached to the returned data frame,
+#' providing the same metrics evaluated over the entire input dataset.
+#'
+#' @returns data.frame. Aggregated speed profile metrics with one row per group.
+#'
+#' @examples
+#' \donttest{
+#' # Get GTFS-RT data collection and route geometries
+#' rt_collect_file <- system.file(
+#'   "extdata/samples", "gtfs_rt_sample_tcb_4_4-CS-TERM.csv", package = "GTFShift"
+#' )
+#' rt_collection <- read.csv(rt_collect_file) |>
+#'   sf::st_as_sf(coords = c("longitude", "latitude"), crs = 4326) |>
+#'   dplyr::select(-speed)
+#'
+#' osm_routes <- sf::st_read(
+#'   system.file("extdata/samples", "osm_routes_tcb.gpkg", package = "GTFShift"),
+#'   quiet = TRUE
+#' ) |>
+#'   dplyr::filter(route_id %in% rt_collection$route_id) |>
+#'   dplyr::mutate(geom = GTFShift::multiline_to_sorted_linestring(geom, metric_crs = 3763))
+#'
+#' # Compute speeds with rt_average_speed()
+#' speed <- GTFShift::rt_average_speed(
+#'   rt_collection = rt_collection,
+#'   trips_geometries = osm_routes,
+#'   rt_collection_trips_geometries_match_col = "route_id",
+#'   metric_crs = 3763
+#' )
+#'
+#' # Compute trip speed profile
+#' profile <- GTFShift::get_trip_speed_profile(speed)
+#' head(profile)
+#'
+#' # Customize aggregation (e.g. by route and day)
+#' route_profile <- GTFShift::get_trip_speed_profile(speed, by = c("route_id", "day"))
+#' head(route_profile)
+#' }
+#'
+#' @seealso \code{GTFShift::rt_average_speed()}
+#'
+#' @import dplyr
+#' @importFrom stats median sd var quantile IQR
+#' @importFrom rlang .data
+#'
+#' @export
+get_trip_speed_profile <- function(
+  rt_speed,
+  by = c("trip_id", "route_id", "day"),
+  dist_col = "distance_along_geometry",
+  speed_col = "speed_kmh",
+  time_col = "timestamp"
+) {
+  # 0. Input validations
+  if (!is.data.frame(rt_speed)) {
+    stop("rt_speed must be a data.frame or sf object")
+  }
+
+  # Drop sf geometry if present for aggregation performance
+  if (inherits(rt_speed, "sf")) {
+    rt_speed <- sf::st_drop_geometry(rt_speed)
+  }
+  rt_speed <- dplyr::ungroup(rt_speed)
+
+  # Check required columns
+  required_cols <- c(dist_col, time_col)
+  missing_cols <- setdiff(required_cols, colnames(rt_speed))
+  if (length(missing_cols) > 0) {
+    stop(paste("rt_speed is missing required column(s):", paste(missing_cols, collapse = ", ")))
+  }
+
+  # 1. Handle "day" column derivation if included in 'by' but absent in rt_speed
+  if (!is.null(by) && "day" %in% by && !"day" %in% colnames(rt_speed)) {
+    if (is.numeric(rt_speed[[time_col]])) {
+      rt_speed$day <- as.Date(as.POSIXct(rt_speed[[time_col]], origin = "1970-01-01", tz = "UTC"))
+    } else {
+      rt_speed$day <- as.Date(rt_speed[[time_col]])
+    }
+  }
+
+  # 2. Validate / resolve grouping columns
+  if (!is.null(by) && length(by) > 0) {
+    if (missing(by)) {
+      # When default 'by' is used, keep only columns that exist in rt_speed
+      by <- intersect(by, colnames(rt_speed))
+      if (length(by) == 0) {
+        stop("None of the default grouping columns ('trip_id', 'route_id', 'day') were found in rt_speed.")
+      }
+    } else {
+      missing_by <- setdiff(by, colnames(rt_speed))
+      if (length(missing_by) > 0) {
+        stop(paste("The following grouping columns specified in 'by' were not found in rt_speed:",
+                   paste(missing_by, collapse = ", ")))
+      }
+    }
+  }
+
+  # Helper: compute time difference in seconds between two timestamps
+  get_time_diff_sec <- function(t_start, t_end) {
+    if (is.na(t_start) || is.na(t_end)) return(NA_real_)
+    if (inherits(t_start, "POSIXt") || inherits(t_start, "Date")) {
+      return(as.numeric(difftime(t_end, t_start, units = "secs")))
+    }
+    if (is.numeric(t_start) && is.numeric(t_end)) {
+      return(as.numeric(t_end - t_start))
+    }
+    t1 <- suppressWarnings(as.POSIXct(t_start))
+    t2 <- suppressWarnings(as.POSIXct(t_end))
+    if (!is.na(t1) && !is.na(t2)) {
+      return(as.numeric(difftime(t2, t1, units = "secs")))
+    }
+    return(as.numeric(t_end - t_start))
+  }
+
+  # Helper: compute commercial speed in km/h
+  compute_speed_kmh <- function(dist_m, time_sec) {
+    if (is.na(dist_m) || is.na(time_sec) || time_sec <= 0) {
+      return(NA_real_)
+    }
+    round((dist_m / 1000) / (time_sec / 3600), 2)
+  }
+
+  # Helper: compute commercial speeds for a single sequence of updates
+  compute_single_trip_commercial_speeds <- function(trip_df) {
+    trip_df <- trip_df[order(trip_df[[time_col]]), , drop = FALSE]
+    n_trip <- nrow(trip_df)
+
+    # First and last updates
+    if (n_trip >= 2) {
+      dt <- get_time_diff_sec(trip_df[[time_col]][1], trip_df[[time_col]][n_trip])
+      dd <- abs(trip_df[[dist_col]][n_trip] - trip_df[[dist_col]][1])
+      comm_speed <- compute_speed_kmh(dd, dt)
+    } else {
+      comm_speed <- NA_real_
+    }
+
+    # 2nd and penultimate updates
+    if (n_trip >= 4) {
+      dt_alt <- get_time_diff_sec(trip_df[[time_col]][2], trip_df[[time_col]][n_trip - 1])
+      dd_alt <- abs(trip_df[[dist_col]][n_trip - 1] - trip_df[[dist_col]][2])
+      comm_speed_alt <- compute_speed_kmh(dd_alt, dt_alt)
+    } else {
+      comm_speed_alt <- NA_real_
+    }
+
+    list(commercial_speed = comm_speed, commercial_speed_alt = comm_speed_alt)
+  }
+
+  # Helper: compute full profile for a group
+  compute_profile <- function(group_df) {
+    group_df <- group_df[order(group_df[[time_col]]), , drop = FALSE]
+    n <- nrow(group_df)
+
+    # Time boundaries to locate trip/group in time
+    timestamp_min <- if (n > 0) group_df[[time_col]][1] else NA
+    timestamp_max <- if (n > 0) group_df[[time_col]][n] else NA
+
+    # Commercial speed calculations
+    if ("trip_id" %in% colnames(group_df) && length(unique(group_df$trip_id)) > 1) {
+      # If group contains multiple trips (e.g. grouped by route_id/day without trip_id),
+      # compute commercial speeds per trip and take average
+      trip_list <- split(group_df, group_df$trip_id)
+      trip_res <- lapply(trip_list, compute_single_trip_commercial_speeds)
+      comm_vec <- vapply(trip_res, function(x) x$commercial_speed, numeric(1))
+      comm_alt_vec <- vapply(trip_res, function(x) x$commercial_speed_alt, numeric(1))
+
+      commercial_speed <- if (any(!is.na(comm_vec))) round(mean(comm_vec, na.rm = TRUE), 2) else NA_real_
+      commercial_speed_alt <- if (any(!is.na(comm_alt_vec))) round(mean(comm_alt_vec, na.rm = TRUE), 2) else NA_real_
+    } else {
+      # Single trip within group (standard case)
+      trip_res <- compute_single_trip_commercial_speeds(group_df)
+      commercial_speed <- trip_res$commercial_speed
+      commercial_speed_alt <- trip_res$commercial_speed_alt
+    }
+
+    # Measures of centrality and spread for speed observations
+    if (speed_col %in% colnames(group_df)) {
+      raw_speeds <- group_df[[speed_col]]
+      valid_speeds <- raw_speeds[!is.na(raw_speeds) & is.finite(raw_speeds)]
+    } else {
+      valid_speeds <- numeric(0)
+    }
+
+    if (length(valid_speeds) > 0) {
+      speed_avg <- round(mean(valid_speeds), 2)
+      speed_median <- round(stats::median(valid_speeds), 2)
+      speed_sd <- if (length(valid_speeds) > 1) round(stats::sd(valid_speeds), 2) else NA_real_
+      speed_var <- if (length(valid_speeds) > 1) round(stats::var(valid_speeds), 2) else NA_real_
+      speed_min <- round(min(valid_speeds), 2)
+      speed_max <- round(max(valid_speeds), 2)
+      speed_p15 <- round(as.numeric(stats::quantile(valid_speeds, probs = 0.15, names = FALSE)), 2)
+      speed_p25 <- round(as.numeric(stats::quantile(valid_speeds, probs = 0.25, names = FALSE)), 2)
+      speed_p75 <- round(as.numeric(stats::quantile(valid_speeds, probs = 0.75, names = FALSE)), 2)
+      speed_p85 <- round(as.numeric(stats::quantile(valid_speeds, probs = 0.85, names = FALSE)), 2)
+      speed_iqr <- round(stats::IQR(valid_speeds), 2)
+      speed_count <- length(valid_speeds)
+    } else {
+      speed_avg <- NA_real_
+      speed_median <- NA_real_
+      speed_sd <- NA_real_
+      speed_var <- NA_real_
+      speed_min <- NA_real_
+      speed_max <- NA_real_
+      speed_p15 <- NA_real_
+      speed_p25 <- NA_real_
+      speed_p75 <- NA_real_
+      speed_p85 <- NA_real_
+      speed_iqr <- NA_real_
+      speed_count <- 0L
+    }
+
+    data.frame(
+      timestamp_min = timestamp_min,
+      timestamp_max = timestamp_max,
+      commercial_speed = commercial_speed,
+      commercial_speed_alt = commercial_speed_alt,
+      speed_avg = speed_avg,
+      speed_median = speed_median,
+      speed_sd = speed_sd,
+      speed_var = speed_var,
+      speed_min = speed_min,
+      speed_max = speed_max,
+      speed_p15 = speed_p15,
+      speed_p25 = speed_p25,
+      speed_p75 = speed_p75,
+      speed_p85 = speed_p85,
+      speed_iqr = speed_iqr,
+      speed_count = speed_count,
+      n_updates = n,
+      stringsAsFactors = FALSE
+    )
+  }
+
+  # 3. Compute speed profile
+  if (is.null(by) || length(by) == 0) {
+    result <- compute_profile(rt_speed)
+  } else {
+    result <- rt_speed |>
+      dplyr::group_by(dplyr::across(dplyr::all_of(by))) |>
+      dplyr::group_modify(function(group_df, keys) {
+        compute_profile(group_df)
+      }) |>
+      dplyr::ungroup()
+  }
+
+  # Attach global summary across the entire dataset as attribute
+  attr(result, "global_summary") <- compute_profile(rt_speed)
+
+  return(result)
+}
diff --git a/R/rt_extend_prioritisation.R b/R/rt_extend_prioritisation.R
index 6c67d4b4..c50446c9 100644
--- a/R/rt_extend_prioritisation.R
+++ b/R/rt_extend_prioritisation.R
@@ -20,8 +20,10 @@
 #' \describe{
 #'   \item{speed_avg}{The average speed of the vehicles on the way.}
 #'   \item{speed_median}{The median speed of the vehicles on the way.}
+#'   \item{speed_p15}{The 15th percentile speed of the vehicles on the way.}
 #'   \item{speed_p25}{The 25th percentile speed of the vehicles on the way.}
 #'   \item{speed_p75}{The 75th percentile speed of the vehicles on the way.}
+#'   \item{speed_p85}{The 85th percentile speed of the vehicles on the way.}
 #'   \item{speed_count}{The number of speed observations on the way.}
 #' }
 #'
@@ -172,8 +174,10 @@ rt_extend_prioritisation <- function(
       dplyr::summarise(
         speed_avg = mean(.data[[rt_attr_speed]], na.rm = TRUE),
         speed_median = stats::median(.data[[rt_attr_speed]], na.rm = TRUE),
+        speed_p15 = stats::quantile(.data[[rt_attr_speed]], probs = 0.15, na.rm = TRUE),
         speed_p25 = stats::quantile(.data[[rt_attr_speed]], probs = 0.25, na.rm = TRUE),
         speed_p75 = stats::quantile(.data[[rt_attr_speed]], probs = 0.75, na.rm = TRUE),
+        speed_p85 = stats::quantile(.data[[rt_attr_speed]], probs = 0.85, na.rm = TRUE),
         speed_count = dplyr::n()
       ) |>
       dplyr::ungroup())
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 2a98f6df..6b9f9cda 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -8,7 +8,7 @@ template:
 
 navbar:
   structure:
-   left:  [intro, reference, articles, dashboard]
+   left:  [intro, reference, articles, news, dashboard]
 
   components:
     dashboard:
diff --git a/inst/extdata/samples/osmextract_lisbon_highways_sample.gpkg b/inst/extdata/samples/osmextract_lisbon_highways_sample.gpkg
new file mode 100644
index 00000000..c9c15286
Binary files /dev/null and b/inst/extdata/samples/osmextract_lisbon_highways_sample.gpkg differ
diff --git a/inst/extdata/samples/osmextract_tcb_network.gpkg b/inst/extdata/samples/osmextract_tcb_network.gpkg
new file mode 100644
index 00000000..e43dd40a
Binary files /dev/null and b/inst/extdata/samples/osmextract_tcb_network.gpkg differ
diff --git a/man/get_trip_speed_profile.Rd b/man/get_trip_speed_profile.Rd
new file mode 100644
index 00000000..dcba7192
--- /dev/null
+++ b/man/get_trip_speed_profile.Rd
@@ -0,0 +1,130 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/get_trip_speed_profile.R
+\name{get_trip_speed_profile}
+\alias{get_trip_speed_profile}
+\title{Get trip speed profile from GTFS-RT speed estimates}
+\usage{
+get_trip_speed_profile(
+  rt_speed,
+  by = c("trip_id", "route_id", "day"),
+  dist_col = "distance_along_geometry",
+  speed_col = "speed_kmh",
+  time_col = "timestamp"
+)
+}
+\arguments{
+\item{rt_speed}{data.frame or sf data.frame. The result of
+\code{GTFShift::rt_average_speed()}. Must contain at least the columns
+\code{distance_along_geometry} and \code{timestamp}.}
+
+\item{by}{Character vector (Default \code{c("trip_id", "route_id", "day")}).
+Columns to aggregate by. If \code{"day"} is included in \code{by} but not present
+in \code{rt_speed}, it is automatically derived from the \code{timestamp} column.
+Set to \code{NULL} or \code{character(0)} to compute metrics across the entire dataset.}
+
+\item{dist_col}{Character (Default \code{"distance_along_geometry"}). Column name
+present in \code{rt_speed} representing cumulative distance along geometry (in meters).}
+
+\item{speed_col}{Character (Default \code{"speed_kmh"}). Column name present in
+\code{rt_speed} representing estimated speed between consecutive updates (in km/h).}
+
+\item{time_col}{Character (Default \code{"timestamp"}). Column name present in
+\code{rt_speed} representing update timestamps (numeric epoch, POSIXct, or Date).}
+}
+\value{
+data.frame. Aggregated speed profile metrics with one row per group.
+}
+\description{
+Computes aggregated speed metrics for trips or groups from the output of
+\code{GTFShift::rt_average_speed()}. Aggregates by trip, route, and day
+(customizable via the \code{by} parameter), calculating commercial speed
+(distance between first and last updates along geometry divided by elapsed time),
+alternative commercial speed (considering the 2nd and penultimate updates to avoid
+terminal wait biases), and measures of centrality and spread for speed observations.
+}
+\details{
+For each group defined by \code{by} (by default, each unique combination of
+\code{trip_id}, \code{route_id}, and \code{day}), observations are ordered
+chronologically by \code{timestamp}.
+
+Let \eqn{\{(d_i, t_i)\}_{i=1}^n} denote the ordered sequence of updates, where
+\eqn{d_i} is the distance along geometry (meters) and \eqn{t_i} is the timestamp (seconds).
+
+\strong{Commercial speed} is calculated as the total distance traveled between the
+first and last updates divided by the elapsed time:
+\deqn{v_{\mathrm{commercial}} = \frac{|d_n - d_1|}{1000} \div \frac{t_n - t_1}{3600}}
+If \eqn{n < 2} or \eqn{t_n \le t_1}, \code{commercial_speed} is \code{NA}.
+
+\strong{Alternative commercial speed} (\code{commercial_speed_alt}) uses the 2nd and
+penultimate (\eqn{n-1}) observations to eliminate potential dwell times or layovers at
+the terminal stops:
+\deqn{v_{\mathrm{commercial\_alt}} = \frac{|d_{n-1} - d_2|}{1000} \div \frac{t_{n-1} - t_2}{3600}}
+If \eqn{n < 4} or \eqn{t_{n-1} \le t_2}, \code{commercial_speed_alt} is \code{NA}.
+
+If \code{by} does not include \code{trip_id} (e.g., aggregating at route or day level),
+\code{commercial_speed} and \code{commercial_speed_alt} are computed per trip and averaged
+across trips within each group.
+
+\strong{Measures of centrality and spread} are calculated from all valid (non-NA, finite)
+speed observations in \code{speed_kmh} for each group:
+\describe{
+  \item{timestamp_min}{Earliest timestamp in the trip/group.}
+  \item{timestamp_max}{Latest timestamp in the trip/group.}
+  \item{commercial_speed}{Commercial speed between first and last update (km/h).}
+  \item{commercial_speed_alt}{Alternative commercial speed between 2nd and penultimate update (km/h).}
+  \item{speed_avg}{Arithmetic mean speed (km/h).}
+  \item{speed_median}{Median speed (km/h).}
+  \item{speed_sd}{Standard deviation of speeds (km/h).}
+  \item{speed_var}{Variance of speeds (km/h)^2.}
+  \item{speed_min}{Minimum observed speed (km/h).}
+  \item{speed_max}{Maximum observed speed (km/h).}
+  \item{speed_p15}{15th percentile speed (km/h).}
+  \item{speed_p25}{25th percentile (1st quartile) speed (km/h).}
+  \item{speed_p75}{75th percentile (3rd quartile) speed (km/h).}
+  \item{speed_p85}{85th percentile speed (km/h).}
+  \item{speed_iqr}{Interquartile range (speed_p75 - speed_p25) (km/h).}
+  \item{speed_count}{Count of valid speed observations.}
+  \item{n_updates}{Total number of position updates in the group.}
+}
+
+In addition, an attribute \code{"global_summary"} is attached to the returned data frame,
+providing the same metrics evaluated over the entire input dataset.
+}
+\examples{
+\donttest{
+# Get GTFS-RT data collection and route geometries
+rt_collect_file <- system.file(
+  "extdata/samples", "gtfs_rt_sample_tcb_4_4-CS-TERM.csv", package = "GTFShift"
+)
+rt_collection <- read.csv(rt_collect_file) |>
+  sf::st_as_sf(coords = c("longitude", "latitude"), crs = 4326) |>
+  dplyr::select(-speed)
+
+osm_routes <- sf::st_read(
+  system.file("extdata/samples", "osm_routes_tcb.gpkg", package = "GTFShift"),
+  quiet = TRUE
+) |>
+  dplyr::filter(route_id \%in\% rt_collection$route_id) |>
+  dplyr::mutate(geom = GTFShift::multiline_to_sorted_linestring(geom, metric_crs = 3763))
+
+# Compute speeds with rt_average_speed()
+speed <- GTFShift::rt_average_speed(
+  rt_collection = rt_collection,
+  trips_geometries = osm_routes,
+  rt_collection_trips_geometries_match_col = "route_id",
+  metric_crs = 3763
+)
+
+# Compute trip speed profile
+profile <- GTFShift::get_trip_speed_profile(speed)
+head(profile)
+
+# Customize aggregation (e.g. by route and day)
+route_profile <- GTFShift::get_trip_speed_profile(speed, by = c("route_id", "day"))
+head(route_profile)
+}
+
+}
+\seealso{
+\code{GTFShift::rt_average_speed()}
+}
diff --git a/man/rt_collect_protobuf.Rd b/man/rt_collect_protobuf.Rd
index f0e41543..bab46fb3 100644
--- a/man/rt_collect_protobuf.Rd
+++ b/man/rt_collect_protobuf.Rd
@@ -61,4 +61,5 @@ head(
     dplyr::select("vehicle.trip.trip_id", "vehicle.position.latitude", "vehicle.position.longitude")
 )
 }
+
 }
diff --git a/man/rt_extend_prioritisation.Rd b/man/rt_extend_prioritisation.Rd
index a1508ff3..3df5ac0b 100644
--- a/man/rt_extend_prioritisation.Rd
+++ b/man/rt_extend_prioritisation.Rd
@@ -28,8 +28,10 @@ sf data.frame. Extended lane prioritisation with the following columns:
 \describe{
   \item{speed_avg}{The average speed of the vehicles on the way.}
   \item{speed_median}{The median speed of the vehicles on the way.}
+  \item{speed_p15}{The 15th percentile speed of the vehicles on the way.}
   \item{speed_p25}{The 25th percentile speed of the vehicles on the way.}
   \item{speed_p75}{The 75th percentile speed of the vehicles on the way.}
+  \item{speed_p85}{The 85th percentile speed of the vehicles on the way.}
   \item{speed_count}{The number of speed observations on the way.}
 }
 }
diff --git a/tests/testthat/test-get_trip_speed_profile.R b/tests/testthat/test-get_trip_speed_profile.R
new file mode 100644
index 00000000..9ad84c91
--- /dev/null
+++ b/tests/testthat/test-get_trip_speed_profile.R
@@ -0,0 +1,493 @@
+library(testthat)
+library(sf)
+
+test_that("get_trip_speed_profile calculates speed profile with default parameters", {
+    # 2 trips, same route, 5 updates each
+    # Trip T1: 1500m over 240s -> commercial_speed = (1.5 km) / (240 / 3600 h) = 22.5 km/h
+    # Trip T1 alt: (1200 - 300)m over (1180 - 1060)s = 900m / 120s -> 27.0 km/h
+    # Trip T2: 1800m over 200s -> commercial_speed = (1.8 km) / (200 / 3600 h) = 32.4 km/h
+    # Trip T2 alt: (1300 - 400)m over (2150 - 2050)s = 900m / 100s -> 32.4 km/h
+    df <- data.frame(
+        trip_id = rep(c("T1", "T2"), each = 5),
+        route_id = "R1",
+        timestamp = c(1000, 1060, 1120, 1180, 1240, 2000, 2050, 2100, 2150, 2200),
+        distance_along_geometry = c(0, 300, 700, 1200, 1500, 0, 400, 800, 1300, 1800),
+        speed_kmh = c(NA, 18, 24, 30, 18, NA, 28.8, 28.8, 36, 36)
+    )
+
+    res <- GTFShift::get_trip_speed_profile(df)
+
+    expect_s3_class(res, "data.frame")
+    expect_equal(nrow(res), 2)
+    expected_cols <- c(
+        "trip_id", "route_id", "day",
+        "timestamp_min", "timestamp_max",
+        "commercial_speed", "commercial_speed_alt",
+        "speed_avg", "speed_median", "speed_sd", "speed_var",
+        "speed_min", "speed_max",
+        "speed_p15", "speed_p25", "speed_p75", "speed_p85",
+        "speed_iqr", "speed_count", "n_updates"
+    )
+    expect_contains(names(res), expected_cols)
+
+    # Trip T1 checks
+    r_t1 <- res[res$trip_id == "T1", ]
+    expect_equal(r_t1$route_id, "R1")
+    expect_equal(r_t1$day, as.Date("1970-01-01"))
+    expect_equal(r_t1$timestamp_min, 1000)
+    expect_equal(r_t1$timestamp_max, 1240)
+    expect_equal(r_t1$commercial_speed, 22.5)
+    expect_equal(r_t1$commercial_speed_alt, 27.0)
+    # T1 valid speeds: c(18, 24, 30, 18)
+    expect_equal(r_t1$speed_avg, 22.5)
+    expect_equal(r_t1$speed_median, 21.0)
+    expect_equal(r_t1$speed_min, 18.0)
+    expect_equal(r_t1$speed_max, 30.0)
+    expect_equal(r_t1$speed_p15, 18.0)
+    expect_equal(r_t1$speed_p25, 18.0)
+    expect_equal(r_t1$speed_p75, 25.5)
+    expect_equal(r_t1$speed_p85, 27.3)
+    expect_equal(r_t1$speed_iqr, 7.5)
+    expect_equal(r_t1$speed_count, 4)
+    expect_equal(r_t1$n_updates, 5)
+
+    # Monotonicity of percentiles
+    expect_true(r_t1$speed_p15 <= r_t1$speed_p25)
+    expect_true(r_t1$speed_p25 <= r_t1$speed_median)
+    expect_true(r_t1$speed_median <= r_t1$speed_p75)
+    expect_true(r_t1$speed_p75 <= r_t1$speed_p85)
+
+    # Trip T2 checks
+    r_t2 <- res[res$trip_id == "T2", ]
+    expect_equal(r_t2$timestamp_min, 2000)
+    expect_equal(r_t2$timestamp_max, 2200)
+    expect_equal(r_t2$commercial_speed, 32.4)
+    expect_equal(r_t2$commercial_speed_alt, 32.4)
+    # T2 valid speeds: c(28.8, 28.8, 36, 36)
+    expect_equal(r_t2$speed_avg, 32.4)
+    expect_equal(r_t2$speed_median, 32.4)
+    expect_equal(r_t2$speed_min, 28.8)
+    expect_equal(r_t2$speed_max, 36.0)
+    expect_equal(r_t2$speed_count, 4)
+    expect_equal(r_t2$n_updates, 5)
+
+    # Check global_summary attribute
+    global_sum <- attr(res, "global_summary")
+    expect_s3_class(global_sum, "data.frame")
+    expect_equal(nrow(global_sum), 1)
+    expect_equal(global_sum$timestamp_min, 1000)
+    expect_equal(global_sum$timestamp_max, 2200)
+    expect_equal(global_sum$commercial_speed, 27.45) # mean of 22.5 and 32.4
+    expect_equal(global_sum$commercial_speed_alt, 29.7) # mean of 27.0 and 32.4
+    expect_equal(global_sum$speed_count, 8)
+    expect_equal(global_sum$n_updates, 10)
+})
+
+test_that("get_trip_speed_profile supports custom aggregation columns in 'by'", {
+    df <- data.frame(
+        trip_id = rep(c("T1", "T2"), each = 5),
+        route_id = "R1",
+        timestamp = c(1000, 1060, 1120, 1180, 1240, 2000, 2050, 2100, 2150, 2200),
+        distance_along_geometry = c(0, 300, 700, 1200, 1500, 0, 400, 800, 1300, 1800),
+        speed_kmh = c(NA, 18, 24, 30, 18, NA, 28.8, 28.8, 36, 36)
+    )
+
+    # 1. Grouping by route_id and day (collapsing multiple trips into one group)
+    res_route <- GTFShift::get_trip_speed_profile(df, by = c("route_id", "day"))
+    expect_equal(nrow(res_route), 1)
+    expect_equal(res_route$route_id, "R1")
+    # Commercial speeds should be averaged across the two trips
+    expect_equal(res_route$commercial_speed, 27.45)
+    expect_equal(res_route$commercial_speed_alt, 29.7)
+    expect_equal(res_route$n_updates, 10)
+    expect_equal(res_route$speed_count, 8)
+
+    # 2. Grouping solely by trip_id
+    res_trip <- GTFShift::get_trip_speed_profile(df, by = "trip_id")
+    expect_equal(nrow(res_trip), 2)
+    expect_false("route_id" %in% names(res_trip))
+    expect_false("day" %in% names(res_trip))
+    expect_equal(res_trip$commercial_speed, c(22.5, 32.4))
+
+    # 3. by = NULL (whole dataset evaluated without grouping)
+    res_null <- GTFShift::get_trip_speed_profile(df, by = NULL)
+    expect_equal(nrow(res_null), 1)
+    expect_equal(res_null$commercial_speed, 27.45)
+    expect_equal(res_null$n_updates, 10)
+
+    # 4. by = character(0)
+    res_char0 <- GTFShift::get_trip_speed_profile(df, by = character(0))
+    expect_equal(nrow(res_char0), 1)
+    expect_equal(res_char0$commercial_speed, 27.45)
+
+    # 5. Default by when route_id is missing: keeps available columns (trip_id and derived day)
+    df_no_route <- df[, c("trip_id", "timestamp", "distance_along_geometry", "speed_kmh")]
+    res_default_partial <- GTFShift::get_trip_speed_profile(df_no_route)
+    expect_equal(nrow(res_default_partial), 2)
+    expect_contains(names(res_default_partial), c("trip_id", "day"))
+    expect_false("route_id" %in% names(res_default_partial))
+})
+
+test_that("get_trip_speed_profile derives day column from different timestamp types", {
+    # 1. POSIXct timestamp
+    t_posix <- as.POSIXct(c("2024-05-01 10:00:00", "2024-05-01 10:02:00", "2024-05-01 10:04:00", "2024-05-01 10:06:00"), tz = "UTC")
+    df_posix <- data.frame(
+        trip_id = "T1",
+        route_id = "R1",
+        timestamp = t_posix,
+        distance_along_geometry = c(0, 500, 1000, 1500),
+        speed_kmh = c(NA, 15, 15, 15)
+    )
+    res_posix <- GTFShift::get_trip_speed_profile(df_posix)
+    expect_equal(res_posix$day, as.Date("2024-05-01"))
+    expect_equal(res_posix$commercial_speed, 15.0)
+    expect_equal(res_posix$commercial_speed_alt, 15.0)
+
+    # 2. Date timestamp
+    df_date <- data.frame(
+        trip_id = "T1",
+        route_id = "R1",
+        timestamp = as.Date(c("2024-06-01", "2024-06-02", "2024-06-03", "2024-06-04")),
+        distance_along_geometry = c(0, 100000, 200000, 300000),
+        speed_kmh = c(NA, 4.17, 4.17, 4.17)
+    )
+    res_date <- GTFShift::get_trip_speed_profile(df_date)
+    expect_equal(res_date$day, as.Date(c("2024-06-01", "2024-06-02", "2024-06-03", "2024-06-04")))
+
+    res_date_trip <- GTFShift::get_trip_speed_profile(df_date, by = "trip_id")
+    expect_equal(res_date_trip$commercial_speed, 4.17)
+
+    # 3. Pre-existing day column is preserved
+    df_custom_day <- df_posix
+    df_custom_day$day <- as.Date("2020-01-01")
+    res_custom_day <- GTFShift::get_trip_speed_profile(df_custom_day)
+    expect_equal(res_custom_day$day, as.Date("2020-01-01"))
+})
+
+test_that("get_trip_speed_profile drops geometry from sf objects", {
+    df_sf <- sf::st_sf(
+        trip_id = rep("T1", 5),
+        route_id = "R1",
+        timestamp = c(1000, 1060, 1120, 1180, 1240),
+        distance_along_geometry = c(0, 300, 700, 1200, 1500),
+        speed_kmh = c(NA, 18, 24, 30, 18),
+        geometry = sf::st_sfc(lapply(1:5, function(i) sf::st_point(c(i, i))))
+    )
+
+    res_sf <- GTFShift::get_trip_speed_profile(df_sf)
+    expect_false(inherits(res_sf, "sf"))
+    expect_false("geometry" %in% names(res_sf))
+    expect_equal(res_sf$commercial_speed, 22.5)
+})
+
+test_that("get_trip_speed_profile handles update count edge cases (n < 2, 2 <= n < 4, n >= 4)", {
+    # 1 update: both commercial speeds NA, variance/sd NA
+    df_n1 <- data.frame(
+        trip_id = "T1",
+        timestamp = 1000,
+        distance_along_geometry = 0,
+        speed_kmh = 25
+    )
+    res_n1 <- GTFShift::get_trip_speed_profile(df_n1, by = "trip_id")
+    expect_true(is.na(res_n1$commercial_speed))
+    expect_true(is.na(res_n1$commercial_speed_alt))
+    expect_true(is.na(res_n1$speed_sd))
+    expect_true(is.na(res_n1$speed_var))
+    expect_equal(res_n1$speed_avg, 25.0)
+    expect_equal(res_n1$speed_count, 1)
+    expect_equal(res_n1$n_updates, 1)
+
+    # 2 updates: commercial_speed valid, commercial_speed_alt NA
+    df_n2 <- data.frame(
+        trip_id = "T1",
+        timestamp = c(1000, 1060),
+        distance_along_geometry = c(0, 500),
+        speed_kmh = c(NA, 30)
+    )
+    res_n2 <- GTFShift::get_trip_speed_profile(df_n2, by = "trip_id")
+    expect_equal(res_n2$commercial_speed, 30.0)
+    expect_true(is.na(res_n2$commercial_speed_alt))
+
+    # 3 updates: commercial_speed valid, commercial_speed_alt NA
+    df_n3 <- data.frame(
+        trip_id = "T1",
+        timestamp = c(1000, 1060, 1120),
+        distance_along_geometry = c(0, 500, 1000),
+        speed_kmh = c(NA, 30, 30)
+    )
+    res_n3 <- GTFShift::get_trip_speed_profile(df_n3, by = "trip_id")
+    expect_equal(res_n3$commercial_speed, 30.0)
+    expect_true(is.na(res_n3$commercial_speed_alt))
+
+    # 4 updates: both valid
+    df_n4 <- data.frame(
+        trip_id = "T1",
+        timestamp = c(1000, 1060, 1120, 1180),
+        distance_along_geometry = c(0, 500, 1000, 1500),
+        speed_kmh = c(NA, 30, 30, 30)
+    )
+    res_n4 <- GTFShift::get_trip_speed_profile(df_n4, by = "trip_id")
+    expect_equal(res_n4$commercial_speed, 30.0)
+    expect_equal(res_n4$commercial_speed_alt, 30.0)
+
+    # Multi-trip group where all trips have n < 2
+    df_multi_n1 <- data.frame(
+        route_id = "R1",
+        trip_id = c("T1", "T2"),
+        timestamp = c(1000, 2000),
+        distance_along_geometry = c(0, 0),
+        speed_kmh = c(20, 30)
+    )
+    res_multi_n1 <- GTFShift::get_trip_speed_profile(df_multi_n1, by = "route_id")
+    expect_true(is.na(res_multi_n1$commercial_speed))
+    expect_true(is.na(res_multi_n1$commercial_speed_alt))
+
+    # Multi-trip group where all trips have n < 4
+    df_multi_n2 <- data.frame(
+        route_id = "R1",
+        trip_id = c("T1", "T1", "T2", "T2"),
+        timestamp = c(1000, 1060, 2000, 2060),
+        distance_along_geometry = c(0, 500, 0, 600),
+        speed_kmh = c(NA, 30, NA, 36)
+    )
+    res_multi_n2 <- GTFShift::get_trip_speed_profile(df_multi_n2, by = "route_id")
+    expect_equal(res_multi_n2$commercial_speed, 33.0) # mean of 30 and 36
+    expect_true(is.na(res_multi_n2$commercial_speed_alt))
+
+    # Multi-trip group where one trip has 4 updates and one has 1 update
+    df_multi_mixed <- data.frame(
+        route_id = "R1",
+        trip_id = c("T1", "T1", "T1", "T1", "T2"),
+        timestamp = c(1000, 1060, 1120, 1180, 2000),
+        distance_along_geometry = c(0, 500, 1000, 1500, 0),
+        speed_kmh = c(NA, 30, 30, 30, 20)
+    )
+    res_multi_mixed <- GTFShift::get_trip_speed_profile(df_multi_mixed, by = "route_id")
+    expect_equal(res_multi_mixed$commercial_speed, 30.0) # only T1 has valid commercial speed
+    expect_equal(res_multi_mixed$commercial_speed_alt, 30.0)
+})
+
+test_that("get_trip_speed_profile handles unordered updates, zero/negative elapsed time, and missing values", {
+    # 1. Unordered updates are sorted chronologically
+    df_sorted <- data.frame(
+        trip_id = "T1",
+        timestamp = c(1000, 1060, 1120, 1180),
+        distance_along_geometry = c(0, 500, 1000, 1500),
+        speed_kmh = c(NA, 30, 30, 30)
+    )
+    df_scrambled <- df_sorted[c(3, 1, 4, 2), ]
+
+    res_sorted <- GTFShift::get_trip_speed_profile(df_sorted, by = "trip_id")
+    res_scrambled <- GTFShift::get_trip_speed_profile(df_scrambled, by = "trip_id")
+    expect_equal(res_scrambled$timestamp_min, res_sorted$timestamp_min)
+    expect_equal(res_scrambled$timestamp_max, res_sorted$timestamp_max)
+    expect_equal(res_scrambled$commercial_speed, res_sorted$commercial_speed)
+    expect_equal(res_scrambled$commercial_speed_alt, res_sorted$commercial_speed_alt)
+
+    # 2. Zero or negative time difference yields NA commercial speed
+    df_zero_time <- data.frame(
+        trip_id = "T1",
+        timestamp = c(1000, 1000),
+        distance_along_geometry = c(0, 500),
+        speed_kmh = c(NA, 30)
+    )
+    res_zero <- GTFShift::get_trip_speed_profile(df_zero_time, by = "trip_id")
+    expect_true(is.na(res_zero$commercial_speed))
+
+    # 3. NA timestamp or NA distance yields NA commercial speed
+    df_na_time <- data.frame(
+        trip_id = "T1",
+        timestamp = c(NA, 1060),
+        distance_along_geometry = c(0, 500),
+        speed_kmh = c(NA, 30)
+    )
+    res_na_time <- GTFShift::get_trip_speed_profile(df_na_time, by = "trip_id")
+    expect_true(is.na(res_na_time$commercial_speed))
+
+    df_na_dist <- data.frame(
+        trip_id = "T1",
+        timestamp = c(1000, 1060),
+        distance_along_geometry = c(NA, 500),
+        speed_kmh = c(NA, 30)
+    )
+    res_na_dist <- GTFShift::get_trip_speed_profile(df_na_dist, by = "trip_id")
+    expect_true(is.na(res_na_dist$commercial_speed))
+
+    # 4. Backward distance progression uses absolute difference
+    df_backward <- data.frame(
+        trip_id = "T1",
+        timestamp = c(1000, 1060, 1120, 1180),
+        distance_along_geometry = c(1500, 1000, 500, 0),
+        speed_kmh = c(NA, 30, 30, 30)
+    )
+    res_backward <- GTFShift::get_trip_speed_profile(df_backward, by = "trip_id")
+    expect_equal(res_backward$commercial_speed, 30.0)
+    expect_equal(res_backward$commercial_speed_alt, 30.0)
+
+    # 5. Character timestamps parseable as POSIXct
+    df_char <- data.frame(
+        trip_id = "T1",
+        timestamp = c("2024-05-01 10:00:00", "2024-05-01 10:02:00", "2024-05-01 10:04:00", "2024-05-01 10:06:00"),
+        distance_along_geometry = c(0, 500, 1000, 1500),
+        speed_kmh = c(NA, 15, 15, 15)
+    )
+    res_char <- GTFShift::get_trip_speed_profile(df_char, by = "trip_id")
+    expect_equal(res_char$commercial_speed, 15.0)
+    expect_equal(res_char$commercial_speed_alt, 15.0)
+})
+
+test_that("get_trip_speed_profile handles speed column variations and edge cases", {
+    # 1. Non-finite speed values (Inf, -Inf) are ignored
+    df_inf <- data.frame(
+        trip_id = "T1",
+        timestamp = c(1000, 1060, 1120, 1180),
+        distance_along_geometry = c(0, 500, 1000, 1500),
+        speed_kmh = c(20, Inf, -Inf, 40)
+    )
+    res_inf <- GTFShift::get_trip_speed_profile(df_inf, by = "trip_id")
+    expect_equal(res_inf$speed_avg, 30.0)
+    expect_equal(res_inf$speed_min, 20.0)
+    expect_equal(res_inf$speed_max, 40.0)
+    expect_equal(res_inf$speed_count, 2)
+
+    # 2. All speeds NA
+    df_all_na <- data.frame(
+        trip_id = "T1",
+        timestamp = c(1000, 1060),
+        distance_along_geometry = c(0, 500),
+        speed_kmh = c(NA_real_, NA_real_)
+    )
+    res_all_na <- GTFShift::get_trip_speed_profile(df_all_na, by = "trip_id")
+    expect_equal(res_all_na$speed_count, 0)
+    expect_true(is.na(res_all_na$speed_avg))
+    expect_true(is.na(res_all_na$speed_median))
+    expect_true(is.na(res_all_na$speed_sd))
+    expect_true(is.na(res_all_na$speed_var))
+    expect_true(is.na(res_all_na$speed_min))
+    expect_true(is.na(res_all_na$speed_max))
+    expect_true(is.na(res_all_na$speed_p15))
+    expect_true(is.na(res_all_na$speed_p25))
+    expect_true(is.na(res_all_na$speed_p75))
+    expect_true(is.na(res_all_na$speed_p85))
+    expect_true(is.na(res_all_na$speed_iqr))
+    expect_equal(res_all_na$commercial_speed, 30.0)
+
+    # 3. Missing speed_col completely
+    df_no_speed <- data.frame(
+        trip_id = "T1",
+        timestamp = c(1000, 1060),
+        distance_along_geometry = c(0, 500)
+    )
+    res_no_speed <- GTFShift::get_trip_speed_profile(df_no_speed, by = "trip_id")
+    expect_equal(res_no_speed$speed_count, 0)
+    expect_true(is.na(res_no_speed$speed_avg))
+    expect_equal(res_no_speed$commercial_speed, 30.0)
+
+    # 4. Custom column names
+    df_custom <- data.frame(
+        my_trip = "T1",
+        my_time = c(1000, 1060, 1120, 1180),
+        my_dist = c(0, 500, 1000, 1500),
+        my_spd = c(NA, 30, 30, 30)
+    )
+    res_custom <- GTFShift::get_trip_speed_profile(
+        df_custom,
+        by = "my_trip",
+        dist_col = "my_dist",
+        speed_col = "my_spd",
+        time_col = "my_time"
+    )
+    expect_equal(res_custom$commercial_speed, 30.0)
+    expect_equal(res_custom$speed_avg, 30.0)
+    expect_equal(res_custom$speed_count, 3)
+
+    # 5. Grouped data.frame input is ungrouped properly
+    df_grouped <- dplyr::group_by(df_inf, trip_id)
+    res_grouped <- GTFShift::get_trip_speed_profile(df_grouped, by = "trip_id")
+    expect_equal(res_grouped$commercial_speed, 30.0)
+})
+
+test_that("get_trip_speed_profile stops with informative error on invalid inputs", {
+    # 1. Non-dataframe input
+    expect_error(
+        GTFShift::get_trip_speed_profile(123),
+        "rt_speed must be a data.frame or sf object"
+    )
+    expect_error(
+        GTFShift::get_trip_speed_profile(list(a = 1)),
+        "rt_speed must be a data.frame or sf object"
+    )
+
+    # 2. Missing required columns
+    expect_error(
+        GTFShift::get_trip_speed_profile(data.frame(x = 1)),
+        "rt_speed is missing required column\\(s\\): distance_along_geometry, timestamp"
+    )
+    expect_error(
+        GTFShift::get_trip_speed_profile(
+            data.frame(x = 1),
+            dist_col = "my_dist",
+            time_col = "my_time"
+        ),
+        "rt_speed is missing required column\\(s\\): my_dist, my_time"
+    )
+
+    # 3. Specified grouping column in by is missing
+    df_valid <- data.frame(
+        trip_id = "T1",
+        timestamp = 1000,
+        distance_along_geometry = 0
+    )
+    expect_error(
+        GTFShift::get_trip_speed_profile(df_valid, by = c("non_existent_col")),
+        "The following grouping columns specified in 'by' were not found in rt_speed: non_existent_col"
+    )
+})
+
+test_that("get_trip_speed_profile integrates seamlessly with rt_average_speed output", {
+    trip_geom <- sf::st_sf(
+        trip_id = c("T1", "T2"),
+        geometry = sf::st_sfc(
+            sf::st_linestring(matrix(c(0, 0, 1000, 0), ncol = 2, byrow = TRUE)),
+            sf::st_linestring(matrix(c(0, 100, 1000, 100), ncol = 2, byrow = TRUE)),
+            crs = 3857
+        )
+    )
+
+    rt_updates <- sf::st_sf(
+        trip_id = c("T1", "T1", "T1", "T1", "T2", "T2", "T2", "T2"),
+        route_id = "R1",
+        timestamp = c(1000, 1060, 1120, 1180, 2000, 2060, 2120, 2180),
+        geometry = sf::st_sfc(
+            sf::st_point(c(0, 0)),
+            sf::st_point(c(250, 0)),
+            sf::st_point(c(600, 0)),
+            sf::st_point(c(900, 0)),
+            sf::st_point(c(0, 100)),
+            sf::st_point(c(300, 100)),
+            sf::st_point(c(650, 100)),
+            sf::st_point(c(1000, 100)),
+            crs = 3857
+        )
+    )
+
+    speed_sf <- GTFShift::rt_average_speed(rt_updates, trip_geom, metric_crs = 3857)
+    profile <- GTFShift::get_trip_speed_profile(speed_sf)
+
+    expect_s3_class(profile, "data.frame")
+    expect_false(inherits(profile, "sf"))
+    expect_equal(nrow(profile), 2)
+    expect_equal(profile$trip_id, c("T1", "T2"))
+    expect_equal(profile$route_id, c("R1", "R1"))
+    expect_true(all(!is.na(profile$commercial_speed)))
+    expect_true(all(!is.na(profile$commercial_speed_alt)))
+    expect_true(all(!is.na(profile$speed_avg)))
+    expect_true(all(!is.na(profile$speed_p15)))
+    expect_true(all(!is.na(profile$speed_p85)))
+
+    # Global summary check
+    global_summary <- attr(profile, "global_summary")
+    expect_s3_class(global_summary, "data.frame")
+    expect_equal(nrow(global_summary), 1)
+    expect_equal(global_summary$n_updates, 8)
+})
diff --git a/tests/testthat/test-rt_extend_prioritisation.R b/tests/testthat/test-rt_extend_prioritisation.R
index 270f069d..14e5edb4 100644
--- a/tests/testthat/test-rt_extend_prioritisation.R
+++ b/tests/testthat/test-rt_extend_prioritisation.R
@@ -25,20 +25,80 @@ test_that("rt_extend_prioritisation extends lane prioritisation with speed metri
 
     res <- GTFShift::rt_extend_prioritisation(lanes_sf, rt_points, metric_crs = 3857)
     expect_s3_class(res, "sf")
-    expect_contains(names(res), c("speed_avg", "speed_median", "speed_p25", "speed_p75", "speed_count"))
+    expect_contains(names(res), c("speed_avg", "speed_median", "speed_p15", "speed_p25", "speed_p75", "speed_p85", "speed_count"))
 
     # For w1: speeds 10 and 20 (point 3 at 80 is STOPPED_AT so filtered out)
     expect_equal(res$speed_avg[res$way_osm_id == "w1"], 15)
     expect_equal(res$speed_median[res$way_osm_id == "w1"], 15)
+    expect_equal(unname(res$speed_p15[res$way_osm_id == "w1"]), 11.5)
     expect_equal(unname(res$speed_p25[res$way_osm_id == "w1"]), 12.5)
     expect_equal(unname(res$speed_p75[res$way_osm_id == "w1"]), 17.5)
+    expect_equal(unname(res$speed_p85[res$way_osm_id == "w1"]), 18.5)
     expect_equal(res$speed_count[res$way_osm_id == "w1"], 2)
 
     # For w2: speed 40
     expect_equal(res$speed_avg[res$way_osm_id == "w2"], 40)
+    expect_equal(res$speed_median[res$way_osm_id == "w2"], 40)
+    expect_equal(unname(res$speed_p15[res$way_osm_id == "w2"]), 40)
+    expect_equal(unname(res$speed_p25[res$way_osm_id == "w2"]), 40)
+    expect_equal(unname(res$speed_p75[res$way_osm_id == "w2"]), 40)
+    expect_equal(unname(res$speed_p85[res$way_osm_id == "w2"]), 40)
     expect_equal(res$speed_count[res$way_osm_id == "w2"], 1)
 })
 
+test_that("rt_extend_prioritisation accurately computes speed_p15 and speed_p85 with multiple observations and NA values", {
+    lanes_sf <- st_sf(
+        way_osm_id = c("w1", "w2", "w_empty"),
+        geometry = st_sfc(
+            st_linestring(matrix(c(-9.150, 38.700, -9.140, 38.700), ncol = 2, byrow = TRUE)),
+            st_linestring(matrix(c(-9.130, 38.700, -9.120, 38.700), ncol = 2, byrow = TRUE)),
+            st_linestring(matrix(c(-9.100, 38.700, -9.090, 38.700), ncol = 2, byrow = TRUE)),
+            crs = 4326
+        )
+    )
+
+    # w1: 10 points with speeds 10, 20, ..., 100
+    w1_lons <- seq(-9.149, -9.141, length.out = 10)
+    w1_points <- lapply(w1_lons, function(x) st_point(c(x, 38.700)))
+    w1_speeds <- seq(10, 100, by = 10)
+
+    # w2: 4 points, one with NA speed: speeds 20, 40, NA, 60
+    w2_lons <- seq(-9.128, -9.122, length.out = 4)
+    w2_points <- lapply(w2_lons, function(x) st_point(c(x, 38.700)))
+    w2_speeds <- c(20, 40, NA, 60)
+
+    rt_points <- st_sf(
+        speed = c(w1_speeds, w2_speeds),
+        current_status = "IN_TRANSIT_TO",
+        geometry = st_sfc(c(w1_points, w2_points), crs = 4326)
+    )
+
+    res <- GTFShift::rt_extend_prioritisation(lanes_sf, rt_points, metric_crs = 3763)
+
+    # w1 percentiles: 15% is 23.5, 25% is 32.5, 75% is 77.5, 85% is 86.5
+    expect_equal(unname(res$speed_p15[res$way_osm_id == "w1"]), 23.5)
+    expect_equal(unname(res$speed_p25[res$way_osm_id == "w1"]), 32.5)
+    expect_equal(unname(res$speed_p75[res$way_osm_id == "w1"]), 77.5)
+    expect_equal(unname(res$speed_p85[res$way_osm_id == "w1"]), 86.5)
+    expect_equal(res$speed_count[res$way_osm_id == "w1"], 10)
+
+    # Verify monotonic percentile relationships: p15 <= p25 <= median <= p75 <= p85
+    expect_true(res$speed_p15[res$way_osm_id == "w1"] <= res$speed_p25[res$way_osm_id == "w1"])
+    expect_true(res$speed_p25[res$way_osm_id == "w1"] <= res$speed_median[res$way_osm_id == "w1"])
+    expect_true(res$speed_median[res$way_osm_id == "w1"] <= res$speed_p75[res$way_osm_id == "w1"])
+    expect_true(res$speed_p75[res$way_osm_id == "w1"] <= res$speed_p85[res$way_osm_id == "w1"])
+
+    # w2 percentiles with NA: speeds c(20, 40, 60) -> 15% is 26, 85% is 54
+    expect_equal(unname(res$speed_p15[res$way_osm_id == "w2"]), 26.0)
+    expect_equal(unname(res$speed_p85[res$way_osm_id == "w2"]), 54.0)
+    expect_equal(res$speed_count[res$way_osm_id == "w2"], 4)
+
+    # w_empty: no points captured -> speed_p15 and speed_p85 are NA
+    expect_true(is.na(res$speed_p15[res$way_osm_id == "w_empty"]))
+    expect_true(is.na(res$speed_p85[res$way_osm_id == "w_empty"]))
+    expect_true(is.na(res$speed_count[res$way_osm_id == "w_empty"]))
+})
+
 test_that("rt_extend_prioritisation raises warning when metric_crs is default", {
     lanes_sf <- st_sf(
         way_osm_id = "w1",