1. Data Ingestion & Configuration

We configure file paths relative to the R script location (located one level up from the data/ directory). New test CSV files can be added directly to the file_paths vector.

# Define CSV data files (relative to R/ directory)
file_paths <- c(
  "data/Another-Sample-Aristotle-Report-universal-oat.csv",
  "data/Sample-Aristotle-Report-universal-oat.csv",
  "data/Still-Another-Sample-Aristotle-Report-universal-oat.csv",
  "data/Theriome-Aristotle-Plus-6-2026-universal-oat.csv"
)

# Function to read and clean a single Theriome CSV export
read_theriome_csv <- function(path) {
  # Verify file exists
  if (!file.exists(path)) {
    warning(paste("File not found:", path))
    return(NULL)
  }
  
  # Read CSV skipping first 30 metadata rows (row 31 is column header)
  # Analyte data spans rows 31 through 327 (header and 296 analytes)
  raw_df <- read_csv(path, skip = 30, show_col_types = FALSE)
  
  # Filter to analyte rows (using the category column for robustness)
  analyte_df <- raw_df %>%
    filter(category == "Analytes") %>%
    select(
      canonicalAnalyteId,
      sourceAnalyteName,
      metric,
      value,
      refLow,
      refHigh
    ) %>%
    mutate(
      source_file = basename(path),
      z_score = metric / 100.0,  # Convert metric (Z * 100) to actual Z-score
      value = as.numeric(value),
      refHigh = as.numeric(refHigh),
      refLow = as.numeric(refLow)
    )
  
  return(analyte_df)
}

# Ingest and pool all available test files
raw_combined <- map_dfr(file_paths, read_theriome_csv)

# Resolve name inconsistencies using Canonical ID.
# Count frequencies of each name per canonical ID.
name_freqs <- raw_combined %>%
  count(canonicalAnalyteId, sourceAnalyteName) %>%
  filter(!is.na(canonicalAnalyteId))

# Determine the resolved name for each canonical ID, handling overrides
resolved_names <- name_freqs %>%
  group_by(canonicalAnalyteId) %>%
  arrange(desc(n), sourceAnalyteName, .by_group = TRUE) %>%
  slice(1) %>%
  ungroup() %>%
  select(canonicalAnalyteId, resolved_name = sourceAnalyteName) %>%
  mutate(
    resolved_name = case_when(
      canonicalAnalyteId == 7165 ~ "Oxidized Glutahione",
      canonicalAnalyteId == 7176 ~ "Prenenolone sulfate",
      TRUE ~ resolved_name
    )
  )

# Replace sourceAnalyteName with resolved_name in raw_combined
raw_combined <- raw_combined %>%
  left_join(resolved_names, by = "canonicalAnalyteId") %>%
  mutate(sourceAnalyteName = ifelse(!is.na(resolved_name), resolved_name, sourceAnalyteName)) %>%
  select(-resolved_name)

Below is the successfully loaded test reports:

Summary of Loaded Test Reports
Source File Observations Min Value Max Value
Another-Sample-Aristotle-Report-universal-oat.csv 296 0.00 17366942
Sample-Aristotle-Report-universal-oat.csv 296 0.00 43433943
Still-Another-Sample-Aristotle-Report-universal-oat.csv 296 2.44 18985247
Theriome-Aristotle-Plus-6-2026-universal-oat.csv 296 0.00 9091720

A total of 4 test reports were successfully ingested, pooling 1184 total observations across 296 distinct analytes. Parameters

2. Mathematical Derivation of Reference Parameters

Even though population means are omitted from the export CSVs, we mathematically recover both the population mean (\(\mu\)) and standard deviation (\(\sigma\)) for each test using two complementary methods:

  1. Single-Test Algebraic Solver: Since \(Z = \frac{X - \mu}{\sigma}\) and \(X_{\text{high}} = \mu + 3\sigma\), we solve for \(\sigma\) and \(\mu\) directly:

