我想使用 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 步。从模型系数创建数据框:
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
第 3 步。合并两个数据框。现在所有引用因子水平都存在，我们有一列指定术语:
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()
第 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"))
关于r - ggplot 中的森林图，引用水平来自回归模型，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/49690498/