Как нарисовать гистограмму, разделенную по уровням переменных, контролируя другие переменные с помощью множественной регрессии?
Как я могу нарисовать гистограмму для средних значений, контролируя другие переменные с помощью регрессии - по принципу разделения столбцов по переменным?
Моя общая проблема
Я провожу исследование, чтобы выяснить, какой фрукт мне больше нравится: манго, банан или яблоко. С этой целью я произвольно выбираю 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графику. Однако, как видно, мой метод длинный (т.е. не сжатый) и повторяющийся. Особенно по сравнению с простотой получения гистограммы для неисправленных средних значений, как показано в начале. Я буду очень рад, если у кого-то появятся идеи, как сделать код короче и проще.
Ответы
Я не уверен, что вы получаете нужные статистические величины при подборе модели на подмножествах данных. Лучший способ задать вопросы, которые вы хотите задать, - это использовать более полную модель (включая слепоту в модели), а затем вычислить контрасты моделей для различий в средних оценках между каждой группой.
При этом вот код, который делает то, что вы хотите.
- Сначала мы создаем
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. Это будет рассматриваться как признак того, что мы делаем немного лучше, чем модель с подгруппами. , так как мы знаем основную истину в отношении выборочных распределений.