\[\sigma_{\text{derived}} = \frac{X_{\text{high}} - X}{3 - Z} \quad \text{and} \quad \mu_{\text{derived}} = X_{\text{high}} - 3\sigma_{\text{derived}}\]

  1. Theoretical \(Z\)-Score Floor (\(Z_{\text{min}}\)): Calculated as \(Z_{\text{min}} = -\frac{\mu}{\sigma} = -\frac{1}{\text{CV}}\).
# Calculate derived parameters per observation
df_processed <- raw_combined %>%
  mutate(
    # Single-test algebraic recovery of sigma and mu
    sigma_derived = ifelse(abs(3 - z_score) > 1e-4, (refHigh - value) / (3 - z_score), NA_real_),
    mu_derived = refHigh - (3 * sigma_derived),
    cv_derived = ifelse(mu_derived > 0, sigma_derived / mu_derived, NA_real_),
    z_floor_derived = ifelse(!is.na(cv_derived) & cv_derived > 0, -1 / cv_derived, NA_real_)
  )

Derived Reference Parameters Preview

Below is a preview of the derived standard deviation (\(\sigma_{\text{derived}}\)), mean (\(\mu_{\text{derived}}\)), coefficient of variation (CV), and theoretical lower Z-score floor (\(Z_{\text{min}}\)) for the first 10 observations:

Derived Reference Parameters Sample
Analyte Name Raw Value Ref High Z-Score Derived Std Dev (σ) Derived Mean (μ) Derived CV Z Floor (Z_min)
4-Imidazoleacetic acid 839 82246.2 -0.65 22303.342 15336.173 1.454 -0.688
Acetylcholine 15145 97154.9 -0.10 26454.806 17790.481 1.487 -0.672
Cystamine 55 1323.0 -0.39 374.041 200.876 1.862 -0.537
Dopamine 128 19914.0 -0.41 5802.346 2506.962 2.314 -0.432
Epinephrine 361 9274.0 -0.45 2583.478 1523.565 1.696 -0.590
Histamine 402 176996.4 -0.32 53191.084 17423.147 3.053 -0.328
Imidazole 2296 70427.0 -0.79 17976.517 16497.449 1.090 -0.918
Metanephrine 1641 340947.1 -0.26 104081.626 28702.223 3.626 -0.276
Methylguanidine 27070 71451.2 0.61 18569.540 15742.581 1.180 -0.848
Methylhistamine 2973 16526.2 -0.56 3807.079 5104.964 0.746 -1.341

For each observation, the underlying population mean (\(\mu\)) and standard deviation (\(\sigma\)) are algebraically derived using the upper limit of the reference range (\(X_{\text{high}}\)) and the reported Z-score. This allows us to reconstruct the coefficient of variation (CV) and compute the theoretical minimum possible Z-score floor (\(Z_{\text{min}} = -1/\text{CV}\)) for every analyte.

3. Regression Fits & Model Comparison (\(R^2_{\text{raw}}\) vs \(R^2_{\text{log}}\))

For each analyte with multiple test observations, we evaluate whether \(Z\) scales linearly with raw concentration \(X\) or logarithmic concentration \(\ln(X)\).

# Group by analyte and perform model fitting
analyte_summary <- df_processed %>%
  group_by(canonicalAnalyteId, sourceAnalyteName) %>%
  summarise(
    n_obs = n(),
    mean_value = mean(value, na.rm = TRUE),
    mean_z = mean(z_score, na.rm = TRUE),
    avg_refHigh = mean(refHigh, na.rm = TRUE),
    avg_mu_derived = mean(mu_derived, na.rm = TRUE),
    avg_sigma_derived = mean(sigma_derived, na.rm = TRUE),
    avg_cv = mean(cv_derived, na.rm = TRUE),
    avg_z_floor = mean(z_floor_derived, na.rm = TRUE),
    
    # Regression Fits (Raw vs Log) - robustly handled to avoid NA errors
    fit_raw_r2 = tryCatch({
      df_temp <- data.frame(z = z_score, val = value) %>% 
        filter(!is.na(z), !is.na(val))
      if (nrow(df_temp) >= 2 && n_distinct(df_temp$val) >= 2) {
        summary(lm(z ~ val, data = df_temp))$r.squared
      } else NA_real_
    }, error = function(e) NA_real_),
    
    fit_log_r2 = tryCatch({
      df_temp <- data.frame(z = z_score, val = value) %>% 
        filter(!is.na(z), !is.na(val), val > 0)
      if (nrow(df_temp) >= 2 && n_distinct(df_temp$val) >= 2) {
        summary(lm(z ~ log(val), data = df_temp))$r.squared
      } else NA_real_
    }, error = function(e) NA_real_),
    
    .groups = "drop"
  ) %>%
  mutate(
    raw_fit_dominates = ifelse(!is.na(fit_raw_r2) & !is.na(fit_log_r2), fit_raw_r2 >= fit_log_r2, TRUE)
  )

