如何在lm_robust回归输出中自动显示因子变量的参考水平？

Got it! Here's a straightforward way to automatically add reference level rows for factor variables to your lm_robust output—perfect for grabbing those reference group names and values for plotting. I’ve put together a custom function that handles this, and I’ll walk you through using it with your example data.

Step 1: Define the helper function #

This function will scan your model for factor variables, identify their reference levels, and insert the reference row into the coefficient table:

library(estimatr)
library(tibble)
library(dplyr)

add_reference_levels <- function(model) {
  # Pull model metadata to analyze variables
  model_terms <- terms(model)
  model_frame <- model.frame(model)
  
  # Convert the model's coefficient table to a tibble for easy manipulation
  coef_table <- as_tibble(model$coefficients, rownames = "term")
  
  # Loop through each term to check for factor variables
  for (term in attr(model_terms, "term.labels")) {
    var_name <- all.vars(as.formula(paste0("~", term)))[1]
    if (is.factor(model_frame[[var_name]])) {
      # Get the reference level (R's default is the first factor level)
      ref_level <- levels(model_frame[[var_name]])[1]
      ref_term_label <- paste0(var_name, ref_level)
      
      # Skip if the reference term is already in the table (shouldn't happen)
      if (!ref_term_label %in% coef_table$term) {
        # Build the reference row: Estimate = 0, other metrics set to NA
        ref_row <- tibble(
          term = ref_term_label,
          Estimate = 0,
          `Std. Error` = NA_real_,
          `t value` = NA_real_,
          `Pr(>|t|)` = NA_real_,
          `CI Lower` = NA_real_,
          `CI Upper` = NA_real_,
          DF = model$df_residual
        )
        
        # Insert the reference row right before the first non-reference term for this variable
        term_positions <- grep(paste0("^", var_name), coef_table$term)
        if (length(term_positions) > 0) {
          coef_table <- coef_table %>%
            add_row(!!!ref_row, .before = term_positions[1])
        } else {
          coef_table <- bind_rows(coef_table, ref_row)
        }
      }
    }
  }
  
  return(coef_table)
}

Step 2: Test with your example data #

Let’s run this with your sample code (I added set.seed() for reproducibility):

set.seed(123)
N = 20000
x = rbinom(N, 1, prob = 0.4)
y = 0.4*x + rnorm(N)
df <- data.frame(x,y)
df$x <- factor(df$x)

# Fit the robust linear model
lm_model <- lm_robust(df, formula = y ~ x)

# Add reference levels to the output
lm_model_with_ref <- add_reference_levels(lm_model)

# Print the full result
print(lm_model_with_ref, n = Inf)

Expected Output #

You’ll get a coefficient table that includes the reference level (x0 here) with an estimate of 0, making it easy to use for plotting:

# A tibble: 3 × 8
  term        Estimate `Std. Error` `t value` `Pr(>|t|)` `CI Lower` `CI Upper`    DF
  <chr>          <dbl>        <dbl>     <dbl>       <dbl>      <dbl>      <dbl> <dbl>
1 (Intercept)   0.0032        0.0092      0.35       0.73     -0.015       0.021 19998
2 x0             0           NA          NA          NA        NA          NA    19998
3 x1             0.397        0.0145     27.4        0        0.369       0.425 19998

Quick Notes #

• The function uses the first level of your factor as the reference (R’s default). If you changed the reference level with relevel(), it will pick up that updated level automatically.
• If you want to label the reference row differently (e.g., replace NA with "Reference"), you can adjust the ref_row definition—just be mindful of data types (you might need to convert columns to character if mixing strings and numbers).

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