Come disegnare un grafico a barre diviso per livelli variabili, mentre si controlla per altre variabili tramite regressione multipla?
Come posso disegnare un grafico a barre per i mezzi, mentre controllo per altre variabili attraverso la regressione, in modo diviso per barre?
Il mio problema generale
Conduco una ricerca per capire quale frutto è più gradevole: mango, banana o mela. A tal fine, vado avanti e assaggio 100 persone a caso. Chiedo loro di valutare, su una scala da 1 a 5, il grado di gradimento di ciascuno dei frutti. Raccolgo anche alcune informazioni demografiche su di loro: sesso, età, livello di istruzione e se sono daltonici o meno perché penso che la visione dei colori potrebbe alterare i risultati. Ma il mio problema è che dopo la raccolta dei dati, mi rendo conto che il mio campione potrebbe non rappresentare bene la popolazione generale. Ho l'80% di maschi mentre nella popolazione il sesso è diviso in modo più equo. Il livello di istruzione nel mio campione è piuttosto uniforme, anche se nella popolazione è più comune possedere solo un diploma di scuola superiore che un dottorato di ricerca. Anche l'età non è rappresentativa.
Pertanto, è probabile che il solo calcolo dei mezzi per il gradimento della frutta basato sul mio campione sia limitato in termini di conclusioni generalizzate a livello di popolazione. Un modo per affrontare questo problema è eseguire una regressione multipla per controllare i dati demografici distorti.
Voglio tracciare i risultati delle regressioni in un grafico a barre, dove divido le barre (fianco a fianco) in base ai livelli di visione dei colori (daltonici o no).
I miei dati
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
Se voglio solo ottenere i valori medi per ogni livello di gradimento della frutta
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")
Ma cosa succede se voglio correggere questi mezzi per rappresentare meglio la popolazione?
Posso usare la regressione multipla
Correggerò l'età media della popolazione che è di 45 anni
Correggerò la corretta divisione 50-50 per il sesso
Correggerò per il livello di istruzione comune che è il liceo (codificato
2nei miei dati)Ho anche un motivo per credere che l'età influenzi il gusto della frutta in modo non lineare, quindi ne terrò conto.
lm(fruit ~ I(age - 45) + I((age - 45)^2) + I(is_male - 0.5) + I(education_level - 2)
Eseguirò i dati dei tre frutti (mela, banana, mango) attraverso lo stesso modello, estrarrò l'intercetta e la considererò come la media corretta dopo aver controllato i dati demografici.
Per prima cosa, eseguirò le regressioni sui dati solo con persone daltoniche
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
Quindi traccia il grafico a barre corretto solo per i daltonici
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")
In secondo luogo, ripeterò lo stesso processo di regressione sui dati solo con la visione dei colori
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")
Quello che cerco alla fine è combinare le trame corrette
Nota finale
Questa è principalmente una domanda su ggplote grafica. Tuttavia, come si può vedere, il mio metodo è lungo (cioè non conciso) e ripetitivo. Soprattutto rispetto alla semplicità di ottenere solo barplot per mezzi non corretti, come dimostrato all'inizio. Sarei molto felice se qualcuno avesse anche idee su come rendere il codice più breve e più semplice.
Risposte
Non sono convinto che stiate ottenendo le quantità statistiche desiderate quando adattate il modello ai sottoinsiemi di dati. Un modo migliore per porre le domande che si desidera porre sarebbe con un modello più completo (includere la cecità nel modello) e quindi calcolare i contrasti del modello per le differenze nel punteggio medio tra ciascun gruppo.
Detto questo, ecco un codice che fa quello che vuoi.
- Per prima cosa abbiamo
pivot_longerle colonne della frutta in modo che i tuoi dati siano in formato lungo. - Quindi abbiamo
group_byil tipo di frutta e le variabili di cecità e chiamiamonestche ci fornisce set di dati separati per ciascun tipo di frutta e categorie di cecità. - Quindi utilizziamo
purrr::mapper adattare un modello a ciascuno di questi set di dati. broom::tidyebroom::confint_tidyforniscici le statistiche che vogliamo per i modelli.- Quindi dobbiamo eliminare i riepiloghi del modello e filtrare in modo specifico le righe che corrispondono all'intercetta.
- Ora abbiamo i dati necessari per creare la figura, lascio a te il 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))
MODIFICARE
Nel caso in cui preferisci adattare un singolo modello (aumentando così la dimensione del campione e riducendo i valori delle tue stime). Puoi inserire is_colorblind nel modello come file 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))
Dovresti quindi ottenere previsioni per due osservazioni, la "persona media daltonica" e la "persona media che non è daltonica":
new_data <- expand_grid(age = 45, is_male = .5,
education_level = 2.5, is_colorblinded = c(0,1))
Potresti quindi fare come prima, adattando il nuovo modello con un po 'di programmazione funzionale, ma group_by(name)invece di namee 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 questo apporteresti una piccola modifica al vecchio codice di plottaggio:
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))
Quando confronti i due grafici, raggruppati e sottogruppi, noterai che gli intervalli di confidenza si restringono e le stime per le medie si avvicinano per lo più a 3. Questo sarebbe visto come un segno che stiamo andando un po 'meglio del modello sottogruppo , poiché conosciamo la verità fondamentale per quanto riguarda le distribuzioni campionate.