Regression Model Fits and Analysis

We evaluate both raw concentration (\(X\)) and logarithmic concentration (\(\ln(X)\)) models for each analyte. If an analyte’s Z-score fits a linear model of raw concentration better than log concentration, raw_fit_dominates is set to TRUE.

Below is the summary of the model fits:

Comparison of Raw vs Logarithmic Regression Model Fits
Total Evaluated Analytes Raw Fit Better or Equal Log Fit Better Average R² (Raw) Average R² (Log)
296 217 79 0.8097 0.7582

Top Analytes Where Logarithmic Fit Dominates

The table below lists the top 10 analytes where a logarithmic concentration model fits the reported Z-score significantly better than a raw concentration model:

Top Log-dominant Analytes
Analyte Name R² Raw Model R² Log Model R² Improvement
2-Deoxyadenosine 0.2482 0.4458 0.1977
5-Hydroxytryptophan 0.2247 0.4128 0.1881
Pantothenic acid 0.2495 0.4368 0.1873
Betaine 0.1547 0.3284 0.1738
Phosphocreatine 0.4592 0.6323 0.1731
Glucose 1,6-bisphosphate 0.1804 0.3489 0.1685
Lauric acid 0.0011 0.1530 0.1519
Naproxen 0.4155 0.5569 0.1415
Capric acid 0.1637 0.3023 0.1385
Indole-3-lactic acid 0.4461 0.5806 0.1346

4. Analyte Explorer

This interactive visualization tool allows you to examine individual analytes in detail.

How to Use

  1. Select an Analyte from the dropdown menu in the upper-left of the plot to see its specific distribution and test results.
  2. Interpret the Subplots:
    • Left Subplot (Linear Scale) displays the Raw Concentration values of the tests, the derived Raw model normal distribution (blue), and the fitted Log model log-normal distribution (orange).
    • Right Subplot (Log Scale) displays the Log-transformed Concentration values, the raw model distribution mapped to log space (blue), and the log model normal distribution (orange).
  3. Analyze the Markers:
    • Test Results are marked as grey dashed vertical lines.
    • Inferred Mean (\(\mu\)) is represented by a solid green vertical line.
    • Reference Limit (\(X_{\text{high}}\)) is marked as a solid red vertical line.

This has been disabled because it is extremely resource intensive.

library(plotly)

# Unique analytes sorted alphabetically by name
analytes_list <- unique(df_processed %>% select(canonicalAnalyteId, sourceAnalyteName)) %>%
  arrange(sourceAnalyteName)

# Total number of analytes
M <- nrow(analytes_list)

p_explorer <- plot_ly()

