如何绘制按变量水平分割的条形图，同时通过多元回归控制其他变量？
How to draw a barplot split by variable levels, while controlling for other variables via multiple regression?
如何绘制均值条形图，同时通过回归控制其他变量——以按变量拆分条形的方式？
我的一般问题
我进行了一项研究，以确定哪种水果更受欢迎：芒果、香蕉或苹果。为此，我继续随机抽样 100 人。我要求他们以 1-5 的等级对每种水果的喜爱程度进行评分。我还收集了一些关于他们的人口统计信息：性别、年龄、教育水平，以及他们是否是色盲，因为我认为色觉可能会改变结果。 但我的问题是，在收集数据后，我意识到我的样本可能不能很好地代表一般人群。我有 80% 的男性，而在人群中，性别分布更为平均。我的样本中的教育水平非常统一，尽管在人口中只持有高中文凭比拥有博士学位更常见。年龄也不具有代表性。
因此，仅根据我的样本计算水果喜好的均值在将结论推广到人口水平方面可能会受到限制。解决此问题的一种方法是 运行 进行多元回归以控制有偏差的人口统计数据。
我想在条形图中绘制回归结果，我根据色觉水平（色盲与否）拆分条形图（并排）。
我的数据
library(tidyverse)
set.seed(123)
fruit_liking_df <-
data.frame(
id = 1:100,
i_love_apple = sample(c(1:5), 100, replace = TRUE),
i_love_banana = sample(c(1:5), 100, replace = TRUE),
i_love_mango = sample(c(1:5), 100, replace = TRUE),
age = sample(c(20:70), 100, replace = TRUE),
is_male = sample(c(0, 1), 100, prob = c(0.2, 0.8), replace = TRUE),
education_level = sample(c(1:4), 100, replace = TRUE),
is_colorblinded = sample(c(0, 1), 100, replace = TRUE)
)
> as_tibble(fruit_liking_df)
## # A tibble: 100 x 8
## id i_love_apple i_love_banana i_love_mango age is_male education_level is_colorblinded
## <int> <int> <int> <int> <int> <dbl> <int> <dbl>
## 1 1 3 5 2 50 1 2 0
## 2 2 3 3 1 49 1 1 0
## 3 3 2 1 5 70 1 1 1
## 4 4 2 2 5 41 1 3 1
## 5 5 3 1 1 49 1 4 0
## 6 6 5 2 1 29 0 1 0
## 7 7 4 5 5 35 1 3 0
## 8 8 1 3 5 24 0 3 0
## 9 9 2 4 2 55 1 2 0
## 10 10 3 4 2 69 1 4 0
## # ... with 90 more rows
如果我只想得到每个水果喜欢程度的平均值
fruit_liking_df_for_barplot <-
fruit_liking_df %>%
pivot_longer(.,
cols = c(i_love_apple, i_love_banana, i_love_mango),
names_to = "fruit",
values_to = "rating") %>%
select(id, fruit, rating, everything())
ggplot(fruit_liking_df_for_barplot, aes(fruit, rating, fill = as_factor(is_colorblinded))) +
stat_summary(fun = mean,
geom = "bar",
position = "dodge") +
## errorbars
stat_summary(fun.data = mean_se,
geom = "errorbar",
position = "dodge") +
## bar labels
stat_summary(
aes(label = round(..y.., 2)),
fun = mean,
geom = "text",
position = position_dodge(width = 1),
vjust = 2,
color = "white") +
scale_fill_discrete(name = "is colorblind?",
labels = c("not colorblind", "colorblind")) +
ggtitle("liking fruits, without correcting for demographics")
但是，如果我想更正这些方法以更好地代表总体呢？
我可以使用多元回归
我将校正人口的平均年龄 45
我会纠正正确的 50-50 性别比例
我将更正普通教育水平，即高中（在我的数据中编码为 2）
我也有理由相信年龄会以非线性方式影响对水果的喜爱，所以我也会考虑到这一点。
lm(fruit ~ I(age - 45) + I((age - 45)^2) + I(is_male - 0.5) + I(education_level - 2)
我将运行三种水果数据（苹果、香蕉、芒果）通过相同的模型，提取截距，并认为这是控制人口统计数据后的校正均值。
首先，我将 运行 仅对色盲人群进行数据回归
library(broom)
dep_vars <- c("i_love_apple",
"i_love_banana",
"i_love_mango")
regresults_only_colorblind <-
lapply(dep_vars, function(dv) {
tmplm <-
lm(
get(dv) ~ I(age - 45) + I((age - 45)^2) + I(is_male - 0.5) + I(education_level - 2),
data = filter(fruit_liking_df, is_colorblinded == 1)
)
broom::tidy(tmplm) %>%
slice(1) %>%
select(estimate, std.error)
})
data_for_corrected_barplot_only_colorblind <-
regresults_only_colorblind %>%
bind_rows %>%
rename(intercept = estimate) %>%
add_column(dep_vars, .before = c("intercept", "std.error"))
## # A tibble: 3 x 3
## dep_vars intercept std.error
## <chr> <dbl> <dbl>
## 1 i_love_apple 3.07 0.411
## 2 i_love_banana 2.97 0.533
## 3 i_love_mango 3.30 0.423
然后仅针对色盲绘制校正后的条形图
ggplot(data_for_corrected_barplot_only_colorblind,
aes(x = dep_vars, y = intercept)) +
geom_bar(stat = "identity", width = 0.7, fill = "firebrick3") +
geom_errorbar(aes(ymin = intercept - std.error, ymax = intercept + std.error),
width = 0.2) +
geom_text(aes(label=round(intercept, 2)), vjust=1.6, color="white", size=3.5) +
ggtitle("liking fruits after correction for demogrpahics \n colorblind subset only")
其次，我将对仅具有色觉的数据重复相同的回归过程
dep_vars <- c("i_love_apple",
"i_love_banana",
"i_love_mango")
regresults_only_colorvision <-
lapply(dep_vars, function(dv) {
tmplm <-
lm(
get(dv) ~ I(age - 45) + I((age - 45)^2) + I(is_male - 0.5) + I(education_level - 2),
data = filter(fruit_liking_df, is_colorblinded == 0) ## <- this is the important change here
)
broom::tidy(tmplm) %>%
slice(1) %>%
select(estimate, std.error)
})
data_for_corrected_barplot_only_colorvision <-
regresults_only_colorvision %>%
bind_rows %>%
rename(intercept = estimate) %>%
add_column(dep_vars, .before = c("intercept", "std.error"))
ggplot(data_for_corrected_barplot_only_colorvision,
aes(x = dep_vars, y = intercept)) +
geom_bar(stat = "identity", width = 0.7, fill = "orchid3") +
geom_errorbar(aes(ymin = intercept - std.error, ymax = intercept + std.error),
width = 0.2) +
geom_text(aes(label=round(intercept, 2)), vjust=1.6, color="white", size=3.5) +
ggtitle("liking fruits after correction for demogrpahics \n colorvision subset only")
我最终要寻找的是合并更正后的图
最后的注释
这主要是关于 ggplot 和图形的问题。但是，可以看出，我的方法很长（即不简洁）且重复。尤其是相对于仅获取未校正均值的条形图的简单性而言，如开头所示。如果有人对如何使代码更短更简单也有想法，我将非常高兴。
我不相信你在将模型拟合到数据子集时得到了你想要的统计量。提出您想问的问题的更好方法是使用更完整的模型（包括模型中的失明），然后计算 model contrasts 每组之间平均分的差异。
也就是说，这里有一些代码可以满足您的需求。
- 首先我们
pivot_longer 水果列，以便您的数据采用长格式。
- 然后我们
group_by 水果类型和盲目变量并调用 nest，这为我们提供了每种水果类型和盲目类别的单独数据集。
- 然后我们使用
purrr::map 为每个数据集拟合一个模型。
broom::tidy 和 broom::confint_tidy 为我们提供了我们想要的模型统计信息。- 然后我们需要解除模型摘要的嵌套，并专门向下过滤到与截距相对应的行。
- 我们现在已经有了创建图形所需的数据，剩下的就交给你了。
library(tidyverse)
set.seed(123)
fruit_liking_df <-
data.frame(
id = 1:100,
i_love_apple = sample(c(1:5), 100, replace = TRUE),
i_love_banana = sample(c(1:5), 100, replace = TRUE),
i_love_mango = sample(c(1:5), 100, replace = TRUE),
age = sample(c(20:70), 100, replace = TRUE),
is_male = sample(c(0, 1), 100, prob = c(0.2, 0.8), replace = TRUE),
education_level = sample(c(1:4), 100, replace = TRUE),
is_colorblinded = sample(c(0, 1), 100, replace = TRUE)
)
model_fits <- fruit_liking_df %>%
pivot_longer(starts_with("i_love"), values_to = "fruit") %>%
group_by(name, is_colorblinded) %>%
nest() %>%
mutate(model_fit = map(data, ~ lm(data = .x, fruit ~ I(age - 45) +
I((age - 45)^2) +
I(is_male - 0.5) +
I(education_level - 2))),
model_summary = map(model_fit, ~ bind_cols(broom::tidy(.x), broom::confint_tidy(.x))))
model_fits %>%
unnest(model_summary) %>%
filter(term == "(Intercept)") %>%
ggplot(aes(x = name, y = estimate, group = is_colorblinded,
fill = as_factor(is_colorblinded), colour = as_factor(is_colorblinded))) +
geom_bar(stat = "identity", position = position_dodge(width = .95)) +
geom_errorbar(stat = "identity", aes(ymin = conf.low, ymax = conf.high),
colour = "black", width = .15, position = position_dodge(width = .95))
编辑
在您更愿意拟合单个模型的情况下（从而增加样本量并减少估计值）。您可以将 is_colorblind 作为 factor.
拉入模型
lm(data = .x, fruit ~ I(age - 45) +
I((age - 45)^2) + I(is_male - 0.5) +
I(education_level - 2) +
as.factor(is_colorblind))
然后您可能想要获得对“色盲的普通人”和“非色盲的普通人”两个观察结果的预测：
new_data <- expand_grid(age = 45, is_male = .5,
education_level = 2.5, is_colorblinded = c(0,1))
然后您可以像以前一样，用一些函数式编程来拟合新模型，但是 group_by(name) 而不是 name 和 is_colorblind。
model_fits_ungrouped <- fruit_liking_df %>%
pivot_longer(starts_with("i_love"), values_to = "fruit") %>%
group_by(name) %>%
tidyr::nest() %>%
mutate(model_fit = map(data, ~ lm(data = .x, fruit ~ I(age - 45) +
I((age - 45)^2) +
I(is_male - .5) +
I(education_level - 2) +
as.factor(is_colorblinded))),
predicted_values = map(model_fit, ~ bind_cols(new_data,
as.data.frame(predict(newdata = new_data, .x,
type = "response", se.fit = T))) %>%
rowwise() %>%
mutate(estimate = fit,
conf.low = fit - qt(.975, df) * se.fit,
conf.high = fit + qt(.975, df) * se.fit)))
用这个你可以对旧的绘图代码做一个小的改动：
model_fits_ungrouped %>%
unnest(predicted_values) %>%
ggplot(aes(x = name, y = estimate, group = is_colorblinded,
fill = as_factor(is_colorblinded), colour = as_factor(is_colorblinded))) +
geom_bar(stat = "identity", position = position_dodge(width = .95)) +
geom_errorbar(stat = "identity", aes(ymin = conf.low, ymax = conf.high),
colour = "black", width = .15, position = position_dodge(width = .95))
当您比较分组和子分组的两个图时，您会注意到置信区间缩小并且均值的估计值大多接近 3。这将被视为我们正在做一些工作的标志比分组模型更好，因为我们知道关于采样分布的基本事实。
如何绘制均值条形图，同时通过回归控制其他变量——以按变量拆分条形的方式？
我的一般问题
我进行了一项研究，以确定哪种水果更受欢迎：芒果、香蕉或苹果。为此，我继续随机抽样 100 人。我要求他们以 1-5 的等级对每种水果的喜爱程度进行评分。我还收集了一些关于他们的人口统计信息：性别、年龄、教育水平，以及他们是否是色盲，因为我认为色觉可能会改变结果。 但我的问题是，在收集数据后，我意识到我的样本可能不能很好地代表一般人群。我有 80% 的男性，而在人群中，性别分布更为平均。我的样本中的教育水平非常统一，尽管在人口中只持有高中文凭比拥有博士学位更常见。年龄也不具有代表性。
因此，仅根据我的样本计算水果喜好的均值在将结论推广到人口水平方面可能会受到限制。解决此问题的一种方法是 运行 进行多元回归以控制有偏差的人口统计数据。
我想在条形图中绘制回归结果，我根据色觉水平（色盲与否）拆分条形图（并排）。
我的数据
library(tidyverse)
set.seed(123)
fruit_liking_df <-
data.frame(
id = 1:100,
i_love_apple = sample(c(1:5), 100, replace = TRUE),
i_love_banana = sample(c(1:5), 100, replace = TRUE),
i_love_mango = sample(c(1:5), 100, replace = TRUE),
age = sample(c(20:70), 100, replace = TRUE),
is_male = sample(c(0, 1), 100, prob = c(0.2, 0.8), replace = TRUE),
education_level = sample(c(1:4), 100, replace = TRUE),
is_colorblinded = sample(c(0, 1), 100, replace = TRUE)
)
> as_tibble(fruit_liking_df)
## # A tibble: 100 x 8
## id i_love_apple i_love_banana i_love_mango age is_male education_level is_colorblinded
## <int> <int> <int> <int> <int> <dbl> <int> <dbl>
## 1 1 3 5 2 50 1 2 0
## 2 2 3 3 1 49 1 1 0
## 3 3 2 1 5 70 1 1 1
## 4 4 2 2 5 41 1 3 1
## 5 5 3 1 1 49 1 4 0
## 6 6 5 2 1 29 0 1 0
## 7 7 4 5 5 35 1 3 0
## 8 8 1 3 5 24 0 3 0
## 9 9 2 4 2 55 1 2 0
## 10 10 3 4 2 69 1 4 0
## # ... with 90 more rows
如果我只想得到每个水果喜欢程度的平均值
fruit_liking_df_for_barplot <-
fruit_liking_df %>%
pivot_longer(.,
cols = c(i_love_apple, i_love_banana, i_love_mango),
names_to = "fruit",
values_to = "rating") %>%
select(id, fruit, rating, everything())
ggplot(fruit_liking_df_for_barplot, aes(fruit, rating, fill = as_factor(is_colorblinded))) +
stat_summary(fun = mean,
geom = "bar",
position = "dodge") +
## errorbars
stat_summary(fun.data = mean_se,
geom = "errorbar",
position = "dodge") +
## bar labels
stat_summary(
aes(label = round(..y.., 2)),
fun = mean,
geom = "text",
position = position_dodge(width = 1),
vjust = 2,
color = "white") +
scale_fill_discrete(name = "is colorblind?",
labels = c("not colorblind", "colorblind")) +
ggtitle("liking fruits, without correcting for demographics")
但是，如果我想更正这些方法以更好地代表总体呢？
我可以使用多元回归
我将校正人口的平均年龄 45
我会纠正正确的 50-50 性别比例
我将更正普通教育水平，即高中（在我的数据中编码为
2）我也有理由相信年龄会以非线性方式影响对水果的喜爱，所以我也会考虑到这一点。
lm(fruit ~ I(age - 45) + I((age - 45)^2) + I(is_male - 0.5) + I(education_level - 2)
我将运行三种水果数据（苹果、香蕉、芒果）通过相同的模型，提取截距，并认为这是控制人口统计数据后的校正均值。
首先，我将 运行 仅对色盲人群进行数据回归
library(broom)
dep_vars <- c("i_love_apple",
"i_love_banana",
"i_love_mango")
regresults_only_colorblind <-
lapply(dep_vars, function(dv) {
tmplm <-
lm(
get(dv) ~ I(age - 45) + I((age - 45)^2) + I(is_male - 0.5) + I(education_level - 2),
data = filter(fruit_liking_df, is_colorblinded == 1)
)
broom::tidy(tmplm) %>%
slice(1) %>%
select(estimate, std.error)
})
data_for_corrected_barplot_only_colorblind <-
regresults_only_colorblind %>%
bind_rows %>%
rename(intercept = estimate) %>%
add_column(dep_vars, .before = c("intercept", "std.error"))
## # A tibble: 3 x 3
## dep_vars intercept std.error
## <chr> <dbl> <dbl>
## 1 i_love_apple 3.07 0.411
## 2 i_love_banana 2.97 0.533
## 3 i_love_mango 3.30 0.423
然后仅针对色盲绘制校正后的条形图
ggplot(data_for_corrected_barplot_only_colorblind,
aes(x = dep_vars, y = intercept)) +
geom_bar(stat = "identity", width = 0.7, fill = "firebrick3") +
geom_errorbar(aes(ymin = intercept - std.error, ymax = intercept + std.error),
width = 0.2) +
geom_text(aes(label=round(intercept, 2)), vjust=1.6, color="white", size=3.5) +
ggtitle("liking fruits after correction for demogrpahics \n colorblind subset only")
其次，我将对仅具有色觉的数据重复相同的回归过程
dep_vars <- c("i_love_apple",
"i_love_banana",
"i_love_mango")
regresults_only_colorvision <-
lapply(dep_vars, function(dv) {
tmplm <-
lm(
get(dv) ~ I(age - 45) + I((age - 45)^2) + I(is_male - 0.5) + I(education_level - 2),
data = filter(fruit_liking_df, is_colorblinded == 0) ## <- this is the important change here
)
broom::tidy(tmplm) %>%
slice(1) %>%
select(estimate, std.error)
})
data_for_corrected_barplot_only_colorvision <-
regresults_only_colorvision %>%
bind_rows %>%
rename(intercept = estimate) %>%
add_column(dep_vars, .before = c("intercept", "std.error"))
ggplot(data_for_corrected_barplot_only_colorvision,
aes(x = dep_vars, y = intercept)) +
geom_bar(stat = "identity", width = 0.7, fill = "orchid3") +
geom_errorbar(aes(ymin = intercept - std.error, ymax = intercept + std.error),
width = 0.2) +
geom_text(aes(label=round(intercept, 2)), vjust=1.6, color="white", size=3.5) +
ggtitle("liking fruits after correction for demogrpahics \n colorvision subset only")
我最终要寻找的是合并更正后的图
最后的注释
这主要是关于 ggplot 和图形的问题。但是，可以看出，我的方法很长（即不简洁）且重复。尤其是相对于仅获取未校正均值的条形图的简单性而言，如开头所示。如果有人对如何使代码更短更简单也有想法，我将非常高兴。
我不相信你在将模型拟合到数据子集时得到了你想要的统计量。提出您想问的问题的更好方法是使用更完整的模型（包括模型中的失明），然后计算 model contrasts 每组之间平均分的差异。
也就是说，这里有一些代码可以满足您的需求。
- 首先我们
pivot_longer水果列，以便您的数据采用长格式。 - 然后我们
group_by水果类型和盲目变量并调用nest，这为我们提供了每种水果类型和盲目类别的单独数据集。 - 然后我们使用
purrr::map为每个数据集拟合一个模型。 broom::tidy和broom::confint_tidy为我们提供了我们想要的模型统计信息。- 然后我们需要解除模型摘要的嵌套，并专门向下过滤到与截距相对应的行。
- 我们现在已经有了创建图形所需的数据，剩下的就交给你了。
library(tidyverse)
set.seed(123)
fruit_liking_df <-
data.frame(
id = 1:100,
i_love_apple = sample(c(1:5), 100, replace = TRUE),
i_love_banana = sample(c(1:5), 100, replace = TRUE),
i_love_mango = sample(c(1:5), 100, replace = TRUE),
age = sample(c(20:70), 100, replace = TRUE),
is_male = sample(c(0, 1), 100, prob = c(0.2, 0.8), replace = TRUE),
education_level = sample(c(1:4), 100, replace = TRUE),
is_colorblinded = sample(c(0, 1), 100, replace = TRUE)
)
model_fits <- fruit_liking_df %>%
pivot_longer(starts_with("i_love"), values_to = "fruit") %>%
group_by(name, is_colorblinded) %>%
nest() %>%
mutate(model_fit = map(data, ~ lm(data = .x, fruit ~ I(age - 45) +
I((age - 45)^2) +
I(is_male - 0.5) +
I(education_level - 2))),
model_summary = map(model_fit, ~ bind_cols(broom::tidy(.x), broom::confint_tidy(.x))))
model_fits %>%
unnest(model_summary) %>%
filter(term == "(Intercept)") %>%
ggplot(aes(x = name, y = estimate, group = is_colorblinded,
fill = as_factor(is_colorblinded), colour = as_factor(is_colorblinded))) +
geom_bar(stat = "identity", position = position_dodge(width = .95)) +
geom_errorbar(stat = "identity", aes(ymin = conf.low, ymax = conf.high),
colour = "black", width = .15, position = position_dodge(width = .95))
编辑
在您更愿意拟合单个模型的情况下（从而增加样本量并减少估计值）。您可以将 is_colorblind 作为 factor.
lm(data = .x, fruit ~ I(age - 45) +
I((age - 45)^2) + I(is_male - 0.5) +
I(education_level - 2) +
as.factor(is_colorblind))
然后您可能想要获得对“色盲的普通人”和“非色盲的普通人”两个观察结果的预测：
new_data <- expand_grid(age = 45, is_male = .5,
education_level = 2.5, is_colorblinded = c(0,1))
然后您可以像以前一样，用一些函数式编程来拟合新模型，但是 group_by(name) 而不是 name 和 is_colorblind。
model_fits_ungrouped <- fruit_liking_df %>%
pivot_longer(starts_with("i_love"), values_to = "fruit") %>%
group_by(name) %>%
tidyr::nest() %>%
mutate(model_fit = map(data, ~ lm(data = .x, fruit ~ I(age - 45) +
I((age - 45)^2) +
I(is_male - .5) +
I(education_level - 2) +
as.factor(is_colorblinded))),
predicted_values = map(model_fit, ~ bind_cols(new_data,
as.data.frame(predict(newdata = new_data, .x,
type = "response", se.fit = T))) %>%
rowwise() %>%
mutate(estimate = fit,
conf.low = fit - qt(.975, df) * se.fit,
conf.high = fit + qt(.975, df) * se.fit)))
用这个你可以对旧的绘图代码做一个小的改动：
model_fits_ungrouped %>%
unnest(predicted_values) %>%
ggplot(aes(x = name, y = estimate, group = is_colorblinded,
fill = as_factor(is_colorblinded), colour = as_factor(is_colorblinded))) +
geom_bar(stat = "identity", position = position_dodge(width = .95)) +
geom_errorbar(stat = "identity", aes(ymin = conf.low, ymax = conf.high),
colour = "black", width = .15, position = position_dodge(width = .95))
当您比较分组和子分组的两个图时，您会注意到置信区间缩小并且均值的估计值大多接近 3。这将被视为我们正在做一些工作的标志比分组模型更好，因为我们知道关于采样分布的基本事实。