¿Cómo dibujar una gráfica de barras dividida por niveles de variable, mientras se controlan otras variables a través de regresión múltiple?
¿Cómo puedo dibujar una gráfica de barras para las medias, mientras controlo otras variables a través de la regresión, en forma de barras divididas por vars?
Mi problema general
Realizo una investigación para averiguar qué fruta es más agradable: mango, plátano o manzana. Con este fin, sigo adelante y hago una muestra de 100 personas al azar. Les pido que califiquen, en una escala del 1 al 5, el grado de agrado de cada una de las frutas. También recopilo información demográfica sobre ellos: género, edad, nivel de educación y si son daltónicos o no porque creo que la visión del color podría alterar los resultados. Pero mi problema es que después de la recopilación de datos, me doy cuenta de que mi muestra podría no representar bien a la población general. Tengo un 80% de hombres, mientras que en la población el sexo está dividido de manera más uniforme. El nivel de educación en mi muestra es bastante uniforme, aunque en la población es más común tener solo un diploma de escuela secundaria que un doctorado. La edad tampoco es representativa.
Por lo tanto, es probable que el simple cálculo de las medias para el gusto por la fruta basado en mi muestra sea limitado en términos de generalizar las conclusiones al nivel de la población. Una forma de abordar este problema es ejecutar una regresión múltiple para controlar los datos demográficos sesgados.
Quiero trazar los resultados de la (s) regresión (es) en un diagrama de barras, donde divido las barras (una al lado de la otra) de acuerdo con los niveles de visión del color (daltónico o no).
Mis datos
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
Si solo quiero obtener los valores medios para cada nivel de gusto por la fruta
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")
Pero, ¿y si quiero corregir estos medios para representar mejor a la población?
Puedo usar regresión múltiple
Corregiré por la edad promedio de la población que es de 45
Corregiré la división correcta 50-50 para el sexo.
Corregiré el nivel de educación común que es la escuela secundaria (codificado
2en mis datos)También tengo una razón para creer que la edad afecta el gusto por las frutas de una manera no lineal, así que también lo explicaré.
lm(fruit ~ I(age - 45) + I((age - 45)^2) + I(is_male - 0.5) + I(education_level - 2)
Ejecutaré los datos de las tres frutas (manzana, plátano, mango) a través del mismo modelo, extraeré la intersección y lo consideraré como la media corregida después de controlar los datos demográficos.
Primero, ejecutaré las regresiones en datos solo con personas daltónicas
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
Luego trazar el diagrama de barras corregido solo para daltónicos
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")
En segundo lugar, repetiré el mismo proceso de regresión en datos con visión de color únicamente.
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")
Lo que busco en última instancia es combinar las parcelas corregidas
Nota final
Esta es principalmente una pregunta sobre ggplotgráficos. Sin embargo, como puede verse, mi método es largo (es decir, no conciso) y repetitivo. Especialmente en relación con la simplicidad de simplemente obtener un gráfico de barras para medios no corregidos, como se demostró al principio. Estaré muy feliz si alguien también tiene ideas sobre cómo hacer que el código sea más corto y simple.
Respuestas
No estoy convencido de que esté obteniendo las cantidades estadísticas que desea cuando ajusta el modelo a los subconjuntos de datos. Una mejor manera de hacer las preguntas que desea hacer sería con un modelo más completo (incluir ceguera en el modelo) y luego calcular los contrastes del modelo para las diferencias en la puntuación media entre cada grupo.
Dicho esto, aquí hay un código que hace lo que quieres.
- Primero ponemos
pivot_longerlas columnas de frutas para que sus datos estén en formato largo. - Luego, seleccionamos
group_byel tipo de fruta y las variables de ceguera y llamamos, lonestque nos da conjuntos de datos separados para cada tipo de fruta y categorías de ceguera. - Luego usamos
purrr::mappara ajustar un modelo a cada uno de esos conjuntos de datos. broom::tidyybroom::confint_tidydanos las estadísticas que queremos para los modelos.- Luego, necesitamos desanidar los resúmenes del modelo y filtrar específicamente a las filas que corresponden a la intersección.
- Ahora tenemos los datos que necesitamos para crear la figura, te dejo el resto.
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))
EDITAR
En el caso de que prefiera ajustar un solo modelo (aumentando así el tamaño de la muestra y reduciendo el conjunto de sus estimaciones). Puede introducir is_colorblind en el modelo como 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))
Luego, querrá obtener predicciones para dos observaciones, la "persona promedio que es daltónica" y la "persona promedio que no es daltónica":
new_data <- expand_grid(age = 45, is_male = .5,
education_level = 2.5, is_colorblinded = c(0,1))
Luego puede hacer lo mismo que antes, ajustando el nuevo modelo con algo de programación funcional, pero en group_by(name)lugar de namey 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)))
Con esto, haría un cambio menor al antiguo código de trazado:
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))
Cuando compare los dos gráficos, agrupados y subgrupos, notará que los intervalos de confianza se reducen y las estimaciones de las medias se acercan en su mayoría a 3. Esto se vería como una señal de que lo estamos haciendo un poco mejor que el modelo subgrupado , ya que conocemos la verdad básica con respecto a las distribuciones muestreadas.