# Loop through each analyte and add its traces
for (i in 1:M) {
  # Get data for this analyte
  id <- analytes_list$canonicalAnalyteId[i]
  name <- analytes_list$sourceAnalyteName[i]
  obs <- df_processed %>% filter(canonicalAnalyteId == id)
  
  ref_high <- mean(obs$refHigh, na.rm = TRUE)
  if (is.na(ref_high) || ref_high <= 0) ref_high <- 100 # safe default
  
  mu_raw <- mean(obs$mu_derived, na.rm = TRUE)
  sigma_raw <- mean(obs$sigma_derived, na.rm = TRUE)
  
  # Fallback for raw parameters if they cannot be derived
  if (is.na(mu_raw) || is.na(sigma_raw) || sigma_raw <= 0) {
    mu_raw <- mean(obs$value, na.rm = TRUE)
    sigma_raw <- sd(obs$value, na.rm = TRUE)
    if (is.na(sigma_raw) || sigma_raw <= 0) {
      sigma_raw <- ref_high / 3
      if (is.na(mu_raw)) mu_raw <- ref_high / 2
    }
  }
  
  # Fit log model for log-normal parameters
  df_fit <- obs %>% filter(!is.na(value), !is.na(z_score), value > 0)
  fit_log_ok <- FALSE
  if (nrow(df_fit) >= 2 && n_distinct(df_fit$value) >= 2) {
    fit <- tryCatch(lm(z_score ~ log(value), data = df_fit), error = function(e) NULL)
    if (!is.null(fit)) {
      coefs <- coef(fit)
      beta <- coefs[1]
      alpha <- coefs[2]
      if (!is.na(alpha) && !is.na(beta) && abs(alpha) > 1e-5) {
        mu_log <- -beta / alpha
        sigma_log <- 1 / abs(alpha)
        fit_log_ok <- TRUE
      }
    }
  }
  
  # Fallback for log parameters
  if (!fit_log_ok) {
    if (mu_raw > 0) {
      mu_log <- log(mu_raw^2 / sqrt(mu_raw^2 + sigma_raw^2))
      sigma_log <- sqrt(log(1 + (sigma_raw / mu_raw)^2))
    } else {
      mu_log <- log(ref_high / 2)
      sigma_log <- 1
    }
  }
  
  # ----------------------------------------------------
  # 1. Left Plot (Linear scale)
  # ----------------------------------------------------
  max_val <- max(ref_high, max(obs$value, na.rm = TRUE), mu_raw + 3 * sigma_raw, na.rm = TRUE)
  x_grid_lin <- seq(0, max_val * 1.2, length.out = 100)
  
  y_raw_lin <- dnorm(x_grid_lin, mean = mu_raw, sd = sigma_raw)
  y_log_lin <- dlnorm(x_grid_lin, meanlog = mu_log, sdlog = sigma_log)
  
  max_y_lin <- max(y_raw_lin, y_log_lin, na.rm = TRUE)
  if (is.na(max_y_lin) || max_y_lin <= 0) max_y_lin <- 1.0
  
  # Trace 1: Raw density (linear)
  p_explorer <- p_explorer %>% add_trace(
    x = x_grid_lin, y = y_raw_lin, type = "scatter", mode = "lines",
    line = list(color = "#3182ce", width = 2), name = "Raw Model Density",
    xaxis = "x", yaxis = "y", visible = (i == 1),
    showlegend = (i == 1), legendgroup = "raw_dens",
    hoverinfo = "text", text = paste("Raw Model Density:", round(y_raw_lin, 5))
  )
  
  # Trace 2: Log density (linear)
  p_explorer <- p_explorer %>% add_trace(
    x = x_grid_lin, y = y_log_lin, type = "scatter", mode = "lines",
    line = list(color = "#dd6b20", width = 2), name = "Log Model Density",
    xaxis = "x", yaxis = "y", visible = (i == 1),
    showlegend = (i == 1), legendgroup = "log_dens",
    hoverinfo = "text", text = paste("Log Model Density:", round(y_log_lin, 5))
  )
  
  # Trace 3: Test results (linear vertical lines)
  x_tests_lin <- rep(obs$value, each = 3) * c(1, 1, NA)
  y_tests_lin <- rep(c(0, max_y_lin * 1.05, NA), nrow(obs))
  p_explorer <- p_explorer %>% add_trace(
    x = x_tests_lin, y = y_tests_lin, type = "scatter", mode = "lines",
    line = list(color = "#718096", dash = "dash", width = 1), name = "Test Results",
    xaxis = "x", yaxis = "y", visible = (i == 1),
    showlegend = (i == 1), legendgroup = "test_res",
    hoverinfo = "none"
  )
  
  # Trace 4: Inferred mean (linear)
  p_explorer <- p_explorer %>% add_trace(
    x = c(mu_raw, mu_raw), y = c(0, max_y_lin * 1.05), type = "scatter", mode = "lines",
    line = list(color = "#2f855a", width = 2), name = "Inferred Mean (\u03bc)",
    xaxis = "x", yaxis = "y", visible = (i == 1),
    showlegend = (i == 1), legendgroup = "inf_mean",
    hoverinfo = "text", text = paste("Inferred Mean (\u03bc):", round(mu_raw, 2))
  )
  
  # Trace 5: Reference limit (linear)
  p_explorer <- p_explorer %>% add_trace(
    x = c(ref_high, ref_high), y = c(0, max_y_lin * 1.05), type = "scatter", mode = "lines",
    line = list(color = "#e53e3e", width = 2), name = "Reference Limit (X_high)",
    xaxis = "x", yaxis = "y", visible = (i == 1),
    showlegend = (i == 1), legendgroup = "ref_limit",
    hoverinfo = "text", text = paste("Reference Limit:", round(ref_high, 2))
  )
  
  # ----------------------------------------------------
  # 2. Right Plot (Log scale)
  # ----------------------------------------------------
  max_val_log <- log(max_val + 0.1)
  y_grid_log <- seq(-1.5, max_val_log * 1.2, length.out = 100)
  
  # Density values in log space
  dens_log_log <- dnorm(y_grid_log, mean = mu_log, sd = sigma_log)
  
  # Raw model transformed to log scale
  exp_y <- exp(y_grid_log)
  dens_raw_log <- dnorm(exp_y, mean = mu_raw, sd = sigma_raw) * exp_y
  
  max_y_log <- max(dens_log_log, dens_raw_log, na.rm = TRUE)
  if (is.na(max_y_log) || max_y_log <= 0) max_y_log <- 1.0
  
  # Trace 6: Raw density (log space)
  p_explorer <- p_explorer %>% add_trace(
    x = y_grid_log, y = dens_raw_log, type = "scatter", mode = "lines",
    line = list(color = "#3182ce", width = 2), name = "Raw Model Density",
    xaxis = "x2", yaxis = "y2", visible = (i == 1),
    showlegend = FALSE, legendgroup = "raw_dens",
    hoverinfo = "text", text = paste("Raw Density (Log-Space):", round(dens_raw_log, 5))
  )
  
  # Trace 7: Log density (log space)
  p_explorer <- p_explorer %>% add_trace(
    x = y_grid_log, y = dens_log_log, type = "scatter", mode = "lines",
    line = list(color = "#dd6b20", width = 2), name = "Log Model Density",
    xaxis = "x2", yaxis = "y2", visible = (i == 1),
    showlegend = FALSE, legendgroup = "log_dens",
    hoverinfo = "text", text = paste("Log Density (Log-Space):", round(dens_log_log, 5))
  )
  
  # Trace 8: Test results (log space vertical lines)
  log_vals <- ifelse(obs$value == 0, -1, log(obs$value))
  x_tests_log <- rep(log_vals, each = 3) * c(1, 1, NA)
  y_tests_log <- rep(c(0, max_y_log * 1.05, NA), nrow(obs))
  p_explorer <- p_explorer %>% add_trace(
    x = x_tests_log, y = y_tests_log, type = "scatter", mode = "lines",
    line = list(color = "#718096", dash = "dash", width = 1), name = "Test Results",
    xaxis = "x2", yaxis = "y2", visible = (i == 1),
    showlegend = FALSE, legendgroup = "test_res",
    hoverinfo = "none"
  )
  
  # Trace 9: Inferred mean (log space)
  log_mu_raw <- ifelse(mu_raw <= 0, -1, log(mu_raw))
  p_explorer <- p_explorer %>% add_trace(
    x = c(log_mu_raw, log_mu_raw), y = c(0, max_y_log * 1.05), type = "scatter", mode = "lines",
    line = list(color = "#2f855a", width = 2), name = "Inferred Mean (\u03bc)",
    xaxis = "x2", yaxis = "y2", visible = (i == 1),
    showlegend = FALSE, legendgroup = "inf_mean",
    hoverinfo = "text", text = paste("Inferred Mean (Log-Space):", round(log_mu_raw, 2))
  )
  
  # Trace 10: Reference limit (log space)
  log_ref_high <- log(ref_high)
  p_explorer <- p_explorer %>% add_trace(
    x = c(log_ref_high, log_ref_high), y = c(0, max_y_log * 1.05), type = "scatter", mode = "lines",
    line = list(color = "#e53e3e", width = 2), name = "Reference Limit (X_high)",
    xaxis = "x2", yaxis = "y2", visible = (i == 1),
    showlegend = FALSE, legendgroup = "ref_limit",
    hoverinfo = "text", text = paste("Reference Limit (Log-Space):", round(log_ref_high, 2))
  )
}

