基于group_by()实现分机台Passing-Bablock回归分析需求

Got it, let's streamline your Passing-Bablock regression workflow so you don't have to manually filter each machine one by one. Here's how to use dplyr grouping and purrr to run these regressions in batches for both parameter pairs across all 4 machines:

Step 1: Clean Up Data Preparation #

First, let's simplify your data setup—no need for repeated mutate_at calls since we can define the data as numeric from the start:

# Load required packages
require(mcr)
require(dplyr)
require(purrr)

# Create and format the dataset
df_reg <- tibble(
  MACHINE = rep(c("M1", "M2", "M3", "M4"), each = 10),
  P1 = c(2.09,3.71,4.71,4.30,4.45,3.29,1.96,3.01,3.33,3.30,
         2.06,3.81,4.53,4.55,4.62,3.51,2.01,3.08,3.53,3.48,
         2.06,3.74,4.60,4.41,4.46,3.37,1.95,3.00,3.25,3.32,
         2.06,3.78,4.55,4.29,4.47,3.29,1.91,2.92,3.23,3.22),
  P2 = c(6.16,4.28,4.62,1.21,7.09,7.82,10.84,3.50,2.98,2.70,
         6.15,4.45,4.77,1.90,7.94,8.07,11.13,3.65,3.33,2.78,
         6.05,4.34,4.60,1.85,7.89,8.10,11.72,3.79,3.20,2.91,
         6.16,4.47,4.66,1.85,7.78,7.96,11.16,3.63,3.24,2.88),
  P1_ref = rep(c(1.935,3.735,4.765,4.505,4.680,3.385,1.885,2.935,3.350,3.290), 4),
  P2_ref = rep(c(6.255,4.700,4.945,2.055,8.430,9.130,12.155,3.990,3.715,3.285), 4)
)

Step 2: Build a Reusable Regression Function #

Create a helper function that runs the Passing-Bablock regression and extracts key summary stats (you can tweak this to include more details if needed):

run_passing_bablock <- function(data, ref_col, test_col) {
  # Run the regression with bootstrap CIs
  model <- mcreg(
    x = data[[ref_col]],
    y = data[[test_col]],
    method.reg = "PaBa",
    method.ci = "bootstrap",
    boot.n = 1000,
    digits = 3
  )
  
  # Extract critical results into a tidy tibble
  tibble(
    parameter_pair = paste(test_col, "vs", ref_col),
    intercept = model$para[1],
    intercept_ci_lower = model$para.ci[1,1],
    intercept_ci_upper = model$para.ci[1,2],
    slope = model$para[2],
    slope_ci_lower = model$para.ci[2,1],
    slope_ci_upper = model$para.ci[2,2],
    r_squared = model$rsq
  )
}

Step 3: Batch Process with Grouping #

Use group_by() + nest() to group data by machine, then map() to apply the regression to both parameter pairs for each group:

# Run all regressions in one go
batch_results <- df_reg %>%
  group_by(MACHINE) %>%
  nest() %>% # Nest data into groups per machine
  mutate(
    # Run regression for P1 vs P1_ref
    p1_results = map(data, ~run_passing_bablock(.x, ref_col = "P1_ref", test_col = "P1")),
    # Run regression for P2 vs P2_ref
    p2_results = map(data, ~run_passing_bablock(.x, ref_col = "P2_ref", test_col = "P2"))
  ) %>%
  unnest(c(p1_results, p2_results)) %>% # Unnest results into a flat table
  select(-data) %>% # Remove nested data column
  ungroup()

# View full results
print(batch_results, n = Inf)

Optional: Keep Full Model Objects for Plotting #

If you want to generate diagnostic plots for each regression, modify the code to store the full mcreg objects:

# Store full model objects instead of just summaries
batch_models <- df_reg %>%
  group_by(MACHINE) %>%
  nest() %>%
  mutate(
    p1_model = map(data, ~mcreg(x = .x$P1_ref, y = .x$P1, method.reg = "PaBa")),
    p2_model = map(data, ~mcreg(x = .x$P2_ref, y = .x$P2, method.reg = "PaBa"))
  )

# Example: Plot M1's P1 vs P1_ref regression
plot(batch_models$p1_model[[1]], main = "M1: P1 vs Reference P1 (Passing-Bablock)")

Key Tips: #

• The bootstrap method for confidence intervals uses 1000 iterations here—adjust boot.n based on your computational power and precision needs.
• The helper function can be extended to include additional outputs from the mcreg object (like p-values or bootstrap summary stats) if you need them.

内容的提问来源于stack exchange，提问作者user4686711