Bagaimana cara menggambar barplot yang dipisahkan oleh level variabel, sambil mengontrol variabel lain melalui regresi berganda?
Bagaimana cara menggambar barplot untuk sarana, sambil mengontrol variabel lain melalui regresi - dengan cara split-bar-by-vars?
Masalah umum saya
Saya melakukan penelitian untuk mencari tahu buah mana yang lebih disukai: mangga, pisang, atau apel. Untuk tujuan ini, saya melanjutkan dan mengambil sampel 100 orang secara acak. Saya meminta mereka untuk menilai, pada skala 1-5, tingkat kesukaan setiap buah. Saya juga mengumpulkan beberapa informasi demografis tentang mereka: jenis kelamin, usia, tingkat pendidikan, dan apakah mereka buta warna atau tidak karena menurut saya penglihatan warna dapat mengubah hasil. Tetapi masalah saya adalah setelah pengumpulan data, saya menyadari bahwa sampel saya mungkin tidak mewakili populasi umum dengan baik. Saya memiliki 80% laki-laki sementara dalam populasi, jenis kelamin terbagi lebih merata. Tingkat pendidikan dalam sampel saya cukup seragam, meskipun dalam populasi lebih umum hanya memiliki ijazah sekolah menengah daripada memiliki gelar PhD. Umur juga tidak mewakili.
Oleh karena itu, hanya cara menghitung untuk menyukai buah berdasarkan sampel saya kemungkinan besar akan dibatasi dalam hal menggeneralisasikan kesimpulan ke tingkat populasi. Salah satu cara untuk mengatasi masalah ini adalah dengan menjalankan regresi berganda untuk mengontrol data demografi yang bias.
Saya ingin memplot hasil regresi dalam barplot, di mana saya membagi batang (berdampingan) sesuai dengan tingkat penglihatan warna (buta warna atau tidak).
Data saya
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
Jika saya hanya ingin mendapatkan nilai rata-rata untuk setiap tingkat kesukaan buah
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")
Tetapi bagaimana jika saya ingin mengoreksi cara-cara ini untuk mewakili populasi dengan lebih baik?
Saya bisa menggunakan regresi berganda
Saya akan mengoreksi usia rata-rata dalam populasi yaitu 45
Saya akan mengoreksi pembagian 50-50 yang benar untuk seks
Saya akan mengoreksi untuk tingkat pendidikan umum yaitu sekolah menengah (dikodekan
2dalam data saya)Saya juga memiliki alasan untuk percaya bahwa usia memengaruhi rasa suka buah-buahan secara non-linear, jadi saya akan menjelaskannya juga.
lm(fruit ~ I(age - 45) + I((age - 45)^2) + I(is_male - 0.5) + I(education_level - 2)
Saya akan menjalankan data tiga buah (apel, pisang, mangga) melalui model yang sama, mengekstrak intersep, dan menganggapnya sebagai mean yang dikoreksi setelah mengontrol data demografi.
Pertama, saya akan menjalankan regresi pada data hanya dengan orang buta warna
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
Kemudian plot barplot yang dikoreksi hanya untuk buta warna
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")
Kedua, saya akan mengulangi proses regresi yang sama pada data dengan penglihatan warna saja
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")
Yang akhirnya saya cari adalah menggabungkan plot yang dikoreksi
Catatan terakhir
Ini terutama pertanyaan tentang ggplotdan grafik. Namun, seperti yang bisa dilihat, metode saya panjang (tidak ringkas) dan berulang. Terutama terkait dengan kesederhanaan hanya mendapatkan barplot untuk cara yang tidak dikoreksi, seperti yang ditunjukkan di awal. Saya akan sangat senang jika seseorang juga memiliki ide tentang cara membuat kode lebih pendek dan sederhana.
Jawaban
Saya tidak yakin Anda mendapatkan kuantitas statistik yang Anda inginkan saat menyesuaikan model pada subset data. Cara yang lebih baik untuk mengajukan pertanyaan yang ingin Anda tanyakan adalah dengan model yang lebih lengkap (termasuk kebutaan dalam model) dan kemudian menghitung kontras model untuk perbedaan skor rata-rata antara masing-masing kelompok.
Karena itu, berikut adalah beberapa kode yang melakukan apa yang Anda inginkan.
- Pertama kita
pivot_longerkolom buah agar data Anda dalam format panjang. - Kemudian kita
group_byvariabel jenis buah dan kebutaan dan panggilannestyang memberi kita kumpulan data terpisah untuk setiap jenis buah dan kategori kebutaan. - Kami kemudian menggunakan
purrr::mapuntuk menyesuaikan model dengan masing-masing set data tersebut. broom::tidydanbroom::confint_tidyberikan statistik yang kami inginkan untuk model tersebut.- Kemudian kita perlu menghapus ringkasan model dan memfilter secara khusus ke baris yang sesuai dengan intersep.
- Kami sekarang memiliki data yang kami butuhkan untuk membuat gambar, saya serahkan sisanya kepada Anda.
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))
EDIT
Dalam kasus di mana Anda lebih suka menyesuaikan model tunggal (sehingga meningkatkan ukuran sampel dan mengurangi perkiraan Anda). Anda dapat menarik is_colorblind ke dalam model sebagai 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))
Anda kemudian ingin mendapatkan prediksi untuk dua pengamatan, "orang biasa yang buta warna" dan "orang biasa yang tidak buta warna":
new_data <- expand_grid(age = 45, is_male = .5,
education_level = 2.5, is_colorblinded = c(0,1))
Anda kemudian dapat melakukan seperti sebelumnya, menyesuaikan model baru dengan beberapa pemrograman fungsional, tetapi group_by(name)alih-alih namedan 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)))
Dengan ini, Anda akan membuat perubahan kecil pada kode plot lama:
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))
Saat Anda membandingkan dua plot, yang dikelompokkan dan disubkelompokkan, Anda akan melihat bahwa interval kepercayaan menyusut dan perkiraan rata-rata mendekati 3. Ini akan terlihat sebagai tanda bahwa kami melakukan sedikit lebih baik daripada model subkelompok , karena kita mengetahui kebenaran dasar sehubungan dengan distribusi sampel.