# Construct the dropdown buttons
buttons <- list()
for (k in 1:M) {
  visible_vector <- rep(FALSE, 10 * M)
  visible_vector[(10 * (k - 1) + 1):(10 * k)] <- TRUE
  
  buttons[[k]] <- list(
    method = "restyle",
    args = list("visible", as.list(visible_vector)),
    label = analytes_list$sourceAnalyteName[k]
  )
}

# Define tick positions for custom log axis
raw_ticks <- c(0, 0.1, 0.5, 1, 5, 10, 50, 100, 500, 1000, 5000, 10000, 50000, 100000, 500000)
tickvals <- ifelse(raw_ticks == 0, -1, log(raw_ticks))
ticktext <- as.character(raw_ticks)

# Set layouts and configure dropdown
p_explorer <- p_explorer %>% layout(
  xaxis = list(domain = c(0, 0.47), title = "Raw Concentration (Linear Scale)"),
  yaxis = list(title = "Density"),
  xaxis2 = list(
    domain = c(0.53, 1), 
    title = "Log-transformed Concentration (Raw Value Ticks; 0 mapped to -1)",
    anchor = "x2",
    tickvals = tickvals,
    ticktext = ticktext
  ),
  yaxis2 = list(title = "Density", anchor = "x2"),
  updatemenus = list(
    list(
      buttons = buttons,
      direction = "down",
      showactive = TRUE,
      x = 0.05, y = 1.25,
      xanchor = "left",
      yanchor = "top"
    )
  ),
  margin = list(t = 100),
  legend = list(orientation = "h", y = -0.2)
)

