使用R语言metafor包绘制森林图：指数化率显示为0的问题求助

在使用R语言metafor包绘制森林图时，对数转换后的发病率（含合并效应量）可正常展示，但尝试将效应量转换为实际发病率（指数化后）时，图中所有数值均显示为0。实际发病率范围为11至500，需按地区分组展示不同年龄组。

模型拟合代码 #

使用对数转换后的发病率拟合随机效应模型：

result1 <- rma(yi = log_rate, sei = se_log_rate, data = df_aggregated_year_grp1, method = "REML")

此模型对应的对数尺度森林图可正常显示。

尝试的森林图代码 #

第一次尝试设置指数转换与x轴范围：

forest(result1,
       atransf = exp, # 对模型估计值应用指数转换
       refline = log(1), # 无效应参考线
       xlab = "Incidence Rate per 100,000", # x轴标签
       xlim = c(0, max(df_sub$exp_upper_ci) * 1.1), # 包含最大置信上限的x轴范围
       slab = paste(df_sub$study_ID_years, df_sub$age_group, sep = ": ") # 研究标签
)

调整后问题仍存在，第二次尝试代码：

forest(result1, 
       atransf=exp,  # 应用指数转换
       xlab="Incidence Rate per 100,000",  # 更新x轴标签
       xlim=c(0, max(df_aggregated_year_grp1$exp_upper_ci)),  # 适配指数化后数值范围的x轴限制
       slab=paste(df_aggregated_year_grp1$study_ID_years, df_aggregated_year_grp1$age_group, sep=": "),
       cex=0.7) 

指数化后数据范围 #

summary(df_aggregated_year_grp1$exp_rate)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  11.44   45.10   82.70  144.61  152.78  539.53 

summary(df_aggregated_year_grp1$exp_lower_ci)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  7.689  33.510  49.464  93.501  93.739 387.612 

summary(df_aggregated_year_grp1$exp_upper_ci)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  17.02   72.32  133.07  236.24  341.38  750.99 

解决思路 #

• 修正x轴范围计算：确保xlim使用指数化后的置信上限最大值，同时处理可能存在的NA值，修改为：

  xlim = c(0, max(df_aggregated_year_grp1$exp_upper_ci, na.rm = TRUE) * 1.1)
  

• 手动指定刻度点：metafor默认刻度基于原对数尺度，需手动设置对数刻度点，转换后显示实际发病率：

  forest(result1, 
         atransf = exp,
         at = log(c(10, 50, 100, 200, 500, 750)), # 对数尺度的刻度点
         xlab = "Incidence Rate per 100,000",
         xlim = c(0, max(df_aggregated_year_grp1$exp_upper_ci, na.rm = TRUE) * 1.1),
         slab = paste(df_aggregated_year_grp1$study_ID_years, df_aggregated_year_grp1$age_group, sep = ": "),
         cex = 0.7)
  

• 手动替换轴标签：先绘制对数尺度森林图，再覆盖x轴标签：

  # 绘制对数尺度森林图，隐藏原x轴标签
  forest(result1, 
         refline = log(1),
         xlab = "", # 先清空标签
         xlim = log(c(1, max(df_aggregated_year_grp1$exp_upper_ci, na.rm = TRUE) * 1.1)),
         slab = paste(df_aggregated_year_grp1$study_ID_years, df_aggregated_year_grp1$age_group, sep = ": "),
         cex = 0.7)
  # 移除原x轴刻度
  axis(side = 1, labels = FALSE)
  # 添加自定义刻度与实际发病率标签
  at_log <- log(c(10, 50, 100, 200, 500, 750))
  at_labels <- c("10", "50", "100", "200", "500", "750")
  axis(side = 1, at = at_log, labels = at_labels)
  # 添加x轴标题
  mtext("Incidence Rate per 100,000", side = 1, line = 2.5, cex = 0.8)
  

• 验证模型结果：通过print(result1)查看模型输出，确认指数化后的合并效应量（exp(result1$b)）在数据范围内，排除模型拟合错误。

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