使用回归模型中的参考水平在ggplot中制作森林图

3
我想使用ggplot2包制作森林图,并且我对我的输出感到满意(见下面的森林图)。
该图显示了回归模型中给定变量的水平(比率比和置信区间),以及参考水平。
问题是生成图需要大量的手动劳动。
第一个问题是,我希望参考水平在图中跟随给定变量的其他水平,所以我手动输入了每个这样的参考水平(请参见下表)。为使ggplot2正常工作,我输入了任意负的比率比和置信区间值作为参考水平,然后将绘图限制设置为从零到较大的正数范围。
第二个问题是,因为我的原始变量位于单个列中,所以我手动输入颜色,这很耗时。
是否有更简单的方法生成这样的图?任何帮助都将不胜感激。
# DATA 
mtcars
mtcars$gear <- as.factor(mtcars$gear)
mtcars$carb <- as.factor(mtcars$carb)

# PREPARE ODDS RATIO & CONFIDENCE INTERVALS DATA FRAME 
model = lm(mpg ~ gear + carb + disp, data = mtcars ) # make regression model
forest_table = data.frame(
  or= round(exp(coef(model)),2), 
  round(exp(confint(model, level = 0.95)),2), 
  check.names = F) # make a table with odds ratio and confidence intervals
names(forest_table) = c("or", "ci_lb", "ci_ub") # give columns clear names
library(data.table)
setDT(forest_table, keep.rownames = TRUE)[] # turn row names into a column
forest_table <- as.data.frame(forest_table) # turn table into a data frame
forest_table <- forest_table[-1, ] # get rid of the intercept row

# ADD ROWS WITH REFERENCE LEVELS TO PREPARED DATA FRAME
r <- 2 # row after which new row is to be inserted
newrow <- c("3 reference", -10.00, -9.00, -11.00) # row to be inserted 
forest_table <- rbind(forest_table[1:r, ], newrow, forest_table[-(1:r), ]) # insert row
r <- 8 # row after which new row is to be inserted
newrow <- c("1 reference", -10.00, -9.00, -11.00) # row to be inserted 
forest_table <- rbind(forest_table[1:r, ], newrow, forest_table[-(1:r), ]) # insert row

# FIX CLASSES IN PREPARED DATA FRAME 
forest_table$or <- as.numeric(forest_table$or)
forest_table$ci_lb <- as.numeric(forest_table$ci_lb)
forest_table$ci_ub <- as.numeric(forest_table$ci_ub)

# ADD DUMMY VARIABLE TO CONTROL ORDER IN PLOT 
forest_table$order <- as.factor(rep(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10))) # create dummy variable 
forest_table$order <- factor(forest_table$order, 
                             levels = rev(levels(forest_table$order))) 
# use dummy variable to counteract ggplot2 default of reversing the order of levels in 
# the prepared data frame when plotting  

# PLOT
library(ggplot2)
forestplot <- ggplot(forest_table, aes(or, order)) + 
  geom_point(size = 5, shape = 18, aes(colour = order)) + # data points
  geom_errorbarh(aes(xmax = ci_ub, xmin = ci_lb, colour = order), 
                 height = 0.15) + # error bars
  geom_vline(xintercept = 1, linetype = "longdash") + # line marking 0 on x axis
  scale_x_continuous(breaks = seq(0, 40000, 10000), 
                     labels = seq(0, 40000, 10000),
                     limits = c(0, 50000)) + # x axis scale and labels
 scale_colour_manual(values = c("blue", "red", "red", "red", "red", "red", "red", 
                                "green", "green", "green")) # manually set one colour per variable 
1个回答

0

您可以在表格中添加一个额外的列,指示模型中的实际术语,并根据此进行颜色分配。以下是一个示例:

步骤1。从模型系数创建数据框:

new.table <- data.frame(
  coef = names(coef(model)),
  or = round(exp(coef(model)), 2),
  ci_lb = round(exp(confint(model, level = 0.95)), 2)[, 1],
  ci_ub = round(exp(confint(model, level = 0.95)), 2)[, 2],
  stringsAsFactors = FALSE, row.names = NULL
)

> new.table
         coef           or        ci_lb        ci_ub
1 (Intercept) 1.226831e+11 249767693.18 6.026058e+13
2       gear4 5.396000e+01         0.31 9.403190e+03
3       gear5 2.193800e+02         1.03 4.662360e+04
4       carb2 1.400000e-01         0.00 4.340000e+00
5       carb3 2.000000e-02         0.00 1.280000e+00
6       carb4 0.000000e+00         0.00 2.000000e-01
7       carb6 0.000000e+00         0.00 3.700000e-01
8       carb8 0.000000e+00         0.00 2.100000e-01
9        disp 9.800000e-01         0.96 1.000000e+00