p_explorer

5. Distribution of Model Fits (\(R^2\) Histograms)

To understand how well the reported Z-scores fit raw vs. logarithmic models across all analytes, we plot the distribution of their \(R^2\) values.

# Prepare data for R2 histograms
r2_data <- analyte_summary %>%
  filter(!is.na(fit_raw_r2) | !is.na(fit_log_r2)) %>%
  select(canonicalAnalyteId, sourceAnalyteName, fit_raw_r2, fit_log_r2) %>%
  pivot_longer(
    cols = c(fit_raw_r2, fit_log_r2),
    names_to = "model_type",
    values_to = "r2_value"
  ) %>%
  mutate(
    model_type = ifelse(model_type == "fit_raw_r2", "Raw Linear Model", "Log-Linear Model")
  )

# Plot interactive histograms using ggplot and plotly
p_hist <- ggplot(r2_data, aes(x = r2_value, fill = model_type)) +
  geom_histogram(alpha = 0.6, binwidth = 0.05, position = "identity", color = "white") +
  scale_fill_manual(values = c("Raw Linear Model" = "#3182ce", "Log-Linear Model" = "#dd6b20")) +
  theme_minimal() +
  labs(
    title = "Distribution of R² Fits Across Analytes",
    x = "R² Value",
    y = "Count of Analytes",
    fill = "Model Type"
  )

ggplotly(p_hist)

