使用sjPlot::tab_model将对数变换的线性混合效应模型系数转为百分比变化

问题背景 #

我拟合了对数变换的线性混合效应模型，希望用sjPlot::tab_model生成汇总表时，将对数尺度的线性系数转换为百分比变化（而非保留原对数尺度）。手动转换系数的逻辑已经实现，但使用tab_model的transform参数时触发错误：

Error in if (fam.info$is_linear) transform <- NULL else transform <- "exp" : argument is of length zero

模拟数据代码 #

library(lme4)
library(lmerTest)
library(dplyr)

set.seed(1234)
dat_short <- data.frame(
  dv = c(
    # t1 ctrl组数据
    rnorm(mean=0.8, sd=0.1, n=6), #Long Healthy
    rnorm(mean=0.7, sd=0.1, n=6), #Lat Healthy
    rnorm(mean=0.6, sd=0.1, n=4),  #Long Damaged
    rnorm(mean=0.5, sd=0.1, n=4),  #Lat Damaged
    # t2 ctrl组数据
    rnorm(mean=0.7, sd=0.1, n=6), #Long Healthy
    rnorm(mean=0.6, sd=0.1, n=6), #Lat Healthy
    rnorm(mean=0.5, sd=0.1, n=4),  #Long Damaged
    rnorm(mean=0.4, sd=0.1, n=4),  #Lat Damaged
    # t1 trt组数据
    rnorm(mean=0.8, sd=0.1, n=6), #Long Healthy
    rnorm(mean=0.7, sd=0.1, n=6), #Lat Healthy
    rnorm(mean=0.6, sd=0.15, n=4), #Long Damaged
    rnorm(mean=0.5, sd=0.15, n=4), #Lat Damaged
    # t2 trt组数据
    rnorm(mean=0.7, sd=0.1, n=6), #Long Healthy
    rnorm(mean=0.6, sd=0.1, n=6), #Lat Healthy
    rnorm(mean=0.65, sd=0.15, n=4),#Long Damaged
    rnorm(mean=0.55, sd=0.15, n=4) #Lat Damaged
  ),
  id=c(rep(c("subj_1", "subj_2"), times=c(40, 40))),
  intervention=c(rep(c("ctrl", "trt"), times=c(40, 40))), 
  timepoint=c(rep(rep(c("t1", "t2"), times=c(20, 20)),2)), 
  direction=c(rep(rep(c("long", "lat", "long", "lat"), times=c(6, 6, 4, 4)),4)),
  region=c(rep(rep(c("healthy", "damaged"), times=c(12, 8)),4))
)  |>
  mutate(dv = case_when(
         id == "subj_1" ~ dv + runif(1, min = 0.01, max = 0.2),
         id == "subj_2" ~ dv))

speed_measures <-data.frame(
  n_speed = c(
    # t1 ctrl组速度数据
    round(runif(8, min = 3, max =10),0)
  ),
  id=c(rep(c("subj_1", "subj_2"), times=c(4, 4))),
  timepoint=c(rep(rep(c("t1", "t2"), times=c(2, 2)),2)), 
  direction=c(rep(rep(c("long", "lat"), times=c(1, 1)))
))  

dat_short_combined <- speed_measures |> left_join(dat_short) |> slice_sample(n = 70)

模型拟合代码 #

lmm_1_short <- lmer(dv ~ intervention*timepoint*region + direction + (1|id), data=dat_short)

手动转换系数代码 #

# 将回归系数转换为百分比变化
lmm_1_short_perc_summary <- coef(lmm_1_short)$id |>
  summarise(across(.fns = ~100*( exp(.x)-1) ) ) |>
  summarise(across(.fns = ~mean(.x))) |>
  select(-"(Intercept)") |>
  tidyr::pivot_longer(cols = everything(), names_to="Coefficient", values_to = "% Change") |>
  # 加入回归结果中的p值等信息
  cbind(coef(summary(lmm_1_short))[-1,-1]) |>
  # 计算置信区间并转换为百分比变化
  cbind(data.frame(confint(lmm_1_short))[4:11,]) |>
  transmute(Coefficient, across(where(is.numeric), .fns = ~round(.x,3))) |>
  rename("2.5% CI" = `X2.5..`, "97.5% CI" = `X97.5..`) 

row.names(lmm_1_short_perc_summary) <- 1:nrow(lmm_1_short_perc_summary)
lmm_1_short_perc_summary

错误原因与正确实现方法 #

错误原因 #

sjPlot默认将lmer拟合的线性混合模型识别为线性模型，会自动忽略transform参数（将其设为NULL），导致自定义转换函数无法执行，从而触发上述错误。

正确实现方式 #

通过force.transform = TRUE强制开启转换功能，同时定义百分比转换函数传递给transform参数，tab_model会自动转换系数和置信区间，直接生成sjPlot风格的汇总表：

# 定义百分比转换函数：对数系数转百分比变化
perc_transform <- function(x) 100 * (exp(x) - 1)

# 生成带百分比转换的sjPlot汇总表
sjPlot::tab_model(lmm_1_short,
                  show.intercept = FALSE,
                  show.reflvl = TRUE,
                  show.obs = TRUE,
                  df.method = "satterthwaite",
                  transform = perc_transform,  # 应用转换函数
                  force.transform = TRUE,      # 强制线性模型执行转换
                  # 可选：自定义变量标签提升可读性
                  pred.labels = c("干预：trt", "时间点：t2", "区域：damaged", 
                                  "干预×时间点", "干预×区域", "时间点×区域", 
                                  "干预×时间点×区域", "方向：lat"),
                  dv.labels = "因变量"
)

这个方法无需手动合并结果，直接生成符合要求的格式化表格，同时保留sjPlot默认的样式和统计信息展示。

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