步骤2。请注意,在“lm”模型中,model$xlevels包含包括多个因子水平的术语的信息。

> model$xlevels
$gear
[1] "3" "4" "5"

$carb
[1] "1" "2" "3" "4" "6" "8"

这可以用于创建包含所有因子水平的参考数据框:

library(dplyr)
library(data.table)

terms.with.levels <- names(model$xlevels)
df.with.levels <- lapply(terms.with.levels, 
       function(x) data.frame(term = x,
                              coef = paste0(x, model$xlevels[[x]]),
                              stringsAsFactors = FALSE)) %>%
  rbindlist()

> df.with.levels
   term  coef
1: gear gear3
2: gear gear4
3: gear gear5
4: carb carb1
5: carb carb2
6: carb carb3
7: carb carb4
8: carb carb6
9: carb carb8

步骤三。合并这两个数据框。现在所有的参考因子水平都已经存在,并且我们有一列指定了术语:

new.table <- merge(new.table, df.with.levels, all = TRUE)

> new.table
          coef           or        ci_lb        ci_ub term
1  (Intercept) 1.226831e+11 249767693.18 6.026058e+13 <NA>
2        carb1           NA           NA           NA carb
3        carb2 1.400000e-01         0.00 4.340000e+00 carb
4        carb3 2.000000e-02         0.00 1.280000e+00 carb
5        carb4 0.000000e+00         0.00 2.000000e-01 carb
6        carb6 0.000000e+00         0.00 3.700000e-01 carb
7        carb8 0.000000e+00         0.00 2.100000e-01 carb
8         disp 9.800000e-01         0.96 1.000000e+00 <NA>
9        gear3           NA           NA           NA gear
10       gear4 5.396000e+01         0.31 9.403190e+03 gear
11       gear5 2.193800e+02         1.03 4.662360e+04 gear

步骤4。进一步修改数据框:

new.table <- new.table %>%

  # drop intercept
  filter(coef != "(Intercept)") %>%

  # indicate whether each row is for a reference level
  mutate(is.reference = is.na(or)) %>%

  # for non-factor term in the model (e.g. disp) which
  # have a single coefficient, term == coef
  mutate(term = ifelse(is.na(term), coef, term)) %>%

  # set reference levels' x values to 0
  mutate_at(vars(or, ci_lb, ci_ub),
            funs(ifelse(is.reference, 0, .))) %>%

  # order terms according to the model specifications
  mutate(term = factor(term,
                       levels = attr(model$terms, "term.labels")))

> new.table
    coef     or ci_lb    ci_ub term is.reference
1  carb1   0.00  0.00     0.00 carb         TRUE
2  carb2   0.14  0.00     4.34 carb        FALSE
3  carb3   0.02  0.00     1.28 carb        FALSE
4  carb4   0.00  0.00     0.20 carb        FALSE
5  carb6   0.00  0.00     0.37 carb        FALSE
6  carb8   0.00  0.00     0.21 carb        FALSE
7   disp   0.98  0.96     1.00 disp        FALSE
8  gear3   0.00  0.00     0.00 gear         TRUE
9  gear4  53.96  0.31  9403.19 gear        FALSE
10 gear5 219.38  1.03 46623.60 gear        FALSE

步骤5。创建图表。通过将参考水平的alpha设置为0(即100%透明度)可以隐藏它们,并且系数的顺序由外观的顺序控制:

p <- ggplot(new.table,
       aes(x = or, xmin = ci_lb, xmax = ci_ub,
           y = coef, color = term, alpha = !is.reference)) +
  geom_vline(xintercept = 1, linetype = "longdash") +
  geom_errorbarh(height = 0.15) +
  geom_point(size = 5, shape = 18) +
  facet_grid(term~., scales = "free_y", space = "free_y") +
  scale_alpha_identity()

basic plot

步骤6。如有需要,进一步调整图表:

p +
  # specify colour for each term in the model
  scale_color_manual(values = c("gear" = "green",
                                "carb" = "red",
                                "disp" = "blue")) +

  # hide facet labels
  theme(strip.background = element_blank(),
        strip.text = element_blank()) +

  # remove spacing between facets
  theme(panel.spacing = unit(0, "pt"))

tweaked plot


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接