Analysis of Low \(R^2\) Values and Possible Explanations

In this platform audit, we notice that while many analytes have very high \(R^2\) values (close to 1.0, indicating a perfect linear or logarithmic mapping), some analytes exhibit low \(R^2\) fits for both models.

Here are the primary analytical explanations for low \(R^2\) values:

  1. Quantization and Precision Limits: For analytes with narrow reference ranges or low concentrations, the exported raw concentration values (\(X\)) might be rounded (quantized) to integers or low-decimal values. Since the Z-scores are exported with higher precision (or calculated using unrounded floats), plotting rounded \(X\) against precise \(Z\) introduces step-like patterns that degrade the correlation, producing a low \(R^2\).
  2. Standardization Engine Bugs/Inconsistencies: The platform’s scoring engine might compute Z-scores using reference intervals (\(\mu\) and \(\sigma\)) that are not constant across all samples. If different reference intervals were used for different patients (e.g. based on age, gender, or batch-specific values) but pooled together, a single regression line will not fit the combined data, resulting in a low \(R^2\).
  3. Outliers and Censored Data: Severe clinical outliers or values that are flagged as “below detection limit” (which might be filled with a constant or arbitrary value) disrupt the linear regression line, pulling down the \(R^2\) metric.
  4. Non-linear scoring profiles: Some analytes might use multi-segment linear interpolation, curvilinear ranges, or clinical threshold-based Z-score step functions rather than standard Gaussian mapping.

6. Platform Audit Summary Metrics

total_analytes <- nrow(analyte_summary)
evaluatable <- sum(!is.na(analyte_summary$fit_raw_r2))
raw_omitted_count <- sum(analyte_summary$raw_fit_dominates, na.rm = TRUE)
pct_omitted <- round((raw_omitted_count / total_analytes) * 100, 1)

median_floor <- median(analyte_summary$avg_z_floor, na.rm = TRUE)

Below are the key platform audit summary metrics recovered from the ingested data:

Executive Audit Summary Metrics
Metric Value
Total Analytes Evaluated 296
Analytes with Sufficient Data for Fit (N >= 2) 296
Analytes Fitting Raw Linear Model Better 217
Percentage of Analytes fitting Raw Model Better 73.3%
Median Theoretical Lower Z-Score Floor -0.724

Key Findings: - Raw Linear Dominance: 73.3 % of the evaluatable analytes have a better fit to a raw linear model than to a logarithmic model. This indicates the platform’s standardization engine is applying a standard linear Z-score mapping on raw concentration metrics. - Diagnostic Floor: The median theoretical lower floor is -0.724. At zero concentration, patients will have a suppressed lower boundary on their reported Z-score, demonstrating why standard raw-concentration standardization fails to provide fine metabolic deficiency resolution.

7. Interactive Concentration vs. Z-Score Explorer

Hover over any data point to inspect the analyte name, raw value, reported Z-score, derived population mean, upper reference range, and computed minimum Z-score floor.

# Prepare dataframe for plotting
plot_data <- df_processed %>%
  left_join(
    analyte_summary %>% select(canonicalAnalyteId, fit_raw_r2, fit_log_r2, raw_fit_dominates),
    by = "canonicalAnalyteId"
  ) %>%
  mutate(
    original_z = z_score,
    # Winsorize Z-score to [-4, +4]
    z_score = pmin(pmax(original_z, -4), 4),
    winsorized = !is.na(original_z) & (original_z < -4 | original_z > 4),
    
    # Log transform x axis. Transform 0 to -1.
    x_trans = ifelse(!is.na(value) & value == 0, -1, ifelse(!is.na(value) & value > 0, log(value), NA_real_)),
    
    hover_info = paste0(
      "<b>Analyte:</b> ", sourceAnalyteName, "<br>",
      "<b>Canonical ID:</b> ", canonicalAnalyteId, "<br>",
      "<b>Report File:</b> ", source_file, "<br>",
      "<b>Raw Value (X):</b> ", round(value, 3), "<br>",
      "<b>Log-transformed X:</b> ", round(x_trans, 3), "<br>",
      "<b>Reported Z-Score:</b> ", round(original_z, 2), 
        ifelse(winsorized, paste0(" (Winsorized to ", round(z_score, 2), ")"), ""), "<br>",
      "<b>Ref Range:</b> [0 - ", round(refHigh, 1), "]<br>",
      "<b>Derived Pop Mean (&mu;):</b> ", round(mu_derived, 1), "<br>",
      "<b>Derived Pop Std (&sigma;):</b> ", round(sigma_derived, 1), "<br>",
      "<b>Theoretical Floor (Z_min):</b> ", round(z_floor_derived, 2)
    ),
    z_score_status = ifelse(!is.na(z_score) & z_score < 0, "Negative", "Positive")
  )

