重回帰を介して他の変数を制御しながら、変数レベルで分割されたバープロットを描画するにはどうすればよいですか？
Aug 20 2020
回帰によって他の変数を制御しながら、平均のバープロットを描画するにはどうすればよいですか？バーごとに分割されますか？
私の一般的な問題
マンゴー、バナナ、リンゴのどれがより好感が持てるのかを調べるために調査を行っています。この目的のために、私は先に進み、ランダムに100人をサンプリングします。それぞれの果物の好みの程度を1〜5のスケールで評価してもらいます。また、性別、年齢、教育レベル、色覚異常が結果を変える可能性があるため、色覚異常かどうかなど、それらに関する人口統計情報も収集します。しかし、私の問題は、データ収集後、私のサンプルが一般的な母集団をうまく表していない可能性があることに気付いたことです。私は80％の男性がいますが、人口の中で性別はより均等に分かれています。私のサンプルの教育レベルはかなり均一ですが、人口では博士号を取得するよりも高校の卒業証書のみを取得する方が一般的です。年齢も代表的ではありません。
したがって、私のサンプルに基づいて果物の好みの平均を計算するだけでは、結論を人口レベルに一般化するという点で制限される可能性があります。この問題に対処する1つの方法は、重回帰を実行して、偏った人口統計データを制御することです。
回帰の結果をバープロットにプロットしたいと思います。ここでは、色覚レベル（色覚異常かどうか）に従ってバーを（並べて）分割します。
私のデータ
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)
3つの果物データ（リンゴ、バナナ、マンゴー）を同じモデルで実行し、切片を抽出し、人口統計データを制御した後の修正平均と見なします。
まず、色覚異常の人だけでデータの回帰を実行します
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グラフィックスに関する質問です。ただし、ご覧のとおり、私の方法は長く（つまり、簡潔ではありません）、反復的です。特に、最初に示したように、修正されていない平均のバープロットを取得するという単純さに関連します。コードを短く簡単にする方法について誰かがアイデアを持っていれば、私はとても幸せです。
回答
1 BrianLang Aug 20 2020 at 16:37
モデルをデータサブセットに適合させたときに必要な統計量が得られるとは確信していません。質問したい質問をするためのより良い方法は、より完全なモデル（モデルに失明を含む）を使用してから、各グループ間の平均スコアの違いについてモデルのコントラストを計算することです。
そうは言っても、ここにあなたが望むことをするいくつかのコードがあります。
- まず
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))
編集
単一のモデルに適合させたい場合（したがって、サンプルサイズを増やし、推定値のseを減らします）。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))
次に、「色覚異常の平均的な人」と「色覚異常のない平均的な人」の2つの観測値の予測を取得する必要があります。
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))
グループ化とサブグループ化された2つのプロットを比較すると、信頼区間が縮小し、平均の推定値がほぼ3に近づくことがわかります。これは、サブグループ化されたモデルよりも少し優れていることを示しています。 、サンプリングされた分布に関するグラウンドトゥルースがわかっているため。