如何将GLM模型系数正确绑定至汇总数据表？

Got it, let's walk through exactly how to do this with the mtcars dataset—this approach works for any GLM model and summary dataset you're working with. Here's a step-by-step breakdown:

Step 1: Build Your GLM Model #

First, let's create a sample GLM (we'll use linear regression here, a special case of GLM, for simplicity) using am (transmission type: 0 = automatic, 1 = manual) as our categorical predictor and mpg as the response variable. This gives us a clear reference group (am=0) to handle.

# Load required library for data manipulation
library(dplyr)

# Build the GLM model
model <- glm(mpg ~ factor(am), data = mtcars)

# Check the model summary to confirm coefficients
summary(model)

Step 2: Extract Model Metrics & Prepare Benchmark Values #

Next, we'll pull the coefficients, standard errors, and p-values from the model. We also need to set up the reference group (am=0) with a coefficient of 1, using the intercept's standard error and p-value as you requested.

# Extract coefficient table from model summary
model_coef_table <- as.data.frame(summary(model)$coefficients)
colnames(model_coef_table) <- c("Estimate", "Std.Error", "t_stat", "p_value")

# Grab intercept's SE and p-value for the reference group
intercept_se <- model_coef_table["(Intercept)", "Std.Error"]
intercept_p <- model_coef_table["(Intercept)", "p_value"]

# Create a row for the reference group (am=0)
reference_row <- data.frame(
  Estimate = 1,
  Std.Error = intercept_se,
  t_stat = NA,  # No t-stat for reference group, mark as NA
  p_value = intercept_p,
  am = 0  # Match the reference group's value
)

# Clean up the model's non-intercept coefficients and add the am value
model_coef_clean <- model_coef_table[-1, ] %>%
  mutate(am = as.integer(sub("factor\\(am\\)(\\d)", "\\1", rownames(.)))) %>%
  rownames_to_column(var = "predictor") %>%
  select(am, Estimate, Std.Error, t_stat, p_value)

# Combine reference group and model coefficients
full_coef_data <- bind_rows(reference_row, model_coef_clean)

Step 3: Create Your Summary Dataset #

Let's make a sample summary dataset for mtcars—we'll group by am and calculate average horsepower and weight, but you can replace this with your actual summary data.

# Build summary dataset (customize this to match your data!)
summary_dataset <- mtcars %>%
  group_by(am) %>%
  summarise(
    avg_hp = mean(hp),
    avg_weight = mean(wt),
    .groups = "drop"
  )

Step 4: Merge Everything Together #

Finally, we'll join the summary dataset with our coefficient data to get the final merged table.

# Merge summary data with coefficient metrics
final_merged_data <- left_join(summary_dataset, full_coef_data, by = "am") %>%
  rename(
    coefficient = Estimate,
    standard_error = Std.Error
  ) %>%
  select(am, avg_hp, avg_weight, coefficient, standard_error, p_value)

# View the final result
print(final_merged_data)

What the Final Output Looks Like #

You'll end up with a table where:

• For the reference group (am=0), coefficient is 1, standard_error and p_value match the model's intercept values
• For the non-reference group (am=1), you get the actual coefficient, SE, and p-value from the model
• All metrics are tied directly to your summary dataset's groups

If you're working with a more complex model (like logistic regression, or multiple categorical predictors), you just need to adjust how you extract and map the reference groups—this core framework stays the same.

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