p <- ggplot(plot_data, aes(x = x_trans, y = z_score, text = hover_info, color = z_score_status)) +
  geom_point(alpha = 0.7, size = 2.2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "#718096") +
  geom_hline(yintercept = median_floor, linetype = "dotted", color = "#e53e3e") +
  scale_color_manual(
    values = c("Negative" = "#e53e3e", "Positive" = "#3182ce"),
    labels = c("Negative" = "Negative Z-Score (Suppressed)", "Positive" = "Positive Z-Score")
  ) +
  theme_minimal() +
  labs(
    title = "Theriome Aristotle+: Log Concentration vs. Winsorized Z-Score",
    x = "Log-transformed Concentration Value (ln(X); 0 mapped to -1)",
    y = "Winsorized Z-Score (bounded at [-4, +4])",
    color = "Z-Score Status"
  )

ggplotly(p, tooltip = "text")

Winsorized Z-Score Log

The interactive table below records all observations that had their Z-scores truncated because they fell outside the [-4, +4] range:

8. Comprehensive 296-Analyte Reference & Floor Table

Sort and filter across all 296 analytes to inspect derived population metrics and theoretical lower-bound floors.

table_data <- analyte_summary %>%
  select(
    canonicalAnalyteId,
    sourceAnalyteName,
    n_obs,
    avg_refHigh,
    avg_mu_derived,
    avg_sigma_derived,
    avg_cv,
    avg_z_floor,
    fit_raw_r2,
    fit_log_r2
  ) %>%
  mutate(
    # Convert canonical ID to character to render text search filter in DT
    canonicalAnalyteId = as.character(canonicalAnalyteId),
    avg_refHigh = round(avg_refHigh, 1),
    avg_mu_derived = round(avg_mu_derived, 1),
    avg_sigma_derived = round(avg_sigma_derived, 1),
    avg_cv = round(avg_cv, 3),
    avg_z_floor = round(avg_z_floor, 2),
    fit_raw_r2 = round(fit_raw_r2, 4),
    fit_log_r2 = round(fit_log_r2, 4)
  )

datatable(
  table_data,
  options = list(
    pageLength = 15,
    searchHighlight = TRUE,
    autoWidth = TRUE
  ),
  filter = "top",
  colnames = c(
    "Canonical ID", "Analyte Name", "N Tests", "Ref High", 
    "Derived Mean (\u03bc)", "Derived Std (\u03c3)", "Derived CV", 
    "Z Floor (Z_min)", "R\u00b2 Raw", "R\u00b2 Log"
  )
)

9. Conclusions & Analytical Recommendations

  1. Universal Raw Gaussian Standardization: Across all tested files, reported Z-scores strictly obey \(Z = \frac{X - \mu}{\sigma}\) on raw concentrations, omitting log-normal transformation.
  2. Artificial Floor Bottleneck: Because concentrations cannot drop below zero (\(X \ge 0\)), every analyte is bounded by \(Z_{\text{min}} = -1/\text{CV}\). The median lower floor across the panel sits at \(Z = -0.72\).
  3. Actionable Fix: Transitioning the scoring engine to Inverse Normal Rank Transformation (INRT) or proper log-normal standardization will instantly restore diagnostic resolution for lower-tail metabolic deficiencies.