如何在循环中从lmerTest模型提取交互项的p值？

问题背景 #

使用lmerTest包拟合多个线性混合模型后，尝试通过循环提取genderfemale:time（二阶交互）和genderfemale:time:[解释变量]（三阶交互）的p值时，出现错误：

Error in coef(summary(model))[, 5] : subscript out of bounds

原循环代码及其中一个模型的summary输出如下：

原代码 #

library(lmerTest)

# Create a list of models with interaction terms to loop over
models <- list(
  mixed_age_interaction,
  mixed_tnfi_year_interaction,
  mixed_crp_interaction
)

# Create a list of explanatory variables to loop over
explanatoryVariables <- list(
  "age_at_diagnosis",
  "bio_drug_start_year",
  "crp"
)

loop_function <- function(models, explanatoryVariables) {
  # Create an empty data frame to store the results
  coef_df <- data.frame(adj_coef_gender_sex = numeric(), coef_interaction_term = numeric(), explanatory_variable = character(), adj_coef_pvalue = numeric())
  
  # Loop over the models and explanatory variables
  for (i in seq_along(models)) {
    model <- models[[i]]
    explanatoryVariable <- explanatoryVariables[[i]]
    
    # Extract the adjusted coefficients for the gender*time interaction
    adj_coef <- fixef(model)["genderfemale:time"]
    
    # Extract the fixed effect of the interaction term
    interaction_coef <- fixef(model)[paste0("genderfemale:time:", explanatoryVariable)]
    
    # Extract the p-value for the adjusted coefficient for gender*time
    adj_coef_pvalue <- coef(summary(model))[,5]["genderfemale:time"]
    
    # Add a row to the data frame with the results for this model
    coef_df <- bind_rows(coef_df, data.frame(adj_coef_gender_sex = adj_coef, coef_interaction_term = interaction_coef, explanatory_variable = explanatoryVariable, adj_coef_pvalue = adj_coef_pvalue))
  }
  return(coef_df)
}

# Loop over the models and extract the fixed effects
coef_df <- loop_function(models, explanatoryVariables)
coef_df

示例模型summary #

Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: basdai ~ 1 + gender + time + age_at_diagnosis + gender * time +  
    time * age_at_diagnosis + gender * age_at_diagnosis + gender *  
    time * age_at_diagnosis + (1 | ID) + (1 | country)
   Data: dat

      AIC       BIC    logLik  deviance  df.resid 
 254340.9  254431.8 -127159.5  254318.9     28557 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.3170 -0.6463 -0.0233  0.6092  4.3180 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 154.62   12.434  
 country  (Intercept)  32.44    5.695  
 Residual             316.74   17.797  
Number of obs: 28568, groups:  ID, 11207; country, 13

Fixed effects:
                                     Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)                         4.669e+01  1.792e+00  2.082e+01  26.048  < 2e-16 ***
genderfemale                        2.368e+00  1.308e+00  1.999e+04   1.810   0.0703 .  
time                               -1.451e+01  4.220e-01  2.164e+04 -34.382  < 2e-16 ***
age_at_diagnosis                    9.907e-02  2.220e-02  1.963e+04   4.463 8.12e-06 ***
genderfemale:time                   1.431e-01  7.391e-01  2.262e+04   0.194   0.8464    
time:age_at_diagnosis               8.188e-02  1.172e-02  2.185e+04   6.986 2.90e-12 ***
genderfemale:age_at_diagnosis       8.547e-02  3.453e-02  2.006e+04   2.476   0.0133 *  
genderfemale:time:age_at_diagnosis  4.852e-03  1.967e-02  2.274e+04   0.247   0.8052    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) gndrfm time   ag_t_d gndrf: tm:g__ gnd:__
genderfemal -0.280                                          
time        -0.241  0.331                                   
age_t_dgnss -0.434  0.587  0.511                            
gendrfml:tm  0.139 -0.519 -0.570 -0.293                     
tm:g_t_dgns  0.228 -0.313 -0.951 -0.533  0.543              
gndrfml:g__  0.276 -0.953 -0.329 -0.639  0.495  0.343        
gndrfml::__ -0.137  0.491  0.567  0.319 -0.954 -0.596 -0.516

错误原因 #

1. 列索引不稳定：coef(summary(model))的列数/列位置可能因模型拟合方式（如是否使用Satterthwaite自由度）变化，直接用[,5]取p值列不可靠。
2. 项存在性未检查：若某个模型中不存在目标交互项，直接索引会导致越界错误。

修正方案 #

关键修改点 #

• 用列名Pr(>|t|)替代固定索引提取p值，保证兼容性
• 先检查目标项是否存在于模型固定效应中，避免报错
• 扩展结果数据框，同时存储二阶和三阶交互的p值

修正后的完整代码 #

library(lmerTest)
library(dplyr)

# 模型列表
models <- list(
  mixed_age_interaction,
  mixed_tnfi_year_interaction,
  mixed_crp_interaction
)

# 解释变量列表
explanatoryVariables <- list(
  "age_at_diagnosis",
  "bio_drug_start_year",
  "crp"
)

loop_function <- function(models, explanatoryVariables) {
  # 初始化结果数据框，新增两个p值列
  coef_df <- data.frame(
    adj_coef_gender_sex = numeric(),
    coef_interaction_term = numeric(),
    explanatory_variable = character(),
    gender_time_pvalue = numeric(),
    three_way_pvalue = numeric(),
    stringsAsFactors = FALSE
  )
  
  for (i in seq_along(models)) {
    model <- models[[i]]
    expl_var <- explanatoryVariables[[i]]
    
    # 获取模型的固定效应系数表（含p值）
    model_coef_summary <- coef(summary(model))
    # 获取所有固定效应项名
    fixed_terms <- rownames(model_coef_summary)
    
    # 定义目标项名
    two_way_term <- "genderfemale:time"
    three_way_term <- paste0("genderfemale:time:", expl_var)
    
    # 提取二阶交互的系数和p值（不存在则设为NA）
    adj_coef <- if(two_way_term %in% fixed_terms) fixef(model)[two_way_term] else NA
    two_way_p <- if(two_way_term %in% fixed_terms) model_coef_summary[two_way_term, "Pr(>|t|)"] else NA
    
    # 提取三阶交互的系数和p值（不存在则设为NA）
    three_way_coef <- if(three_way_term %in% fixed_terms) fixef(model)[three_way_term] else NA
    three_way_p <- if(three_way_term %in% fixed_terms) model_coef_summary[three_way_term, "Pr(>|t|)"] else NA
    
    # 新增行到结果框
    coef_df <- bind_rows(coef_df, data.frame(
      adj_coef_gender_sex = adj_coef,
      coef_interaction_term = three_way_coef,
      explanatory_variable = expl_var,
      gender_time_pvalue = two_way_p,
      three_way_pvalue = three_way_p,
      stringsAsFactors = FALSE
    ))
  }
  return(coef_df)
}

# 运行循环提取结果
coef_df <- loop_function(models, explanatoryVariables)
print(coef_df)

说明 #

• 代码中加入了项存在性检查，若模型中没有目标交互项，对应值会设为NA，避免报错
• 使用列名Pr(>|t|)提取p值，无论模型输出列顺序如何都能正确定位
• 结果数据框新增了gender_time_pvalue（二阶交互p值）和three_way_pvalue（三阶交互p值）两列，完整保存所需信息

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